Numerical instability in the Eady example

[zhihua-zheng]

Hi,

I want to reproduce Figure 1 in this recent paper introducing the OceanBioME. By the way, Oceananigans + biogeochemistry is so cool! Thank you all for the hard work!

I followed the script here to set up a simulation for Eady turbulence. However, I noticed some weird behaviors at the beginning of the simulation (see animation 1). The near-surface slumping of isopycnal immediately after initialization does not seem physical. This (numerical instability) is potentially related to the generous parameters in the time step wizard.

I have tried a smaller CFL number and more limiting parameters in the time step wizard:

wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_Δt = 2minutes, min_change = 0.5, max_change = 1.5)

This gives a smaller vertical velocity and there is no near-surface slumping (see animation 2), which are expected (?) on the first day after initialization.

This behavior is easy to overlook, because using the original code, I can still see the submesoscale eddy develop later in the simulation. The model does not report any warnings or errors. A related question on that is:

Does the initial (artificial) near-surface stratification affect the dynamics of mixed-layer instability later on?

animation1.mp4

animation2.mp4

[johnryantaylor]

Hi @zhihua-zheng, thanks for spotting this! We had some issues with the initial time step for this example before (the model was blowing up for other parameter choices), but I hadn't seen this. I agree that we should use a timestep that doesn't cause any numerical problems. A CFL number of 0.2 is fairly restrictive. Since the problem is only arising at the srtart, I wonder if we could get around it by reducing the maximum timestep and leaving the CFL number at 0.75 (or perhaps 0.5). Have you tried that?

[zhihua-zheng]

Yes, indeed. In the case above, I have reduced both the CFL number and max_Δt. I think it is possible to just reduce max_Δt and keep the CFL number around 0.75, with some trial and error. Or maybe it is even better to start with a small max_Δt and change it to a larger value after some transition time (CliMA/Oceananigans.jl#1361 (reply in thread)).

[jagoosw]

Hi @zhihua-zheng,

Thanks for pointing this out!

Is the only difference between the two plots the CFL number? It seems odd how different the results are.

[zhihua-zheng]

Yes, the only difference is the parameters used in the time step wizard, including cfl, max_Δt, etc, as described in the post. I have seen this kind of behavior in other setups as well. But how different the submesoscale fields are in the end is another question.