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:

ModuleLocationTendency 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, Oceanosticsjulia> 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)  ĝₓ bjulia> 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_anomalyBuoyancy term in TendencyBudget 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), 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.

U-momentum

Oceanostics.UMomentumEquation.AdvectionType
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, Oceanosticsjulia> 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)
source
Oceanostics.UMomentumEquation.BuoyancyAccelerationType
BuoyancyAcceleration(model, buoyancy, tracers; location)

Calculate the buoyancy acceleration in the x-direction as

BUOY = ĝₓ b

where ĝₓ is the x-component of the gravitational unit vector and b is the buoyancy.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b);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
source
Oceanostics.UMomentumEquation.CoriolisAccelerationType
CoriolisAcceleration(model, coriolis, velocities; location)

Calculate the Coriolis acceleration in the x-direction as

COR = f × u

where f is the Coriolis parameter vector and u is the velocity vector.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid; coriolis=FPlane(f=1e-4));julia> COR = UMomentumEquation.CoriolisAcceleration(model)UCoriolisAcceleration KernelFunctionOperation at (Face, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: x_f_cross_U (generic function with 11 methods)└── arguments: ("FPlane", "NamedTuple")└── computes: Coriolis acceleration (x)  (f⃗ × u⃗)ₓ
source
Oceanostics.UMomentumEquation.PressureGradientType
PressureGradient(model, hydrostatic_pressure; location)

Calculate the hydrostatic pressure gradient force in the x-direction as

PRES = ∂p/∂x

where p is the hydrostatic pressure anomaly.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> PRES = UMomentumEquation.PressureGradient(model)UPressureGradient KernelFunctionOperation at (Face, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: hydrostatic_pressure_gradient_x (generic function with 2 methods)└── arguments: ("Nothing",)└── computes: hydrostatic pressure gradient (x)  ∂p/∂x
source
Oceanostics.UMomentumEquation.BarotropicPressureGradientType
BarotropicPressureGradient(model, free_surface; location)

Calculate the explicit barotropic pressure gradient in the x-direction. For HydrostaticFreeSurfaceModel this is the contribution of the free surface $η$ to the horizontal pressure gradient that's treated explicitly during time-stepping (i.e., not solved implicitly). For ExplicitFreeSurface it equals g ∂x η; for ImplicitFreeSurface or SplitExplicitFreeSurface the contribution is handled inside the pressure solve so this kernel returns zero. For NonhydrostaticModel (free_surface == nothing) it also returns zero. This term completes the u-momentum budget closure on HFS.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> BARO = UMomentumEquation.BarotropicPressureGradient(model)UBarotropicPressureGradient KernelFunctionOperation at (Face, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: explicit_barotropic_pressure_x_gradient (generic function with 4 methods)└── arguments: ("Nothing",)└── computes: barotropic pressure gradient (x)  g ∂η/∂x
source
Oceanostics.UMomentumEquation.ViscousDissipationType
ViscousDissipation(
    model,
    closure,
    diffusivities,
    clock,
    model_fields,
    buoyancy;
    location
)

Calculate the viscous dissipation term (excluding immersed boundaries) as

VISC = ∂ⱼ τ₁ⱼ,

where τ₁ⱼ is the viscous stress tensor for the x-momentum equation.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> VISC = UMomentumEquation.ViscousDissipation(model)UViscousDissipation KernelFunctionOperation at (Face, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: ∂ⱼ_τ₁ⱼ (generic function with 10 methods)└── arguments: ("Nothing", "Nothing", "Clock", "NamedTuple", "Nothing")└── computes: viscous term (interior, x)  ∂ⱼτ₁ⱼ
source
Oceanostics.UMomentumEquation.ImmersedViscousDissipationType
ImmersedViscousDissipation(
    model,
    velocities,
    u_immersed_bc,
    closure,
    diffusivities,
    clock,
    model_fields;
    location
)

Calculate the viscous dissipation term due to immersed boundaries as

VISC = ∂ⱼ τ₁ⱼ,

where τ₁ⱼ is the immersed boundary viscous stress tensor for the x-momentum equation.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> VISC = UMomentumEquation.ImmersedViscousDissipation(model)UImmersedViscousDissipation KernelFunctionOperation at (Face, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: immersed_∂ⱼ_τ₁ⱼ (generic function with 2 methods)└── arguments: ("NamedTuple", "Nothing", "Nothing", "Nothing", "Clock", "NamedTuple")└── computes: viscous term through immersed boundaries (x)  ∂ⱼτ₁ⱼ
source
Oceanostics.UMomentumEquation.TotalViscousDissipationType
TotalViscousDissipation(
    model,
    velocities,
    u_immersed_bc,
    closure,
    diffusivities,
    clock,
    model_fields,
    buoyancy;
    location
)

Calculate the total viscous dissipation term as

VISC = ∂ⱼ τ₁ⱼ + ∂ⱼ τ₁ⱼ_immersed,

where τ₁ⱼ is the interior viscous stress tensor and τ₁ⱼ_immersed is the immersed boundary viscous stress tensor for the x-momentum equation.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> VISC = UMomentumEquation.TotalViscousDissipation(model)UTotalViscousDissipation KernelFunctionOperation at (Face, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: total_∂ⱼ_τ₁ⱼ (generic function with 1 method)└── arguments: ("NamedTuple", "Nothing", "Nothing", "Nothing", "Clock", "NamedTuple", "Nothing")└── computes: total viscous term (interior + immersed, x)  ∂ⱼτ₁ⱼ
source
Oceanostics.UMomentumEquation.StokesShearType
StokesShear(model, stokes_drift, velocities, time; location)

Calculate the Stokes shear term as

STOKES_SHEAR = (∇ × uˢ) × u

where uˢ is the Stokes drift velocity and u is the velocity vector.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> STOKES = UMomentumEquation.StokesShear(model)UStokesShear KernelFunctionOperation at (Face, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: x_curl_Uˢ_cross_U (generic function with 3 methods)└── arguments: ("Nothing", "NamedTuple", "Float64")└── computes: Stokes shear forcing (x)  ((∇ × u⃗ˢ) × u⃗)ₓ
source
Oceanostics.UMomentumEquation.StokesTendencyType
StokesTendency(model, stokes_drift, time; location)

Calculate the Stokes tendency term as

STOKES_TEND = ∂uˢ/∂t

where uˢ is the Stokes drift velocity.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> STOKES = UMomentumEquation.StokesTendency(model)UStokesTendency KernelFunctionOperation at (Face, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: ∂t_uˢ (generic function with 4 methods)└── arguments: ("Nothing", "Float64")└── computes: Stokes drift tendency (x)  ∂uˢ/∂t
source
Oceanostics.UMomentumEquation.ForcingType
Forcing(model, forcing, clock, model_fields; location)

Calculate the forcing term Fᵘ on the x-momentum equation for model.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> FORC = UMomentumEquation.Forcing(model)UForcing KernelFunctionOperation at (Face, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: forcing_fcc (generic function with 1 method)└── arguments: ("Returns", "Clock", "NamedTuple")└── computes: momentum forcing (x)  Fᵘ
source
Oceanostics.UMomentumEquation.TendencyType
Tendency(
    model,
    advection_scheme,
    coriolis,
    closure,
    u_immersed_bc,
    velocities,
    free_surface,
    tracers,
    buoyancy,
    closure_fields,
    hydrostatic_pressure_anomaly,
    auxiliary_fields,
    vertical_coordinate,
    clock,
    forcing_func;
    location
)

Calculate the u-momentum tendency Gᵘ as computed by Oceananigans, where Oceananigans writes the momentum equation as ∂_t u = Gᵘ - ∂ₓ p_n. Gᵘ is the sum of every term in u_velocity_tendency (NH) or hydrostatic_free_surface_u_velocity_tendency (HFS): advection, buoyancy (NH only — absorbed into the hydrostatic pressure for HFS), Coriolis, hydrostatic pressure gradient, viscous and immersed-viscous stress divergence, Stokes shear and tendency (NH only), and forcing. The terms exposed by UMomentumEquation reconstruct Gᵘ to machine precision.

Gᵘ is not the full time derivative ∂_t u. The remaining piece -∂ₓ p_n — the nonhydrostatic pressure gradient for NonhydrostaticModel, or the implicit barotropic free-surface correction for HydrostaticFreeSurfaceModel with ImplicitFreeSurface — is applied separately by the pressure-projection / barotropic-correction step of the time-stepper and is not part of Gᵘ or of any diagnostic in this module.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> TEND = UMomentumEquation.Tendency(model)UTendency KernelFunctionOperation at (Face, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: u_velocity_tendency (generic function with 1 method)└── arguments: ("Centered", "Nothing", "Nothing", "Nothing", "Nothing", "Nothing", "Oceananigans.Models.NonhydrostaticModels.BackgroundFields", "NamedTuple", "NamedTuple", "NamedTuple", "Nothing", "Nothing", "Clock", "Returns")└── computes: total tendency of the u-momentum equation
source

V-momentum

Oceanostics.VMomentumEquation.AdvectionType
Advection(model, velocities, advection_scheme; location)

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

ADV = ∂ⱼ (uⱼ v)

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

ADV = (u · ∇) v
julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> ADV = VMomentumEquation.Advection(model)VAdvection KernelFunctionOperation at (Center, Face, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: div_𝐯v (generic function with 4 methods)└── arguments: ("Centered", "NamedTuple", "Field")└── computes: advection of v-momentum  ∂ⱼ(uⱼv)
source
Oceanostics.VMomentumEquation.BuoyancyAccelerationType
BuoyancyAcceleration(model, buoyancy, tracers; location)

Calculate the buoyancy acceleration in the y-direction as

BUOY = ĝᵧ b

where ĝᵧ is the y-component of the gravitational unit vector and b is the buoyancy.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b);julia> BUOY = VMomentumEquation.BuoyancyAcceleration(model)VBuoyancyAcceleration KernelFunctionOperation at (Center, Face, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: y_dot_g_bᶜᶠᶜ (generic function with 3 methods)└── arguments: ("BuoyancyForce", "NamedTuple")└── computes: buoyancy acceleration (y)  ĝ_y b
source
Oceanostics.VMomentumEquation.CoriolisAccelerationType
CoriolisAcceleration(model, coriolis, velocities; location)

Calculate the Coriolis acceleration in the y-direction as

COR = (f × u)ᵧ

where f is the Coriolis parameter vector and u is the velocity vector.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid; coriolis=FPlane(f=1e-4));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
source
Oceanostics.VMomentumEquation.PressureGradientType
PressureGradient(model, hydrostatic_pressure; location)

Calculate the hydrostatic pressure gradient force in the y-direction as

PRES = ∂p/∂y

where p is the hydrostatic pressure anomaly.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> PRES = VMomentumEquation.PressureGradient(model)VPressureGradient KernelFunctionOperation at (Center, Face, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: hydrostatic_pressure_gradient_y (generic function with 2 methods)└── arguments: ("Nothing",)└── computes: hydrostatic pressure gradient (y)  ∂p/∂y
source
Oceanostics.VMomentumEquation.BarotropicPressureGradientType
BarotropicPressureGradient(model, free_surface; location)

Calculate the explicit barotropic pressure gradient in the y-direction. For HydrostaticFreeSurfaceModel this is the contribution of the free surface $η$ to the horizontal pressure gradient that's treated explicitly during time-stepping (i.e., not solved implicitly). For ExplicitFreeSurface it equals g ∂y η; for ImplicitFreeSurface or SplitExplicitFreeSurface the contribution is handled inside the pressure solve so this kernel returns zero. For NonhydrostaticModel (free_surface == nothing) it also returns zero. This term completes the v-momentum budget closure on HFS.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> BARO = VMomentumEquation.BarotropicPressureGradient(model)VBarotropicPressureGradient KernelFunctionOperation at (Center, Face, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: explicit_barotropic_pressure_y_gradient (generic function with 4 methods)└── arguments: ("Nothing",)└── computes: barotropic pressure gradient (y)  g ∂η/∂y
source
Oceanostics.VMomentumEquation.ViscousDissipationType
ViscousDissipation(
    model,
    closure,
    diffusivities,
    clock,
    model_fields,
    buoyancy;
    location
)

Calculate the viscous dissipation term (excluding immersed boundaries) as

VISC = ∂ⱼ τ₂ⱼ,

where τ₂ⱼ is the viscous stress tensor for the y-momentum equation.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> VISC = VMomentumEquation.ViscousDissipation(model)VViscousDissipation KernelFunctionOperation at (Center, Face, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: ∂ⱼ_τ₂ⱼ (generic function with 10 methods)└── arguments: ("Nothing", "Nothing", "Clock", "NamedTuple", "Nothing")└── computes: viscous term (interior, y)  ∂ⱼτ₂ⱼ
source
Oceanostics.VMomentumEquation.ImmersedViscousDissipationType
ImmersedViscousDissipation(
    model,
    velocities,
    v_immersed_bc,
    closure,
    diffusivities,
    clock,
    model_fields;
    location
)

Calculate the viscous dissipation term due to immersed boundaries as

VISC = ∂ⱼ τ₂ⱼ,

where τ₂ⱼ is the immersed boundary viscous stress tensor for the y-momentum equation.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> VISC = VMomentumEquation.ImmersedViscousDissipation(model)VImmersedViscousDissipation KernelFunctionOperation at (Center, Face, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: immersed_∂ⱼ_τ₂ⱼ (generic function with 2 methods)└── arguments: ("NamedTuple", "Nothing", "Nothing", "Nothing", "Clock", "NamedTuple")└── computes: viscous term through immersed boundaries (y)  ∂ⱼτ₂ⱼ
source
Oceanostics.VMomentumEquation.TotalViscousDissipationType
TotalViscousDissipation(
    model,
    velocities,
    v_immersed_bc,
    closure,
    diffusivities,
    clock,
    model_fields,
    buoyancy;
    location
)

Calculate the total viscous dissipation term as

VISC = ∂ⱼ τ₂ⱼ + ∂ⱼ τ₂ⱼ_immersed,

where τ₂ⱼ is the interior viscous stress tensor and τ₂ⱼ_immersed is the immersed boundary viscous stress tensor for the y-momentum equation.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> VISC = VMomentumEquation.TotalViscousDissipation(model)VTotalViscousDissipation KernelFunctionOperation at (Center, Face, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: total_∂ⱼ_τ₂ⱼ (generic function with 1 method)└── arguments: ("NamedTuple", "Nothing", "Nothing", "Nothing", "Clock", "NamedTuple", "Nothing")└── computes: total viscous term (interior + immersed, y)  ∂ⱼτ₂ⱼ
source
Oceanostics.VMomentumEquation.StokesShearType
StokesShear(model, stokes_drift, velocities, time; location)

Calculate the Stokes shear term as

STOKES_SHEAR = ((∇ × uˢ) × u)ᵧ

where uˢ is the Stokes drift velocity and u is the velocity vector.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> STOKES = VMomentumEquation.StokesShear(model)VStokesShear KernelFunctionOperation at (Center, Face, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: y_curl_Uˢ_cross_U (generic function with 3 methods)└── arguments: ("Nothing", "NamedTuple", "Float64")└── computes: Stokes shear forcing (y)  ((∇ × u⃗ˢ) × u⃗)_y
source
Oceanostics.VMomentumEquation.StokesTendencyType
StokesTendency(model, stokes_drift, time; location)

Calculate the Stokes tendency term as

STOKES_TEND = ∂vˢ/∂t

where vˢ is the y-component of the Stokes drift velocity.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> STOKES = VMomentumEquation.StokesTendency(model)VStokesTendency KernelFunctionOperation at (Center, Face, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: ∂t_vˢ (generic function with 4 methods)└── arguments: ("Nothing", "Float64")└── computes: Stokes drift tendency (y)  ∂vˢ/∂t
source
Oceanostics.VMomentumEquation.ForcingType
Forcing(model, forcing, clock, model_fields; location)

Calculate the forcing term Fᵛ on the y-momentum equation for model.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> FORC = VMomentumEquation.Forcing(model)VForcing KernelFunctionOperation at (Center, Face, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: forcing_cfc (generic function with 1 method)└── arguments: ("Returns", "Clock", "NamedTuple")└── computes: momentum forcing (y)  Fᵛ
source
Oceanostics.VMomentumEquation.TendencyType
Tendency(
    model,
    advection_scheme,
    coriolis,
    closure,
    v_immersed_bc,
    velocities,
    free_surface,
    tracers,
    buoyancy,
    closure_fields,
    hydrostatic_pressure_anomaly,
    auxiliary_fields,
    vertical_coordinate,
    clock,
    forcing_func;
    location
)

Calculate the v-momentum tendency Gᵛ as computed by Oceananigans, where Oceananigans writes the momentum equation as ∂_t v = Gᵛ - ∂ᵧ p_n. Gᵛ is the sum of every term in v_velocity_tendency (NH) or hydrostatic_free_surface_v_velocity_tendency (HFS): advection, buoyancy (NH only — absorbed into the hydrostatic pressure for HFS), Coriolis, hydrostatic pressure gradient, viscous and immersed-viscous stress divergence, Stokes shear and tendency (NH only), and forcing. The terms exposed by VMomentumEquation reconstruct Gᵛ to machine precision.

Gᵛ is not the full time derivative ∂_t v. The remaining piece -∂ᵧ p_n — the nonhydrostatic pressure gradient for NonhydrostaticModel, or the implicit barotropic free-surface correction for HydrostaticFreeSurfaceModel with ImplicitFreeSurface — is applied separately by the pressure-projection / barotropic-correction step of the time-stepper and is not part of Gᵛ or of any diagnostic in this module.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> TEND = VMomentumEquation.Tendency(model)VTendency KernelFunctionOperation at (Center, Face, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: v_velocity_tendency (generic function with 1 method)└── arguments: ("Centered", "Nothing", "Nothing", "Nothing", "Nothing", "Nothing", "Oceananigans.Models.NonhydrostaticModels.BackgroundFields", "NamedTuple", "NamedTuple", "NamedTuple", "Nothing", "Nothing", "Clock", "Returns")└── computes: total tendency of the v-momentum equation
source

W-momentum

Oceanostics.WMomentumEquation.AdvectionType
Advection(model, velocities, advection_scheme; location)

Calculate the advection of w-momentum as

ADV = ∂ⱼ (uⱼ w)

using Oceananigans' kernel div_𝐯w.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> ADV = WMomentumEquation.Advection(model)WAdvection KernelFunctionOperation at (Center, Center, Face)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: div_𝐯w (generic function with 4 methods)└── arguments: ("Centered", "NamedTuple", "Field")└── computes: advection of w-momentum  ∂ⱼ(uⱼw)
source
Oceanostics.WMomentumEquation.BuoyancyAccelerationType
BuoyancyAcceleration(model, buoyancy, tracers; location)

Calculate the buoyancy acceleration in the z-direction as

BUOY = ĝ_z b

where ĝ_z is the z-component of the gravitational unit vector and b is the buoyancy.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b);julia> BUOY = WMomentumEquation.BuoyancyAcceleration(model)WBuoyancyAcceleration KernelFunctionOperation at (Center, Center, Face)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: z_dot_g_bᶜᶜᶠ (generic function with 3 methods)└── arguments: ("BuoyancyForce", "NamedTuple")└── computes: buoyancy acceleration (z)  ĝ_z b
source
Oceanostics.WMomentumEquation.CoriolisAccelerationType
CoriolisAcceleration(model, coriolis, velocities; location)

Calculate the Coriolis acceleration in the z-direction as

COR = (f × u)_z

where f is the Coriolis parameter vector and u is the velocity vector.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid; coriolis=FPlane(f=1e-4));julia> COR = WMomentumEquation.CoriolisAcceleration(model)WCoriolisAcceleration KernelFunctionOperation at (Center, Center, Face)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: z_f_cross_U (generic function with 6 methods)└── arguments: ("FPlane", "NamedTuple")└── computes: Coriolis acceleration (z)  (f⃗ × u⃗)_z
source
Oceanostics.WMomentumEquation.ViscousDissipationType
ViscousDissipation(
    model,
    closure,
    diffusivities,
    clock,
    model_fields,
    buoyancy;
    location
)

Calculate the viscous dissipation term (excluding immersed boundaries) as

VISC = ∂ⱼ τ₃ⱼ,

where τ₃ⱼ is the viscous stress tensor for the z-momentum equation.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> VISC = WMomentumEquation.ViscousDissipation(model)WViscousDissipation KernelFunctionOperation at (Center, Center, Face)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: ∂ⱼ_τ₃ⱼ (generic function with 8 methods)└── arguments: ("Nothing", "Nothing", "Clock", "NamedTuple", "Nothing")└── computes: viscous term (interior, z)  ∂ⱼτ₃ⱼ
source
Oceanostics.WMomentumEquation.ImmersedViscousDissipationType
ImmersedViscousDissipation(
    model,
    velocities,
    w_immersed_bc,
    closure,
    diffusivities,
    clock,
    model_fields;
    location
)

Calculate the viscous dissipation term due to immersed boundaries as

VISC = ∂ⱼ τ₃ⱼ,

where τ₃ⱼ is the immersed boundary viscous stress tensor for the z-momentum equation.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> VISC = WMomentumEquation.ImmersedViscousDissipation(model)WImmersedViscousDissipation KernelFunctionOperation at (Center, Center, Face)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: immersed_∂ⱼ_τ₃ⱼ (generic function with 2 methods)└── arguments: ("NamedTuple", "Nothing", "Nothing", "Nothing", "Clock", "NamedTuple")└── computes: viscous term through immersed boundaries (z)  ∂ⱼτ₃ⱼ
source
Oceanostics.WMomentumEquation.TotalViscousDissipationType
TotalViscousDissipation(
    model,
    velocities,
    w_immersed_bc,
    closure,
    diffusivities,
    clock,
    model_fields,
    buoyancy;
    location
)

Calculate the total viscous dissipation term as

VISC = ∂ⱼ τ₃ⱼ + ∂ⱼ τ₃ⱼ_immersed,

where τ₃ⱼ is the interior viscous stress tensor and τ₃ⱼ_immersed is the immersed boundary viscous stress tensor for the z-momentum equation.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> VISC = WMomentumEquation.TotalViscousDissipation(model)WTotalViscousDissipation KernelFunctionOperation at (Center, Center, Face)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: total_∂ⱼ_τ₃ⱼ (generic function with 1 method)└── arguments: ("NamedTuple", "Nothing", "Nothing", "Nothing", "Clock", "NamedTuple", "Nothing")└── computes: total viscous term (interior + immersed, z)  ∂ⱼτ₃ⱼ
source
Oceanostics.WMomentumEquation.StokesShearType
StokesShear(model, stokes_drift, velocities, time; location)

Calculate the Stokes shear term as

STOKES_SHEAR = ((∇ × uˢ) × u)_z

where uˢ is the Stokes drift velocity and u is the velocity vector.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> STOKES = WMomentumEquation.StokesShear(model)WStokesShear KernelFunctionOperation at (Center, Center, Face)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: z_curl_Uˢ_cross_U (generic function with 3 methods)└── arguments: ("Nothing", "NamedTuple", "Float64")└── computes: Stokes shear forcing (z)  ((∇ × u⃗ˢ) × u⃗)_z
source
Oceanostics.WMomentumEquation.StokesTendencyType
StokesTendency(model, stokes_drift, time; location)

Calculate the Stokes tendency term as

STOKES_TEND = ∂wˢ/∂t

where wˢ is the z-component of the Stokes drift velocity.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> STOKES = WMomentumEquation.StokesTendency(model)WStokesTendency KernelFunctionOperation at (Center, Center, Face)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: ∂t_wˢ (generic function with 4 methods)└── arguments: ("Nothing", "Float64")└── computes: Stokes drift tendency (z)  ∂wˢ/∂t
source
Oceanostics.WMomentumEquation.ForcingType
Forcing(model, forcing, clock, model_fields; location)

Calculate the forcing term on the z-momentum equation for model.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> FORC = WMomentumEquation.Forcing(model)WForcing KernelFunctionOperation at (Center, Center, Face)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: forcing_ccf (generic function with 1 method)└── arguments: ("Returns", "Clock", "NamedTuple")└── computes: momentum forcing (z)  Fʷ
source
Oceanostics.WMomentumEquation.TendencyType
Tendency(
    model,
    advection_scheme,
    coriolis,
    stokes_drift,
    closure,
    w_immersed_bc,
    buoyancy,
    background_fields,
    velocities,
    tracers,
    auxiliary_fields,
    diffusivities,
    hydrostatic_pressure,
    clock,
    forcing_func;
    location
)

Calculate the w-momentum tendency as computed by Oceananigans, where Oceananigans writes the vertical momentum equation as ∂_t w = Gʷ - ∂_z p_NHS. is the sum of every term in w_velocity_tendency: advection, buoyancy (only when hydrostatic_pressure_anomaly = nothing; otherwise absorbed via the maybe_z_dot_g_bᶜᶜᶠ branch), Coriolis, viscous and immersed-viscous stress divergence, Stokes shear and tendency, and forcing. The terms exposed by WMomentumEquation reconstruct to machine precision in either pressure-splitting mode.

is not the full time derivative ∂_t w. The remaining piece -∂_z p_NHS (the nonhydrostatic pressure gradient) is applied separately by the pressure-projection step of the time-stepper and is not part of or of any diagnostic in this module.

HydrostaticFreeSurfaceModel does not have a prognostic w-momentum equation (w is diagnosed from continuity), so Tendency is only defined for NonhydrostaticModel.

julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> TEND = WMomentumEquation.Tendency(model)WTendency KernelFunctionOperation at (Center, Center, Face)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: w_velocity_tendency (generic function with 1 method)└── arguments: ("Centered", "Nothing", "Nothing", "Nothing", "Nothing", "Nothing", "Oceananigans.Models.NonhydrostaticModels.BackgroundFields", "NamedTuple", "NamedTuple", "NamedTuple", "Nothing", "Nothing", "Clock", "Returns")└── computes: total tendency of the w-momentum equation
source