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 TEOS10EquationOfState.
• 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, SeawaterPolynomials, CUDA"

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

using Oceananigans
using Oceananigans.Units

using CUDA
using Random
using Printf
using CairoMakie
using SeawaterPolynomials.TEOS10: TEOS10EquationOfState

The grid

We use 128²×64 grid points with 1 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 = 128    # number of points in each of horizontal directions
Nz = 64          # number of points in the vertical direction

Lx = Ly = 128    # (m) domain horizontal extents
Lz = 64          # (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_interfaces(k) = Lz * (ζ₀(k) * Σ(k) - 1)

grid = RectilinearGrid(GPU(),
                       size = (Nx, Nx, Nz),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = z_interfaces)

128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CUDAGPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 128.0) regularly spaced with Δx=1.0
├── Periodic y ∈ [0.0, 128.0) regularly spaced with Δy=1.0
└── Bounded  z ∈ [-64.0, 0.0] variably spaced with min(Δz)=0.833413, max(Δz)=1.96618

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

fig = Figure(size=(1200, 800))
ax = Axis(fig[1, 1], ylabel = "z (m)", xlabel = "Vertical spacing (m)")

lines!(ax, zspacings(grid, Center()))
scatter!(ax, zspacings(grid, Center()))

fig

Buoyancy that depends on temperature and salinity

We use the SeawaterBuoyancy model with the TEOS10 equation of state,

ρₒ = 1026 # kg m⁻³, average density at the surface of the world ocean
equation_of_state = TEOS10EquationOfState(reference_density=ρₒ)
buoyancy = SeawaterBuoyancy(; equation_of_state)

SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: BoussinesqEquationOfState{Float64}

Boundary conditions

We calculate the surface temperature flux associated with surface cooling of 200 W m⁻², reference density ρₒ, and heat capacity cᴾ,

Q = 200   # W m⁻², surface _heat_ flux
cᴾ = 3991 # J K⁻¹ kg⁻¹, typical heat capacity for seawater

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

4.884283985946938e-5

Finally, we impose a temperature gradient dTdz both initially (see "Initial conditions" section below) and at the bottom of the domain, culminating in the boundary conditions on temperature,

dTdz = 0.01 # K m⁻¹

T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Jᵀ),
                                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ᴰ = 2e-3 # dimensionless drag coefficient
ρₐ = 1.2  # kg m⁻³, approximate average density of air at sea-level
τx = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²

-0.00023391812865497074

The boundary conditions on u are thus

u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τx))

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.000233918
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

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

@inline Jˢ(x, y, t, S, evaporation_rate) = - evaporation_rate * S # [salinity unit] m s⁻¹

where S is salinity. We use an evaporation 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 Jˢ, indicating that Jˢ depends on salinity S and passing the parameter evaporation_rate,

evaporation_bc = FluxBoundaryCondition(Jˢ, field_dependencies=:S, parameters=evaporation_rate)

FluxBoundaryCondition: ContinuousBoundaryFunction Jˢ 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 Jˢ at (Nothing, Nothing, Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

Model instantiation

We fill in the final details of the model here, i.e., 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 are not explicitly resolved.

model = NonhydrostaticModel(; grid, buoyancy,
                            tracers = (:T, :S),
                            coriolis = FPlane(f=1e-4),
                            closure = AnisotropicMinimumDissipation(),
                            boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs))

NonhydrostaticModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CUDAGPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: (T, S)
├── closure: AnisotropicMinimumDissipation{ExplicitTimeDiscretization, @NamedTuple{T::Float64, S::Float64}, Float64, Nothing}
├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
└── coriolis: FPlane{Float64}(f=0.0001)

Note: To use the Smagorinsky-Lilly turbulence closure (with a constant model coefficient) rather than AnisotropicMinimumDissipation, use closure = SmagorinskyLilly() in the model 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(τx)) * 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 2 hours, with adaptive time-stepping and progress printing.

simulation = Simulation(model, Δt=10, stop_time=2hours)

Simulation of NonhydrostaticModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 10 seconds
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 2 hours
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── 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

The TimeStepWizard helps ensure stable time-stepping with a Courant-Freidrichs-Lewy (CFL) number of 1.0.

wizard = TimeStepWizard(cfl=1, 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(40))

We then set up the simulation:

Output

We use the JLD2Writer to save $x, z$ slices of the velocity fields, tracer fields, and eddy diffusivities. The prefix keyword argument to JLD2Writer indicates that output will be saved in ocean_wind_mixing_and_convection.jld2.

# Create a NamedTuple with eddy viscosity
eddy_viscosity = (; νₑ = model.closure_fields.νₑ)

filename = "ocean_wind_mixing_and_convection"

simulation.output_writers[:slices] =
    JLD2Writer(model, merge(model.velocities, model.tracers, eddy_viscosity),
               filename = filename * ".jld2",
               indices = (:, grid.Ny/2, :),
               schedule = TimeInterval(1minute),
               overwrite_existing = true)

JLD2Writer 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]
├── file_splitting: NoFileSplitting
└── file size: 63.8 KiB

We're ready:

run!(simulation)

[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 11 seconds, max(|w|) = 1.3e-05 ms⁻¹, wall time: 0 seconds
[ Info:     ... simulation initialization complete (7.721 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (2.805 seconds).
Iteration: 0040, time: 7.344 minutes, Δt: 7.912 seconds, max(|w|) = 4.2e-05 ms⁻¹, wall time: 11.123 seconds
Iteration: 0080, time: 11.441 minutes, Δt: 4.724 seconds, max(|w|) = 8.8e-03 ms⁻¹, wall time: 11.642 seconds
Iteration: 0120, time: 14.202 minutes, Δt: 4.056 seconds, max(|w|) = 2.3e-02 ms⁻¹, wall time: 12.011 seconds
Iteration: 0160, time: 16.894 minutes, Δt: 4.082 seconds, max(|w|) = 2.5e-02 ms⁻¹, wall time: 12.429 seconds
Iteration: 0200, time: 19.542 minutes, Δt: 3.989 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 12.800 seconds
Iteration: 0240, time: 22.134 minutes, Δt: 3.900 seconds, max(|w|) = 2.5e-02 ms⁻¹, wall time: 13.209 seconds
Iteration: 0280, time: 24.603 minutes, Δt: 3.579 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 13.592 seconds
Iteration: 0320, time: 26.984 minutes, Δt: 3.558 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 13.971 seconds
Iteration: 0360, time: 29.171 minutes, Δt: 3.494 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 14.339 seconds
Iteration: 0400, time: 31.389 minutes, Δt: 3.212 seconds, max(|w|) = 2.6e-02 ms⁻¹, wall time: 14.711 seconds
Iteration: 0440, time: 33.489 minutes, Δt: 3.178 seconds, max(|w|) = 2.8e-02 ms⁻¹, wall time: 15.085 seconds
Iteration: 0480, time: 35.578 minutes, Δt: 3.211 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 15.478 seconds
Iteration: 0520, time: 37.690 minutes, Δt: 3.116 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 15.859 seconds
Iteration: 0560, time: 39.722 minutes, Δt: 3.110 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 16.223 seconds
Iteration: 0600, time: 41.693 minutes, Δt: 2.727 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 16.606 seconds
Iteration: 0640, time: 43.547 minutes, Δt: 3.047 seconds, max(|w|) = 3.8e-02 ms⁻¹, wall time: 16.988 seconds
Iteration: 0680, time: 45.560 minutes, Δt: 2.923 seconds, max(|w|) = 3.4e-02 ms⁻¹, wall time: 17.359 seconds
Iteration: 0720, time: 47.432 minutes, Δt: 2.858 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 17.737 seconds
Iteration: 0760, time: 49.285 minutes, Δt: 2.790 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 18.115 seconds
Iteration: 0800, time: 51.145 minutes, Δt: 2.909 seconds, max(|w|) = 3.7e-02 ms⁻¹, wall time: 18.485 seconds
Iteration: 0840, time: 53 minutes, Δt: 2.834 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 18.865 seconds
Iteration: 0880, time: 54.850 minutes, Δt: 2.835 seconds, max(|w|) = 4.2e-02 ms⁻¹, wall time: 19.249 seconds
Iteration: 0920, time: 56.635 minutes, Δt: 2.705 seconds, max(|w|) = 4.5e-02 ms⁻¹, wall time: 19.641 seconds
Iteration: 0960, time: 58.356 minutes, Δt: 2.657 seconds, max(|w|) = 4.1e-02 ms⁻¹, wall time: 20.025 seconds
Iteration: 1000, time: 1.001 hours, Δt: 2.633 seconds, max(|w|) = 4.1e-02 ms⁻¹, wall time: 20.391 seconds
Iteration: 1040, time: 1.031 hours, Δt: 2.638 seconds, max(|w|) = 4.2e-02 ms⁻¹, wall time: 20.764 seconds
Iteration: 1080, time: 1.060 hours, Δt: 2.571 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 21.148 seconds
Iteration: 1120, time: 1.088 hours, Δt: 2.629 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 21.540 seconds
Iteration: 1160, time: 1.117 hours, Δt: 2.655 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 21.912 seconds
Iteration: 1200, time: 1.145 hours, Δt: 2.711 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 22.283 seconds
Iteration: 1240, time: 1.175 hours, Δt: 2.708 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 22.661 seconds
Iteration: 1280, time: 1.203 hours, Δt: 2.648 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 23.047 seconds
Iteration: 1320, time: 1.232 hours, Δt: 2.524 seconds, max(|w|) = 4.4e-02 ms⁻¹, wall time: 23.430 seconds
Iteration: 1360, time: 1.259 hours, Δt: 2.693 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 23.814 seconds
Iteration: 1400, time: 1.288 hours, Δt: 2.495 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 24.219 seconds
Iteration: 1440, time: 1.315 hours, Δt: 2.545 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 24.577 seconds
Iteration: 1480, time: 1.343 hours, Δt: 2.475 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 24.961 seconds
Iteration: 1520, time: 1.370 hours, Δt: 2.491 seconds, max(|w|) = 6.4e-02 ms⁻¹, wall time: 25.347 seconds
Iteration: 1560, time: 1.398 hours, Δt: 2.490 seconds, max(|w|) = 6.6e-02 ms⁻¹, wall time: 25.729 seconds
Iteration: 1600, time: 1.425 hours, Δt: 2.536 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 26.119 seconds
Iteration: 1640, time: 1.453 hours, Δt: 2.432 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 26.503 seconds
Iteration: 1680, time: 1.480 hours, Δt: 2.501 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 26.890 seconds
Iteration: 1720, time: 1.508 hours, Δt: 2.594 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 27.289 seconds
Iteration: 1760, time: 1.535 hours, Δt: 2.424 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 27.729 seconds
Iteration: 1800, time: 1.562 hours, Δt: 2.487 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 28.237 seconds
Iteration: 1840, time: 1.589 hours, Δt: 2.482 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 28.679 seconds
Iteration: 1880, time: 1.616 hours, Δt: 2.520 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 29.149 seconds
Iteration: 1920, time: 1.644 hours, Δt: 2.381 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 29.577 seconds
Iteration: 1960, time: 1.670 hours, Δt: 2.344 seconds, max(|w|) = 6.0e-02 ms⁻¹, wall time: 29.968 seconds
Iteration: 2000, time: 1.696 hours, Δt: 2.443 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 30.356 seconds
Iteration: 2040, time: 1.723 hours, Δt: 2.442 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 30.751 seconds
Iteration: 2080, time: 1.750 hours, Δt: 2.444 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 31.144 seconds
Iteration: 2120, time: 1.776 hours, Δt: 2.434 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 31.555 seconds
Iteration: 2160, time: 1.803 hours, Δt: 2.470 seconds, max(|w|) = 5.9e-02 ms⁻¹, wall time: 31.999 seconds
Iteration: 2200, time: 1.829 hours, Δt: 2.438 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 32.400 seconds
Iteration: 2240, time: 1.856 hours, Δt: 2.436 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 32.775 seconds
Iteration: 2280, time: 1.883 hours, Δt: 2.417 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 33.159 seconds
Iteration: 2320, time: 1.909 hours, Δt: 2.405 seconds, max(|w|) = 5.5e-02 ms⁻¹, wall time: 33.587 seconds
Iteration: 2360, time: 1.935 hours, Δt: 2.438 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 33.988 seconds
Iteration: 2400, time: 1.962 hours, Δt: 2.334 seconds, max(|w|) = 5.9e-02 ms⁻¹, wall time: 34.394 seconds
Iteration: 2440, time: 1.987 hours, Δt: 2.337 seconds, max(|w|) = 6.8e-02 ms⁻¹, wall time: 34.796 seconds
[ Info: Simulation is stopping after running for 34.988 seconds.
[ Info: Simulation time 2 hours equals or exceeds stop time 2 hours.

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

(w = 128×1×65×121 FieldTimeSeries{InMemory} located at (Center, Center, Face) of w at ocean_wind_mixing_and_convection.jld2
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: w
└── data: 134×1×71×121 OffsetArray(::Array{Float64, 4}, -2:131, 64:64, -2:68, 1:121) with eltype Float64 with indices -2:131×64:64×-2:68×1:121
    └── max=0.0381107, min=-0.0584959, mean=2.63258e-5, T = 128×1×64×121 FieldTimeSeries{InMemory} located at (Center, Center, Center) of T at ocean_wind_mixing_and_convection.jld2
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: T
└── data: 134×1×70×121 OffsetArray(::Array{Float64, 4}, -2:131, 64:64, -2:67, 1:121) with eltype Float64 with indices -2:131×64:64×-2:67×1:121
    └── max=19.9958, min=0.0, mean=18.5885, S = 128×1×64×121 FieldTimeSeries{InMemory} located at (Center, Center, Center) of S at ocean_wind_mixing_and_convection.jld2
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: S
└── data: 134×1×70×121 OffsetArray(::Array{Float64, 4}, -2:131, 64:64, -2:67, 1:121) with eltype Float64 with indices -2:131×64:64×-2:67×1:121
    └── max=35.0137, min=0.0, mean=33.0007, νₑ = 128×1×64×121 FieldTimeSeries{InMemory} located at (Center, Center, Center) of νₑ at ocean_wind_mixing_and_convection.jld2
├── grid: 128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: νₑ
└── data: 134×1×70×121 OffsetArray(::Array{Float64, 4}, -2:131, 64:64, -2:67, 1:121) with eltype Float64 with indices -2:131×64:64×-2:67×1:121
    └── max=0.0191491, min=0.0, mean=0.000425782)

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 Observables work we refer to Makie.jl's Documentation.

n = Observable(intro)

 wₙ = @lift time_series.w[$n]
 Tₙ = @lift time_series.T[$n]
 Sₙ = @lift time_series.S[$n]
νₑₙ = @lift time_series.νₑ[$n]

fig = Figure(size = (1800, 900))

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, wₙ; colormap = :balance, colorrange = wlims)
Colorbar(fig[2, 2], hm_w; label = "m s⁻¹")

hm_T = heatmap!(ax_T, Tₙ; colormap = :thermal, colorrange = Tlims)
Colorbar(fig[2, 4], hm_T; label = "ᵒC")

hm_S = heatmap!(ax_S, Sₙ; colormap = :haline, colorrange = Slims)
Colorbar(fig[3, 2], hm_S; label = "g / kg")

hm_νₑ = heatmap!(ax_νₑ, νₑₙ; 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..."

CairoMakie.record(fig, filename * ".mp4", frames, framerate=8) do i
    n[] = i
end

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

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.12.2
Commit ca9b6662be4 (2025-11-20 16:25 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 128 × AMD EPYC 9374F 32-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
Environment:
  LD_LIBRARY_PATH = 
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia-oceananigans
  JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-27500/docs/
  JULIA_VERSION = 1.12.2
  JULIA_LOAD_PATH = @:@v#.#:@stdlib
  JULIA_VERSION_ENZYME = 1.10.10
  JULIA_PYTHONCALL_EXE = /var/lib/buildkite-agent/Oceananigans.jl-27500/docs/.CondaPkg/.pixi/envs/default/bin/python
  JULIA_DEBUG = Literate

These were the top-level packages installed in the environment:

import Pkg
Pkg.status()

Status `~/Oceananigans.jl-27500/docs/Project.toml`
  [79e6a3ab] Adapt v4.4.0
  [052768ef] CUDA v5.9.5
  [13f3f980] CairoMakie v0.15.8
  [e30172f5] Documenter v1.16.1
  [daee34ce] DocumenterCitations v1.4.1
  [033835bb] JLD2 v0.6.3
  [98b081ad] Literate v2.21.0
  [da04e1cc] MPI v0.20.23
  [85f8d34a] NCDatasets v0.14.10
  [9e8cae18] Oceananigans v0.102.4 `~/Oceananigans.jl-27500`
  [f27b6e38] Polynomials v4.1.0
  [6038ab10] Rotations v1.7.1
  [d496a93d] SeawaterPolynomials v0.3.10
  [09ab397b] StructArrays v0.7.2
  [bdfc003b] TimesDates v0.3.3
  [2e0b0046] XESMF v0.1.6
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.12.0
  [44cfe95a] Pkg v1.12.0

This page was generated using Literate.jl.