Tilted bottom boundary layer
This example simulates a two-dimensional oceanic bottom boundary layer in a domain that's tilted with respect to gravity. We simulate the perturbation away from a constant along-slope (y-direction) velocity and a constant density stratification. This perturbation develops into a turbulent bottom boundary layer due to momentum loss at the bottom boundary modeled with a quadratic drag law.
This example illustrates
changing the direction of gravitational acceleration in the buoyancy model;
changing the axis of rotation for Coriolis forces.
Install dependencies
First let's make sure we have all required packages installed.
using Pkg
pkg"add Oceananigans, NCDatasets, CairoMakie"The domain
We create a grid with finer resolution near the bottom,
using Oceananigans
using Oceananigans.Units
using Random
Random.seed!(42) # for reproducible results
Lx = 200meters
Lz = 100meters
Nx = 64
Nz = 64
refinement = 1.8 # controls spacing near the bottom (higher means finer spaced)
stretching = 10 # controls rate of stretching away from the bottom
h(k) = (Nz + 1 - k) / Nz
ζ(k) = 1 + (h(k) - 1) / refinement
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))
z_faces(k) = - Lz * (ζ(k) * Σ(k) - 1)
grid = RectilinearGrid(topology = (Periodic, Flat, Bounded),
size = (Nx, Nz),
x = (0, Lx),
z = z_faces)64×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── Periodic x ∈ [0.0, 200.0) regularly spaced with Δx=3.125
├── Flat y
└── Bounded z ∈ [-0.0, 100.0] variably spaced with min(Δz)=0.868817, max(Δz)=6.55496Let's make sure the grid spacing is both finer and near-uniform at the bottom,
using CairoMakie
scatterlines(zspacings(grid, Center()),
axis = (ylabel = "Depth (m)",
xlabel = "Vertical spacing (m)"))
Tilting the domain
We use a domain that's tilted with respect to gravity by
θ = 3 # degrees3so that
ẑ = (sind(θ), 0, cosd(θ))(0.052335956242943835, 0, 0.9986295347545738)Changing the vertical direction impacts both the gravity_unit_vector for BuoyancyForce as well as the rotation_axis for Coriolis forces,
buoyancy = BuoyancyForce(BuoyancyTracer(), gravity_unit_vector = .-ẑ)
coriolis = ConstantCartesianCoriolis(f = 1e-4, rotation_axis = ẑ)ConstantCartesianCoriolis{Float64}: fx = 5.23e-06, fy = 0.00e+00, fz = 9.99e-05where above we used a constant Coriolis parameter
@inline constant_stratification(x, z, t, p) = p.N² * (x * p.ẑ[1] + z * p.ẑ[3])constant_stratification (generic function with 1 method)is not periodic in BackgroundField,
N² = 1e-5 # s⁻²
B∞_field = BackgroundField(constant_stratification, parameters=(; ẑ, N²))BackgroundField{typeof(Main.var"##8349".constant_stratification), @NamedTuple{ẑ::Tuple{Float64, Int64, Float64}, N²::Float64}}
├── func: constant_stratification (generic function with 1 method)
└── parameters: (ẑ = (0.052335956242943835, 0, 0.9986295347545738), N² = 1.0e-5)We choose to impose a bottom boundary condition of zero total diffusive buoyancy flux across the seafloor,
This shows that to impose a no-flux boundary condition on the total buoyancy field
b_bcs = FieldBoundaryConditions(bottom = GradientBoundaryCondition(-N² * cosd(θ)))Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── bottom: GradientBoundaryCondition: -9.9863e-6
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)Bottom drag and along-slope interior velocity
We impose bottom drag that follows Monin–Obukhov theory. We use BulkDrag to create the drag boundary conditions, which computes a quadratic drag proportional to the total velocity (including the background velocity):
V∞ = 0.1 # m s⁻¹
ℓ = 0.1 # roughness length (m)
ϰ = 0.4 # von Kármán constant
z₁ = first(znodes(grid, Center())) # height of the grid center closest to the bottom
cᴰ = (ϰ / log(z₁ / ℓ))^2
drag_bc = BulkDrag(coefficient=cᴰ, background_velocities=(0, V∞, 0))
u_bcs = FieldBoundaryConditions(bottom=drag_bc)
v_bcs = FieldBoundaryConditions(bottom=drag_bc)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.07416276026775863)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)Note that, similar to the buoyancy boundary conditions, we had to include the background flow in the drag calculation.
Let us also create a BackgroundField for the along-slope interior velocity:
V∞_field = BackgroundField(V∞)BackgroundField{Float64, Nothing}
├── func: 0.1
└── parameters: nothingCreate the NonhydrostaticModel
We are now ready to create the model. We create a NonhydrostaticModel with a fifth-order UpwindBiased advection scheme and a constant viscosity and diffusivity. Here we use a smallish value of
closure = ScalarDiffusivity(ν=1e-4, κ=1e-4)
model = NonhydrostaticModel(grid; buoyancy, coriolis, closure,
advection = UpwindBiased(order=5),
tracers = :b,
boundary_conditions = (u=u_bcs, v=v_bcs, b=b_bcs),
background_fields = (; b=B∞_field, v=V∞_field))NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 64×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme:
│ ├── momentum: UpwindBiased(order=5)
│ └── b: UpwindBiased(order=5)
├── tracers: b
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0001, κ=(b=0.0001,))
├── buoyancy: BuoyancyTracer with ĝ = (-0.052336, 0.0, -0.99863)
└── coriolis: ConstantCartesianCoriolis{Float64}Let's introduce a bit of random noise at the bottom of the domain to speed up the onset of turbulence:
noise(x, z) = 1e-3 * randn() * exp(-(10z)^2 / grid.Lz^2)
set!(model, u=noise, w=noise)Create and run a simulation
We are now ready to create the simulation. We begin by setting the initial time step conservatively, based on the smallest horizontal grid spacing and the interior velocity.
Δt₀ = 0.5 * minimum_xspacing(grid) / V∞
simulation = Simulation(model, Δt = Δt₀, stop_time = 1day)Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 15.625 seconds
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 1 day
├── 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 entriesWe use a TimeStepWizard to adapt our time-step,
conjure_time_step_wizard!(simulation, IterationInterval(4), max_change=1.1, cfl=0.7)and also we add another callback to print a progress message,
using Printf
progress_message(sim) =
@printf("Iteration: %04d, time: %s, Δt: %s, max|w|: %.1e m s⁻¹, wall time: %s\n",
iteration(sim), prettytime(sim), prettytime(sim.Δt),
maximum(abs, sim.model.velocities.w), prettytime(sim.run_wall_time))
add_callback!(simulation, progress_message, IterationInterval(200))Add outputs to the simulation
We output the total buoyancy NetCDFWriter, which needs NCDatasets to be loaded:
u, v, w = model.velocities
b = model.tracers.b
B∞ = model.background_fields.tracers.b
B = b + B∞
V = v + V∞
ωy = ∂z(u) - ∂x(w)
outputs = (; u, V, w, B, ωy)
using NCDatasets
filename = joinpath(@__DIR__, "tilted_bottom_boundary_layer.nc")
simulation.output_writers[:fields] = NetCDFWriter(model, outputs; filename,
schedule = TimeInterval(20minutes),
overwrite_files = true)NetCDFWriter scheduled on TimeInterval(20 minutes):
├── filepath: tilted_bottom_boundary_layer.nc
├── dimensions: time(0), x_faa(64), x_caa(64), z_aaf(65), z_aac(64)
├── 5 outputs: (u, V, w, B, ωy)
├── array_type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: (file not yet created)Now we just run it!
run!(simulation)[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 17.188 seconds, max|w|: 1.0e-03 m s⁻¹, wall time: 0 seconds
[ Info: ... simulation initialization complete (6.843 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (1.044 seconds).
Iteration: 0200, time: 4.523 hours, Δt: 1.637 minutes, max|w|: 5.2e-03 m s⁻¹, wall time: 8.033 seconds
Iteration: 0400, time: 9.325 hours, Δt: 1.647 minutes, max|w|: 5.1e-03 m s⁻¹, wall time: 8.184 seconds
Iteration: 0600, time: 13.834 hours, Δt: 1.087 minutes, max|w|: 6.3e-03 m s⁻¹, wall time: 8.337 seconds
Iteration: 0800, time: 17.714 hours, Δt: 1.427 minutes, max|w|: 4.5e-03 m s⁻¹, wall time: 8.486 seconds
Iteration: 1000, time: 21.144 hours, Δt: 1.216 minutes, max|w|: 6.6e-03 m s⁻¹, wall time: 8.628 seconds
[ Info: Simulation is stopping after running for 8.739 seconds.
[ Info: Simulation time 1 day equals or exceeds stop time 1 day.Visualize the results
We load the output as FieldTimeSeries and create an animation showing the
fig = Figure(size = (800, 600))
axis_kwargs = (xlabel = "Across-slope distance (m)",
ylabel = "Slope-normal\ndistance (m)",
limits = ((0, Lx), (0, Lz)))
ax_ω = Axis(fig[2, 1]; title = "Along-slope vorticity", axis_kwargs...)
ax_v = Axis(fig[3, 1]; title = "Along-slope velocity (v)", axis_kwargs...)
n = Observable(1)
ωy_timeseries = FieldTimeSeries(filename, "ωy")
B_timeseries = FieldTimeSeries(filename, "B")
V_timeseries = FieldTimeSeries(filename, "V")
ωyn = @lift ωy_timeseries[$n]
Bn = @lift B_timeseries[$n]
Vn = @lift V_timeseries[$n]
buoyancy_levels = -1e-3:5e-5:1e-3
hm_ω = heatmap!(ax_ω, ωyn, colorrange = (-0.015, +0.015), colormap = :balance)
Colorbar(fig[2, 2], hm_ω; label = "s⁻¹")
contour!(ax_ω, Bn, levels=buoyancy_levels, color=:black)
hm_v = heatmap!(ax_v, Vn, colorrange = (-V∞, +V∞), colormap = :balance)
Colorbar(fig[3, 2], hm_v; label = "m s⁻¹")
contour!(ax_v, Bn, levels=buoyancy_levels, color=:black)
times = ωy_timeseries.times
title = @lift "t = " * prettytime(times[$n])
fig[1, :] = Label(fig, title, fontsize=20, tellwidth=false)
fig
Finally, we record a movie.
frames = 1:length(times)
record(fig, "tilted_bottom_boundary_layer.mp4", frames, framerate=12) do i
n[] = i
endJulia version and environment information
This example was executed with the following version of 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-34814/docs/
JULIA_DEBUG = LiterateThese were the top-level packages installed in the environment:
import Pkg
Pkg.status()Status `~/Oceananigans.jl-34814/docs/Project.toml`
[79e6a3ab] Adapt v4.7.4
⌃ [052768ef] CUDA v6.2.2
[13f3f980] CairoMakie v0.15.15
[8e50ac2c] ConservativeRegridding v0.2.16
⌅ [e30172f5] Documenter v1.17.0
[daee34ce] DocumenterCitations v1.5.0
[4710194d] DocumenterVitepress v0.3.7
[7da242da] Enzyme v0.13.216
[033835bb] JLD2 v0.6.7
⌅ [63c18a36] KernelAbstractions v0.9.45
[98b081ad] Literate v2.21.0
[da04e1cc] MPI v0.20.28
[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.3
[bdfc003b] TimesDates v0.3.3
⌃ [0a941bbe] Zarr v0.10.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.