Conceptual bug in Redi slope triad calculation

[vanroekel]

Currently the slope triad computations (the critical ingredient of the redi parameterization) are defined to do horizontal gradients in model index space, see here

https://github.com/E3SM-Project/E3SM/blob/master/components/mpas-ocean/src/shared/mpas_ocn_gm.F#L304-L322

Given this quantity should be the isopycnal slope in physical space, I believe we need to redefine

dTdx = (activeTracers(indexTemperature, k, cell2) &
                          - activeTracers(indexTemperature, k, cell1)) &
                         *dcEdgeInv

as

dTdx = (activeTracers(indexTemperature, k, cell2) &
                          - activeTracers(indexTemperature, k, cell1)) &
                         *dcEdgeInv + dTdz*dz/dx

where dTdz and dzdx would be edge based quantities (similar to how the PGF zmid is calculated). Through much of the global ocean this won't matter, but will matter near bathymetry and at the ice cavity transition. A similar change is needed for salinity (dSdx) I think.

[vanroekel]

pinging @mark-petersen in particular for thoughts.

[vanroekel]

also noting that the Griffies 1998 paper on this does say 'isoneutral diffusion in z-coordinate models', so I think for z* we need an alteration and for PBCs too I think

[mark-petersen]

I'm not sure about this. I called the variable dTdx for simplicity, but if it the coordinates are tilted, this is the derivative in the direction of the horizontal coordinate, so we could have called it dTds where z=s(x) is the surface following the vertical middle of the cells. I need to think about this more.

[vanroekel]

I don't quite understand what you are saying. My thinking is the slope needs to be in z-space and now it is in model coordinate space (or middle of the cell space). The diffusion should be along the isopycnal and I think you'd want the horizontal derivative to be on a fixed z not between model levels.

I guess a simple test would be to modify the parabolic bowl to have large horizontal variations in thickness with steep isopycnals and see if the diffusion remains isopycnal

[mark-petersen]

@vanroekel I thought more about this, and I agree with your code alteration.

I was tripped up by your wording previously. You said

Given this quantity should be the isopycnal slope in physical space ... dT/dx = ...

I thought you meant that x should be an along-isopycnal derivative here. But I think you meant was simply

Given this quantity should be the temperature slope in physical space ... dT/dx = ...

OK, so let's start at the beginning. Let's call $(x,z)$ the physical, orthogonal coordinate system. Since we have tilted coordinates, we will use zeta: $\zeta(x)$ for the vertical location of a layer. So now our 3D arrays, like temperature, are $T(x,\zeta(x))$ and a derivative in the physical horizontal direction $x$ is:

$\frac{d}{dx} T(x,\zeta(x)) = \frac{\partial T}{\partial x} + \frac{\partial T}{\partial\zeta} \frac{d\zeta}{dx}$
The second term corresponds with your addition + dTdz*dz/dx above, so this makes sense. We should make a note that in the code dTdz and dz/dx are in array space, so that this z is really the zeta in the expansion above.

Looking back at the Griffies et al. 1998 paper, our dT/dx and dT/dz terms are used here:

where in the compact notation m can either be x or z. The $\theta$ is temperature and $s$ is salinity. This is then used to compute the isopycnal slope here:

The Griffies paper only considers a z-level model, as those were the models they were writing for at the time (MOM and POP).

Note that the term dTdz*dz/dx will be small due to the aspect ratio of the grid cells. Even in strongly tilted layers near bathymetry and ice shelves, dz/dx might be at most 100m/10km = 0.01. But dT/dz could be much larger than dT/dx in a stratified ocean, so perhaps the correction could be important in locations with strong tilting and stratification.

[vanroekel]

Thanks @mark-petersen I see the confusion now. Sorry for that. I think we are fully on the same page. In thinking about this, the one place I think it might matter most is near ice shelf cavities. Let me put the proposed changes into master and do a G-case with and without the changes as a first step.