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 Pkgpkg"add Oceananigans, Oceanostics, CairoMakie, NCDatasets"using Oceananigansusing Oceananigans.UnitsParameters
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 = 1kilometersLy = 2kilometersH_target = 70.0 # [m] requested depth; the stretched grid below lands a little deeperh = 30.0 # [m] mixed-layer depthf = 1e-4 # [s⁻¹] Coriolis frequencyM² = 3e-8 # [s⁻²] cross-front buoyancy gradient inside each frontN²_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 front250.0Grid
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 frontN_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 wavelengthRi_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: Smagorinskyclosure = Smagorinsky(coefficient=0.3)Smagorinsky closure with
├── coefficient = 0.3
└── Pr = 1.0Model
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/4y₂ = +Ly/4double_ramp(y, w) = ramp(y - y₁, w) - ramp(y - y₂, w)# the stratification integrated from the surface down, continuous across the mixed-layer baseb_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 RandomRandom.seed!(8675309)ϵb = 1e-2 * Δb # noise amplitude, a percent of the front's buoyancy jumpbᵢ(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 boundsimulation = Simulation(model, Δt = max_Δt, stop_time = 8days)conjure_time_step_wizard!(simulation, IterationInterval(5), cfl = 0.7, max_Δt = max_Δt)using Oceanosticsadd_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.0202361For 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.velocitiesb = 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.000655031Output
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 NCDatasetsfilename = 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: (file not yet created)Run the simulation
run!(simulation)[ Info: Initializing simulation...
┌ Info: iter = 0, [000.00%] time = 0 seconds, Δt = 8.450 minutes, walltime = 43.535 seconds, 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 (56.328 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (5.472 seconds).
┌ Info: iter = 100, [006.25%] time = 12 hours, Δt = 6.232 minutes, walltime = 2.148 minutes, walltime / timestep = 853.385 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.643 minutes, walltime / timestep = 296.619 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 = 3.137 minutes, walltime / timestep = 296.874 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 = 3.635 minutes, walltime / timestep = 298.907 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 = 4.136 minutes, walltime / timestep = 300.237 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 = 4.631 minutes, walltime / timestep = 297.372 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 = 5.085 minutes, walltime / timestep = 272.321 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 = 5.606 minutes, walltime / timestep = 312.590 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 = 6.086 minutes, walltime / timestep = 287.718 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 = 6.570 minutes, walltime / timestep = 290.272 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 = 7.057 minutes, walltime / timestep = 292.514 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 = 7.556 minutes, walltime / timestep = 299.245 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 = 8.123 minutes, walltime / timestep = 340.199 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 = 8.630 minutes, walltime / timestep = 304.090 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 = 9.133 minutes, walltime / timestep = 302.138 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 8.812 minutes.
[ Info: Simulation time 8 days equals or exceeds stop time 8 days.
Process the results
using CairoMakieds = 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.0055754795f0The 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) - 1idx2 = 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_pairdeₖdt = (eₖ_bud[idx2] .- eₖ_bud[idx1]) ./ Δt_pairpair_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 + Φ_paireₖ_resid = @. -deₖdt + wb_pair - εₖ_pairPlotting
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[] = iend[ Info: Animating...
This page was generated using Literate.jl.