Problems with instability timescale in submesoscale front

[rafacmartins]

Hi all,
I've been trying to run a simulation of a submesoscale front that is in thermal-wind balance and that goes unstable when adding some white noise, but the time to develop the instabilities are much longer than we expected. To test if the numerical setup was the issue, I tried to replicate the results of the work by Boccaletti et al 2007 (Mixed Layer Instabilities and Restratification) using his initial conditions and the numerical setup of my simulations, and I found some similar results:

The pertubation KE evolution looks the same, and even the features developed on his simulations (done with MITgcm) appear quite similar, at least in terms of the time of development of the instabilities.

So, we started to look for some diagnoses of my initial conditions, in terms of some non-dimensional numbers, as Ri, Fr and Bu.

As a statement, my initial conditions looks like this:

For the calculation of the Ri I've been using the Oceanostics.jl package, and for the Fr and Bu I've used the following formulas:

$Fr = \frac{U}{NH}$
$Bu = (\frac{Ro}{Fr})^{2}$

where
$N = \sqrt{\frac{\partial b}{\partial z}}$

The problem is that when I calculate those properties I have some really weird results.
For the Ri, I have:

For the PV:

The code that I use to run the simulations is the follow:

#Importing useful packages
using Oceananigans
using Oceananigans.Advection
using Oceananigans.Units
using NCDatasets
using CUDA
using SpecialFunctions
using Oceanostics

## Open data from where I obtain the value of some parameters that I use to build the Initial condition
ds = NCDataset("Rho_front_8_adjusted_rho0_corrected.nc")
## Setting general constants
f0 = FPlane(latitude=37).f #s-¹
g=9.81 #m/s²
for (name, value) in ds.attrib
    @eval $(Symbol(name)) = $value
end
c0 = -sqrt(pi) * (f0 * U0 * l * rho0) / (sqrt(2) * g * delta0)
c1 = sqrt(pi) * delta2 / 2
y0 = (77e3/2)-(Lx/2)+(x0)
Ly = 77e3meters
Lx = 177e3meters
Lz = float(Lz)meters
z0 = Lz - 2*delta0
z1 = z0

#--------------------Creating the grid
Nx = 300 
Ny = 300
Nz = 50
grid = RectilinearGrid(GPU(),size = (Nx, Ny, Nz),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = Lz * (cos.(LinRange(π/2, 0, Nz+1)) .- 1), #stretched on vertical
                       topology = (Periodic, Bounded, Bounded),
                        halo=(5,5,5))

ϵb = 1e-8 #noise amplitude
uᵢ(x,y,z) = U0 * exp((z) / delta0) * exp(-(y0 - y)^2 / (2 * l^2)) #+ ϵb * randn() #adding white noise
uz(x,y,z) = U0/delta0 * exp((z)/delta0) * exp(-(y0 - y)^2 / (2 * l^2))
RHO_star(z) = rho_ml + drho*(1 - tanh((z+2*delta0)/delta1))*((-z+c1*erf(z/delta2) + c1*erf((z+Lz)/delta2))/(2*Lz))
RHO_prime(y,z) = c0 * exp((z)/delta0) * (erf((y0-y)/(sqrt(2)*l))+erf((y0)/(sqrt(2)*l)))
RHO(y,z) = RHO_prime(y,z) + RHO_star(z)
bᵢ(x,y,z) = -g * ((RHO(y,z)-rho0)/ rho0) #+ ϵb * randn() #turning into buoyancy and adding the white noise

#--------------- Sponges

const sponge_size = 2.5kilometers #width of the lateral sponge
const slope = 0.75kilometers #slope of the tanh, the smaller it gets the sharper gets the slope
x, y, z = nodes(grid, (Center(), Center(), Center())) #getting the grid dimensions 
const ymin = minimum(y)
const ymax = maximum(y)
#Creating a mask func to know where to apply the sponge layer. It is 0 outside of the sponge layer domain and 1 inside
@inline mask_func(x,y,z) = ((
     tanh((y-(ymax-sponge_size))/slope)
    *tanh((y-(ymin+sponge_size))/slope)
)+1)/2

bottom_mask = GaussianMask{:z}(center=-Lz, width=15meters)#Creating the mask for the bottom sponge layer
north_sponge = Relaxation(rate=1/30minutes, mask=mask_func, target=0) #in fact creating the sponge layer for the model
bottom_sponge = Relaxation(rate=1/30minutes, mask=bottom_mask, target=0) #in fact creating the sponge layer for the model

mom_forcings =  (north_sponge, bottom_sponge)
forcings = (; u = mom_forcings, v = mom_forcings)


#Setting the hydrostatic model within f-plane configuration, centered in 37°N
model = HydrostaticFreeSurfaceModel(; grid,
                                    coriolis = FPlane(latitude = 37),
                                    buoyancy = BuoyancyTracer(),
                                    tracers = (:b,:e),
                                    forcing = forcings, #appling sponge layers
                                    momentum_advection = WENO(),
                                    tracer_advection = WENO())

set!(model, b=bᵢ,u = uᵢ) #Appling the initial conditions to the model
simulation = Simulation(model, Δt = 50seconds, stop_time = 30days) #defining the time step and stop time

#------------- Preparing the output
b, = simulation.model.tracers
u, v, w = simulation.model.velocities
ζ = ∂x(v) - ∂y(u)
div = ∂x(u) + ∂y(v)
Ri = Oceanostics.RichardsonNumber(model, model.velocities..., b)
Ro = Oceanostics.RossbyNumber(model)
PV = Oceanostics.ErtelPotentialVorticity(model, model.velocities..., b, model.coriolis)

#interpolating to the center of the cells to make it easier to analyze
extra_outputs = (;
    u=@at((Center, Center, Center), u),
    v=@at((Center, Center, Center), v), 
    w=@at((Center, Center, Center), w), 
    ζ=@at((Center, Center, Center),ζ),
    div = @at((Center, Center, Center),div),
    Ri = @at((Center, Center, Center),Ri),
    Ro = @at((Center, Center, Center),Ro),
    PV = @at((Center, Center, Center),PV),
)
outputs = merge(model.tracers,extra_outputs);

name = "control_simulation_without_noise"
base_dir = "/scratch/home3/rmartins"
output_folder = joinpath(base_dir, name)
mkpath(output_folder)
println("Folder created: $output_folder")

# Setting the NetCDFOutputWriter
simulation.output_writers[:fields] = NetCDFWriter(
    model,
    outputs,
    overwrite_existing=true,
    filename=name*".nc",
    dir = output_folder,
    schedule=TimeInterval(1hour),
    global_attributes = Dict("U0" => U0 , "delta0" => delta0, "delta1" => delta1, "delta2" => delta2, "l" => l, "rho_ml" => rho_ml, "drho" => drho, "rho0" => rho0,"z0" => z0, "z1" => z1, "c0" => c0, "c1" => c1, "y0" => y0, "ϵb"=>ϵb, "slope"=>slope, "sponge_size"=>sponge_size)
)

#Creating a callback function to print the progress of the simulation and do the diagnosis while still running
using Printf

function print_progress(simulation)
    u, v, w = simulation.model.velocities

    # Print a progress message
    msg = @sprintf("i: %04d, t: %s, Δt: %s, umax = (%.1e, %.1e, %.1e) ms⁻¹, wall time: %s\n",
                   iteration(simulation),
                   prettytime(time(simulation)),
                   prettytime(simulation.Δt),
                   maximum(abs, u), maximum(abs, v), maximum(abs, w),
                   prettytime(simulation.run_wall_time))

    @info msg

    return nothing
end
simulation.callbacks[:progress] = Callback(print_progress, TimeInterval(4hours))

#Running the simulation
run!(simulation)

This version of the code does not include white noise (so instabilities do not develop), since I only want to analyze the initial conditions.

From #4485, I know that my initial conditions won't be in perfect thermal-wind balance due to the way that the equations are discretized in Oceananigans.
Still, I would like to hear your thoughts on why the perturbations take so long to develop. If we estimate the theoretical timescale of baroclinic instability using
$T = \frac{1}{0.31} \frac{\sqrt{Ri}}{|f|}$)
I find something close to $T = 5.1 hrs$ when considering $Ri = 0.25$, and about $T=10hr$ for $Ri = 0.95$. This is far from the timescale observed in my outputs, where instabilities only start to appear after about 10–12 days.

In addition, there are the results shown in the figures for Ri and PV, which may be more of a discussion related to Oceanostics.jl (and I can move the discussion there if needed).

Thank you in advance!

[glwagner]

Is it possible that there is a bug in your initial condition expressions? If this is possible, I'd suggest rewriting them directly in terms of b, and also expressing them in terms of the (bulk) non-dimensional numbers that you want to prescribe. This will help you diagnose whether there is a problem with the initial conditions, or with the diagnostics you are using to compute the point-wise non-dimensional numbers.