Baroclinic adjustment and the potential energy budget

In this example we spin up a pair of submesoscale fronts in a doubly-periodic channel, let them go baroclinically unstable, and close the volume-integrated potential energy budget against the kinetic energy one.

Before starting, make sure you have the required packages installed for this example, which can be done with

using Pkg
pkg"add Oceananigans, Oceanostics, CairoMakie, NCDatasets"

using Oceananigans
using Oceananigans.Units

Parameters

A 1 km by 2 km channel, 70 m deep, with a 30 m mixed layer over a stratified interior. The two fronts are 250 m wide and carry a cross-front buoyancy gradient M²:

Lx = 1kilometers
Ly = 2kilometers
H_target = 70.0      # [m] requested depth; the stretched grid below lands a little deeper
h  = 30.0            # [m] mixed-layer depth

f       = 1e-4       # [s⁻¹] Coriolis frequency
M²      = 3e-8       # [s⁻²] cross-front buoyancy gradient inside each front
N²_ml   = 9e-8       # [s⁻²] mixed-layer stratification (z > -h)
N²_int  = 1.8e-6     # [s⁻²] interior stratification (z < -h)
w_front = 250.0      # [m]   width of each front

250.0

Grid

Doubly periodic in the horizontal and bounded in the vertical. Periodicity in y is what the double front buys us, and it is what makes the transport terms of both budgets integrate to zero. The vertical coordinate is surface-intensified: a constant 2 m spacing over the top 32 m resolves the mixed layer the instability lives in, and stretches below it.

The cells come out roughly 16 m by 16 m by 2 to 8 m, within an order of magnitude of isotropic. That matters for the closure below: Smagorinsky builds its filter width from the cell volume, so it is only meaningful on a grid that is not wildly stretched in one direction.

z = ReferenceToStretchedDiscretization(extent = H_target,
                                       bias = :right, bias_edge = 0,   # fine spacing at the surface
                                       constant_spacing = 2,
                                       constant_spacing_extent = 32,
                                       stretching = PowerLawStretching(1.08))

grid = RectilinearGrid(size = (64, 128, length(z)),
                       x = (0, Lx), y = (-Ly/2, Ly/2), z = z,
                       topology = (Periodic, Periodic, Bounded))

H = grid.Lz   # [m] the depth the grid actually has

@info "Grid: $(size(grid)), depth $(round(H, digits=1)) m, Δx = $(round(minimum_xspacing(grid), digits=1)) m, Δz from $(round(minimum_zspacing(grid), digits=2)) m to $(round(maximum(zspacings(grid, Center())), digits=2)) m"

[ Info: Grid: (64, 128, 28), depth 72.9 m, Δx = 15.6 m, Δz from 2.0 m to 5.66 m

Derived quantities

The mixed layer sets the scale of the instability: its deformation radius Ld = N h / f, the fastest-growing wavelength λ ≈ 3.9 Ld, and the growth rate σ ≈ 0.31 M²/N. The balanced Richardson number Ri = N²/α² is 1 in the mixed layer, which is what lets the resolved strain reach the stratification, and so what lets the closure below do anything at all.

α     = M² / f              # [s⁻¹]   thermal-wind shear inside a front
N_ml  = √N²_ml              # [s⁻¹]
N_int = √N²_int             # [s⁻¹]
Ld    = N_ml * h / f        # [m]     mixed-layer deformation radius
λ     = 3.9 * Ld            # [m]     fastest-growing wavelength
Ri_ml = N²_ml / α^2         # []      balanced Richardson number in the mixed layer
Δb    = M² * w_front        # [m s⁻²] buoyancy jump across each front
Ū     = α * H               # [m s⁻¹] thermal-wind velocity scale
σ     = 0.31 * M² / N_ml    # [s⁻¹]   growth rate of the mixed-layer mode

@info "Ld = $(round(Ld, digits=1)) m, λ = $(round(λ, digits=1)) m in a $(round(Int, Ly)) m channel, Ri = $(round(Ri_ml, digits=2))"
@info "Δb = $(round(Δb, sigdigits=3)) m s⁻², Ū = $(round(100Ū, digits=2)) cm/s, growth time 1/σ = $(prettytime(1/σ))"

[ Info: Ld = 90.0 m, λ = 351.0 m in a 2000 m channel, Ri = 1.0
[ Info: Δb = 7.5e-6 m s⁻², Ū = 2.19 cm/s, growth time 1/σ = 8.961 hours

Closure

We use a Smagorinsky closure with a constant coefficient to model turbulence stresses. We use a high coefficient to make sure the instability is well-resolved in this very coarse example:

using Oceananigans.TurbulenceClosures.Smagorinskys: Smagorinsky
closure = Smagorinsky(coefficient=0.3)

Smagorinsky closure with
├── coefficient = 0.3
└── Pr = 1.0

Model

The advection scheme is Centered, which adds no dissipation of its own, so every sink in the budgets below is one we write down and compute. We set up an f-plane Coriolis:

model = NonhydrostaticModel(grid;
                            coriolis = FPlane(; f),
                            buoyancy = BuoyancyTracer(),
                            tracers = :b,
                            timestepper = :RungeKutta3,
                            advection = Centered(order=4),
                            closure = closure)

NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 64×128×28 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=4)
├── tracers: b
├── closure: Smagorinsky with coefficient = Float64, Pr=(b = 1.0,)
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
└── coriolis: FPlane{Float64}(f=0.0001)

Initial condition

The buoyancy is the two-layer stratification plus two ramps of width w_front, one rising at y = -Ly/4 and one falling at y = +Ly/4. Their difference returns to zero at the edges of the domain, so b matches across the periodic boundary:

ramp(y, w) = min(max(0, y/w + 1/2), 1)

y₁ = -Ly/4
y₂ = +Ly/4

double_ramp(y, w) = ramp(y - y₁, w) - ramp(y - y₂, w)

# the stratification integrated from the surface down, continuous across the mixed-layer base
b_strat(z) = ifelse(z ≥ -h, N²_ml * z, -N²_ml * h + N²_int * (z + h))

b_strat (generic function with 1 method)

The velocity starts in thermal-wind balance with those fronts, referenced to the bottom, and a little noise on the buoyancy seeds the instability:

using Random
Random.seed!(8675309)

ϵb = 1e-2 * Δb    # noise amplitude, a percent of the front's buoyancy jump

bᵢ(x, y, z) = b_strat(z) + Δb * double_ramp(y, w_front) + ϵb * randn()

ramp_prime(y, w) = (-w/2 < y < w/2) ? 1/w : zero(y)
double_ramp_prime(y, w) = ramp_prime(y - y₁, w) - ramp_prime(y - y₂, w)

uᵢ(x, y, z) = -(Δb * double_ramp_prime(y, w_front) / f) * (z + H)

set!(model, u=uᵢ, b=bᵢ)

Simulation

max_Δt = min(minimum_xspacing(grid) / Ū, 1 / N_int)   # Smagorinsky sets ν from the flow, so no fixed diffusive bound

simulation = Simulation(model, Δt = max_Δt, stop_time = 8days)
conjure_time_step_wizard!(simulation, IterationInterval(5), cfl = 0.7, max_Δt = max_Δt)

using Oceanostics
add_callback!(simulation, ProgressMessengers.TimedMessenger(), IterationInterval(100))

The potential energy budget

Given eₚ = -bz, we multiply the budget equation for b by -z. Pulling z inside each derivative splits terms into a transport and a conversion, and over this domain every transport integrates away, leaving just the two conversions:

\[\frac{d}{dt}\int e_p\, dV = -\int wb\, dV + \int \Phi\, dV ,\]

with $wb$ (PotentialToKineticEnergyConversion) the exchange with the kinetic energy and $\Phi = \kappa\,\partial b/\partial z$ (PotentialEnergyDiffusiveVerticalBuoyancyFlux) the work diffusion does against gravity.

eₚ = PotentialEnergy(model)
wb = PotentialToKineticEnergyConversion(model)
Φ  = PotentialEnergyDiffusiveVerticalBuoyancyFlux(model)

PotentialEnergyDiffusiveVerticalBuoyancyFlux KernelFunctionOperation at (Center, Center, Center)
├── grid: 64×128×28 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: diffusive_buoyancy_flux_ccc (generic function with 1 method)
└── arguments: ("Oceananigans.TurbulenceClosures.Smagorinskys.Smagorinsky", "NamedTuple", "Val", "Field", "Clock", "NamedTuple", "BuoyancyForce")
└── computes: diffusive vertical buoyancy flux  Φ = κ ∂b/∂z = -q₃

The kinetic energy budget

The other side of the exchange, written with eₖ = ½uᵢuᵢ to match eₚ. wb is the same term in both and carries opposite signs, so it cancels from their sum:

\[\frac{d}{dt}\int e_k\, dV = \int wb\, dV - \int \varepsilon_k\, dV .\]

eₖ = KineticEnergy(model)
εₖ = KineticEnergyDissipationRate(model)

∫eₚ = ∫dV(eₚ)
∫Φ  = ∫dV(Φ)
∫eₖ = ∫dV(eₖ)
∫wb = ∫dV(wb)
∫εₖ = ∫dV(εₖ)

1×1×1 Field{Nothing, Nothing, Nothing} reduced over dims = (1, 2, 3) on RectilinearGrid on CPU
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (1, 1, 1)
├── grid: 64×128×28 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── operand: Integral of BinaryOperation at (Center, Center, Center) over dims (1, 2, 3)
├── status: time=0.0
└── data: 1×1×1 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, 1:1) with eltype Float64 with indices 1:1×1:1×1:1
    └── max=0.0202361, min=0.0202361, mean=0.0202361

For the movie we keep the surface vorticity, buoyancy, and the closure's own eddy viscosity, which shows where Smagorinsky is actually acting:

u, v, w = model.velocities
b = model.tracers.b
ζ = ∂x(v) - ∂y(u)
νₑ = viscosity(model)

64×128×28 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 64×128×28 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 70×134×34 OffsetArray(::Array{Float64, 3}, -2:67, -2:131, -2:31) with eltype Float64 with indices -2:67×-2:131×-2:31
    └── max=0.00557548, min=0.0, mean=0.000655031

Output

A snapshot writer for the surface maps and a budget writer for the volume integrals, the latter on ConsecutiveIterations, which takes a second sample one model step after each output time so we can finite-difference d/dt across that step.

using NCDatasets
filename = joinpath(@__DIR__, "baroclinic_adjustment")

simulation.output_writers[:fields] =
    NetCDFWriter(model, (; ζ, b, νₑ),
                 filename = filename,
                 schedule = TimeInterval(3hours),
                 indices = (:, :, grid.Nz),
                 overwrite_existing = true)

simulation.output_writers[:budget] =
    NetCDFWriter(model, (; ∫eₚ, ∫Φ, ∫eₖ, ∫wb, ∫εₖ),
                 filename = filename * "_budget",
                 schedule = ConsecutiveIterations(TimeInterval(3hours)),
                 overwrite_existing = true)

NetCDFWriter scheduled on ConsecutiveIterations(TimeInterval(3 hours), 1):
├── filepath: baroclinic_adjustment_budget.nc
├── dimensions: time(0), y_afa(128), x_faa(64), x_caa(64), y_aca(128), z_aaf(29), z_aac(28)
├── 5 outputs: (∫Φ, ∫εₖ, ∫wb, ∫eₚ, ∫eₖ)
├── array_type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: 35.1 KiB

Run the simulation

run!(simulation)

[ Info: Initializing simulation...
┌ Info: iter =      0,  [000.00%] time = 0 seconds,  Δt = 8.450 minutes,  walltime = 1.166 minutes,  walltime / timestep = 0 seconds
└       |u⃗|ₘₐₓ = [2.16e-02,  0.00e+00,  0.00e+00] m/s,  advective CFL = 0.7,  diffusive CFL = 0.71,  νₘₐₓ = 0.0056 m²/s
[ Info:     ... simulation initialization complete (14.271 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.544 seconds).
┌ Info: iter =    100,  [006.25%] time = 12 hours,  Δt = 6.232 minutes,  walltime = 1.771 minutes,  walltime / timestep = 363.377 ms
└       |u⃗|ₘₐₓ = [2.45e-02,  5.37e-03,  5.28e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.43,  νₘₐₓ = 0.0046 m²/s
┌ Info: iter =    200,  [011.70%] time = 22.464 hours,  Δt = 6.822 minutes,  walltime = 2.142 minutes,  walltime / timestep = 222.217 ms
└       |u⃗|ₘₐₓ = [2.45e-02,  2.07e-03,  1.88e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.42,  νₘₐₓ = 0.0041 m²/s
┌ Info: iter =    300,  [017.84%] time = 1.427 days,  Δt = 7.578 minutes,  walltime = 2.505 minutes,  walltime / timestep = 218.151 ms
└       |u⃗|ₘₐₓ = [2.26e-02,  1.45e-03,  2.03e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.41,  νₘₐₓ = 0.0036 m²/s
┌ Info: iter =    400,  [024.36%] time = 1.949 days,  Δt = 7.568 minutes,  walltime = 2.870 minutes,  walltime / timestep = 219.017 ms
└       |u⃗|ₘₐₓ = [2.24e-02,  2.55e-03,  3.95e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.39,  νₘₐₓ = 0.0035 m²/s
┌ Info: iter =    500,  [030.51%] time = 2.441 days,  Δt = 6.717 minutes,  walltime = 3.233 minutes,  walltime / timestep = 217.772 ms
└       |u⃗|ₘₐₓ = [2.28e-02,  6.27e-03,  8.33e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.36,  νₘₐₓ = 0.0036 m²/s
┌ Info: iter =    600,  [036.18%] time = 2.894 days,  Δt = 6.958 minutes,  walltime = 3.597 minutes,  walltime / timestep = 218.100 ms
└       |u⃗|ₘₐₓ = [2.26e-02,  8.12e-03,  7.60e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.33,  νₘₐₓ = 0.0031 m²/s
┌ Info: iter =    700,  [042.25%] time = 3.380 days,  Δt = 7.265 minutes,  walltime = 3.934 minutes,  walltime / timestep = 202.097 ms
└       |u⃗|ₘₐₓ = [2.23e-02,  7.98e-03,  4.67e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.32,  νₘₐₓ = 0.0029 m²/s
┌ Info: iter =    800,  [048.64%] time = 3.891 days,  Δt = 7.725 minutes,  walltime = 4.340 minutes,  walltime / timestep = 243.618 ms
└       |u⃗|ₘₐₓ = [2.17e-02,  7.73e-03,  5.60e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.33,  νₘₐₓ = 0.0028 m²/s
┌ Info: iter =    900,  [055.01%] time = 4.401 days,  Δt = 7.476 minutes,  walltime = 4.725 minutes,  walltime / timestep = 231.022 ms
└       |u⃗|ₘₐₓ = [2.04e-02,  9.94e-03,  5.78e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.32,  νₘₐₓ = 0.0028 m²/s
┌ Info: iter =   1000,  [061.40%] time = 4.912 days,  Δt = 7.576 minutes,  walltime = 5.108 minutes,  walltime / timestep = 229.834 ms
└       |u⃗|ₘₐₓ = [1.94e-02,  1.10e-02,  1.01e-03] m/s,  advective CFL = 0.7,  diffusive CFL = 0.32,  νₘₐₓ = 0.0028 m²/s
┌ Info: iter =   1100,  [067.84%] time = 5.427 days,  Δt = 7.568 minutes,  walltime = 5.492 minutes,  walltime / timestep = 230.626 ms
└       |u⃗|ₘₐₓ = [1.81e-02,  1.11e-02,  1.76e-03] m/s,  advective CFL = 0.7,  diffusive CFL = 0.38,  νₘₐₓ = 0.0033 m²/s
┌ Info: iter =   1200,  [074.42%] time = 5.954 days,  Δt = 8.441 minutes,  walltime = 5.876 minutes,  walltime / timestep = 230.185 ms
└       |u⃗|ₘₐₓ = [1.64e-02,  8.95e-03,  1.39e-03] m/s,  advective CFL = 0.7,  diffusive CFL = 0.43,  νₘₐₓ = 0.0034 m²/s
┌ Info: iter =   1300,  [081.39%] time = 6.511 days,  Δt = 7.920 minutes,  walltime = 6.319 minutes,  walltime / timestep = 265.703 ms
└       |u⃗|ₘₐₓ = [1.44e-02,  9.23e-03,  8.71e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.3,  νₘₐₓ = 0.0025 m²/s
┌ Info: iter =   1400,  [088.21%] time = 7.057 days,  Δt = 8.188 minutes,  walltime = 6.703 minutes,  walltime / timestep = 230.426 ms
└       |u⃗|ₘₐₓ = [1.35e-02,  9.97e-03,  5.04e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.36,  νₘₐₓ = 0.0029 m²/s
┌ Info: iter =   1500,  [095.14%] time = 7.611 days,  Δt = 8.420 minutes,  walltime = 7.087 minutes,  walltime / timestep = 230.861 ms
└       |u⃗|ₘₐₓ = [1.37e-02,  9.25e-03,  7.05e-04] m/s,  advective CFL = 0.7,  diffusive CFL = 0.36,  νₘₐₓ = 0.0028 m²/s
[ Info: Simulation is stopping after running for 6.225 minutes.
[ Info: Simulation time 8 days equals or exceeds stop time 8 days.

Process the results

using CairoMakie

ds = NCDataset(simulation.output_writers[:fields].filepath)
times = ds["time"][:]
x_faa = ds["x_faa"][:]; y_afa = ds["y_afa"][:]
x_caa = ds["x_caa"][:]; y_aca = ds["y_aca"][:]
ζ_arr = ds["ζ"][:, :, 1, :]
b_arr = ds["b"][:, :, 1, :]
ν_arr = ds["νₑ"][:, :, 1, :]
close(ds)

ζlim = maximum(abs, ζ_arr)
νlim = maximum(ν_arr)

0.0055754795f0

The budget scalars come in consecutive-iteration pairs (2k-1, 2k); a one-step finite difference inside each pair gives the tendencies, and every source term is averaged over the same pair.

ds_b = NCDataset(simulation.output_writers[:budget].filepath)
t_bud    = ds_b["time"][:]
eₚ_bud = ds_b["∫eₚ"][:]
Φ_bud  = ds_b["∫Φ"][:]
eₖ_bud = ds_b["∫eₖ"][:]
wb_bud = ds_b["∫wb"][:]
εₖ_bud = ds_b["∫εₖ"][:]
close(ds_b)

idx1 = 1:2:length(t_bud) - 1
idx2 = 2:2:length(t_bud)

Δt_pair = t_bud[idx2] .- t_bud[idx1]
t_pair  = @. 0.5 * (t_bud[idx1] + t_bud[idx2])

deₚdt = (eₚ_bud[idx2] .- eₚ_bud[idx1]) ./ Δt_pair
deₖdt = (eₖ_bud[idx2] .- eₖ_bud[idx1]) ./ Δt_pair

pair_mean(x) = @. 0.5 * (x[idx1] + x[idx2])

Φ_pair  = pair_mean(Φ_bud)
wb_pair = pair_mean(wb_bud)
εₖ_pair = pair_mean(εₖ_bud);

Both budgets in sum-to-zero form: every curve is plotted with the sign it carries here, so each panel below adds up to its residual.

eₚ_resid = @. -deₚdt - wb_pair + Φ_pair
eₖ_resid = @. -deₖdt + wb_pair - εₖ_pair

Plotting

set_theme!(Theme(fontsize = 18))
fig = Figure(size = (1500, 950))

n = Observable(1)

panel_kwargs = (xlabel = "x [m]", ylabel = "y [m]", aspect = DataAspect(), height = 240)
axζ = Axis(fig[2, 1]; title = "vertical vorticity, ζ",  panel_kwargs...)
axb = Axis(fig[2, 3]; title = "surface buoyancy, b",    panel_kwargs...)
axν = Axis(fig[2, 5]; title = "eddy viscosity, νₑ",     panel_kwargs...)

ζₙ = @lift ζ_arr[:, :, $n]
bₙ = @lift b_arr[:, :, $n]
νₙ = @lift ν_arr[:, :, $n]

hmζ = heatmap!(axζ, x_faa, y_afa, ζₙ; colormap = :balance, colorrange = (-ζlim, ζlim))
Colorbar(fig[2, 2], hmζ)

hmb = heatmap!(axb, x_caa, y_aca, bₙ; colormap = :thermal)
Colorbar(fig[2, 4], hmb)

hmν = heatmap!(axν, x_caa, y_aca, νₙ; colormap = :tempo, colorrange = (0, νlim))
Colorbar(fig[2, 6], hmν)

budget_kwargs = (xlabel = "time [days]", ylabel = "[m⁵ s⁻³]")

ax_p = Axis(fig[3, 1:6]; title = "Volume-integrated potential energy budget", budget_kwargs...)
lines!(ax_p, t_pair ./ day, -deₚdt,    label = "-d(∫eₚ)/dt")
lines!(ax_p, t_pair ./ day, -wb_pair,  label = "-∫wb dV")
lines!(ax_p, t_pair ./ day,  Φ_pair,   label = "∫Φ dV")
lines!(ax_p, t_pair ./ day,  eₚ_resid, label = "residual", color = :black, linestyle = :dash)
axislegend(ax_p; position = :rt, labelsize = 10, nbanks = 2)

ax_k = Axis(fig[4, 1:6]; title = "Volume-integrated kinetic energy budget", budget_kwargs...)
lines!(ax_k, t_pair ./ day, -deₖdt,    label = "-d(∫eₖ)/dt")
lines!(ax_k, t_pair ./ day,  wb_pair,  label = "∫wb dV")
lines!(ax_k, t_pair ./ day, -εₖ_pair,  label = "-∫εₖ dV")
lines!(ax_k, t_pair ./ day,  eₖ_resid, label = "residual", color = :black, linestyle = :dash)
axislegend(ax_k; position = :rt, labelsize = 10, nbanks = 2)

vlines!(ax_p, @lift(times[$n] / day), color = :black, linestyle = :dot)
vlines!(ax_k, @lift(times[$n] / day), color = :black, linestyle = :dot)

title = @lift "Baroclinic adjustment, t = " * prettytime(times[$n])
fig[1, 1:6] = Label(fig, title, fontsize = 22, tellwidth = false)

@info "Animating..."
record(fig, "baroclinic_adjustment.mp4", 1:length(times), framerate = 8) do i
    n[] = i
end

[ Info: Animating...

This page was generated using Literate.jl.