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,

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 = CenterField(grid),
         hydrostatic_pressure_anomaly = DefaultHydrostaticPressureAnomaly(),
                   diffusivity_fields = nothing,
                      pressure_solver = nothing,
                     auxiliary_fields = NamedTuple())

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
├── run_wall_time: 3.768 seconds
├── run_wall_time / iteration: 376.812 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,

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

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