Mistake in stretching matrix in multilayerqg.jl

[mpudig]

Hi,

We (Shafer Smith's group at NYU) noticed a possible typo in how the reduced gravities are defined in multilayerqg (https://github.com/FourierFlows/GeophysicalFlows.jl/blob/main/src/multilayerqg.jl#L355).

It is our understanding that in the multi-layer QG equations, one should divide by a reference density in the reduced gravities (see Vallis, 2017 p.186), rather than by the density of each layer as is currently done in multilayerqg (as well as in its documentation).

@navidcy, could you comment on whether this is in fact a typo or if there was a reason behind this choice for the reduced gravities?

Thanks!

[navidcy]

Hm... trying to remember.

Does the division with reference density implies that the Boussinesq approximation holds?

[navidcy]

The discussion #325 is probably relevant. The discussion was pointing out two things: a sign error (we fixed in #329) but also pointing out that having the denominator in the reduced gravities vary with layer brings problems in setting up an Eady problem.

We should resume the discussion. Perhaps we close this issue and continue chatting in that discussion and when we converge on what to do we should open a PR and fix it?

[mpudig]

Sorry @navidcy, I hadn't seen this extensive discussion! We originally picked up this mistake in pyqg as Matt Lobo seems to have as well - ha! Happy to chat on that discussion. I believe that yes the reference density comes from Boussinesq.

[navidcy]

No worries!

[navidcy]

hi @mpudig, don't we still need to resolve whether we divide with $\rho_1$ or $\rho_n$?

[mpudig]

Hi @navidcy, sorry for the slow uptake!

Is the conclusion from discussion #325 that we should do away with accepting $\rho$ as an input field for multilayerqg and only accept buoyancy? That way, if the stratification is constant then the reduced gravity will be constant layer-to-layer as we would expect.

multilayerqg will then not need any definition of $g$ or $\rho$ (nor a reference density) and will just require a buoyancy field input, $b$. I have had a tentative go at this change (see my PR).

If instead we do want to keep $\rho$ as an input, then yes I think we should simply divide by $\rho_1$ (the first layer value) to form the reduced gravity array (in order to adhere to the Boussinesq approximation and to avoid spurious internal PV gradients as Matt Lobo noted in the discussion).

Let me know what you think and happy to make changes as you see fit!