Equation error in GLS model coded into gls_corstep.F

[harcourt]

Hi.

I've been working on putting different turbulence closures into the COAWST version of ROMS and into GOTM, and I think I have found a longstanding problem in the published account of an aspect of the GLS turbulence closure model, and in the resulting implementation in ROMS.

The problem involves a smoothing of the transition to a limiter on the nondimensional buoyancy forcing Gh in the ROMS/Nonlinear/gls_corstep.F subroutine. In a current or recent version of this routine with svn line

!svn $Id: gls_corstep.F 830 2017-01-24 21:21:11Z arango $

the code with the indicated line numbers reads:

Code: Select all

1058             Gh=MIN(Gh,Gh-(Gh-gls_Ghcri)**2/                             &
1059      &                    (Gh+gls_Gh0-2.0_r8*gls_Ghcri))

and this is indeed consistent with the published equations for this aspect of the application of limiters to Gh with the original statement in Eq. (19) of Burchard and Peterson (1999) and its reiteration in the Burchard et al. (1999) GOTM report to the EC. This equation was altered erroneously in Eq. (33) of Warner et al (2005) but this is a typographical error (placing 'Gh - ' into the numerator rather than in front of the fraction). That typo has no direct bearing here on its longstanding implementation in gls_corstep.F, which is faithful to the preceding mathematical statements of how this smoothed limiter should be applied to Gh.

The problem instead is that the original mathematical statement of the smoothed Gh-limiter in Burchard and Peterson (1999) or in Burchard et al. (1999) appear to be incorrect, given the effect it is reported to have in the accompanying Figure 1 of Burchard and Peterson (1999) or Figure 2.3 of Burchard et al (1999). In order to obtain the behavior reported in those figures -- and which makes much more sense than the currently coded limiter -- the lines of code in gls_corstep.F need to be changed to read

Code: Select all

1058            Gh=MIN(Gh,Gh-MAX(0.0_r8,Gh-gls_Ghcri)**2/                   &
1059      &                    (Gh+gls_Gh0-2.0_r8*gls_Ghcri))

and the Equation (33) in Warner et al (2005) should read something like:

Gh = Gh_{unlimit} - MAX(0,Gh_{unlimit} - Gh_{crit})^2 / (Gh_{unlimit} + Gh0 + 2 Gh_{crit})

or

Gh = Gh_{unlimit} - (Gh_{unlimit} - Gh_{crit})^2 / (Gh_{unlimit} + Gh0 + 2 Gh_{crit}), for Gh_{unlimit} > Gh_{crit}

Because Burchard et al (1999) and Burchard and Peterson (1999) apply this limiter to a differently-defined nondimensional buoyancy forcing \alpha_N with the opposite sign, the equations (2.90 & 19, resp.) in those publications should have read something like

\alpha_N = MAX{ \tilde{\alpha}_N , \tilde{\alpha}_N - MIN(0,\tilde{\alpha}_N - \alpha_c)^2/(\tilde{\alpha}_N + \alpha_{min} - 2 \alpha_c)

or as is but appended:

\alpha_N = MAX{ \tilde{\alpha}_N , \tilde{\alpha}_N - (\tilde{\alpha}_N - \alpha_c)^2/(\tilde{\alpha}_N + \alpha_{min} - 2 \alpha_c), for \tilde{\alpha}_N < \alpha_c

.

That this is indeed the case can be verified by plotting the different resulting behaviors of Gh or \alpha_N and comparing it to those figures.

The effect of this model error is that when unstable buoyancy forcing is very strong, vertical diffusivity gets unintentionally shut down over the range of Gh approaching the critical value Ghcri, before rising back up to unlimited Gh at Ghcri and transitioning smoothly beyond that as intended, to a maximum value. The smoother is working right for Gh>Ghcri, but puts a notch in Gh down to zero just below Gh=Ghcri, and that notch is wider for some gls closure options than others.

In evaluating the impact of this, bear in mind this is also the regime where local single-point closures are expected to perform poorly due to their typical lack of a nonlocal buoyancy flux. Having the vertical diffusivity turn off when even the correctly limited diffusivity underpredicts buoyancy flux exacerbates this poor performance. Dysfunctionality due to the limiter error has probably been attributed to this fundamental model feature, so that may explain why the erroneous limiter has stuck around so long.

Other implementations such as current/recent versions of GOTM appear to have unceremoniously abandoned this limiter smoothing in favor of a simpler hard upper limit on Gh -- I could not find corresponding code in GOTM. A current or recent implementation in NEMO does show vestiges of an earlier presence in the definition of an unused parameter rghcri in their NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 corresponding to Ghcri.

If anyone is aware of other models where this issue of the incorrectly determined smooth Gh limiter has propagated, please add a comment or direct attention here.

Ramsey Harcourt
University of Washington, APL
harcourt@uw.edu

References

Burchard, H., Bolding, K., Villarreal, M.R., 1999. GOTM––a general ocean turbulence model. Theory, applications
and test cases. Technical Report EUR 18745 EN, European Commission.

Burchard, H., Petersen, O., 1999. Models of turbulence in the marine environment. A comparative study of twoequation
turbulence models. J. Marine Syst. 21, 29–53.

Warner, J. C., C. R. Sherwood, H. G. Arango and R. P. Signell, 2005: Performance of four turbulence closure models implemented using a generic length scale method, Ocean Modelling, 8, 81-113.

[harcourt]

Here is a Matlab program (with apologies to open-source folk) to demonstrate the problem with the published & ROMS-coded equation for the smoothed Gh limiter, based on 3 different closure options that invoke it in ROMS within gls_corstep.F . High resolution of Gh is needed to see the problem in discretized plot, as the error is only significant over a small interval below Ghcri, but the behavior looks like it could trip up closure model performance.

Code: Select all

closure(1).name = 'CANUTO A';
closure(1).ROMSComment = 'gls_Gh0 = 0.0329_r8 ! 0.0329 GOTM, 0.0673 Burchard';
closure(1).gls_Gh0 = 0.0329;
closure(1).gls_Ghcri = 0.03;

closure(2).name = 'CANUTO B';
closure(2).ROMScodeComment = 'gls_Gh0 = 0.0444_r8 ! 0.044 GOTM, 0.0673 Burchard';
closure(2).gls_Gh0 = 0.0444;
closure(2).gls_Ghcri = 0.0414;

closure(3).name='Default MY25/KC94/G88';
closure(3).ROMScodeComment = 'None, this is just the default if not Canuto A or B';
closure(3).gls_Gh0 = 0.028;
closure(3).gls_Ghcri = 0.02;

gls_Ghmin = -0.28;

for ic=1:length(closure),
gls_Gh0=closure(ic).gls_Gh0;
gls_Ghcri=closure(ic).gls_Ghcri;

% Gh_ul is unlimited Gh
% Gh_li_eqn is limited Gh per published eqn. and implemented in ROMS code
% Gh_li_fig is limited Gh consistent with figures accompanying eqns.

Gh_ul=-.3:0.0001:0.1;

Gh_li_eqn=min(gls_Gh0,Gh_ul); 
Gh_li_eqn=min(Gh_li_eqn,Gh_li_eqn-(Gh_li_eqn-gls_Ghcri).^2./(Gh_li_eqn+gls_Gh0-2*gls_Ghcri));
Gh_li_eqn=max(Gh_li_eqn,gls_Ghmin);

Gh_li_fig=min(gls_Gh0,Gh_ul); 
Gh_li_fig=min(Gh_li_fig,Gh_li_fig-max(0,Gh_li_fig-gls_Ghcri).^2./(Gh_li_fig+gls_Gh0-2*gls_Ghcri));
Gh_li_fig=max(Gh_li_fig,gls_Ghmin);

fh(ic)=figure;
plot(Gh_ul,Gh_li_fig,'g',Gh_ul,Gh_li_eqn,'r--')

title(char(closure(ic).name))
xlabel('Unlimited Gh')
ylabel('Limited Gh')
legend('Limiter matching published fig.','Limited per published eqn.','location','NorthWest'), legend boxoff
end

[arango]

Thank you for looking this carefully and for providing a Matlab script to show the difference between the formulation currently coded in ROMS and the correction to the limiter (green curve). Obviously, the coded formulation has a discontinuity spike just before it remains constant. It is clearly wrong. I don't know what the repercussions of this error in the solutions are. I don't have the time to explore it. I will update the code with your limiter suggestion.

We should start giving certificates to users like you that provides us such a clear information and solution to an algorithm error.

[jcwarner]

Hey Hernan. Before you update the code, how about you let me take a look at it and get back to you next week. I spent a good bit of time to code this in, and in a separate email chain Hans is looking at this as well. Give the people who coded it a chance to look it over first.

thanks,
john

[arango]

I haven't updated the code yet. I think that Ramsey Harcourt is correct. The shape functions used in mostly all vertical mixing parameterization should be smooth and continuous. We should stay away from sharp discontinuities that lead to eye bullet values in vertical diffusion and vertical viscosity. I get frequently such bullets in my simulations with the GLS. It is good that you and Hans are looking at it. Ramsey field of expertise is LES and turbulence, and have published good papers about it.

[harcourt]

I had some exchange with Hans Burchard on this with CC to John Warner -- ought to have cc'd you too, Hernan. I summarize that exchange here:

Burchard:

Around 2005, we had a major update in the turbulence routines coordinated by Lars, and I cannot find any pre-2005 GOTM code where we use the limiting method by Burchard & Petersen (1999). As you say, without a great ceremony we have moved to a hard limiter as

an = max (an, 0.5 an_min)

and that works fine, it seems.

I do also think that our old implementation

an = max (an, an-(an-an_c)**2 / (an + an_min - 2 * an_c) )

worked, since probably, we had an if-condition (if (an.lt.an_c) then ...) in front of it. With that included, the above equation should be identical with what you propose in your post:

an = max (an, an-min(0,an-an_c)**2 / (an + an_min - 2 * an_c) )

It is all really long ago that I dealt with things like this, but isn't this problem only relevant for quasi-equilibrium stability functions, which depend on an only? As far as I know it, most people working with GOTM use the non-equilibrium versions. Isn't that also be true for the ROMS community?

Me (Harcourt):

... the current implementation of this [smooth] limiter in ROMS should be changed, either as I suggested, to correspond accurately to the original intention, or to a hard limiter, which should be located at anLimitFact=0.5 times the value (Ghmax or anMin) of Gh or an where the equilibrium state would force Shear production Sm*as to become negative. The former change would be quick and easy, following the discussion ... on the board, while the latter may require translating the functional determination of (2x) the limit anMin on 'an' in GOTM to one on Gh in ROMS....

More info/comment:

The impact of correcting the existing code to the smooth limiter originally intended by Burchard and Petersen (1999), following proposed code change in my original post, may well be minor for the 'weak-equilibrium' closures, where Gm or as are not just dictated as functions of Gh or as in the equilibrium state, as they are for quasi-equilibrium closures of Kantha & Clayson, (1994); Galperin et al., (1988) and Harcourt (2013; 2015). I'm uncertain of the scale of problems it could be responsible for in quasi-equilibrium closures, but it is likely to be more significant because this limiter (right or wrong...) is normally 'active' for convectively unstable profiles in those closures, so fixing the code to remove the discontinuities just below Ghcri does seem like the prudent thing to do.

Updating the ROMS GLS weak-equilibrium closures to be consistent with the more recent formulations and limiters of Umlauf & Burchard (2005 CSR Review) might be a more extensive endeavor, as the newer limiters on Gh (an) and Gm (as) in current GOTM implementation (see src/turbulence/cmue_c.F90 of GOTM code) look significantly different. In addition to a different formulation for the hard limiter on Gh (an), there appear to be dependencies on higher powers of Gh (an) for the limiter on Gm (as) in weak-equilibrium closures. It's not clear to me if these limiters would then actually be different for the weak-equilibrium closures implemented in ROMS, or if those extra terms involved would end up just being zero for Canuto A,B.