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

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,
                            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: 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-10/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|: 8.9e-04 m s⁻¹, wall time: 20.858 seconds
[ Info:     ... simulation initialization complete (7.723 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (14.284 seconds).
Iteration: 0200, time: 26.883 minutes, Δt: 7.728 seconds, max|w|: 3.7e-04 m s⁻¹, wall time: 42.404 seconds
Iteration: 0400, time: 52.210 minutes, Δt: 7.530 seconds, max|w|: 3.4e-04 m s⁻¹, wall time: 45.822 seconds
Iteration: 0600, time: 1.284 hours, Δt: 7.443 seconds, max|w|: 2.8e-04 m s⁻¹, wall time: 49.255 seconds
Iteration: 0800, time: 1.693 hours, Δt: 7.394 seconds, max|w|: 2.5e-04 m s⁻¹, wall time: 52.838 seconds
Iteration: 1000, time: 2.102 hours, Δt: 7.361 seconds, max|w|: 2.2e-04 m s⁻¹, wall time: 56.317 seconds
Iteration: 1200, time: 2.509 hours, Δt: 7.342 seconds, max|w|: 2.8e-04 m s⁻¹, wall time: 59.814 seconds
Iteration: 1400, time: 2.915 hours, Δt: 7.325 seconds, max|w|: 6.3e-04 m s⁻¹, wall time: 1.062 minutes
Iteration: 1600, time: 3.321 hours, Δt: 7.312 seconds, max|w|: 2.2e-03 m s⁻¹, wall time: 1.126 minutes
Iteration: 1800, time: 3.726 hours, Δt: 7.305 seconds, max|w|: 5.3e-03 m s⁻¹, wall time: 1.189 minutes
Iteration: 2000, time: 4.136 hours, Δt: 7.560 seconds, max|w|: 5.6e-03 m s⁻¹, wall time: 1.255 minutes
Iteration: 2200, time: 4.552 hours, Δt: 7.504 seconds, max|w|: 5.0e-03 m s⁻¹, wall time: 1.333 minutes
Iteration: 2400, time: 4.964 hours, Δt: 7.391 seconds, max|w|: 4.6e-03 m s⁻¹, wall time: 1.412 minutes
Iteration: 2600, time: 5.370 hours, Δt: 7.365 seconds, max|w|: 4.3e-03 m s⁻¹, wall time: 1.483 minutes
Iteration: 2800, time: 5.774 hours, Δt: 7.240 seconds, max|w|: 4.2e-03 m s⁻¹, wall time: 1.548 minutes
Iteration: 3000, time: 6.175 hours, Δt: 7.263 seconds, max|w|: 4.4e-03 m s⁻¹, wall time: 1.614 minutes
Iteration: 3200, time: 6.576 hours, Δt: 7.278 seconds, max|w|: 4.9e-03 m s⁻¹, wall time: 1.683 minutes
Iteration: 3400, time: 6.979 hours, Δt: 7.250 seconds, max|w|: 3.6e-03 m s⁻¹, wall time: 1.763 minutes
Iteration: 3600, time: 7.380 hours, Δt: 7.256 seconds, max|w|: 3.9e-03 m s⁻¹, wall time: 1.835 minutes
Iteration: 3800, time: 7.781 hours, Δt: 7.260 seconds, max|w|: 4.3e-03 m s⁻¹, wall time: 1.905 minutes
Iteration: 4000, time: 8.181 hours, Δt: 7.239 seconds, max|w|: 5.5e-03 m s⁻¹, wall time: 1.983 minutes
Iteration: 4200, time: 8.583 hours, Δt: 7.257 seconds, max|w|: 6.5e-03 m s⁻¹, wall time: 2.063 minutes
Iteration: 4400, time: 8.986 hours, Δt: 7.210 seconds, max|w|: 6.6e-03 m s⁻¹, wall time: 2.142 minutes
Iteration: 4600, time: 9.388 hours, Δt: 7.312 seconds, max|w|: 8.1e-03 m s⁻¹, wall time: 2.222 minutes
Iteration: 4800, time: 9.793 hours, Δt: 7.322 seconds, max|w|: 6.9e-03 m s⁻¹, wall time: 2.297 minutes
Iteration: 5000, time: 10.203 hours, Δt: 7.355 seconds, max|w|: 6.1e-03 m s⁻¹, wall time: 2.367 minutes
Iteration: 5200, time: 10.609 hours, Δt: 7.247 seconds, max|w|: 4.8e-03 m s⁻¹, wall time: 2.454 minutes
Iteration: 5400, time: 11.014 hours, Δt: 7.360 seconds, max|w|: 5.4e-03 m s⁻¹, wall time: 2.537 minutes
Iteration: 5600, time: 11.421 hours, Δt: 7.274 seconds, max|w|: 4.3e-03 m s⁻¹, wall time: 2.613 minutes
Iteration: 5800, time: 11.821 hours, Δt: 7.234 seconds, max|w|: 4.4e-03 m s⁻¹, wall time: 2.686 minutes
Iteration: 6000, time: 12.223 hours, Δt: 7.221 seconds, max|w|: 4.0e-03 m s⁻¹, wall time: 2.761 minutes
Iteration: 6200, time: 12.624 hours, Δt: 7.198 seconds, max|w|: 4.1e-03 m s⁻¹, wall time: 2.836 minutes
Iteration: 6400, time: 13.022 hours, Δt: 7.202 seconds, max|w|: 3.6e-03 m s⁻¹, wall time: 2.915 minutes
Iteration: 6600, time: 13.421 hours, Δt: 7.162 seconds, max|w|: 4.2e-03 m s⁻¹, wall time: 2.991 minutes
Iteration: 6800, time: 13.818 hours, Δt: 7.197 seconds, max|w|: 6.8e-03 m s⁻¹, wall time: 3.067 minutes
Iteration: 7000, time: 14.219 hours, Δt: 7.212 seconds, max|w|: 8.7e-03 m s⁻¹, wall time: 3.147 minutes
Iteration: 7200, time: 14.617 hours, Δt: 7.195 seconds, max|w|: 7.7e-03 m s⁻¹, wall time: 3.226 minutes
Iteration: 7400, time: 15.018 hours, Δt: 7.254 seconds, max|w|: 6.0e-03 m s⁻¹, wall time: 3.304 minutes
Iteration: 7600, time: 15.420 hours, Δt: 7.250 seconds, max|w|: 5.6e-03 m s⁻¹, wall time: 3.383 minutes
Iteration: 7800, time: 15.819 hours, Δt: 7.178 seconds, max|w|: 6.1e-03 m s⁻¹, wall time: 3.456 minutes
Iteration: 8000, time: 16.214 hours, Δt: 7.116 seconds, max|w|: 5.6e-03 m s⁻¹, wall time: 3.528 minutes
Iteration: 8200, time: 16.607 hours, Δt: 7.066 seconds, max|w|: 5.4e-03 m s⁻¹, wall time: 3.607 minutes
Iteration: 8400, time: 16.998 hours, Δt: 7.101 seconds, max|w|: 6.2e-03 m s⁻¹, wall time: 3.693 minutes
Iteration: 8600, time: 17.395 hours, Δt: 7.151 seconds, max|w|: 4.9e-03 m s⁻¹, wall time: 3.780 minutes
Iteration: 8800, time: 17.790 hours, Δt: 7.147 seconds, max|w|: 7.0e-03 m s⁻¹, wall time: 3.859 minutes
Iteration: 9000, time: 18.187 hours, Δt: 7.133 seconds, max|w|: 8.7e-03 m s⁻¹, wall time: 3.938 minutes
Iteration: 9200, time: 18.579 hours, Δt: 7.099 seconds, max|w|: 5.9e-03 m s⁻¹, wall time: 4.021 minutes
Iteration: 9400, time: 18.971 hours, Δt: 7.020 seconds, max|w|: 5.8e-03 m s⁻¹, wall time: 4.099 minutes
Iteration: 9600, time: 19.363 hours, Δt: 7.123 seconds, max|w|: 6.4e-03 m s⁻¹, wall time: 4.172 minutes
Iteration: 9800, time: 19.759 hours, Δt: 7.206 seconds, max|w|: 5.1e-03 m s⁻¹, wall time: 4.245 minutes
Iteration: 10000, time: 20.158 hours, Δt: 7.231 seconds, max|w|: 7.0e-03 m s⁻¹, wall time: 4.320 minutes
Iteration: 10200, time: 20.557 hours, Δt: 7.159 seconds, max|w|: 8.7e-03 m s⁻¹, wall time: 4.394 minutes
Iteration: 10400, time: 20.951 hours, Δt: 7.172 seconds, max|w|: 8.7e-03 m s⁻¹, wall time: 4.439 minutes
Iteration: 10600, time: 21.349 hours, Δt: 7.190 seconds, max|w|: 8.9e-03 m s⁻¹, wall time: 4.492 minutes
Iteration: 10800, time: 21.746 hours, Δt: 7.180 seconds, max|w|: 9.0e-03 m s⁻¹, wall time: 4.531 minutes
Iteration: 11000, time: 22.146 hours, Δt: 7.216 seconds, max|w|: 9.0e-03 m s⁻¹, wall time: 4.569 minutes
Iteration: 11200, time: 22.546 hours, Δt: 7.221 seconds, max|w|: 7.8e-03 m s⁻¹, wall time: 4.607 minutes
Iteration: 11400, time: 22.947 hours, Δt: 7.207 seconds, max|w|: 8.4e-03 m s⁻¹, wall time: 4.646 minutes
Iteration: 11600, time: 23.345 hours, Δt: 7.195 seconds, max|w|: 7.5e-03 m s⁻¹, wall time: 4.689 minutes
Iteration: 11800, time: 23.744 hours, Δt: 7.179 seconds, max|w|: 7.6e-03 m s⁻¹, wall time: 4.735 minutes
[ Info: Simulation is stopping after running for 4.479 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

This page was generated using Literate.jl.