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.

Install dependencies

First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, JLD2, Plots"

  Resolving package versions...
No Changes to `/storage7/buildkite-agent/builds/tartarus-mit-edu-6/clima/oceananigans/docs/Project.toml`
No Changes to `/storage7/buildkite-agent/builds/tartarus-mit-edu-6/clima/oceananigans/docs/Manifest.toml`

We start by importing all of the packages and functions that we'll need for this example.

using Random
using Printf
using Plots
using JLD2

using Oceananigans
using Oceananigans.Utils

using Oceananigans.Grids: nodes
using Oceananigans.Advection: UpwindBiasedFifthOrder
using Oceananigans.Diagnostics: FieldMaximum
using Oceananigans.OutputWriters: JLD2OutputWriter, FieldSlicer, TimeInterval

The grid

We use 32³ grid points with 2 m grid spacing in the horizontal and 1 m spacing in the vertical,

grid = RegularCartesianGrid(size=(32, 32, 32), extent=(64, 64, 32))

RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}
                   domain: x ∈ [0.0, 64.0], y ∈ [0.0, 64.0], z ∈ [-32.0, 0.0]
                 topology: (Periodic, Periodic, Bounded)
  resolution (Nx, Ny, Nz): (32, 32, 32)
   halo size (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (2.0, 2.0, 1.0)

Buoyancy that depends on temperature and salinity

We use the SeawaterBuoyancy model with a linear equation of state,

buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=2e-4, β=8e-4))

SeawaterBuoyancy{Float64}: g = 9.80665
└── equation of state: LinearEquationOfState{Float64}: α = 2.00e-04, β = 8.00e-04

where $α$ and $β$ are the thermal expansion and haline contraction coefficients for temperature and salinity.

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⁻¹ s⁻¹, typical heat capacity for seawater

Qᵀ = Qʰ / (ρₒ * cᴾ) # K m⁻¹ s⁻¹, surface _temperature_ flux

4.884283985946938e-5

Finally, 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 = TracerBoundaryConditions(grid,
                                 top = BoundaryCondition(Flux, Qᵀ),
                                 bottom = BoundaryCondition(Gradient, dTdz))

Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions
├── x: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
├── y: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
└── z: CoordinateBoundaryConditions{BoundaryCondition{Gradient,Float64},BoundaryCondition{Flux,Float64}}

Note 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.0002984892787524367

The boundary conditions on u are thus

u_bcs = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ))

Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions
├── x: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
├── y: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
└── z: CoordinateBoundaryConditions{BoundaryCondition{Flux,Nothing},BoundaryCondition{Flux,Float64}}

For salinity, S, we impose an evaporative flux of the form

@inline Qˢ(x, y, t, S, evaporation_rate) = - evaporation_rate * S

where S is salinity. We use an evporation rate of 1 millimeter per hour,

evaporation_rate = 1e-3 / hour

2.7777777777777776e-7

We build the Flux evaporation BoundaryCondition with the function Qˢ, indicating that Qˢ depends on salinity S and passing the parameter evaporation_rate,

evaporation_bc = BoundaryCondition(Flux, Qˢ, field_dependencies=:S, parameters=evaporation_rate)

BoundaryCondition: type=Flux, condition=Qˢ(x, y, t, S, evaporation_rate) in Main.ex-ocean_wind_mixing_and_convection at none:1

The full salinity boundary conditions are

S_bcs = TracerBoundaryConditions(grid, top=evaporation_bc)

Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions
├── x: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
├── y: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
└── z: CoordinateBoundaryConditions{BoundaryCondition{Flux,Nothing},BoundaryCondition{Flux,Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Nothing,Nothing,Nothing,Nothing,typeof(Main.ex-ocean_wind_mixing_and_convection.Qˢ),Float64,Tuple{Symbol},Nothing,Nothing}}}

Model 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 = IncompressibleModel(architecture = CPU(),
                            advection = UpwindBiasedFifthOrder(),
                            timestepper = :RungeKutta3,
                            grid = grid,
                            coriolis = FPlane(f=1e-4),
                            buoyancy = buoyancy,
                            closure = AnisotropicMinimumDissipation(),
                            boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs))

IncompressibleModel{CPU, Float64}(time = 0 seconds, iteration = 0) 
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=32, Ny=32, Nz=32)
├── tracers: (:T, :S)
├── closure: VerstappenAnisotropicMinimumDissipation{Float64,NamedTuple{(:T, :S),Tuple{Float64,Float64}},Float64,NamedTuple{(:T, :S),Tuple{Float64,Float64}}}
├── buoyancy: SeawaterBuoyancy{Float64,LinearEquationOfState{Float64},Nothing,Nothing}
└── coriolis: FPlane{Float64}

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 architecture = CPU() 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 + 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 first build a TimeStepWizard to ensure stable time-stepping with a Courant-Freidrichs-Lewy (CFL) number of 1.0.

wizard = TimeStepWizard(cfl=1.0, Δt=10.0, max_change=1.1, max_Δt=1minute)

TimeStepWizard{Float64}(1.0, Inf, 1.1, 0.5, 60.0, 0.0, 10.0)

Nice progress messaging is helpful:

wmax = FieldMaximum(abs, model.velocities.w)

start_time = time_ns() # so we can print the total elapsed wall time

# Print a progress message
progress_message(sim) =
    @printf("i: %04d, t: %s, Δt: %s, wmax = %.1e ms⁻¹, wall time: %s\n",
            sim.model.clock.iteration, prettytime(model.clock.time),
            prettytime(wizard.Δt), wmax(sim.model),
            prettytime((time_ns() - start_time) * 1e-9))

progress_message (generic function with 1 method)

We then set up the simulation:

simulation = Simulation(model, Δt=wizard, stop_time=40minutes, iteration_interval=10,
                        progress=progress_message)

Simulation{IncompressibleModel{CPU, Float64}}
├── Model clock: time = 0 seconds, iteration = 0 
├── Next time step (TimeStepWizard{Float64}): 10 seconds 
├── Iteration interval: 10
├── Stop criteria: Any[Oceananigans.Simulations.iteration_limit_exceeded, Oceananigans.Simulations.stop_time_exceeded, Oceananigans.Simulations.wall_time_limit_exceeded]
├── Run time: 0 seconds, wall time limit: Inf
├── Stop time: 40 minutes, stop iteration: Inf
├── Diagnostics: OrderedCollections.OrderedDict with 1 entry:
│   └── nan_checker => NaNChecker
└── Output writers: OrderedCollections.OrderedDict with no entries

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.diffusivities.νₑ,)

simulation.output_writers[:slices] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers, eddy_viscosity),
                           prefix = "ocean_wind_mixing_and_convection",
                     field_slicer = FieldSlicer(j=Int(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, :νₑ)
├── field slicer: FieldSlicer(:, 16, :, with_halos=false)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB

We're ready:

run!(simulation)

i: 0010, t: 1.667 minutes, Δt: 10 seconds, wmax = 9.0e-06 ms⁻¹, wall time: 16.252 seconds
i: 0020, t: 3 minutes, Δt: 11 seconds, wmax = 8.1e-06 ms⁻¹, wall time: 17.127 seconds
i: 0030, t: 5 minutes, Δt: 12.100 seconds, wmax = 6.1e-06 ms⁻¹, wall time: 18.010 seconds
i: 0040, t: 7 minutes, Δt: 13.310 seconds, wmax = 5.6e-06 ms⁻¹, wall time: 18.910 seconds
i: 0050, t: 9 minutes, Δt: 14.641 seconds, wmax = 7.5e-06 ms⁻¹, wall time: 19.810 seconds
i: 0060, t: 11 minutes, Δt: 12.417 seconds, wmax = 8.9e-06 ms⁻¹, wall time: 20.718 seconds
i: 0070, t: 12.677 minutes, Δt: 10.162 seconds, wmax = 2.8e-05 ms⁻¹, wall time: 21.636 seconds
i: 0080, t: 14 minutes, Δt: 8.818 seconds, wmax = 8.8e-05 ms⁻¹, wall time: 22.519 seconds
i: 0090, t: 15.266 minutes, Δt: 7.983 seconds, wmax = 2.8e-04 ms⁻¹, wall time: 23.343 seconds
i: 0100, t: 16.366 minutes, Δt: 7.314 seconds, wmax = 8.6e-04 ms⁻¹, wall time: 24.193 seconds
i: 0110, t: 17.453 minutes, Δt: 6.795 seconds, wmax = 2.7e-03 ms⁻¹, wall time: 25.012 seconds
i: 0120, t: 18.419 minutes, Δt: 6.290 seconds, wmax = 7.4e-03 ms⁻¹, wall time: 25.838 seconds
i: 0130, t: 19.389 minutes, Δt: 5.830 seconds, wmax = 1.6e-02 ms⁻¹, wall time: 26.688 seconds
i: 0140, t: 20.268 minutes, Δt: 5.356 seconds, wmax = 2.9e-02 ms⁻¹, wall time: 27.537 seconds
i: 0150, t: 21.084 minutes, Δt: 5.018 seconds, wmax = 4.2e-02 ms⁻¹, wall time: 28.372 seconds
i: 0160, t: 21.926 minutes, Δt: 5.051 seconds, wmax = 5.9e-02 ms⁻¹, wall time: 29.285 seconds
i: 0170, t: 22.787 minutes, Δt: 5.248 seconds, wmax = 7.5e-02 ms⁻¹, wall time: 30.124 seconds
i: 0180, t: 23.617 minutes, Δt: 5.292 seconds, wmax = 6.7e-02 ms⁻¹, wall time: 30.960 seconds
i: 0190, t: 24.457 minutes, Δt: 5.481 seconds, wmax = 7.6e-02 ms⁻¹, wall time: 31.796 seconds
i: 0200, t: 25.272 minutes, Δt: 5.440 seconds, wmax = 7.2e-02 ms⁻¹, wall time: 32.635 seconds
i: 0210, t: 26.091 minutes, Δt: 5.433 seconds, wmax = 5.8e-02 ms⁻¹, wall time: 33.515 seconds
i: 0220, t: 26.953 minutes, Δt: 5.177 seconds, wmax = 7.0e-02 ms⁻¹, wall time: 34.308 seconds
i: 0230, t: 27.854 minutes, Δt: 5.694 seconds, wmax = 7.4e-02 ms⁻¹, wall time: 35.156 seconds
i: 0240, t: 28.827 minutes, Δt: 6.206 seconds, wmax = 7.9e-02 ms⁻¹, wall time: 36.011 seconds
i: 0250, t: 29.910 minutes, Δt: 6.826 seconds, wmax = 6.9e-02 ms⁻¹, wall time: 36.972 seconds
i: 0260, t: 31 minutes, Δt: 7.395 seconds, wmax = 6.6e-02 ms⁻¹, wall time: 37.859 seconds
i: 0270, t: 32.125 minutes, Δt: 7.473 seconds, wmax = 6.2e-02 ms⁻¹, wall time: 38.697 seconds
i: 0280, t: 33.254 minutes, Δt: 7.620 seconds, wmax = 5.6e-02 ms⁻¹, wall time: 39.545 seconds
i: 0290, t: 34.525 minutes, Δt: 7.878 seconds, wmax = 6.1e-02 ms⁻¹, wall time: 40.386 seconds
i: 0300, t: 35.749 minutes, Δt: 7.486 seconds, wmax = 6.2e-02 ms⁻¹, wall time: 41.230 seconds
i: 0310, t: 36.875 minutes, Δt: 7.502 seconds, wmax = 6.0e-02 ms⁻¹, wall time: 42.071 seconds
i: 0320, t: 37.914 minutes, Δt: 6.857 seconds, wmax = 5.5e-02 ms⁻¹, wall time: 42.917 seconds
i: 0330, t: 39 minutes, Δt: 7.125 seconds, wmax = 5.2e-02 ms⁻¹, wall time: 43.809 seconds
i: 0339, t: 40 minutes, Δt: 7.112 seconds, wmax = 5.9e-02 ms⁻¹, wall time: 44.584 seconds
[ 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 creating coordinate arrays, opening the file, building a vector of the iterations that we saved data at, and defining functions for computing colorbar limits:

# Coordinate arrays
xw, yw, zw = nodes(model.velocities.w)
xT, yT, zT = nodes(model.tracers.T)

# Open the file with our data
file = jldopen(simulation.output_writers[:slices].filepath)

# Extract a vector of iterations
iterations = parse.(Int, keys(file["timeseries/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
end

Main.ex-ocean_wind_mixing_and_convection.sequential_levels

We start the animation at t = 10minutes since things are pretty boring till then:

times = [file["timeseries/t/$iter"] for iter in iterations]
intro = searchsortedfirst(times, 10minutes)

anim = @animate for (i, iter) in enumerate(iterations[intro:end])

    @info "Drawing frame $i from iteration $iter..."

    t = file["timeseries/t/$iter"]
    w = file["timeseries/w/$iter"][:, 1, :]
    T = file["timeseries/T/$iter"][:, 1, :]
    S = file["timeseries/S/$iter"][:, 1, :]
    νₑ = file["timeseries/νₑ/$iter"][:, 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])

    iter == iterations[end] && close(file)
end

This page was generated using Literate.jl.