Implementation of GL90 vertical viscosity parameterization

[NoraLoose]

The parameterization "GL90''

This PR implements a vertical viscosity parameterization a la Greatbatch and Lamb (1990), Ferreira & Marshall (2006) and Zhao & Vallis (2008), hereafter referred to as the GL90 vertical viscosity parameterization. This vertical viscosity scheme redistributes momentum in the vertical, and is the equivalent of the Gent & McWilliams (1990) parameterization, but in a TWA (thickness-weighted averaged) set of equations. The vertical viscosity coefficient nu is computed from kappa_GM via thermal wind balance, and the following relation:

nu = kappa_GM * f^2 / N^2,

where f is the Coriolis parameter, and N^2 is the buoyancy frequency. kappa_GM can vary horizontally and vertically, but for now the implementation assumes horizontally and vertically constant kappa_GM. The vertical viscosity del_z ( nu del_z u) is applied to the momentum equation with stress-free boundary conditions at the top and bottom.

Using GL90 with NeverWorld2: first results

I have run NeverWorld2 0.5 degree with GL90 switched on (and GM90 switched off). This is the only switch that has to be done in order to run adiabatic NeverWorld2 in the TWA equations. The implementation is working as expected:

• Figure 1 and 2 below show that the viscosity (diagnosed online via Kv_u, Kv_v) is enhanced where stratification is weak (because nu ~ 1/N^2), and zero at the equator (because nu ~ f^2).
• The kinetic energy budget closes, and behaves as expected (see Figure 3 and 4 below). That is, when substituting GM90 by GL90, the main KE sink term is not conversion to PE anymore; instead, the main sink is the GL90 work which acts directly on the KE budget (consistent with the Bleck energy cycle).

Next steps

I want to work on a few more things before we can consider merging this. I am hoping I can get input from @Hallberg-NOAA on the following questions / roadmap:

• [ ] KE budget diagnostics: I would like to implement a diagnostic for the GL_work (similar to the GM_work), i.e., the energetic effect of GL90 on the KE budget. In the plot below, the blue line is the difference KE_visc - KE_stress, and includes both the GL_work and the energetic effects of bottom drag and vertical viscosity not associated with GL90 (e.g., through background viscosity or enhanced viscosity in the mixed layer). Diagnosing GL_work would require a new diagnostic in MOM_vert_friction alongside with d[uv]_dt_visc and d[uv]_dt_str. I think to do this, we have to compute the coupling coefficients a_u, a_v in vertvisc_coef separately for the Kv_gl90_[uv] contributions. Would you agree?
• [ ] Code efficiency & structure: I am computing the GL90 viscosity as part of a new subroutine in MOM_vert_friction. This subroutine is hacked together choosing the path of least resistance. For example, I am calling subroutines find_eta and calc_isoneutral_slopes to get at N^2. I'm sure this can be done more efficiently. Ultimately, we probably want to move the whole GL90 routine into its own module, because it will grow as we allow kappa_GM to vary vertically and horizontally. Thoughts?

@Hallberg-NOAA, or other people interested, would you have time for a meeting to discuss the code questions above? Thanks!

[NoraLoose]

Total vertical viscosity along 3 zonal sections, diagnosed online (1) using GM with kappa_GM = 1000 m^2/s (first row) and (2) using GL with kappa_GM = 1000 m^2 (second row). When GL is switched on, the vertical viscosity is enhanced by a few orders of magnitudes, especially where stratification is weak (in the deep ocean)

[NoraLoose]

As previous figure, but now vertical viscosity is shown along 3 meridional sections. These figures emphasize how viscosity goes to zero as the equator is approached.

[NoraLoose]

KE budget using (a) GM with kappa_GM = 1000m^2/s, (b) GM with kappa_GM = 2000m^2/s, (c) GL with kappa_GL = 1000m^2/s, (b) GL with kappa_GM = 2000m^2/s.

• In the GM framework (a, b), KE gets sent to PE (red line), from where it can be taken out by GM;
• In the GL framework (c, d), KE gets taken out by vertical viscosity (blue line); this "GL_work" acts directly on the KE reservoir and parameterizes the effect of baroclinic instability in the TWA framework. Conversion to PE (red line) integrates to ~zero over the domain.

[NoraLoose]

Final figure: The domain integral of the previous figure, which just underlines the previous point.

[NoraLoose]

cc @MFJansen @sdbachman @gustavo-marques @StephenGriffies @adcroft

[adcroft]

It is remarkable how similar the solutions are between the two schemes.

I agree that making a new module makes sense (eventually). In the meantime, the efficiency issue could be addressed if we modularized (broke apart) calc_isoneutral_slopes() so that you can access N2_u and N2_v more readily.

Yes, to diagnose the conversion we need to keep a copy of the GL viscous coupling coefficients.

[StephenGriffies]

@adcroft "It is remarkable how similar the solutions are between the two schemes." I did not see any figures with fields to show the solutions are similar. The energetics quite distinct. Did I miss something?

[adcroft]

I did not see any figures with fields to show the solutions are similar.

@NoraLoose posted 12 panels showing the hydrography. The only major differences I see are near the surface,

[StephenGriffies]

I only see viscosity and energy line plots. Must be working on a different planet.

[StephenGriffies]

Anyhow, the similarities are expected.

[adcroft]

I only see viscosity and energy line plots. Must be working on a different planet.

View the PR on GitHub. https://github.com/ocean-eddy-cpt/MOM6/pull/23#issuecomment-1077038807

[NoraLoose]

@adcroft and @StephenGriffies: I think you are looking at the same plots. The plots are indeed showing viscosity, but Alistair may be referring to the interface height positions (shown by the black lines) as "hydrography".

[NoraLoose]

And yes, so far I haven't seen big differences in the solutions between the two schemes, except for the energetics. I have been only looking at 500-day averages and the 0.5 degree simulation, though. (That's what the above plots are for.)

[NoraLoose]

In the meantime, the efficiency issue could be addressed if we modularized (broke apart) calc_isoneutral_slopes() so that you can access N2_u and N2_v more readily.

@adcroft I agree, this would be a good thing to do. I searched around in the MOM6 code to see how computation of N2 is handled in other places. I found a very simple solution (4 lines of code instead of calling 2 subroutines):

https://github.com/ocean-eddy-cpt/MOM6/blob/d8fc2ef6a7309e239d3edcf217b3d524594668fa/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90#L941-L944

This solution would work for NeverWorld2 but not for more complex (non-adiabatic) MOM6 setups, correct?

[adcroft]

Correct: that fn does not give N^2 if density is variable along layers or different from that implied by g' . It also has a "fix" to handle vanished layers that might not be appropriate for this application. In your case you want 1/N^2 so it might be better to think of it as h/g' in which case vanished layers are not the problem but instead a vanishing density difference (g') leads to large numbers. We should get Bob in a call to discuss options (as you asked higher up).

[NoraLoose]

Good point about the peculiarities of computing 1/N^2 (rather than N^2). This implies that we should probably write a new subroutine for this particular application (rather than calling subroutines that are specialized to compute N^2, not 1/N^2).

[iangrooms]

The Ferrari/Griffies/Nurser/Vallis approach could potentially be used here both to set the stress to 0 at the top and bottom boundaries, and also to deal with 'interpolating' through regions of low or negative N^2. Specifically, the vertical stress in this parameterization is

tau = (kappa f^2/N^2)du/dz

(possibly up to a minus sign - I can't recall the convention for the definition of stress). If you instead define tau to be the solution of

c^2 d^2tau/dz^2 - N^2 tau = -kappa f^2du/dz

and truncate negative values of N^2, then you get a stress that has the desired boundary conditions and also handles small or negative N^2 gracefully. It also allows re-purposing the existing FGNV10 code from GM to GL. I think it would require re-tooling the implicit vertical viscosity solver, but it might end up just requiring 2 tridiagonal solves instead of 1.

[MFJansen]

I think it would be good to take a step back to the SSW theory to think about what exactly this "1/N^2" term here should be, rather than taking the z-coordinate expression and then trying to approximate it in SSW.

I also agree with Ian about FGNV being potentially useful here - in fact this issue of infinite viscosity for vanishing stratification is basically the same as the issue of infinite streamfunction for vanishing stratification in GM, so all the things that people have been doing to deal with that should apply.

[adcroft]

FGNV is an interesting idea but the stress is already being solved implicitly-in-time leading to an elliptic equation. This would be a second equation for stress but I think the issue is we need to calculate a bounded viscosity.

[MFJansen]

Following up on my point above, if I’m getting this right, I think we ultimately want the form stress at each interface to be given by fK_GM\grad\eta. Using Margules, this is f^2/g’ K_GM\Delta u_g, with \Delta u_g the (geostrophic) velocity contrast at that interface (and to get the whole thing to look like a regular vertical viscosity we need to approximate u_g as u). So following this logic, 1/N^2 basically just gets dz/g' where the dz should just cancel the dz that appears in the finite difference to get the velocity shear, no?

[MFJansen]

I'm also going to repeat a suggestion I made in the chat during the call: Would it be worth trying to just implement the form stress as fK_GM\grad\eta, which in fact can be written as f \Psi_GM, with Psi_GM computed as before (thus allowing us to reuse all that code that already exists - including tapering schemes etc.)?

[NoraLoose]

@MFJansen, your derivation above matches exactly the one we did in Section 6.2 in the overleaf CPT: TWA & Stacked Shallow Water. And yes, in SSW stratification 1/N^2 is exactly equal to h/g'. (This is also what calc_isoneutral_slopes, the subroutine I am calling right now, computes if you are in SSW mode.)

[iangrooms]

I like @MFJansen's idea to try the version that is more directly connected to GM and does not require a geostrophic assumption. That said, I'm worried about numerical stability. If the vertical viscosity version is a good approximation of the interface height gradient version and we treat it explicitly in time, it might lead to instability.

If you only want to deal with viscosity instead of stress, you can do FGNV for the viscosity. Normally nu = kappa f^2/N^2, so instead let

c^2 d^2nu/dz^2 - N^2 nu = -kappa f^2

If you want nu to go to 0 on the boundaries then it's exactly the same idea as for isopycnal slopes. Not sure if we want nu=0 on the boundary or not. In my comments about FGNV I'm using expressions with N^2 as shorthand, but I agree with @MFJansen that we should ultimately use the correct expression for whatever vertical coordinate and discretization we're using.

[MFJansen]

And yes, in SSW stratification 1/N^2 is exactly equal to h/g'. (This is also what calc_isoneutral_slopes, the subroutine I am calling right now, computes if you are in SSW mode.)

I guess the question is what exactly is "h" here. I was arguing that it would be perhaps most consistent with the SSW equation to effectively define it such that h (du/dz)= \Delta u, where (du/dz) is the discrete version of the vertical shear as it enters in the model's strain tensor. So h is effectively the "Delta z" used in the model's strain computation (although if that's higher order then it gets weird... ). Maybe that is already how it's done? Sorry If I'm stating the obvious....

[MFJansen]

Ian raises an interesting question about the boundary condition. The bottom boundary condition for the stress is really the SGS bottom form stress. Setting that to zero is consistent with what current GM implementations do, so it's a good starting point, but there's an opportunity here to very naturally include a parameterization of sgs form drag. (All I'm saying is that there's a lot more fun to be had with this... :)

[Hallberg-NOAA]

I am joining this conversation late (too much scheduling all day long), so my comments may be a bit scattered across the last few hours of conversations.

First, the vertical viscosity is implemented an implicit tridiagonal solver for the velocity. It is incredibly robust because it does not actually need to calculate any differences (so roundoff is no issue), and it works very well even in the limit where the effective coupling between layers (measured as a length scale calculated from viscosity*timestep/vertical grid spacing) is very large compared with the layer thicknesses. Whatever discretization we choose these are not properties we should sacrifice.

We should never calculate or use N^-2 if we can help it, if for no other reason than that it can be infinite in easily physically realizable situations. For example, in the energetically-constrained vertical diffusivity, the natural inclination would be to use the Osborn relationship (\kappa = \Gamma \epsilon / N^2) to specify the diffusivity, but we are able to avoid this by rewriting the energetic constraint as the vertical integral of the potential energy change integrated over a timestep as a function of the local diffusivity. (My Ocean Sciences talk goes over this carefully - see the second talk under https://osm2022.secure-platform.com/a/solicitations/3/sessiongallery/schedule/items/181 if you still have your Ocean Sciences credentials.) This is one of the big reasons why our ePBL boundary layer mixing scheme works so well.

I suspect that we could do something similar for GL90. I am intrigued by the idea of calculating the GL90 viscosity with an elliptic solver. A zero boundary condition might work for the GL90 viscosity, but for the total viscosity it is absolutely the wrong thing to do. We use the implicit vertical viscosity with a stress boundary condition to get the wind stresses into the ocean and as a part of the implementation of a quadratic bottom drag law via a no-slip boundary condition and an appropriate viscosity profile. Both of these require nonzero (if not necessarily large) vertical viscosities at the ocean top and bottom, and they are critical for handling the pressure gradient calculations in massless, steeply sloping boundary layers over topography.

It should be possible to introduce explicit stresses (as in one of Malte's suggestions to base the stresses on the thermal wind shears, rather than the actual shears), but we would have to be careful to design the formulation to avoid excessively large accelerations, and to include enough implicit vertical viscosity to dominate the velocities when the layers get thin or close to the boundaries. Bill Large has been working on implementing explicit near-surface wind stresses that are at an angle to the model's resolved gradients, and this can work as long as there are also implicit down gradient stresses. An elliptic solver for the stresses might be a good approach for managing the accelerations.

[NoraLoose]

@Hallberg-NOAA: thanks for the pointer to your OSM talk, I really enjoyed it. Question: Are ocean models in SSW configuration also "Murphy's Law machines"?

In SSW, we have 1/N^2 = h/g', so this expression won't grow to infinity (since g' is fixed). There is still the problem that we compute N^2 at time step n, even though we would really want N^2 at time step n+1 to be consistent with the implicit solver. But we could argue that the relation kappa * f^2 / N^2 is only an approximation itself (obtained via thermal wind balance), so being off by a time step for N^2 may be fine.

In short, SSW NeverWorld2 could be a nice test bed to see whether GL90 improves the flow and energetics compared to GM90 at all. In SSW we don't need to do a lot of the sophisticated things that are mentioned above (I think?!). If we see that GL90 is better than GM90 in NeverWorld2, we could take the next step and think about elliptic solvers etc., but maybe we will realize that this is not even worth the effort. What do people think? Maybe I am missing something.

[Hallberg-NOAA]

All ocean models are Murphy's Law Machines, but there are fewer things that can go wrong in SSW mode.

[codecov-commenter]

Codecov Report

Merging #23 (e1e2429) into dev/cpt (d8fc2ef) will increase coverage by 0.01%. The diff coverage is 46.09%.

@@             Coverage Diff             @@
##           dev/cpt      #23      +/-   ##
===========================================
+ Coverage    29.04%   29.06%   +0.01%     
===========================================
  Files          244      244              
  Lines        71849    71953     +104     
===========================================
+ Hits         20869    20913      +44     
- Misses       50980    51040      +60     

Impacted Files	Coverage Δ
src/core/MOM_variables.F90	28.57% <ø> (ø)
...c/parameterizations/vertical/MOM_vert_friction.F90	41.78% <34.61%> (-0.76%)	:arrow_down:
src/diagnostics/MOM_diagnostics.F90	75.85% <93.33%> (+0.31%)	:arrow_up:
src/core/MOM_dynamics_split_RK2.F90	61.80% <100.00%> (ø)
src/core/MOM_dynamics_unsplit.F90	64.42% <100.00%> (ø)
src/core/MOM_dynamics_unsplit_RK2.F90	63.13% <100.00%> (ø)
src/framework/MOM_restart.F90	25.68% <0.00%> (-0.27%)	:arrow_down:
config_src/drivers/solo_driver/MOM_driver.F90	67.88% <0.00%> (-0.14%)	:arrow_down:
src/ocean_data_assim/MOM_oda_driver.F90	0.00% <0.00%> (ø)
src/user/baroclinic_zone_initialization.F90	0.00% <0.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update d8fc2ef...e1e2429. Read the comment docs.

[NoraLoose]

As implemented now, the GL90 vertical viscosity feels the effect of the bottom boundary layer because in MOM_vert_friction

1. visc_gl90_[uv] is first added to the total vertical viscosity Kv_tot
2. and then Kv_tot is modified if it is within the influence of the BBL,

see this code snippet, which is part of this PR.

The fact that the GL90 vertical viscosity feels the BBL can be seen particularly clearly when looking at the blue line in the following figure.

The figure shows the vertical profile of the diagnosed GL90 vertical viscosity at the location (20E, 50S) for a suite of experiments. The blue line is from an experiment where I specify a depth-independent vertical GL90 viscosity via

with alpha_GL = 3.e8 m^2 s and f the Coriolis parameter, via the following (newly implemented) code option:

https://github.com/NoraLoose/MOM6/blob/df52a0d980ed0aed0e2810104ec9af20303a6193/src/parameterizations/vertical/MOM_vert_friction.F90#L263-L264

We expect a vertically constant blue line, but because of the BBL the blue line has a sharp peak toward the bottom.

Question: How should the GL90 vertical viscosity be handled? Should we let it feel or not let it feel the BBL?

cc @Hallberg-NOAA @iangrooms @MFJansen

[NoraLoose]

Hi @Hallberg-NOAA, @adcroft, @MFJansen, @iangrooms, @sdbachman, @gustavo-marques et al.,

I have some issues with the GL90 vertical viscosity scheme in grid cells that are horizontally adjacent to vanished layers. To illustrate the problem, the following figure shows two meridional sections (sketched in green in subpanel on the left), close to the continental slope. Shown is the 1 degree NW2 simulation - so the two sections are for neighboring grid cells.

First row: Section along 3E, where the lowermost layers are all vanished because they are riding over the topographic slope. Second row: Section along 4E, where the same layers have thickness O(100m). As you can see, these layers have large spurious zonal velocities (left and middle panel).

The following sketch may be helpful to understand the problem better:

The coupling coefficients a_GL that go into the implicit vertical viscosity scheme live on interfaces and velocity points. (I am setting these coupling coefficients in the subroutine find_coupling_coef_gl90, which is newly added as part of this PR.) In SSW mode, they are (on paper) given by

a_GL(I,K) = kappa f^2 / g'(K)

But we probably want these coupling coefficients to not be active when we are adjacent to vanished layers. So I have tried to multiply them by h / (h + epsilon) (green term in sketch), where h is an upwind-biased thickness estimate if you are close to topography and next to vanished layers. So we should have h~0 as well as (green term) ~ 0 for the case illustrated above. However, this did not eliminate the problem.

Any recommendations for what to try next? Thanks for your help!

[adcroft]

We need to figure out if this is coming from the vertical grid used for the viscosity or the viscosity itself. OOI, what happens if you set KV to a large constant value of order the GL value? Is HARMONIC_VISC true of false (which affects the other viscosities)? Perhaps time to meet about this...

[NoraLoose]

@adcroft, that’s a good point, thanks! HARMONIC_VISC = False and KV = 1E-4 m2 s-1 for the simulations shown above.

There doesn’t seem to be a problem if I use (first column) a large KV = 1 m2 s-1, HARMONIC_VISC = False or (second column) GM instead of GL, with KV = 1E-4 m2 s-1, HARMONIC_VISC = False. In this figure, the colorbar range is [-0.01, 0.01] (compared to [-0.1, 0.1] above) to make sure we are not missing any spurious velocities. But it seems there aren't any!

So the problem must be the a_cpl_GL90 coefficient itself. Here are some hints: The problem

• does not arise if I set

        a_cpl_gl90(I,K) = 2.0 * f2 * CS%alpha_gl90 / (hvel(I,k) + hvel(I,k-1) + h_neglect)

  which corresponds to depth-independent nu_GL (see third column in first figure);

• does arise if I set

        a_cpl_gl90(I,K) = f2 * CS%kappa_gl90 / GV%g_prime(K)

  or

  
        a_cpl_gl90(I,K) = f2 * CS%kappa_gl90 / GV%g_prime(K) * (hvel(I,k) + hvel(I,k-1)) / (hvel(I,k) + hvel(I,k-1) + h_neglect)

[NoraLoose]

I have looked at some more diagnostics to get to the bottom of the problem.

This figure compares the vertical profile of Hu_visc for GM (green) vs. GL (pink) in a 1 degree simulation. Hu_visc is the thickness at u-points that is used in the vertical viscosity scheme (essentially hvel in the above schematic). The 4 columns show the vertical profiles as ones moves grid-point-wise across the continental shelf: 2E, 3E, 4E, 5E.

The second row is a copy of the first row, but each panel zooms into the leftmost part of the panel above (the first 10m). What is striking is that GL allows layers of Hu_visc = O(1m) (circled in black), while with GM the vertical transition from vanished to non-vanished layers seems to be much sharper: from Hu_visc = 0m to Hu_visc = O(100m) for the layer above.

For reference, here are the coupling coefficients used in the vertical viscosity scheme, diagnosed online for the same time step as the figure above. Total coupling coefficient au_visc (first row) and the GL90 coupling coefficient au_gl90_visc (second row). Note that the scale of the x-axis varies by orders of magnitude across subpanels.

The GL90 coupling coefficient does not do anything crazy, but it is always > 0. So maybe it somehow interferes with the peculiar vertical shape of Hu_visc?

cc @adcroft @Hallberg-NOAA

[NoraLoose]

Closing this PR because this code has been merged into the NOAA-GFDL MOM6 branch via https://github.com/NOAA-GFDL/MOM6/pull/268 and https://github.com/NOAA-GFDL/MOM6/pull/293.