How to add a "background" field within the Hydrostatic model configuration

[rafacmartins]

Hi all,
I'm currently running some simulations of a submesoscale front in thermal wind balance that goes unstable when I add noise to it. My current configuration is zonally homogeneous, so the instabilities develop everywhere simultaneously. I was thinking of trying to change this. To control where the instabilities may develop I thought about using a mesoscale strain field, that is positive from 0 to Lx/2 and is negative from Lx/2 to Lx. This, at least in my conception, would induce frontogenesis in the first half and frontolysis in the second half. The problem is that I don't have an idea on how to add that in my current setup, specially considering that I'm using the HydrostaticFreeSurfaceModel, which does not allow the usage of BackgroundField.

Bellow is a minimum example of my current setup.

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

#Setting useful constants 
f0 = FPlane(latitude=37).f #s-¹
g = 9.81 #m/s²

Ly = 77e3meters
Lx = 177e3meters
Lz = 200meters

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

#Setting the values for the parameters used on the initial condition
delta0 = 25
rhoml = 1023
deltarho = 1.3
rho0 = 1025
l = 2.5e3
U0 = 1.0
y0 = Ly/2
z0 = -25
mld = 50
deltasm = 15
H = Lz

#Setting the initial condition in thermal-wind balance
ϵb = 1e-8
uᵢ(x,y,z) = U0/2  * (1+tanh((z-z0)/delta0))*exp(-1/2 * ((y-y0)^2/(l^2)))  + ϵb * randn()
F(z) = 1/ cosh((z-z0) / delta0)^2
c0 = (f0*U0*l*rho0)/(2*g*delta0) * ((sqrt(2*pi))/2)
B(y) = erf((y-y0)/(sqrt(2)*l))+erf(y0/(sqrt(2)*l))
RHO_prime(y,z) = c0*F(z)*B(y)
RHO_star(z) = (rhoml + 
    (deltarho/2)*(1-tanh((z+mld)/deltasm))*
    (1-(z+mld-deltasm)/(H-mld)))
RHO(y,z) = RHO_star(z) + RHO_prime(y,z)
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. 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,
                                    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

I did not include the section on output definitions, since it is not directly relevant to the question I’m asking.

Any ideas on how to apply this "mesoscale strain" field would be appreciate! Thanks in advance!

[glwagner]

Unfortunately the hydrostatic model does not support background fields right now!

On the other hand, you are not using a vertical mixing closure and you have a rectilinear grid. Therefore, you should be able to use the NonhydrostaticModel with very similar results as the HydrostaticFreeSurfaceModel.