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-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: 17.340 seconds
[ Info:     ... simulation initialization complete (5.077 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (12.998 seconds).
Iteration: 0200, time: 26.730 minutes, Δt: 7.679 seconds, max|w|: 4.9e-04 m s⁻¹, wall time: 37.205 seconds
Iteration: 0400, time: 51.951 minutes, Δt: 7.524 seconds, max|w|: 3.7e-04 m s⁻¹, wall time: 41.143 seconds
Iteration: 0600, time: 1.280 hours, Δt: 7.435 seconds, max|w|: 2.6e-04 m s⁻¹, wall time: 45.393 seconds
Iteration: 0800, time: 1.691 hours, Δt: 7.389 seconds, max|w|: 3.0e-04 m s⁻¹, wall time: 49.263 seconds
Iteration: 1000, time: 2.100 hours, Δt: 7.354 seconds, max|w|: 2.6e-04 m s⁻¹, wall time: 52.316 seconds
Iteration: 1200, time: 2.507 hours, Δt: 7.340 seconds, max|w|: 2.3e-04 m s⁻¹, wall time: 55.558 seconds
Iteration: 1400, time: 2.913 hours, Δt: 7.329 seconds, max|w|: 4.3e-04 m s⁻¹, wall time: 1.000 minutes
Iteration: 1600, time: 3.319 hours, Δt: 7.316 seconds, max|w|: 1.4e-03 m s⁻¹, wall time: 1.059 minutes
Iteration: 1800, time: 3.724 hours, Δt: 7.343 seconds, max|w|: 4.0e-03 m s⁻¹, wall time: 1.121 minutes
Iteration: 2000, time: 4.131 hours, Δt: 7.372 seconds, max|w|: 5.9e-03 m s⁻¹, wall time: 1.190 minutes
Iteration: 2200, time: 4.545 hours, Δt: 7.523 seconds, max|w|: 4.9e-03 m s⁻¹, wall time: 1.268 minutes
Iteration: 2400, time: 4.961 hours, Δt: 7.501 seconds, max|w|: 5.1e-03 m s⁻¹, wall time: 1.342 minutes
Iteration: 2600, time: 5.372 hours, Δt: 7.353 seconds, max|w|: 5.2e-03 m s⁻¹, wall time: 1.410 minutes
Iteration: 2800, time: 5.779 hours, Δt: 7.334 seconds, max|w|: 3.8e-03 m s⁻¹, wall time: 1.459 minutes
Iteration: 3000, time: 6.182 hours, Δt: 7.249 seconds, max|w|: 4.1e-03 m s⁻¹, wall time: 1.510 minutes
Iteration: 3200, time: 6.583 hours, Δt: 7.256 seconds, max|w|: 4.1e-03 m s⁻¹, wall time: 1.556 minutes
Iteration: 3400, time: 6.986 hours, Δt: 7.258 seconds, max|w|: 4.0e-03 m s⁻¹, wall time: 1.603 minutes
Iteration: 3600, time: 7.386 hours, Δt: 7.240 seconds, max|w|: 3.8e-03 m s⁻¹, wall time: 1.653 minutes
Iteration: 3800, time: 7.787 hours, Δt: 7.231 seconds, max|w|: 3.6e-03 m s⁻¹, wall time: 1.712 minutes
Iteration: 4000, time: 8.189 hours, Δt: 7.231 seconds, max|w|: 4.4e-03 m s⁻¹, wall time: 1.791 minutes
Iteration: 4200, time: 8.593 hours, Δt: 7.287 seconds, max|w|: 5.6e-03 m s⁻¹, wall time: 1.827 minutes
Iteration: 4400, time: 8.997 hours, Δt: 7.282 seconds, max|w|: 6.5e-03 m s⁻¹, wall time: 1.866 minutes
Iteration: 4600, time: 9.400 hours, Δt: 7.303 seconds, max|w|: 8.4e-03 m s⁻¹, wall time: 1.899 minutes
Iteration: 4800, time: 9.809 hours, Δt: 7.401 seconds, max|w|: 5.8e-03 m s⁻¹, wall time: 1.932 minutes
Iteration: 5000, time: 10.221 hours, Δt: 7.410 seconds, max|w|: 6.3e-03 m s⁻¹, wall time: 1.965 minutes
Iteration: 5200, time: 10.629 hours, Δt: 7.298 seconds, max|w|: 5.1e-03 m s⁻¹, wall time: 1.998 minutes
Iteration: 5400, time: 11.035 hours, Δt: 7.308 seconds, max|w|: 5.4e-03 m s⁻¹, wall time: 2.031 minutes
Iteration: 5600, time: 11.438 hours, Δt: 7.283 seconds, max|w|: 3.9e-03 m s⁻¹, wall time: 2.064 minutes
Iteration: 5800, time: 11.843 hours, Δt: 7.327 seconds, max|w|: 5.4e-03 m s⁻¹, wall time: 2.097 minutes
Iteration: 6000, time: 12.247 hours, Δt: 7.272 seconds, max|w|: 3.9e-03 m s⁻¹, wall time: 2.130 minutes
Iteration: 6200, time: 12.649 hours, Δt: 7.228 seconds, max|w|: 3.5e-03 m s⁻¹, wall time: 2.162 minutes
Iteration: 6400, time: 13.046 hours, Δt: 7.179 seconds, max|w|: 4.3e-03 m s⁻¹, wall time: 2.202 minutes
Iteration: 6600, time: 13.446 hours, Δt: 7.256 seconds, max|w|: 3.3e-03 m s⁻¹, wall time: 2.262 minutes
Iteration: 6800, time: 13.847 hours, Δt: 7.187 seconds, max|w|: 4.6e-03 m s⁻¹, wall time: 2.302 minutes
Iteration: 7000, time: 14.247 hours, Δt: 7.257 seconds, max|w|: 5.5e-03 m s⁻¹, wall time: 2.354 minutes
Iteration: 7200, time: 14.650 hours, Δt: 7.249 seconds, max|w|: 4.3e-03 m s⁻¹, wall time: 2.396 minutes
Iteration: 7400, time: 15.050 hours, Δt: 7.266 seconds, max|w|: 6.2e-03 m s⁻¹, wall time: 2.439 minutes
Iteration: 7600, time: 15.452 hours, Δt: 7.284 seconds, max|w|: 6.8e-03 m s⁻¹, wall time: 2.494 minutes
Iteration: 7800, time: 15.854 hours, Δt: 7.238 seconds, max|w|: 7.1e-03 m s⁻¹, wall time: 2.547 minutes
Iteration: 8000, time: 16.256 hours, Δt: 7.234 seconds, max|w|: 6.9e-03 m s⁻¹, wall time: 2.600 minutes
Iteration: 8200, time: 16.657 hours, Δt: 7.255 seconds, max|w|: 6.8e-03 m s⁻¹, wall time: 2.655 minutes
Iteration: 8400, time: 17.058 hours, Δt: 7.181 seconds, max|w|: 7.1e-03 m s⁻¹, wall time: 2.721 minutes
Iteration: 8600, time: 17.457 hours, Δt: 7.144 seconds, max|w|: 4.7e-03 m s⁻¹, wall time: 2.771 minutes
Iteration: 8800, time: 17.853 hours, Δt: 7.206 seconds, max|w|: 5.5e-03 m s⁻¹, wall time: 2.845 minutes
Iteration: 9000, time: 18.253 hours, Δt: 7.242 seconds, max|w|: 4.7e-03 m s⁻¹, wall time: 2.911 minutes
Iteration: 9200, time: 18.654 hours, Δt: 7.201 seconds, max|w|: 4.8e-03 m s⁻¹, wall time: 2.973 minutes
Iteration: 9400, time: 19.052 hours, Δt: 7.211 seconds, max|w|: 5.7e-03 m s⁻¹, wall time: 3.025 minutes
Iteration: 9600, time: 19.452 hours, Δt: 7.199 seconds, max|w|: 6.6e-03 m s⁻¹, wall time: 3.074 minutes
Iteration: 9800, time: 19.847 hours, Δt: 7.191 seconds, max|w|: 9.4e-03 m s⁻¹, wall time: 3.138 minutes
Iteration: 10000, time: 20.245 hours, Δt: 7.147 seconds, max|w|: 6.4e-03 m s⁻¹, wall time: 3.184 minutes
Iteration: 10200, time: 20.642 hours, Δt: 7.163 seconds, max|w|: 8.0e-03 m s⁻¹, wall time: 3.231 minutes
Iteration: 10400, time: 21.039 hours, Δt: 6.979 seconds, max|w|: 8.2e-03 m s⁻¹, wall time: 3.293 minutes
Iteration: 10600, time: 21.433 hours, Δt: 7.210 seconds, max|w|: 7.2e-03 m s⁻¹, wall time: 3.342 minutes
Iteration: 10800, time: 21.832 hours, Δt: 7.103 seconds, max|w|: 7.7e-03 m s⁻¹, wall time: 3.391 minutes
Iteration: 11000, time: 22.230 hours, Δt: 7.152 seconds, max|w|: 7.5e-03 m s⁻¹, wall time: 3.439 minutes
Iteration: 11200, time: 22.628 hours, Δt: 7.188 seconds, max|w|: 6.0e-03 m s⁻¹, wall time: 3.498 minutes
Iteration: 11400, time: 23.026 hours, Δt: 7.197 seconds, max|w|: 6.0e-03 m s⁻¹, wall time: 3.568 minutes
Iteration: 11600, time: 23.425 hours, Δt: 7.147 seconds, max|w|: 5.4e-03 m s⁻¹, wall time: 3.634 minutes
Iteration: 11800, time: 23.822 hours, Δt: 7.150 seconds, max|w|: 7.3e-03 m s⁻¹, wall time: 3.709 minutes
[ Info: Simulation is stopping after running for 3.497 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.