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. A single (Eady) front would need a uniform background buoyancy gradient, which makes the domain's potential energy infinite. 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:
│ ├── momentum: Centered(order=4)
│ └── b: 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\, \mathrm{d}V = -\int wb\, \mathrm{d}V + \int \Phi\, \mathrm{d}V ,\]
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\, \mathrm{d}V = \int wb\, \mathrm{d}V - \int \varepsilon_k\, \mathrm{d}V .\]
eₖ = KineticEnergy(model)εₖ = KineticEnergyDissipationRate(model)∫eₚ = ∫dV(eₚ)∫Φ = ∫dV(Φ)∫eₖ = ∫dV(eₖ)∫wb = ∫dV(wb)∫εₖ = ∫dV(εₖ)∂ₜ∫eₚ = TimeDerivative(∫eₚ)∂ₜ∫eₖ = TimeDerivative(∫eₖ)TimeDerivative of 1×1×1 Field{Nothing, Nothing, Nothing} reduced over dims = (1, 2, 3) on RectilinearGrid on CPUFor the movie we keep the surface vorticity, buoyancy, the closure's own eddy viscosity (which shows where Smagorinsky is actually acting), and the mixed layer depth. Its criterion threshold is set to the initial buoyancy jump across the mixed layer, N²_ml * h, since the default threshold assumes realistic ocean stratification, well above what this coarse example reaches:
u, v, w = model.velocitiesb = model.tracers.bζ = ∂x(v) - ∂y(u)νₑ = viscosity(model)MLD = MixedLayerDepth(grid, model.buoyancy.formulation, (; b); criterion = BuoyancyAnomalyCriterion(threshold = -N²_ml * h))MixedLayerDepth KernelFunctionOperation at (Center, Center, ⋅)
├── grid: 64×128×28 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: Oceanostics.FlowDiagnostics.MixedLayerDepthKernel{BuoyancyAnomalyCriterion{Float64}}
└── arguments: ("BuoyancyTracer", "NamedTuple")
└── computes: mixed layer depth (shallowest depth where the criterion is met)Output
A snapshot writer for the surface maps and a budget writer for the volume integrals, both on TimeInterval(3hours).
using NCDatasetsfilename = joinpath(@__DIR__, "baroclinic_adjustment")simulation.output_writers[:fields] = NetCDFWriter(model, (; ζ, b, νₑ, MLD), filename = filename, schedule = TimeInterval(3hours), indices = (:, :, grid.Nz), overwrite_files = true)simulation.output_writers[:budget] = NetCDFWriter(model, (; ∂ₜ∫eₚ, ∫Φ, ∂ₜ∫eₖ, ∫wb, ∫εₖ), filename = filename * "_budget", schedule = TimeInterval(3hours), overwrite_files = true)NetCDFWriter scheduled on TimeInterval(3 hours):
├── filepath: baroclinic_adjustment_budget.nc
├── dimensions: time(0), x_faa(64), x_caa(64), y_afa(128), y_aca(128), z_aaf(29), z_aac(28)
├── 5 outputs: (∂ₜ∫eₚ, ∫Φ, ∂ₜ∫eₖ, ∫wb, ∫εₖ)
├── 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 = 39.187 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 (24.795 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (635.818 ms).
┌ Info: iter = 100, [006.25%] time = 12 hours, Δt = 6.235 minutes, walltime = 57.342 seconds, walltime / timestep = 181.407 ms
└ |u⃗|ₘₐₓ = [2.44e-02, 5.37e-03, 5.35e-04] m/s, advective CFL = 0.7, diffusive CFL = 0.43, νₘₐₓ = 0.0046 m²/s
┌ Info: iter = 200, [011.70%] time = 22.465 hours, Δt = 6.827 minutes, walltime = 1.151 minutes, walltime / timestep = 116.990 ms
└ |u⃗|ₘₐₓ = [2.45e-02, 2.10e-03, 1.75e-04] m/s, advective CFL = 0.7, diffusive CFL = 0.42, νₘₐₓ = 0.0041 m²/s
┌ Info: iter = 300, [017.85%] time = 1.428 days, Δt = 7.606 minutes, walltime = 1.359 minutes, walltime / timestep = 124.995 ms
└ |u⃗|ₘₐₓ = [2.26e-02, 1.47e-03, 1.93e-04] m/s, advective CFL = 0.7, diffusive CFL = 0.41, νₘₐₓ = 0.0036 m²/s
┌ Info: iter = 400, [024.35%] time = 1.948 days, Δt = 7.454 minutes, walltime = 1.554 minutes, walltime / timestep = 116.967 ms
└ |u⃗|ₘₐₓ = [2.24e-02, 2.95e-03, 3.39e-04] m/s, advective CFL = 0.7, diffusive CFL = 0.39, νₘₐₓ = 0.0035 m²/s
┌ Info: iter = 500, [030.40%] time = 2.432 days, Δt = 6.880 minutes, walltime = 1.766 minutes, walltime / timestep = 127.032 ms
└ |u⃗|ₘₐₓ = [2.28e-02, 5.99e-03, 1.29e-03] m/s, advective CFL = 0.7, diffusive CFL = 0.4, νₘₐₓ = 0.0039 m²/s
┌ Info: iter = 600, [036.24%] time = 2.899 days, Δt = 7.041 minutes, walltime = 1.984 minutes, walltime / timestep = 131.066 ms
└ |u⃗|ₘₐₓ = [2.06e-02, 7.86e-03, 4.75e-04] m/s, advective CFL = 0.7, diffusive CFL = 0.32, νₘₐₓ = 0.003 m²/s
┌ Info: iter = 700, [042.44%] time = 3.396 days, Δt = 7.366 minutes, walltime = 2.189 minutes, walltime / timestep = 123.184 ms
└ |u⃗|ₘₐₓ = [2.12e-02, 6.10e-03, 2.86e-04] m/s, advective CFL = 0.7, diffusive CFL = 0.32, νₘₐₓ = 0.0029 m²/s
┌ Info: iter = 800, [049.17%] time = 3.934 days, Δt = 8.445 minutes, walltime = 2.393 minutes, walltime / timestep = 122.435 ms
└ |u⃗|ₘₐₓ = [1.96e-02, 5.31e-03, 3.65e-04] m/s, advective CFL = 0.7, diffusive CFL = 0.35, νₘₐₓ = 0.0027 m²/s
┌ Info: iter = 900, [056.25%] time = 4.500 days, Δt = 8.343 minutes, walltime = 2.597 minutes, walltime / timestep = 121.823 ms
└ |u⃗|ₘₐₓ = [1.95e-02, 6.36e-03, 2.85e-04] m/s, advective CFL = 0.7, diffusive CFL = 0.33, νₘₐₓ = 0.0026 m²/s
┌ Info: iter = 1000, [063.20%] time = 5.056 days, Δt = 8.087 minutes, walltime = 2.828 minutes, walltime / timestep = 138.951 ms
└ |u⃗|ₘₐₓ = [1.88e-02, 7.83e-03, 3.45e-04] m/s, advective CFL = 0.7, diffusive CFL = 0.32, νₘₐₓ = 0.0027 m²/s
┌ Info: iter = 1100, [070.01%] time = 5.601 days, Δt = 8.115 minutes, walltime = 3.031 minutes, walltime / timestep = 121.556 ms
└ |u⃗|ₘₐₓ = [1.84e-02, 1.01e-02, 6.39e-04] m/s, advective CFL = 0.7, diffusive CFL = 0.33, νₘₐₓ = 0.0027 m²/s
┌ Info: iter = 1200, [076.77%] time = 6.141 days, Δt = 8.002 minutes, walltime = 3.267 minutes, walltime / timestep = 141.625 ms
└ |u⃗|ₘₐₓ = [1.79e-02, 1.19e-02, 1.36e-03] m/s, advective CFL = 0.7, diffusive CFL = 0.39, νₘₐₓ = 0.0032 m²/s
┌ Info: iter = 1300, [083.25%] time = 6.660 days, Δt = 7.497 minutes, walltime = 3.464 minutes, walltime / timestep = 118.627 ms
└ |u⃗|ₘₐₓ = [1.64e-02, 1.16e-02, 1.66e-03] m/s, advective CFL = 0.7, diffusive CFL = 0.38, νₘₐₓ = 0.0034 m²/s
┌ Info: iter = 1400, [089.43%] time = 7.155 days, Δt = 7.734 minutes, walltime = 3.660 minutes, walltime / timestep = 117.432 ms
└ |u⃗|ₘₐₓ = [1.53e-02, 1.12e-02, 1.68e-03] m/s, advective CFL = 0.7, diffusive CFL = 0.38, νₘₐₓ = 0.0032 m²/s
┌ Info: iter = 1500, [096.18%] time = 7.694 days, Δt = 7.983 minutes, walltime = 3.867 minutes, walltime / timestep = 124.145 ms
└ |u⃗|ₘₐₓ = [1.38e-02, 1.11e-02, 1.55e-03] m/s, advective CFL = 0.7, diffusive CFL = 0.42, νₘₐₓ = 0.0035 m²/s
[ Info: Simulation is stopping after running for 3.621 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, :]MLD_arr = ds["MLD"][:, :, :] # (Center, Center, Nothing)-located: no z dimension to sliceclose(ds)ζlim = maximum(abs, ζ_arr)νlim = maximum(ν_arr)0.0055754795f0Every budget record carries both tendencies and every source term at the same time. The first record has no earlier state to difference against and is written as zero, so the budgets start from the second.
ds_b = NCDataset(simulation.output_writers[:budget].filepath)nb = 2:length(ds_b["time"])t_bud = ds_b["time"][nb]deₚdt = ds_b["∂ₜ∫eₚ"][nb]deₖdt = ds_b["∂ₜ∫eₖ"][nb]Φ_bud = ds_b["∫Φ"][nb]wb_bud = ds_b["∫wb"][nb]εₖ_bud = ds_b["∫εₖ"][nb]close(ds_b)closed DatasetBoth 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_bud + Φ_budeₖ_resid = @. -deₖdt + wb_bud - εₖ_budPlotting
set_theme!(Theme(fontsize = 18))fig = Figure(size = (1000, 1800))n = Observable(1)panel_kwargs = (xlabel = "x [m]", ylabel = "y [m]", aspect = DataAspect())axζ = Axis(fig[2, 1][1, 1]; title = "vertical vorticity, ζ", panel_kwargs...)axb = Axis(fig[2, 2][1, 1]; title = "surface buoyancy, b", panel_kwargs...)axν = Axis(fig[3, 1][1, 1]; title = "eddy viscosity, νₑ", panel_kwargs...)axh = Axis(fig[3, 2][1, 1]; title = "mixed layer depth, MLD", panel_kwargs...)ζₙ = @lift ζ_arr[:, :, $n]bₙ = @lift b_arr[:, :, $n]νₙ = @lift ν_arr[:, :, $n]MLDₙ = @lift MLD_arr[:, :, $n]hmζ = heatmap!(axζ, x_faa, y_afa, ζₙ; colormap = :balance, colorrange = (-ζlim, ζlim))Colorbar(fig[2, 1][1, 2], hmζ)hmb = heatmap!(axb, x_caa, y_aca, bₙ; colormap = :thermal)Colorbar(fig[2, 2][1, 2], hmb)hmν = heatmap!(axν, x_caa, y_aca, νₙ; colormap = :tempo, colorrange = (0, νlim))Colorbar(fig[3, 1][1, 2], hmν)hmh = heatmap!(axh, x_caa, y_aca, MLDₙ; colormap = :deep, colorrange = (-h, 0))Colorbar(fig[3, 2][1, 2], hmh)budget_kwargs = (xlabel = "time [days]", ylabel = "[m⁵ s⁻³]")ax_p = Axis(fig[4, 1:2]; title = "Volume-integrated potential energy budget", budget_kwargs...)lines!(ax_p, t_bud ./ day, -deₚdt, label = "-d(∫eₚ)/dt")lines!(ax_p, t_bud ./ day, -wb_bud, label = "-∫wb dV")lines!(ax_p, t_bud ./ day, Φ_bud, label = "∫Φ dV")lines!(ax_p, t_bud ./ day, eₚ_resid, label = "residual", color = :black, linestyle = :dash)axislegend(ax_p; position = :rt, labelsize = 16, nbanks = 2)ax_k = Axis(fig[5, 1:2]; title = "Volume-integrated kinetic energy budget", budget_kwargs...)lines!(ax_k, t_bud ./ day, -deₖdt, label = "-d(∫eₖ)/dt")lines!(ax_k, t_bud ./ day, wb_bud, label = "∫wb dV")lines!(ax_k, t_bud ./ day, -εₖ_bud, label = "-∫εₖ dV")lines!(ax_k, t_bud ./ day, eₖ_resid, label = "residual", color = :black, linestyle = :dash)axislegend(ax_k; position = :rt, labelsize = 16, 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:2] = Label(fig, title, fontsize = 22, tellwidth = false)rowsize!(fig.layout, 1, Fixed(50))rowsize!(fig.layout, 2, Fixed(500))rowsize!(fig.layout, 3, Fixed(500))rowsize!(fig.layout, 4, Fixed(350))rowsize!(fig.layout, 5, Fixed(350))colsize!(fig.layout, 1, Fixed(450))colsize!(fig.layout, 2, Fixed(450))resize_to_layout!(fig)@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.