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

Turbulence closures

We prescribe the values of vertical viscocity and diffusivity according to the ratio of the vertical and lateral grid spacing.

Δx, Δz = Lx/Nx, Lz/Nz
𝒜 = Δz/Δx # Grid cell aspect ratio.

κh = 0.1    # [m² s⁻¹] horizontal diffusivity
νh = 0.1    # [m² s⁻¹] horizontal viscosity
κz = 𝒜 * κh # [m² s⁻¹] vertical diffusivity
νz = 𝒜 * νh # [m² s⁻¹] vertical viscosity

horizontal_diffusive_closure = HorizontalScalarDiffusivity(ν = νh, κ = κh)

vertical_diffusive_closure = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization();
                                                       ν = νz, κ = κz)

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,
                                    closure = (vertical_diffusive_closure, horizontal_diffusive_closure),
                                    momentum_advection = WENO(),
                                    tracer_advection = WENO(),
                                    free_surface = ImplicitFreeSurface())
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: Tuple with 2 closures:
│   ├── VerticalScalarDiffusivity{VerticallyImplicitTimeDiscretization}(ν=0.00016, κ=(b=0.00016,))
│   └── HorizontalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.1, κ=(b=0.1,))
├── 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 siulation'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$) averages of buoyancy and zonal velocity $u$.

u, v, w = model.velocities

B = Field(Average(b, dims=1))
U = Field(Average(u, 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),
                                                       indices)
end

simulation.output_writers[:zonal] = JLD2OutputWriter(model, (b=B,);
                                                     schedule = TimeInterval(save_fields_interval),
                                                     filename = filename * "_zonal_average")
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: 30.055 seconds, max(u): (0.000e+00, 0.000e+00, 0.000e+00) m/s, next Δt: 5.500 minutes
[ Info:     ... simulation initialization complete (14.555 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (2.166 minutes).
[01.17%] i: 100, t: 11.193 hours, wall time: 2.617 minutes, max(u): (6.988e-01, 3.326e-01, 3.363e-03) m/s, next Δt: 8.858 minutes
[03.03%] i: 200, t: 1.213 days, wall time: 18.331 seconds, max(u): (6.061e-01, 3.520e-01, 2.524e-03) m/s, next Δt: 14.266 minutes
[05.99%] i: 300, t: 2.396 days, wall time: 16.667 seconds, max(u): (5.584e-01, 2.879e-01, 2.026e-03) m/s, next Δt: 20 minutes
[09.44%] i: 400, t: 3.778 days, wall time: 16.629 seconds, max(u): (4.729e-01, 2.783e-01, 1.811e-03) m/s, next Δt: 20 minutes
[12.92%] i: 500, t: 5.167 days, wall time: 14.779 seconds, max(u): (4.612e-01, 2.459e-01, 1.641e-03) m/s, next Δt: 20 minutes
[16.39%] i: 600, t: 6.556 days, wall time: 14.745 seconds, max(u): (4.581e-01, 2.598e-01, 1.913e-03) m/s, next Δt: 20 minutes
[19.86%] i: 700, t: 7.944 days, wall time: 14.803 seconds, max(u): (4.867e-01, 4.474e-01, 2.172e-03) m/s, next Δt: 20 minutes
[23.14%] i: 800, t: 9.257 days, wall time: 14.900 seconds, max(u): (5.575e-01, 6.053e-01, 2.478e-03) m/s, next Δt: 18.286 minutes
[26.05%] i: 900, t: 10.419 days, wall time: 14.649 seconds, max(u): (7.626e-01, 7.252e-01, 5.117e-03) m/s, next Δt: 15.446 minutes
[28.84%] i: 1000, t: 11.537 days, wall time: 14.726 seconds, max(u): (7.789e-01, 7.582e-01, 4.354e-03) m/s, next Δt: 19.141 minutes
[31.71%] i: 1100, t: 12.684 days, wall time: 14.630 seconds, max(u): (8.256e-01, 8.570e-01, 4.517e-03) m/s, next Δt: 15.376 minutes
[34.53%] i: 1200, t: 13.811 days, wall time: 14.796 seconds, max(u): (7.795e-01, 8.646e-01, 4.680e-03) m/s, next Δt: 17.807 minutes
[37.53%] i: 1300, t: 15.012 days, wall time: 14.468 seconds, max(u): (8.063e-01, 8.427e-01, 4.690e-03) m/s, next Δt: 17.768 minutes
[40.86%] i: 1400, t: 16.345 days, wall time: 14.861 seconds, max(u): (8.867e-01, 8.449e-01, 4.127e-03) m/s, next Δt: 20 minutes
[44.31%] i: 1500, t: 17.722 days, wall time: 15.507 seconds, max(u): (8.768e-01, 8.053e-01, 3.708e-03) m/s, next Δt: 20 minutes
[47.78%] i: 1600, t: 19.111 days, wall time: 15.103 seconds, max(u): (7.566e-01, 9.701e-01, 3.179e-03) m/s, next Δt: 20 minutes
[51.25%] i: 1700, t: 20.500 days, wall time: 14.714 seconds, max(u): (8.296e-01, 9.522e-01, 2.683e-03) m/s, next Δt: 20 minutes
[54.72%] i: 1800, t: 21.889 days, wall time: 14.746 seconds, max(u): (8.599e-01, 9.753e-01, 2.753e-03) m/s, next Δt: 20 minutes
[58.19%] i: 1900, t: 23.278 days, wall time: 14.710 seconds, max(u): (9.487e-01, 8.012e-01, 2.621e-03) m/s, next Δt: 20 minutes
[61.67%] i: 2000, t: 24.667 days, wall time: 14.513 seconds, max(u): (9.170e-01, 7.861e-01, 2.705e-03) m/s, next Δt: 20 minutes
[65.14%] i: 2100, t: 26.056 days, wall time: 14.508 seconds, max(u): (9.553e-01, 8.331e-01, 2.521e-03) m/s, next Δt: 20 minutes
[68.61%] i: 2200, t: 27.444 days, wall time: 14.091 seconds, max(u): (8.980e-01, 8.257e-01, 2.346e-03) m/s, next Δt: 20 minutes
[72.08%] i: 2300, t: 28.833 days, wall time: 15.478 seconds, max(u): (8.570e-01, 7.846e-01, 2.270e-03) m/s, next Δt: 20 minutes
[75.56%] i: 2400, t: 30.222 days, wall time: 14.534 seconds, max(u): (8.801e-01, 8.854e-01, 2.685e-03) m/s, next Δt: 20 minutes
[79.03%] i: 2500, t: 31.611 days, wall time: 14.133 seconds, max(u): (8.467e-01, 8.814e-01, 2.843e-03) m/s, next Δt: 20 minutes
[82.50%] i: 2600, t: 33 days, wall time: 14.664 seconds, max(u): (9.581e-01, 9.291e-01, 3.393e-03) m/s, next Δt: 20 minutes
[85.97%] i: 2700, t: 34.389 days, wall time: 14.092 seconds, max(u): (9.091e-01, 8.803e-01, 2.717e-03) m/s, next Δt: 20 minutes
[89.44%] i: 2800, t: 35.778 days, wall time: 14.437 seconds, max(u): (8.803e-01, 8.275e-01, 2.423e-03) m/s, next Δt: 20 minutes
[92.92%] i: 2900, t: 37.167 days, wall time: 14.371 seconds, max(u): (8.637e-01, 8.134e-01, 2.062e-03) m/s, next Δt: 20 minutes
[96.39%] i: 3000, t: 38.556 days, wall time: 14.656 seconds, max(u): (8.912e-01, 9.689e-01, 1.788e-03) m/s, next Δt: 20 minutes
[99.86%] i: 3100, t: 39.944 days, wall time: 14.314 seconds, max(u): (7.771e-01, 1.008e+00, 1.803e-03) m/s, next Δt: 20 minutes
[ Info: Simulation is stopping. Model time 40 days has hit or exceeded simulation stop time 40 days.
[ Info: Simulation completed in 10.176 minutes

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{Int64} with 0 listeners. Value:
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; textsize = 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.