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.

julia

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.

julia

using Oceananigans
using Oceananigans.Units

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

Random.seed!(1969) # for reproducible results

Random.TaskLocalRNG()

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:

julia

Nx = Ny = 128
Nz = 64

Lx = Ly = 128meters
Lz = 64meters

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, Ny, Nz),
                       halo = (4, 4, 4),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = z_interfaces)

128×128×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CUDAGPU with 4×4×4 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:

julia

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,

julia

ρₒ = 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ᴾ,

julia

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 the "Initial conditions" section below) and at the bottom of the domain, culminating in the boundary conditions on temperature,

julia

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:

julia

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

julia

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

julia

@inline Jˢ(x, y, t, S, evaporation_rate) = - evaporation_rate * S

We use an evaporation rate of 1 millimeter per hour,

julia

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,

julia

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

FluxBoundaryCondition: ContinuousBoundaryFunction Jˢ at (Nothing, Nothing, Nothing)

The full salinity boundary conditions are

julia

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, advection scheme, and use the DynamicSmagorinsky closure for large eddy simulation. The effect of both the WENO advection scheme and DynamicSmagorinsky turbulence closure is to dissipate variance at the grid scale. In the context of large eddy simulation, this dissipation may be interpreted as approximating a forward cascade of kinetic energy from resolved motions into motions that are smaller than the grid scale and not explicitly resolved. Note that dissipation of grid-scale variance can be achieved by the WENO advection scheme alone with closure = nothing. Typically, using an explicit closure = DynamicSmagorinsky() will produce stronger dissipation of kinetic energy; whether or not this leads to a higher quality numerical solution depends on the context. An explicit closure = DynamicSmagorinsky() is also useful for diagnosing the kinetic energy dissipation rate (the dissipation rate associated with WENO advection can be computed in principle, but is challenging and relatively computationally intensive).

In the DynamicSmagorinsky closure used below, a dynamic, multi-scale method is used to estimate the Smagorinsky coefficient at every point in time and space. Specifically, we use an algorithm that assumes the coefficient does not depend on the spatial scale of the implicit filter that separates the "true", underlying, and unresolved flow from the "filtered", or computed flow. This implementation corresponds to the "scale-invariant" formulation described by .

julia

model = NonhydrostaticModel(grid; buoyancy,
                            advection = WENO(order=7),
                            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 4×4×4 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme:
│   ├── momentum: WENO{4, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=7)
│   ├── T: WENO{4, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=7)
│   └── S: WENO{4, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=7)
├── 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 ĝ = 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 top and bottom, while our initial condition for velocity consists only of random noise scaled by the friction velocity .

julia

Ξ(z) = randn() * z / Lz * (1 + z / Lz)

Tᵢ(x, y, z) = 20 + dTdz * z + dTdz * Lz * 2e-6 * Ξ(z)
uᵢ(x, y, z) = sqrt(abs(τx)) * 1e-3 * Ξ(z)

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.

julia

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-Friedrichs-Lewy (CFL) number of 0.7.

julia

conjure_time_step_wizard!(simulation, cfl=0.7)

Nice progress messaging is helpful:

julia

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(200))

Output

We use the JLD2Writer to save slices of the velocity fields, tracer fields, and eddy viscosity in ocean_wind_mixing_and_convection.jld2.

julia

eddy_viscosity = (; νₑ = model.closure_fields.νₑ)

filename = "ocean_wind_mixing_and_convection.jld2"

simulation.output_writers[:slices] =
    JLD2Writer(model, merge(model.velocities, model.tracers, eddy_viscosity);
               filename,
               indices = (:, Ny ÷ 2, :),
               schedule = TimeInterval(1minute),
               overwrite_files = 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: [:coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

We're ready:

julia

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 (19.221 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (5.064 seconds).
Iteration: 0200, time: 15.999 minutes, Δt: 2.483 seconds, max(|w|) = 7.6e-06 ms⁻¹, wall time: 33.152 seconds
Iteration: 0400, time: 22.704 minutes, Δt: 1.717 seconds, max(|w|) = 4.7e-05 ms⁻¹, wall time: 40.914 seconds
Iteration: 0600, time: 27.730 minutes, Δt: 1.317 seconds, max(|w|) = 2.9e-03 ms⁻¹, wall time: 46.921 seconds
Iteration: 0800, time: 31.106 minutes, Δt: 782.504 ms, max(|w|) = 1.6e-01 ms⁻¹, wall time: 52.759 seconds
Iteration: 1000, time: 33.730 minutes, Δt: 994.684 ms, max(|w|) = 1.8e-01 ms⁻¹, wall time: 58.081 seconds
Iteration: 1200, time: 38.291 minutes, Δt: 1.718 seconds, max(|w|) = 1.3e-01 ms⁻¹, wall time: 1.067 minutes
Iteration: 1400, time: 44.407 minutes, Δt: 1.865 seconds, max(|w|) = 7.5e-02 ms⁻¹, wall time: 1.171 minutes
Iteration: 1600, time: 50.590 minutes, Δt: 1.847 seconds, max(|w|) = 7.1e-02 ms⁻¹, wall time: 1.271 minutes
Iteration: 1800, time: 56.513 minutes, Δt: 1.768 seconds, max(|w|) = 7.2e-02 ms⁻¹, wall time: 1.372 minutes
Iteration: 2000, time: 1.038 hours, Δt: 1.669 seconds, max(|w|) = 7.9e-02 ms⁻¹, wall time: 1.472 minutes
Iteration: 2200, time: 1.130 hours, Δt: 1.648 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 1.564 minutes
Iteration: 2400, time: 1.222 hours, Δt: 1.718 seconds, max(|w|) = 7.9e-02 ms⁻¹, wall time: 1.663 minutes
Iteration: 2600, time: 1.312 hours, Δt: 1.529 seconds, max(|w|) = 7.0e-02 ms⁻¹, wall time: 1.761 minutes
Iteration: 2800, time: 1.400 hours, Δt: 1.663 seconds, max(|w|) = 7.3e-02 ms⁻¹, wall time: 1.877 minutes
Iteration: 3000, time: 1.487 hours, Δt: 1.492 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 1.976 minutes
Iteration: 3200, time: 1.574 hours, Δt: 1.578 seconds, max(|w|) = 6.9e-02 ms⁻¹, wall time: 2.077 minutes
Iteration: 3400, time: 1.660 hours, Δt: 1.612 seconds, max(|w|) = 6.6e-02 ms⁻¹, wall time: 2.214 minutes
Iteration: 3600, time: 1.745 hours, Δt: 1.510 seconds, max(|w|) = 8.3e-02 ms⁻¹, wall time: 2.322 minutes
Iteration: 3800, time: 1.827 hours, Δt: 1.407 seconds, max(|w|) = 6.8e-02 ms⁻¹, wall time: 2.422 minutes
Iteration: 4000, time: 1.906 hours, Δt: 1.423 seconds, max(|w|) = 7.5e-02 ms⁻¹, wall time: 2.521 minutes
Iteration: 4200, time: 1.988 hours, Δt: 1.571 seconds, max(|w|) = 7.7e-02 ms⁻¹, wall time: 2.637 minutes
[ Info: Simulation is stopping after running for 2.655 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.

julia

time_series = (w = FieldTimeSeries(filename, "w"),
               T = FieldTimeSeries(filename, "T"),
               S = FieldTimeSeries(filename, "S"),
               νₑ = FieldTimeSeries(filename, "νₑ"))

(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 4×4×4 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: w
└── data: 136×1×73×121 StaticStridesOffsetArray(::Array{Float64, 4}, -3:132, 64:64, -3:69, 1:121) with eltype Float64 with indices -3:132×64:64×-3:69×1:121
    └── max=0.0851124, min=-0.212471, mean=-1.41625e-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 4×4×4 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: T
└── data: 136×1×72×121 StaticStridesOffsetArray(::Array{Float64, 4}, -3:132, 64:64, -3:68, 1:121) with eltype Float64 with indices -3:132×64:64×-3:68×1:121
    └── max=20.0007, min=0.0, mean=18.072, 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 4×4×4 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: S
└── data: 136×1×72×121 StaticStridesOffsetArray(::Array{Float64, 4}, -3:132, 64:64, -3:68, 1:121) with eltype Float64 with indices -3:132×64:64×-3:68×1:121
    └── max=35.033, min=0.0, mean=32.084, νₑ = 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 4×4×4 halo
├── indices: (:, 64:64, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: ocean_wind_mixing_and_convection.jld2
├── name: νₑ
└── data: 136×1×72×121 StaticStridesOffsetArray(::Array{Float64, 4}, -3:132, 64:64, -3:68, 1:121) with eltype Float64 with indices -3:132×64:64×-3:68×1:121
    └── max=0.0106354, min=0.0, mean=5.31883e-5)

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.

julia

times = time_series.w.times
n = Observable(length(times))

 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(Lx / Lz),
               limits = ((0, Lx), (-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 viscosity", 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. We start the animation at minutes since things are pretty boring till then:

julia

intro = searchsortedfirst(times, 10minutes)
frames = intro:length(times)

@info "Making a motion picture of ocean wind mixing and convection..."

CairoMakie.record(fig, "ocean_wind_mixing_and_convection.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:

julia

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.13.2
Commit 45b6d15e34b (2026-10-14 20:33 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-20.1.8 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
Environment:
  JULIA_LOAD_PATH = @:@v#.#:@stdlib
  JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.2+0.x64.linux.gnu/local/share/julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.2+0.x64.linux.gnu/share/julia
  JULIA_VERSION_ENZYME = 1.13.1
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  LD_LIBRARY_PATH = 
  JULIA_MAX_NUM_PRECOMPILE_FILES = 24
  JULIA_VERSION = 1.13.2
  JULIA_CUDA_USE_COMPAT = false
  JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-35948/docs/
  JULIA_DEBUG = Literate

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

julia

import Pkg
Pkg.status()

Status `~/Oceananigans.jl-35948/docs/Project.toml`
  [79e6a3ab] Adapt v4.8.0
⌃ [052768ef] CUDA v6.2.2
  [13f3f980] CairoMakie v0.15.17
  [8e50ac2c] ConservativeRegridding v0.2.20
  [e30172f5] Documenter v1.19.0
  [daee34ce] DocumenterCitations v1.5.0
  [4710194d] DocumenterVitepress v0.3.9
  [7da242da] Enzyme v0.13.218
  [033835bb] JLD2 v0.6.7
⌅ [63c18a36] KernelAbstractions v0.9.45
  [98b081ad] Literate v2.21.0
  [da04e1cc] MPI v0.20.29
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.117.1 `..`
  [429524aa] Optim v2.3.2
  [f27b6e38] Polynomials v4.1.6
  [3c362404] Reactant v0.2.305
  [6038ab10] Rotations v1.7.1
  [d496a93d] SeawaterPolynomials v0.3.11
  [09ab397b] StructArrays v0.7.4
  [bdfc003b] TimesDates v0.3.3
  [0a941bbe] Zarr v0.11.2
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.13.0
  [44cfe95a] Pkg v1.13.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.

`NonhydrostaticModel`

`HydrostaticFreeSurfaceModel`

`ShallowWaterModel`

Wind- and convection-driven mixing in an ocean surface boundary layer

Install dependencies

The grid

Buoyancy that depends on temperature and salinity

Boundary conditions

Model instantiation

Initial conditions

Setting up a simulation

Output

Turbulence visualization

Julia version and environment information

Wind- and convection-driven mixing in an ocean surface boundary layer ​

Install dependencies ​

The grid ​

Buoyancy that depends on temperature and salinity ​

Boundary conditions ​

Model instantiation ​

Initial conditions ​

Setting up a simulation ​

Output ​

Turbulence visualization ​

Julia version and environment information ​

Wind- and convection-driven mixing in an ocean surface boundary layer

Install dependencies

The grid

Buoyancy that depends on temperature and salinity

Boundary conditions

Model instantiation

Initial conditions

Setting up a simulation

Output

Turbulence visualization

Julia version and environment information