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 allows 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 Type

NonhydrostaticModel(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 = nothing,
                    hydrostatic_pressure_anomaly = DefaultHydrostaticPressureAnomaly(),
                    closure_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.

Arguments

• grid: (required) The resolution and discrete geometry on which the model is solved. The architecture (CPU/GPU) that the model is solved on is inferred from the architecture of the grid. Note that the grid needs to be regularly spaced in the horizontal dimensions, and .

Keyword arguments

• advection: The scheme that advects velocities and tracers. See Oceananigans.Advection.

• buoyancy: The buoyancy model. See Oceananigans.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 for model. See Oceananigans.TurbulenceClosures.

• boundary_conditions: NamedTuple containing field boundary conditions.

• tracers: A tuple of symbols defining the names of the modeled tracers, or a NamedTuple of preallocated CenterFields.

• timestepper: A symbol or a TimeStepper object that specifies the time-stepping method. Supported symbols include :QuasiAdamsBashforth2, :RungeKutta3. Default: :RungeKutta3.

• 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 for tracers.

• velocities: The model velocities. Default: nothing.

• nonhydrostatic_pressure: The nonhydrostatic pressure field. Default: nothing.

• hydrostatic_pressure_anomaly: An optional field that stores the part of the nonhydrostatic pressure in hydrostatic balance with the buoyancy field. If CenterField(grid), the anomaly is precomputed by vertically integrating the buoyancy field. In this case, the nonhydrostatic_pressure represents only the part of pressure that deviates from the hydrostatic anomaly. If nothing (default), the anomaly is not computed. Note: for grids with periodic vertical topology, the hydrostatic pressure anomaly is set to nothing by default.

• closure_fields: Diffusivity fields. Default: nothing.

• pressure_solver: Pressure solver to be used in the model. If nothing (default), the model constructor chooses the default based on the grid provide.

• auxiliary_fields: NamedTuple of auxiliary fields. Default: nothing

source

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.00513064, min=-0.00475115, mean=2.87991e-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
├── run_wall_time: 1.973 seconds
├── run_wall_time / iteration: 197.259 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 => 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

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 Type

HydrostaticFreeSurfaceModel(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,
                            closure_fields = nothing,
                            auxiliary_fields = NamedTuple(),
                            vertical_coordinate = default_vertical_coordinate(grid))

Construct a hydrostatic model with a free surface on grid.

Arguments

• grid: (required) The resolution and discrete geometry on which model is solved. The architecture (CPU/GPU) that the model is solved is inferred from the architecture of the grid.

Keyword arguments

• momentum_advection: The scheme that advects velocities. See Oceananigans.Advection.

• tracer_advection: The scheme that advects tracers. See Oceananigans.Advection.

• buoyancy: The buoyancy model. See Oceananigans.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 the grid. If the grid is a RectilinearGrid that is regularly spaced in the horizontal the default is an ImplicitFreeSurface solver with solver_method = :FFTBasedPoissonSolver. In all other cases, the default is a SplitExplicitFreeSurface.

• tracers: A tuple of symbols defining the names of the modeled tracers, or a NamedTuple of preallocated CenterFields.

• forcing: NamedTuple of user-defined forcing functions that contribute to solution tendencies.

• closure: The turbulence closure for model. See Oceananigans.TurbulenceClosures.

• timestepper: A symbol or a TimeStepper object that specifies the time-stepping method. Supported symbols include :QuasiAdamsBashforth2, :SplitRungeKutta2, :SplitRungeKutta3, :SplitRungeKutta4, :SplitRungeKutta5. Default: :QuasiAdamsBashforth2.

• boundary_conditions: NamedTuple containing field boundary conditions.

• particles: Lagrangian particles to be advected with the flow. Default: nothing.

• biogeochemistry: Biogeochemical model for tracers.

• velocities: The model velocities. Default: nothing.

• pressure: Hydrostatic pressure field. Default: nothing.

• closure_fields: Closure fields. Default: nothing.

• auxiliary_fields: NamedTuple of auxiliary fields. Default: nothing.

• vertical_coordinate: Algorithm for grid evolution: ZStarCoordinate() or ZCoordinate(grid). Default: default_vertical_coordinate(grid), which returns ZStarCoordinate(grid) for grids with MutableVerticalDiscretization otherwise returns ZCoordinate().

source

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

HydrostaticFreeSurfaceModel{CPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0)
├── grid: 180×80×10 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 6×6×3 halo
├── 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, Oceananigans.Coriolis.HydrostaticFormulation}

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
├── 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.8866, min=16.1633, 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 and HydrostaticFreeSurfaceModel.