Baroclinic adjustment

In this example, we simulate the evolution and equilibration of a baroclinically unstable front.

Install dependencies

First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, CairoMakie"

using Oceananigans
using Oceananigans.Units

Grid

We use a three-dimensional channel that is periodic in the x direction:

Lx = 1000kilometers # east-west extent [m]
Ly = 1000kilometers # north-south extent [m]
Lz = 1kilometers    # depth [m]

grid = RectilinearGrid(size = (48, 48, 8),
                       x = (0, Lx),
                       y = (-Ly/2, Ly/2),
                       z = (-Lz, 0),
                       topology = (Periodic, Bounded, Bounded))

48×48×8 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0e6)          regularly spaced with Δx=20833.3
├── Bounded  y ∈ [-500000.0, 500000.0] regularly spaced with Δy=20833.3
└── Bounded  z ∈ [-1000.0, 0.0]        regularly spaced with Δz=125.0

Model

We built a HydrostaticFreeSurfaceModel with an ImplicitFreeSurface solver. Regarding Coriolis, we use a beta-plane centered at 45° South.

model = HydrostaticFreeSurfaceModel(grid;
                                    coriolis = BetaPlane(latitude = -45),
                                    buoyancy = BuoyancyTracer(),
                                    tracers = :b,
                                    momentum_advection = WENO(),
                                    tracer_advection = WENO())

HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 48×48×8 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── solver: FFTImplicitFreeSurfaceSolver
├── advection scheme: 
│   ├── momentum: WENO{3, Float64, Float32}(order=5)
│   └── b: WENO{3, Float64, Float32}(order=5)
├── vertical_coordinate: ZCoordinate
└── coriolis: BetaPlane{Float64}

We start our simulation from rest with a baroclinically unstable buoyancy distribution. We use ramp(y, Δy), defined below, to specify a front with width Δy and horizontal buoyancy gradient M². We impose the front on top of a vertical buoyancy gradient N² and a bit of noise.

"""
    ramp(y, Δy)

Linear ramp from 0 to 1 between -Δy/2 and +Δy/2.

For example:
```
            y < -Δy/2 => ramp = 0
    -Δy/2 < y < -Δy/2 => ramp = y / Δy
            y >  Δy/2 => ramp = 1
```
"""
ramp(y, Δy) = min(max(0, y/Δy + 1/2), 1)

N² = 1e-5 # [s⁻²] buoyancy frequency / stratification
M² = 1e-7 # [s⁻²] horizontal buoyancy gradient

Δy = 100kilometers # width of the region of the front
Δb = Δy * M²       # buoyancy jump associated with the front
ϵb = 1e-2 * Δb     # noise amplitude

bᵢ(x, y, z) = N² * z + Δb * ramp(y, Δy) + ϵb * randn()

set!(model, b=bᵢ)

Let's visualize the initial buoyancy distribution.

using CairoMakie
set_theme!(Theme(fontsize = 20))

# Build coordinates with units of kilometers
x, y, z = 1e-3 .* nodes(grid, (Center(), Center(), Center()))

b = model.tracers.b

fig, ax, hm = heatmap(view(b, 1, :, :),
                      colormap = :deep,
                      axis = (xlabel = "y [km]",
                              ylabel = "z [km]",
                              title = "b(x=0, y, z, t=0)",
                              titlesize = 24))

Colorbar(fig[1, 2], hm, label = "[m s⁻²]")

fig

Simulation

Now let's build a Simulation.

simulation = Simulation(model, Δt=20minutes, stop_time=20days)

Simulation of HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 20 minutes
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 20 days
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: OrderedDict with no entries

We add a TimeStepWizard callback to adapt the simulation's time-step,

conjure_time_step_wizard!(simulation, IterationInterval(20), cfl=0.2, max_Δt=20minutes)

Also, we add a callback to print a message about how the simulation is going,

using Printf

wall_clock = Ref(time_ns())

function print_progress(sim)
    u, v, w = model.velocities
    progress = 100 * (time(sim) / sim.stop_time)
    elapsed = (time_ns() - wall_clock[]) / 1e9

    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            progress, iteration(sim), prettytime(sim), prettytime(elapsed),
            maximum(abs, u), maximum(abs, v), maximum(abs, w), prettytime(sim.Δt))

    wall_clock[] = time_ns()

    return nothing
end

add_callback!(simulation, print_progress, IterationInterval(100))

Diagnostics/Output

Here, we save the buoyancy, $b$, at the edges of our domain as well as the zonal ($x$) average of buoyancy.

u, v, w = model.velocities
ζ = ∂x(v) - ∂y(u)
B = Average(b, dims=1)
U = Average(u, dims=1)
V = Average(v, dims=1)

filename = "baroclinic_adjustment"
save_fields_interval = 0.5day

slicers = (east = (grid.Nx, :, :),
           north = (:, grid.Ny, :),
           bottom = (:, :, 1),
           top = (:, :, grid.Nz))

for side in keys(slicers)
    indices = slicers[side]

    simulation.output_writers[side] = JLD2Writer(model, (; b, ζ);
                                                 filename = filename * "_$(side)_slice",
                                                 schedule = TimeInterval(save_fields_interval),
                                                 overwrite_existing = true,
                                                 indices)
end

simulation.output_writers[:zonal] = JLD2Writer(model, (; b=B, u=U, v=V);
                                               filename = filename * "_zonal_average",
                                               schedule = TimeInterval(save_fields_interval),
                                               overwrite_existing = true)

JLD2Writer scheduled on TimeInterval(12 hours):
├── filepath: baroclinic_adjustment_zonal_average.jld2
├── 3 outputs: (b, u, v)
├── array_type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

Now we're ready to run.

@info "Running the simulation..."

run!(simulation)

@info "Simulation completed in " * prettytime(simulation.run_wall_time)

[ Info: Running the simulation...
[ Info: Initializing simulation...
[00.00%] i: 0, t: 0 seconds, wall time: 10.006 seconds, max(u): (0.000e+00, 0.000e+00, 0.000e+00) m/s, next Δt: 20 minutes
[ Info:     ... simulation initialization complete (13.386 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (3.738 seconds).
[06.94%] i: 100, t: 1.389 days, wall time: 15.696 seconds, max(u): (1.203e-01, 1.199e-01, 1.697e-03) m/s, next Δt: 20 minutes
[13.89%] i: 200, t: 2.778 days, wall time: 1.224 seconds, max(u): (2.115e-01, 2.018e-01, 1.975e-03) m/s, next Δt: 20 minutes
[20.83%] i: 300, t: 4.167 days, wall time: 955.523 ms, max(u): (2.929e-01, 2.864e-01, 1.929e-03) m/s, next Δt: 20 minutes
[27.78%] i: 400, t: 5.556 days, wall time: 862.593 ms, max(u): (3.865e-01, 4.059e-01, 1.908e-03) m/s, next Δt: 20 minutes
[34.72%] i: 500, t: 6.944 days, wall time: 827.374 ms, max(u): (4.905e-01, 5.763e-01, 2.241e-03) m/s, next Δt: 20 minutes
[41.67%] i: 600, t: 8.333 days, wall time: 859.111 ms, max(u): (7.433e-01, 1.002e+00, 3.352e-03) m/s, next Δt: 20 minutes
[48.61%] i: 700, t: 9.722 days, wall time: 1.065 seconds, max(u): (1.230e+00, 1.174e+00, 4.743e-03) m/s, next Δt: 20 minutes
[55.56%] i: 800, t: 11.111 days, wall time: 890.884 ms, max(u): (1.445e+00, 1.192e+00, 4.867e-03) m/s, next Δt: 20 minutes
[62.50%] i: 900, t: 12.500 days, wall time: 893.758 ms, max(u): (1.356e+00, 1.200e+00, 4.645e-03) m/s, next Δt: 20 minutes
[69.44%] i: 1000, t: 13.889 days, wall time: 865.962 ms, max(u): (1.318e+00, 1.168e+00, 5.551e-03) m/s, next Δt: 20 minutes
[76.39%] i: 1100, t: 15.278 days, wall time: 927.860 ms, max(u): (1.336e+00, 1.183e+00, 4.386e-03) m/s, next Δt: 20 minutes
[83.33%] i: 1200, t: 16.667 days, wall time: 930.854 ms, max(u): (1.373e+00, 1.123e+00, 2.855e-03) m/s, next Δt: 20 minutes
[90.28%] i: 1300, t: 18.056 days, wall time: 879.198 ms, max(u): (1.575e+00, 1.059e+00, 2.518e-03) m/s, next Δt: 20 minutes
[97.22%] i: 1400, t: 19.444 days, wall time: 924.097 ms, max(u): (1.284e+00, 1.123e+00, 2.962e-03) m/s, next Δt: 20 minutes
[ Info: Simulation is stopping after running for 30.600 seconds.
[ Info: Simulation time 20 days equals or exceeds stop time 20 days.
[ Info: Simulation completed in 30.642 seconds

Visualization

All that's left is to make a pretty movie. Actually, we make two visualizations here. First, we illustrate how to make a 3D visualization with Makie's Axis3 and Makie.surface. Then we make a movie in 2D. We use CairoMakie in this example, but note that using GLMakie is more convenient on a system with OpenGL, as figures will be displayed on the screen.

using CairoMakie

Three-dimensional visualization

We load the saved buoyancy output on the top, north, and east surface as FieldTimeSerieses.

filename = "baroclinic_adjustment"

sides = keys(slicers)

slice_filenames = NamedTuple(side => filename * "_$(side)_slice.jld2" for side in sides)

b_timeserieses = (east   = FieldTimeSeries(slice_filenames.east, "b"),
                  north  = FieldTimeSeries(slice_filenames.north, "b"),
                  top    = FieldTimeSeries(slice_filenames.top, "b"))

B_timeseries = FieldTimeSeries(filename * "_zonal_average.jld2", "b")

times = B_timeseries.times
grid = B_timeseries.grid

48×48×8 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0e6)          regularly spaced with Δx=20833.3
├── Bounded  y ∈ [-500000.0, 500000.0] regularly spaced with Δy=20833.3
└── Bounded  z ∈ [-1000.0, 0.0]        regularly spaced with Δz=125.0

We build the coordinates. We rescale horizontal coordinates to kilometers.

xb, yb, zb = nodes(b_timeserieses.east)

xb = xb ./ 1e3 # convert m -> km
yb = yb ./ 1e3 # convert m -> km

Nx, Ny, Nz = size(grid)

x_xz = repeat(x, 1, Nz)
y_xz_north = y[end] * ones(Nx, Nz)
z_xz = repeat(reshape(z, 1, Nz), Nx, 1)

x_yz_east = x[end] * ones(Ny, Nz)
y_yz = repeat(y, 1, Nz)
z_yz = repeat(reshape(z, 1, Nz), grid.Ny, 1)

x_xy = x
y_xy = y
z_xy_top = z[end] * ones(grid.Nx, grid.Ny)

Then we create a 3D axis. We use zonal_slice_displacement to control where the plot of the instantaneous zonal average flow is located.

fig = Figure(size = (1600, 800))

zonal_slice_displacement = 1.2

ax = Axis3(fig[2, 1],
           aspect=(1, 1, 1/5),
           xlabel = "x (km)",
           ylabel = "y (km)",
           zlabel = "z (m)",
           xlabeloffset = 100,
           ylabeloffset = 100,
           zlabeloffset = 100,
           limits = ((x[1], zonal_slice_displacement * x[end]), (y[1], y[end]), (z[1], z[end])),
           elevation = 0.45,
           azimuth = 6.8,
           xspinesvisible = false,
           zgridvisible = false,
           protrusions = 40,
           perspectiveness = 0.7)

Axis3 with 12 plots:
 ┣━ Poly{Tuple{GeometryBasics.Polygon{2, Float64}}}
 ┣━ Poly{Tuple{GeometryBasics.Polygon{2, Float64}}}
 ┣━ Poly{Tuple{GeometryBasics.Polygon{2, Float64}}}
 ┣━ LineSegments{Tuple{Base.ReinterpretArray{Point{3, Float64}, 1, Tuple{Point{3, Float64}, Point{3, Float64}}, Vector{Tuple{Point{3, Float64}, Point{3, Float64}}}, false}}}
 ┣━ LineSegments{Tuple{Base.ReinterpretArray{Point{3, Float64}, 1, Tuple{Point{3, Float64}, Point{3, Float64}}, Vector{Tuple{Point{3, Float64}, Point{3, Float64}}}, false}}}
 ┣━ LineSegments{Tuple{Vector{Point{3, Float64}}}}
 ┣━ LineSegments{Tuple{Base.ReinterpretArray{Point{3, Float64}, 1, Tuple{Point{3, Float64}, Point{3, Float64}}, Vector{Tuple{Point{3, Float64}, Point{3, Float64}}}, false}}}
 ┣━ LineSegments{Tuple{Base.ReinterpretArray{Point{3, Float64}, 1, Tuple{Point{3, Float64}, Point{3, Float64}}, Vector{Tuple{Point{3, Float64}, Point{3, Float64}}}, false}}}
 ┣━ LineSegments{Tuple{Vector{Point{3, Float64}}}}
 ┣━ LineSegments{Tuple{Base.ReinterpretArray{Point{3, Float64}, 1, Tuple{Point{3, Float64}, Point{3, Float64}}, Vector{Tuple{Point{3, Float64}, Point{3, Float64}}}, false}}}
 ┣━ LineSegments{Tuple{Base.ReinterpretArray{Point{3, Float64}, 1, Tuple{Point{3, Float64}, Point{3, Float64}}, Vector{Tuple{Point{3, Float64}, Point{3, Float64}}}, false}}}
 ┗━ LineSegments{Tuple{Vector{Point{3, Float64}}}}

We use data from the final savepoint for the 3D plot. Note that this plot can easily be animated by using Makie's Observable. To dive into Observables, check out Makie.jl's Documentation.

n = length(times)

41

Now let's make a 3D plot of the buoyancy and in front of it we'll use the zonally-averaged output to plot the instantaneous zonal-average of the buoyancy.

b_slices = (east   = interior(b_timeserieses.east[n], 1, :, :),
            north  = interior(b_timeserieses.north[n], :, 1, :),
            top    = interior(b_timeserieses.top[n], :, :, 1))

# Zonally-averaged buoyancy
B = interior(B_timeseries[n], 1, :, :)

clims = 1.1 .* extrema(b_timeserieses.top[n][:])

kwargs = (colorrange=clims, colormap=:deep, shading=NoShading)

surface!(ax, x_yz_east, y_yz, z_yz;  color = b_slices.east, kwargs...)
surface!(ax, x_xz, y_xz_north, z_xz; color = b_slices.north, kwargs...)
surface!(ax, x_xy, y_xy, z_xy_top;   color = b_slices.top, kwargs...)

sf = surface!(ax, zonal_slice_displacement .* x_yz_east, y_yz, z_yz; color = B, kwargs...)

contour!(ax, y, z, B; transformation = (:yz, zonal_slice_displacement * x[end]),
         levels = 15, linewidth = 2, color = :black)

Colorbar(fig[2, 2], sf, label = "m s⁻²", height = Relative(0.4), tellheight=false)

title = "Buoyancy at t = " * string(round(times[n] / day, digits=1)) * " days"
fig[1, 1:2] = Label(fig, title; fontsize = 24, tellwidth = false, padding = (0, 0, -120, 0))

rowgap!(fig.layout, 1, Relative(-0.2))
colgap!(fig.layout, 1, Relative(-0.1))

save("baroclinic_adjustment_3d.png", fig)

Two-dimensional movie

We make a 2D movie that shows buoyancy $b$ and vertical vorticity $ζ$ at the surface, as well as the zonally-averaged zonal and meridional velocities $U$ and $V$ in the $(y, z)$ plane. First we load the FieldTimeSeries and extract the additional coordinates we'll need for plotting

ζ_timeseries = FieldTimeSeries(slice_filenames.top, "ζ")
U_timeseries = FieldTimeSeries(filename * "_zonal_average.jld2", "u")
B_timeseries = FieldTimeSeries(filename * "_zonal_average.jld2", "b")
V_timeseries = FieldTimeSeries(filename * "_zonal_average.jld2", "v")

xζ, yζ, zζ = nodes(ζ_timeseries)
yv = ynodes(V_timeseries)

xζ = xζ ./ 1e3 # convert m -> km
yζ = yζ ./ 1e3 # convert m -> km
yv = yv ./ 1e3 # convert m -> km

-500.0:20.833333333333332:500.0

Next, we set up a plot with 4 panels. The top panels are large and square, while the bottom panels get a reduced aspect ratio through rowsize!.

fig = Figure(size=(1800, 1000))

axb = Axis(fig[1, 2], xlabel="x (km)", ylabel="y (km)", aspect=1)
axζ = Axis(fig[1, 3], xlabel="x (km)", ylabel="y (km)", aspect=1, yaxisposition=:right)

axu = Axis(fig[2, 2], xlabel="y (km)", ylabel="z (m)")
axv = Axis(fig[2, 3], xlabel="y (km)", ylabel="z (m)", yaxisposition=:right)

rowsize!(fig.layout, 2, Relative(0.3))

To prepare a plot for animation, we index the timeseries with an Observable,

n = Observable(1)

b_top = @lift interior(b_timeserieses.top[$n], :, :, 1)
ζ_top = @lift interior(ζ_timeseries[$n], :, :, 1)
U = @lift interior(U_timeseries[$n], 1, :, :)
V = @lift interior(V_timeseries[$n], 1, :, :)
B = @lift interior(B_timeseries[$n], 1, :, :)

Observable([-0.009375827386975288 -0.008115331642329693 -0.006863572634756565 -0.005603509489446878 -0.0043804459273815155 -0.003082681680098176 -0.0018720811931416392 -0.0006198332412168384; -0.009389028884470463 -0.008108125999569893 -0.006888588424772024 -0.005608034320175648 -0.004391692113131285 -0.0031215734779834747 -0.0018882937729358673 -0.0006314754718914628; -0.00937504693865776 -0.00812914315611124 -0.006900915876030922 -0.005610805004835129 -0.0043731979094445705 -0.0031205632258206606 -0.0018719263607636094 -0.0005882501718588173; -0.009362530894577503 -0.00813062023371458 -0.00686543807387352 -0.0056110708974301815 -0.004375075455754995 -0.0031367165502160788 -0.0018835797673091292 -0.0006136768497526646; -0.00939097162336111 -0.008145217783749104 -0.006880767177790403 -0.005610702559351921 -0.004367046058177948 -0.0031253276392817497 -0.0018593142740428448 -0.0006105767679400742; -0.009387333877384663 -0.008131183683872223 -0.006868558004498482 -0.00562877394258976 -0.004388961475342512 -0.0031285472214221954 -0.0018581830663606524 -0.0006302017136476934; -0.009390502236783504 -0.00813261978328228 -0.006878410466015339 -0.005656544119119644 -0.0043508149683475494 -0.0031609931029379368 -0.0018707929411903024 -0.0006238670321181417; -0.009367218241095543 -0.00812871940433979 -0.006863442715257406 -0.005648050457239151 -0.004363846033811569 -0.0031042692717164755 -0.0018723832909017801 -0.000642432423774153; -0.009364893659949303 -0.008124938234686852 -0.00687911082059145 -0.005637975409626961 -0.004398196469992399 -0.003155940445140004 -0.00188451015856117 -0.0006222432712092996; -0.009365012869238853 -0.008109252899885178 -0.006872318685054779 -0.00560746667906642 -0.004393348004668951 -0.0031277681700885296 -0.001892683794721961 -0.0006040071602910757; -0.009354833513498306 -0.008124666288495064 -0.006891452707350254 -0.005621962249279022 -0.004370904993265867 -0.003124233800917864 -0.0018820477416738868 -0.0006421724683605134; -0.009382533840835094 -0.008131904527544975 -0.00688491016626358 -0.005637418944388628 -0.004381454084068537 -0.0031363628804683685 -0.0018650087295100093 -0.0006447449559345841; -0.009364292025566101 -0.008116628043353558 -0.006895150989294052 -0.005619329400360584 -0.0043585943058133125 -0.0031167089473456144 -0.0018730968004092574 -0.0006326051079668105; -0.009386088699102402 -0.008101114071905613 -0.006887619849294424 -0.005645467434078455 -0.00439302297309041 -0.0031352490186691284 -0.0018653955776244402 -0.0006170714041218162; -0.009373744949698448 -0.008105888031423092 -0.00687916437163949 -0.005615980830043554 -0.004398385528475046 -0.0031093035358935595 -0.001874720212072134 -0.0006396498065441847; -0.009378956630825996 -0.00812617875635624 -0.0068449475802481174 -0.005638033617287874 -0.0043646241538226604 -0.0031073507852852345 -0.0018473061500117183 -0.0006225585821084678; -0.009391129948198795 -0.008097817189991474 -0.00688178138807416 -0.005640231538563967 -0.004388230387121439 -0.003119709901511669 -0.0018697087652981281 -0.0006211657309904695; -0.009363553486764431 -0.008144947700202465 -0.006876219995319843 -0.00561447162181139 -0.004380031954497099 -0.0031223997939378023 -0.0018853064393624663 -0.0006196764297783375; -0.009381498210132122 -0.008124117739498615 -0.006876965053379536 -0.005622352473437786 -0.004378326237201691 -0.0031231811735779047 -0.0018886730540543795 -0.0006228406564332545; -0.009371857158839703 -0.008123648352921009 -0.006906024180352688 -0.005607191007584333 -0.0043906886130571365 -0.003120406297966838 -0.0018841439159587026 -0.0006344981375150383; -0.00937062967568636 -0.008120059035718441 -0.00688138185068965 -0.005622204393148422 -0.00439249025657773 -0.003132396610453725 -0.0018874832894653082 -0.0006402337457984686; -0.009351012296974659 -0.008144540712237358 -0.0068702735006809235 -0.005651468876749277 -0.0044003804214298725 -0.0030905886087566614 -0.001887109479866922 -0.0006336977821774781; -0.007523529697209597 -0.006239814683794975 -0.00498583447188139 -0.0037630884908139706 -0.002511846600100398 -0.0012456688564270735 1.8546379578765482e-5 0.0012333592167124152; -0.005399083252996206 -0.004136588424444199 -0.0029300199821591377 -0.0016687465831637383 -0.0004145460552535951 0.0008268540259450674 0.002056600758805871 0.003345222445204854; -0.0033382968977093697 -0.002077592769637704 -0.0008231117390096188 0.0003866321058012545 0.0016703458968549967 0.00291971443220973 0.004178056493401527 0.005410141311585903; -0.0012398229446262121 1.2650294593186118e-5 0.0012564744101837277 0.0024854247458279133 0.003729530842974782 0.004984708968549967 0.006237974856048822 0.007489944342523813; 0.0006297442014329135 0.001901932992041111 0.0031278699170798063 0.004353807773441076 0.005604709964245558 0.006875153630971909 0.008128521032631397 0.009381299838423729; 0.0006115239229984581 0.0018665243405848742 0.003136987565085292 0.004372523166239262 0.005630094092339277 0.006868950091302395 0.008093606680631638 0.009398231282830238; 0.0006319331587292254 0.0018894035601988435 0.0031201739329844713 0.004360498394817114 0.00562030915170908 0.006854303181171417 0.00813127402216196 0.0093600545078516; 0.0005870038294233382 0.0018695889739319682 0.003108419245108962 0.004389649257063866 0.00560262193903327 0.006903393659740686 0.008126487955451012 0.009367507882416248; 0.0006445850594900548 0.0018698772182688117 0.0031389149371534586 0.004377532750368118 0.005592058878391981 0.006852137856185436 0.008110196329653263 0.009377486072480679; 0.0006269316654652357 0.0018871191423386335 0.0031300412956625223 0.00434942077845335 0.005637639202177525 0.006885149981826544 0.00812543649226427 0.009362134151160717; 0.0006144155049696565 0.0018579558236524463 0.003084383672103286 0.004353346768766642 0.005620869342237711 0.006893177516758442 0.008125100284814835 0.009389102458953857; 0.000646879430860281 0.0018549795495346189 0.0031409556977450848 0.004361655563116074 0.005640025716274977 0.006856212392449379 0.00810568779706955 0.0093765240162611; 0.0006160188931971788 0.0018892448861151934 0.003129462944343686 0.004358119796961546 0.005627302452921867 0.006873182021081448 0.00813320092856884 0.00937605556100607; 0.0006203322554938495 0.0018723472021520138 0.003124258480966091 0.004370648879557848 0.005630529019981623 0.006857578177005053 0.00813406053930521 0.009374463930726051; 0.0006082246545702219 0.0018709335708990693 0.0031114777084439993 0.004376907832920551 0.005587043706327677 0.006889630109071732 0.008139167912304401 0.009375883266329765; 0.0006314062047749758 0.0018740821396932006 0.003125484101474285 0.004367649555206299 0.005627807695418596 0.006866486743092537 0.008147566579282284 0.009383738972246647; 0.0006103495834395289 0.0018583409255370498 0.003128924872726202 0.004364558961242437 0.005653732921928167 0.006897184532135725 0.008123865351080894 0.009354389272630215; 0.0006553507409989834 0.0018758268561214209 0.003111141500994563 0.004381299018859863 0.00560422707349062 0.006862929090857506 0.008147510699927807 0.009353441186249256; 0.0006391274509951472 0.0018756308127194643 0.00312030129134655 0.0043699308298528194 0.005600322969257832 0.00688052736222744 0.008110073395073414 0.009392245672643185; 0.0006393387448042631 0.0018752857577055693 0.0031192402821034193 0.004365513566881418 0.0056174201890826225 0.006891972851008177 0.008161993697285652 0.009370330721139908; 0.0006243929965421557 0.0018807636806741357 0.0031122458167374134 0.0043747685849666595 0.005635312292724848 0.0068741594441235065 0.008128662593662739 0.009397156536579132; 0.0006463482277467847 0.0018762322142720222 0.003137442981824279 0.004377619829028845 0.0056273238733410835 0.006855629850178957 0.008130747824907303 0.009353269822895527; 0.0006259552901610732 0.001894457614980638 0.0031009891536086798 0.004384420812129974 0.005630373489111662 0.0068703447468578815 0.008140373975038528 0.009372220374643803; 0.0006673220777884126 0.0018733944743871689 0.003123084781691432 0.004375624004751444 0.005658276379108429 0.006854484789073467 0.008122681640088558 0.009356718510389328; 0.0006122922641225159 0.0019074927549809217 0.003130801720544696 0.004377448000013828 0.005631558131426573 0.006878209766000509 0.008112089708447456 0.009359347634017467; 0.0006096081342548132 0.0018842784920707345 0.0031412714160978794 0.004380751401185989 0.005640337243676186 0.006879671476781368 0.00811453815549612 0.009354242123663425])

and then build our plot:

hm = heatmap!(axb, xb, yb, b_top, colorrange=(0, Δb), colormap=:thermal)
Colorbar(fig[1, 1], hm, flipaxis=false, label="Surface b(x, y) (m s⁻²)")

hm = heatmap!(axζ, xζ, yζ, ζ_top, colorrange=(-5e-5, 5e-5), colormap=:balance)
Colorbar(fig[1, 4], hm, label="Surface ζ(x, y) (s⁻¹)")

hm = heatmap!(axu, yb, zb, U; colorrange=(-5e-1, 5e-1), colormap=:balance)
Colorbar(fig[2, 1], hm, flipaxis=false, label="Zonally-averaged U(y, z) (m s⁻¹)")
contour!(axu, yb, zb, B; levels=15, color=:black)

hm = heatmap!(axv, yv, zb, V; colorrange=(-1e-1, 1e-1), colormap=:balance)
Colorbar(fig[2, 4], hm, label="Zonally-averaged V(y, z) (m s⁻¹)")
contour!(axv, yb, zb, B; levels=15, color=:black)

Finally, we're ready to record the movie.

frames = 1:length(times)

record(fig, filename * ".mp4", frames, framerate=8) do i
    n[] = i
end

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.12.4
Commit 01a2eadb047 (2026-01-06 16:56 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 128 × AMD EPYC 9374F 32-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
Environment:
  LD_LIBRARY_PATH = 
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia-oceananigans
  JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-29180/docs/
  JULIA_VERSION = 1.12.4
  JULIA_LOAD_PATH = @:@v#.#:@stdlib
  JULIA_VERSION_ENZYME = 1.10.10
  JULIA_PYTHONCALL_EXE = /var/lib/buildkite-agent/Oceananigans.jl-29180/docs/.CondaPkg/.pixi/envs/default/bin/python
  JULIA_DEBUG = Literate

These were the top-level packages installed in the environment:

import Pkg
Pkg.status()

Status `~/Oceananigans.jl-29180/docs/Project.toml`
  [79e6a3ab] Adapt v4.4.0
  [052768ef] CUDA v5.9.6
  [13f3f980] CairoMakie v0.15.8
  [e30172f5] Documenter v1.16.1
  [daee34ce] DocumenterCitations v1.4.1
  [033835bb] JLD2 v0.6.3
  [63c18a36] KernelAbstractions v0.9.39
  [98b081ad] Literate v2.21.0
  [da04e1cc] MPI v0.20.23
  [85f8d34a] NCDatasets v0.14.11
  [9e8cae18] Oceananigans v0.104.3 `..`
  [f27b6e38] Polynomials v4.1.0
  [6038ab10] Rotations v1.7.1
  [d496a93d] SeawaterPolynomials v0.3.10
  [09ab397b] StructArrays v0.7.2
  [bdfc003b] TimesDates v0.3.3
  [2e0b0046] XESMF v0.1.6
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.12.0
  [44cfe95a] Pkg v1.12.1

This page was generated using Literate.jl.