Why is the Kruskal Wallis test statistic approximately chi-square distributed? According to Wikipedia, the Kruskal-Wallis test statistic H is approximately Chi Square distributed under the null hypothesis. However I can’t find any resources explaining why that’s the case. 
I’m curious what the connection between the sum of squared normal distributions and H is. 
 A: I'll give a somewhat intuitive explanation of why we might expect a chi-squared approximation in large samples (including illustrating a connection to sums of squared normals) and give a couple of references along the way. 
Let us begin by starting with an ordinary F statistic for one-way ANOVA.
This F-distribution of the test statistic arises from having independent chi-squared variates divided by their respective df in the numerator and denominator when the null is true. The numerator is a variance among group means, while the denominator is the variance of observations around their group means (which shouldn't be impacted by the truth or falsity of the null).
The size of the numerator is related to the underlying $\sigma^2$, so the denominator is needed in there so that we have some way of judging when the numerator is "too big".
Imagine changing from measuring in $\textit{mm}$ to measuring in $m$; the size of the numerator in the F would reduce by a factor of a million.
 It needs to be scaled to something that is not a function of the variance (more specifically, we need a pivotal quantity).
If you replaced the data by their ranks, the numerator is not affected by the scale of the original observations (again imagine changing from measuring in $\textit{mm}$ to measuring in $m$ -- the ranks would be unchanged). We can consider the numerator of an F-statistic applied to the ranks without regard to the variance of the original data (but we will eventually have to deal with variation in the ranks). For simplicity consider sampling from continuous populations (no ties).
Now given the null and under the assumption of normality, the numerator of an ordinary ANOVA F-statistic will have a chi-squared distribution divided by its df. Let's just consider the sum of squares term instead and look at what happens when we switch to ranks. 
That is, we're looking at something like this:
$\text{SST} =  \sum_{i=1}^g n_i[\bar{R_i}-(N+1)/2]^2$
where $g$ is the number of groups, $n_i$ the number of observations in the $i$-th group, $\bar{R}_i$ is the average rank in the $i$-th group and $N$ is the overall number of observations.
Each term in that sum is $n_i$ times the square of a deviation between a group mean rank and an overall mean rank. The overall mean rank is of course simply the average of the integer from $1$ to $N$. The group mean ranks are averages of ranks to each group (under the null these will be sampled without replacement from the full set of ranks). 
In large samples, those average ranks $\bar{R}_i$ will be approximately normally distributed. If we subtract the overall mean we would have a zero-mean, approximately normal random variable with variance related to the variance of the available ranks and inversely related to $n_i$. So let's multiply that by $\sqrt{n_i}$ and divide by $\sqrt{(N^2-1)/12}$ (the standard deviation of a uniform over the available ranks):
$Z_i = \frac{\sqrt{n_i}(\bar{R_i}-(N+1)/2)}{\sqrt{(N^2-1)/12}}$
These $Z_i$ are approximately standard normal. Now if we square them and add them up we get a "standardized" version of the above numerator $\text{SST}/V(N)$.
However, we haven't accounted for the a couple of things there, including fact that these $Z_i$ values are not independent (if you know all but one of them the last one is determined).
Accounting for this, we can write $H=\frac{N-1}{N}\sum_i Z_i^2$, where $H$ should be approximately chi-squared with $g-1$ d.f.
Now simple algebra can be used to transform 
$H = \frac{(N-1)}{N} \sum_{i=1}^g \frac{n_i[\bar{R_i}-(N+1)/2]^2}{(N^2-1)/12}$
into the more usual form
$H = \frac{12}{N(N+1)}\sum_{i=1}^g n_i \bar{R}_{i}^2 -\ 3(N+1)$
or you can just proceed from considering $\sum_i n_i[\bar{R_i}-(N+1)/2]^2$, suitably scaling that and then expanding and simplifying.
Hopefully this gives both some motivation for the form of the statistic and some indication for why we'd expect to get a chi-squared approximation by starting from the SS term in the numerator of an F applied to the ranks (and then suitably scaling it).
The second edition of Conover's Practical Nonparametric Statistics gives a little more detail for the argument (see p235, including explicit discussion of the effect of sampling without replacement), and a couple of pages later shows that $H$ is monotonic in the F-statistic applied to the ranks. That is, if you were able to produce suitable tables for $F$ applied to the ranks (they won't quite be distributed as F) and tables for $H$ (or directly produced the permutation distribution for both) they should reject in the same circumstances.
(Some other details can be gleaned from the original Kruskal and Wallis paper, but I think Conover's explanation provides a clear summary.)
Kruskal and Wallis offer several other approximations in their 1952 paper.
