Wind- and convection-driven mixing in an ocean surface boundary layer
This example simulates mixing by three-dimensional turbulence in an ocean surface boundary layer driven by atmospheric winds and convection. It demonstrates:
- How to set-up a grid with varying spacing in the vertical direction
- How to use the
SeawaterBuoyancymodel for buoyancy with a linear equation of state. - How to use a turbulence closure for large eddy simulation.
- How to use a function to impose a boundary condition.
Install dependencies
First let's make sure we have all required packages installed.
using Pkg
pkg"add Oceananigans, Plots"We start by importing all of the packages and functions that we'll need for this example.
using Random
using Printf
using Plots
using Oceananigans
using Oceananigans.Units: minute, minutes, hourThe grid
We use 32²×24 grid points with 2 m grid spacing in the horizontal and varying spacing in the vertical, with higher resolution closer to the surface. Here we use a stretching function for the vertical nodes that maintains relatively constant vertical spacing in the mixed layer, which is desirable from a numerical standpoint:
Nz = 24 # number of points in the vertical direction
Lz = 32 # (m) domain depth
refinement = 1.2 # controls spacing near surface (higher means finer spaced)
stretching = 12 # controls rate of stretching at bottom
# Normalized height ranging from 0 to 1
h(k) = (k - 1) / 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(size = (32, 32, Nz),
x = (0, 64),
y = (0, 64),
z = z_faces)32×32×24 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
├── Periodic x ∈ [0.0, 64.0) regularly spaced with Δx=2.0
├── Periodic y ∈ [0.0, 64.0) regularly spaced with Δy=2.0
└── Bounded z ∈ [-32.0, 0.0] variably spaced with min(Δz)=1.11123, max(Δz)=2.53571We plot vertical spacing versus depth to inspect the prescribed grid stretching:
plot(grid.Δzᵃᵃᶜ[1:grid.Nz], grid.zᵃᵃᶜ[1:grid.Nz],
marker = :circle,
ylabel = "Depth (m)",
xlabel = "Vertical spacing (m)",
legend = nothing)Buoyancy that depends on temperature and salinity
We use the SeawaterBuoyancy model with a linear equation of state,
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion = 2e-4,
haline_contraction = 8e-4))SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation of state: LinearEquationOfState(thermal_expansion=0.0002, haline_contraction=0.0008)Boundary conditions
We calculate the surface temperature flux associated with surface heating of 200 W m⁻², reference density ρₒ, and heat capacity cᴾ,
Qʰ = 200 # W m⁻², surface _heat_ flux
ρₒ = 1026 # kg m⁻³, average density at the surface of the world ocean
cᴾ = 3991 # J K⁻¹ kg⁻¹, typical heat capacity for seawater
Qᵀ = Qʰ / (ρₒ * cᴾ) # K m s⁻¹, surface _temperature_ flux4.884283985946938e-5Finally, we impose a temperature gradient dTdz both initially and at the bottom of the domain, culminating in the boundary conditions on temperature,
dTdz = 0.01 # K m⁻¹
T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵀ),
bottom = GradientBoundaryCondition(dTdz))Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition
├── east: DefaultBoundaryCondition
├── south: DefaultBoundaryCondition
├── north: DefaultBoundaryCondition
├── bottom: GradientBoundaryCondition: 0.01
├── top: FluxBoundaryCondition: 4.88428e-5
└── immersed: FluxBoundaryCondition: NothingNote that a positive temperature flux at the surface of the ocean implies cooling. This is because a positive temperature flux implies that temperature is fluxed upwards, out of the ocean.
For the velocity field, we imagine a wind blowing over the ocean surface with an average velocity at 10 meters u₁₀, and use a drag coefficient cᴰ to estimate the kinematic stress (that is, stress divided by density) exerted by the wind on the ocean:
u₁₀ = 10 # m s⁻¹, average wind velocity 10 meters above the ocean
cᴰ = 2.5e-3 # dimensionless drag coefficient
ρₐ = 1.225 # kg m⁻³, average density of air at sea-level
Qᵘ = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²-0.0002984892787524367The boundary conditions on u are thus
u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition
├── east: DefaultBoundaryCondition
├── south: DefaultBoundaryCondition
├── north: DefaultBoundaryCondition
├── bottom: DefaultBoundaryCondition
├── top: FluxBoundaryCondition: -0.000298489
└── immersed: FluxBoundaryCondition: NothingFor salinity, S, we impose an evaporative flux of the form
@inline Qˢ(x, y, t, S, evaporation_rate) = - evaporation_rate * S # [salinity unit] m s⁻¹where S is salinity. We use an evporation rate of 1 millimeter per hour,
evaporation_rate = 1e-3 / hour # m s⁻¹2.7777777777777776e-7We build the Flux evaporation BoundaryCondition with the function Qˢ, indicating that Qˢ depends on salinity S and passing the parameter evaporation_rate,
evaporation_bc = FluxBoundaryCondition(Qˢ, field_dependencies=:S, parameters=evaporation_rate)FluxBoundaryCondition: ContinuousBoundaryFunction Qˢ at (Nothing, Nothing, Nothing)The full salinity boundary conditions are
S_bcs = FieldBoundaryConditions(top=evaporation_bc)Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition
├── east: DefaultBoundaryCondition
├── south: DefaultBoundaryCondition
├── north: DefaultBoundaryCondition
├── bottom: DefaultBoundaryCondition
├── top: FluxBoundaryCondition: ContinuousBoundaryFunction Qˢ at (Nothing, Nothing, Nothing)
└── immersed: FluxBoundaryCondition: NothingModel instantiation
We fill in the final details of the model here: upwind-biased 5th-order advection for momentum and tracers, 3rd-order Runge-Kutta time-stepping, Coriolis forces, and the AnisotropicMinimumDissipation closure for large eddy simulation to model the effect of turbulent motions at scales smaller than the grid scale that we cannot explicitly resolve.
model = NonhydrostaticModel(advection = UpwindBiasedFifthOrder(),
timestepper = :RungeKutta3,
grid = grid,
tracers = (:T, :S),
coriolis = FPlane(f=1e-4),
buoyancy = buoyancy,
closure = AnisotropicMinimumDissipation(),
boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs))NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 32×32×24 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── tracers: (T, S)
├── closure: AnisotropicMinimumDissipation{ExplicitTimeDiscretization, NamedTuple{(:T, :S), Tuple{Float64, Float64}}, Float64, Nothing}
├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.0002, haline_contraction=0.0008) with -ĝ = ZDirection
└── coriolis: FPlane{Float64}(f=0.0001)Notes:
To use the Smagorinsky-Lilly turbulence closure (with a constant model coefficient) rather than
AnisotropicMinimumDissipation, useclosure = SmagorinskyLilly()in the model constructor.To change the
architecturetoGPU, replacearchitecture = CPU()witharchitecture = GPU().
Initial conditions
Our initial condition for temperature consists of a linear stratification superposed with random noise damped at the walls, while our initial condition for velocity consists only of random noise.
# Random noise damped at top and bottom
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise
# Temperature initial condition: a stable density gradient with random noise superposed.
Tᵢ(x, y, z) = 20 + dTdz * z + dTdz * model.grid.Lz * 1e-6 * Ξ(z)
# Velocity initial condition: random noise scaled by the friction velocity.
uᵢ(x, y, z) = sqrt(abs(Qᵘ)) * 1e-3 * Ξ(z)
# `set!` the `model` fields using functions or constants:
set!(model, u=uᵢ, w=uᵢ, T=Tᵢ, S=35)Setting up a simulation
We set-up a simulation with an initial time-step of 10 seconds that stops at 40 minutes, with adaptive time-stepping and progress printing.
simulation = Simulation(model, Δt=10.0, stop_time=40minutes)Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 10 seconds
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN years
├── Stop time: 40 minutes
├── 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 entriesThe TimeStepWizard helps ensure stable time-stepping with a Courant-Freidrichs-Lewy (CFL) number of 1.0.
wizard = TimeStepWizard(cfl=1.0, max_change=1.1, max_Δt=1minute)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))Callback of TimeStepWizard(cfl=1.0, max_Δt=60.0, min_Δt=0.0) on IterationInterval(10)Nice progress messaging is helpful:
# Print a progress message
progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, max(|w|) = %.1e ms⁻¹, wall time: %s\n",
iteration(sim),
prettytime(sim),
prettytime(sim.Δt),
maximum(abs, sim.model.velocities.w),
prettytime(sim.run_wall_time))
simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(10))Callback of progress_message on IterationInterval(10)We then set up the simulation:
Output
We use the JLD2OutputWriter to save $x, z$ slices of the velocity fields, tracer fields, and eddy diffusivities. The prefix keyword argument to JLD2OutputWriter indicates that output will be saved in ocean_wind_mixing_and_convection.jld2.
# Create a NamedTuple with eddy viscosity
eddy_viscosity = (; νₑ = model.diffusivity_fields.νₑ)
simulation.output_writers[:slices] =
JLD2OutputWriter(model, merge(model.velocities, model.tracers, eddy_viscosity),
prefix = "ocean_wind_mixing_and_convection",
indices = (:, grid.Ny/2, :),
schedule = TimeInterval(1minute),
force = true)JLD2OutputWriter scheduled on TimeInterval(1 minute):
├── filepath: ./ocean_wind_mixing_and_convection.jld2
├── 6 outputs: (u, v, w, T, S, νₑ)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiBWe're ready:
run!(simulation)[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 11 seconds, max(|w|) = 9.9e-06 ms⁻¹, wall time: 0 seconds
[ Info: ... simulation initialization complete (4.781 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (1.604 minutes).
Iteration: 0010, time: 1.733 minutes, Δt: 12.100 seconds, max(|w|) = 9.4e-06 ms⁻¹, wall time: 1.695 minutes
Iteration: 0020, time: 3.403 minutes, Δt: 13.310 seconds, max(|w|) = 8.1e-06 ms⁻¹, wall time: 1.710 minutes
Iteration: 0030, time: 5.444 minutes, Δt: 14.641 seconds, max(|w|) = 6.6e-06 ms⁻¹, wall time: 1.722 minutes
Iteration: 0040, time: 7.488 minutes, Δt: 16.105 seconds, max(|w|) = 6.5e-06 ms⁻¹, wall time: 1.732 minutes
Iteration: 0050, time: 10 minutes, Δt: 12.417 seconds, max(|w|) = 8.2e-06 ms⁻¹, wall time: 1.745 minutes
Iteration: 0060, time: 12 minutes, Δt: 10.350 seconds, max(|w|) = 9.6e-06 ms⁻¹, wall time: 1.755 minutes
Iteration: 0070, time: 13.517 minutes, Δt: 9.189 seconds, max(|w|) = 1.5e-05 ms⁻¹, wall time: 1.766 minutes
Iteration: 0080, time: 14.919 minutes, Δt: 8.326 seconds, max(|w|) = 4.2e-05 ms⁻¹, wall time: 1.776 minutes
Iteration: 0090, time: 16.139 minutes, Δt: 7.694 seconds, max(|w|) = 1.1e-04 ms⁻¹, wall time: 1.789 minutes
Iteration: 0100, time: 17.385 minutes, Δt: 7.134 seconds, max(|w|) = 3.1e-04 ms⁻¹, wall time: 1.799 minutes
Iteration: 0110, time: 18.476 minutes, Δt: 6.691 seconds, max(|w|) = 9.1e-04 ms⁻¹, wall time: 1.810 minutes
Iteration: 0120, time: 19.558 minutes, Δt: 6.252 seconds, max(|w|) = 2.6e-03 ms⁻¹, wall time: 1.822 minutes
Iteration: 0130, time: 20.521 minutes, Δt: 5.812 seconds, max(|w|) = 6.6e-03 ms⁻¹, wall time: 1.832 minutes
Iteration: 0140, time: 21.484 minutes, Δt: 5.314 seconds, max(|w|) = 1.8e-02 ms⁻¹, wall time: 1.843 minutes
Iteration: 0150, time: 22.354 minutes, Δt: 5.129 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 1.853 minutes
Iteration: 0160, time: 23.171 minutes, Δt: 4.907 seconds, max(|w|) = 4.6e-02 ms⁻¹, wall time: 1.865 minutes
Iteration: 0170, time: 23.989 minutes, Δt: 5.015 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 1.877 minutes
Iteration: 0180, time: 24.752 minutes, Δt: 4.652 seconds, max(|w|) = 7.4e-02 ms⁻¹, wall time: 1.889 minutes
Iteration: 0190, time: 25.465 minutes, Δt: 4.772 seconds, max(|w|) = 6.3e-02 ms⁻¹, wall time: 1.901 minutes
Iteration: 0200, time: 26.239 minutes, Δt: 5.249 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 1.911 minutes
Iteration: 0210, time: 27.087 minutes, Δt: 5.774 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 1.922 minutes
Iteration: 0220, time: 28 minutes, Δt: 6.351 seconds, max(|w|) = 8.1e-02 ms⁻¹, wall time: 1.932 minutes
Iteration: 0230, time: 29 minutes, Δt: 6.881 seconds, max(|w|) = 9.5e-02 ms⁻¹, wall time: 1.944 minutes
Iteration: 0240, time: 30.115 minutes, Δt: 7.519 seconds, max(|w|) = 8.0e-02 ms⁻¹, wall time: 1.954 minutes
Iteration: 0250, time: 31.251 minutes, Δt: 7.677 seconds, max(|w|) = 8.2e-02 ms⁻¹, wall time: 1.965 minutes
Iteration: 0260, time: 32.512 minutes, Δt: 7.851 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 1.977 minutes
Iteration: 0270, time: 33.785 minutes, Δt: 7.753 seconds, max(|w|) = 7.5e-02 ms⁻¹, wall time: 1.989 minutes
Iteration: 0280, time: 35 minutes, Δt: 7.296 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 2.003 minutes
Iteration: 0290, time: 36.122 minutes, Δt: 7.799 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 2.014 minutes
Iteration: 0300, time: 37.390 minutes, Δt: 7.526 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 2.027 minutes
Iteration: 0310, time: 38.627 minutes, Δt: 7.524 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 2.038 minutes
Iteration: 0320, time: 39.878 minutes, Δt: 7.228 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 2.049 minutes
[ Info: Simulation is stopping. Model time 40 minutes has hit or exceeded simulation stop time 40 minutes.Turbulence visualization
We animate the data saved in ocean_wind_mixing_and_convection.jld2. We prepare for animating the flow by loading the data into FieldTimeSeries and defining functions for computing colorbar limits.
filepath = "ocean_wind_mixing_and_convection.jld2"
time_series = (w = FieldTimeSeries(filepath, "w"),
T = FieldTimeSeries(filepath, "T"),
S = FieldTimeSeries(filepath, "S"),
νₑ = FieldTimeSeries(filepath, "νₑ"))
# Coordinate arrays
xw, yw, zw = nodes(time_series.w)
xT, yT, zT = nodes(time_series.T)
""" Returns colorbar levels equispaced between `(-clim, clim)` and encompassing the extrema of `c`. """
function divergent_levels(c, clim, nlevels=21)
cmax = maximum(abs, c)
levels = clim > cmax ? range(-clim, stop=clim, length=nlevels) : range(-cmax, stop=cmax, length=nlevels)
return (levels[1], levels[end]), levels
end
""" Returns colorbar levels equispaced between `clims` and encompassing the extrema of `c`."""
function sequential_levels(c, clims, nlevels=20)
levels = range(clims[1], stop=clims[2], length=nlevels)
cmin, cmax = minimum(c), maximum(c)
cmin < clims[1] && (levels = vcat([cmin], levels))
cmax > clims[2] && (levels = vcat(levels, [cmax]))
return clims, levels
endWe start the animation at t = 10minutes since things are pretty boring till then:
times = time_series.w.times
intro = searchsortedfirst(times, 10minutes)
anim = @animate for (i, t) in enumerate(times[intro:end])
@info "Drawing frame $i of $(length(times[intro:end]))..."
w = interior(time_series.w[i], :, 1, :)
T = interior(time_series.T[i], :, 1, :)
S = interior(time_series.S[i], :, 1, :)
νₑ = interior(time_series.νₑ[i], :, 1, :)
wlims, wlevels = divergent_levels(w, 2e-2)
Tlims, Tlevels = sequential_levels(T, (19.7, 19.99))
Slims, Slevels = sequential_levels(S, (35, 35.005))
νlims, νlevels = sequential_levels(νₑ, (1e-6, 5e-3))
kwargs = (linewidth=0, xlabel="x (m)", ylabel="z (m)", aspectratio=1,
xlims=(0, grid.Lx), ylims=(-grid.Lz, 0))
w_plot = contourf(xw, zw, w'; color=:balance, clims=wlims, levels=wlevels, kwargs...)
T_plot = contourf(xT, zT, T'; color=:thermal, clims=Tlims, levels=Tlevels, kwargs...)
S_plot = contourf(xT, zT, S'; color=:haline, clims=Slims, levels=Slevels, kwargs...)
# We use a heatmap for the eddy viscosity to observe how it varies on the grid scale.
ν_plot = heatmap(xT, zT, νₑ'; color=:thermal, clims=νlims, levels=νlevels, kwargs...)
w_title = @sprintf("vertical velocity (m s⁻¹), t = %s", prettytime(t))
T_title = "temperature (ᵒC)"
S_title = "salinity (g kg⁻¹)"
ν_title = "eddy viscosity (m² s⁻¹)"
# Arrange the plots side-by-side.
plot(w_plot, T_plot, S_plot, ν_plot, layout=(2, 2), size=(1200, 600),
title=[w_title T_title S_title ν_title])
endThis page was generated using Literate.jl.