# 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
`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, CairoMakie"
```

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

```
using Random
using Printf
using CairoMakie
using Oceananigans
using Oceananigans.Units: minute, minutes, hour
```

## The 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:

```
Nx = Ny = 32 # number of points in each of horizontal directions
Nz = 24 # number of points in the vertical direction
Lx = Ly = 64 # (m) domain horizontal extents
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 = (Nx, Nx, Nz),
x = (0, Lx),
y = (0, Ly),
z = z_faces)
```

```
32×32×24 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 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.53571
```

We plot vertical spacing versus depth to inspect the prescribed grid stretching:

```
fig = Figure(size=(1200, 800))
ax = Axis(fig[1, 1], ylabel = "Depth (m)", xlabel = "Vertical spacing (m)")
lines!(ax, zspacings(grid, Center()), znodes(grid, Center()))
scatter!(ax, zspacings(grid, Center()), znodes(grid, Center()))
fig
```

## 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 cooling of 200 W m⁻², reference density `ρₒ`

, and heat capacity `cᴾ`

,

```
Qʰ = 200.0 # W m⁻², surface _heat_ flux
ρₒ = 1026.0 # kg m⁻³, average density at the surface of the world ocean
cᴾ = 3991.0 # J K⁻¹ kg⁻¹, 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 = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵀ),
bottom = GradientBoundaryCondition(dTdz))
```

```
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: GradientBoundaryCondition: 0.01
├── top: FluxBoundaryCondition: 4.88428e-5
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
```

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 = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))`

```
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: -0.000298489
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
```

For 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-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 = 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 (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: ContinuousBoundaryFunction Qˢ at (Nothing, Nothing, Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: 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 = NonhydrostaticModel(; grid, buoyancy,
advection = UpwindBiasedFifthOrder(),
timestepper = :RungeKutta3,
tracers = (:T, :S),
coriolis = FPlane(f=1e-4),
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 ĝ = NegativeZDirection()
└── coriolis: FPlane{Float64}(f=0.0001)
```

Notes:

To use the Smagorinsky-Lilly turbulence closure (with a constant model coefficient) rather than

`AnisotropicMinimumDissipation`

, use`closure = SmagorinskyLilly()`

in the model constructor.To change the architecture to

`GPU`

, replace`CPU()`

with`GPU()`

inside the`grid`

constructor.

## 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 days
├── 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 entries
```

The `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))
add_callback!(simulation, progress_message, IterationInterval(20))
```

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.νₑ)
filename = "ocean_wind_mixing_and_convection"
simulation.output_writers[:slices] =
JLD2OutputWriter(model, merge(model.velocities, model.tracers, eddy_viscosity),
filename = filename * ".jld2",
indices = (:, grid.Ny/2, :),
schedule = TimeInterval(1minute),
overwrite_existing = true)
```

```
JLD2OutputWriter scheduled on TimeInterval(1 minute):
├── filepath: ./ocean_wind_mixing_and_convection.jld2
├── 6 outputs: (u, v, w, T, S, νₑ)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB
```

We're ready:

`run!(simulation)`

```
[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 11 seconds, max(|w|) = 1.1e-05 ms⁻¹, wall time: 0 seconds
[ Info: ... simulation initialization complete (5.651 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (7.828 seconds).
Iteration: 0020, time: 3.403 minutes, Δt: 13.310 seconds, max(|w|) = 8.9e-06 ms⁻¹, wall time: 14.840 seconds
Iteration: 0040, time: 7.488 minutes, Δt: 16.105 seconds, max(|w|) = 5.9e-06 ms⁻¹, wall time: 16.481 seconds
Iteration: 0060, time: 12 minutes, Δt: 9.990 seconds, max(|w|) = 9.0e-06 ms⁻¹, wall time: 18.103 seconds
Iteration: 0080, time: 14.884 minutes, Δt: 7.990 seconds, max(|w|) = 2.8e-05 ms⁻¹, wall time: 19.746 seconds
Iteration: 0100, time: 17.122 minutes, Δt: 6.903 seconds, max(|w|) = 9.5e-05 ms⁻¹, wall time: 21.364 seconds
Iteration: 0120, time: 19.215 minutes, Δt: 6.104 seconds, max(|w|) = 4.8e-04 ms⁻¹, wall time: 22.937 seconds
Iteration: 0140, time: 21 minutes, Δt: 5.515 seconds, max(|w|) = 2.2e-03 ms⁻¹, wall time: 24.557 seconds
Iteration: 0160, time: 22.781 minutes, Δt: 4.908 seconds, max(|w|) = 9.7e-03 ms⁻¹, wall time: 26.148 seconds
Iteration: 0180, time: 24.322 minutes, Δt: 4.666 seconds, max(|w|) = 2.4e-02 ms⁻¹, wall time: 27.761 seconds
Iteration: 0200, time: 25.830 minutes, Δt: 4.407 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 29.368 seconds
Iteration: 0220, time: 27.228 minutes, Δt: 4.805 seconds, max(|w|) = 6.8e-02 ms⁻¹, wall time: 31.007 seconds
Iteration: 0240, time: 28.797 minutes, Δt: 5.246 seconds, max(|w|) = 1.0e-01 ms⁻¹, wall time: 32.621 seconds
Iteration: 0260, time: 30.479 minutes, Δt: 5.727 seconds, max(|w|) = 1.1e-01 ms⁻¹, wall time: 34.236 seconds
Iteration: 0280, time: 32.280 minutes, Δt: 6.133 seconds, max(|w|) = 8.7e-02 ms⁻¹, wall time: 35.872 seconds
Iteration: 0300, time: 34.215 minutes, Δt: 5.981 seconds, max(|w|) = 7.7e-02 ms⁻¹, wall time: 37.536 seconds
Iteration: 0320, time: 36.202 minutes, Δt: 6.654 seconds, max(|w|) = 7.6e-02 ms⁻¹, wall time: 39.195 seconds
Iteration: 0340, time: 38.340 minutes, Δt: 6.376 seconds, max(|w|) = 6.7e-02 ms⁻¹, wall time: 40.825 seconds
[ Info: Simulation is stopping after running for 42.248 seconds.
[ Info: Simulation time 40 minutes equals or exceeds 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 = filename * ".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)
```

`([1.0, 3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0, 23.0, 25.0, 27.0, 29.0, 31.0, 33.0, 35.0, 37.0, 39.0, 41.0, 43.0, 45.0, 47.0, 49.0, 51.0, 53.0, 55.0, 57.0, 59.0, 61.0, 63.0], [1.0, 3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0, 23.0, 25.0, 27.0, 29.0, 31.0, 33.0, 35.0, 37.0, 39.0, 41.0, 43.0, 45.0, 47.0, 49.0, 51.0, 53.0, 55.0, 57.0, 59.0, 61.0, 63.0], [-30.73214655800215, -28.344120885446316, -26.245517670529452, -24.406267643406927, -22.77515636448867, -21.301119444298777, -19.9410533865803, -18.66145654635607, -17.43754626408272, -16.251584220508228, -15.09116977761991, -13.947785764082138, -12.815662268579008, -11.690933327108011, -10.57103192002591, -9.45426630377412, -8.339528643431164, -7.226097424980429, -6.1135049573687645, -5.001449332843377, -3.889736372710594, -2.7782415919634484, -1.666885417004453, -0.5556171156552487])`

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

```
times = time_series.w.times
intro = searchsortedfirst(times, 10minutes)
```

`11`

We are now ready to animate using Makie. We use Makie's `Observable`

to animate the data. To dive into how `Observable`

s work we refer to Makie.jl's Documentation.

```
n = Observable(intro)
wₙ = @lift interior(time_series.w[$n], :, 1, :)
Tₙ = @lift interior(time_series.T[$n], :, 1, :)
Sₙ = @lift interior(time_series.S[$n], :, 1, :)
νₑₙ = @lift interior(time_series.νₑ[$n], :, 1, :)
fig = Figure(size = (1000, 500))
axis_kwargs = (xlabel="x (m)",
ylabel="z (m)",
aspect = AxisAspect(grid.Lx/grid.Lz),
limits = ((0, grid.Lx), (-grid.Lz, 0)))
ax_w = Axis(fig[2, 1]; title = "Vertical velocity", axis_kwargs...)
ax_T = Axis(fig[2, 3]; title = "Temperature", axis_kwargs...)
ax_S = Axis(fig[3, 1]; title = "Salinity", axis_kwargs...)
ax_νₑ = Axis(fig[3, 3]; title = "Eddy viscocity", axis_kwargs...)
title = @lift @sprintf("t = %s", prettytime(times[$n]))
wlims = (-0.05, 0.05)
Tlims = (19.7, 19.99)
Slims = (35, 35.005)
νₑlims = (1e-6, 5e-3)
hm_w = heatmap!(ax_w, xw, zw, wₙ; colormap = :balance, colorrange = wlims)
Colorbar(fig[2, 2], hm_w; label = "m s⁻¹")
hm_T = heatmap!(ax_T, xT, zT, Tₙ; colormap = :thermal, colorrange = Tlims)
Colorbar(fig[2, 4], hm_T; label = "ᵒC")
hm_S = heatmap!(ax_S, xT, zT, Sₙ; colormap = :haline, colorrange = Slims)
Colorbar(fig[3, 2], hm_S; label = "g / kg")
hm_νₑ = heatmap!(ax_νₑ, xT, zT, νₑₙ; colormap = :thermal, colorrange = νₑlims)
Colorbar(fig[3, 4], hm_νₑ; label = "m s⁻²")
fig[1, 1:4] = Label(fig, title, fontsize=24, tellwidth=false)
fig
```

And now record a movie.

```
frames = intro:length(times)
@info "Making a motion picture of ocean wind mixing and convection..."
record(fig, filename * ".mp4", frames, framerate=8) do i
n[] = i
end
```

```
[ Info: Making a motion picture of ocean wind mixing and convection...
```

*This page was generated using Literate.jl.*