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 (8.579 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (3.651 seconds).
Iteration: 0040, time: 7.344 minutes, Δt: 7.911 seconds, max(|w|) = 4.3e-05 ms⁻¹, wall time: 12.947 seconds
Iteration: 0080, time: 11.439 minutes, Δt: 4.718 seconds, max(|w|) = 9.3e-03 ms⁻¹, wall time: 13.407 seconds
Iteration: 0120, time: 14.209 minutes, Δt: 4.133 seconds, max(|w|) = 2.1e-02 ms⁻¹, wall time: 14.461 seconds
Iteration: 0160, time: 16.924 minutes, Δt: 4.097 seconds, max(|w|) = 2.7e-02 ms⁻¹, wall time: 15.721 seconds
Iteration: 0200, time: 19.524 minutes, Δt: 3.987 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 16.629 seconds
Iteration: 0240, time: 22 minutes, Δt: 3.912 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 16.995 seconds
Iteration: 0280, time: 24.416 minutes, Δt: 3.630 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 17.415 seconds
Iteration: 0320, time: 26.805 minutes, Δt: 3.725 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 17.821 seconds
Iteration: 0360, time: 29.059 minutes, Δt: 3.588 seconds, max(|w|) = 3.5e-02 ms⁻¹, wall time: 18.188 seconds
Iteration: 0400, time: 31.339 minutes, Δt: 3.402 seconds, max(|w|) = 3.4e-02 ms⁻¹, wall time: 18.608 seconds
Iteration: 0440, time: 33.567 minutes, Δt: 3.450 seconds, max(|w|) = 3.4e-02 ms⁻¹, wall time: 19.062 seconds
Iteration: 0480, time: 35.753 minutes, Δt: 3.168 seconds, max(|w|) = 2.9e-02 ms⁻¹, wall time: 19.428 seconds
Iteration: 0520, time: 37.851 minutes, Δt: 3.208 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 19.831 seconds
Iteration: 0560, time: 39.953 minutes, Δt: 3.204 seconds, max(|w|) = 3.2e-02 ms⁻¹, wall time: 20.224 seconds
Iteration: 0600, time: 42 minutes, Δt: 3.229 seconds, max(|w|) = 3.0e-02 ms⁻¹, wall time: 20.668 seconds
Iteration: 0640, time: 44.050 minutes, Δt: 2.880 seconds, max(|w|) = 3.1e-02 ms⁻¹, wall time: 21.065 seconds
Iteration: 0680, time: 45.962 minutes, Δt: 2.953 seconds, max(|w|) = 3.2e-02 ms⁻¹, wall time: 21.443 seconds
Iteration: 0720, time: 47.840 minutes, Δt: 2.912 seconds, max(|w|) = 3.3e-02 ms⁻¹, wall time: 21.808 seconds
Iteration: 0760, time: 49.714 minutes, Δt: 2.893 seconds, max(|w|) = 4.2e-02 ms⁻¹, wall time: 22.210 seconds
Iteration: 0800, time: 51.641 minutes, Δt: 2.979 seconds, max(|w|) = 4.6e-02 ms⁻¹, wall time: 22.591 seconds
Iteration: 0840, time: 53.517 minutes, Δt: 2.755 seconds, max(|w|) = 4.6e-02 ms⁻¹, wall time: 22.969 seconds
Iteration: 0880, time: 55.316 minutes, Δt: 2.525 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 23.360 seconds
Iteration: 0920, time: 56.993 minutes, Δt: 2.591 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 23.809 seconds
Iteration: 0960, time: 58.714 minutes, Δt: 2.603 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 24.210 seconds
Iteration: 1000, time: 1.007 hours, Δt: 2.604 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 24.616 seconds
Iteration: 1040, time: 1.035 hours, Δt: 2.770 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 25.034 seconds
Iteration: 1080, time: 1.066 hours, Δt: 2.776 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 25.426 seconds
Iteration: 1120, time: 1.095 hours, Δt: 2.738 seconds, max(|w|) = 4.7e-02 ms⁻¹, wall time: 25.855 seconds
Iteration: 1160, time: 1.126 hours, Δt: 2.821 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 26.319 seconds
Iteration: 1200, time: 1.156 hours, Δt: 2.643 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 26.732 seconds
Iteration: 1240, time: 1.184 hours, Δt: 2.701 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 27.102 seconds
Iteration: 1280, time: 1.214 hours, Δt: 2.801 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 27.528 seconds
Iteration: 1320, time: 1.244 hours, Δt: 2.737 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 27.987 seconds
Iteration: 1360, time: 1.273 hours, Δt: 2.611 seconds, max(|w|) = 5.0e-02 ms⁻¹, wall time: 28.408 seconds
Iteration: 1400, time: 1.301 hours, Δt: 2.443 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 28.826 seconds
Iteration: 1440, time: 1.328 hours, Δt: 2.463 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 29.200 seconds
Iteration: 1480, time: 1.356 hours, Δt: 2.599 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 29.584 seconds
Iteration: 1520, time: 1.384 hours, Δt: 2.504 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 29.963 seconds
Iteration: 1560, time: 1.412 hours, Δt: 2.653 seconds, max(|w|) = 5.7e-02 ms⁻¹, wall time: 30.380 seconds
Iteration: 1600, time: 1.440 hours, Δt: 2.513 seconds, max(|w|) = 4.2e-02 ms⁻¹, wall time: 30.879 seconds
Iteration: 1640, time: 1.467 hours, Δt: 2.487 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 31.269 seconds
Iteration: 1680, time: 1.494 hours, Δt: 2.473 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 31.691 seconds
Iteration: 1720, time: 1.521 hours, Δt: 2.496 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 32.167 seconds
Iteration: 1760, time: 1.548 hours, Δt: 2.267 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 32.544 seconds
Iteration: 1800, time: 1.572 hours, Δt: 2.317 seconds, max(|w|) = 5.8e-02 ms⁻¹, wall time: 32.958 seconds
Iteration: 1840, time: 1.598 hours, Δt: 2.242 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 33.386 seconds
Iteration: 1880, time: 1.622 hours, Δt: 2.214 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 33.862 seconds
Iteration: 1920, time: 1.646 hours, Δt: 2.326 seconds, max(|w|) = 5.1e-02 ms⁻¹, wall time: 34.389 seconds
Iteration: 1960, time: 1.670 hours, Δt: 2.211 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 34.917 seconds
Iteration: 2000, time: 1.695 hours, Δt: 2.279 seconds, max(|w|) = 4.3e-02 ms⁻¹, wall time: 35.501 seconds
Iteration: 2040, time: 1.719 hours, Δt: 2.259 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 35.900 seconds
Iteration: 2080, time: 1.744 hours, Δt: 2.268 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 36.338 seconds
Iteration: 2120, time: 1.768 hours, Δt: 2.272 seconds, max(|w|) = 4.8e-02 ms⁻¹, wall time: 36.756 seconds
Iteration: 2160, time: 1.793 hours, Δt: 2.261 seconds, max(|w|) = 4.9e-02 ms⁻¹, wall time: 37.134 seconds
Iteration: 2200, time: 1.817 hours, Δt: 2.299 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 37.546 seconds
Iteration: 2240, time: 1.843 hours, Δt: 2.326 seconds, max(|w|) = 5.2e-02 ms⁻¹, wall time: 38.078 seconds
Iteration: 2280, time: 1.868 hours, Δt: 2.233 seconds, max(|w|) = 5.6e-02 ms⁻¹, wall time: 38.481 seconds
Iteration: 2320, time: 1.893 hours, Δt: 2.329 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 39.058 seconds
Iteration: 2360, time: 1.918 hours, Δt: 2.164 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 39.488 seconds
Iteration: 2400, time: 1.942 hours, Δt: 2.198 seconds, max(|w|) = 5.3e-02 ms⁻¹, wall time: 39.915 seconds
Iteration: 2440, time: 1.967 hours, Δt: 2.346 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 40.323 seconds
Iteration: 2480, time: 1.992 hours, Δt: 2.254 seconds, max(|w|) = 5.4e-02 ms⁻¹, wall time: 40.744 seconds
[ Info: Simulation is stopping after running for 40.879 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.0437704, min=-0.0431349, mean=2.36976e-6, 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.5884, 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.0136, 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.0197221, min=0.0, mean=0.000397864)

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

This page was generated using Literate.jl.