Help with initial velocity profile

[sam12396]

Hi all!

I am having an issue where ~30 minutes into my run, the pycnocline spontaneously mixes and diffuses. This is concomitant with an increase in eddy viscosity and the interface.

I am using the AnisotropicMinimumDissipation closure.

There isn't an instability in my initial density profile so I expect that the issue is related to the velocity noise I initialize with. Although I have tried a number of different initial velocity profiles, the general rule of thumb I have been following is noise in the surface that decays before it gets to the pycnocline.

I am open to any tips but my specific questions are:
1) What order of magnitude should the noise be?
I have mostly been using very low amplitude noise O(10^-4) at the surface and it decays to O(10^-9) at the pycnocline.

2) Can this be an issue of shear in the noise profile?
What do you recommend for a surface to pycnocline velocity noise difference?

As a note, the issue stopped when I was using the ConstantSmagorinsky closure

Thanks for the help!

[glwagner]

The AnisotropicMinimumDissipation closure can only produce non-zero diffusivities if velocity gradients are present:

https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/#Anisotropic-minimum-dissipation-(AMD)-turbulence-closure

What produces velocity gradients in your simulation?

This setup is similar, though I think the "pycnocline" could be thicker than in your simulation:

https://github.com/CliMA/LESbrary.jl/blob/master/examples/three_layer_constant_fluxes.jl

The noise is

noise = 1e-6 * Ξ(z) * dθdz_surface_layer * grid.Lz

where

Ξ(z) = rand() * exp(z / 8)

The amount of noise you "need" depends on your resolution, and the type of simulation you'd like to perform. With higher resolutions you need less noise (and your eddy diffusivity should decrease).

I think the best way to find out is to run a series of simulations with different noise levels, including "no noise" (which may not work at all for you, depending on how the simulation is forced).

Can you describe in more detail the physical scenario that you are trying to model?

[sam12396]

Thank you so much for all the help!

That looks like it is doing the trick, I even tested running for an hour just to be sure! I also did a test with max_Δt=1minute and it worked fine. Most of the steps were on the order of seconds so it was probably still being kept low by the cfl conditions.

I will run some tests tomorrow with forcing turned on.

A few questions so that I don't have to bother you again if this happens:

1. So in the end, the problem was a combination of the initial velocity profile, vertical resolution, and max_Δt?

• If the vertical resolution was an issue, do you think we hit a sweet spot or can I go finer? I was planning on 128 cells over an 85m extent.

2. What is the using Oceananigans.BuoyancyModels: BuoyancyField line doing? It was giving me an error and I saw that BuoyancyField was unused so I commented the line out and everything worked fine.

[glwagner]

Of course! You have a cool problem!

So in the end, the problem was a combination of the initial velocity profile, vertical resolution, and max_Δt?

The experiments I ran suggest the problem was due to the time-step. Your initial condition has an extremely strong stratification. I can't remember the numbers but I think N^2 = 0.01, roughly? This means that the time-scale of a buoyancy frequency internal wave is 1 / N ~ 10 seconds. Perhaps that was simply the issue? (I think I made a mistake when I did this calculation initially --- it would be worth redoing). The time-step restriction imposed by the buoyancy oscillations is independent of resolution (the dispersion relation for buoyancy oscillations is trivial!) so resolution may not play a role. I think the specific time-step restriction, and whether or not it is affected by resolution, is something you are in a better position to answer via systematic experimentation.

I strongly recommend running a suite of simulations to discover the answer to these questions for your own edification, and also to become comfortable with performing systematic parameter sweeps and numerical experiments for the purpose of discovering scientific knowledge! All you have to do is wrap your setup in a function that takes in the time-step and vertical resolution as arguments, and then call it a few times...

I also don't think the velocity profile mattered much. Your simulation was exhibiting a numerical instability that probably would have been triggered by any non-trivial initial condition... (experimentation would answer this more conclusively).

If the vertical resolution was an issue, do you think we hit a sweet spot or can I go finer? I was planning on 128 cells over an 85m extent.

I think we have only figured out how to "initialize" your simulation stably. Your scientific application has different requirements. I don't think we know what you will need once you start simulating turbulence --- which, happily, you are now prepared to do because you know what initial time-step to use! Such a strong stratification may require fairly high resolutions (0.25 m or finer?) in both horizontal and vertical resolution, because strong stratification reduces the size of turbulent eddies. The key length scale is the Ozmidov length scale, which is an emergent parameter in your simulatinos and thus can only be diagnosed after the results are in. The Ozmidov length is bigger if your forcing is stronger or your stratification is weaker. You'll need resolution around the Ozmidov length scale for the LES closures to work properly. Thinking all of this through, it implies that with strong stratification, you'd also need very strong forcing to get away with coarse resolution. It's not really possible to make a priori calculations to predict the Ozmidov scale a priori (unless you had run similar experiments before, and could use those as a guide...). Hopefully you can start running simulations soon and make some discoveries about what resolution you'll need!

What is the using Oceananigans.BuoyancyModels: BuoyancyField line doing? It was giving me an error and I saw that BuoyancyField was unused so I commented the line out and everything worked fine.

b = BuoyancyField(model) returns an object that represents buoyancy in your system. It works no matter what buoyancy model you use (linear or nonlinear equation of state, etc). Calling compute!(b) computes the buoyancy field. You can diagnose the buoyancy frequency of your stratification by calculating the vertical derivative, N² = ComputedField(∂z(b)). You can use this to plot the vertical profile of N², and also calculate the time-step restriction imposed by the buoyancy frequency peak in the middle of your pycnocline.

[glwagner]

Ah, referring to the figure I made above with the buoyancy frequency profile, it looks like you have a buoyancy time-scale of 1 / sqrt(0.006) \approx 12.9 seconds. I bet this was the problem (doh!) I think you need to limit your time-step to 1/2 or 1/4 of that time-scale for stable simulations. Please run some experiments and figure it out --- I'd love to know!

[sam12396]

The buoyancy frequency, time step stuff definitely makes sense and that is really interesting. I will definitely do a few tests on a matrix of max_Δt and dz options. I am pretty new to Julia but, from examples I have seen, it doesn't seem that hard to write the function wrapper so I will take your script from this thread and make some of those w*2 animations for different cases. Hopefully, we won't see any sensitivity to dz.

b = BuoyancyField(model) returns an object that represents buoyancy in your system. It works no matter what buoyancy model you use (linear or nonlinear equation of state, etc).

That's interesting, my machine was telling me that Oceananigans.BuoyancyModels didn't exist but maybe I am running a slightly older version. I have had some trouble when switching versions because I am not the most apt coder even in languages I work in more frequently, but I am learning the Julia basics ,slow and steady.

It's not really possible to make a priori calculations to predict the Ozmidov scale a priori (unless you had run similar experiments before, and could use those as a guide...)

Just so I am sure I understand, the idea here is to run a simulation at some coarse resolution, calculate the Ozmidov length scale, run at that new resolution, and rinse and repeat until there is less sensitivity? I'll be honest, I am definitely concerned about having to run at 0.25m resolution because I would like to keep my 100m x 100m x 85m domain extent and that would push me well over 1 million points. If we don't run at a resolution near the Ozmidov length scale, we are leaving more turbulence up to the closure correct? Sorry if this is an LES 101 question but if that were a class I would take it haha.

[glwagner]

No worries about being new to Julia! I hope exercises like the one I've proposed to investigate a stable time-step will be useful. The stable time-step I found (1 second) is 1/13th of the crucial time scale. I bet you can get away with a (slightly) larger time-step. As turbulence develops in your forced simulation, the time-step will probably be limited by CFL rather than the buoyancy frequency...

That's interesting, my machine was telling me that Oceananigans.BuoyancyModels didn't exist but maybe I am running a slightly older version.

In older versions the object BuoyancyField was exported by the module Oceananigans.Buoyancy.

Just so I am sure I understand, the idea here is to run a simulation at some coarse resolution, calculate the Ozmidov length scale, run at that new resolution, and rinse and repeat until there is less sensitivity?

I apologize if I made things more complicated than they need to be. The basic approach is to run two simulations at different resolutions (perhaps a factor of two different). If the results are the same, you might be reasonably convinced your results are "independent from grid resolution". A practical approach is to run simulations at the highest resolution you can afford, and then run a second simulation at 1/2 or 2/3 that resolution and compare the results. If you're running lots and lots of simulations (so you'd like to use a relatively coarse resolution), you might run a single test case at twice the resolution to test whether your results depend on resolution.

This is just an empirical test to answer the question: "are my results independent of grid resolution". If your results depend on grid resolution there is reason to suspect that they aren't realistic / reflective of real physics (instead, they reflect the behavior of some LES model).

The stuff I was saying about the Ozmidov scale is a theoretical explanation on why LES results might depend on resolution. Roughly speaking, LES closures seem to work well when motions with scales close to the grid scale are relatively isotropic and thus not strongly affected by stratification. This could help explain why finer resolution results are more likely to be unaffected by resolution, and why your results will be affected by grid resolution when the grid spacing exceeds a threshold. This is an empirical result --- you could probably write a paper about it! LES of stratified fluids is not that easy; there's much to be discovered. Suffice to say that you should simply be careful about interpreting your results, and make some effort to show that your results are not the arbitrary result of a model that happened to use 0.83 m resolution (and would be different if your grid resolution was 0.64 m instead), but rather that your results pertain to real fluid physics and ocean phenomena.

As for learning, a good place to start learning about LES might be to read Vreugdenhil and Taylor 2018, and probably also the references contained there. There might be some good textbooks too (my knowledge mostly comes from the literature...)