Why does $R^2$ grow when more predictor variables are added to a model? I do understand that $ R^2 = \frac{\text{SSR}}{\text{SST}}= 1- \frac{SSE}{SST}$, however, I don't understand what changes when more predictor variables are added and how $R^2$ is affected accordingly. Can someone clarify? 
 A: It is the result of the fitting process that takes place in the OLS regression. Each variable is regressed against all others, and what is left unexplained (residuals) is carried over. In a way, the regression process looks for explanations in the variance in the data, but it doesn't really excel at telling what is signal and what is noise.
In fact, if you were to just include variables composed of random noise, you could still see how there would be progressive overfitting of this noise in a misleading attempt at explaining the variability in the "dependent" variable. 
I did this test in here, and plotted the resulting effect on the $RSS$ as the number of non-sensical variables increased:

This is why it is advisable to use adjusted $R^2$ instead of $R^2$ to judge whether it is a good idea to include more variables in a model.
A: Let's suppose that we've got two models:
$$
Y = \beta_0 + \beta_1 X_1 + \varepsilon \tag{M1}
$$
and
$$
Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \varepsilon \tag{M2}
$$
This means that we have
$$
RSS_1 = \sum_{i=1}^n (Y_i - \hat \beta_0 - \hat \beta_1^{(1)} X_1)^2
$$
and
$$
RSS_2 = \sum_{i=1}^n (Y_i - \hat \beta_0 - \hat \beta_1^{(2)} X_1 - \hat \beta_2 X_2)^2.
$$
Model $M2$ contains model $M1$ as a special case, so there is no way that $RSS_1 < RSS_2$: we can just set $\hat \beta_2 = 0$ and $\hat \beta_1^{(1)} = \hat \beta_1^{(2)}$ in order to get $RSS_1 = RSS_2$. Much more likely is that $RSS_2 < RSS_1$ because we have an extra parameter so we can fit the data more closely.
This reveals the big problem with the unadjusted $R^2$: there is no penalty for model complexity. A more complicated model will almost always fit the data better so $R^2$ will prefer this model, even if the extra complexity is just modeling noise. That's why other methods like the adjusted $R^2$ (as mentioned in Antoni Parellada's answer) and $AIC$ are popular, since these take into account both the fit of the model to the data while also penalizing model complexity.
