excess mixing in k-epsilon GLS mixing; bug or misconfigured?

[jpringle]

Dear all --

I think I have stumbled across evidence of a bug in the GLS vertical mixing scheme in ROMS, at least for the k-epsilon scheme. It remains possible that I am just doing something stupid.

I find that under the k-epsilon scheme vertical stratification diffuses away much more quickly than under MY25 or LMD schemes, and faster than I would anticipate. I have also compared ROMS k-epsilon to a k-epsilon solution in GOTM, and also find that ROMS mixes away density gradients more quickly than GOTM with the same mixing scheme.

To test this, I have a run with linear equation of state and a linear innitial temperature stratification equivalent to a rho_z=-0.01 kg/m^4 in a periodic channel with no horizontal gradients in density or wind forcing; the only thing that should happen is diffusion of stratification.

I am eager to know if anyone knows of any way I could have misconfigured k-epsilon to get this excess mixing? And if not, where would one start to debug this code? Who is the expert on the GLS schemes in ROMS? This is a common scheme in coastal applications, and it would be nice to know that it is correctly implemented.

I have attached all files needed to replicate the run below; however, briefly, my three mixing schemes are activated with the following options:

Code: Select all

#define GLS_MIXING
#if defined GLS_MIXING 
# undef KANTHA_CLAYSON
# define CANUTO_A
# define N2S2_HORAVG
# define RI_SPLINES
#endif

#undef MY25_MIXING
#if defined MY25_MIXING 
# define KANTHA_CLAYSON
# define N2S2_HORAVG
# define RI_SPLINES
#endif

#undef LMD_MIXING 
#ifdef LMD_MIXING
# define ANA_SRFLUX
# define LMD_CONVEC
# define LMD_DDMIX
# define LMD_SKPP
# define LMD_BKPP
# define LMD_NONLOCAL
# define LMD_RIMIX
#endif

and my vertical mixing parameters for the GLS model from the *.in file are

Code: Select all

! Turbulent closure parameters.

     AKK_BAK == 5.0d-6                          ! m2/s
     AKP_BAK == 5.0d-6                          ! m2/s
      TKENU2 == 0.0d0                           ! m2/s
      TKENU4 == 0.0d0                           ! m4/s

! Generic length-scale turbulence closure parameters.

       GLS_P == 3.0d0                           ! K-epsilon
       GLS_M == 1.5d0
       GLS_N == -1.0d0
    GLS_Kmin == 7.6d-6
    GLS_Pmin == 1.0d-12

    GLS_CMU0 == 0.5477d0
      GLS_C1 == 1.44d0
      GLS_C2 == 1.92d0
     GLS_C3M == -0.4d0
     GLS_C3P == 1.0d0
    GLS_SIGK == 1.0d0
    GLS_SIGP == 1.30d0

In the following pictures, you can see profiles of temperature, log10(tracer diffusivity) and horizontal velocity (which is zero) for the mixing schemes in ROMS and the temperature solution for GOTM k-epsilon. After one day of running, the mixing levels for k-epsilon in ROMS are 6.7 times greater than those for LMD, and 63 times greater than for MY25.



By day 45 you can see that the temperature stratification is almost gone in ROMS k-epsilon, but largely intact in all the other mixing schemes, including k-epsilon in GOTM.



An animation of the evolution can be found at http://oxbow.sr.unh.edu/data/bump02_compareGOTM30.mp4

As I said above, I am eager to know if anyone knows of any way I could have misconfigured k-epsilon to get this excess mixing. And if not, where would one start to debug this code? Who is the expert on the GLS schemes in ROMS?

Thanks,
Jamie Pringle

[jpringle]

I was asked offline why I expected the diffusion of temperature to be slower than ROMS's k-epsilon in the ocean -- why did I think the GOTM k-epsilon, ROMS MY25 and ROMS LMD was closer to the truth?

In this stratified ocean, with no horizontal gradients of density and no forcing, there is no flow. Finite density gradient and zero velocity shear should result in an infinite gradient Richardson number and no turbulence production. Diffusivity should be near molecular. (in real life this is changed somewhat by breaking IWs from the unresolved GM spectra and other things that must be parameterized in the model). Whatever -- the mixing should not be that far from molecular in this case. And it is far. Happy to help look for the issue with someone who understands the mixing code (and it may be the result of my stupidity in the end).

Jamie

[pmaccc]

Very interesting. What was AKT_BAK for your simulations?

[jpringle]

ALT_BAK is 1.0d-6 for both T and S, so much less than observed... ¯\_(ツ)_/¯

[wilkin]

Do you know about this post regarding minimum TKE ...
viewtopic.php?f=17&t=2310&p=8441&hilit= ... ully#p8458

[jpringle]

Both in John Wilkin's reply and in some private emails, I have received some good hints related to minimum turbulence level parameters in the mixing scheme. I shall investigate and compare to the primary literature for this scheme in the next few days, and will post a summary here soon.

Thanks for the help everyone!
Jamie

[jcwarner]

i have been out of touch for a awhile with some family losses but trying to get back to work.
When I coded GLS into ROMS in 2003 (or whenever) I spent a considerable effort to compare my code to GOTM. I do remember some small differences to GOTM, but essentially our implementation gave the same results as GOTM (at that time). If you look in our GLS paper OM 2005, there are several tests, one that shows mixed layer deepening as compared to the Price solution, and the results were similar to the analytical solution. So i dont think we have any issues with GLS in ROMS, but i could be wrong.

When i saw your cpp defs, i noticed you used CANUTO
#if defined GLS_MIXING
# undef KANTHA_CLAYSON
# define CANUTO_A
# define N2S2_HORAVG
# define RI_SPLINES
#endif

CANUTO will allow mixing up to a Ri of about 1, but Kantha Clayson stops at about 0.2. Suggest you try GLS with KANTHA_CLAYSON.

If you run MY25, you only have the choice of KANTHA_CLAYSON or GALPERIN (if KC is not defined). So it is not really a fair test to compare GLS with CANUTO to MY25 with KANTHA_CLAYSON.

I do not have a lot of experience with LMD so i will leave others to comment on that.

-john

[jpringle]

Thanks to John Wilkin, John Warner, Nirnimesh Kumar, the excellent Warner et al. 2005 Ocean Modeling paper on mixing schemes in ROMS and the excellent note by Malcolm Scully (viewtopic.php?f=17&t=2310&p=8441&hilit= ... ully#p8458 ), I have run this issue down.

In summary:

• Since Warner et al. (2005) the default value of k_min in GOTM (at least in their entrainment test case) has decreased from 7.6e-6 to 1.0e-10. For the reason givien immediately below, it is not useful to ask what the "right" value is.
• As pointed out by Nirni Kumar and Malcolm Scully, when the water is stratified and there is no local turbulence production (i.e. Ri>>0.25 or 1), the vertical mixing of T and S limits to AKt ~ 0.05*TKE*N^-1 . When there is no turbulence production, background turbulence kinetic energy (TKE) for the GLS turbulence closure schemes is given by GLS_Kmin from the input parameter list.
• What is the "right" value of GLS_Kmin? For an actual ROMS run it will be determined by unresolved sources of TKE so is application dependent. E.g. in the stratified coastal ocean I would worry about un-captured internal wave breaking events. For my application I would tend towards the smaller value of 1.0e-10 because I don't think the background stratification will go away as quickly (in 10s of days) as it does with the larger (7.6e-6) value of GLS_Kmin.
• John Warner points out in his message above that the default ROMS parameters agree with the default parameters in GOTM for the wind driven entrainment experiment. Leaving aside changes in the GOTM value of GLS_Kmin (I find references to it being 7.6e-6, 1e-9 and 1e-10 in various places; it is 1e-10 now), the evolution of the mixed layer depth in this experiment is not very sensitive to GLS_Kmin. I have estabilished this experimentally in ROMS and GOTM. Why? Because the mixed layer is a region of TKE generation, so GLS_Kmin is not important in this wind-driven test.
• Fundamentally, this entire thread is not a question of turbulence parameterization schemes. It is a discussion of what happens when there is no turbulence generation (Ri big).

Hernan -- I do think that the default value of GKS_Kmin should be smaller. In the absance of forcing, do we expect stratification to erode away from 50m of water in 30 days? I think the background diffusivity should be close to molecular or an order of magnitude bigger to fit most observed values. It is application dependent, but the current default seems very big.

Jamie

p.s. for reference, here are plots of the evolution of the stratification of the initial condition described above using k-epsilon from GOTM, ROMS k-epsilon with default parameters (labeled k-epsilon on figure below), ROMS k-epsilon with GLS_Kmin=1e-10, and LMD. Images are shown for 5, 10, 20 and 40 days.

[bperfect]

Thanks for the great discussion, everyone. I have been attempting to quantify the turbulent dissipation for an application involving topographic wakes and I recently noticed that the suggested Kmin value for the GLS implementation of k-epsilon was resulting in extremely large background dissipation values. I didn't notice any severe distortion of the background stratification like Jamie, but that could be because of differences in the domain dimensions.

I was hesitant to simply reduce the Kmin parameter because my understanding was that very small values of TKE could cause numerical problems whenever the reciprocal of k is needed. However, reducing Kmin to 1e-10 would likely fix my problem (in the sense that the signal that I'm interested in would dominate the background dissipation).