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.55496

Note 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; # degrees

so 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

ĝ = [sind(θ), 0, cosd(θ)]

3-element Vector{Float64}:
 0.08715574274765818
 0.0
 0.9961946980917455

Changing 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 = -ĝ)

f₀ = 1e-4/second
coriolis = ConstantCartesianCoriolis(f = f₀, rotation_axis = ĝ)

ConstantCartesianCoriolis{Float64}: fx = 8.72e-06, fy = 0.00e+00, fz = 9.96e-05

The 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.ĝ[1] + z * p.ĝ[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=(; ĝ, 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 => 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 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: ("64×1×64 Field{Face, Center, Center} on RectilinearGrid on CPU", "64×1×64 Field{Center, Face, Center} on RectilinearGrid on CPU", "64×1×65 Field{Center, Center, Face} on RectilinearGrid on CPU", "BinaryOperation at (Center, Center, Center)", "8.71557e-6", "0.0", "9.96195e-5")

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.8 KiB

Run 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.05e-03,  0.00e+00,  1.36e-03] m/s,  advective CFL = 0.0087,  walltime / timestep = 0 seconds
[ Info:     ... simulation initialization complete (10.772 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.292 seconds).
[ Info: 082.86%,  time = 9.943 hours,  Δt = 1.069 minutes,  |u⃗|ₘₐₓ = [1.42e-02,  9.32e-02,  8.73e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 25.576 ms
[ Info: Simulation is stopping after running for 0 seconds.
[ Info: Simulation time 12 hours equals or exceeds stop time 12 hours.

Now we'll read the results and plot an animation

using Rasters

ds = RasterStack(simulation.output_writers[:nc].filepath)

┌ 64×64×65×64×37 RasterStack ┐
├────────────────────────────┴─────────────────────────────────────────── dims ┐
  ↓ x_faa Sampled{Float32} [0.0, 3.125, …, 193.75, 196.875] ForwardOrdered Regular Points,
  → x_caa Sampled{Float32} [1.5625, 4.6875, …, 195.3125, 198.4375] ForwardOrdered Regular Points,
  ↗ z_aaf Sampled{Float32} [-0.0, 0.86881673, …, 93.44504, 100.0] ForwardOrdered Irregular Points,
  ⬔ z_aac Sampled{Float32} [0.43440837, 1.3032823, …, 90.52517, 96.72252] ForwardOrdered Irregular Points,
  ◩ Ti    Sampled{Float64} [0.0, 1200.0, …, 42000.0, 43200.0] ForwardOrdered Regular Points
├────────────────────────────────────────────────────────────────────── layers ┤
  :Δx_caa eltype: Union{Missing, Float32} dims: x_caa size: 64
  :Δx_faa eltype: Union{Missing, Float32} dims: x_faa size: 64
  :Δz_aac eltype: Union{Missing, Float32} dims: z_aac size: 64
  :Δz_aaf eltype: Union{Missing, Float32} dims: z_aaf size: 65
  :PV     eltype: Union{Missing, Float32} dims: x_faa, z_aaf, Ti size: 64×65×37
  :Ri     eltype: Union{Missing, Float32} dims: x_caa, z_aaf, Ti size: 64×65×37
  :Ro     eltype: Union{Missing, Float32} dims: x_faa, z_aaf, Ti size: 64×65×37
  :b      eltype: Union{Missing, Float32} dims: x_caa, z_aac, Ti size: 64×64×37
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.NCDsource} of Dict{String, Any} with 6 entries:
  "interval"             => 1200.0
  "Oceananigans"         => "This file was generated using "
  "Julia"                => "This file was generated using Julia Version 1.10.9…
  "output time interval" => "Output was saved every 20 minutes."
  "date"                 => "This file was generated on 2025-03-19T23:10:38.676…
  "schedule"             => "TimeInterval"
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(x_faa = (0.0f0, 196.875f0), x_caa = (1.5625f0, 198.4375f0), z_aaf = (-0.0f0, 100.0f0), z_aac = (0.43440837f0, 96.72252f0), Ti = (0.0, 43200.0))
└──────────────────────────────────────────────────────────────────────────────┘

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 = Array(dims(ds, :x_caa))
x_faa = Array(dims(ds, :x_faa))
z_aac = Array(dims(ds, :z_aac))
z_aaf = Array(dims(ds, :z_aaf))

bₙ = @lift Array(ds.b[Ti=$n, y_aca=Near(0)])

Riₙ = @lift Array(ds.Ri[Ti=$n, y_aca=Near(0)])
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 Array(ds.Ro[Ti=$n, y_afa=Near(0)])
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 Array(ds.PV[Ti=$n, y_afa=Near(0)])
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);

┌ Warning: (DimensionalData.Dimensions.Dim{:y_aca},) dims were not found in object.
└ @ DimensionalData.Dimensions ~/.julia/packages/DimensionalData/a0k29/src/Dimensions/primitives.jl:849
┌ Warning: (DimensionalData.Dimensions.Dim{:y_aca},) dims were not found in object.
└ @ DimensionalData.Dimensions ~/.julia/packages/DimensionalData/a0k29/src/Dimensions/primitives.jl:849
┌ Warning: (DimensionalData.Dimensions.Dim{:y_afa},) dims were not found in object.
└ @ DimensionalData.Dimensions ~/.julia/packages/DimensionalData/a0k29/src/Dimensions/primitives.jl:849
┌ Warning: (DimensionalData.Dimensions.Dim{:y_afa},) dims were not found in object.
└ @ DimensionalData.Dimensions ~/.julia/packages/DimensionalData/a0k29/src/Dimensions/primitives.jl:849

Now we mark the time by placing a vertical line in the bottom panel and adding a helpful title

times = dims(ds, :Ti)
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

"tilted_bottom_boundary_layer.mp4"

The 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.