Transition to small w^2 regimes

[vanroekel]

The model currently struggles massively when w2 becomes small. This is mostly due to tumd terms which are defined as tumd = wt/(Mc+1e-12). When w2 goes to zero Mc goes to zero, this causes tumd to go to infinity and the model becomes unstable. It makes me wonder how we should handle this case. I see a couple solutions

1. we enforce a minimum in w3 (say 1e-8 or 1e-7)
2. we somehow transition to traditional high order closure forms for terms with tumd and sumd in them. w2t and w2s are fine as sigma=0.5 sets those terms to zero.

I'm a bit uncomfortable with (1) as it is non conservative in energy, but I guess we could think of this as an artificial background diffusivity, where kappa = w2*tau

I'm curious what others think.

[vanroekel]

pinging @BrodiePearson, @qingli411,

[BrodiePearson]

@vanroekel Where does the model typically encounter these small w2 values, is it at the entraining interface, at the surface, or somewhere else? And does your first suggestion refer to setting a minimum value for w2 or for w3?

If the problem is at the interfaces then I might interpret that as the regions where local (sub-plume-scale) processes dominate mixing so the plume formulation is not suitable. Perhaps something similar to your first point would be appropriate with a a coupling between the background diffusivity and the sub-plume scale TKE. Thinking about the flux being either non-local (plume-driven) or local (sub-plume driven):
$\overline{w'\theta'}=max(Mc*\theta_{umd},-\kappa\partial_z\theta)$
where $\kappa \propto h\sqrt{e_{sps}}$ with $h$ the boundary layer depth or distance from the surface, and $e_{sps}$ the sub-plume scale TKE. If the first (plume) term is bigger then you diagnose $\theta_{umd}$ as you already do, but if the second (diffusive) term dominates you set $\theta_{umd}=0$ at that level?

[vanroekel]

@BrodiePearson interesting way to think about this. It makes sense to me. To clarify, the min value is on w2 and is indeed at the entraining interface.

This makes me wonder though if we should for terms like

w2tend4(k,iCell) = (2.0_RKIND - 4.0_RKIND/3.0_RKIND*C_2)*Mc(k,iCell)* &
                  (grav*alphaT(k,iCell)*tumd(k,iCell) - grav*betaS(k,iCell)*sumd(k,iCell))

should also include a SPS contribution, i.e.

w2tend4(k,iCell) = (2.0_RKIND - 4.0_RKIND/3.0_RKIND*C_2)*Mc(k,iCell)* &
                   (grav*alphaT(k,iCell)*tumd(k,iCell) - grav*betaS(k,iCell)*sumd(k,iCell)) + &
                   g*alpha*areaFraction(k,iCell)*wt_spsD(k,iCell) + (1.0_RKIND - areaFraction(k,ICell))*wt_spsU(k,iCell) + g*beta*salinity

and in similar places in u2 and v2

[BrodiePearson]

@vanroekel Is the code breaking because of the tumd (and sumd) source terms in the w2 equation that you posted or because of the wt budget terms which depend on tumd^2?

Since this is a plume-scale budget, I'm not convinced that sub-plume scale fluxes belong (but I'm still thinking about this). I'll take a more thorough look through the equations over the next few days and see if I spot anything.

I also realized that the limiter I suggested cannot be imposed practically, so perhaps some transition that applies under small $\frac{\overline{w'w'}}{\frac{2}{3}e_{sps}}$ values could be applied (i.e. plume formulation is not applicable)

[vanroekel]

@BrodiePearson I'm not actually sure what the fundamental issue is with the code. It has been very difficult to tease out. But my feeling is that the model is struggling to handle situations were w3 is transitioning from slight positive values to negative values (as entrainment happens). As it goes through zero, something is causing w2 to be zero at some layers, causing blow ups. What I tend to see first though is a somewhat 2dz pattern in w2 and w3. I'm still trying to dig into the data.

Thinking on this more, I think you are right that the SPS terms should not appear. Fundamentally, it seems when we are strongly turbulent we expect w't' from the plume scale budget to far exceed the SPS contribution and vice versa when turbulence is weak. It may boil down to finding the right condition for when to set tumd = 0.

I agree as well on the w2t limiter. I don't think it is able to be imposed, and in fact the code is more stable with out it. I've added a different realizability condition on w3, where when we have sigma get clipped (greater than 0.99 or less than 0.01) w3 is limited such that it is consistent with the clipped sigma value. This has helped some, but hasn't fixed the fundamental problem of model stability.