Advection on the C-grid: one term, three ways to compute it
The advection term moves a quantity along with the flow. In the continuous equations it can be written three ways that are all equal. On a grid they stop being equal: each one puts the approximation in a different place and keeps a different set of quantities exactly. This page goes from the PDE to the stencils, in finitevolX's conventions, and runs all three on the same test.
Same term, different approximation
- Flux form \(\nabla\cdot(\vect{v} q)\): move \(q\) to the faces where \(u, v\) live, multiply, take a difference. The difference is exact; all the error is in the face value of \(q\). Conserves \(\sum q\) for any face rule. This is where upwind and WENO live.
- Advective form \(\vect{v}\cdot\nabla q\): take the gradient of \(q\) at the centre and move \(u, v\) there. The error is in the gradient and in moving the velocity. Conserves nothing by construction. This is the slot a spectral gradient fits.
- Jacobian form \(J(\psi, q)\), for QG: combine \(\psi\) and \(q\) directly. Arakawa's average of three versions conserves \(\sum q\), \(\sum q^2\) and energy exactly, but adds no dissipation, so fronts ring unless viscosity is added.
- Centred schemes never damp; upwind-biased ones always do. On a sharp front that is the whole difference: the centred schemes keep the variance and ring, upwind smears, WENO5 does neither.
What the advection term says
Follow a parcel of water. If nothing creates or destroys \(q\) inside it, the parcel keeps its value, and the only reason \(q\) changes at a fixed point is that different water arrives. That is the material derivative being zero:
\[ \frac{\mathrm{D} q}{\mathrm{D} t} = \underbrace{\frac{\partial q}{\partial t}}_{\text{change at a fixed point}} + \underbrace{\vect{v}\cdot\nabla q}_{\text{advection}} = 0 . \]Here \(q(\vect{s}, t)\), with \(\vect{s} = (x, y)\), is any carried quantity (a tracer, layer thickness, potential vorticity), \(\vect{v} = (u, v)\) the velocity, and \(\nabla q := \nabla_{\vect{s}} q = (\partial_x q, \partial_y q)\). In one dimension with constant \(u\), the solution is the initial profile sliding along unchanged:
\[ \partial_t q + u \partial_x q = 0 \Longrightarrow q(x, t) = q_0(x - u t). \]So an exact scheme would move a sharp front without smearing it or adding wiggles. Every numerical scheme falls short of that in one of two ways: it smears (numerical diffusion) or it rings (dispersion, which makes over- and undershoots). The running example is the Gulf Stream PV front of the other pages: a jump in \(q\) a few grid cells wide, carried by eddies.
Where the term appears
| Model | Equation | Advected quantity | Velocity |
|---|---|---|---|
| Tracer | \(\partial_t c + \nabla\cdot(\vect{v} c) = \kappa\nabla^2 c\) | concentration \(c\) | given |
| Shallow water, mass | \(\partial_t h + \nabla\cdot(h\vect{v}) = 0\) | layer thickness \(h\) | prognostic |
| Shallow water, momentum | \(\partial_t\vect{v} + (\vect{v}\cdot\nabla)\vect{v} + f \hat{\vect{z}}\times\vect{v} = -g\nabla\eta\) | the velocity itself | prognostic |
| Navier–Stokes | \(\partial_t\vect{v} + (\vect{v}\cdot\nabla)\vect{v} = -\nabla p/\rho_0 + \nu\nabla^2\vect{v},\ \ \nabla\cdot\vect{v} = 0\) | the velocity itself | prognostic, divergence-free |
| QG | \(\partial_t q + J(\psi, q) = F + D\) | potential vorticity \(q\) | \(u = -\partial_y\psi,\ v = \partial_x\psi\), from \(q\) by inversion |
In QG, \(q = \nabla^2\psi - \psi/L_d^2 + \beta y\), so the velocity is computed from \(q\) itself through the Helmholtz inversion. That makes the problem nonlinear even though each step looks like advection by a given velocity. Momentum advection is treated at the end; the scalar case comes first because everything about it carries over.
Three ways to write the same term
The product rule gives the link between the advective and the flux form:
\[ \nabla\cdot(\vect{v} q) = \vect{v}\cdot\nabla q + q \nabla\cdot\vect{v} . \]When the flow is non-divergent (\(\nabla\cdot\vect{v} = 0\): incompressible flow, and every QG flow) the two are equal. With a streamfunction, \(u = -\partial_y\psi\) and \(v = \partial_x\psi\), the same term is a Jacobian, and the Jacobian itself has three equivalent forms:
\[ J(\psi, q) = \partial_x\psi \partial_y q - \partial_y\psi \partial_x q = \partial_x(\psi \partial_y q) - \partial_y(\psi \partial_x q) = \partial_y(q \partial_x\psi) - \partial_x(q \partial_y\psi) = \vect{v}\cdot\nabla q . \]In a closed or periodic domain with \(\nabla\cdot\vect{v} = 0\), the continuous term conserves three things, each by an integration by parts:
\[ \frac{\mathrm{d}}{\mathrm{d}t} \int q \mathrm{d}A = 0, \] \[ \frac{\mathrm{d}}{\mathrm{d}t} \int \frac{1}{2} q^2 \mathrm{d}A = 0, \] \[ \int \psi J(\psi, q) \mathrm{d}A = 0 . \]The first is mass (or total tracer), the second variance (enstrophy, for PV), and the third gives energy conservation in QG. On a grid the integrations by parts become sums, and they only work if the stencils are built for them. That is the real difference between the three forms.
The grid
finitevolX uses the Arakawa C-grid with same-index colocation: every array has the same shape [Ny, Nx], and the index [j, i] means a different physical point depending on which kind of array it is.
V[j,i] ──── X[j,i] ──── V[j,i+1]
│ │
│ T[j,i] │ T cell centre (j, i)
│ │ U east face (j, i+1/2)
V[j-1,i] ── X[j-1,i] ── V[j-1,i+1] V north face (j+1/2, i)
X NE corner (j+1/2, i+1/2)
U[j,i-1] U[j,i] U[j,i+1]
In one dimension, along a row:
T[i-1] T[i] T[i+1]
|------o------|------o------|------o------|
U[i-1] U[i]
(i-1/2) (i+1/2)
q lives at o (cell averages), u lives at | (faces)
| Point | Holds | Why there |
|---|---|---|
| T, centre | \(q\), \(h\), tracers, pressure | cell averages: what a finite-volume cell stores |
| U, east face | \(u\), x-fluxes | a flux is what crosses a face |
| V, north face | \(v\), y-fluxes | same, in y |
| X, corner | \(\psi\), relative vorticity | differences of \(\psi\) across a face give the velocity on that face |
Arrays carry a one-cell ghost ring: the interior is [1:-1, 1:-1]. A neighbour is a shifted slice of the same shape as the interior: q[1:-1, 2:] is q[j, i+1], q[2:, 1:-1] is q[j+1, i]. finitevolX's raw stencils (diff_x_fwd, avg_y_bwd, …) return the interior, and interior(values, like) pads the result back to [Ny, Nx]. Every stencil line carries its half-index formula as a comment, and every array its shape and point type.
Putting \(\psi\) at the corners gives velocities on the faces that are divergence-free to round-off, because the four corner values cancel in pairs:
import finitevolx as fvx import jax.numpy as jnp # Shapes: Ny, Nx = 130 × 130 total, 128 × 128 interior plus a one-cell ghost ring. # Same-index colocation: # T[j, i] → cell centre (j, i) U[j, i] → east face (j, i+1/2) # V[j, i] → north face (j+1/2, i) X[j, i] → NE corner (j+1/2, i+1/2) # today (finitevolx) grid = fvx.CartesianGrid2D.from_interior(128, 128, 1.024e6, 1.024e6) dx, dy = grid.dx, grid.dy # q: [Ny, Nx] - advected quantity (PV, a tracer) at T-points # psi: [Ny, Nx] - streamfunction at X-points (corners) # u: [Ny, Nx] - x-velocity at U-points (east faces) # v: [Ny, Nx] - y-velocity at V-points (north faces) # today (finitevolx) # u[j, i+1/2] = -(psi[j+1/2, i+1/2] - psi[j-1/2, i+1/2]) / dy u = fvx.interior(-fvx.diff_y_bwd(psi) / dy, psi) # [Ny, Nx] U-points # today (finitevolx) # v[j+1/2, i] = (psi[j+1/2, i+1/2] - psi[j+1/2, i-1/2]) / dx v = fvx.interior(fvx.diff_x_bwd(psi) / dx, psi) # [Ny, Nx] V-points # today (finitevolx) # div[j, i] = (u[j, i+1/2] - u[j, i-1/2]) / dx + (v[j+1/2, i] - v[j-1/2, i]) / dy # zero to round-off: the corner values cancel pairwise div = fvx.interior(fvx.diff_x_bwd(u) / dx + fvx.diff_y_bwd(v) / dy, u) # [Ny, Nx] T
Markers: today (finitevolx) names ship on finitevolX main (checked 24 Oct 2026, 831c69f). Shape column: total shape and point type, as in finitevolX.
Form 1 · Flux form (finite volume)
Integrate \(\partial_t q + \nabla\cdot(\vect{v} q) = 0\) over cell \((j, i)\) and apply the divergence theorem. The cell average changes only through what crosses its four faces:
\[ \frac{\mathrm{d}q_{j,i}}{\mathrm{d}t} = -\frac{F^x_{j, i+\frac12} - F^x_{j, i-\frac12}}{\Delta x} - \frac{F^y_{j+\frac12, i} - F^y_{j-\frac12, i}}{\Delta y}, \] \[ F^x_{j, i+\frac12} = u_{j, i+\frac12} \hat q_{j, i+\frac12} . \]This is where the C-grid earns its keep: \(u\) is already on the face. The one missing piece is \(\hat q_{j, i+\frac12}\), the value of \(q\) on the face, which has to be estimated from the cell averages on either side. The final difference is exact. So in flux form, choosing a scheme means choosing how to estimate \(\hat q\), and nothing else:
| Face rule | Face value for u ≥ 0 | Effect |
|---|---|---|
| centred | \(\frac{1}{2}(q_i + q_{i+1})\) | second order, no damping, rings at fronts |
| upwind 1 | \(q_i\) | first order, strong damping, never rings |
| WENO5 | a smoothness-weighted blend of three 3-cell fits | fifth order where smooth, upwind-like at fronts |
# --- 1. face values: the only approximation in flux form --- # today (finitevolx) # centred: q[j, i+1/2] = (q[j, i] + q[j, i+1]) / 2 q_east = fvx.interior(fvx.avg_x_fwd(q), q) # [Ny, Nx] U-points # today (finitevolx) # centred: q[j+1/2, i] = (q[j, i] + q[j+1, i]) / 2 q_north = fvx.interior(fvx.avg_y_fwd(q), q) # [Ny, Nx] V-points # upwind1 instead: take the cell the flow comes from # q[j, i+1/2] = q[j, i] if u[j, i+1/2] >= 0 else q[j, i+1] q_here = q[..., 1:-1, 1:-1] # [Ny-2, Nx-2] T (j, i) q_east_cell = q[..., 1:-1, 2:] # [Ny-2, Nx-2] T (j, i+1) u_face = u[..., 1:-1, 1:-1] # [Ny-2, Nx-2] U (j, i+1/2) q_up = jnp.where(u_face >= 0.0, q_here, q_east_cell) # [Ny-2, Nx-2] U q_east_up = fvx.interior(q_up, q) # [Ny, Nx] U-points # --- 2. fluxes at the faces --- # fe[j, i+1/2] = u[j, i+1/2] * q[j, i+1/2] flux_east = u * q_east # [Ny, Nx] U-points # fn[j+1/2, i] = v[j+1/2, i] * q[j+1/2, i] flux_north = v * q_north # [Ny, Nx] V-points # --- 3. divergence: exact, telescoping --- # today (finitevolx) # dfe[j, i] = (fe[j, i+1/2] - fe[j, i-1/2]) / dx div_x = fvx.diff_x_bwd(flux_east) / dx # [Ny-2, Nx-2] T # today (finitevolx) # dfn[j, i] = (fn[j+1/2, i] - fn[j-1/2, i]) / dy div_y = fvx.diff_y_bwd(flux_north) / dy # [Ny-2, Nx-2] T # dq[j, i] = -(dfe[j, i] + dfn[j, i]) tendency = fvx.interior(-(div_x + div_y), q) # [Ny, Nx] T-points # the same pipeline, with WENO5 faces, masks and closed walls: # today (finitevolx) tendency_weno = fvx.Advection2D(grid)(q, u, v, method="weno5") # [Ny, Nx] T
What it conserves
Mass, for any face rule. The east flux of cell \(i\) is the west flux of cell \(i+1\), so summing the tendency over all cells, every interior flux cancels (it telescopes):
\[ \sum_i \frac{\mathrm{d}q_i}{\mathrm{d}t} \Delta x = -\sum_i(F_{i+\frac12} - F_{i-\frac12}) = -(F_{\text{east wall}} - F_{\text{west wall}}) = 0 . \]Variance, for the centred rule only, when the discrete velocity is divergence-free. Summing by parts and using \(F_{i+\frac12} = u_{i+\frac12}\frac{1}{2}(q_i + q_{i+1})\):
\[ \frac{\mathrm{d}}{\mathrm{d}t}\sum_i\frac{1}{2} q_i^2 = \sum_i F_{i+\frac12} (q_{i+1} - q_i) = \sum_i \frac{1}{2} u_{i+\frac12}(q_{i+1}^2 - q_i^2) = -\frac{1}{2}\sum_i q_i^2 (u_{i+\frac12} - u_{i-\frac12}) = 0 . \](per unit \(\Delta x\); in 2-D the last bracket becomes the discrete divergence, which is zero by construction.) Any upwind bias adds the term \(-\vert u \vert \kappa (q_{i+1}-q_i)\) to the flux, which removes variance: that is the dissipation discussed on the learned WENO page.
Form 2 · Advective form
Compute the gradient of \(q\) at the cell centre and multiply by the velocity there:
\[ \frac{\mathrm{d}q_{j,i}}{\mathrm{d}t} = -\left( \bar u_{j,i} \frac{q_{j,i+1} - q_{j,i-1}}{2\Delta x} + \bar v_{j,i} \frac{q_{j+1,i} - q_{j-1,i}}{2\Delta y} \right), \] \[ \bar u_{j,i} = \frac{1}{2}(u_{j, i-\frac12} + u_{j, i+\frac12}). \]Now \(q\) stays where it is and the velocity is moved, from the faces to the centre. The gradient itself is a choice too: here the centred difference, which in finitevolX is a forward difference to the faces followed by an average back to the centre.
# --- 1. gradient of q: forward to the faces, then back to the centre --- # today (finitevolx) # dq_dx[j, i+1/2] = (q[j, i+1] - q[j, i]) / dx dqdx_east = fvx.interior(fvx.diff_x_fwd(q) / dx, q) # [Ny, Nx] U-points # today (finitevolx) # dq_dx[j, i] = (dq_dx[j, i-1/2] + dq_dx[j, i+1/2]) / 2 # = (q[j, i+1] - q[j, i-1]) / (2 dx) dqdx = fvx.interior(fvx.avg_x_bwd(dqdx_east), q) # [Ny, Nx] T-points # today (finitevolx) # dq_dy[j+1/2, i] = (q[j+1, i] - q[j, i]) / dy dqdy_north = fvx.interior(fvx.diff_y_fwd(q) / dy, q) # [Ny, Nx] V-points # today (finitevolx) # dq_dy[j, i] = (q[j+1, i] - q[j-1, i]) / (2 dy) dqdy = fvx.interior(fvx.avg_y_bwd(dqdy_north), q) # [Ny, Nx] T-points # --- 2. velocity: move u and v from the faces to the centre --- # today (finitevolx) # u[j, i] = (u[j, i-1/2] + u[j, i+1/2]) / 2 u_center = fvx.interior(fvx.avg_x_bwd(u), u) # [Ny, Nx] T-points # today (finitevolx) # v[j, i] = (v[j-1/2, i] + v[j+1/2, i]) / 2 v_center = fvx.interior(fvx.avg_y_bwd(v), v) # [Ny, Nx] T-points # --- 3. product at the centre --- # u ∂ₓq advect_x = u_center * dqdx # [Ny, Nx] T-points # v ∂ᵧq advect_y = v_center * dqdy # [Ny, Nx] T-points # dq[j, i] = -(u ∂ₓq + v ∂ᵧq)[j, i] tendency = -(advect_x + advect_y) # [Ny, Nx] T-points
What it conserves
Nothing, by construction. \(\sum_i \bar u_i (q_{i+1} - q_{i-1})\) does not telescope, because \(\bar u\) changes from cell to cell, so neither \(\sum q\) nor \(\sum q^2\) is guaranteed. For special flows the errors cancel; for a general flow they do not (see the table below).
The advective form is still the natural one in two places. In semi-Lagrangian schemes, which trace parcels back along \(\vect{v}\). And in spectral methods, where the gradient is computed in Fourier space: that is the slot a spectral derivative fills, and why a "perfect gradient" does not help the flux form, whose difficulty is the face value, not the gradient. For momentum in Navier–Stokes codes, a common compromise is the skew-symmetric form, half flux plus half advective, which conserves kinetic energy with centred stencils.
Form 3 · Jacobian form (Arakawa)
QG writes the term with \(\psi\) and \(q\) directly. Each of the three continuous forms of \(J\) gives a centred stencil (\(J^{++}\), \(J^{+\times}\), \(J^{\times+}\)); each alone conserves only some of the invariants. Arakawa (1966) showed that their average conserves all of them:
\[ J_{\text{A}} = \frac{1}{3} (J^{++} + J^{+\times} + J^{\times+}), \] \[ \sum J_{\text{A}} = 0, \] \[ \sum q J_{\text{A}}(\psi, q) = 0, \] \[ \sum \psi J_{\text{A}}(\psi, q) = 0, \] \[ J_{\text{A}}(\psi, \psi) = 0 . \]The three forms, as finitevolX writes them (on a collocated grid, f is \(\psi\) and g is \(q\)):
Jpp[j,i] = ( (f[j,i+1] - f[j,i-1]) * (g[j+1,i] - g[j-1,i])
- (f[j+1,i] - f[j-1,i]) * (g[j,i+1] - g[j,i-1]) ) / (4 dx dy)
Jpx[j,i] = ( f[j,i+1] * (g[j+1,i+1] - g[j-1,i+1])
- f[j,i-1] * (g[j+1,i-1] - g[j-1,i-1])
- f[j+1,i] * (g[j+1,i+1] - g[j+1,i-1])
+ f[j-1,i] * (g[j-1,i+1] - g[j-1,i-1]) ) / (4 dx dy)
Jxp[j,i] = ( g[j+1,i] * (f[j+1,i+1] - f[j+1,i-1])
- g[j-1,i] * (f[j-1,i+1] - f[j-1,i-1])
- g[j,i+1] * (f[j+1,i+1] - f[j-1,i+1])
+ g[j,i-1] * (f[j+1,i-1] - f[j-1,i-1]) ) / (4 dx dy)
\(J^{++}\) discretises \(\partial_x\psi \partial_y q - \partial_y\psi \partial_x q\); \(J^{+\times}\) discretises \(\partial_x(\psi \partial_y q) - \partial_y(\psi \partial_x q)\); \(J^{\times+}\) discretises \(\partial_y(q \partial_x\psi) - \partial_x(q \partial_y\psi)\). In code, \(\psi\) is first averaged from the corners to the centres so it sits with \(q\):
# psi_t: [Ny, Nx] - streamfunction collocated with q at T-points # today (finitevolx) # psi[j, i] = (psi[j-1/2, i-1/2] + psi[j-1/2, i+1/2] # + psi[j+1/2, i-1/2] + psi[j+1/2, i+1/2]) / 4 psi_t = fvx.interior(fvx.avg_xy_bwd(psi), psi) # [Ny, Nx] T-points # --- one of the three forms, written out: J⁺⁺, the plain centred one --- psi_e, psi_w = psi_t[..., 1:-1, 2:], psi_t[..., 1:-1, :-2] # [Ny-2, Nx-2] (j, i±1) psi_n, psi_s = psi_t[..., 2:, 1:-1], psi_t[..., :-2, 1:-1] # [Ny-2, Nx-2] (j±1, i) q_e, q_w = q[..., 1:-1, 2:], q[..., 1:-1, :-2] # [Ny-2, Nx-2] (j, i±1) q_n, q_s = q[..., 2:, 1:-1], q[..., :-2, 1:-1] # [Ny-2, Nx-2] (j±1, i) # Jpp[j, i] = ((psi[j, i+1] - psi[j, i-1]) * (q[j+1, i] - q[j-1, i]) # - (psi[j+1, i] - psi[j-1, i]) * (q[j, i+1] - q[j, i-1])) / (4 dx dy) dpsi_x = psi_e - psi_w # [Ny-2, Nx-2] dq_y = q_n - q_s # [Ny-2, Nx-2] dpsi_y = psi_n - psi_s # [Ny-2, Nx-2] dq_x = q_e - q_w # [Ny-2, Nx-2] jacobian_pp = (dpsi_x * dq_y - dpsi_y * dq_x) / (4.0 * dx * dy) # [Ny-2, Nx-2] T # --- the full Arakawa average J = (J⁺⁺ + J⁺ˣ + Jˣ⁺) / 3 --- # today (finitevolx) jacobian = fvx.arakawa_jacobian(psi_t, q, dx, dy) # [Ny, Nx] → [Ny-2, Nx-2] T # dq[j, i] = -J(psi, q)[j, i] tendency = fvx.interior(-jacobian, q) # [Ny, Nx] T-points
The price of exact conservation is zero numerical dissipation. Variance that cascades to the grid scale has nowhere to go, so a QG model with Arakawa advection needs explicit dissipation, the \(\nu\nabla^6\psi\) term of the other pages, or a learned closure in its place.
All three, on the same test
A doubly periodic box of 128 × 128 cells, a steady cellular flow \(\psi = \sin x \sin y\) placed at the corners (so the discrete divergence is \(4\times10^{-15}\)), and two initial fields centred just off the middle of one flow cell: a smooth Gaussian and a sharp disc. Each scheme runs one revolution of the cell centre with classical RK4 at Courant number 0.25.
Conservation, from the stencils alone
First the semi-discrete check: the instantaneous tendency of each invariant, for a general divergence-free flow (16 random Fourier modes in \(\psi\)) and a smooth \(q\). Zero to round-off means conserved by the stencils.
| Scheme | d/dt Σq | d/dt Σq² | Σψ · dq/dt (energy) |
|---|---|---|---|
| flux, centred | ✓ 0 | ✓ 0 | ✗ −3.3 × 10⁻² |
| flux, upwind 1 | ✓ 0 | −8.2 (damps) | ✗ −1.1 × 10⁻¹ |
| flux, WENO5 | ✓ 0 | −3.7 × 10⁻³ (damps) | ✗ −1.2 × 10⁻² |
| advective, centred | ✗ −2.7 × 10⁻⁴ | ✗ +9.0 × 10⁻³ | ✗ −1.7 × 10⁻² |
| Jacobian, Arakawa | ✓ 0 | ✓ 0 | ✓ 0 |
"0" is below \(10^{-15}\). The energy column uses \(\psi\) averaged to the cell centres. The advective form gains variance here: with nothing to stop it, that is how a non-conservative scheme can drift or blow up over long runs.
After one revolution
| Scheme | Disc: min | Disc: max | Disc: change in Σq² | Gaussian: change in Σq² |
|---|---|---|---|---|
| exact | 0.00 | 1.00 | 0 | 0 |
| flux, centred | −0.48 | 1.53 | −2 × 10⁻⁶ | −1 × 10⁻⁸ |
| flux, upwind 1 | 0.00 | 0.67 | −66 % | −70 % |
| flux, WENO5 | 0.00 | 1.00 | −20.5 % | −1.8 % |
| advective, centred | −0.54 | 1.54 | −0.8 % | +0.7 % |
| Jacobian, Arakawa | −0.61 | 1.54 | −6 × 10⁻⁷ | −1 × 10⁻⁸ |
The tiny changes for the conservative centred schemes are RK4's time-stepping error, not the stencils. Read the two tables together: the centred schemes keep the variance and ring; upwind 1 never rings and loses two thirds of the variance; WENO5 loses variance only near the front (1.8 % for the smooth Gaussian, 20.5 % for the disc) and stays in bounds. A front that is only a few cells wide cannot be both kept sharp and kept smooth by a linear scheme: that is Godunov's theorem, from the learned WENO page, in a picture.
Momentum advection, briefly
In shallow water and Navier–Stokes the velocity advects itself. On the C-grid the usual route is the vector-invariant form, which splits \((\vect{v}\cdot\nabla)\vect{v}\) into a rotational and a gradient part:
\[ (\vect{v}\cdot\nabla)\vect{v} = \zeta \hat{\vect{z}}\times\vect{v} + \nabla \left( \frac{1}{2} \Vert \vect{v} \Vert^2 \right), \] \[ \zeta = \partial_x v - \partial_y u . \]| Piece | Lives at | Needs moving |
|---|---|---|
| \(\zeta\) (or PV \(\zeta/h\)) | X, corners | to U and V, to multiply \(v\) and \(u\) |
| \(\frac{1}{2} \Vert \vect{v} \Vert^2\) | T, centres (average of \(u^2\), \(v^2\)) | no; its gradient lands on U and V |
| \(\zeta \hat{\vect{z}}\times\vect{v} = (-\zeta v,\ \zeta u)\) | U and V | \(v\) to U-points and \(u\) to V-points |
The same question returns: how to average \(\zeta\) and the velocities onto the faces. Different averages conserve energy, enstrophy, or both (Sadourny's and Arakawa–Lamb's schemes); finitevolX exposes them as MomentumAdvection2D and the PV-flux helpers such as pv_flux_arakawa_lamb. The structure of the choice is the same as for a scalar.
Back to the question
| Flux | Advective | Jacobian (Arakawa) | |
|---|---|---|---|
| What is moved between points | \(q\), to the faces | \(u, v\), to the centres | \(\psi\) to the centres; then averaged combinations |
| Where the error lives | the face value \(\hat q\) only | the gradient and \(\bar u\) | the centred differences |
| \(\sum q\) | ✓ any face rule | ✗ | ✓ |
| \(\sum q^2\) | ✓ centred; damped if upwind-biased | ✗ | ✓ |
| Energy (QG) | ✗ | ✗ | ✓ |
| Dissipation | chosen by the face rule | none (centred) | none; add it explicitly |
| Fronts | handled by upwind / WENO faces | ring | ring |
| Natural for | finite volume, tracers, shallow-water mass | spectral, semi-Lagrangian | QG with explicit closure |
So: yes, in flux form the face value of \(q\) is where everything happens, and that is exactly the decision WENO makes. The staggering does not cause it (finite volume needs face values on any grid); it only means \(u\) is already where the flux needs it. A spectral gradient would improve the advective form's gradient, but not the flux form's face value, and neither form gets dissipation for free.