Tilted bottom boundary layer example
This example is based on the similar Oceananigans example and 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 constant density stratification. This perturbation develops into a turbulent bottom boundary layer due to momentum loss at the bottom boundary.
First let's make sure we have all required packages installed.
using Pkg
pkg"add Oceananigans, Oceanostics, Rasters, CairoMakie"Grid
We start by creating a $x, z$ grid with 64² cells and finer resolution near the bottom:
using Oceananigans
using Oceananigans.Units
Lx = 200meters
Lz = 100meters
Nx = 64
Nz = 64
refinement = 1.8 # controls spacing near surface (higher means finer spaced)
stretching = 10 # controls rate of stretching at 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.55496Note that, with the z faces defined as above, the spacings near the bottom are approximately constant, becoming progressively coarser moving up.
Tilting the domain
We use a domain that's tilted with respect to gravity by
θ = 5; # degreesso that $x$ is the along-slope direction, $z$ is the across-sloce direction that is perpendicular to the bottom, and the unit vector anti-aligned with gravity is
ĝ = [sind(θ), 0, cosd(θ)]3-element Vector{Float64}:
0.08715574274765818
0.0
0.9961946980917455Changing the vertical direction impacts both the gravity_unit_vector for Buoyancy as well as the rotation_axis for Coriolis forces,
buoyancy = BuoyancyForce(BuoyancyTracer(), gravity_unit_vector = -ĝ)
f₀ = 1e-4/second
coriolis = ConstantCartesianCoriolis(f = f₀, rotation_axis = ĝ)ConstantCartesianCoriolis{Float64}: fx = 8.72e-06, fy = 0.00e+00, fz = 9.96e-05The tilting also affects the kind of density stratified flows we can model. The simulate an environment that's uniformly stratified, with a stratification frequency
N² = 1e-5/second^2;In a tilted coordinate, this can be achieved with
@inline constant_stratification(x, z, t, p) = p.N² * (x * p.ĝ[1] + z * p.ĝ[3]);However, this distribution is not periodic in $x$ and can't be explicitly modelled on an $x$-periodic grid such as the one used here. Instead, we simulate periodic perturbations away from the constant density stratification by imposing a constant stratification as a BackgroundField,
B_field = BackgroundField(constant_stratification, parameters=(; ĝ, N²))BackgroundField{typeof(Main.constant_stratification), @NamedTuple{ĝ::Vector{Float64}, N²::Float64}}
├── func: constant_stratification (generic function with 1 method)
└── parameters: (ĝ = [0.08715574274765818, 0.0, 0.9961946980917455], N² = 1.0e-5)Bottom drag
We impose bottom drag that follows Monin-Obukhov theory and include the background flow in the drag calculation, which is the only effect the background flow has on the problem
V∞ = 0.1meters/second
z₀ = 0.1meters # (roughness length)
κ = 0.4 # von Karman constant
z₁ = znodes(grid, Center())[1] # Closest grid center to the bottom
cᴰ = (κ / log(z₁ / z₀))^2 # Drag coefficient
@inline drag_u(x, t, u, v, p) = - p.cᴰ * √(u^2 + (v + p.V∞)^2) * u
@inline drag_v(x, t, u, v, p) = - p.cᴰ * √(u^2 + (v + p.V∞)^2) * (v + p.V∞)
drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=(; cᴰ, V∞))
drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=(; cᴰ, V∞))
u_bcs = FieldBoundaryConditions(bottom = drag_bc_u)
v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: FluxBoundaryCondition: ContinuousBoundaryFunction drag_v at (Nothing, Nothing, Nothing)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)Create model and simulation
We are now ready to create the model. We create a NonhydrostaticModel with a 5th UpwindBiased advection scheme, a RungeKutta3 timestepper, and a constant viscosity and diffusivity.
closure = ScalarDiffusivity(ν=2e-4, κ=2e-4)
model = NonhydrostaticModel(; grid, buoyancy, coriolis, closure,
timestepper = :RungeKutta3,
advection = UpwindBiased(order=5),
tracers = :b,
boundary_conditions = (u=u_bcs, v=v_bcs),
background_fields = (; b=B_field))
noise(x, z) = 1e-3 * randn() * exp(-(10z)^2/grid.Lz^2)
set!(model, u=noise, w=noise)The bottom-intensified noise above should accelerate the emergence of turbulence close to the wall.
We are now ready to create the simulation. We begin by setting the initial time step conservatively, based on the smallest grid size of our domain and set-up a
using Oceananigans.Units
simulation = Simulation(model, Δt = 0.5 * minimum_zspacing(grid) / V∞, stop_time = 12hours)Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 4.344 seconds
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: 12 hours
├── Stop iteration: Inf
├── Wall time limit: Inf
├── Minimum relative step: 0.0
├── Callbacks: OrderedDict with 4 entries:
│ ├── stop_time_exceeded => 4
│ ├── stop_iteration_exceeded => -
│ ├── wall_time_limit_exceeded => e
│ └── nan_checker => }
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entriesWe use TimeStepWizard to maximize Δt
wizard = TimeStepWizard(max_change=1.1, cfl=0.7)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(4))Callback of TimeStepWizard(cfl=0.7, max_Δt=Inf, min_Δt=0.0) on IterationInterval(4)Model diagnostics
We set-up a custom progress messenger using Oceanostics.ProgressMessengers, which allows us to combine different ProgressMessengers into one:
using Oceanostics.ProgressMessengers
walltime_per_timestep = StepDuration() # This needs to instantiated here, and not in the function below
progress(simulation) = @info (PercentageProgress(with_prefix=false, with_units=false) + SimulationTime() + TimeStep() + MaxVelocities() + AdvectiveCFLNumber() + walltime_per_timestep)(simulation)
simulation.callbacks[:progress] = Callback(progress, IterationInterval(400))Callback of progress on IterationInterval(400)We now define some useful diagnostics for the flow. Namely, we define RichardsonNumber, RossbyNumber and ErtelPotentialVorticity:
using Oceanostics
b = model.tracers.b + model.background_fields.tracers.b
Ri = RichardsonNumber(model, model.velocities..., b)
Ro = RossbyNumber(model)
PV = ErtelPotentialVorticity(model, model.velocities..., b, model.coriolis)KernelFunctionOperation at (Face, Face, Face)
├── grid: 64×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── kernel_function: ertel_potential_vorticity_fff (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Oceananigans.AbstractOperations.BinaryOperation", "Float64", "Float64", "Float64")Note that the calculation of these quantities depends on the alignment with the true (geophysical) vertical and the rotation axis. Oceanostics already takes that into consideration by using model.buoyancy and model.coriolis, making their calculation much easier. Furthermore, passing the flag add_background=true automatically adds the model's BackgroundFields to the resolved perturbations, which is important in our case for the correct calculation of $\nabla b$ with the background stratification.
Now we write these quantities to a NetCDF file:
output_fields = (; Ri, Ro, PV, b)
using NCDatasets
filename = "tilted_bottom_boundary_layer"
simulation.output_writers[:nc] = NetCDFWriter(model, output_fields,
filename = joinpath(@__DIR__, filename),
schedule = TimeInterval(20minutes),
overwrite_existing = true)NetCDFWriter scheduled on TimeInterval(20 minutes):
├── filepath: build/generated/tilted_bottom_boundary_layer.nc
├── dimensions: time(0), x_faa(64), x_caa(64), z_aaf(65), z_aac(64)
├── 4 outputs: (Ri, b, PV, Ro)
├── array_type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: 22.4 KiBRun the simulation and process results
To run the simulation:
run!(simulation)[ Info: Initializing simulation...
[ Info: 000.00%, time = 0 seconds, Δt = 4.778 seconds, |u⃗|ₘₐₓ = [2.40e-03, 0.00e+00, 1.11e-03] m/s, advective CFL = 0.0084, walltime / timestep = 0 seconds
[ Info: ... simulation initialization complete (9.552 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (3.966 seconds).
[ Info: 086.11%, time = 10.333 hours, Δt = 1.396 minutes, |u⃗|ₘₐₓ = [1.64e-02, 9.35e-02, 5.52e-03] m/s, advective CFL = 0.7, walltime / timestep = 23.242 ms
[ Info: Simulation is stopping after running for 16.564 seconds.
[ Info: Simulation time 12 hours equals or exceeds stop time 12 hours.Now we'll read the results and plot an animation
ds = NCDataset(simulation.output_writers[:nc].filepath)Dataset: /home/runner/work/Oceanostics.jl/Oceanostics.jl/docs/build/generated/tilted_bottom_boundary_layer.nc
Group: /
Dimensions
time = 37
x_faa = 64
x_caa = 64
z_aaf = 65
z_aac = 64
Variables
time (37)
Datatype: Float64 (Float64)
Dimensions: time
Attributes:
units = seconds
long_name = Time
x_faa (64)
Datatype: Float32 (Float32)
Dimensions: x_faa
Attributes:
units = m
long_name = Cell face locations in the x-direction.
x_caa (64)
Datatype: Float32 (Float32)
Dimensions: x_caa
Attributes:
units = m
long_name = Cell center locations in the x-direction.
z_aaf (65)
Datatype: Float32 (Float32)
Dimensions: z_aaf
Attributes:
units = m
long_name = Cell face locations in the z-direction.
z_aac (64)
Datatype: Float32 (Float32)
Dimensions: z_aac
Attributes:
units = m
long_name = Cell center locations in the z-direction.
Δx_caa (64)
Datatype: Float32 (Float32)
Dimensions: x_caa
Attributes:
units = m
long_name = Spacings between cell faces (located at the cell centers) in the x-direction.
Δx_faa (64)
Datatype: Float32 (Float32)
Dimensions: x_faa
Attributes:
units = m
long_name = Spacings between cell centers (located at the cell faces) in the x-direction.
Δz_aac (64)
Datatype: Float32 (Float32)
Dimensions: z_aac
Attributes:
units = m
long_name = Spacings between cell faces (located at cell centers) in the z-direction.
Δz_aaf (65)
Datatype: Float32 (Float32)
Dimensions: z_aaf
Attributes:
units = m
long_name = Spacings between cell centers (located at cell faces) in the z-direction.
PV (64 × 65 × 37)
Datatype: Float32 (Float32)
Dimensions: x_faa × z_aaf × time
Ri (64 × 65 × 37)
Datatype: Float32 (Float32)
Dimensions: x_caa × z_aaf × time
Ro (64 × 65 × 37)
Datatype: Float32 (Float32)
Dimensions: x_faa × z_aaf × time
b (64 × 64 × 37)
Datatype: Float32 (Float32)
Dimensions: x_caa × z_aac × time
Attributes:
units = m/s²
long_name = Buoyancy
Global attributes
Julia = This file was generated using
Oceananigans = This file was generated using
date = This file was generated on 2025-07-23T20:59:27.728 local time (2025-07-23T20:59:27.728 UTC).
interval = 1200.0
output time interval = Output was saved every 20 minutes.
schedule = TimeInterval
Groups
Dataset: /home/runner/work/Oceanostics.jl/Oceanostics.jl/docs/build/generated/tilted_bottom_boundary_layer.nc
Group: grid_reconstruction
Dimensions
x_f = 65
z_f = 65
Variables
x_f (65)
Datatype: Float32 (Float32)
Dimensions: x_f
z_f (65)
Datatype: Float32 (Float32)
Dimensions: z_f
Global attributes
Hx = 3
Hy = 0
Hz = 3
Nx = 64
Ny = 1
Nz = 64
TX = Periodic
TY = Flat
TZ = Bounded
eltype = Float64
type = RectilinearGrid
x_spacing = regular
y_spacing = flat
z_spacing = irregular
We now use Makie to create the figure and its axes
using CairoMakie
set_theme!(Theme(fontsize = 20))
fig = Figure()
kwargs = (xlabel="x [m]", ylabel="z [m]", height=150, width=250)
ax1 = Axis(fig[2, 1]; title = "Ri", kwargs...)
ax2 = Axis(fig[2, 2]; title = "Ro", kwargs...)
ax3 = Axis(fig[2, 3]; title = "PV", kwargs...);Next we an Observable to lift the values at each specific time and plot heatmaps, along with their colorbars, with buoyancy contours on top
n = Observable(1)
x_caa = ds["x_caa"][:]
x_faa = ds["x_faa"][:]
z_aac = ds["z_aac"][:]
z_aaf = ds["z_aaf"][:]
bₙ = @lift ds["b"][:, :, $n]
Riₙ = @lift ds["Ri"][:, :, $n]
hm1 = heatmap!(ax1, x_caa, z_aaf, Riₙ; colormap = :coolwarm, colorrange = (-1, +1))
contour!(ax1, x_caa, z_aac, bₙ; levels=10, color=:white, linestyle=:dash, linewidth=0.5)
Colorbar(fig[3, 1], hm1, vertical=false, height=8, ticklabelsize=14)
Roₙ = @lift ds["Ro"][:, :, $n]
hm2 = heatmap!(ax2, x_faa, z_aaf, Roₙ; colormap = :balance, colorrange = (-10, +10))
contour!(ax2, x_caa, z_aac, bₙ; levels=10, color=:black, linestyle=:dash, linewidth=0.5)
Colorbar(fig[3, 2], hm2, vertical=false, height=8, ticklabelsize=14)
PVₙ = @lift ds["PV"][:, :, $n]
hm3 = heatmap!(ax3, x_faa, z_aaf, PVₙ; colormap = :coolwarm, colorrange = N²*f₀.*(-1.5, +1.5))
contour!(ax3, x_caa, z_aac, bₙ; levels=10, color=:white, linestyle=:dash, linewidth=0.5)
Colorbar(fig[3, 3], hm3, vertical=false, height=8, ticklabelsize=14);Now we mark the time by placing a vertical line in the bottom panel and adding a helpful title
times = ds["time"][:]
title = @lift "Time = " * string(prettytime(times[$n]))
fig[1, 1:3] = Label(fig, title, fontsize=24, tellwidth=false);Finally, we adjust the figure dimensions to fit all the panels and record a movie
resize_to_layout!(fig)
@info "Animating..."
record(fig, filename * ".mp4", 1:length(times), framerate=10) do i
n[] = i
end
close(ds)closed DatasetThe animation shows negative PV being produced at the bottom due to drag, which leads to the emergence of centrifulgal-symmetric instabilities, which become turbulent and erode stratification (as can be seen by inspecting $Ri$). Note that there are some boundary effects on the upper boundary, likely caused by interaction internal waves that are produced by the bottom turbulence. These effects are, to some degree, expected, and a sponge/relaxation layer at the top is needed to minimize them in a production-ready code.
This page was generated using Literate.jl.