How to calculate turbulent shear rate / energy dissipation

[pwallhead]

Hi

I would like to calculate (and output) the turbulent shear rate [s-1], which I think means first calculating the turbulent energy dissipation (then dividing by kinematic viscosity and taking square root). I can see how this is done in the COAWST sediment dynamics code, using Warner et al. (2005) Eqn (12) (see sed_flocs.F). But how should I do it in the case of Mellor-Yamada scheme (MY25_MIXING)? It seems like the COAWST flocculating sediment code requires GLS_MIXING and would not work this MY25_MIXING --- unless I missed something?

Any tips greatly appreciated,
Cheers
Phil

[jcwarner]

the turbulence and the sediment should be separate.
Both MY25 and GLS will output tke and gls.
IF you use MY25, then tke is 0.5q^2 (i think, it has been a while) and gls will be the q^2l.
IF you use GLS , then tke is tke and gls will be the set depending on the p m n parameters (gls can be dissipation, or omega, etc).
does that help?

I am not sure if the floc code requires gls mixing. i dont think it would.
-j

[wilkin]

Phil,

The GLS k-kl closure is effectively Mellor-Yamada 2.5.

Here is my Python code for calculating dissipation and turbulence length scale from the output, reading the GLS p,m,n parameters from the netcdf file (so this is generic for any chosen closure). ds is a DataArray opened with xarray (I was actually using xroms, which has convenient features for )

Code: Select all

# Calculate dissipation rate and turbulence length scale
# GLS parameters p, m, n and c_mu0 are saved in output file so are part of ds
# Dissipation and length scale are given by Warner et al. equations (12) and (15) 
#   where k is variable tke in ROMS output
#   and psi is variable gls in ROMS output
# The variable gls takes on different definitions depending on the closure scheme ...
#   k-epsilon, k-kl, etc.
p = ds.gls_p
m = ds.gls_m
n = ds.gls_n

# epsilon = C^(3+p/n) k^(3/2+m/n) psi^(-1/n)
dissipation = np.power(ds.gls_cmu0,3+p/n)*np.power(ds.tke,3/2+m/n)*np.power(ds.gls,-1/n)
ds["dissipation"] = dissipation
ds.dissipation.attrs["long_name"] = "dissipation"
ds.dissipation.attrs["units"] = "meter2 second-3"
# ds.dissipation.attrs["grid"] = ds.gls.grid
# ds.dissipation.attrs["time"] = ds.gls.time

ds["logdissipation"] = np.log10(ds.dissipation)
ds.logdissipation.attrs["long_name"] = "log10 dissipation"

# l = C^3 k^(3/2) epsilon^-1
lengthscale = np.power(ds.gls_cmu0,3)*np.power(ds.tke,3/2)*np.power(dissipation,-1)
ds["lengthscale"] = lengthscale
ds.lengthscale.attrs["long_name"] = "length scale"
ds.lengthscale.attrs["units"] = "meter"

[pwallhead]

Thanks John and John!

I checked the ROMS code to compare how the turbulent energy dissipation is calculated under MY25_MIXING (my25_corstep.F) vs. GLS_MIXING (gls_corstep.F). It seems that the MY25_MIXING code is indeed consistent with the GLS_MIXING formulation for epsilon (Warner et al. 2005 Eqn 12), for the special case of (k-kl) scheme (p=0, m=1, n=1), *assuming* a fixed (hard-coded) value of gls_cmu0 = 0.4939. My rationale is below --- it would be great if someone could cross-validate...

In both files, the implicit Euler step uses 'BCK' to set the denominator for the tke update. In gls_corstep.F the Warner et al. 2005 formula enters one of the terms contributing to BCK (losing a factor of tke due to the implicit update). In my25_corstep.F the corresponding term is 2.0_r8*Qdiss, where Qdiss includes the dependence on tke and gls. To make this latter consistent with gls_corstep.F I think we must have (2/my_B1) = (2/16.6) = gls_cmu0^3, hence gls_cmu0 = 0.4939. This makes me think there might have been a mistake in the code for my25_corstep.F, because a more canonical value of gls_cmu0 is (2^1.5/16.6)^(1/3) = 0.5544 (Warner et al., Table 2). Hence I wonder if the my25_corstep.F code should read 2.0_r8*SQRT(2.0_r8)*Qdiss rather than 2.0_r8*Qdiss?

In any case I think the code in sed_floc.F where Warner et al. Eqn 12 is applied (under FLOC_TURB_DISS) will not work in the case of MY25_MIXING because it uses the constant gls_cmu0(ng), which is set only if GLS_MIXING is defined. Unless I missed something ...

Phil