Possible inconsistency in the eddy definition of diffusivities/closures in the docs

[tomchor]

Maybe I'm missing something, but when I look at the equations in the closures page it seems that ν_eddy does not include the molecular viscosity ν. However looking at this closure page in the model physics it seems like the molecular viscosity is included in the Smagorisnky and AMD eddy viscosities.

Looking at model.diffusivities.ν_e it seems like that the lowest value is always the molecular (background) viscosity, so I'm guessing that the convention used by the model is the latter one. Should those be consistent with each other?

Thanks!

[glwagner]

You're right that the molecular viscosity is added to νₑ:

https://github.com/CliMA/Oceananigans.jl/blob/81db22f4a26396142e8cd5b5a4c50c75790c1d50/src/TurbulenceClosures/turbulence_closure_implementations/verstappen_anisotropic_minimum_dissipation.jl#L122-L136

It seems nicer to me to keep these separate in the docs, but it's less expensive / convenient to combine them in the code.

We could define νˢᵍˢ and νₑ = νˢᵍˢ + ν in the docs to solve this issue.

It's also been pointed out that making the default value of ν to be a particular number for a particular temperature water is weird. We should probably have default ν = 0.

[glwagner]

@tomchor I think you should feel free to edit the docs (and also make other appropriate clarifying changes if they are appropriate). I think we can call νₑ the "effective diffusivity" and νˢᵍˢ the subgrid or subfilter diffusivity.

[tomchor]

We could define νˢᵍˢ and νₑ = νˢᵍˢ + ν in the docs to solve this issue.

On the surface I definitely agree with this. Although I prefer to keep calling ν_e the eddy diffusivity, since that's what everyone calls it and it's supposed to model the action of eddies. The caveat here is that this doesn't match the convention that the code uses, right? Which I think can be pretty confusing.

It's also been pointed out that making the default value of ν to be a particular number for a particular temperature water is weird. We should probably have default ν = 0.

This seems like a matter of preference, but I kinda like that the default is a reasonable value. Saves me the trouble of looking up what the molecular values are if I'm simulating something :)

[glwagner]

This seems like a matter of preference, but I kinda like that the default is a reasonable value.

Ok, let's think about it. It also doesn't make sense if what you're doing is non-dimensional...

I've found that AMD doesn't actually work that well when turbulence is weak and the subfilter diffusivity is "close" to the molecular value. It might be more symbolic than anything to have a non-zero molecular diffusivity, since setting it to zero might not have an observable effect on results ... ? What do you think? Have you found that it changes your results?

One downside of specifying a non-zero yet negligible number (in other words, the difference between this number and zero cannot be observed) is that it implies that non-zeroness is important when it isn't. So if it really doesn't matter, I'm worried that it distracts from the communication of scientific results...

The default for IsotropicDiffusivity is also this "20 degree water" number. Maybe that's more important to change, since simulations with IsotropicDiffusivity in Oceananigans are probably idealized rather than DNS of a miniscule box.

[glwagner]

Although I prefer to keep calling ν_e the eddy diffusivity, since that's what everyone calls it and it's supposed to model the action of eddies. The caveat here is that this doesn't match the convention that the code uses, right? Which I think can be pretty confusing.

You've pointed out here that the docs are inconsistent, so they must be changed. But we could change notation in both the code and the docs if we think it will make either one more clear.

We could also change the closure implementation so that the contribution of the molecular diffusivity to flux divergences is added separately, rather than adding the molecular diffusivity to νₑ in calculate_diffusivities!.

[tomchor]

You've pointed out here that the docs are inconsistent, so they must be changed. But we could change notation in both the code and the docs if we think it will make either one more clear.

I'm not sure what the best approach is for the code. Also, I'm happy to change the docs to be consistent after we agree on what to do, but I'd rather not be the one changing the code right now (if we do decide to change it). That would require me to spend a considerable amount of time learning how to contribute to Oceananigans given that I've never started a project in Julia or contributed to a repo that wasn't my own).

[tomchor]

You've pointed out here that the docs are inconsistent, so they must be changed. But we could change notation in both the code and the docs if we think it will make either one more clear.

So, I thought about this a little more, and I do think that the notation used in the code and the docs could change to something more accurate and consistent. Right now, if I understand correctly, the code considers ν_sgs to model the action of eddies, v to be molecular diffusion, and v_e to be the sum of both:

https://github.com/CliMA/Oceananigans.jl/blob/81db22f4a26396142e8cd5b5a4c50c75790c1d50/src/TurbulenceClosures/turbulence_closure_implementations/verstappen_anisotropic_minimum_dissipation.jl#L135

However, a better way of naming thing (which both the docs and code could follow) would be

• take ν_e (or ν_eddy if you prefer) to be the turbulent closure alone. Because that is meant to model exclusively the action of eddies, so the name makes sense.
• ν could continue being the molecular diffusivity (we could even change it to ν_m or ν_molec to be more explicit)
• ν_sgs is then ν_eddy + ν_molec, since this is generally what turbulence literature uses. And it makes sense, since both unresolved eddies and molecular diffusion are processes that happen at the sub-grid scales.

This would require changing nomenclature in both the code and the docs, unfortunately, but the upside is that is doesn't require any change in the programming itself (it would have to change if we separated eddy and molecular viscosities).

Thoughts?

[francispoulin]

I have not used eddy diffusivities yet so I can't say that I have had any experience with this myself, but based on commments from @tomchor, I will see that I like using the names ν_eddy, ν_molecular and ν_sgs or ν_total to be the eddy, molecular and total diffusivities.

Having said that, I don't imagine that for large scale flows people will ever need to use the total as the eddy viscosities would dominate.

[glwagner]

take ν_e (or ν_eddy if you prefer) to be the turbulent closure alone. Because that is meant to model exclusively the action of eddies, so the name makes sense.

I think the e could be "effective" rather than "eddy". I'm not sure it makes sense to refer to the molecular diffusivity as "sub-grid scale". Shouldn't we refer to the molecular contribution to the flux as "resolved" physics?

[glwagner]

Having said that, I don't imagine that for large scale flows people will ever need to use the total as the eddy viscosities would dominate.

Eddy viscosities dominate in almost all ocean problems, including large eddy simulation with grids on meter scales --- the exceptions being "resolved LES" of the bottom boundary layer, where in some circumstances (very small problems) the molecular contribution can be significant. For typical ocean conditions we have found that the molecular contribution is less than a percent in most cases (which may be less than numerical diffusion...)

[jklymak]

I think it is important to differentiate whether the eddy diffusivity has the "molecular" as a lower bound or if it is additive. Folks run with a-physical molecular diffusivities and viscocities all the time (ie low Reynolds number) I know MITGcm will usually choose the max of different physics rather than add them.

[glwagner]

The main LES closure we use (Anisotropic Minimum Dissipation) is supposedly designed so that the "eddy" or "sub-grid scale" component vanishes to 0 in laminar flow. Thus when using AMD in transitional problems where laminar flow and turbulent flow with unresolved scales coexist, it'd be appropriate to add the molecular diffusivity to the nonlinear sub-grid scale contribution. But in practice I think LES in transitional flows is pretty tricky... so I think this consideration might be more of a theoretical one rather than something of pressing practical importance...

The other abstraction we have designed for using multiple turbulence closures in one problem is also additive. However, we could easily design an alternative abstraction that uses the "maximum" (or minimum). Perhaps a turbulence closure called MultipleTurbulenceClosures(combination_method, closure1, closure2, ...) would be useful, where combination_method might be Sum(), Maximum(), Minimum(), etc.

[tomchor]

I think the e could be "effective" rather than "eddy". I'm not sure it makes sense to refer to the molecular diffusivity as "sub-grid scale". Shouldn't we refer to the molecular contribution to the flux as "resolved" physics?

I come from an engineering background, so my views might be different, but I don't remember ever seeing something called "effective" diffusivity. Also, to me at least it's kind of ambiguous what process go into some "effective" diffusivity, which I think can confuse users (I certainly would be confused).

Plus, one could argue that molecular diffusivity itself isn't technically "resolved physics", since we're not resolving the molecules! It's just the same down-gradient closure as the eddy diffusivity. It works really well in this case because there's a huge scale separation between the molecule size and the smallest motion size, which is the Kolmogorov scale. While this scale separation is natural for molecular diffusivity, it's completely artificial (non-existent really) for the eddies, which is why it doesn't work nearly as well. (Corrsin published a really nice paper about this in 1976.)

So I still think that ν_sgs = ν_eddy + ν_molec makes sense and is unambiguous. I think being unambiguous is the most important part here, though, since we want people to understand the code easy to be able to contribute, no?

[tomchor]

Perhaps a turbulence closure called MultipleTurbulenceClosures(combination_method, closure1, closure2, ...) would be useful, where combination_method might be Sum(), Maximum(), Minimum(), etc.

I like this idea!

Also, like people pointed out, the treatment of molecular viscosity won't really make any quantitative changes in 99% of the simulations run I think. But I think it's good to get it "right" because (i) we should be striving for consistency in the physics no matter what, (ii) it may make some differences if a precise closure of the TKE equation is necessary, and (iii) Oceananigans is so general that it is possible to simulate things at very small scales (small enough that the transition between turbulent/laminar flow may be important).

[glwagner]

I agree that "eddy viscosity" is common and "effective viscosity" is not, and abusing the letter "e" could lead to confusion.

However, I also think that ν_sgs is synonymous with "eddy viscosity" and that the molecular fluxes are not "sub grid" or "sub filter".

The most important question is whether or not constant isotropic component of viscosity / diffusivity should be included in the arrays saved as model.diffusivities. What are people's opinions about that?

If would like to save the total viscosity there then I think we should call it just that (ν_total), or find a comparably simple name.

If we would like to include just the "eddy" or nonlinear viscosity in model.diffusivities, then we can keep the current name in the code νₑ, and change the physics docs so that they are consistent with the code. The physics docs are currently consistent with the code -- its just that we don't want to call "sgs + molecular" as νₑ or "eddy".

As pointed out in the original post, the numerical implementation docs are out of sync with both the physics docs and the code. I think this is because they were essentially written as a stand-alone tutorial to large eddy simulation rather than something that pertains specifically to Oceananigans. I think a lot of the numerical implementation docs are a bit dated. It might be better to nuke this section than confuse future users. But an alternative is to spend some time sprucing it up. More than just the large eddy simulation section needs to be improved.

[francispoulin]

Lots of interesting discussion here. Nomenclature can be a pain as people use a lot of different words for the same things and sometimes use the same word for different things. I wonder whether finding a good review article and citing that would be a good way to establish our foundation? I imagine there is a lot to choose from out there but don't have a specific one in mind.

Below are a few thoughts.

1) I would argue that molecular viscosity is real (resolved) physics, in these continuum models, is derived from physical principles and you should not change the coefficients based on the grid you are using. These are measured in experiments and are what the community to believe as true when doing lab experiments.

2) For large-scale flows, as many of us consider, molecular viscosity is insignificant and should probably be ignored as it has no real meaning on the large scales. That is why summing molecular and eddy viscosities does not seem like such a good idea to me. However, adding a very small number to a large number will nto make much of a difference for all intensive purposes. I would think that the user has the onous of responsibility of picking the viscosity that is appropriate for their problem.

3) For large-scale flows where you need some parameterization, I think that using eddy or sgs are accepted by large groups in the community. Effective seems a little less desirable for me as it does not say where it comes from.

[tomchor]

The discussion about molecular viscosity being resolved/modeled physics is a deep one. Probably too deep for a github issue. After many discussions with Jim McWilliams he managed to convince me that the only difference between the modeling of eddy versus molecular diffusion is scale separation between processes.

The issue about always keeping a more rigorous set of equations and including the molecular diffusivity by default (even though it's probably not important in the vast majority of cases) is one of many modeling philosophies. I highly recommend the paper by Kerry Emanuel. It's only an opinion piece, so obviously not the objective truth, but I quite agree with it. This passage stuck with me:

Sometimes the quest for better simulations subordinates even simple physics. About 20 years ago I pointed out that most models of that era neglected to turn dissipated kinetic energy back into heat. For most atmospheric phenomena, this is indeed a small term in the thermodynamic energy budget (though technically required to close any net energy budget), but in strong windstorms like hurricanes, it becomes important. Moreover, no substantial computational benefit accrues from neglecting it. A few weeks later, a researcher came to me to report that he had added this term to his model and found that it made simulated hurricanes too intense, so he took it out again. [...] For this researcher, getting the “right answer” was the goal, even if it is obtained for the wrong reasons. Still today, the conversion of dissipated kinetic energy back into heat remains an optional switch (whose default position is “off”) in a state‐of‐the‐art hurricane prediction model.

I'm definitely not accusing anyone of doing this here! But I just wanted to explain a bit where I'm coming from.

That said, I agree with @glwagner that maybe we should focus on the more pressing issues at hand:

• 1 - Should Oceananigans keep adding the molecular viscosity as a background to the eddy viscosity? My vote is yes since it's more inline with the physics, which are additive, and we don't gain significant computational efficiency from neglecting it. But that's just my vote. (I also like calling it ν_total, as it is unambiguous.)
• 2 - What should we do about the numerical implementation docs being out of sync with both the physics docs and the code? I'm happy to change the docs myself, but I think we haven't converged on how to do that yet. I'm okay with nuking one of the docs sections. I don't see why two separate sections are needed and the physics / turbulence closures section is more than enough to explain to the user what's going on IMO.

[jklymak]

For large-scale flows, as many of us consider, molecular viscosity is insignificant and should probably be ignored as it has no real meaning on the large scales. That is why summing molecular and eddy viscosities does not seem like such a good idea to me. However, adding a very small number to a large number will nto make much of a difference for all intensive purposes. I would think that the user has the onous of responsibility of picking the viscosity that is appropriate for their problem.

I agree with this, however, there are whole classes of flows where you use an eddy diffusivity from some sort of scheme based on resolved physics and have a higher background diffusivity to represent diffusivity due to processes that are not resolved. The only way you can get away from doing that is if you are going to claim you are always resolving the outer length scales of whatever turbulence may be in your flow, and that is really not practical a lot of the time.

The classic situation in the ocean is when you are resolving km-scale flows, where Smagorinski may be a perfectly reasonable lateral closure, but you still have an internal wave continuum that you are not even close to resolving. The typical kludge is to realize that you have that mixing in the flow, and increase the "background" diffusivity.

[glwagner]

Should Oceananigans keep adding the molecular viscosity as a background to the eddy viscosity? My vote is yes since it's more inline with the physics, which are additive, and we don't gain significant computational efficiency from neglecting it.

I think my suggestion is being misinterpreted --- I was referring only to changes to software and data structures; not changes to the partial differential equations being solved.

If we do not "add" the constant background viscosity component to model.diffusivities.νₑ, we would instead add the contribution of the constant background viscosity to the total diffusive flux divergence separately from the contribution associated with the nonlinear eddy viscosity (similar to how the diffusive flux divergence for constant background viscosity is currently calculated without the auxiliary field model.diffusivities.νₑ).

Thus my suggestion only changes the content of model.diffusivities.νₑ and otherwise preserves both the user interface and the default behavior for Oceananigans.TurbulenceClosures.AnisotropicMinimumDissipation.

We could also completely remove "constant background" parameter from AnisotropicMinimumDissipation and preserve existing behavior, since users can use multiple turbulence closures to achieve the same effect, as in

closure = (AnisotropicMinimumDissipation(), IsotropicViscosity(ν=1.05e-6))

for example. This might be a better user interface than what we have now.

As a side note, since there's no constant background viscosity that's appropriate for all scenarios, it should probably be zero by default if we continue to allow users to set it as a parameter in AnisotropicMinimumDissipation (eg #1278).

[glwagner]

I'm thinking about introducing a new container for "auxiliary fields" called model.auxiliaries, getting rid of model.diffusivities, and putting references to the "eddy diffusivity" in the LES closure parameter structs.

Are we okay with removing the constant background viscosities and diffusivities from the LES closures, and using closure tuples instead to implement this behavior? Much is simplified this way I think, because we can now have closure.νₑ straightforwardly. I think scripts are also clearer because users would write:

closure = (AnisotropicMinimumDissipation(), IsotropicDiffusivity(ν=1.05e-6))

To implement an LES closure in addition to a constant background isotropic viscosity.

[tomchor]

I agree!