What is the primary reason that someone would apply the square root transformation to their data? I always observe that doing this always increases the $R^2$. However, this is probably just due to centering the data. Any thoughts are appreciated!
In general, parametric regression / GLM assume that the relationship between the $Y$ variable and each $X$ variable is linear, that the residuals once you've fitted the model follow a normal distribution and that the size of the residuals stays about the same all the way along your fitted line(s). When your data don't conform to these assumptions, transformations can help.
It should be intuitive that if $Y$ is proportional to $X^2$ then square-rooting $Y$ linearises this relationship, leading to a model that better fits the assumptions and that explains more variance (has higher $R^2$). Square rooting $Y$ also helps when you have the problem that the size of your residuals progressively increases as your values of $X$ increase (i.e. the scatter of data points around the fitted line gets more marked as you move along it). Think of the shape of a square root function: it increases steeply at first but then saturates. So applying a square root transform inflates smaller numbers but stabilises bigger ones. So you can think of it as pushing small residuals at low $X$ values away from the fitted line and squishing large residuals at high $X$ values towards the line. (This is mental shorthand not proper maths!)
As Dmitrij and ocram say, this is just one possible transformation which will help in certain circumstances, and tools like the Box-Cox formula can help you to pick the most useful one. I would advise getting into the habit of always looking at a plots of residuals against fitted values (and also a normal probability plot or histogram of residuals) when you fit a model. You'll find you'll often end up being able to see from these what sort of transformation will help.
The square-root transformation is just a special case of Box-Cox power transformation (a nice overview by Pengfi Li, could be useful reading and is found here), with $\lambda = 0.5$ and omitting some centering.
The aim of the Box-Cox transformations is to ensure the usual assumptions for Linear Model hold. That is, $y\sim N(X\beta, \sigma^2 I_n)$.
However this a priori fixed value could be (and probably is) not optimal. In R you may consider a function from
powerTransform that helps to estimate an optimal value for Box-Cox transformations for each of the variables participated in linear regression or any data you work with (see the
example(powerTransform) for further details).