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 = Buoyancy(model = 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, y, 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{(:ĝ, :N²), Tuple{Vector{Float64}, 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, y, t, u, v, p) = - p.cᴰ * √(u^2 + (v + p.V∞)^2) * u
@inline drag_v(x, y, 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 an UpwindBiasedFifthOrder 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 = UpwindBiasedFifthOrder(),
                            tracers = :b,
                            boundary_conditions = (u=u_bcs, v=v_bcs),
                            background_fields = (; b=B_field))

noise(x, y, 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
├── 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

Ri = RichardsonNumber(model, add_background=true)
Ro = RossbyNumber(model)
PV = ErtelPotentialVorticity(model, add_background=true)

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 = model.tracers.b + model.background_fields.tracers.b)

filename = "tilted_bottom_boundary_layer"
simulation.output_writers[:nc] = NetCDFOutputWriter(model, output_fields,
                                                    filename = joinpath(@__DIR__, filename),
                                                    schedule = TimeInterval(20minutes),
                                                    overwrite_existing = true)

NetCDFOutputWriter scheduled on TimeInterval(20 minutes):
├── filepath: /home/runner/work/Oceanostics.jl/Oceanostics.jl/docs/build/generated/tilted_bottom_boundary_layer.nc
├── dimensions: zC(64), zF(65), xC(64), yF(1), xF(64), yC(1), time(0)
├── 4 outputs: (Ri, b, PV, Ro)
└── array type: Array{Float64}

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⃗|ₘₐₓ = [3.91e-03,  0.00e+00,  1.33e-03] m/s,  advective CFL = 0.011,  walltime / timestep = 0 seconds
[ Info:     ... simulation initialization complete (10.073 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (7.173 seconds).
[ Info: 007.44%,  time = 53.557 minutes,  Δt = 7.802 seconds,  |u⃗|ₘₐₓ = [6.51e-03,  8.78e-02,  3.25e-04] m/s,  advective CFL = 0.7,  walltime / timestep = 42.096 ms
[ Info: 014.54%,  time = 1.745 hours,  Δt = 7.629 seconds,  |u⃗|ₘₐₓ = [1.05e-02,  8.96e-02,  1.95e-04] m/s,  advective CFL = 0.7,  walltime / timestep = 13.128 ms
[ Info: 021.57%,  time = 2.588 hours,  Δt = 7.588 seconds,  |u⃗|ₘₐₓ = [1.07e-02,  9.03e-02,  1.54e-04] m/s,  advective CFL = 0.7,  walltime / timestep = 13.190 ms
[ Info: 028.55%,  time = 3.426 hours,  Δt = 7.560 seconds,  |u⃗|ₘₐₓ = [9.12e-03,  9.13e-02,  1.03e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 13.192 ms
[ Info: 035.49%,  time = 4.259 hours,  Δt = 7.527 seconds,  |u⃗|ₘₐₓ = [9.14e-03,  9.29e-02,  6.08e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 13.097 ms
[ Info: 042.51%,  time = 5.102 hours,  Δt = 7.628 seconds,  |u⃗|ₘₐₓ = [8.19e-03,  9.14e-02,  1.97e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 13.138 ms
[ Info: 049.44%,  time = 5.933 hours,  Δt = 7.404 seconds,  |u⃗|ₘₐₓ = [9.02e-03,  9.39e-02,  1.28e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 13.019 ms
[ Info: 056.25%,  time = 6.750 hours,  Δt = 7.363 seconds,  |u⃗|ₘₐₓ = [6.22e-03,  9.42e-02,  1.47e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 13.023 ms
[ Info: 063.03%,  time = 7.563 hours,  Δt = 7.329 seconds,  |u⃗|ₘₐₓ = [7.85e-03,  9.48e-02,  1.23e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 12.919 ms
[ Info: 069.78%,  time = 8.374 hours,  Δt = 7.284 seconds,  |u⃗|ₘₐₓ = [9.38e-03,  9.56e-02,  4.02e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 13.024 ms
[ Info: 076.50%,  time = 9.180 hours,  Δt = 7.266 seconds,  |u⃗|ₘₐₓ = [9.58e-03,  9.56e-02,  7.34e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 13.044 ms
[ Info: 083.31%,  time = 9.997 hours,  Δt = 7.493 seconds,  |u⃗|ₘₐₓ = [1.29e-02,  9.22e-02,  4.64e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 13.106 ms
[ Info: 090.16%,  time = 10.819 hours,  Δt = 7.410 seconds,  |u⃗|ₘₐₓ = [1.17e-02,  9.38e-02,  3.12e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 13.096 ms
[ Info: 096.96%,  time = 11.635 hours,  Δt = 7.314 seconds,  |u⃗|ₘₐₓ = [1.23e-02,  9.51e-02,  2.71e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 13.058 ms
[ Info: Simulation is stopping after running for 1.549 minutes.
[ 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)

RasterStack with dimensions: 
  Dim{:xC} Sampled{Float64} Float64[1.5625, 4.6875, …, 195.3125, 198.4375] ForwardOrdered Regular Points,
  Dim{:yC} Sampled{Float64} Float64[1.0] ForwardOrdered Regular Points,
  Dim{:zF} Sampled{Float64} Float64[-0.0, 0.8688167217695386, …, 93.44503754878271, 100.0] ForwardOrdered Irregular Points,
  Ti Sampled{Float64} Float64[0.0, 1200.0000000000002, …, 42000.0, 43200.0] ForwardOrdered Regular Points,
  Dim{:zC} Sampled{Float64} Float64[0.4344083608847693, 1.303282217470314, …, 90.52516773625328, 96.72251877439135] ForwardOrdered Irregular Points,
  Dim{:xF} Sampled{Float64} Float64[0.0, 3.125, …, 193.75, 196.875] ForwardOrdered Regular Points,
  Dim{:yF} Sampled{Float64} Float64[1.0] ForwardOrdered Regular Points
and 4 layers:
  :Ri Float64 dims: Dim{:xC}, Dim{:yC}, Dim{:zF}, Ti (64×1×65×37)
  :b  Float64 dims: Dim{:xC}, Dim{:yC}, Dim{:zC}, Ti (64×1×64×37)
  :PV Float64 dims: Dim{:xF}, Dim{:yF}, Dim{:zF}, Ti (64×1×65×37)
  :Ro Float64 dims: Dim{:xF}, Dim{:yF}, Dim{:zF}, Ti (64×1×65×37)

with 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.9.3\…
  "output time interval" => "Output was saved every 20 minutes."
  "date"                 => "This file was generated on 2023-10-06T16:15:17.534…
  "schedule"             => "TimeInterval"

We now use Makie to create the figure and its axes

using CairoMakie

set_theme!(Theme(fontsize = 20))
fig = Figure()

kwargs = (xlabel="x", ylabel="z", 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)

bₙ = @lift set(ds.b[Ti=$n, yC=Near(0)], :xC => X, :zC => Z)

Riₙ = @lift set(ds.Ri[Ti=$n, yC=Near(0)], :xC => X, :zF => Z)
hm1 = heatmap!(ax1, Riₙ; colormap = :coolwarm, colorrange = (-1, +1))
contour!(ax1, bₙ; levels=10, color=:white, linestyle=:dash, linewidth=0.5)
Colorbar(fig[3, 1], hm1, vertical=false, height=8, ticklabelsize=14)

Roₙ = @lift set(ds.Ro[Ti=$n, yF=Near(0)], :xF => X, :zF => Z)
hm2 = heatmap!(ax2, Roₙ; colormap = :balance, colorrange = (-10, +10))
contour!(ax2, bₙ; levels=10, color=:black, linestyle=:dash, linewidth=0.5)
Colorbar(fig[3, 2], hm2, vertical=false, height=8, ticklabelsize=14)

PVₙ = @lift set(ds.PV[Ti=$n, yF=Near(0)], :xF => X, :zF => Z)
hm3 = heatmap!(ax3, PVₙ; colormap = :coolwarm, colorrange = N²*f₀.*(-1.5, +1.5))
contour!(ax3, 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 = 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.