Implementation of GL90 vertical viscosity parameterization #23
Conversation
|
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) |
|
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. |
|
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.
|
|
Final figure: The domain integral of the previous figure, which just underlines the previous point. |
|
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) Yes, to diagnose the conversion we need to keep a copy of the GL viscous coupling coefficients. |
|
@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? |
@NoraLoose posted 12 panels showing the hydrography. The only major differences I see are near the surface, |
|
I only see viscosity and energy line plots. Must be working on a different planet. |
|
Anyhow, the similarities are expected. |
View the PR on GitHub. #23 (comment) |
|
@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". |
|
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.) |
@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): MOM6/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 Lines 941 to 944 in d8fc2ef This solution would work for NeverWorld2 but not for more complex (non-adiabatic) MOM6 setups, correct? |
|
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). |
|
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). |
|
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. |
|
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. |
|
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. |
|
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? |
|
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.)? |
- don't let the GL90 coupling coefficient feel BBL anymore - effectively use hvel (which close to bottom has upwind biased thickness) to avoid spurious effects close to almost vanished layers and topography - GL90 coupling coefficient is now computed more computationally efficient: only at one j-index at a time, consistent with other routines - remove dependency on tv, which may have to re-added in the future for computing N^2. but not needed in SSW because here N^2 = g'/h.
|
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.
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 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 Any recommendations for what to try next? Thanks for your help! |
|
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 |
|
@adcroft, that’s a good point, thanks! There doesn’t seem to be a problem if I use (first column) a large
So the problem must be the
which corresponds to depth-independent nu_GL (see third column in first figure);
or |
|
I have looked at some more diagnostics to get to the bottom of the problem. This figure compares the vertical profile of 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
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
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 |
* switch around computation of ADp%du_dt_str and ADp%du_dt_visc_gl90 * add more comments in ADp%du_dt_visc_gl90 computation
* KE_u = u * h * du_dt_visc_gl90, where u should be the velocity before the update; this is why KE_[uv] had to be moved up within MOM_vert_friction, to ensure the original velocity can be used
- necessary otherwise model crashes if neither d[uv]_dt_visc_gl90 nor KE_visc_gl90 is asked for by user
|
Closing this PR because this code has been merged into the NOAA-GFDL MOM6 branch via NOAA-GFDL#268 and |
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:
Kv_u,Kv_v) is enhanced where stratification is weak (because nu ~ 1/N^2), and zero at the equator (because nu ~ f^2).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:
GL_work(similar to theGM_work), i.e., the energetic effect of GL90 on the KE budget. In the plot below, the blue line is the differenceKE_visc-KE_stress, and includes both theGL_workand 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). DiagnosingGL_workwould require a new diagnostic inMOM_vert_frictionalongside withd[uv]_dt_viscandd[uv]_dt_str. I think to do this, we have to compute the coupling coefficientsa_u,a_vinvertvisc_coefseparately for theKv_gl90_[uv]contributions. Would you agree?MOM_vert_friction. This subroutine is hacked together choosing the path of least resistance. For example, I am calling subroutinesfind_etaandcalc_isoneutral_slopesto 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!