Models
In general in Oceananigans, the "model" object serves two main purposes:
- models store the configuration of a set of discrete equations. The discrete equations imply rules for evolving prognostic variables, and computing diagnostic variables from the prognostic state.
- models provide a container for the prognostic and diagnostic state of those discrete equations at a particular time.
Two Oceananigans models for ocean simulations
In addition to defining the abstract concept of a "model" that can be used with Simulation, Oceananigans provides two mature model implementations for simulating ocean-flavored fluid dynamics. Both of these integrate the Navier-Stokes equations within the Boussinesq approximation (we call these the "Boussinesq equations" for short): the NonhydrostaticModel and the HydrostaticFreeSurfaceModel.
The NonhydrostaticModel integrates the Boussinesq equations without making the hydrostatic approximation, and therefore possessing a prognostic vertical momentum equation. The NonhydrostaticModel is useful for simulations that resolve three-dimensional turbulence, such as large eddy simulations on RectilinearGrid with grid spacings of O(1 m), as well as direct numerical simulation. The NonhydrostaticModel may also be used for idealized classroom problems, as in the two-dimensional turbulence example.
The HydrostaticFreeSurfaceModel integrates the hydrostatic or "primitive" Boussinesq equations with a free surface on its top boundary. The hydrostatic approximation allosw the HydrostaticFreeSurfaceModel to achieve much higher efficiency in simulations on curvilinear grids used for large-scale regional or global simulations such as LatitudeLongitudeGrid, TripolarGrid, ConformalCubedSphereGrid, and other OrthogonalSphericalShellGrids such as RotatedLatitudeLongitudeGrid. Because they span larger domains, simulations with the HydrostaticFreeSurfaceModel also usually involve coarser grid spacings of O(30 m) up to O(100 km). Such coarse-grained simulations are usually paired with more elaborate turbulence closures or "parameterizations" than small-scale simulations with NonhydrostaticModel, such as the vertical mixing schemes CATKEVerticalDiffusivity, RiBasedVerticalDiffusivity, and TKEDissipationVerticalDiffusivity, and the mesoscale turbulence closure IsopycnalSkewSymmetricDiffusivity (a.k.a. "Gent-McWilliams plus Redi").
A third, experimental ShallowWaterModel solves the shallow water equations.
Configuring NonhydrostaticModel with keyword arguments
To illustrate the specification of discrete equations, consider first the docstring for NonhydrostaticModel,
Oceananigans.Models.NonhydrostaticModels.NonhydrostaticModel
— TypeNonhydrostaticModel(; grid,
clock = Clock{eltype(grid)}(time = 0),
advection = Centered(),
buoyancy = nothing,
coriolis = nothing,
stokes_drift = nothing,
forcing::NamedTuple = NamedTuple(),
closure = nothing,
boundary_conditions::NamedTuple = NamedTuple(),
tracers = (),
timestepper = :RungeKutta3,
background_fields::NamedTuple = NamedTuple(),
particles::ParticlesOrNothing = nothing,
biogeochemistry::AbstractBGCOrNothing = nothing,
velocities = nothing,
nonhydrostatic_pressure = CenterField(grid),
hydrostatic_pressure_anomaly = DefaultHydrostaticPressureAnomaly(),
diffusivity_fields = nothing,
pressure_solver = nothing,
auxiliary_fields = NamedTuple())
Construct a model for a non-hydrostatic, incompressible fluid on grid
, using the Boussinesq approximation when buoyancy != nothing
. By default, all Bounded directions are rigid and impenetrable.
Keyword arguments
grid
: (required) The resolution and discrete geometry on which themodel
is solved. The architecture (CPU/GPU) that the model is solved on is inferred from the architecture of thegrid
. Note that the grid needs to be regularly spaced in the horizontal dimensions, $x$ and $y$.advection
: The scheme that advects velocities and tracers. SeeOceananigans.Advection
.buoyancy
: The buoyancy model. SeeOceananigans.BuoyancyFormulations
.coriolis
: Parameters for the background rotation rate of the model.stokes_drift
: Parameters for Stokes drift fields associated with surface waves. Default:nothing
.forcing
:NamedTuple
of user-defined forcing functions that contribute to solution tendencies.closure
: The turbulence closure formodel
. SeeOceananigans.TurbulenceClosures
.boundary_conditions
:NamedTuple
containing field boundary conditions.tracers
: A tuple of symbols defining the names of the modeled tracers, or aNamedTuple
of preallocatedCenterField
s.timestepper
: A symbol that specifies the time-stepping method. Either:QuasiAdamsBashforth2
or:RungeKutta3
(default).background_fields
:NamedTuple
with background fields (e.g., background flow). Default:nothing
.particles
: Lagrangian particles to be advected with the flow. Default:nothing
.biogeochemistry
: Biogeochemical model fortracers
.velocities
: The model velocities. Default:nothing
.nonhydrostatic_pressure
: The nonhydrostatic pressure field. Default:CenterField(grid)
.hydrostatic_pressure_anomaly
: An optional field that stores the part of the nonhydrostatic pressure in hydrostatic balance with the buoyancy field. IfCenterField(grid)
(default), the anomaly is precomputed by vertically integrating the buoyancy field. In this case, thenonhydrostatic_pressure
represents only the part of pressure that deviates from the hydrostatic anomaly. Ifnothing
, the anomaly is not computed.diffusivity_fields
: Diffusivity fields. Default:nothing
.pressure_solver
: Pressure solver to be used in the model. Ifnothing
(default), the model constructor chooses the default based on thegrid
provide.auxiliary_fields
:NamedTuple
of auxiliary fields. Default:nothing
The configuration operations for NonhydrostaticModel include "discretization choices", such as the advection
scheme, as well as aspects of the continuous underlying equations, such as the formulation of the buoyancy force.
For our first example, we build the default NonhydrostaticModel
(which is quite boring):
using Oceananigans
grid = RectilinearGrid(size=(8, 8, 8), extent=(8, 8, 8))
nh = NonhydrostaticModel(; grid)
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: ()
├── closure: Nothing
├── buoyancy: Nothing
└── coriolis: Nothing
The default NonhydrostaticModel
has no tracers, no buoyancy force, no Coriolis force, and a second-order advection scheme. We next consider a slightly more exciting NonhydrostaticModel configured with a WENO advection scheme, the temperature/salinity-based SeawaterBuoyancy, a boundary condition on the zonal momentum, and a passive tracer forced by a cooked-up surface flux called "c":
using Oceananigans
grid = RectilinearGrid(size=(128, 128), halo=(5, 5), x=(0, 256), z=(-128, 0),
topology = (Periodic, Flat, Bounded))
# Numerical method and physics choices
advection = WENO(order=9) # ninth‑order upwind for momentum and tracers
buoyancy = BuoyancyTracer()
coriolis = FPlane(f=1e-4)
τx = - 8e-5 # m² s⁻² (τₓ < 0 ⟹ eastward wind stress)
u_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(τx))
@inline Jc(x, t, Lx) = cos(2π / Lx * x)
c_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(Jc, parameters=grid.Lx))
model = NonhydrostaticModel(; grid, advection, buoyancy, coriolis,
tracers = (:b, :c),
boundary_conditions = (; u=u_bcs, c=c_bcs))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 5×0×5 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: WENO{5, Float64, Float32}(order=9)
├── tracers: (b, c)
├── closure: Nothing
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
└── coriolis: FPlane{Float64}(f=0.0001)
Mutation of the model state
In addition to providing an interface for configuring equations, models also store the prognostic and diagnostic state associated with the solution to those equations. Models thus also provide an interface for "setting" or fixing the prognostic state, which is typically invoked to determine the initial conditions of a simulation. To illustrate this we consider setting the above model to a stably-stratified and noisy condition:
N² = 1e-5
bᵢ(x, z) = N² * z + 1e-6 * randn()
uᵢ(x, z) = 1e-3 * randn()
set!(model, b=bᵢ, u=uᵢ, w=uᵢ)
model.tracers.b
128×1×128 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 128×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 5×0×5 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 138×1×138 OffsetArray(::Array{Float64, 3}, -4:133, 1:1, -4:133) with eltype Float64 with indices -4:133×1:1×-4:133
└── max=-2.59705e-6, min=-0.00127806, mean=-0.000640007
Invoking set!
above determine the model tracer b
and the velocity components u
and w
. set!
also computes the diagnostic state of a model, which in the case of NonhydrostaticModel
includes the nonhydrostatic component of pressure,
model.pressures.pNHS
128×1×128 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 128×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 5×0×5 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 138×1×138 OffsetArray(::Array{Float64, 3}, -4:133, 1:1, -4:133) with eltype Float64 with indices -4:133×1:1×-4:133
└── max=0.00492612, min=-0.00468784, mean=3.7693e-20
Evolving models in time
Model may be integrated or "stepped forward" in time by calling time_step!(model, Δt)
, where Δt
is the time step and thus advancing the model.clock
:
time_step!(model, 1)
model.clock
Clock{Float64, Float64}(time=1 second, iteration=1, last_Δt=1 second)
├── stage: 1
└── last_stage_Δt: 333.333 ms
However, users are strongly encouraged to use the Simulation
interface to manage time-stepping along with other activities, like monitoring progress, writing output to disk, and more.
simulation = Simulation(model, Δt=1, stop_iteration=10)
run!(simulation)
simulation
Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 10 seconds, iteration = 10)
├── Next time step: 1 second
├── Elapsed wall time: 2.950 seconds
├── Wall time per iteration: 294.967 ms
├── Stop time: Inf days
├── Stop iteration: 10.0
├── Wall time limit: Inf
├── Minimum relative step: 0.0
├── Callbacks: OrderedDict with 4 entries:
│ ├── stop_time_exceeded => 4
│ ├── stop_iteration_exceeded => -
│ ├── wall_time_limit_exceeded => e
│ └── nan_checker => }
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entries
Using the HydrostaticFreeSurfaceModel
The HydrostaticFreeSurfaceModel has a similar interface as the NonhydrostaticModel,
using Oceananigans
grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1))
model = HydrostaticFreeSurfaceModel(; grid) # default free surface, no tracers
HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: ()
├── closure: Nothing
├── buoyancy: Nothing
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│ └── solver: FFTImplicitFreeSurfaceSolver
├── advection scheme:
│ └── momentum: VectorInvariant
├── vertical_coordinate: ZCoordinate
└── coriolis: Nothing
The full array of keyword arguments used to configure a HydrostaticFreeSurfaceModel are detailed in the docstring for HydrostaticFreeSurfaceModel,
Oceananigans.Models.HydrostaticFreeSurfaceModels.HydrostaticFreeSurfaceModel
— TypeHydrostaticFreeSurfaceModel(; grid,
clock = Clock{Float64}(time = 0),
momentum_advection = VectorInvariant(),
tracer_advection = Centered(),
buoyancy = SeawaterBuoyancy(eltype(grid)),
coriolis = nothing,
free_surface = [default_free_surface],
forcing::NamedTuple = NamedTuple(),
closure = nothing,
timestepper = :QuasiAdamsBashforth2,
boundary_conditions::NamedTuple = NamedTuple(),
tracers = (:T, :S),
particles::ParticlesOrNothing = nothing,
biogeochemistry::AbstractBGCOrNothing = nothing,
velocities = nothing,
pressure = nothing,
diffusivity_fields = nothing,
auxiliary_fields = NamedTuple(),
vertical_coordinate = default_vertical_coordinate(grid))
Construct a hydrostatic model with a free surface on grid
.
Keyword arguments
grid
: (required) The resolution and discrete geometry on whichmodel
is solved. The architecture (CPU/GPU) that the model is solved is inferred from the architecture of thegrid
.momentum_advection
: The scheme that advects velocities. SeeOceananigans.Advection
.tracer_advection
: The scheme that advects tracers. SeeOceananigans.Advection
.buoyancy
: The buoyancy model. SeeOceananigans.BuoyancyFormulations
.coriolis
: Parameters for the background rotation rate of the model.free_surface
: The free surface model. The default free-surface solver depends on the geometry of thegrid
. If thegrid
is aRectilinearGrid
that is regularly spaced in the horizontal the default is anImplicitFreeSurface
solver withsolver_method = :FFTBasedPoissonSolver
. In all other cases, the default is aSplitExplicitFreeSurface
.tracers
: A tuple of symbols defining the names of the modeled tracers, or aNamedTuple
of preallocatedCenterField
s.forcing
:NamedTuple
of user-defined forcing functions that contribute to solution tendencies.closure
: The turbulence closure formodel
. SeeOceananigans.TurbulenceClosures
.timestepper
: A symbol that specifies the time-stepping method. Either:QuasiAdamsBashforth2
(default) or:SplitRungeKutta3
.boundary_conditions
:NamedTuple
containing field boundary conditions.particles
: Lagrangian particles to be advected with the flow. Default:nothing
.biogeochemistry
: Biogeochemical model fortracers
.velocities
: The model velocities. Default:nothing
.pressure
: Hydrostatic pressure field. Default:nothing
.diffusivity_fields
: Diffusivity fields. Default:nothing
.auxiliary_fields
:NamedTuple
of auxiliary fields. Default:nothing
.vertical_coordinate
: Algorithm for grid evolution:ZStarCoordinate()
orZCoordinate(grid)
. Default:default_vertical_coordinate(grid)
, which returnsZStarCoordinate(grid)
for grids withMutableVerticalDiscretization
otherwise returnsZCoordinate()
.
A bit more involved HydrostaticFreeSurfaceModel example:
using Oceananigans
using SeawaterPolynomials: TEOS10EquationOfState
grid = LatitudeLongitudeGrid(size = (180, 80, 10),
longitude = (0, 360),
latitude = (-80, 80),
z = (-1000, 0),
halo = (6, 6, 3))
momentum_advection = WENOVectorInvariant()
coriolis = HydrostaticSphericalCoriolis()
equation_of_state = TEOS10EquationOfState()
buoyancy = SeawaterBuoyancy(; equation_of_state)
closure = CATKEVerticalDiffusivity()
# Generate a zonal wind stress that mimics Earth's mean winds
# with westerlies in mid-latitudes and easterlies near equator and poles
function zonal_wind_stress(λ, φ, t)
# Parameters
τ₀ = 1e-4 # Maximum wind stress magnitude (N/m²)
φ₀ = 30 # Latitude of maximum westerlies (degrees)
dφ = 10
# Approximate wind stress pattern
return - τ₀ * (+ exp(-(φ - φ₀)^2 / 2dφ^2)
- exp(-(φ + φ₀)^2 / 2dφ^2)
- 0.3 * exp(-φ^2 / dφ^2))
end
u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(zonal_wind_stress))
model = HydrostaticFreeSurfaceModel(; grid, momentum_advection, coriolis, closure, buoyancy,
boundary_conditions = (; u=u_bcs), tracers=(:T, :S, :e))
HydrostaticFreeSurfaceModel{CPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0)
├── grid: 180×80×10 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 6×6×3 halo and with precomputed metrics
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: (T, S, e)
├── closure: CATKEVerticalDiffusivity{VerticallyImplicitTimeDiscretization}
├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
├── free surface: SplitExplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│ └── substepping: FixedTimeStepSize(4.911 minutes)
├── advection scheme:
│ ├── momentum: WENOVectorInvariant{5, Float64, Float32}(vorticity_order=9, vertical_order=5)
│ ├── T: Centered(order=2)
│ ├── S: Centered(order=2)
│ └── e: Centered(order=2)
├── vertical_coordinate: ZCoordinate
└── coriolis: HydrostaticSphericalCoriolis{Oceananigans.Advection.EnstrophyConserving{Float64}, Float64}
Mutating the state of the HydrostaticFreeSurfaceModel works similarly as for the NonhydrostaticModel –- except that the vertical velocity cannot be set!
, because vertical velocity is not prognostic in the hydrostatic equations.
using SeawaterPolynomials
N² = 1e-5
T₀ = 20
S₀ = 35
eos = model.buoyancy.formulation.equation_of_state
α = SeawaterPolynomials.thermal_expansion(T₀, S₀, 0, eos)
g = model.buoyancy.formulation.gravitational_acceleration
dTdz = N² / (α * g)
Tᵢ(λ, φ, z) = T₀ + dTdz * z + 1e-3 * T₀ * randn()
uᵢ(λ, φ, z) = 1e-3 * randn()
set!(model, T=Tᵢ, S=S₀, u=uᵢ, v=uᵢ)
model.tracers.T
180×80×10 Field{Center, Center, Center} on LatitudeLongitudeGrid on CPU
├── grid: 180×80×10 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 6×6×3 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: ZeroFlux, north: ZeroFlux, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 192×92×16 OffsetArray(::Array{Float64, 3}, -5:186, -5:86, -2:13) with eltype Float64 with indices -5:186×-5:86×-2:13
└── max=19.8764, min=16.1734, mean=18.0249
Where to go next
- See Quick start for a compact, end-to-end workflow
- See the Examples gallery for longer tutorials covering specific cases, including large eddy simulation, Kelvin–Helmholtz instability and baroclinic instability.
- Other pages in the Models section: in‑depth pages on buoyancy, forcing, boundary conditions, closures, diagnostics, and output.
- Physics: governing equations and numerical forms for
NonhydrostaticModel
andHydrostaticFreeSurfaceModel
.