Problematic behavior of the RiBasedVerticalDiffusivity

[simone-silvestri]

with @sandreza we found out that the new RiBasedVerticalDiffusivity has some stability issues.
These issues are related to the entrainment diffusivity added to the tracer diffusivity.
In particular, the change is

 κᵉⁿ = ifelse(entraining, Cᵉⁿ, zero(grid))

to

κᵉⁿ = ifelse(entraining, Cᵉⁿ * Qᵇ / N², zero(grid))

where

entraining = (N²_above < 0) & (!convecting) & (Qᵇ > 0) 

Previously, a maximum of Cᵉⁿ was added (a constant free parameter with a value of 0.1 calibrated from LES).
When substituting this with a fraction with N² at the denominator, the diffusivity has no upper bound, exploding towards very unphysical values where stratification is low.
This, compounded with the averaging in time, has the effect of having an ever-growing diffusivity which eventually leads to N² = 0 and consequently κᵉⁿ == NaN

There are to possible fix to this:

1. revert to a constant "entraining" diffusivity
2. specify a κᵉⁿ_max and κᵉⁿ = min(κᵉⁿ_max, ifelse(entraining, Cᵉⁿ * Qᵇ / N², zero(grid)))

[glwagner]

Interesting! Good find.

This not the only thing to worry about regarding the entrainment diffusivity. It's effect also seems to depend on vertical resolution, and it only seems useful at low resolutions.

In general, I recommend switching to CATKEVerticalDiffusivity, where we have put a lot of effort into solving those problems.

I'm going to suggest a different solution that setting some maximum capping diffusivity. I think setting a maximum diffusivity could work, but that there are better solutions that more directly address the underlying issue.

I think the issue here is that "entrainment" should not occur across low stratification regions. Thus this model for "entrainment" assumes the stratification is "strong".

More generally, we are dividing the boundary layer into a convecting region and an entraining region. We identify the convecting region by N² < 0. However, in reality the convecting region should extend further down into regions of stable, but still weak stratification below the unstable region. The "entraining" layer is more properly where stratification significantly increases over the weakly stratified lower part of the convecting region. (This analysis is based on LES.) RiBasedVerticalDiffusivity's formulation --- apart from being overly simple --- simply cannot cope with the complexity of real boundary layer structure because it is not truly non-local. Really addressing this issue (eg implementing a more physically accurate model of boundary layer structure) might require using a vertical integral instead of a 3D kernel. (We decided that because CATKE seems like a decent parameterization, it's not worth it to pursue those accuracy improvements).

One way to fix the issue but retaining a 3D kernel is to change the way we estimate "convecting" and "entraining" regions. Those criteria are

Oceananigans.jl/src/TurbulenceClosures/turbulence_closure_implementations/ri_based_vertical_diffusivity.jl

Lines 247 to 248 in 94e99e1

	convecting = N² < 0 # applies regardless of Qᵇ
	entraining = (N²_above < 0) & (!convecting) & (Qᵇ > 0)

If we introduce a new parameter, say minimum_entrainment_buoyancy_frequency, then we can modify this criteria to read

convecting = N² < 0 # applies regardless of Qᵇ 

N²_entrainment = minimum_entrainment_buoyancy_frequency
entraining = (N²_above < 0) & (N² > N²_entrainment) (Qᵇ > 0)

Now the user can determine minimum_entrainment_buoyancy_frequency to stabilize their model. Probably some default like minimum_entrainment_buoyancy_frequency = 1e-7 would work for Earth-like situations.

[glwagner]

And note I think a more satisfying model would assess the buoyancy jump across the entrainment region relative to the "convective potential" in the column above. Then we could say that we're convecting until we reach a grid cell with a significant stable buoyancy jump that's an appreciable multiple of the convective potential above. That would introduce a non-dimensional parameter rather than a dimensional one. But it would also require a vertical integral.