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]

Nx = 64
Ny = 64
Nz = 40

grid = RectilinearGrid(size = (Nx, Ny, Nz),
                       x = (0, Lx),
                       y = (-Ly/2, Ly/2),
                       z = (-Lz, 0),
                       topology = (Periodic, Bounded, Bounded))
64×64×40 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0e6)          regularly spaced with Δx=15625.0
├── Bounded  y ∈ [-500000.0, 500000.0] regularly spaced with Δy=15625.0
└── Bounded  z ∈ [-1000.0, 0.0]        regularly spaced with Δz=25.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: 64×64×40 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with -ĝ = ZDirection
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── solver: FFTImplicitFreeSurfaceSolver
└── coriolis: BetaPlane{Float64}

We want to initialize our model with a baroclinically unstable front plus some small-amplitude 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)

We then use ramp(y, Δy) to construct an initial buoyancy configuration of a baroclinically unstable front. The front has a buoyancy jump Δb over a latitudinal width Δy.

N² = 4e-6 # [s⁻²] buoyancy frequency / stratification
M² = 8e-8 # [s⁻²] horizontal buoyancy gradient

Δy = 50kilometers # 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

x, y, z = 1e-3 .* nodes((Center, Center, Center), grid) # convert m -> km

b = model.tracers.b

fig, ax, hm = heatmap(y, z, interior(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⁻²]")

Now let's built a Simulation.

Δt₀ = 5minutes
stop_time = 40days

simulation = Simulation(model, Δt=Δt₀, stop_time=stop_time)
Simulation of HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 5 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN years
├── Stop time: 40 days
├── Stop iteration : Inf
├── Wall time limit: Inf
├── 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
└── Diagnostics: OrderedDict with no entries

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

wizard = TimeStepWizard(cfl=0.2, max_change=1.1, max_Δt=20minutes)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(20))
Callback of TimeStepWizard(cfl=0.2, max_Δt=1200.0, min_Δt=0.0) on IterationInterval(20)

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

using Printf

wall_clock = [time_ns()]

function print_progress(sim)
    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            100 * (sim.model.clock.time / sim.stop_time),
            sim.model.clock.iteration,
            prettytime(sim.model.clock.time),
            prettytime(1e-9 * (time_ns() - wall_clock[1])),
            maximum(abs, sim.model.velocities.u),
            maximum(abs, sim.model.velocities.v),
            maximum(abs, sim.model.velocities.w),
            prettytime(sim.Δt))

    wall_clock[1] = time_ns()

    return nothing
end

simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(100))
Callback of print_progress on IterationInterval(100)

Diagnostics/Output

Add some diagnostics. 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

B = Field(Average(b, dims=1))

filename = "baroclinic_adjustment"
save_fields_interval = 0.5day

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

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

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

simulation.output_writers[:zonal] = JLD2OutputWriter(model, (b=B,);
                                                     filename = filename * "_zonal_average",
                                                     schedule = TimeInterval(save_fields_interval),
                                                     overwrite_existing = true)
JLD2OutputWriter scheduled on TimeInterval(12 hours):
├── filepath: ./baroclinic_adjustment_zonal_average.jld2
├── 1 outputs: b
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB

Now let's 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: 12.078 seconds, max(u): (0.000e+00, 0.000e+00, 0.000e+00) m/s, next Δt: 5.500 minutes
[ Info:     ... simulation initialization complete (11.787 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (1.461 minutes).
[01.17%] i: 100, t: 11.193 hours, wall time: 4.384 minutes, max(u): (6.972e-01, 3.336e-01, 3.438e-03) m/s, next Δt: 8.858 minutes
[03.03%] i: 200, t: 1.213 days, wall time: 2.798 minutes, max(u): (6.074e-01, 3.530e-01, 2.463e-03) m/s, next Δt: 14.266 minutes
[05.99%] i: 300, t: 2.396 days, wall time: 2.853 minutes, max(u): (5.596e-01, 2.916e-01, 2.090e-03) m/s, next Δt: 20 minutes
[09.44%] i: 400, t: 3.778 days, wall time: 2.891 minutes, max(u): (4.668e-01, 2.641e-01, 1.793e-03) m/s, next Δt: 20 minutes
[12.92%] i: 500, t: 5.167 days, wall time: 2.870 minutes, max(u): (4.459e-01, 2.486e-01, 1.533e-03) m/s, next Δt: 20 minutes
[16.39%] i: 600, t: 6.556 days, wall time: 2.795 minutes, max(u): (4.442e-01, 2.517e-01, 1.673e-03) m/s, next Δt: 20 minutes
[19.86%] i: 700, t: 7.944 days, wall time: 2.800 minutes, max(u): (4.592e-01, 3.615e-01, 2.140e-03) m/s, next Δt: 20 minutes
[23.33%] i: 800, t: 9.333 days, wall time: 2.801 minutes, max(u): (5.028e-01, 5.446e-01, 2.924e-03) m/s, next Δt: 20 minutes
[26.76%] i: 900, t: 10.704 days, wall time: 2.777 minutes, max(u): (5.834e-01, 6.706e-01, 4.024e-03) m/s, next Δt: 20 minutes
[29.84%] i: 1000, t: 11.937 days, wall time: 2.789 minutes, max(u): (7.672e-01, 7.477e-01, 4.965e-03) m/s, next Δt: 16.783 minutes
[32.39%] i: 1100, t: 12.956 days, wall time: 2.801 minutes, max(u): (8.439e-01, 8.187e-01, 6.074e-03) m/s, next Δt: 13.721 minutes
[34.80%] i: 1200, t: 13.921 days, wall time: 2.788 minutes, max(u): (9.442e-01, 8.045e-01, 5.485e-03) m/s, next Δt: 15.192 minutes
[37.74%] i: 1300, t: 15.096 days, wall time: 2.782 minutes, max(u): (8.119e-01, 9.187e-01, 4.123e-03) m/s, next Δt: 20 minutes
[41.18%] i: 1400, t: 16.472 days, wall time: 2.779 minutes, max(u): (9.490e-01, 8.187e-01, 4.513e-03) m/s, next Δt: 18.464 minutes
[44.55%] i: 1500, t: 17.819 days, wall time: 2.749 minutes, max(u): (1.012e+00, 8.352e-01, 3.783e-03) m/s, next Δt: 20 minutes
[48.02%] i: 1600, t: 19.208 days, wall time: 2.739 minutes, max(u): (8.846e-01, 9.608e-01, 3.701e-03) m/s, next Δt: 20 minutes
[51.49%] i: 1700, t: 20.597 days, wall time: 2.729 minutes, max(u): (8.097e-01, 9.447e-01, 3.862e-03) m/s, next Δt: 20 minutes
[54.97%] i: 1800, t: 21.986 days, wall time: 2.742 minutes, max(u): (8.443e-01, 8.793e-01, 3.152e-03) m/s, next Δt: 20 minutes
[58.44%] i: 1900, t: 23.375 days, wall time: 2.729 minutes, max(u): (8.467e-01, 7.550e-01, 2.807e-03) m/s, next Δt: 20 minutes
[61.91%] i: 2000, t: 24.764 days, wall time: 2.759 minutes, max(u): (8.410e-01, 8.492e-01, 2.356e-03) m/s, next Δt: 20 minutes
[65.38%] i: 2100, t: 26.153 days, wall time: 2.738 minutes, max(u): (7.638e-01, 8.596e-01, 2.728e-03) m/s, next Δt: 20 minutes
[68.85%] i: 2200, t: 27.542 days, wall time: 2.693 minutes, max(u): (8.790e-01, 8.182e-01, 2.616e-03) m/s, next Δt: 20 minutes
[72.33%] i: 2300, t: 28.931 days, wall time: 2.709 minutes, max(u): (8.934e-01, 8.540e-01, 2.853e-03) m/s, next Δt: 20 minutes
[75.80%] i: 2400, t: 30.319 days, wall time: 2.432 minutes, max(u): (9.855e-01, 8.855e-01, 2.471e-03) m/s, next Δt: 20 minutes
[79.27%] i: 2500, t: 31.708 days, wall time: 2.447 minutes, max(u): (1.133e+00, 9.073e-01, 2.510e-03) m/s, next Δt: 20 minutes
[82.74%] i: 2600, t: 33.097 days, wall time: 2.420 minutes, max(u): (9.845e-01, 9.360e-01, 2.425e-03) m/s, next Δt: 20 minutes
[86.22%] i: 2700, t: 34.486 days, wall time: 2.415 minutes, max(u): (9.758e-01, 8.624e-01, 2.588e-03) m/s, next Δt: 20 minutes
[89.69%] i: 2800, t: 35.875 days, wall time: 2.426 minutes, max(u): (9.523e-01, 9.636e-01, 1.928e-03) m/s, next Δt: 20 minutes
[93.16%] i: 2900, t: 37.264 days, wall time: 2.432 minutes, max(u): (8.566e-01, 1.021e+00, 2.106e-03) m/s, next Δt: 20 minutes
[96.63%] i: 3000, t: 38.653 days, wall time: 2.409 minutes, max(u): (8.816e-01, 1.004e+00, 2.466e-03) m/s, next Δt: 20 minutes
[ Info: Simulation is stopping after running for 1.414 hours.
[ Info: Simulation time 40 days equals or exceeds stop time 40 days.
[ Info: Simulation completed in 1.415 hours

Visualization

Now we are ready to visualize our resutls! We use CairoMakie in this example. On a system with OpenGL using GLMakie is more convenient as figures will be displayed on the screen.

using CairoMakie

We load the saved buoyancy output on the top, bottom, 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"),
                  bottom = FieldTimeSeries(slice_filenames.bottom, "b"),
                  top    = FieldTimeSeries(slice_filenames.top, "b"))

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

We build the coordinates. We rescale horizontal coordinates so that they correspond to kilometers.

x, y, z = nodes(b_timeserieses.east)

x = x .* 1e-3 # convert m -> km
y = y .* 1e-3 # convert m -> km

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)
z_xy_bottom = z[1] * 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(resolution = (900, 520))

zonal_slice_displacement = 1.2

ax = Axis3(fig[2, 1], aspect=(1, 1, 1/5),
           xlabel="x (km)", ylabel="y (km)", zlabel="z (m)",
           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)

We use Makie's Observable to animate the data. To dive into how Observables work we refer to Makie.jl's Documentation.

n = Observable(1)
Observable(1)

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   = @lift(interior(b_timeserieses.east[$n], 1, :, :)),
            north  = @lift(interior(b_timeserieses.north[$n], :, 1, :)),
            bottom = @lift(interior(b_timeserieses.bottom[$n], :, :, 1)),
            top    = @lift(interior(b_timeserieses.top[$n], :, :, 1)))

avg_b = @lift interior(avg_b_timeseries[$n], 1, :, :)

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

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

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_bottom ; color = b_slices.bottom, 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 = avg_b, kwargs...)

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

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

Finally, we add a figure title with the time of the snapshot and then record a movie.

times = avg_b_timeseries.times

title = @lift "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))

frames = 1:length(times)

record(fig, filename * ".mp4", frames, framerate=8) do i
    msg = string("Plotting frame ", i, " of ", frames[end])
    print(msg * " \r")
    n[] = i
end
Plotting frame 1 of 81 
Plotting frame 2 of 81 
Plotting frame 3 of 81 
Plotting frame 4 of 81 
Plotting frame 5 of 81 
Plotting frame 6 of 81 
Plotting frame 7 of 81 
Plotting frame 8 of 81 
Plotting frame 9 of 81 
Plotting frame 10 of 81 
Plotting frame 11 of 81 
Plotting frame 12 of 81 
Plotting frame 13 of 81 
Plotting frame 14 of 81 
Plotting frame 15 of 81 
Plotting frame 16 of 81 
Plotting frame 17 of 81 
Plotting frame 18 of 81 
Plotting frame 19 of 81 
Plotting frame 20 of 81 
Plotting frame 21 of 81 
Plotting frame 22 of 81 
Plotting frame 23 of 81 
Plotting frame 24 of 81 
Plotting frame 25 of 81 
Plotting frame 26 of 81 
Plotting frame 27 of 81 
Plotting frame 28 of 81 
Plotting frame 29 of 81 
Plotting frame 30 of 81 
Plotting frame 31 of 81 
Plotting frame 32 of 81 
Plotting frame 33 of 81 
Plotting frame 34 of 81 
Plotting frame 35 of 81 
Plotting frame 36 of 81 
Plotting frame 37 of 81 
Plotting frame 38 of 81 
Plotting frame 39 of 81 
Plotting frame 40 of 81 
Plotting frame 41 of 81 
Plotting frame 42 of 81 
Plotting frame 43 of 81 
Plotting frame 44 of 81 
Plotting frame 45 of 81 
Plotting frame 46 of 81 
Plotting frame 47 of 81 
Plotting frame 48 of 81 
Plotting frame 49 of 81 
Plotting frame 50 of 81 
Plotting frame 51 of 81 
Plotting frame 52 of 81 
Plotting frame 53 of 81 
Plotting frame 54 of 81 
Plotting frame 55 of 81 
Plotting frame 56 of 81 
Plotting frame 57 of 81 
Plotting frame 58 of 81 
Plotting frame 59 of 81 
Plotting frame 60 of 81 
Plotting frame 61 of 81 
Plotting frame 62 of 81 
Plotting frame 63 of 81 
Plotting frame 64 of 81 
Plotting frame 65 of 81 
Plotting frame 66 of 81 
Plotting frame 67 of 81 
Plotting frame 68 of 81 
Plotting frame 69 of 81 
Plotting frame 70 of 81 
Plotting frame 71 of 81 
Plotting frame 72 of 81 
Plotting frame 73 of 81 
Plotting frame 74 of 81 
Plotting frame 75 of 81 
Plotting frame 76 of 81 
Plotting frame 77 of 81 
Plotting frame 78 of 81 
Plotting frame 79 of 81 
Plotting frame 80 of 81 
Plotting frame 81 of 81


This page was generated using Literate.jl.