pkg/kl10 density calculation bug?

[loisbaker]

I've found some funny behaviour of the kl10 parameterisation. The package is supposed to enhance mixing and diffusivity when there are density overturns, but it's turning on in my simulation for areas that are stably stratified. I think it could be due to how density profiles are calculated and sorted (then overturn length scale inferred):

In-situ density is used to find the unsorted density profile (call to FIND_RHO_SCALAR on line 93 of pkg/kl10/kl10_calc.F). I don't think that this should give any overturns at all due to the added vertical density gradient from using in-situ density rather than neutral density. However, kl10 does find overturns because the call to FIND_RHO_SCALAR (line 93 of pkg/kl10/kl10_calc.F) gives totPhiHyd as the pressure argument, rather than totPhiHyd*rhoConst. So density is being found as if it were at (say) 2 m depth rather than 2000 m depth. This gives a diffusivity/viscosity output that looks mostly sensible in terms of overturn location, but can turn on very large areas of high diffusivity/viscosity where it thinks that the (stable) density profile is unstable. Interested to know if my interpretation is correct!

[jklymak]

Hi @loisbaker its possible that kl10 was never tested on a nonlinear equation of state? If so that would be a pretty large error for anyone trying to use kl10 that way. I don't find how to calculate density in the mitgcm to be particularly clear but I thought I copied the routine from elsewhere in the code, so your hypothesis should be tested. It's also the case that density calculation changed from when I originally wrote this so a bug could have crept in.

[loisbaker]

Hi @jklymak - thanks for the reply, that makes much more sense if it's only supposed to be used with a linear equation of state! Otherwise, some neutral density would need to be used instead of in situ density. What's the best way to make this clear to future users?

It looks to me like other calls of FIND_RHO_SCALAR in the code use a pressure rather than a pressure potential (p/rho), I could be missing something though.

[jklymak]

Does the mitgcm have a way to calculate a potential or neutral density? Potential density is probably ok for an overturn I'd expect though there may be small errors.

Certainly if this is a practical limitation that we can't overcome, the documentation should be fixed and/or a check out into the code to warn folks.

[loisbaker]

The vertical gradient of iso-neutral density is available, so I'm having a go at editing the function to use that instead of in-situ density profiles. I'll update with any progress!

[loisbaker]

I've just had a chance to come back to this @jklymak. I think KL10 can easily be made to apply to a nonlinear EOS by replacing lines 93-95 pkg/kl10/kl10_calc.F:

CALL FIND_RHO_SCALAR(theta(I,J,K,bi,bj), salt(I,J,K,bi,bj),totPhiHyd(I,J,K,bi,bj),rhot,myThid )
rhoS(K)=rhot

with:
rhoS(K)= rhoS(K-1) + rkSign*drC(K)*SigmaR(I,J,K)

Where SigmaR is the vertical gradient of iso-neutral density.

I've checked that this amounts to the same thing in the linear EOS case, but it does result in a quite different (although qualitatively similar) output in the internal_wave verification example after a few hundred timesteps. I'm fairly sure this is due to errors in the density calculation that start at 1e-12 then get progressively larger - if this sounds like expected behaviour, I'm happy to create a PR.

[jklymak]

@loisbaker any help would be appreciated. I don't run simulations with non linear equation of state, but this should be made to work with them.