Partially answered in comments:
There are two ways to think about it. One is to think of it as an interaction with itself. When we use interactions we generally do include the main effects and the interaction term: e.g. on ucla site. Another is to think of it as a polynomial. There is no reason why the relationship should be linear. If a polynomial fits better then why not use it? The Data Generating Process might be non-linear.
– Toby
BTW in this case think of your lungs as two cylinders that are a function of someone's size. The volume of a cylinder is $πr^2h$. People who are taller probably are also larger in other dimensions, they tend to be right?, so the $r$ will be larger. And since the cylinder is a nonlinear function of $r$ it will be a nonlinear function of size. It seems to make sense to me to include height squared. Write it in logs and you have $\log(πr^2 h) = \log(\pi) + \log(r^2) + \log(h)$. If $h$ and $r$ are both a function of size for which height is a reasonable proxy then this is your model.
– Toby
( This argument seems to lead to log transforms of both response and height )
Using height and height^2 is fitting a quadratic in height, other predictors aside. It could be just a way of catching some curvature in the data. There are even cases, e.g. trajectory of a projectile, in which a quadratic is exactly the shape you expect.
– Nick Cox
That's very helpful. I can't see a biological reason to have height#height^2 as an interaction term? I do, however, see why you'd square height in this specific example. My take on his question is, more generally when we transform non-normal covariates, whether we should include both the transformed and original variables in the model?
– bobmcpop
Toby's argument in comments can be extended this way: often $r$ wil be only approximately proportional to $h$, $r$ wil grow slower (and slower and slower) as $h$ grows, so maybe $r\propto h^\alpha$ for some $0<\alpha<1$ is a better model. That leads to $y (\text{volume})\propto h^{1+2\alpha}$ and taking logarithms on both sides leads to a model in only $\log h$, where we in addition have an theoretical expectation for the value of slope: $1+2\alpha$, so we can see from the fit if this geometrical speculations holds for this case. And no, in this model we do not need $h$ untransformed as an additional term.