Large vertical velocities at early time steps in frontal instability simulation

[rafacmartins]

Hi all,

To give some context to my problem: I’ve been working on setting up a simulation of a submesoscale front in thermal-wind balance, which becomes unstable a few days later due to some some noise added to the buoyancy and zonal velocity fields.

The issue is that at the very beginning of the simulation (just a few time steps in), I observe relatively large vertical velocity values, and I don’t see any corresponding features in the horizontal Rossby number field, which I usually use to track the evolution of the simulation visually.

In the left panel of the video, I’ve plotted the horizontal Rossby number field. I can’t spot any obvious issues there; the frontal jet appears to stay in balance for a few days, until the instabilities begin to grow (see Figure 1 for reference). In the right panel, I show a vertical cross-section of vertical velocity along the front at the center of the domain. This panel shows relatively large vertical velocities in the very first time steps of the simulation.

combined_plot.mp4

Figure 1.

Below is the KE plot of the simulation (Figure 2) and a plot of the time variance of the RMS vertical velocities (Figure 3), where I can see a spike at the very beginning.

Figure 2

Figure 3.

Here is the code I used to run this simulation:

#Importing useful packages
using Oceananigans
using Oceananigans.Advection
using Oceananigans.Units
using NCDatasets
using CUDA
using SpecialFunctions
### Definying some parameters used to create the initial conditions
g = 9.81             # m/s²
f0 = 8.7e-5          # 1/s
Lx = 150000meters # east-west extent [m]
Ly = 150000meters # north-south extent [m]
Lz = 600meters    # depth [m]
U0 = 0.4             # m/s
delta0 = 50          # m
Nx = 300  
Ny = 200 
Nz = 30 
grid = RectilinearGrid(GPU(),size = (Nx, Ny, Nz),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = Lz * (cos.(LinRange(π/2, 0, Nz+1)) .- 1),
                       topology = (Periodic, Bounded, Bounded),
                        halo=(5,5,5))
uᵢ(x,y,z) = U0 * exp((z) / delta0) * exp(-(y0 - y)^2 / (2 * l^2)) + ϵb * randn() 
RHOstar(z) = rho_ml + drho * (1 - tanh((z+Lz - (Lz-2*delta0)) / delta1)) * (-z + c1 * erf(z / delta2) + c1 * erf((z + Lz) / delta2)) / (2 * Lz)
RHOprime(y,z) = c0 * exp(z / delta0) * (erf((y0 - y) / (sqrt(2) * l)) + erf(y0 / (sqrt(2) * l)))
RHO(y,z) = RHOprime(y,z) + RHOstar(z)
ϵb = 1e-8 #Noise amplitude
bᵢ(x,y,z) = -g * ((RHO(y,z)-rho0)/ rho0) + ϵb * randn()

#--------------- 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.
@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=100meters)#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)

horizontal_closure = HorizontalScalarDiffusivity()
vertical_closure = ScalarDiffusivity()
closure = (horizontal_closure, vertical_closure)

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

set!(model, b=bᵢ,u = uᵢ)
simulation = Simulation(model, Δt = 50seconds, stop_time = 30days)

#Preparing the output
b, = simulation.model.tracers
u, v, w = simulation.model.velocities
ζ = ∂x(v) - ∂y(u)
div = ∂x(u) + ∂y(v)

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),
)
outputs = merge(model.tracers,extra_outputs);

name = "original_front"
output_folder = name

# creating the folder if it not exists
mkpath(output_folder)

# Setting the NetCDFOutputWriter
simulation.output_writers[:fields] = NetCDFOutputWriter(
    model,
    outputs,
    overwrite_existing=true,
    filename=name*".nc",
    dir = output_folder,
    schedule=TimeInterval(1hour)
)

#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)

I'm wondering if anyone has experienced something similar or has any suggestions on what might be causing this issue. I'm particularly puzzled by the early spike in vertical velocities. Could this be related to something numerical happening at the beginning of the run? Any ideas on diagnostics I could check, or things I should look into to better understand where this behavior is coming from, would be very much appreciated!

[glwagner]

Very quick thought --- you are initializing the velocity and buoyancy field to be in thermal wind balance, correct? However, your initial conditions are expressed as continuous functions. This means that there is a small discrepancy between the analytical thermal wind balance being prescribed, and the "true" discrete thermal wind balance that is obeyed by the model. The small difference between discrete and analytical thermal wind balance could cause a small adjustment to occur at the beginning of the simulation.

Solving this may not be trivial, but is possible. In the first place, you can verify that the effect depends on resolution. The higher the resolution, the smaller the difference between the discrete thermal wind balance obeyed exactly by the model and the analytica thermal wind balance you assume to derive the velocity and buoyancy/density fields.

Another possibility is to set the geostrophic streamfunction first, and then differentiate it to find the initial conditions for velocity and buoyancy. For example,

ψ = Field{Face, Face, Center}(grid)
one = CenterField(grid)
set!(one, 1)
set!(ψ, geostrophic_streamfunction)
uᵢ = + ∂y(ψ)
bᵢ = one * f * ∂z(ψ)

The purpose of one being to correctly interpolate the vertical derivative of the streamfunction to the buoyancy location. I don't think this is discretely exact either because of the interpolation (and we do not invoke the stencil being used to compute the Coriolis force). It'd be interesting to develop a feature that actually computes the velocity field in exact discrete thermal wind balance...