Definition of reduced gravity in MultiLayerQG model

[mjclobo]

With the current definition of reduced gravity in the multi-layer QG model it isn't possible to produce a pure Eady instability. Namely, even with constant stratification, constant vertical shear of horizontal velocity, and zero beta there will still be an interior PV gradient.

This is because the stretching matrix's interior rows won't be symmetric about the center. I.e., using the notation from the documentation, in a pure Eady case we should have $F_{j-1/2,k} = F_{j+1/2,k}$ for rows $k=2,..,n-1$, because we should have $g_{j-1/2}' = g_{j+1/2}'$.

But we don't have equal reduced gravities since their denominators will be different! This leads to an interior PV gradient.

To me this seems problematic if one was interested in, for example, the vertical structure of idealized instabilities.

Is this a common convention in QG models? If so, what are the reasons for defining reduced gravity this way? I noticed that pyqg does the same thing as GeophysicalFlows.jl. Though Vallis, 2017, uses the surface value of density in the denominator of his definition of reduced gravity for a QG model.

Thanks for your time, and thanks so much for this amazing model! I've learned so much just by digging through the code, and hope to use a version of the model for some of my PhD research :)

[navidcy]

Thanks for pointing this out. I'll try to think though and get back to you!

[navidcy]

I wanted to revisit the MultiLayerQG docs to add a derivation from the continuous QG equations to the discrete, layered version. I'll try to do that and in the process think about the denominators of g'. :)

[mjclobo]

Thanks for the reply, Navid. That sounds great. No rush! I look forward to hearing your thoughts.

I actually just found a true bug in the same function. On line 367 of multilayerqg.jl you should have "+ Fp[j]" rather than "- Fp[j]".

I'm going to create a pull request for this issue, since it is more obviously incorrect. For now I won't add the reduced gravity definition correction.

The easiest way to validate these changes is by running a 3+ layer (I'm using 8) model with no beta, constant shear, and constant stratification, and seeing if there is zero interior PV (in params.Qy).

Note that by only changing the more obvious bug (the sign error) you will still have non-zero interior PV. However, the non-zero interior PV will be of much smaller magnitude than when the sign error isn't fixed.

[glwagner]

That's useful @mjclobo . If you want to be complete, you can first open an issue and document the problem (ie link to the wrong line, and show some plots that illustrate how it leads to spuriously large interior PV).

Note that you can paste links directly to lines of code in github, and github will usefully display the code on that line (if its in the same repo), ie

GeophysicalFlows.jl/src/multilayerqg.jl

Line 367 in c0c42ce

CUDA.@allowscalar @views Qy[:, :, j] = @. Qy[:, :, j] - Fp[j] * (U[:, :, j+1] - U[:, :, j]) + Fm[j-1] * (U[:, :, j-1] - U[:, :, j])

[mjclobo]

Thanks for the tip @glwagner. I opened issue #328. Let me know if you'd like any further detail there!

[mjclobo]

For further clarification of this possible issue, I have made a couple of plots. Both plots can be reproduced using my script found in a branch of this repo that I made for this issue. Reduced gravity is defined here in multilayerqg.jl:

GeophysicalFlows.jl/src/multilayerqg.jl

Line 355 in f8fd4b4

g′ = T(g) * (ρ[2:nlayers] - ρ[1:nlayers-1]) ./ ρ[2:nlayers] # reduced gravity at each interface

In the plots "original g'" refers to the current definition used by GeophysicalFlows.jl (and pyqg; shown above) and "constant $\rho$ ref" refers to possible change in the definition of reduced gravity, as given by Vallis, 2017, i.e., $g_i' = g (\rho_{i+1} - \rho_{i})/\rho_{0}$, where $\rho_{0}$ is some reference density, generally taken to be that of the uppermost layer.

The first plot compares results comparing use of the two definitions of reduced gravity under scrutiny, when using low absolute values of density (i.e., $\mathcal{O}(1)$ ), like those given by the Phillips instability example. The left plot shows meridional gradients of background potential vorticity. The right plot shows the vertical structure of the most unstable mode, found using traditional linear stability analysis. The current definition of reduced gravity leads to a spurious interior potential vorticity gradient, which leads to a subsurface baroclinic instability. Also, using unrealistically small density values leads to small boundary PV gradients that aren't able to interact via Eady edge ways to cause an instability.

The second plot compares results using $\mathcal{O}(1000)$ values of density, as in the ocean. Firstly, we see that since $\Delta \rho / \rho = \mathcal{O}(10^{-3})$ there is no longer a difference in the two definitions of reduced gravity. Secondly, we see that using realistic values of density leads to strong enough boundary PV gradients to cause an Eady instability (cf. the "constant $\rho$ ref" vertical structure in the first plot).

One further point is that layer-wise conservation of volume is already an assumption in QG models that uses the Boussinesq approximation (pointed out to me by Stephen Griffies). From these results it is suggested that (a) the denominator in the definition of reduced gravity be changed to a constant reference density, (b) the Phillips BCI example's density values should be changed to $\mathcal{O}(1000)$, and (c) the Boussinesq assumption should be explicitly stated in any future derivation of the layered QG model.

Happy to discuss any of this further!

[mjclobo]

Sorry, I wasn't clear in my assumptions. I meant absolute density, relative to $\mathcal{O}(1)$ values of $\Delta \rho$, or just the scale $\mathcal{O} (\Delta \rho / \rho )$, like you've said.

I'm unsure of what you mean by the Phillips problem being non-dimensionalized. But what matters is that the density values that you pass in (currently $\rho$ = [4.0, 5.0]) are those used in the definition of reduced gravity, which both you and I have shown. It is here that their magnitude, relative to their difference, matters, right? Since this definition of reduced gravity is what is used to calculate the meridional PV gradient.

And yes, what you're pointing out is what I'm suggesting to change. I believe that the denominator should be a reference value, generally taken as $\rho_{1}$.

[glwagner]

I'm unsure of what you mean by the Phillips problem being non-dimensionalized.

I'm referring to the example: Phillips instability example

where f=1, g=1, and L=2pi

[glwagner]

And yes, what you're pointing out is what I'm suggesting to change. I believe that the denominator should be a reference value

I guess we are specifying something like

density = reference_density + perturbation

But with the current interface, it is not enforced that the perturbations are small relative to the reference density. As a result one can easily specify a problem incorrectly.

A different interface would explicitly require defining reference_density and the profile of perturbations.

Does that encapsulate this issue? There is a problem when the reference is not large compared to perturbations? I'm not sure we want to define a reference implicitly. It seems better to require it to be explicitly specified if we are going to make that distinction.

[mjclobo]

Yes, that captures the problem well, though two things are being talked about that are separate, I think.

(1) Technically, any constant reference density will fix the problem that I illustrated above. When you have constant shear and stratification, any constant reference density in the denominator of the reduced gravity definition will give you zero interior PV gradient, regardless of the magnitude of that reference value. This is because $\Delta \rho$ is constant for interior layers, so we just need the density in the denominator to be constant, in order for PV to be constant (i.e., no interior gradient).

(2) Being consistent with the Boussinesq approximation, the densities should satisfy $\mathcal{O}(\Delta \rho / \rho_0 ) \ll 1$.

Satisfying (2) makes (1) unimportant, using an order of magnitude argument (as shown by my second set of plots, above).

So there are two ``degrees of freedom'' with this problem. But what you suggest solves both of these issues. Thanks for your time :)

[glwagner]

I think a while ago we discussed specifying something like the buoyancy of each layer, and dispensing with gravitational acceleration (and a reference density) altogether. Then the reduced gravity is just

g′ = b[2] - b[1]

Not sure if that helps. I thought we had good reasons not to use buoyancy though, but I can't remember them...

Thanks for your time :)

no problem, of course

[Empyreal092]

Hi all,

Ryan Du here, another PhD student of Shafer's at NYU. Matt (@mpudig) and I realized this issue as well. And indeed as @mjclobo points out, pyqg does the same thing. We opened an issue about this in pyqg a while back (pyqg/issues/325).

I also think the choice is problematic. With the simplest constant stratification case (which includes Eady), the stretching matrix is not the discretization of dz(dz()), which should be symmetric away from boundary conditions.

Vallis 17 has derivations of multi-layer QG on p.185 with the definition of $g'$ in the text below (5.88b). He uses a constant reference density. The two-level derivation from Boussinesq is at p.195. For this issue, it is simpler to think about the Boussinesq derivation, where the reference density is in the definition of $b=-g\delta \rho/\rho_0$, where $\rho_0$ is a constant reference density. @glwagner suggestion of taking $b$ as input and defining $g′ = b[2] - b[1]$ should work.

[navidcy]

Sounds good! Let's open a PR to fix this? If you wanna start on this I can help out where needed.

[navidcy]

Perhaps we can include another example with MultiLayerQG that showcases the Eady problem?

[glwagner]

Beautiful! Buoyancy is the Boussinesq Way.

[navidcy]

PR #343 and v0.16 deal with this.