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-1/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: 51.550 seconds
[ Info:     ... simulation initialization complete (9.648 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (1.188 minutes).
Iteration: 0100, time: 30.941 minutes, Δt: 51.774 seconds, max|w|: 2.1e-18 m s⁻¹, wall time: 2.242 minutes
Iteration: 0200, time: 3.779 hours, Δt: 2.247 minutes, max|w|: 3.5e-16 m s⁻¹, wall time: 2.324 minutes
Iteration: 0300, time: 8.667 hours, Δt: 1.325 minutes, max|w|: 1.5e-02 m s⁻¹, wall time: 2.409 minutes
Iteration: 0400, time: 10.600 hours, Δt: 1.286 minutes, max|w|: 7.9e-03 m s⁻¹, wall time: 2.490 minutes
Iteration: 0500, time: 12.378 hours, Δt: 1.337 minutes, max|w|: 7.6e-03 m s⁻¹, wall time: 2.571 minutes
Iteration: 0600, time: 14.503 hours, Δt: 1.159 minutes, max|w|: 8.7e-03 m s⁻¹, wall time: 2.654 minutes
Iteration: 0700, time: 16.435 hours, Δt: 1.223 minutes, max|w|: 8.3e-03 m s⁻¹, wall time: 2.734 minutes
Iteration: 0800, time: 18.173 hours, Δt: 1.106 minutes, max|w|: 9.2e-03 m s⁻¹, wall time: 2.813 minutes
Iteration: 0900, time: 19.876 hours, Δt: 1.128 minutes, max|w|: 9.0e-03 m s⁻¹, wall time: 2.894 minutes
Iteration: 1000, time: 21.439 hours, Δt: 1.112 minutes, max|w|: 9.1e-03 m s⁻¹, wall time: 2.973 minutes
Iteration: 1100, time: 23.044 hours, Δt: 53.628 seconds, max|w|: 1.1e-02 m s⁻¹, wall time: 3.053 minutes
Iteration: 1200, time: 1.028 days, Δt: 53.548 seconds, max|w|: 1.1e-02 m s⁻¹, wall time: 3.135 minutes
Iteration: 1300, time: 1.090 days, Δt: 58.765 seconds, max|w|: 1.0e-02 m s⁻¹, wall time: 3.215 minutes
Iteration: 1400, time: 1.152 days, Δt: 46.371 seconds, max|w|: 1.3e-02 m s⁻¹, wall time: 3.294 minutes
Iteration: 1500, time: 1.217 days, Δt: 48.449 seconds, max|w|: 1.3e-02 m s⁻¹, wall time: 3.376 minutes
Iteration: 1600, time: 1.284 days, Δt: 57.361 seconds, max|w|: 1.1e-02 m s⁻¹, wall time: 3.456 minutes
Iteration: 1700, time: 1.346 days, Δt: 55.337 seconds, max|w|: 1.1e-02 m s⁻¹, wall time: 3.535 minutes
Iteration: 1800, time: 1.404 days, Δt: 55.659 seconds, max|w|: 1.1e-02 m s⁻¹, wall time: 3.617 minutes
Iteration: 1900, time: 1.466 days, Δt: 50.889 seconds, max|w|: 1.2e-02 m s⁻¹, wall time: 3.697 minutes
Iteration: 2000, time: 1.527 days, Δt: 53.665 seconds, max|w|: 1.1e-02 m s⁻¹, wall time: 3.777 minutes
Iteration: 2100, time: 1.587 days, Δt: 59.742 seconds, max|w|: 1.0e-02 m s⁻¹, wall time: 3.858 minutes
Iteration: 2200, time: 1.653 days, Δt: 58.984 seconds, max|w|: 1.0e-02 m s⁻¹, wall time: 3.939 minutes
Iteration: 2300, time: 1.719 days, Δt: 56.479 seconds, max|w|: 1.1e-02 m s⁻¹, wall time: 4.019 minutes
Iteration: 2400, time: 1.785 days, Δt: 1.021 minutes, max|w|: 9.9e-03 m s⁻¹, wall time: 4.101 minutes
Iteration: 2500, time: 1.844 days, Δt: 47.917 seconds, max|w|: 1.3e-02 m s⁻¹, wall time: 4.179 minutes
Iteration: 2600, time: 1.905 days, Δt: 58.335 seconds, max|w|: 1.0e-02 m s⁻¹, wall time: 4.258 minutes
Iteration: 2700, time: 1.969 days, Δt: 57.318 seconds, max|w|: 1.1e-02 m s⁻¹, wall time: 4.340 minutes
[ Info: Simulation is stopping. Model time 2 days has hit or exceeded simulation 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], 1.0:0.0:1.0, [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, textsize=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.