# 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-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.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"##319".constant_stratification), NamedTuple{(:ĝ, :N²), Tuple{Vector{Float64}, 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
├── 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-4/clima/oceananigans/docs/src/generated/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}

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: 11.929 seconds
[ Info:     ... simulation initialization complete (6.454 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (6.813 seconds).
Iteration: 0200, time: 26.871 minutes, Δt: 7.703 seconds, max|w|: 5.4e-04 m s⁻¹, wall time: 25.344 seconds
Iteration: 0400, time: 52.004 minutes, Δt: 7.478 seconds, max|w|: 3.5e-04 m s⁻¹, wall time: 28.273 seconds
Iteration: 0600, time: 1.280 hours, Δt: 7.436 seconds, max|w|: 2.9e-04 m s⁻¹, wall time: 31.172 seconds
Iteration: 0800, time: 1.689 hours, Δt: 7.388 seconds, max|w|: 2.5e-04 m s⁻¹, wall time: 34.137 seconds
Iteration: 1000, time: 2.098 hours, Δt: 7.356 seconds, max|w|: 2.2e-04 m s⁻¹, wall time: 37.068 seconds
Iteration: 1200, time: 2.505 hours, Δt: 7.342 seconds, max|w|: 2.6e-04 m s⁻¹, wall time: 40.024 seconds
Iteration: 1400, time: 2.911 hours, Δt: 7.325 seconds, max|w|: 6.4e-04 m s⁻¹, wall time: 43.012 seconds
Iteration: 1600, time: 3.317 hours, Δt: 7.310 seconds, max|w|: 2.1e-03 m s⁻¹, wall time: 46.351 seconds
Iteration: 1800, time: 3.720 hours, Δt: 7.357 seconds, max|w|: 5.6e-03 m s⁻¹, wall time: 49.325 seconds
Iteration: 2000, time: 4.127 hours, Δt: 7.394 seconds, max|w|: 5.7e-03 m s⁻¹, wall time: 52.255 seconds
Iteration: 2200, time: 4.542 hours, Δt: 7.504 seconds, max|w|: 5.3e-03 m s⁻¹, wall time: 55.222 seconds
Iteration: 2400, time: 4.957 hours, Δt: 7.422 seconds, max|w|: 4.8e-03 m s⁻¹, wall time: 58.153 seconds
Iteration: 2600, time: 5.366 hours, Δt: 7.367 seconds, max|w|: 4.5e-03 m s⁻¹, wall time: 1.019 minutes
Iteration: 2800, time: 5.772 hours, Δt: 7.266 seconds, max|w|: 4.2e-03 m s⁻¹, wall time: 1.068 minutes
Iteration: 3000, time: 6.173 hours, Δt: 7.230 seconds, max|w|: 4.8e-03 m s⁻¹, wall time: 1.118 minutes
Iteration: 3200, time: 6.575 hours, Δt: 7.251 seconds, max|w|: 4.4e-03 m s⁻¹, wall time: 1.167 minutes
Iteration: 3400, time: 6.977 hours, Δt: 7.244 seconds, max|w|: 4.1e-03 m s⁻¹, wall time: 1.217 minutes
Iteration: 3600, time: 7.378 hours, Δt: 7.284 seconds, max|w|: 3.3e-03 m s⁻¹, wall time: 1.267 minutes
Iteration: 3800, time: 7.780 hours, Δt: 7.277 seconds, max|w|: 5.8e-03 m s⁻¹, wall time: 1.317 minutes
Iteration: 4000, time: 8.181 hours, Δt: 7.236 seconds, max|w|: 6.0e-03 m s⁻¹, wall time: 1.366 minutes
Iteration: 4200, time: 8.585 hours, Δt: 7.280 seconds, max|w|: 6.7e-03 m s⁻¹, wall time: 1.416 minutes
Iteration: 4400, time: 8.986 hours, Δt: 7.216 seconds, max|w|: 5.4e-03 m s⁻¹, wall time: 1.466 minutes
Iteration: 4600, time: 9.391 hours, Δt: 7.358 seconds, max|w|: 5.7e-03 m s⁻¹, wall time: 1.516 minutes
Iteration: 4800, time: 9.797 hours, Δt: 7.301 seconds, max|w|: 5.2e-03 m s⁻¹, wall time: 1.566 minutes
Iteration: 5000, time: 10.204 hours, Δt: 7.375 seconds, max|w|: 5.6e-03 m s⁻¹, wall time: 1.615 minutes
Iteration: 5200, time: 10.611 hours, Δt: 7.302 seconds, max|w|: 4.6e-03 m s⁻¹, wall time: 1.665 minutes
Iteration: 5400, time: 11.012 hours, Δt: 7.285 seconds, max|w|: 5.0e-03 m s⁻¹, wall time: 1.714 minutes
Iteration: 5600, time: 11.414 hours, Δt: 7.229 seconds, max|w|: 5.2e-03 m s⁻¹, wall time: 1.763 minutes
Iteration: 5800, time: 11.818 hours, Δt: 7.254 seconds, max|w|: 4.3e-03 m s⁻¹, wall time: 1.812 minutes
Iteration: 6000, time: 12.218 hours, Δt: 7.146 seconds, max|w|: 3.8e-03 m s⁻¹, wall time: 1.861 minutes
Iteration: 6200, time: 12.617 hours, Δt: 7.216 seconds, max|w|: 4.4e-03 m s⁻¹, wall time: 1.910 minutes
Iteration: 6400, time: 13.016 hours, Δt: 7.181 seconds, max|w|: 4.5e-03 m s⁻¹, wall time: 1.959 minutes
Iteration: 6600, time: 13.416 hours, Δt: 7.219 seconds, max|w|: 6.5e-03 m s⁻¹, wall time: 2.007 minutes
Iteration: 6800, time: 13.815 hours, Δt: 7.188 seconds, max|w|: 6.6e-03 m s⁻¹, wall time: 2.056 minutes
Iteration: 7000, time: 14.216 hours, Δt: 7.245 seconds, max|w|: 6.6e-03 m s⁻¹, wall time: 2.105 minutes
Iteration: 7200, time: 14.616 hours, Δt: 7.272 seconds, max|w|: 4.5e-03 m s⁻¹, wall time: 2.154 minutes
Iteration: 7400, time: 15.014 hours, Δt: 7.092 seconds, max|w|: 5.8e-03 m s⁻¹, wall time: 2.203 minutes
Iteration: 7600, time: 15.409 hours, Δt: 7.205 seconds, max|w|: 5.3e-03 m s⁻¹, wall time: 2.252 minutes
Iteration: 7800, time: 15.808 hours, Δt: 7.191 seconds, max|w|: 4.9e-03 m s⁻¹, wall time: 2.301 minutes
Iteration: 8000, time: 16.209 hours, Δt: 7.228 seconds, max|w|: 5.0e-03 m s⁻¹, wall time: 2.351 minutes
Iteration: 8200, time: 16.612 hours, Δt: 7.247 seconds, max|w|: 4.6e-03 m s⁻¹, wall time: 2.401 minutes
Iteration: 8400, time: 17.008 hours, Δt: 7.299 seconds, max|w|: 6.2e-03 m s⁻¹, wall time: 2.450 minutes
Iteration: 8600, time: 17.410 hours, Δt: 7.286 seconds, max|w|: 6.8e-03 m s⁻¹, wall time: 2.499 minutes
Iteration: 8800, time: 17.811 hours, Δt: 7.228 seconds, max|w|: 6.9e-03 m s⁻¹, wall time: 2.548 minutes
Iteration: 9000, time: 18.215 hours, Δt: 7.273 seconds, max|w|: 7.7e-03 m s⁻¹, wall time: 2.598 minutes
Iteration: 9200, time: 18.619 hours, Δt: 7.253 seconds, max|w|: 4.4e-03 m s⁻¹, wall time: 2.647 minutes
Iteration: 9400, time: 19.022 hours, Δt: 7.254 seconds, max|w|: 4.4e-03 m s⁻¹, wall time: 2.696 minutes
Iteration: 9600, time: 19.423 hours, Δt: 7.183 seconds, max|w|: 4.0e-03 m s⁻¹, wall time: 2.746 minutes
Iteration: 9800, time: 19.822 hours, Δt: 7.259 seconds, max|w|: 4.4e-03 m s⁻¹, wall time: 2.795 minutes
Iteration: 10000, time: 20.223 hours, Δt: 7.211 seconds, max|w|: 6.0e-03 m s⁻¹, wall time: 2.844 minutes
Iteration: 10200, time: 20.622 hours, Δt: 7.185 seconds, max|w|: 6.5e-03 m s⁻¹, wall time: 2.893 minutes
Iteration: 10400, time: 21.020 hours, Δt: 7.186 seconds, max|w|: 5.8e-03 m s⁻¹, wall time: 2.942 minutes
Iteration: 10600, time: 21.417 hours, Δt: 7.155 seconds, max|w|: 6.3e-03 m s⁻¹, wall time: 2.991 minutes
Iteration: 10800, time: 21.814 hours, Δt: 7.153 seconds, max|w|: 6.2e-03 m s⁻¹, wall time: 3.040 minutes
Iteration: 11000, time: 22.211 hours, Δt: 7.185 seconds, max|w|: 6.2e-03 m s⁻¹, wall time: 3.088 minutes
Iteration: 11200, time: 22.604 hours, Δt: 7.114 seconds, max|w|: 5.6e-03 m s⁻¹, wall time: 3.138 minutes
Iteration: 11400, time: 23 hours, Δt: 7.153 seconds, max|w|: 5.8e-03 m s⁻¹, wall time: 3.187 minutes
Iteration: 11600, time: 23.399 hours, Δt: 7.209 seconds, max|w|: 5.3e-03 m s⁻¹, wall time: 3.236 minutes
Iteration: 11800, time: 23.795 hours, Δt: 7.103 seconds, max|w|: 6.3e-03 m s⁻¹, wall time: 3.285 minutes
[ Info: Simulation is stopping after running for 3.159 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