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 $400 × 100$ meter $x, z$ grid with $128 × 32$ cells and finer resolution near the bottom,

using Oceananigans

Lx = 400 # m
Lz = 100 # m
Nx = 128
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,
                       halo = (3, 3))
128×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── Periodic x ∈ [0.0, 400.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(grid.Δzᵃᵃᶜ[1:Nz], grid.zᵃᵃᶜ[1:Nz],
      axis = (ylabel = "Depth (m)",
              xlabel = "Vertical spacing (m)"))

scatter!(grid.Δzᵃᵃᶜ[1:Nz], grid.zᵃᵃᶜ[1:Nz])

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

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

where we have used a constant Coriolis parameter $f = 10⁻⁴ \rm{s}⁻¹$. 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, y, z, t, p) = p.N² * (x * p.ĝ[1] + z * p.ĝ[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=(; ĝ, N² = 1e-5))
BackgroundField{typeof(Main.constant_stratification), NamedTuple{(:ĝ, :N²), Tuple{Tuple{Float64, Int64, Float64}, Float64}}}
├── func: constant_stratification (generic function with 1 method)
└── parameters: (ĝ = (0.052335956242943835, 0, 0.9986295347545738), N² = 1.0e-5)

where $N² = 10⁻⁵ \rm{s}⁻¹$ 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₁ = znodes(Center, grid)[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 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.

ν = 1e-4 # m² s⁻¹, small-ish
κ = ν
closure = ScalarDiffusivity(; ν, κ)

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))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×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 -ĝ = Tuple{Float64, Int64, Float64}
└── coriolis: ConstantCartesianCoriolis{Float64}

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 and set-up a

using Oceananigans.Units
using Oceananigans.Grids: min_Δz

simulation = Simulation(model, Δt = 0.5 * min_Δz(grid) / V∞, stop_time = 2days)
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 years
├── Stop time: 2 days
├── 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 adapt our time-step and print a progress message,

using Printf

wizard = TimeStepWizard(max_change=1.1, cfl=0.7)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(4))

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(100))
Callback of progress_message on IterationInterval(100)

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-8/clima/oceananigans/docs/build/generated/tilted_bottom_boundary_layer.nc
├── dimensions: zC(64), zF(65), xC(128), yF(1), xF(128), yC(1), time(0)
├── 5 outputs: (B, w, ωy, V, u)
└── array type: Array{Float32}

Now we just run it!

run!(simulation)
[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 4.778 seconds, max|w|: 0.0e+00 m s⁻¹, wall time: 16.341 seconds
[ Info:     ... simulation initialization complete (6.881 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (1.472 minutes).
Iteration: 0100, time: 30.941 minutes, Δt: 51.774 seconds, max|w|: 1.4e-18 m s⁻¹, wall time: 1.864 minutes
Iteration: 0200, time: 3.779 hours, Δt: 2.247 minutes, max|w|: 5.2e-16 m s⁻¹, wall time: 1.907 minutes
Iteration: 0300, time: 8.667 hours, Δt: 1.104 minutes, max|w|: 1.6e-02 m s⁻¹, wall time: 1.947 minutes
Iteration: 0400, time: 11.151 hours, Δt: 2.951 minutes, max|w|: 2.6e-03 m s⁻¹, wall time: 1.990 minutes
Iteration: 0500, time: 14.307 hours, Δt: 1.299 minutes, max|w|: 7.8e-03 m s⁻¹, wall time: 2.032 minutes
Iteration: 0600, time: 16.359 hours, Δt: 1.516 minutes, max|w|: 6.7e-03 m s⁻¹, wall time: 2.073 minutes
Iteration: 0700, time: 18.476 hours, Δt: 1.425 minutes, max|w|: 7.1e-03 m s⁻¹, wall time: 2.115 minutes
Iteration: 0800, time: 20.843 hours, Δt: 1.103 minutes, max|w|: 9.2e-03 m s⁻¹, wall time: 2.157 minutes
Iteration: 0900, time: 22.660 hours, Δt: 1.307 minutes, max|w|: 7.8e-03 m s⁻¹, wall time: 2.198 minutes
Iteration: 1000, time: 1.030 days, Δt: 1.398 minutes, max|w|: 6.7e-03 m s⁻¹, wall time: 2.237 minutes
Iteration: 1100, time: 1.124 days, Δt: 1.524 minutes, max|w|: 6.5e-03 m s⁻¹, wall time: 2.279 minutes
Iteration: 1200, time: 1.220 days, Δt: 1.310 minutes, max|w|: 7.7e-03 m s⁻¹, wall time: 2.320 minutes
Iteration: 1300, time: 1.301 days, Δt: 1.104 minutes, max|w|: 9.2e-03 m s⁻¹, wall time: 2.362 minutes
Iteration: 1400, time: 1.384 days, Δt: 1.393 minutes, max|w|: 7.3e-03 m s⁻¹, wall time: 2.403 minutes
Iteration: 1500, time: 1.458 days, Δt: 56.262 seconds, max|w|: 1.1e-02 m s⁻¹, wall time: 2.444 minutes
Iteration: 1600, time: 1.527 days, Δt: 53.041 seconds, max|w|: 1.1e-02 m s⁻¹, wall time: 2.485 minutes
Iteration: 1700, time: 1.590 days, Δt: 49.264 seconds, max|w|: 1.2e-02 m s⁻¹, wall time: 2.524 minutes
Iteration: 1800, time: 1.649 days, Δt: 53.475 seconds, max|w|: 1.1e-02 m s⁻¹, wall time: 2.565 minutes
Iteration: 1900, time: 1.712 days, Δt: 1.005 minutes, max|w|: 1.0e-02 m s⁻¹, wall time: 2.606 minutes
Iteration: 2000, time: 1.789 days, Δt: 1.374 minutes, max|w|: 7.4e-03 m s⁻¹, wall time: 2.647 minutes
Iteration: 2100, time: 1.868 days, Δt: 1.149 minutes, max|w|: 8.8e-03 m s⁻¹, wall time: 2.689 minutes
Iteration: 2200, time: 1.944 days, Δt: 1.012 minutes, max|w|: 1.0e-02 m s⁻¹, wall time: 2.732 minutes
[ Info: Simulation is stopping after running for 2.534 minutes.
[ Info: Simulation time 2 days equals or exceeds stop time 2 days.

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  …  370.3125, 373.4375, 376.5625, 379.6875, 382.8125, 385.9375, 389.0625, 392.1875, 395.3125, 398.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  …  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(resolution = (800, 600))

axis_kwargs = (xlabel = "Across-slope distance (x)",
               ylabel = "Slope-normal\ndistance (z)",
               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 = "m 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)
Label()

Finally, we record a movie.

frames = 1:length(times)

record(fig, "tilted_bottom_boundary_layer.mp4", frames, framerate=12) do i
    msg = string("Plotting frame ", i, " of ", frames[end])
    if i%5 == 0 print(msg * " \r") end
    n[] = i
end
Plotting frame 5 of 145 
Plotting frame 10 of 145 
Plotting frame 15 of 145 
Plotting frame 20 of 145 
Plotting frame 25 of 145 
Plotting frame 30 of 145 
Plotting frame 35 of 145 
Plotting frame 40 of 145 
Plotting frame 45 of 145 
Plotting frame 50 of 145 
Plotting frame 55 of 145 
Plotting frame 60 of 145 
Plotting frame 65 of 145 
Plotting frame 70 of 145 
Plotting frame 75 of 145 
Plotting frame 80 of 145 
Plotting frame 85 of 145 
Plotting frame 90 of 145 
Plotting frame 95 of 145 
Plotting frame 100 of 145 
Plotting frame 105 of 145 
Plotting frame 110 of 145 
Plotting frame 115 of 145 
Plotting frame 120 of 145 
Plotting frame 125 of 145 
Plotting frame 130 of 145 
Plotting frame 135 of 145 
Plotting frame 140 of 145 
Plotting frame 145 of 145

Don't forget to close the NetCDF file!

close(ds)
closed NetCDF NCDataset

This page was generated using Literate.jl.