Momentum equation

The UMomentumEquation, VMomentumEquation, and WMomentumEquation modules provide diagnostics for every term on the right-hand side of the momentum equations. The momentum equation for the $i$-th velocity component $u_i$ is

\[\partial_t u_i = \underbrace{-\partial_j(u_j u_i)}_{\text{advection}} + \underbrace{\hat{g}_i \, b}_{\text{buoyancy}} + \underbrace{-\epsilon_{ijk} f_j u_k}_{\text{Coriolis}} + \underbrace{-\partial_i p}_{\text{pressure}} + \underbrace{-\partial_j \tau_{ij}}_{\text{viscous}} + \underbrace{(\nabla \times \mathbf{u}^S) \times \mathbf{u}}_{\text{Stokes shear}} + \underbrace{\partial_t u_i^S}_{\text{Stokes tendency}} + \underbrace{F_{u_i}}_{\text{forcing}}\]

where $\hat{g}_i$ is the $i$-th component of the gravitational unit vector, $b$ is the buoyancy, $f_j$ is the Coriolis frequency, $p$ is the pressure, $\tau_{ij}$ is the viscous/subgrid stress tensor, $\mathbf{u}^S$ is the Stokes drift, and $F_{u_i}$ is the forcing. This decomposition lets the user compute each contribution independently, build diagnostics like budget closure, or analyse the energetics of individual processes.

Each module wraps the corresponding Oceananigans velocity-tendency kernel and provides diagnostics at the natural grid location for that velocity component:

Module	Location	Tendency wrapped
`UMomentumEquation`	`(Face, Center, Center)`	`u_velocity_tendency` (NH) / `hydrostatic_free_surface_u_velocity_tendency` (HFS)
`VMomentumEquation`	`(Center, Face, Center)`	`v_velocity_tendency` (NH) / `hydrostatic_free_surface_v_velocity_tendency` (HFS)
`WMomentumEquation`	`(Center, Center, Face)`	`w_velocity_tendency` (NH only)

Constructors accept either a full Oceananigans model object (the convenience form) or the individual arguments expected by the underlying kernel function. The model-only form is enough for most analysis workflows.

Example

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b, coriolis=FPlane(f=1e-4));

julia> ADV = UMomentumEquation.Advection(model)
UAdvection KernelFunctionOperation at (Face, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: div_𝐯u (generic function with 4 methods)
└── arguments: ("Centered", "NamedTuple", "Field")
└── computes: advection of u-momentum  ∂ⱼ(uⱼu)

julia> BUOY = UMomentumEquation.BuoyancyAcceleration(model)
UBuoyancyAcceleration KernelFunctionOperation at (Face, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: x_dot_g_bᶠᶜᶜ (generic function with 3 methods)
└── arguments: ("BuoyancyForce", "NamedTuple")
└── computes: buoyancy acceleration (x)  ĝₓ b

julia> COR = VMomentumEquation.CoriolisAcceleration(model)
VCoriolisAcceleration KernelFunctionOperation at (Center, Face, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: y_f_cross_U (generic function with 11 methods)
└── arguments: ("FPlane", "NamedTuple")
└── computes: Coriolis acceleration (y)  (f⃗ × u⃗)_y

Budget closure

The combined RHS — Tendency(model) — wraps Oceananigans' internal *_velocity_tendency kernel directly. The sum of the individual term diagnostics, taken with the correct signs from Oceananigans' source, equals Tendency to machine precision.

`NonhydrostaticModel` u/v budget

\[\text{Tendency} = -\text{Advection} + \text{BuoyancyAcceleration} -\text{CoriolisAcceleration} -\text{PressureGradient} -\text{TotalViscousDissipation} +\text{StokesShear} + \text{StokesTendency} + \text{Forcing}\]

`NonhydrostaticModel` w budget

The w-momentum equation has one subtlety: Oceananigans' w_velocity_tendency has no explicit -\partial_z p line. The vertical hydrostatic balance is treated as a property of the pressure-projection step, and whether the buoyancy term appears in the tendency depends on whether the model uses hydrostatic/nonhydrostatic pressure splitting:

`hydrostatic_pressure_anomaly`	Buoyancy term in `Tendency`	Budget formula
`Field` (default, splitting on)	absorbed by $pHY'$	`-ADV - COR - TVISC + SS + ST + FORC`
`nothing` (splitting off)	included as $+b$	`-ADV + BUOY - COR - TVISC + SS + ST + FORC`

`HydrostaticFreeSurfaceModel` u/v budget

HFS uses hydrostatic_free_surface_{u,v}_velocity_tendency, which differs from the NH tendency: there are no Stokes terms, no buoyancy term (absorbed into the hydrostatic pressure anomaly pHY′), and the free-surface contribution to the horizontal pressure gradient appears as a separate BarotropicPressureGradient term:

\[\text{Tendency} = -\text{Advection} -\text{BarotropicPressureGradient} -\text{CoriolisAcceleration} -\text{PressureGradient} -\text{TotalViscousDissipation} +\text{Forcing}\]

For ImplicitFreeSurface (the default) and SplitExplicitFreeSurface, BarotropicPressureGradient returns zero (the contribution is handled inside the pressure solve), so the budget closes without it. For ExplicitFreeSurface it returns g ∂_i η and must be included for closure.

Two caveats worth noting:

Tendency(model) wraps Oceananigans' hydrostatic_free_surface_{u,v}_velocity_tendency, which returns G_u in the decomposition ∂_t u = G_u - ∂_x p_n. The implicit barotropic pressure correction p_n (relevant for ImplicitFreeSurface) is applied separately during time-stepping — not in G_u, and not surfaced as a diagnostic. The budget closure above is therefore against G_u, which is the physical balance of forces; reconstructing the full ∂_t u would additionally require the implicit correction applied by correct_barotropic_mode!.
On horizontally-curvilinear grids (LatitudeLongitudeGrid, OrthogonalSphericalShellGrid) with flux-form advection, hydrostatic_free_surface_*_velocity_tendency includes an extra U_dot_∇*_hydrostatic_metric term that is not surfaced as a diagnostic. This term is identically zero on RectilinearGrid (the configuration the tests use), on any grid with VectorInvariant advection (where the metric is already included in the vorticity/Bernoulli decomposition), and on Nothing advection — so budget closure holds for those configurations.

`HydrostaticFreeSurfaceModel` caveats

WMomentumEquation.Tendency, WMomentumEquation.Forcing(model, Val(:w)), and the U/V/W StokesShear and StokesTendency throw an ArgumentError on HydrostaticFreeSurfaceModel: w is diagnosed from continuity rather than evolved by a prognostic equation, and HFS has no stokes_drift field.
Advection(::HydrostaticFreeSurfaceModel) (and the explicit-args form Advection(::HydrostaticFreeSurfaceModel, u, v, w, scheme)) both wrap U_dot_∇u/U_dot_∇v — the vector-invariant kernel that HFS's tendency actually uses. This works with HFS's default VectorInvariant momentum advection as well as any flux-form scheme. The diagnostic returns a KFO whose kernel function is U_dot_∇u/U_dot_∇v, whereas the NH Advection diagnostic returns a KFO whose kernel function is div_𝐯u/div_𝐯v — both pass the isa Advection type check via the Advection Union alias.

`Forcing` dispatch quirk

Forcing aliases KernelFunctionOperation directly (no narrowing on the kernel function), so its constructor methods are shared across the three modules. To avoid clobbering UMomentumEquation.Forcing(model) when the V and W modules are also loaded, the V/W single-argument convenience form takes an explicit Val tag:

UMomentumEquation.Forcing(model)              # u-momentum forcing
VMomentumEquation.Forcing(model, Val(:v))     # v-momentum forcing
WMomentumEquation.Forcing(model, Val(:w))     # w-momentum forcing

The full explicit form Forcing(model, forcing, clock, model_fields, Val(:u/:v/:w)) is available in all three modules.

U-momentum

Oceanostics.UMomentumEquation.Advection — Type

Advection(model, velocities, advection_scheme; location)

Calculate the advection of u-momentum. For NonhydrostaticModel this wraps the flux-form divergence kernel div_𝐯u as

ADV = ∂ⱼ (uⱼ u)

For HydrostaticFreeSurfaceModel this wraps the vector-invariant kernel U_dot_∇u (which is what hydrostatic_free_surface_u_velocity_tendency actually uses) as

ADV = (u · ∇) u

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(grid);

julia> ADV = UMomentumEquation.Advection(model)
UAdvection KernelFunctionOperation at (Face, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: div_𝐯u (generic function with 4 methods)
└── arguments: ("Centered", "NamedTuple", "Field")
└── computes: advection of u-momentum  ∂ⱼ(uⱼu)