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 use the
SeawaterBuoyancy
model 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;
In addition to Oceananigans.jl
we need Plots
for plotting, Random
for generating random initial conditions, and Printf
for printing progress messages. We also need Oceananigans.OutputWriters
and Oceananigans.Diagnostics
to access some nice features for writing output data to disk.
using Random, Printf, Plots
using Oceananigans, Oceananigans.OutputWriters, Oceananigans.Diagnostics, Oceananigans.Utils,
Oceananigans.BoundaryConditions, Oceananigans.Grids
Model parameters
Here we use an isotropic, cubic grid with Nz
grid points and grid spacing Δz = 1
meter. We specify fluxes of heat, momentum, and salinity via
A temperature flux
Qᵀ
at the top of the domain, which is related to heat flux byQᵀ = Qʰ / (ρ₀ * cᴾ)
, whereQʰ
is the heat flux,ρ₀
is a reference density, andcᴾ
is the heat capacity of seawater. With a reference densityρ₀ = 1026 kg m⁻³
and heat capacitycᴾ = 3991
, our chosen temperature flux ofQᵀ = 5 × 10⁻⁵ K m⁻¹ s⁻¹
corresponds to a heat flux ofQʰ = 204.7 W m⁻²
, a relatively powerful cooling rate.A velocity flux
Qᵘ
at the top of the domain, which is related to thex
momentum fluxτˣ
viaτˣ = ρ₀ * Qᵘ
, whereρ₀
is a reference density. Our chosen value ofQᵘ = -2 × 10⁻⁵ m² s⁻²
roughly corresponds to atmospheric winds ofuᵃ = 2.9 m s⁻¹
in the positivex
-direction, using the parameterizationτ = 0.0025 * |uᵃ| * uᵃ
.An evaporation rate
evaporation = 10⁻⁷ m s⁻¹
, or approximately 0.1 millimeter per hour.
Finally, we use an initial temperature gradient of ∂T/∂z = 0.005 K m⁻¹
, which implies an iniital buoyancy frequency N² = α * g * ∂T/∂z = 9.8 × 10⁻⁶ s⁻²
with a thermal expansion coefficient α = 2 × 10⁻⁴ K⁻¹
and gravitational acceleration g = 9.81 s⁻²
. Note that, by default, the SeawaterBuoyancy
model uses a gravitational acceleration gᴱᵃʳᵗʰ = 9.80665 s⁻²
.
Nz = 32 # Number of grid points in x, y, z
Δz = 1.0 # Grid spacing in x, y, z (meters)
Qᵀ = 5e-5 # Temperature flux at surface
Qᵘ = -2e-5 # Velocity flux at surface
∂T∂z = 0.005 # Initial vertical temperature gradient
evaporation = 1e-7 # Mass-specific evaporation rate [m s⁻¹]
f = 1e-4 # Coriolis parameter
α = 2e-4 # Thermal expansion coefficient
β = 8e-4 # Haline contraction coefficient
Boundary conditions
Here we define Flux
boundary conditions at the surface for u
, T
, and S
, and a Gradient
boundary condition on T
that maintains a constant stratification at the bottom. Our flux boundary condition for salinity uses a function that calculates the salinity flux in terms of the evaporation rate.
grid = RegularCartesianGrid(size=(Nz, Nz, Nz), extent=(Δz*Nz, Δz*Nz, Δz*Nz))
u_bcs = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ))
T_bcs = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵀ),
bottom = BoundaryCondition(Gradient, ∂T∂z))
# Salinity flux: Qˢ = - E * S
@inline Qˢ(i, j, grid, clock, model_fields, evaporation) = @inbounds - evaporation * model_fields.S[i, j, 1]
S_bcs = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, Qˢ, discrete_form=true, parameters=evaporation))
Model instantiation
We instantiate a horizontally-periodic Model
on the CPU with on a RegularCartesianGrid
, using a FPlane
model for rotation (constant rotation rate), a linear equation of state for temperature and salinity, the Anisotropic Minimum Dissipation closure to model the effects of unresolved turbulence, and the previously defined boundary conditions for u
, T
, and S
.
model = IncompressibleModel( architecture = CPU(),
grid = grid,
coriolis = FPlane(f=f),
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=α, β=β)),
closure = AnisotropicMinimumDissipation(),
boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs))
Notes:
To use the Smagorinsky-Lilly turbulence closure (with a constant model coefficient) rather than
AnisotropicMinimumDissipation
, useclosure = ConstantSmagorinsky()
in the model constructor.To change the
architecture
toGPU
, replace thearchitecture
keyword argument 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 + ∂T∂z * z + ∂T∂z * model.grid.Lz * 1e-6 * Ξ(z)
# Velocity initial condition: random noise scaled by the friction velocity.
u₀(x, y, z) = sqrt(abs(Qᵘ)) * 1e-1 * Ξ(z)
set!(model, u=u₀, w=u₀, T=T₀, S=35)
Set up output
We set up an output writer that saves all velocity fields, tracer fields, and the subgrid turbulent diffusivity associated with model.closure
. The prefix
keyword argument to JLD2OutputWriter
indicates that output will be saved in ocean_wind_mixing_and_convection.jld2
.
# Create a NamedTuple containing all the fields to be outputted.
fields_to_output = merge(model.velocities, model.tracers, (νₑ=model.diffusivities.νₑ,))
# Instantiate a JLD2OutputWriter to write fields. We will add it to the simulation before
# running it.
field_writer = JLD2OutputWriter(model, fields_to_output; time_interval=hour/4,
prefix="ocean_wind_mixing_and_convection", force=true)
JLD2OutputWriter{Nothing,Float64,NamedTuple{(:u, :v, :w, :T, :S, :νₑ),Tuple{Field{Face,Cell,Cell,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},NamedTuple{(:x, :y, :z),Tuple{CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{Flux,Nothing},BoundaryCondition{Flux,Float64}}}}},Field{Cell,Face,Cell,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},NamedTuple{(:x, :y, :z),Tuple{CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{Flux,Nothing},BoundaryCondition{Flux,Nothing}}}}},Field{Cell,Cell,Face,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},NamedTuple{(:x, :y, :z),Tuple{CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{NormalFlow,Nothing},BoundaryCondition{NormalFlow,Nothing}}}}},Field{Cell,Cell,Cell,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},NamedTuple{(:x, :y, :z),Tuple{CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{Gradient,Float64},BoundaryCondition{Flux,Float64}}}}},Field{Cell,Cell,Cell,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},NamedTuple{(:x, :y, :z),Tuple{CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{Flux,Nothing},BoundaryCondition{Flux,Oceananigans.BoundaryConditions.ParameterizedDiscreteBoundaryFunction{typeof(Main.ex-ocean_wind_mixing_and_convection.Qˢ),Float64}}}}}},Field{Cell,Cell,Cell,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},NamedTuple{(:x, :y, :z),Tuple{CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}},CoordinateBoundaryConditions{BoundaryCondition{Flux,Nothing},BoundaryCondition{Flux,Nothing}}}}}}},FieldSlicer{Colon,Colon,Colon,Bool},UnionAll,typeof(Oceananigans.OutputWriters.noinit),Array{Symbol,1},Dict{Symbol,Any}}("./ocean_wind_mixing_and_convection.jld2", (u = Field located at (Face, Cell, Cell) ├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (34, 34, 34) ├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=32, Ny=32, Nz=32) └── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=Flux), v = Field located at (Cell, Face, Cell) ├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (34, 34, 34) ├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=32, Ny=32, Nz=32) └── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux), w = Field located at (Cell, Cell, Face) ├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (34, 34, 35) ├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=32, Ny=32, Nz=32) └── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=NormalFlow, top=NormalFlow), T = Field located at (Cell, Cell, Cell) ├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (34, 34, 34) ├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=32, Ny=32, Nz=32) └── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=Gradient, top=Flux), S = Field located at (Cell, Cell, Cell) ├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (34, 34, 34) ├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=32, Ny=32, Nz=32) └── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=Flux), νₑ = Field located at (Cell, Cell, Cell) ├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (34, 34, 34) ├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=32, Ny=32, Nz=32) └── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux)), nothing, 900.0, FieldSlicer{Colon,Colon,Colon,Bool}(Colon(), Colon(), Colon(), false), Array{Float32,N} where N, Oceananigans.OutputWriters.noinit, [:grid, :coriolis, :buoyancy, :closure], 0.0, 1, Inf, true, false, Dict{Symbol,Any}())
Running the simulation
To run the simulation, we instantiate a TimeStepWizard
to ensure stable time-stepping with a Courant-Freidrichs-Lewy (CFL) number of 0.2.
wizard = TimeStepWizard(cfl=0.2, Δt=1.0, max_change=1.1, max_Δt=5.0)
A diagnostic that returns the maximum absolute value of w
by calling wmax(model)
:
wmax = FieldMaximum(abs, model.velocities.w)
Finally, we set up and run the the simulation.
simulation = Simulation(model, Δt=wizard, stop_iteration=0, iteration_interval=10)
simulation.output_writers[:fields] = field_writer
anim = @animate for i in 1:100
# Run the simulation forward
simulation.stop_iteration += 10
walltime = @elapsed run!(simulation)
# Print a progress message
@printf("i: %04d, t: %s, Δt: %s, wmax = %.1e ms⁻¹, wall time: %s\n",
model.clock.iteration, prettytime(model.clock.time), prettytime(wizard.Δt),
wmax(model), prettytime(walltime))
# Coordinate arrays for plotting
xC, zF, zC = xnodes(Cell, grid)[:], znodes(Face, grid)[:], znodes(Cell, grid)[:]
# Slices to plots.
jhalf = ceil(Int, model.grid.Ny/2)
w = Array(interior(model.velocities.w))[:, jhalf, :]
T = Array(interior(model.tracers.T))[:, jhalf, :]
S = Array(interior(model.tracers.S))[:, jhalf, :]
# Plot the slices.
w_plot = heatmap(xC, zF, w', xlabel="x (m)", ylabel="z (m)", color=:balance, clims=(-3e-2, 3e-2))
T_plot = heatmap(xC, zC, T', xlabel="x (m)", ylabel="z (m)", color=:thermal, clims=(19.75, 20.0))
S_plot = heatmap(xC, zC, S', xlabel="x (m)", ylabel="z (m)", color=:haline, clims=(34.99, 35.01))
# Arrange the plots side-by-side.
plot(w_plot, T_plot, S_plot, layout=(1, 3), size=(1600, 400),
title=["vertical velocity (m/s)" "temperature (C)" "salinity (g/kg)"])
end
This page was generated using Literate.jl.