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 by Qᵀ = Qʰ / (ρ₀ * cᴾ), where Qʰ is the heat flux, ρ₀ is a reference density, and cᴾ is the heat capacity of seawater. With a reference density ρ₀ = 1026 kg m⁻³and heat capacity cᴾ = 3991, our chosen temperature flux of Qᵀ = 5 × 10⁻⁵ K m⁻¹ s⁻¹ corresponds to a heat flux of Qʰ = 204.7 W m⁻², a relatively powerful cooling rate.
A velocity flux Qᵘ at the top of the domain, which is related to the x momentum flux τˣ via τˣ = ρ₀ * Qᵘ, where ρ₀ is a reference density. Our chosen value of Qᵘ = -2 × 10⁻⁵ m² s⁻² roughly corresponds to atmospheric winds of uᵃ = 2.9 m s⁻¹ in the positive x-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, use closure = ConstantSmagorinsky() in the model constructor.
To change the architecture to GPU, replace the architecture keyword argument with architecture = 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.