Spatial Operators
This page covers the theory behind the higher-level spatial operators in finitevolX: divergence, vorticity (relative and potential), diffusion, Coriolis, momentum advection (vortex-force form), kinetic energy, and diagnostic quantities.
Overview
The spatial operators build on the primitive difference and interpolation operators to compute physically meaningful quantities on the Arakawa C-grid. The table below lists the main operators and what they compute:
| Operator | Module | Output location | Physical quantity |
|---|---|---|---|
Divergence2D |
divergence | T-points | \(\nabla \cdot (u, v)\) |
Vorticity2D.relative_vorticity |
vorticity | X-points | \(\zeta = \partial v/\partial x - \partial u/\partial y\) |
Vorticity2D.potential_vorticity |
vorticity | X-points | \(q = (\zeta + f)/h\) |
Vorticity2D.pv_flux_energy_conserving |
vorticity | U/V-points | \((q\,Uh, q\,Vh)\) |
Diffusion2D / diffusion_2d |
diffusion | T-points | \(\nabla \cdot (\kappa\,\nabla h)\) |
BiharmonicDiffusion2D |
diffusion | T-points | \(-\kappa\,\nabla^4 h\) |
Coriolis2D |
coriolis | U/V-points | \((+\overline{f_X\bar{v}}, -\overline{f_X\bar{u}})\), energy-conserving |
MomentumAdvection2D |
momentum | U/V-points | Vortex-force form |
Vorticity2D.enstrophy / .potential_enstrophy |
vorticity | X-points | \(\tfrac12\zeta^2\), \(\tfrac12 q^2 h\) |
Energetics2D |
diagnostics | T-points | KE, Bernoulli potential, APE |
Strain2D |
diagnostics | T/X-points | shear & tensor strain, \(\sigma^2\), Okubo–Weiss |
QGPotentialVorticity2D |
diagnostics | T-points | \(q = \nabla^2\psi/f_0 + \beta(y - y_0)/f_0\) (− stretching) |
ArakawaJacobian2D / arakawa_jacobian |
jacobian | T-points / interior | \(J(f, g)\) |
Divergence
Theory
The flux-form divergence of the horizontal velocity field at T-points:
In the shallow-water model the continuity equation is \(\partial h/\partial t + \nabla \cdot (h\,\mathbf{u}) = 0\), so computing the divergence of the mass flux \((hu, hv)\) gives the thickness tendency \(-\delta(hu, hv)\).
For a non-divergent velocity field (\(\delta = 0\)), the mass flux form and
the simpler advective form agree. On a C-grid, the discrete divergence is
exactly zero for any velocity field derived from a streamfunction via
grad_perp — this is the key property that makes QG models mass-conservative.
Usage
from finitevolx import CartesianGrid2D, Divergence2D, divergence_2d
grid = CartesianGrid2D.from_interior(64, 64, 1e6, 1e6)
# Class-based
div_op = Divergence2D(grid=grid)
delta = div_op(u, v) # T-points
# Functional (no grid object needed)
delta = divergence_2d(u, v, dx=grid.dx, dy=grid.dy)
Vorticity
Relative Vorticity
The relative vorticity at X-points (NE corners):
This is the discrete curl of the 2-D velocity field. It is placed at
X-points because this is the natural location for diff_x_V_to_X and
diff_y_U_to_X.
Potential Vorticity (Shallow Water)
The shallow-water potential vorticity at X-points:
where \(f\) and \(h\) are interpolated from T-points to X-points via T_to_X.
Thin-layer singularity
\(q\) is set to NaN where \(h = 0\) to avoid division by zero. Guard
against this with a minimum thickness clamp before calling the PV operator.
PV Fluxes
Vorticity2D offers two ways of combining \(q\) (at X-points) with a pair of
face fields \((F_u, F_v)\) at U/V-points:
| Method | Form | U-point value |
|---|---|---|
pv_flux_energy_conserving |
multiply at the corners, then average | \(\overline{q\,\overline{F_u}^{y}}^{y}\) |
pv_flux_enstrophy_conserving |
average \(q\) to the face, then multiply | \(\overline{q}^{y}\,F_u\) |
pv_flux_blend |
linear blend \(\alpha E + (1-\alpha) Z\) | — |
The names refer to the vortex-force use (Sadourny 1975): when the
cross-face transport is used — \(+q\,V\) in the \(u\)-equation, \(-q\,U\) in the
\(v\)-equation, as MomentumAdvection2D forms it — the corner-product form
conserves total energy and the average-then-multiply form conserves
potential enstrophy for non-divergent flow. Passed the same-component
fields, e.g. (q, u*h, v*h), the methods return the PV transport flux
\((q\,Uh, q\,Vh)\) and neither label applies.
The blend is not the Arakawa–Lamb (1981) scheme
pv_flux_blend is a linear blend of the two Sadourny forms and, for
\(0 < \alpha < 1\), conserves neither invariant exactly. The
Arakawa–Lamb (1981) scheme is a distinct eight-point stencil and is not
implemented. The old name pv_flux_arakawa_lamb (and
scheme="al" below) is a deprecated alias that warns (#283).
The corner products are formed on the full arrays, so the first interior row and column read the south / west ghost corners of the caller's ghost-filled \(q\), \(u\), \(v\): fill them (periodic, free-slip, …) before calling.
Usage
from finitevolx import CartesianGrid2D, Vorticity2D
grid = CartesianGrid2D.from_interior(64, 64, 1e6, 1e6)
vort = Vorticity2D(grid=grid)
# Relative vorticity at X-points
zeta = vort.relative_vorticity(u, v)
# Shallow-water PV at X-points
q = vort.potential_vorticity(u, v, h, f)
# PV transport flux (q*Uh, q*Vh) at U/V-points
qu, qv = vort.pv_flux_energy_conserving(q, u*h, v*h)
Diffusion
Harmonic Diffusion
The harmonic (Laplacian) diffusion computes the tracer tendency \(\partial h/\partial t = \nabla \cdot (\kappa\,\nabla h)\) in flux form using forward-then-backward differences:
-
East-face flux (forward T → U): $$ F_{j,i+\frac{1}{2}}^x = \kappa_{j,i+\frac{1}{2}}\,\frac{h_{j,i+1} - h_{j,i}}{\Delta x} $$
-
North-face flux (forward T → V): $$ F_{j+\frac{1}{2},i}^y = \kappa_{j+\frac{1}{2},i}\,\frac{h_{j+1,i} - h_{j,i}}{\Delta y} $$
-
Tendency (backward U/V → T): $$ \frac{\partial h}{\partial t}\bigg|{j,i} = \frac{F}{2}}^x - F_{j,i-\frac{1}{2}}^x}{\Delta x
- \frac{F_{j+\frac{1}{2},i}^y - F_{j-\frac{1}{2},i}^y}{\Delta y} $$
For uniform \(\kappa\) this reduces to \(\kappa\,\nabla^2 h\). For spatially
varying \(\kappa\) (e.g. near coastlines or mixed-layer parameterisations) the
flux form ensures discrete conservation of the diffused quantity. A
T-point \(\kappa\) field is averaged to each face,
\(\kappa_{j,i+\frac{1}{2}} = \tfrac{1}{2}(\kappa_{j,i} + \kappa_{j,i+1})\) (and
likewise in \(y\)), which keeps the operator symmetric and second-order
accurate; the spherical operators do the same. With wall="open" the wall
faces average the edge cell with the ghost cell of \(\kappa\), so fill the
\(\kappa\) halo consistently with the tracer BC (wrapped for periodic, a copy
of the edge cell for zero gradient).
Boundary conditions by default: Face fluxes at domain walls are zero (no-flux, closed-wall), giving a homogeneous Neumann BC \(\partial h/\partial n = 0\).
Biharmonic Diffusion
Biharmonic diffusion applies the harmonic operator twice:
It damps small scales much more strongly than large scales: for a wavenumber \(k\), harmonic damping scales as \(\kappa k^2\) while biharmonic scales as \(\kappa k^4\). This is used to remove grid-scale noise without excessively damping the resolved mesoscale.
Which diffusion to use
- Harmonic (
Diffusion2D): simple, physically motivated, appropriate for tracer diffusion when scale-selectivity is not critical. - Biharmonic (
BiharmonicDiffusion2D): preferred for velocity diffusion in eddy-resolving models because it acts mainly on unresolved scales. Requires a larger \(\kappa\) to achieve the same total dissipation.
Mask Support
All diffusion classes take an optional mask= (Mask2D / Mask3D) field.
Fluxes through faces with a dry side are zeroed and land-cell tendencies are
zeroed, giving natural no-flux BCs at irregular coastlines without any
additional code. The masked biharmonic operators apply this masked no-flux
Laplacian \(D_m\) twice, \(-\kappa\,D_m(D_m h)\), after zeroing the dry inputs:
the coast is no-flux for both \(h\) and \(\nabla^2 h\), land values (even NaN)
never enter the wet tendency, and the operator conserves \(\sum_\text{wet} h\).
Usage
from finitevolx import Diffusion2D, BiharmonicDiffusion2D, diffusion_2d
from finitevolx import CartesianGrid2D
grid = CartesianGrid2D.from_interior(64, 64, 1e6, 1e6)
# Class-based harmonic diffusion (unmasked)
diff_op = Diffusion2D(grid=grid)
dh_dt = diff_op(h, kappa=100.0)
# With masking — mask is a class field, not a per-call kwarg.
# Diffusion applies the intermediate flux-masking pattern internally
# (flux_x *= mask.u, flux_y *= mask.v, tendency *= mask.h).
diff_op_masked = Diffusion2D(grid=grid, mask=mask)
dh_dt = diff_op_masked(h, kappa=100.0)
# Biharmonic (scale-selective)
biharm_op = BiharmonicDiffusion2D(grid=grid)
dh_dt = biharm_op(h, kappa=1e9)
# Functional API (no grid object, mask-free)
dh_dt = diffusion_2d(h, kappa=100.0, dx=grid.dx, dy=grid.dy)
Coriolis
Theory
The Coriolis term is the vortex force with the planetary vorticity \(f\).
Coriolis2D uses the energy-conserving corner-product form of Sadourny
(1975), the same construction as the 'energy' vortex force of
MomentumAdvection2D: \(f\) (given at T-points) is averaged to the X-points,
multiplied there by the cross velocity averaged to the corner, and the
product is averaged back to the faces:
Both tendencies come from the same corner products, so the kinetic-energy
tendency \(\sum u\,\dot u + \sum v\,\dot v\) vanishes to round-off for any
\(f\), including a \(\beta\)-plane, on a periodic domain or a closed basin (with
the wall BCs of BoundaryConditionSet.closed(point="U") / "V"). For
uniform \(f\) it reduces to the familiar \(+f\,\overline{v}^{xy}\),
\(-f\,\overline{u}^{xy}\) (Interpolation2D.V_to_U / U_to_V). The form
used before #289, \(f_{\text{on-U}}\,\overline{v}_{\text{on-U}}\), conserves
energy only for uniform \(f\).
Usage
from finitevolx import CartesianGrid2D, Coriolis2D
grid = CartesianGrid2D.from_interior(64, 64, 1e6, 1e6)
cor = Coriolis2D(grid=grid)
# f at T-points
f = jnp.full((grid.Ny, grid.Nx), 1e-4)
du_cor, dv_cor = cor(u, v, f)
Momentum Advection (Vortex-Force Form)
Theory
The vector-invariant (vortex-force) form of the horizontal momentum advection separates the advection into a vorticity flux and a kinetic-energy gradient:
where \(\zeta = \partial v/\partial x - \partial u/\partial y\) is the relative vorticity at X-points, \(K = \tfrac{1}{2}(\bar{u}^2 + \bar{v}^2)\) is the kinetic energy at T-points, and \((\zeta\,v)_{\text{on-U}}\) is the vorticity-flux product interpolated to U-points.
This form has important discrete conservation properties depending on how the vorticity-flux product is evaluated. Three schemes are available:
| Scheme | Discrete form of \((\zeta\,v)_{\text{on-U}}\) | Key property | Reference |
|---|---|---|---|
'energy' (default) |
\(\overline{\zeta\,\overline{v}^{xy}}^{y}\): product at the corners, then average | Conserves total KE for non-divergent flow | Sadourny (1975) E-scheme |
'enstrophy' |
\(\overline{\zeta}^{y}\,\overline{v}^{xy}\): average to the face, then multiply | Conserves potential enstrophy | Sadourny (1975) Z-scheme |
'blend' |
\(\tfrac13 E + \tfrac23 Z\) | Conserves neither exactly | linear blend, not Arakawa & Lamb (1981); 'al' is a deprecated alias |
When to use each scheme
'energy': energy-conserving vortex force, the usual choice for long integrations.'enstrophy': preferred when enstrophy conservation matters (e.g. turbulence spindown experiments).'blend': a ⅓–⅔ blend of the two; it is not the Arakawa–Lamb (1981) scheme, which conserves both (not implemented, #283).
Usage
from finitevolx import CartesianGrid2D, MomentumAdvection2D
grid = CartesianGrid2D.from_interior(64, 64, 1e6, 1e6)
# Energy-conserving (default)
madv = MomentumAdvection2D(grid=grid)
du_adv, dv_adv = madv(u, v)
# Enstrophy-conserving
madv_z = MomentumAdvection2D(grid=grid)
du_adv, dv_adv = madv_z(u, v, scheme="enstrophy")
# 1/3 energy + 2/3 enstrophy blend
du_adv, dv_adv = madv(u, v, scheme="blend")
Arakawa Jacobian
Theory
The Arakawa (1966) Jacobian \(J(f, g) = \partial f/\partial x\,\partial g/\partial y - \partial f/\partial y\,\partial g/\partial x\) is the classical advection operator for QG models. Unlike simple centred differences, the three-term Arakawa average
conserves energy, enstrophy, and satisfies \(J(f, f) = 0\) and \(\int J(f, g)\,\mathrm{d}A = 0\) exactly at the discrete level.
The three terms are: - \(J^{++}\): standard centred (advective) form - \(J^{+\times}\): flux form in one direction - \(J^{\times+}\): flux form in the other direction
Conservation under a mask
The discrete identities above hold for the unmasked operator with
suitable boundary conditions (e.g. \(\psi\) constant on the boundary).
ArakawaJacobian2D(grid, mask=mask) zeroes land inputs and dry output
cells, so sums over a masked domain such as \(\sum J\) or
\(\sum q\,J\) are not guaranteed to vanish exactly.
Output shape
The class ArakawaJacobian2D returns the full [..., Ny, Nx] array with
a zero ghost ring (and dry cells zeroed when given a mask), ready for a
time-stepping loop. The functional arakawa_jacobian returns only the
interior [..., Ny−2, Nx−2]; the caller must embed it into the full
array.
Usage
from finitevolx import ArakawaJacobian2D, CartesianGrid2D, arakawa_jacobian
grid = CartesianGrid2D.from_interior(64, 64, 1e6, 1e6)
# QG vorticity advection: J(psi, q)
# psi, q have shape [Ny, Nx] with one ghost cell on each side
jac = ArakawaJacobian2D(grid=grid, mask=mask) # mask optional
Jpsi_q = jac(psi, q)
# Jpsi_q has shape [Ny, Nx]; ghost ring and dry cells are 0
# Functional form: interior only, no mask
Jpsi_q_int = arakawa_jacobian(psi, q, dx=grid.dx, dy=grid.dy)
# Jpsi_q_int has shape [Ny-2, Nx-2]
Decision Guide
What spatial operator do you need?
│
├── Divergence ∇·(u, v)
│ └── Divergence2D or divergence_2d
│
├── Relative vorticity ζ = ∂v/∂x − ∂u/∂y
│ ├── Class-based → Vorticity2D.relative_vorticity
│ └── Functional → Difference2D.curl
│
├── Shallow-water PV q = (ζ+f)/h
│ └── Vorticity2D.potential_vorticity
│
├── Tracer diffusion ∂q/∂t = ∇·(κ∇q)
│ ├── Harmonic → Diffusion2D or diffusion_2d
│ └── Biharmonic → BiharmonicDiffusion2D (scale-selective)
│
├── Coriolis tendency (f×v, −f×u)
│ └── Coriolis2D
│
├── Momentum advection ∂u/∂t|adv
│ ├── Conservative vortex-force form → MomentumAdvection2D
│ └── QG vorticity advection → ArakawaJacobian2D or arakawa_jacobian
│
└── Diagnostics (KE, APE, enstrophy, Okubo-Weiss, …)
├── Energy → Energetics2D
├── Strain / Okubo-Weiss → Strain2D
├── Enstrophy → Vorticity2D.enstrophy / .potential_enstrophy
├── QG PV → QGPotentialVorticity2D
└── See [Diagnostics API reference](api/diagnostics.md)
References
- Sadourny (1975) — C-grid momentum advection, energy- and enstrophy-conserving schemes
- Arakawa & Lamb (1981) — Combined energy/enstrophy conserving scheme
- Arakawa (1966) — Energy/enstrophy-conserving Jacobian operator
- Smagorinsky (1963) — Horizontal diffusion parameterisation
- Haidvogel & Beckmann (1999) — Momentum advection in ocean models