Oceananigans.jl

Appearance

Hydrostatic lock exchange with CATKEVerticalDiffusivity ​

This example simulates a lock exchange problem on a slope. It demonstrates:

  • How to set up a 2D grid with a sloping bottom with an immersed boundary

  • Initializing a hydrostatic free surface model

  • Including variable density initial conditions

  • Applying bottom drag boundary conditions on immersed boundaries using BulkDrag

  • Saving outputs of the simulation

  • Creating an animation with CairoMakie

The Lock Exchange Problem ​

In a lock exchange, two fluids of different densities (due to temperature, salinity, etc.) are separated by a "lock" at time  . Once the lock is removed, the fluids interact and form gravity currents. The lock exchange represents scenarios where waters of different salinities or temperatures meet and form sharp density gradients, for example in estuaries or in the Denmark Strait overflow. Its evolution can be described by the hydrostatic Boussinesq equations. In this example, we use buoyancy as a tracer; see the Boussinesq approximation section for more details on the Boussinesq approximation.

Install dependencies ​

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

julia
using Pkg
pkg"add Oceananigans, CairoMakie"

Import Required Packages ​

julia
using Oceananigans
using Oceananigans.Units
using Printf
using CairoMakie

Set up a 2D rectilinear grid ​

julia
Nx = 128
Nz = 64
L = 8kilometers   # horizontal length
H = 50meters      # depth
50.0

We wrap z into a MutableVerticalDiscretization. This allows the HydrostaticFreeSurfaceModel (which we construct further down) to use a time-evolving, free-surface-following ZStarCoordinate vertical coordinate.

julia
x = (-L/8, 7L/8)
z = MutableVerticalDiscretization((-H, 0))

underlying_grid = RectilinearGrid(; size=(Nx, Nz), x, z, halo=(5, 5),
                                    topology=(Bounded, Flat, Bounded))
128×1×64 RectilinearGrid{Float32, Bounded, Flat, Bounded} on CPU with 5×0×5 halo
├── Bounded  x ∈ [-1000.0, 7000.0] regularly spaced with Δx=62.5
├── Flat y                         
└── Bounded  z ∈ [-50.0, 0.0]      regularly spaced with mutable Δr=0.78125

The bottom slopes upward from    at   to    at  . We describe the sloping bottom with a partial-cell immersed boundary.

julia
slope = H / 2L
bottom(x) = -H + slope * x

grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(bottom))
128×1×64 ImmersedBoundaryGrid{Float32, Bounded, Flat, Bounded} on CPU with 5×0×5 halo:
├── immersed_boundary: PartialCellBottom(mean(zb)=-40.4425, min(zb)=-50.0, max(zb)=-28.2812, ϵ=0.2)
├── underlying_grid: 128×1×64 RectilinearGrid{Float32, Bounded, Flat, Bounded} on CPU with 5×0×5 halo
├── Bounded  x ∈ [-1000.0, 7000.0] regularly spaced with Δx=62.5
├── Flat y                         
└── Bounded  z ∈ [-50.0, 0.0]      regularly spaced with mutable Δr=0.78125

Set up a bottom drag boundary condition on the immersed boundary ​

To apply drag to the flow, we use BulkDrag, which implements quadratic drag proportional to ᴰ, where ᴰ is the drag coefficient and   is the horizontal speed. We use a drag coefficient ᴰ , a reasonable value for seafloor drag.

julia
Cᴰ = 0.002
drag = BulkDrag(coefficient=Cᴰ)
FluxBoundaryCondition: BulkDragFunction(QuadraticFormulation(), Nothing, Cᴰ=0.002)

BulkDrag can be applied both to domain boundaries (like the bottom of a Bounded grid) and to immersed boundaries. Here we apply it to the immersed sloping bottom boundary:

julia
u_bcs = FieldBoundaryConditions(bottom=drag, immersed=drag)
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── bottom: FluxBoundaryCondition: BulkDragFunction(QuadraticFormulation(), Nothing, Cᴰ=0.002)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
└── immersed: FluxBoundaryCondition: BulkDragFunction(QuadraticFormulation(), Nothing, Cᴰ=0.002)

In a 2D simulation with Flat in the -direction, we don't need boundary conditions for .

Initialize the model ​

  • We use a hydrostatic model because the horizontal scale of the flow is much larger than its vertical scale

  • Buoyancy is a tracer that determines the fluid density

  • The vertical closure CATKEVerticalDiffusivity parameterizes small-scale vertical turbulence

  • Weighted Essentially Non-Oscillatory (WENO) advection schemes capture sharp changes in density

julia
model = HydrostaticFreeSurfaceModel(grid;
                                    tracers = :b,
                                    buoyancy = BuoyancyTracer(),
                                    closure = CATKEVerticalDiffusivity(),
                                    momentum_advection = WENO(order=5),
                                    tracer_advection = WENO(order=7),
                                    boundary_conditions = (; u=u_bcs),
                                    free_surface = SplitExplicitFreeSurface(grid; substeps=20))
HydrostaticFreeSurfaceModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×64 ImmersedBoundaryGrid{Float32, Bounded, Flat, Bounded} on CPU with 5×0×5 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: (b, e)
├── closure: CATKEVerticalDiffusivity{VerticallyImplicitTimeDiscretization}
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
├── free surface: SplitExplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── substepping: FixedSubstepNumber(14)
├── advection scheme: 
│   ├── momentum: WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── b: WENO{4, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=7)
│   └── e: WENO{4, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=7)
├── vertical_coordinate: ZStarCoordinate
└── coriolis: Nothing

Set variable density initial conditions ​

The lock at   separates light fluid on the left from denser fluid on the right.

julia
bᵢ(x, z) = x > L/2 ? 0.01 : 0.06
set!(model, b=bᵢ)

Construct a Simulation ​

Fast wave speeds make the equations stiff, so the CFL condition restricts the time step to small values to maintain numerical stability.

julia
simulation = Simulation(model, Δt=1second, stop_time=6hours)
Simulation of HydrostaticFreeSurfaceModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 1 second
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 6 hours
├── 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

The TimeStepWizard is incorporated in the simulation via the conjure_time_step_wizard! helper function and it ensures stable time-stepping with a Courant–Friedrichs–Lewy (CFL) number of 0.3.

julia
conjure_time_step_wizard!(simulation, cfl=0.3)

Track simulation progress ​

We add a callback that prints the simulation progress alongside some flow statistics.

julia
function progress(sim)
    @info @sprintf("Iter: %6d, time: %s, Δt: %s, wall time: %s, max|w| = %6.3e m s⁻¹",
                   iteration(sim),
                   prettytime(sim),
                   prettytime(sim.Δt),
                   prettytime(sim.run_wall_time),
                   maximum(abs, sim.model.velocities.w))
    return nothing
end

add_callback!(simulation, progress, IterationInterval(1000))

Add output writer ​

Here, we construct a JLD2Writer to save buoyancy , turbulent kinetic energy , horizontal velocity , and stratification    every 2 minutes.

julia
b = model.tracers.b
e = model.tracers.e
u = model.velocities.u
N² = ∂z(b)

filename = "hydrostatic_lock_exchange.jld2"
simulation.output_writers[:fields] = JLD2Writer(model, (; b, e, u, N²);
                                                filename,
                                                schedule = TimeInterval(2minutes),
                                                overwrite_files = true)
JLD2Writer scheduled on TimeInterval(2 minutes):
├── filepath: hydrostatic_lock_exchange.jld2
├── 4 outputs: (b, e, u, N²)
├── array_type: Array{Float32}
├── including: [:coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

Run the simulation ​

julia
run!(simulation)

@info "Simulation finished. Output saved to $(filename)"
[ Info: Initializing simulation...
[ Info: Iter:      0, time: 0 seconds, Δt: 1.100 seconds, wall time: 0 seconds, max|w| = 0.000e+00 m s⁻¹
[ Info:     ... simulation initialization complete (8.767 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (17.744 seconds).
[ Info: Iter:   1000, time: 32.353 minutes, Δt: 1.716 seconds, wall time: 34.243 seconds, max|w| = 1.356e-01 m s⁻¹
[ Info: Iter:   2000, time: 1.030 hours, Δt: 1.726 seconds, wall time: 42.286 seconds, max|w| = 1.355e-01 m s⁻¹
[ Info: Iter:   3000, time: 1.480 hours, Δt: 1.837 seconds, wall time: 51.552 seconds, max|w| = 1.261e-01 m s⁻¹
[ Info: Iter:   4000, time: 1.900 hours, Δt: 1.430 seconds, wall time: 1.022 minutes, max|w| = 1.636e-01 m s⁻¹
[ Info: Iter:   5000, time: 2.300 hours, Δt: 832.665 ms, wall time: 1.187 minutes, max|w| = 2.813e-01 m s⁻¹
[ Info: Iter:   6000, time: 3.005 hours, Δt: 3.513 seconds, wall time: 1.358 minutes, max|w| = 6.530e-02 m s⁻¹
[ Info: Iter:   7000, time: 3.917 hours, Δt: 2.231 seconds, wall time: 1.532 minutes, max|w| = 1.028e-01 m s⁻¹
[ Info: Iter:   8000, time: 4.420 hours, Δt: 1.627 seconds, wall time: 1.695 minutes, max|w| = 1.429e-01 m s⁻¹
[ Info: Iter:   9000, time: 4.906 hours, Δt: 2.188 seconds, wall time: 1.860 minutes, max|w| = 1.069e-01 m s⁻¹
[ Info: Iter:  10000, time: 5.556 hours, Δt: 2.664 seconds, wall time: 2.005 minutes, max|w| = 8.704e-02 m s⁻¹
[ Info: Simulation is stopping after running for 2.092 minutes.
[ Info: Simulation time 6 hours equals or exceeds stop time 6 hours.
[ Info: Simulation finished. Output saved to hydrostatic_lock_exchange.jld2

Load the saved time series ​

julia
u_timeseries = FieldTimeSeries(filename, "u")
N²_timeseries = FieldTimeSeries(filename, "N²")
b_timeseries = FieldTimeSeries(filename, "b")
e_timeseries = FieldTimeSeries(filename, "e")
times = b_timeseries.times
0.0:120.0:21600.0

Visualize the simulation output ​

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

julia
n = Observable(1)

title = @lift "t = " * prettytime(times[$n])

uₙ = @lift u_timeseries[$n]
N²ₙ = @lift N²_timeseries[$n]
bₙ = @lift b_timeseries[$n]
eₙ = @lift e_timeseries[$n]

We use the last snapshot to set the color ranges.

julia
umax = maximum(abs, u_timeseries[end])
N²max = maximum(abs, N²_timeseries[end])
bmax = maximum(abs, b_timeseries[end])
emax = maximum(abs, e_timeseries[end])

We visualize , , , and  .

julia
nan_color = :grey
axis_kwargs = (xlabel = "x (m)", ylabel = "z (m)",
               limits = (x, (-H, 0)), titlesize = 18)

fig = Figure(size = (800, 900))
fig[1, :] = Label(fig, title, fontsize = 24, tellwidth = false)

ax_b = Axis(fig[2, 1]; title = "b (buoyancy)", axis_kwargs...)
hm_b = heatmap!(ax_b, bₙ; nan_color, colorrange = (0, bmax), colormap = :thermal)
Colorbar(fig[2, 2], hm_b, label = "m s⁻²")

ax_e = Axis(fig[3, 1]; title = "e (turbulent kinetic energy)", axis_kwargs...)
hm_e = heatmap!(ax_e, eₙ; nan_color, colorrange = (0, emax), colormap = :magma)
Colorbar(fig[3, 2], hm_e, label = "m² s⁻²")

ax_u = Axis(fig[4, 1]; title = "u (horizontal velocity)", axis_kwargs...)
hm_u = heatmap!(ax_u, uₙ; nan_color, colorrange = (-umax, umax), colormap = :balance)
Colorbar(fig[4, 2], hm_u, label = "m s⁻¹")

ax_N² = Axis(fig[5, 1]; title = "N² (stratification)", axis_kwargs...)
hm_N² = heatmap!(ax_N², N²ₙ; nan_color, colorrange = (-N²max/4, N²max), colormap = :haline)
Colorbar(fig[5, 2], hm_N², label = "s⁻²")

record(fig, "hydrostatic_lock_exchange.mp4", 1:length(times); framerate = 8) do i
    n[] = i
end

The visualization shows the time evolution of buoyancy , turbulent kinetic energy , horizontal velocity , and stratification  . Initially, the two water masses are separated horizontally by a sharp density interface. As the flow evolves, gravity currents form and the dense fluid moves beneath the lighter fluid. This shows the characteristic transition from horizontal to vertical density separation for the lock exchange problem.

Julia version and environment information ​

This example was executed with the following version of Julia:

julia
using InteractiveUtils: versioninfo
versioninfo()
Julia Version 1.13.0
Commit d1c37793dd2 (2026-09-09 19:00 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-20.1.8 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
Environment:
  JULIA_LOAD_PATH = @:@v#.#:@stdlib
  JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.0+0.x64.linux.gnu/local/share/julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.0+0.x64.linux.gnu/share/julia
  JULIA_VERSION_ENZYME = 1.13.0
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  LD_LIBRARY_PATH = 
  JULIA_MAX_NUM_PRECOMPILE_FILES = 24
  JULIA_VERSION = 1.13.0
  JULIA_CUDA_USE_COMPAT = false
  JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-34931/docs/
  JULIA_DEBUG = Literate

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

julia
import Pkg
Pkg.status()
Status `~/Oceananigans.jl-34931/docs/Project.toml`
  [79e6a3ab] Adapt v4.7.4
⌃ [052768ef] CUDA v6.2.2
  [13f3f980] CairoMakie v0.15.16
  [8e50ac2c] ConservativeRegridding v0.2.16
⌅ [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.5.0
  [4710194d] DocumenterVitepress v0.3.8
  [7da242da] Enzyme v0.13.216
  [033835bb] JLD2 v0.6.7
⌅ [63c18a36] KernelAbstractions v0.9.45
  [98b081ad] Literate v2.21.0
  [da04e1cc] MPI v0.20.29
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.114.0 `..`
  [429524aa] Optim v2.3.2
  [f27b6e38] Polynomials v4.1.3
⌅ [3c362404] Reactant v0.2.292
  [6038ab10] Rotations v1.7.1
  [d496a93d] SeawaterPolynomials v0.3.11
  [09ab397b] StructArrays v0.7.4
  [bdfc003b] TimesDates v0.3.3
  [0a941bbe] Zarr v0.11.2
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.13.0
  [44cfe95a] Pkg v1.13.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.