Richardson number profile output jumps around and is not smooth

[xkykai]

The Richardson number profile in the strong wind case and free convection case doesn't seem to be correct, https://engaging-web.mit.edu/~alir/lesbrary/6DaySuite/three_layer_constant_fluxes_linear_hr144_Qu0.0e+00_Qb5.0e-08_f1.0e-04_Nh256_Nz128_free_convection/instantaneous_statistics.mp4 https://engaging-web.mit.edu/~alir/lesbrary/6DaySuite/three_layer_constant_fluxes_linear_hr144_Qu7.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_strong_wind/instantaneous_statistics.mp4 The profile looks very jagged in regions were it should be smooth, and jumps around sometimes.

[glwagner]

For free convection the Richardson number isn't well defined (if computed from horizontal averages) because the mean shear is close to zero, right?

[xkykai]

yeah, since the mean shear is close to zero, shouldn't Ri be a huge number in the region? It looks like at certain grid points it drops to below zero. That can be seen in the strong wind case as well.

[glwagner]

True, I guess I would expect it to be uniformly large. Certainly dz(B) is never < 0, so Ri should be strictly positive. Can you screen shot a problematic frame in the movie and paste it here?

[xkykai]

Here you go!

[glwagner]

Right so while Ri attains very large values as expected, it also appears to attain low values below 0 which is unexpected. Where is Ri calculated?

[xkykai]

https://github.com/CliMA/LESbrary.jl/blob/371ce34cbc4352c58ef40358e993b3797acc9d33/examples/three_layer_constant_fluxes.jl#L407 It's here! I have been using this script to run simulations.

[glwagner]

and where is richardson_number_ccf! defined?

It looks like this might be the horizontal average of the local (3D) Richardson number, rather than the ratio between dz(B) and (dz(U)^2 + dz(V)^2), where the capitalized variables are horizontally-averaged. In that case, it is possible that Ri is less than zero because locally Ri can take on values close to - infinity?

[ali-ramadhan]

Ah initially looking at the movie I thought there might be a bug since the Ri(z) goes wonky in between frames, but yeah now that you mention it, it's actually taking the horizontal average of the local Richardson number (which could container some very large positive and negative numbers) which seems like the wrong approach.

Actually even worse, it was using averaged U(z), V(z) but local b(x, y, z) haha.

Should be fixed in https://github.com/CliMA/LESbrary.jl/pull/113

and where is richardson_number_ccf! defined?

https://github.com/tomchor/Oceanostics.jl/blob/c7c6a8bea46d127f243e7c8f610edb7ce03ffcec/src/FlowDiagnostics.jl#L11-L19

[glwagner]

So just to clarify the discussion, I'll mention that there are many definitions of the "Richardson number", and the physical interpretation of this quantity depends on how it is computed.

Most people define a "gradient" Richardson number as

Ri = dz(b) / (dz(u)^2 + dz(v)^2)

There is also a "bulk" Richardson number that can be interpreted as a ratio between potential energy and kinetic energy.

The first thing to mention is that the gradient Richardson number is fundamental to the instability of stratified shear flows. Instabilities often occur when Ri < 1/4. Note that the exact stability criterion depends on the entire flow profile and Ri < 1/4 is strictly only informative for linear shear and stratification. This is a subject of vigorous research in fluid dynamics.

The relative simplicity and predictive power of Ri < 1/4 means that it is often used as a "proxy for mixing" in stratified turbulence studies. Thus in DNS and LES it is common to measure the local Richardson number and interpret turbulent regions / behavior in light of this diagnostic quantity. This is the sense in which the quantity is defined in Oceanostics, where it is defined in three-dimensions, locally, and with the option to include a background flow around which the dynamics might be posed.

In parameterization, we use a quantity that is superficially identical to the local gradient Richardson number. However, this quantity is different than the local gradient Richardson number because it is defined on coarse-grained / horizontally-averaged fields only. For "exact" flows the criterion Ri < 1/4 is paramount. For coarse-grained fields, the specific Ri value is just an indicator the true, underlying subgrid distribution of Ri. Another interpretation of the coarse-grained Ri is the relative role of potential energy dissipation and kinetic energy dissipation in a turbulent flow.

Both the local and coarse-grained Ri are useful quantities. I'm not sure what we want to compute here. However, it is definitely the case that the horizontal average of the local Ri may not be very useful. The issue is that we lose information by horizontal averaging; and it is specifically the local value of Ri in that case which is useful. This is evident in the results here --- apparently, there are very localized parts of the stably stratified region below the boundary layer that experience very small Ri (and are thus susceptible to turbulent mixing). However, despite that these regions of low Ri are localized, when a horizontal average of Ri is taken we see that the horizontal average for the entire slice is less than 0. This shows that we have lost information about the fact that most of the slice is stable and non-turbulent.

Probably what we would like to compute is just dz(B) / (dz(U)^2 + dz(V)^2), where U, V, B are horizontal averages. There's no need to launch a kernel to compute these. We don't have a great way of expressing this computation in the current abstract operations framework. Perhaps this is something to work on.

[ali-ramadhan]

Thanks for the explanation! I agree that Ri = dz(B) / (dz(U)^2 + dz(V)^2 seems to make the most sense in our case (boundary layer parameterizations).

Will reply to your last paragraph in PR #113.