GLwork is not sign-definite

[NoraLoose]

Issue

I find locations along the continental shelf where the GLwork diagnostic (coded up as KE_visc_gl90) is positive (red shading). This is contrary to the expectation that the GLwork should be strictly negative, both in the continuous and discretized world: it describes the energetics of a vertical viscosity. The figure shows the 5-day averaged diagnostic for a 1/2 degree simulation. The problem gets worse with coarser resolution.

Most of these positive values we are seeing are O(10^-13) and smaller, but the largest value of GLwork is O(10^-5). A distribution is shown here:

How to fix it?

The GLwork diagnostic KE_visc_gl90 (together with the GL90 parameterization) are currently part of this branch. The most complex part of computing KE_visc_gl90 is to track the vertical velocity tendency d[uv]_dt_visc_gl90 that arises from the GL90 viscosity through the tridiagonal solver. I followed what is done for d[uv]_dt_visc. Here is the meat of it:

https://github.com/NoraLoose/MOM6/blob/10b4ecd715ba8ec7d5e4dbf4a68a3efc1d7ce513/src/parameterizations/vertical/MOM_vert_friction.F90#L480-L504

@adcroft you mentioned that the problem of indefiniteness could have to do with the barotropic solver. Do you have an idea for what I could try to fix the issue?

[NoraLoose]

I bet @Hallberg-NOAA has some ideas too!

[NoraLoose]

@adcroft and @Hallberg-NOAA, sorry for pinging you again. I am still interested in whether you have ideas for how to fix the problem described above!

[adcroft]

My suggestion was that the line https://github.com/NoraLoose/MOM6/blob/10b4ecd715ba8ec7d5e4dbf4a68a3efc1d7ce513/src/diagnostics/MOM_diagnostics.F90#L1119 is using uh which has been adjusted by the barotropic solver. I wondering if that adjustment has a small shear component due to the BT solver knowing about the shape of the bottom boundary layer. @Hallberg-NOAA is the right person to explain this.

[Hallberg-NOAA]

@adcroft and I think that the problem has to do with the way that the barotropic solver and the vertical viscosity handle their velocities. The velocities in uh and vh are adjusted to give the time-mean barotropic transport, whereas u and v themselves are instantaneous velocities at the end of the timestep. One way to test this would be to run the model with SPLIT=False. This would require the use of a very short dynamic timestep (to resolve the external gravity waves that are usually captured by the barotropic solver), but in this case the transports and velocities are associated with the same instantaneous velocities at the end of the timestep. The corrections are subject to a viscous correction, so it should be possible to reconstruct the full energetic consequences of the GL viscosity, but first let's try to confirm that this is what is going on. For a thorough description of these aspects of the split baroclinic-barotropic timestepping scheme, see Hallberg, R., and A. Adcroft, 2009: Reconciling estimates of the free surface height in Lagrangian vertical coordinate ocean models with mode-split time stepping. Ocean Modelling, 29(1), DOI:10.1016/j.ocemod.2009.02.008.

Another possibility is that this is arising because the time-step-averaged thicknesses at velocity points as used by the continuity solver are not exactly the same as those used by the vertical viscosity code. If this is the problem, changing to an unsplit time stepping scheme could still lead to negative work terms. If this is the case, solutions would require an extensive revision to parts of the dynamic core (which we are already contemplating as an option), but they might be a good idea for a variety of other reasons.

[NoraLoose]

@adcroft @Hallberg-NOAA thanks so much for your input!

I tested your first idea and set SPLIT = False and a very short timestep. I am seeing fewer points with positive values (see right figure), but they are still there. Again, the plots show 5-day averaged GLwork diagnostics.

cc @iangrooms

[adcroft]

I still think we're on the right track. We had two possibilities for why the uh used in the diagnostics differs from the u resulting from the vertical solver; time-averaging and vertical structure. We haven't eliminated either but whichever it is, I think the safest way to make this diagnostic sign definite is to move the KE-rate calculation (currently at MOM_diagnostics.F90#L1115-L1132 ) to where the ADp%du_dt_visc_gl90 is calculated at MOM_vert_friction.F90#L483-L503 and use the actual u that results from the vertical solver at MOM_vert_friction.F90#L477. We might find the full (non-GL) viscosity term has the same problem if we looked at it carefully. Also, moving the diagnostic might mean the overall budget should have a term associated with the barotropic adjustment.

[NoraLoose]

@adcroft thanks for this suggestion! I moved the KE_visc_gl90 diagnostic into the MOM_vert_friction routine, to enable usage of u and h_u that are used in the vertical solver.

This commit reflects the changes I made to arrive at the new KE_visc_gl90 diagnostic.

Unfortunately, KE_visc_gl90 is still sign-indefinite:

I don't really understand why...

[adcroft]

😕 How big are the negatives now?

Double checking we agree on the math. If the vertical viscosity update is summarized as
$$ h u - dt \delta_k \left( \frac{ \nu }{ \bar{h} } \delta_k u \right) = h u^* $$
then
$$ h u^2 - h u u^* = dt u \delta_k \left( \frac{ \nu }{ \overline{h} } \delta_k \right) = h u^* $$
$$ = \delta_k \left( \overline{u} \frac{ dt \nu }{ \bar{h} } \delta_k \right) - \overline{ \frac{ dt \nu }{ \overline{h} } \left( \delta_k \right)^2 } $$
where a_cpl -> $\frac{ dt \nu }{ \bar{h} }$
then the first terms sums to zero (assuming $\nu=0$ at top and bottom), The second term is sign definite IF $\nu$ is single signed.

If we split $\nu = \nu_o + \nu_{GL}$ nothing changes so long as the $\nu_o$ and $\nu_GL$ are each single signed.

I don't see anywhere in the code the boundary values of $\nu$ are being used so believe they are implicitly zero.

Let's (sanity) check that KE_u is sign definite if we substitute du_dt_GL90 with `du_dt_str in your diagnostic.

I am unsure what the "underflow" block does to the sign. Hopefully no damage but "to find out, comment out".

I ask at the top about the magnitude because the code is making me wonder about roundoff. Both the math and code rely on the rearranging of differences and in both cases the residual of the differences is potentially inaccurate when difference large coherent quantities. The term we are indirectly obtaining and that you are plotting is a_cpl_GL90 * ( v(i,J,K-1) - (v(i,J,K) )**2 summed over k=2:nz . This is unequivocally sign-definite (subject to the sign of a_cpl_G90). We are currently diagnosing the left hand side of my first equation, which is a difference and so subject to round off. Calculating the full right hand side also has a difference, but only the first term could give us a roundoff problem.

Aside: I suggest moving block L660-L666 to above block L634-L658 so that c1 is not corrupted by being re-used. I see you were careful to make sure dv_dt_visc still has the original value of v in it. I don't think moving this block fixes the problem.

[NoraLoose]

@adcroft thanks so much for taking a deep dive into my code!

I agree with the math (I just corrected some obvious typos):

I don't see anywhere in the code the boundary values of ν are being used so believe they are implicitly zero.

Correct 👍

Dumb question: Using your notation from above, aren't we computing $h(u*)^2$ rather than $huu^*$ in the energetic calculation? As of now, the velocity update
https://github.com/NoraLoose/MOM6/blob/557ccb99643aab61e0c8c76ef23392522a28e06d/src/parameterizations/vertical/MOM_vert_friction.F90#L500-L504
happens before I compute the diagnostic du_dt_visc_gl90 in
https://github.com/NoraLoose/MOM6/blob/557ccb99643aab61e0c8c76ef23392522a28e06d/src/parameterizations/vertical/MOM_vert_friction.F90#L506-L530

Would it make sense to switch the order?

[NoraLoose]

How big are the negatives now?

The amplitude of the positives are similar to our old diagnostic, when the energetic calculation was done in MOM_diagnostics. The number of positives has slightly decreased, though.

[NoraLoose]

@adcroft I will now perform the other sanity checks that you suggested, including

• using du_dt_str rather than du_dt_visc_gl90 in the diagnostic
• commenting out the "underflow" block (Edit: did not make any difference)

[NoraLoose]

@adcroft @Hallberg-NOAA @iangrooms good news: I finally figured out how to make the KE_visc_gl90 diagnostic sign-definite. My mistake was the following: Even when I calculated KE_visc_gl90 within MOM_vert_friction, I did so at the very end and used the updated velocity u^* for the energetic calculation, i.e.,

KE_visc_gl90 = u^* * h * du_dt_visc_gl90

where I am using Alistair's notation from above.

But what we actually need is

KE_visc_gl90 = u * h * du_dt_visc_gl90,

which uses the original velocity u before the update. This commit, which moves the computation of KE_[uv] up within MOM_vert_friction did the trick (see right most figure):

Thanks everyone for your help with fixing this issue!