Edit on GitHub

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;
  • how to use user-defined model.parameters in a boundary condition function.

In addition to Oceananigans.jl we need PyPlot 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 Oceananigans, Oceananigans.OutputWriters, Oceananigans.Diagnostics
using PyPlot, Random, Printf

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

  1. 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.

  2. 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ᵃ.

  3. 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 = 48       # 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⁻¹]
   end_time = 2hour    # End time for the simulation
          f = 1e-4     # Coriolis parameter
          α = 2e-4     # Thermal expansion coefficient
          β = 8e-4     # Haline contraction coefficient
0.0008

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.

u_bcs = HorizontallyPeriodicBCs(top = BoundaryCondition(Flux, Qᵘ))

T_bcs = HorizontallyPeriodicBCs(   top = BoundaryCondition(Flux, Qᵀ),
                                bottom = BoundaryCondition(Gradient, ∂T∂z))

# Salinity flux: Qˢ = - E * S
@inline Qˢ(i, j, grid, time, iter, U, C, p) = @inbounds -p.evaporation * C.S[i, j, 1]

S_bcs = HorizontallyPeriodicBCs(top = BoundaryCondition(Flux, Qˢ))
Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions
├── x: CoordinateBoundaryConditions{BoundaryCondition{Periodic,Nothing},BoundaryCondition{Periodic,Nothing}}
├── y: CoordinateBoundaryConditions{BoundaryCondition{Periodic,Nothing},BoundaryCondition{Periodic,Nothing}}
└── z: CoordinateBoundaryConditions{BoundaryCondition{Flux,Nothing},BoundaryCondition{Flux,typeof(Main.ex-ocean_wind_mixing_and_convection.Qˢ)}}

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. We also pass the evaporation rate to the container model.parameters for use in the boundary condition function that calculates the salinity flux.

model = Model(
         architecture = CPU(),
                 grid = RegularCartesianGrid(size=(Nz, Nz, Nz), length=(Δz*Nz, Δz*Nz, Δz*Nz)),
             coriolis = FPlane(f=f),
             buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=α, β=β)),
              closure = AnisotropicMinimumDissipation(),
  boundary_conditions = BoundaryConditions(u=u_bcs, T=T_bcs, S=S_bcs),
           parameters = (evaporation = evaporation,)
)
Oceananigans.Model on a CPU architecture (time = 0.000 ns, iteration = 0) 
├── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
├── tracers: (:T, :S)
├── closure: VerstappenAnisotropicMinimumDissipation{Float64,NamedTuple{(:T, :S),Tuple{Float64,Float64}}}
├── buoyancy: SeawaterBuoyancy{Float64,LinearEquationOfState{Float64}}
├── coriolis: FPlane{Float64}
├── output writers: OrderedCollections.OrderedDict{Symbol,Oceananigans.AbstractOutputWriter} with no entries
└── diagnostics: OrderedCollections.OrderedDict{Symbol,Oceananigans.AbstractDiagnostic} with no entries

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()`

Set makeplot = true to live-update a plot of vertical velocity, temperature, and salinity as the simulation runs.

makeplot = false
false

Initial conditions

Out 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 tradient with random noise superposed.
T₀(x, y, z) = 20 + ∂T∂z * z + ∂T∂z * model.grid.Lz * 1e-1 * Ξ(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.
field_writer = JLD2OutputWriter(model, FieldOutputs(fields_to_output); interval=hour/4,
                                prefix="ocean_wind_mixing_and_convection", force=true)

# Add the output writer to the models `output_writers`.
model.output_writers[:fields] = field_writer;
Oceananigans.OutputWriters.JLD2OutputWriter{Nothing,Float64,Dict{Symbol,Oceananigans.OutputWriters.FieldOutput{UnionAll,F} where F},typeof(Oceananigans.OutputWriters.noinit),Array{Symbol,1},Dict{Symbol,Any}}("./ocean_wind_mixing_and_convection.jld2", Dict{Symbol,Oceananigans.OutputWriters.FieldOutput{UnionAll,F} where F}(:T => Oceananigans.OutputWriters.FieldOutput{UnionAll,Field{Oceananigans.Cell,Oceananigans.Cell,Oceananigans.Cell,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}(Array, Field at (Oceananigans.Cell, Oceananigans.Cell, Oceananigans.Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
    ├── size: (48, 48, 48)
    └── domain: x ∈ [0.0, 48.0], y ∈ [0.0, 48.0], z ∈ [-48.0, 0.0]
),:w => Oceananigans.OutputWriters.FieldOutput{UnionAll,Field{Oceananigans.Cell,Oceananigans.Cell,Oceananigans.Face,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}(Array, Field at (Oceananigans.Cell, Oceananigans.Cell, Oceananigans.Face)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
    ├── size: (48, 48, 48)
    └── domain: x ∈ [0.0, 48.0], y ∈ [0.0, 48.0], z ∈ [-48.0, 0.0]
),:v => Oceananigans.OutputWriters.FieldOutput{UnionAll,Field{Oceananigans.Cell,Oceananigans.Face,Oceananigans.Cell,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}(Array, Field at (Oceananigans.Cell, Oceananigans.Face, Oceananigans.Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
    ├── size: (48, 48, 48)
    └── domain: x ∈ [0.0, 48.0], y ∈ [0.0, 48.0], z ∈ [-48.0, 0.0]
),:S => Oceananigans.OutputWriters.FieldOutput{UnionAll,Field{Oceananigans.Cell,Oceananigans.Cell,Oceananigans.Cell,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}(Array, Field at (Oceananigans.Cell, Oceananigans.Cell, Oceananigans.Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
    ├── size: (48, 48, 48)
    └── domain: x ∈ [0.0, 48.0], y ∈ [0.0, 48.0], z ∈ [-48.0, 0.0]
),:u => Oceananigans.OutputWriters.FieldOutput{UnionAll,Field{Oceananigans.Face,Oceananigans.Cell,Oceananigans.Cell,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}(Array, Field at (Oceananigans.Face, Oceananigans.Cell, Oceananigans.Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
    ├── size: (48, 48, 48)
    └── domain: x ∈ [0.0, 48.0], y ∈ [0.0, 48.0], z ∈ [-48.0, 0.0]
),:νₑ => Oceananigans.OutputWriters.FieldOutput{UnionAll,Field{Oceananigans.Cell,Oceananigans.Cell,Oceananigans.Cell,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}(Array, Field at (Oceananigans.Cell, Oceananigans.Cell, Oceananigans.Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
    ├── size: (48, 48, 48)
    └── domain: x ∈ [0.0, 48.0], y ∈ [0.0, 48.0], z ∈ [-48.0, 0.0]
)), 900.0, nothing, Oceananigans.OutputWriters.noinit, Symbol[:grid, :coriolis, :buoyancy, :closure], 0.0, 1, Inf, false, 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)
TimeStepWizard{Float64}(0.2, 0.02, 1.1, 0.5, 5.0, 1.0)

A diagnostic that returns the maximum absolute value of w by calling wmax(model):

wmax = FieldMaximum(abs, model.velocities.w);
Oceananigans.Diagnostics.FieldMaximum{Field{Oceananigans.Cell,Oceananigans.Cell,Oceananigans.Face,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},typeof(abs)}(abs, Field at (Oceananigans.Cell, Oceananigans.Cell, Oceananigans.Face)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
    ├── size: (48, 48, 48)
    └── domain: x ∈ [0.0, 48.0], y ∈ [0.0, 48.0], z ∈ [-48.0, 0.0]
)

We also create a figure and define a plotting function for live plotting of results.

fig, axs = subplots(ncols=3, figsize=(12, 5))

"""
    makeplot!(axs, model)

Make a triptych of x-z slices of vertical velocity, temperature, and salinity
associated with `model` in `axs`.
"""
function makeplot!(axs, model)
    jhalf = floor(Int, model.grid.Nz/2)

    # Coordinate arrays for plotting
    xC = repeat(model.grid.xC, 1, model.grid.Nz)
    zF = repeat(reshape(model.grid.zF[1:end-1], 1, model.grid.Nz), model.grid.Nx, 1)
    zC = repeat(reshape(model.grid.zC, 1, model.grid.Nz), model.grid.Nx, 1)

    sca(axs[1]); cla()
    title("Vertical velocity")
    pcolormesh(xC, zF, Array(interior(model.velocities.w))[:, jhalf, :])
    xlabel("\$ x \$ (m)"); ylabel("\$ z \$ (m)")

    sca(axs[2]); cla()
    title("Temperature")
    pcolormesh(xC, zC, Array(interior(model.tracers.T))[:, jhalf, :])
    xlabel("\$ x \$ (m)")

    sca(axs[3]); cla()
    title("Salinity")
    pcolormesh(xC, zC, Array(interior(model.tracers.S))[:, jhalf, :])
    xlabel("\$ x \$ (m)")

    [ax.set_aspect(1) for ax in axs]
    pause(0.01)

    return nothing
end
Main.ex-ocean_wind_mixing_and_convection.makeplot!

Finally, we run the the model in a while loop.

while model.clock.time < end_time

    # Update the time step associated with `wizard`.
    update_Δt!(wizard, model)

    # Time step the model forward
    walltime = @elapsed time_step!(model, 100, wizard.Δt)

    # 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))

    makeplot && makeplot!(axs, model)
end
i: 0100, t: 1.833 min, Δt: 1.100 s, wmax = 8.4e-03 ms⁻¹, wall time: 16.983 s
i: 0200, t: 3.850 min, Δt: 1.210 s, wmax = 1.3e-02 ms⁻¹, wall time: 9.585 s
i: 0300, t: 6.068 min, Δt: 1.331 s, wmax = 1.6e-02 ms⁻¹, wall time: 9.597 s
i: 0400, t: 8.508 min, Δt: 1.464 s, wmax = 1.5e-02 ms⁻¹, wall time: 9.691 s
i: 0500, t: 11.193 min, Δt: 1.611 s, wmax = 1.4e-02 ms⁻¹, wall time: 9.577 s
i: 0600, t: 14.145 min, Δt: 1.772 s, wmax = 1.5e-02 ms⁻¹, wall time: 9.626 s
i: 0700, t: 17.393 min, Δt: 1.949 s, wmax = 1.3e-02 ms⁻¹, wall time: 9.672 s
i: 0800, t: 20.966 min, Δt: 2.144 s, wmax = 1.9e-02 ms⁻¹, wall time: 9.654 s
i: 0900, t: 24.896 min, Δt: 2.358 s, wmax = 2.4e-02 ms⁻¹, wall time: 9.619 s
i: 1000, t: 29.219 min, Δt: 2.594 s, wmax = 2.3e-02 ms⁻¹, wall time: 9.625 s
i: 1100, t: 33.974 min, Δt: 2.853 s, wmax = 2.3e-02 ms⁻¹, wall time: 9.540 s
i: 1200, t: 39.205 min, Δt: 3.138 s, wmax = 2.9e-02 ms⁻¹, wall time: 9.519 s
i: 1300, t: 44.958 min, Δt: 3.452 s, wmax = 2.6e-02 ms⁻¹, wall time: 9.575 s
i: 1400, t: 51.287 min, Δt: 3.797 s, wmax = 2.5e-02 ms⁻¹, wall time: 9.526 s
i: 1500, t: 58.250 min, Δt: 4.177 s, wmax = 2.7e-02 ms⁻¹, wall time: 9.440 s
i: 1600, t: 1.088 hr, Δt: 4.233 s, wmax = 2.9e-02 ms⁻¹, wall time: 9.457 s
i: 1700, t: 1.208 hr, Δt: 4.312 s, wmax = 2.9e-02 ms⁻¹, wall time: 9.452 s
i: 1800, t: 1.321 hr, Δt: 4.073 s, wmax = 3.5e-02 ms⁻¹, wall time: 9.437 s
i: 1900, t: 1.420 hr, Δt: 3.567 s, wmax = 3.0e-02 ms⁻¹, wall time: 9.374 s
i: 2000, t: 1.529 hr, Δt: 3.923 s, wmax = 3.5e-02 ms⁻¹, wall time: 9.432 s
i: 2100, t: 1.649 hr, Δt: 4.316 s, wmax = 3.6e-02 ms⁻¹, wall time: 9.426 s
i: 2200, t: 1.767 hr, Δt: 4.232 s, wmax = 3.2e-02 ms⁻¹, wall time: 9.447 s
i: 2300, t: 1.885 hr, Δt: 4.253 s, wmax = 3.8e-02 ms⁻¹, wall time: 9.447 s
i: 2400, t: 1.999 hr, Δt: 4.107 s, wmax = 3.5e-02 ms⁻¹, wall time: 9.415 s
i: 2500, t: 2.107 hr, Δt: 3.899 s, wmax = 3.4e-02 ms⁻¹, wall time: 9.400 s

Show the reults in a plot

makeplot!(axs, model)
gcf()

This page was generated using Literate.jl.