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:

using Oceananigans
using Oceananigans.Units

Lx = 200meters
Lz = 100meters
Nx = 64
Nz = 64

grid = RectilinearGrid(topology = (Periodic, Flat, Bounded), size = (Nx, Nz),
                       x = (0, Lx), z = (0, Lz))

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]     regularly spaced with Δz=1.5625

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, 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, 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 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, 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: 7.812 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", "64×1×64 Field{Center, Center, Center} on RectilinearGrid on CPU", "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}
├── file_splitting: NoFileSplitting
└── file size: 18.2 KiB

Run the simulation and process results

To run the simulation:

run!(simulation)

[ Info: Initializing simulation...
[ Info: 000.00%,  time = 0 seconds,  Δt = 8.594 seconds,  |u⃗|ₘₐₓ = [2.37e-03,  0.00e+00,  1.46e-03] m/s,  advective CFL = 0.012,  walltime / timestep = 0 seconds
[ Info:     ... simulation initialization complete (8.208 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (5.638 seconds).
[ Info: 008.54%,  time = 1.025 hours,  Δt = 8.064 seconds,  |u⃗|ₘₐₓ = [1.06e-02,  8.36e-02,  2.79e-04] m/s,  advective CFL = 0.7,  walltime / timestep = 32.458 ms
[ Info: 015.82%,  time = 1.898 hours,  Δt = 7.763 seconds,  |u⃗|ₘₐₓ = [1.16e-02,  8.70e-02,  6.10e-04] m/s,  advective CFL = 0.7,  walltime / timestep = 5.767 ms
[ Info: 022.91%,  time = 2.750 hours,  Δt = 7.648 seconds,  |u⃗|ₘₐₓ = [1.05e-02,  8.99e-02,  2.41e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.785 ms
[ Info: 029.86%,  time = 3.583 hours,  Δt = 7.354 seconds,  |u⃗|ₘₐₓ = [8.88e-03,  9.05e-02,  7.20e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.831 ms
[ Info: 037.13%,  time = 4.456 hours,  Δt = 7.947 seconds,  |u⃗|ₘₐₓ = [1.22e-02,  8.76e-02,  2.52e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.906 ms
[ Info: 044.32%,  time = 5.318 hours,  Δt = 7.677 seconds,  |u⃗|ₘₐₓ = [9.62e-03,  9.06e-02,  2.59e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.933 ms
[ Info: 051.29%,  time = 6.155 hours,  Δt = 7.512 seconds,  |u⃗|ₘₐₓ = [8.27e-03,  9.21e-02,  3.42e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.948 ms
[ Info: 058.17%,  time = 6.981 hours,  Δt = 7.430 seconds,  |u⃗|ₘₐₓ = [1.25e-02,  9.30e-02,  3.24e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.890 ms
[ Info: 064.89%,  time = 7.787 hours,  Δt = 7.301 seconds,  |u⃗|ₘₐₓ = [1.33e-02,  9.31e-02,  7.75e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.945 ms
[ Info: 071.53%,  time = 8.584 hours,  Δt = 7.441 seconds,  |u⃗|ₘₐₓ = [1.30e-02,  9.26e-02,  6.28e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.958 ms
[ Info: 078.57%,  time = 9.428 hours,  Δt = 7.573 seconds,  |u⃗|ₘₐₓ = [1.31e-02,  9.24e-02,  5.87e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.917 ms
[ Info: 085.50%,  time = 10.260 hours,  Δt = 7.225 seconds,  |u⃗|ₘₐₓ = [1.24e-02,  9.45e-02,  3.75e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.916 ms
[ Info: 092.29%,  time = 11.075 hours,  Δt = 7.463 seconds,  |u⃗|ₘₐₓ = [1.21e-02,  9.23e-02,  4.81e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.932 ms
[ Info: 099.14%,  time = 11.897 hours,  Δt = 7.320 seconds,  |u⃗|ₘₐₓ = [1.15e-02,  9.43e-02,  7.29e-03] m/s,  advective CFL = 0.7,  walltime / timestep = 5.840 ms
[ Info: Simulation is stopping after running for 47.150 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)

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, 1.5625, …, 98.4375, 100.0] ForwardOrdered Regular Points,
  Ti Sampled{Float64} Float64[0.0, 1200.0, …, 42000.0, 43200.0] ForwardOrdered Regular Points,
  Dim{:zC} Sampled{Float64} Float64[0.78125, 2.34375, …, 97.65625, 99.21875] ForwardOrdered Regular 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.4\…
  "output time interval" => "Output was saved every 20 minutes."
  "date"                 => "This file was generated on 2024-05-23T11:26:59.885…
  "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.