Front simulation weird patterns at the surface (no forcing)

[iuryt]

Hi all,

I am having some weird patterns for buoyancy at the surface that blows up the model, even that I am not applying any forcing there.

The context

I am currently trying to setup a front simulation to couple with the NP model being developed in #2385 . The idea is to start with an unbalaced front that evolves in time until it reaches its balance. I expect that front will probably radiate some internal waves, but geostrophically adjust after some time.

The initial conditions

I am starting the simulation with a single buoyant front that adds to an initial buoyancy profile.

const g = 9.82 # gravity
const ρₒ = 1026 # reference density

# background density profile based on Argo data
@inline bg(z) = 0.25*tanh(0.0027*(-653.3-z))-6.8*z/1e5+1027.56

# decay function for fronts
@inline decay(z) = (tanh((z+500)/300)+1)/2

@inline front(x, y, z, cy) = tanh((y-cy)/12kilometers)
@inline D(x, y, z) = bg(z) + 0.8*decay(z)*front(x, y, z, 0)/4
@inline B(x, y, z) = -(g/ρₒ)*D(x, y, z)

# setting the initial conditions
set!(model; b=B)

The problem

After 1 timestep it starts to present weird patters at the surface. I am pretty sure I am messing up with something here, can you help me to figure this out?

Full code:

using Oceananigans
using Oceananigans.Units

# define the size and max depth of the simulation
const Ny = 100
const Nz = 48 # number of points in z
const H = 1000 # maximum depth

# create the grid of the model
grid = RectilinearGrid(GPU(),
    size=(Ny, Nz),
    halo=(3,3),
    y=(-(Ny/2)kilometers, (Ny/2)kilometers),
    z=(H * cos.(LinRange(π/2,0,Nz+1)) .- H)meters,
    topology=(Flat, Bounded, Bounded)
)

# define the turbulence closure of the model
horizontal_closure = ScalarDiffusivity(ν=1, κ=1)
vertical_closure = ScalarDiffusivity(ν=1e-5, κ=1e-5)

coriolis = FPlane(latitude=60)

#--------------- Instantiate Model

# create the model
model = NonhydrostaticModel(grid = grid,
                            advection = WENO5(),
                            timestepper = :RungeKutta3,
                            coriolis = coriolis,
                            closure=(horizontal_closure,vertical_closure),
                            tracers = (:b),
                            buoyancy = BuoyancyTracer())

#--------------- Initial Conditions

const g = 9.82 # gravity
const ρₒ = 1026 # reference density

# background density profile based on Argo data
@inline bg(z) = 0.25*tanh(0.0027*(-653.3-z))-6.8*z/1e5+1027.56

# decay function for fronts
@inline decay(z) = (tanh((z+500)/300)+1)/2

@inline front(x, y, z, cy) = tanh((y-cy)/12kilometers)
@inline D(x, y, z) = bg(z) + 0.8*decay(z)*front(x, y, z, 0)/4
@inline B(x, y, z) = -(g/ρₒ)*D(x, y, z)

# setting the initial conditions
set!(model; b=B)

#--------------- Simulation

simulation = Simulation(model, Δt = 10seconds, stop_time = 10minutes)

wizard = TimeStepWizard(cfl=1.0, max_change=1.1, max_Δt=6minutes)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(5))

# merge light and h to the outputs
outputs = merge(model.velocities, model.tracers) # make a NamedTuple with all outputs

# writing the output
simulation.output_writers[:fields] =
    NetCDFOutputWriter(model, outputs, filepath = "data/output.nc",
                     schedule=IterationInterval(1))

# run the simulation
run!(simulation)

[glwagner]

Should we convert to a discussion? There's a nice Q&A format there and you can give people credit for answering your q.

[glwagner]

Might be nice to include the plotting code!

[iuryt]

Should we convert to a discussion? There's a nice Q&A format there and you can give people credit for answering your q.

Of course!

[iuryt]

Can you convert this to a Discussion for me?