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

```
Nz = 24 # number of points in the vertical direction
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(CPU();
size = (32, 32, Nz),
x = (0, 64),
y = (0, 64),
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(resolution=(1200, 800))
ax = Axis(fig[1, 1], ylabel = "Depth (m)", xlabel = "Vertical spacing (m)")
lines!(ax, grid.Δzᵃᵃᶜ[1:grid.Nz], grid.zᵃᵃᶜ[1:grid.Nz])
scatter!(ax, grid.Δzᵃᵃᶜ[1:Nz], grid.zᵃᵃᶜ[1:Nz])
save("ocean_wind_mixing_convection_grid_spacing.svg", 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 heating 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 -ĝ = ZDirection
└── 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 years
├── 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))
simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(20))
```

`Callback of progress_message on 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{Float32}
├── 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|) = 9.9e-06 ms⁻¹, wall time: 0 seconds
[ Info: ... simulation initialization complete (7.177 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (2.106 minutes).
Iteration: 0020, time: 3.403 minutes, Δt: 13.310 seconds, max(|w|) = 8.2e-06 ms⁻¹, wall time: 2.271 minutes
Iteration: 0040, time: 7.488 minutes, Δt: 16.105 seconds, max(|w|) = 6.0e-06 ms⁻¹, wall time: 2.314 minutes
Iteration: 0060, time: 12 minutes, Δt: 10.350 seconds, max(|w|) = 1.2e-05 ms⁻¹, wall time: 2.354 minutes
Iteration: 0080, time: 14.919 minutes, Δt: 8.329 seconds, max(|w|) = 4.0e-05 ms⁻¹, wall time: 2.399 minutes
Iteration: 0100, time: 17.385 minutes, Δt: 7.150 seconds, max(|w|) = 1.4e-04 ms⁻¹, wall time: 2.446 minutes
Iteration: 0120, time: 19.561 minutes, Δt: 6.355 seconds, max(|w|) = 7.1e-04 ms⁻¹, wall time: 2.488 minutes
Iteration: 0140, time: 21.505 minutes, Δt: 5.782 seconds, max(|w|) = 3.6e-03 ms⁻¹, wall time: 2.534 minutes
Iteration: 0160, time: 23.277 minutes, Δt: 5.351 seconds, max(|w|) = 1.2e-02 ms⁻¹, wall time: 2.574 minutes
Iteration: 0180, time: 24.976 minutes, Δt: 5.428 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 2.616 minutes
Iteration: 0200, time: 26.766 minutes, Δt: 6.070 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 2.657 minutes
Iteration: 0220, time: 28.727 minutes, Δt: 6.852 seconds, max(|w|) = 8.3e-02 ms⁻¹, wall time: 2.698 minutes
Iteration: 0240, time: 30.931 minutes, Δt: 7.079 seconds, max(|w|) = 7.8e-02 ms⁻¹, wall time: 2.741 minutes
Iteration: 0260, time: 33.118 minutes, Δt: 7.153 seconds, max(|w|) = 7.2e-02 ms⁻¹, wall time: 2.787 minutes
Iteration: 0280, time: 35.525 minutes, Δt: 8.287 seconds, max(|w|) = 7.7e-02 ms⁻¹, wall time: 2.829 minutes
Iteration: 0300, time: 38 minutes, Δt: 9.409 seconds, max(|w|) = 6.7e-02 ms⁻¹, wall time: 2.873 minutes
[ 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 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 … 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 … 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 … -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 = 10minutes$ 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(resolution = (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, textsize=24, tellwidth=false)
```

`Label()`

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
msg = string("Plotting frame ", i, " of ", frames[end])
print(msg * " \r")
n[] = i
end
```

```
[ Info: Making a motion picture of ocean wind mixing and convection...
Plotting frame 11 of 41
Plotting frame 12 of 41
Plotting frame 13 of 41
Plotting frame 14 of 41
Plotting frame 15 of 41
Plotting frame 16 of 41
Plotting frame 17 of 41
Plotting frame 18 of 41
Plotting frame 19 of 41
Plotting frame 20 of 41
Plotting frame 21 of 41
Plotting frame 22 of 41
Plotting frame 23 of 41
Plotting frame 24 of 41
Plotting frame 25 of 41
Plotting frame 26 of 41
Plotting frame 27 of 41
Plotting frame 28 of 41
Plotting frame 29 of 41
Plotting frame 30 of 41
Plotting frame 31 of 41
Plotting frame 32 of 41
Plotting frame 33 of 41
Plotting frame 34 of 41
Plotting frame 35 of 41
Plotting frame 36 of 41
Plotting frame 37 of 41
Plotting frame 38 of 41
Plotting frame 39 of 41
Plotting frame 40 of 41
Plotting frame 41 of 41
```

*This page was generated using Literate.jl.*