# Tilted bottom boundary layer example

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

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

# Creates a grid with near-constant spacing refinement * Lz / Nz
# near the bottom:
refinement = 1.8 # controls spacing near surface (higher means finer spaced)
stretching = 10  # controls rate of stretching at bottom

# "Warped" height coordinate
h(k) = (Nz + 1 - k) / Nz

# Linear near-surface generator
ζ(k) = 1 + (h(k) - 1) / refinement

# Bottom-intensified stretching function
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))

# Generating function
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

Let's make sure the grid spacing is both finer and near-uniform at the bottom,

using CairoMakie

lines(zspacings(grid, Center()), znodes(grid, Center()),
axis = (ylabel = "Depth (m)",
xlabel = "Vertical spacing (m)"))

scatter!(zspacings(grid, Center()), znodes(grid, Center()))


## Tilting the domain

We use a domain that's tilted with respect to gravity by

θ = 3 # degrees
3

so that $x$ is the along-slope direction, $z$ is the across-slope 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.052335956242943835
0.0
0.9986295347545738

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̂)
coriolis = ConstantCartesianCoriolis(f = 1e-4, rotation_axis = ĝ)
ConstantCartesianCoriolis{Float64}: fx = 5.23e-06, fy = 0.00e+00, fz = 9.99e-05

where above we used a constant Coriolis parameter $f = 10^{-4} \, \rm{s}^{-1}$. The tilting also affects the kind of density stratified flows we can model. In particular, a constant density stratification in the tilted coordinate system

@inline constant_stratification(x, z, t, p) = p.N² * (x * p.ĝ[1] + z * p.ĝ[3])
constant_stratification (generic function with 1 method)

is not periodic in $x$. Thus we cannot explicitly model a constant stratification 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² = 1e-5))
BackgroundField{typeof(Main.var"##269".constant_stratification), @NamedTuple{ĝ::Vector{Float64}, N²::Float64}}
├── func: constant_stratification (generic function with 1 method)
└── parameters: (ĝ = [0.052335956242943835, 0.0, 0.9986295347545738], N² = 1.0e-5)

where $N^2 = 10^{-5} \rm{s}^{-2}$ is the background buoyancy gradient.

## Bottom drag

We impose bottom drag that follows Monin–Obukhov theory. We include the background flow in the drag calculation, which is the only effect the background flow enters the problem,

V∞ = 0.1 # m s⁻¹
z₀ = 0.1 # m (roughness length)
κ = 0.4  # von Karman constant

z₁ = first(znodes(grid, Center())) # 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 the NonhydrostaticModel

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. Here we use a smallish value of $10^{-4} \, \rm{m}^2\, \rm{s}^{-1}$.

closure = ScalarDiffusivity(ν=1e-4, κ=1e-4)

model = NonhydrostaticModel(; grid, buoyancy, coriolis, closure,
timestepper = :RungeKutta3,
tracers = :b,
boundary_conditions = (u=u_bcs, v=v_bcs),
background_fields = (; b=B_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: Upwind Biased reconstruction order 5
├── tracers: b
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0001, κ=(b=0.0001,))
├── buoyancy: BuoyancyTracer with ĝ = Tuple{Float64, Float64, Float64}
└── 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 grid size of our domain.

simulation = Simulation(model, Δt = 0.5 * minimum_zspacing(grid) / V∞, stop_time = 1day)
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: 1 day
├── 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 a TimeStepWizard to adapt our time-step,

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)

and also we add another callback to print a progress message,

using Printf

start_time = time_ns() # so we can print the total elapsed wall time

progress_message(sim) =
@printf("Iteration: %04d, time: %s, Δt: %s, max|w|: %.1e m s⁻¹, wall time: %s\n",
iteration(sim), prettytime(time(sim)),
prettytime(sim.Δt), maximum(abs, sim.model.velocities.w),
prettytime((time_ns() - start_time) * 1e-9))

simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(200))
Callback of progress_message on IterationInterval(200)

## Add outputs to the simulation

We add outputs to our model using the NetCDFOutputWriter,

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)

simulation.output_writers[:fields] = NetCDFOutputWriter(model, outputs;
filename = joinpath(@__DIR__, "tilted_bottom_boundary_layer.nc"),
schedule = TimeInterval(20minutes),
overwrite_existing = true)
NetCDFOutputWriter scheduled on TimeInterval(20 minutes):
├── filepath: /var/lib/buildkite-agent/builds/tartarus-2/clima/oceananigans/docs/src/literated/tilted_bottom_boundary_layer.nc
├── dimensions: zC(64), zF(65), xC(64), yF(1), xF(64), yC(1), time(0)
├── 5 outputs: (B, w, ωy, V, u)
└── array type: Array{Float64}
├── file_splitting: NoFileSplitting
└── file size: 19.2 KiB

Now we just run it!

run!(simulation)
[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 4.778 seconds, max|w|: 1.3e-03 m s⁻¹, wall time: 9.887 seconds
[ Info:     ... simulation initialization complete (3.679 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (6.269 seconds).
Iteration: 0200, time: 26.888 minutes, Δt: 7.720 seconds, max|w|: 4.2e-04 m s⁻¹, wall time: 19.585 seconds
Iteration: 0400, time: 52.234 minutes, Δt: 7.541 seconds, max|w|: 3.6e-04 m s⁻¹, wall time: 21.155 seconds
Iteration: 0600, time: 1.287 hours, Δt: 7.446 seconds, max|w|: 2.6e-04 m s⁻¹, wall time: 22.727 seconds
Iteration: 0800, time: 1.697 hours, Δt: 7.391 seconds, max|w|: 2.7e-04 m s⁻¹, wall time: 24.300 seconds
Iteration: 1000, time: 2.106 hours, Δt: 7.359 seconds, max|w|: 2.1e-04 m s⁻¹, wall time: 25.872 seconds
Iteration: 1200, time: 2.513 hours, Δt: 7.337 seconds, max|w|: 3.2e-04 m s⁻¹, wall time: 27.440 seconds
Iteration: 1400, time: 2.919 hours, Δt: 7.325 seconds, max|w|: 6.7e-04 m s⁻¹, wall time: 29.008 seconds
Iteration: 1600, time: 3.325 hours, Δt: 7.314 seconds, max|w|: 2.1e-03 m s⁻¹, wall time: 30.573 seconds
Iteration: 1800, time: 3.730 hours, Δt: 7.337 seconds, max|w|: 5.2e-03 m s⁻¹, wall time: 32.166 seconds
Iteration: 2000, time: 4.142 hours, Δt: 7.524 seconds, max|w|: 5.8e-03 m s⁻¹, wall time: 33.741 seconds
Iteration: 2200, time: 4.559 hours, Δt: 7.518 seconds, max|w|: 5.1e-03 m s⁻¹, wall time: 35.327 seconds
Iteration: 2400, time: 4.973 hours, Δt: 7.402 seconds, max|w|: 4.1e-03 m s⁻¹, wall time: 36.905 seconds
Iteration: 2600, time: 5.376 hours, Δt: 7.285 seconds, max|w|: 3.8e-03 m s⁻¹, wall time: 38.483 seconds
Iteration: 2800, time: 5.780 hours, Δt: 7.301 seconds, max|w|: 3.8e-03 m s⁻¹, wall time: 40.050 seconds
Iteration: 3000, time: 6.183 hours, Δt: 7.223 seconds, max|w|: 4.1e-03 m s⁻¹, wall time: 41.617 seconds
Iteration: 3200, time: 6.585 hours, Δt: 7.252 seconds, max|w|: 3.9e-03 m s⁻¹, wall time: 43.193 seconds
Iteration: 3400, time: 6.987 hours, Δt: 7.224 seconds, max|w|: 3.5e-03 m s⁻¹, wall time: 44.766 seconds
Iteration: 3600, time: 7.386 hours, Δt: 7.231 seconds, max|w|: 3.2e-03 m s⁻¹, wall time: 46.361 seconds
Iteration: 3800, time: 7.787 hours, Δt: 7.228 seconds, max|w|: 4.0e-03 m s⁻¹, wall time: 47.943 seconds
Iteration: 4000, time: 8.190 hours, Δt: 7.301 seconds, max|w|: 5.4e-03 m s⁻¹, wall time: 49.553 seconds
Iteration: 4200, time: 8.594 hours, Δt: 7.274 seconds, max|w|: 6.5e-03 m s⁻¹, wall time: 51.195 seconds
Iteration: 4400, time: 9 hours, Δt: 7.294 seconds, max|w|: 4.9e-03 m s⁻¹, wall time: 52.798 seconds
Iteration: 4600, time: 9.404 hours, Δt: 7.311 seconds, max|w|: 6.6e-03 m s⁻¹, wall time: 54.381 seconds
Iteration: 4800, time: 9.811 hours, Δt: 7.362 seconds, max|w|: 9.3e-03 m s⁻¹, wall time: 55.960 seconds
Iteration: 5000, time: 10.216 hours, Δt: 7.333 seconds, max|w|: 6.3e-03 m s⁻¹, wall time: 57.545 seconds
Iteration: 5200, time: 10.623 hours, Δt: 7.380 seconds, max|w|: 5.8e-03 m s⁻¹, wall time: 59.130 seconds
Iteration: 5400, time: 11.033 hours, Δt: 7.390 seconds, max|w|: 6.4e-03 m s⁻¹, wall time: 1.012 minutes
Iteration: 5600, time: 11.441 hours, Δt: 7.284 seconds, max|w|: 5.4e-03 m s⁻¹, wall time: 1.039 minutes
Iteration: 5800, time: 11.846 hours, Δt: 7.298 seconds, max|w|: 5.0e-03 m s⁻¹, wall time: 1.066 minutes
Iteration: 6000, time: 12.249 hours, Δt: 7.270 seconds, max|w|: 6.3e-03 m s⁻¹, wall time: 1.103 minutes
Iteration: 6200, time: 12.653 hours, Δt: 7.264 seconds, max|w|: 4.1e-03 m s⁻¹, wall time: 1.133 minutes
Iteration: 6400, time: 13.052 hours, Δt: 7.169 seconds, max|w|: 5.0e-03 m s⁻¹, wall time: 1.160 minutes
Iteration: 6600, time: 13.451 hours, Δt: 7.183 seconds, max|w|: 3.7e-03 m s⁻¹, wall time: 1.187 minutes
Iteration: 6800, time: 13.850 hours, Δt: 7.211 seconds, max|w|: 4.3e-03 m s⁻¹, wall time: 1.213 minutes
Iteration: 7000, time: 14.250 hours, Δt: 7.174 seconds, max|w|: 7.4e-03 m s⁻¹, wall time: 1.239 minutes
Iteration: 7200, time: 14.647 hours, Δt: 7.150 seconds, max|w|: 7.1e-03 m s⁻¹, wall time: 1.267 minutes
Iteration: 7400, time: 15.048 hours, Δt: 7.197 seconds, max|w|: 5.5e-03 m s⁻¹, wall time: 1.294 minutes
Iteration: 7600, time: 15.448 hours, Δt: 7.251 seconds, max|w|: 6.7e-03 m s⁻¹, wall time: 1.321 minutes
Iteration: 7800, time: 15.851 hours, Δt: 7.265 seconds, max|w|: 7.1e-03 m s⁻¹, wall time: 1.348 minutes
Iteration: 8000, time: 16.255 hours, Δt: 7.263 seconds, max|w|: 6.5e-03 m s⁻¹, wall time: 1.375 minutes
Iteration: 8200, time: 16.656 hours, Δt: 7.218 seconds, max|w|: 5.5e-03 m s⁻¹, wall time: 1.402 minutes
Iteration: 8400, time: 17.054 hours, Δt: 7.164 seconds, max|w|: 5.6e-03 m s⁻¹, wall time: 1.429 minutes
Iteration: 8600, time: 17.450 hours, Δt: 7.111 seconds, max|w|: 6.3e-03 m s⁻¹, wall time: 1.456 minutes
Iteration: 8800, time: 17.847 hours, Δt: 7.213 seconds, max|w|: 4.9e-03 m s⁻¹, wall time: 1.483 minutes
Iteration: 9000, time: 18.250 hours, Δt: 7.233 seconds, max|w|: 6.6e-03 m s⁻¹, wall time: 1.510 minutes
Iteration: 9200, time: 18.651 hours, Δt: 7.258 seconds, max|w|: 5.3e-03 m s⁻¹, wall time: 1.537 minutes
Iteration: 9400, time: 19.050 hours, Δt: 7.237 seconds, max|w|: 6.3e-03 m s⁻¹, wall time: 1.564 minutes
Iteration: 9600, time: 19.452 hours, Δt: 7.263 seconds, max|w|: 5.4e-03 m s⁻¹, wall time: 1.597 minutes
Iteration: 9800, time: 19.855 hours, Δt: 7.170 seconds, max|w|: 6.1e-03 m s⁻¹, wall time: 1.636 minutes
Iteration: 10000, time: 20.254 hours, Δt: 7.222 seconds, max|w|: 5.8e-03 m s⁻¹, wall time: 1.672 minutes
Iteration: 10200, time: 20.655 hours, Δt: 7.236 seconds, max|w|: 5.7e-03 m s⁻¹, wall time: 1.713 minutes
Iteration: 10400, time: 21.054 hours, Δt: 7.168 seconds, max|w|: 7.6e-03 m s⁻¹, wall time: 1.754 minutes
Iteration: 10600, time: 21.451 hours, Δt: 7.150 seconds, max|w|: 5.3e-03 m s⁻¹, wall time: 1.795 minutes
Iteration: 10800, time: 21.848 hours, Δt: 7.176 seconds, max|w|: 5.6e-03 m s⁻¹, wall time: 1.836 minutes
Iteration: 11000, time: 22.248 hours, Δt: 7.177 seconds, max|w|: 6.3e-03 m s⁻¹, wall time: 1.879 minutes
Iteration: 11200, time: 22.645 hours, Δt: 7.146 seconds, max|w|: 6.5e-03 m s⁻¹, wall time: 1.919 minutes
Iteration: 11400, time: 23.040 hours, Δt: 7.220 seconds, max|w|: 6.6e-03 m s⁻¹, wall time: 1.961 minutes
Iteration: 11600, time: 23.439 hours, Δt: 7.174 seconds, max|w|: 5.5e-03 m s⁻¹, wall time: 2.002 minutes
Iteration: 11800, time: 23.839 hours, Δt: 7.200 seconds, max|w|: 6.7e-03 m s⁻¹, wall time: 2.043 minutes
[ Info: Simulation is stopping after running for 1.926 minutes.
[ Info: Simulation time 1 day equals or exceeds stop time 1 day.


## Visualize the results

First we load the required package to load NetCDF output files and define the coordinates for plotting using existing objects:

using NCDatasets, CairoMakie

xω, yω, zω = nodes(ωy)
xv, yv, zv = nodes(V)
([1.5625, 4.6875, 7.8125, 10.9375, 14.0625, 17.1875, 20.3125, 23.4375, 26.5625, 29.6875, 32.8125, 35.9375, 39.0625, 42.1875, 45.3125, 48.4375, 51.5625, 54.6875, 57.8125, 60.9375, 64.0625, 67.1875, 70.3125, 73.4375, 76.5625, 79.6875, 82.8125, 85.9375, 89.0625, 92.1875, 95.3125, 98.4375, 101.5625, 104.6875, 107.8125, 110.9375, 114.0625, 117.1875, 120.3125, 123.4375, 126.5625, 129.6875, 132.8125, 135.9375, 139.0625, 142.1875, 145.3125, 148.4375, 151.5625, 154.6875, 157.8125, 160.9375, 164.0625, 167.1875, 170.3125, 173.4375, 176.5625, 179.6875, 182.8125, 185.9375, 189.0625, 192.1875, 195.3125, 198.4375], StepRangeLen(1.0, 0.0, 1), [0.4344083608847693, 1.303282217470314, 2.1722793473188617, 3.041419168936904, 3.9107241431416084, 4.780220246692868, 5.649937519111903, 6.519910693894487, 7.390179927024476, 8.26079163764491, 9.131799477986252, 10.00326545222724, 10.875261206921866, 11.747869519021842, 12.621186011421164, 13.495321130420056, 14.370402424632783, 15.24657717074116, 16.124015398230675, 17.00291337296054, 17.88349760825156, 18.766029482284154, 19.650810552161538, 20.538188668213625, 21.428565007224726, 22.322402160523495, 23.22023343257072, 24.122673528152244, 25.030430831905008, 25.944321513097137, 26.865285721823064, 27.794406180595843, 28.732929518320184, 29.682290742481506, 30.64414130083751, 31.62038124678409, 32.61319609381579, 33.62509902514256, 34.65897921569778, 35.71815712673464, 36.80644774933731, 37.92823290397964, 39.08854385039015, 40.29315562720595, 41.548694726134066, 42.86276191262984, 44.24407223660885, 45.702614534695954, 47.24983301231516, 48.89883381190526, 50.66461982505285, 52.56435739360966, 54.61767896995204, 56.84702627115933, 59.27803896724765, 61.93999449018695, 64.86630513774502, 68.09507927238374, 71.66975407682737, 75.63980801853194, 80.06156188607441, 84.99907797852119, 90.52516773625328, 96.72251877439135])

Read in the simulation's output_writer for the two-dimensional fields and then create an animation showing the $y$-component of vorticity.

ds = NCDataset(simulation.output_writers[:fields].filepath, "r")

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 = @lift ds["ωy"][:, 1, :, $n] hm_ω = heatmap!(ax_ω, xω, zω, ωy, colorrange = (-0.015, +0.015), colormap = :balance) Colorbar(fig[2, 2], hm_ω; label = "s⁻¹") V = @lift ds["V"][:, 1, :,$n]
V_max = @lift maximum(abs, ds["V"][:, 1, :, $n]) hm_v = heatmap!(ax_v, xv, zv, V, colorrange = (-V∞, +V∞), colormap = :balance) Colorbar(fig[3, 2], hm_v; label = "m s⁻¹") times = collect(ds["time"]) title = @lift "t = " * string(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
end

Don't forget to close the NetCDF file!

close(ds)
closed Dataset