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 use the (scale-invariant) DynamicSmagorinsky 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 = DynamicSmagorinsky(),
                            boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs))

NonhydrostaticModel{CUDAGPU, 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: Smagorinsky with coefficient = DynamicCoefficient(averaging = LagrangianAveraging(), schedule = IterationInterval(1, 0)), Pr=(T = 1.0, S = 1.0)
├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
└── coriolis: FPlane{Float64}(f=0.0001)

Note: To use the (constant-coefficient) Smagorinsky-Lilly turbulence closure rather than DynamicSmagorinsky, 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{CUDAGPU, 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: 0 bytes (file not yet created)

We're ready:

run!(simulation)

[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 11 seconds, max(|w|) = 1.5e-05 ms⁻¹, wall time: 0 seconds
[ Info:     ... simulation initialization complete (9.939 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.971 seconds).
Iteration: 0040, time: 7.344 minutes, Δt: 7.914 seconds, max(|w|) = 7.2e-06 ms⁻¹, wall time: 23.252 seconds
Iteration: 0080, time: 11.453 minutes, Δt: 5.016 seconds, max(|w|) = 1.0e-05 ms⁻¹, wall time: 31.432 seconds
Iteration: 0120, time: 14.416 minutes, Δt: 3.953 seconds, max(|w|) = 8.7e-06 ms⁻¹, wall time: 40.556 seconds
Iteration: 0160, time: 16.823 minutes, Δt: 3.365 seconds, max(|w|) = 9.1e-06 ms⁻¹, wall time: 57.543 seconds
Iteration: 0200, time: 18.880 minutes, Δt: 2.981 seconds, max(|w|) = 2.8e-05 ms⁻¹, wall time: 1.267 minutes
Iteration: 0240, time: 20.744 minutes, Δt: 2.698 seconds, max(|w|) = 9.0e-05 ms⁻¹, wall time: 1.524 minutes
Iteration: 0280, time: 22.464 minutes, Δt: 2.473 seconds, max(|w|) = 2.5e-04 ms⁻¹, wall time: 1.819 minutes
Iteration: 0320, time: 24 minutes, Δt: 2.286 seconds, max(|w|) = 7.5e-04 ms⁻¹, wall time: 2.129 minutes
Iteration: 0360, time: 25.467 minutes, Δt: 2.097 seconds, max(|w|) = 2.3e-03 ms⁻¹, wall time: 2.355 minutes
Iteration: 0400, time: 26.783 minutes, Δt: 1.845 seconds, max(|w|) = 5.8e-03 ms⁻¹, wall time: 2.500 minutes
Iteration: 0440, time: 27.901 minutes, Δt: 1.552 seconds, max(|w|) = 2.0e-02 ms⁻¹, wall time: 2.637 minutes
Iteration: 0480, time: 28.860 minutes, Δt: 1.195 seconds, max(|w|) = 1.1e-01 ms⁻¹, wall time: 2.786 minutes
Iteration: 0520, time: 29.556 minutes, Δt: 1.051 seconds, max(|w|) = 2.3e-01 ms⁻¹, wall time: 2.949 minutes
Iteration: 0560, time: 30.206 minutes, Δt: 981.532 ms, max(|w|) = 2.9e-01 ms⁻¹, wall time: 3.250 minutes
Iteration: 0600, time: 30.838 minutes, Δt: 1.017 seconds, max(|w|) = 2.4e-01 ms⁻¹, wall time: 3.551 minutes
Iteration: 0640, time: 31.516 minutes, Δt: 1.194 seconds, max(|w|) = 2.2e-01 ms⁻¹, wall time: 3.891 minutes
Iteration: 0680, time: 32.429 minutes, Δt: 1.716 seconds, max(|w|) = 1.9e-01 ms⁻¹, wall time: 4.197 minutes
Iteration: 0720, time: 33.630 minutes, Δt: 2.015 seconds, max(|w|) = 1.5e-01 ms⁻¹, wall time: 4.493 minutes
Iteration: 0760, time: 35.078 minutes, Δt: 2.411 seconds, max(|w|) = 1.2e-01 ms⁻¹, wall time: 4.736 minutes
Iteration: 0800, time: 36.740 minutes, Δt: 2.728 seconds, max(|w|) = 1.2e-01 ms⁻¹, wall time: 4.886 minutes
Iteration: 0840, time: 38.542 minutes, Δt: 2.711 seconds, max(|w|) = 9.3e-02 ms⁻¹, wall time: 5.022 minutes
Iteration: 0880, time: 40.355 minutes, Δt: 2.930 seconds, max(|w|) = 9.8e-02 ms⁻¹, wall time: 5.157 minutes
Iteration: 0920, time: 42.259 minutes, Δt: 2.984 seconds, max(|w|) = 7.3e-02 ms⁻¹, wall time: 5.293 minutes
Iteration: 0960, time: 44.260 minutes, Δt: 3.069 seconds, max(|w|) = 7.1e-02 ms⁻¹, wall time: 5.430 minutes
Iteration: 1000, time: 46.199 minutes, Δt: 3.085 seconds, max(|w|) = 6.7e-02 ms⁻¹, wall time: 5.564 minutes
Iteration: 1040, time: 48.144 minutes, Δt: 2.977 seconds, max(|w|) = 6.8e-02 ms⁻¹, wall time: 5.701 minutes
Iteration: 1080, time: 50.047 minutes, Δt: 2.871 seconds, max(|w|) = 7.9e-02 ms⁻¹, wall time: 5.835 minutes
Iteration: 1120, time: 51.980 minutes, Δt: 2.759 seconds, max(|w|) = 6.6e-02 ms⁻¹, wall time: 6.119 minutes
Iteration: 1160, time: 53.831 minutes, Δt: 2.817 seconds, max(|w|) = 6.7e-02 ms⁻¹, wall time: 6.421 minutes
Iteration: 1200, time: 55.636 minutes, Δt: 2.770 seconds, max(|w|) = 7.0e-02 ms⁻¹, wall time: 6.720 minutes
Iteration: 1240, time: 57.459 minutes, Δt: 2.695 seconds, max(|w|) = 6.6e-02 ms⁻¹, wall time: 7.021 minutes
Iteration: 1280, time: 59.285 minutes, Δt: 2.700 seconds, max(|w|) = 6.8e-02 ms⁻¹, wall time: 7.167 minutes
Iteration: 1320, time: 1.018 hours, Δt: 2.653 seconds, max(|w|) = 7.1e-02 ms⁻¹, wall time: 7.307 minutes
Iteration: 1360, time: 1.048 hours, Δt: 2.639 seconds, max(|w|) = 7.3e-02 ms⁻¹, wall time: 7.447 minutes
Iteration: 1400, time: 1.077 hours, Δt: 2.746 seconds, max(|w|) = 6.6e-02 ms⁻¹, wall time: 7.586 minutes
Iteration: 1440, time: 1.107 hours, Δt: 2.738 seconds, max(|w|) = 7.4e-02 ms⁻¹, wall time: 7.728 minutes
Iteration: 1480, time: 1.135 hours, Δt: 2.666 seconds, max(|w|) = 7.0e-02 ms⁻¹, wall time: 7.868 minutes
Iteration: 1520, time: 1.163 hours, Δt: 2.510 seconds, max(|w|) = 6.7e-02 ms⁻¹, wall time: 8.062 minutes
Iteration: 1560, time: 1.191 hours, Δt: 2.575 seconds, max(|w|) = 7.1e-02 ms⁻¹, wall time: 8.420 minutes
Iteration: 1600, time: 1.219 hours, Δt: 2.607 seconds, max(|w|) = 7.8e-02 ms⁻¹, wall time: 8.830 minutes
Iteration: 1640, time: 1.246 hours, Δt: 2.661 seconds, max(|w|) = 8.2e-02 ms⁻¹, wall time: 9.217 minutes
Iteration: 1680, time: 1.274 hours, Δt: 2.460 seconds, max(|w|) = 6.8e-02 ms⁻¹, wall time: 9.408 minutes
Iteration: 1720, time: 1.301 hours, Δt: 2.324 seconds, max(|w|) = 7.3e-02 ms⁻¹, wall time: 9.574 minutes
Iteration: 1760, time: 1.328 hours, Δt: 2.453 seconds, max(|w|) = 8.5e-02 ms⁻¹, wall time: 9.725 minutes
Iteration: 1800, time: 1.355 hours, Δt: 2.460 seconds, max(|w|) = 8.0e-02 ms⁻¹, wall time: 9.870 minutes
Iteration: 1840, time: 1.382 hours, Δt: 2.570 seconds, max(|w|) = 8.3e-02 ms⁻¹, wall time: 10.010 minutes
Iteration: 1880, time: 1.408 hours, Δt: 2.394 seconds, max(|w|) = 7.7e-02 ms⁻¹, wall time: 10.151 minutes
Iteration: 1920, time: 1.435 hours, Δt: 2.465 seconds, max(|w|) = 8.0e-02 ms⁻¹, wall time: 10.292 minutes
Iteration: 1960, time: 1.462 hours, Δt: 2.406 seconds, max(|w|) = 7.5e-02 ms⁻¹, wall time: 10.432 minutes
Iteration: 2000, time: 1.487 hours, Δt: 2.596 seconds, max(|w|) = 7.8e-02 ms⁻¹, wall time: 10.576 minutes
Iteration: 2040, time: 1.514 hours, Δt: 2.431 seconds, max(|w|) = 7.8e-02 ms⁻¹, wall time: 10.720 minutes
Iteration: 2080, time: 1.541 hours, Δt: 2.499 seconds, max(|w|) = 9.0e-02 ms⁻¹, wall time: 10.887 minutes
Iteration: 2120, time: 1.567 hours, Δt: 2.326 seconds, max(|w|) = 8.2e-02 ms⁻¹, wall time: 11.058 minutes
Iteration: 2160, time: 1.593 hours, Δt: 2.266 seconds, max(|w|) = 7.7e-02 ms⁻¹, wall time: 11.199 minutes
Iteration: 2200, time: 1.619 hours, Δt: 2.465 seconds, max(|w|) = 8.2e-02 ms⁻¹, wall time: 11.342 minutes
Iteration: 2240, time: 1.645 hours, Δt: 2.498 seconds, max(|w|) = 7.9e-02 ms⁻¹, wall time: 11.485 minutes
Iteration: 2280, time: 1.672 hours, Δt: 2.459 seconds, max(|w|) = 7.6e-02 ms⁻¹, wall time: 11.629 minutes
Iteration: 2320, time: 1.699 hours, Δt: 2.328 seconds, max(|w|) = 8.0e-02 ms⁻¹, wall time: 11.772 minutes
Iteration: 2360, time: 1.724 hours, Δt: 2.290 seconds, max(|w|) = 8.2e-02 ms⁻¹, wall time: 11.917 minutes
Iteration: 2400, time: 1.749 hours, Δt: 2.329 seconds, max(|w|) = 8.5e-02 ms⁻¹, wall time: 12.062 minutes
Iteration: 2440, time: 1.775 hours, Δt: 2.404 seconds, max(|w|) = 7.8e-02 ms⁻¹, wall time: 12.209 minutes
Iteration: 2480, time: 1.801 hours, Δt: 2.403 seconds, max(|w|) = 7.5e-02 ms⁻¹, wall time: 12.353 minutes
Iteration: 2520, time: 1.827 hours, Δt: 2.372 seconds, max(|w|) = 8.2e-02 ms⁻¹, wall time: 12.497 minutes
Iteration: 2560, time: 1.853 hours, Δt: 2.391 seconds, max(|w|) = 7.1e-02 ms⁻¹, wall time: 12.647 minutes
Iteration: 2600, time: 1.878 hours, Δt: 2.354 seconds, max(|w|) = 7.5e-02 ms⁻¹, wall time: 12.786 minutes
Iteration: 2640, time: 1.903 hours, Δt: 2.202 seconds, max(|w|) = 7.7e-02 ms⁻¹, wall time: 12.933 minutes
Iteration: 2680, time: 1.928 hours, Δt: 2.253 seconds, max(|w|) = 8.2e-02 ms⁻¹, wall time: 13.078 minutes
Iteration: 2720, time: 1.953 hours, Δt: 2.197 seconds, max(|w|) = 7.8e-02 ms⁻¹, wall time: 13.222 minutes
Iteration: 2760, time: 1.978 hours, Δt: 2.325 seconds, max(|w|) = 8.5e-02 ms⁻¹, wall time: 13.366 minutes
[ Info: Simulation is stopping after running for 13.509 minutes.
[ 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.0985793, min=-0.17129, mean=-1.46908e-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=20.0173, 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.0278, 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.0125521, min=0.0, mean=9.97815e-5)

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.3
Commit 966d0af0fdf (2025-12-15 11:20 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-28307/docs/
  JULIA_VERSION = 1.12.3
  JULIA_LOAD_PATH = @:@v#.#:@stdlib
  JULIA_VERSION_ENZYME = 1.10.10
  JULIA_PYTHONCALL_EXE = /var/lib/buildkite-agent/Oceananigans.jl-28307/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-28307/docs/Project.toml`
  [79e6a3ab] Adapt v4.4.0
  [052768ef] CUDA v5.9.6
  [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.104.0 `~/Oceananigans.jl-28307`
  [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.1

This page was generated using Literate.jl.