M4E 2026 · Politecnico di Milano · 9–11 November

Polytopal Discontinuous Galerkin Methods for Geophysical Fluid Dynamics

The rotating shallow water equations on polygonal grids

Pedro S. Peixoto

with Felipe A. V. B. Alves, Ivan Fumagalli and Paola F. Antonietti

Instituto de Matemática e Estatística, Universidade de São Paulo

Visiting professor at MOX, Politecnico di Milano

The global atmosphere

Six days of an operational forecast, wind at 10 m: cyclones, jets and trade winds, a thin layer of air on a rotating sphere

Six days of wind at 10 m · ▶ Video: www.ime.usp.br/~pedrosp/talks/2026-polimi/figures/video-weather.mp4

ECMWF open data (CC BY 4.0), IFS forecast of 29 Oct 2026, 00 UTC

(0.5 min) What global models must simulate: the operational ECMWF forecast of 29 October 2026, six days of surface wind (colour: speed; particles: direction, sped up so the flow reads). The mid-latitude cyclones travelling eastwards in both hemispheres, the strong Southern Ocean westerlies, the trade winds, a typhoon in the western Pacific. All of it in a layer of air a few km thick on a rotating sphere: the reason why the shallow water equations are our testbed later. Next: on which grids we solve them.

Global Models Grids

\(\newcommand{\jump}[1]{[\![\,#1\,]\!]}\newcommand{\avg}[1]{\{\!\!\{\,#1\,\}\!\!\}}\)

• Lat–lon: tiny polar cells limit explicit time steps
• Semi-implicit, semi-Lagrangian: long steps, but global solves (e.g. ECMWF IFS model, UKMetOffice)

• Quasi-uniform cells: explicit and local, scales on many cores
• Icosahedral/Voronoi: MPAS, ICON, NICAM, DYNAMICO, GRIST, MONAN

Explicit: local, but the step is set by the smallest cell · Implicit: long steps, global communication

Williamson, J. Meteor. Soc. Japan (2007); Staniforth & Thuburn, Q. J. R. Meteorol. Soc. (2012)

(1 min) Global models solve the equations on the sphere, and the grid and the time stepping go together. On a lat-lon grid the meridians converge: the polar cells are tiny, and an explicit scheme must take tiny time steps (CFL). The classical way out is implicit or semi-implicit time stepping with semi-Lagrangian advection, as in ECMWF’s IFS (spectral transform): long steps, but each step needs global operations (transforms, elliptic solves), so communication grows with the number of cores. Quasi-uniform grids remove the small cells, so explicit, local schemes (in the horizontal: MPAS, ICON and others are horizontally explicit, vertically implicit) become efficient and scale. Icosahedral grids come as triangles (ICON) or as their dual, hexagons (MPAS, NICAM): the third and fourth globes are duals of each other. They also allow smooth local refinement (last globe: refined over South America, as in MONAN).

MPAS and MONAN

C-grid finite volumes (TRSK) on Voronoi cells

Model for Ocean-laNd-Atmosphere predictioN

MONAN, Brazil’s community model built on MPAS, operational at INPE since 3 September 2026: global ≈ 10 km, 11-day forecasts twice a day; next, 5 km over South America

Skamarock et al., Mon. Wea. Rev. (2012); Ringler et al., J. Comput. Phys. (2010); monanadmin.github.io (Sep. 2026)

(1 min) MPAS (NCAR/LANL): nonhydrostatic model on spherical centroidal Voronoi grids with C-grid staggering and the TRSK finite-volume scheme. MONAN, the Brazilian community model led by INPE, uses the MPAS dynamical core, with grids that can be refined over South America (globe in the middle, a JIGSAW mesh like those of MPAS). News: MONAN went operational at INPE on 3 September 2026, global at about 10 km, forecasts up to 11 days, twice a day (00 and 12 UTC), on the new Jaci supercomputer, evaluated better than INPE’s previous model (BAM) and competitive with international centres; next, a 5 km regional configuration over South America. Right: a MONAN regional precipitation forecast.

The catch: low order on irregular cells

Divergence of a solid-body rotation (exactly zero), icosahedral grid, 40 962 cells (≈ 120 km)

TRSK (MPAS) truncation error, momentum equation: no convergence on irregular grids

• Errors follow the grid, not the flow: grid imprinting
• Higher-order finite volumes need wide stencils: less locality

Inconsistency, but with enough accuracy to allow MONAN to be operational!?!

Peixoto & Barros, J. Comput. Phys. (2013); Peixoto, J. Comput. Phys. (2016); Weller, Thuburn & Cotter, Mon. Wea. Rev. (2012)

(1 min) The low-order finite volumes in use are sensitive to the cell geometry. Left: the discrete divergence of a solid-body rotation should vanish, but its error follows the icosahedral grid lines (grid imprinting; glevel 6, 40 962 cells, about 120 km): the strongest lines are the edges of the icosahedron the grid is built from (top left). Right: the TRSK momentum operator is inconsistent on irregular grids, its error does not converge. Raising the order of finite volumes on unstructured grids requires wide reconstruction stencils.

Goal: high order, local in space

High-order discretizations on polygonal grids, robust to cell irregularity,
that keep the data locality of finite volumes

(1 min) What we want: high order, robustness to the cell shape, and locality. A high-order finite-volume reconstruction reaches several rings of cells (left), so the halos and the communication grow with the order. Discontinuous Galerkin keeps a polynomial of degree p inside each cell and couples cells only through face fluxes (middle), so a parallel rank needs a single layer of ghost cells and one exchange per Runge-Kutta stage (right), whatever the order.

Polytopal discontinuous Galerkin (PolyDG)

• DG on arbitrary polygons and polyhedra
• \(hp\)-theory on general meshes, faces of any size
• Exact, quadrature-free integration on polytopes
• Seismic waves, brain, poroelasticity, …
• Not yet for geophysical flows: rotation, balance, vortices

This work: PolyDG for the rotating shallow water equations

open-source C++ PolyDG library

Cangiani, Georgoulis & Houston (2014); Antonietti et al. (2016, 2018, 2021, 2022). Meshes: Corti, Bonizzoni & Antonietti, J. Sci. Comput. (2024)

(0.6 min) Polytopal DG, a strength of the MOX group: DG on general polygons and polyhedra, with a mature hp-convergence theory and quadrature-free integration. It has been used for seismic waves, brain, poroelasticity, but not for geophysical flows, whose discretizations must also respect rotation, geostrophic balance and vortical dynamics. To our knowledge this is the first high-order DG for the rotating shallow water equations on general polygonal grids. Implemented in Vulpes.

Rotating shallow water equations

Vector-invariant form:

\[ \begin{alignedat}{2} \partial_t \eta + \nabla\cdot(\eta\mathbf{u}) &= 0 &&\qquad {\small\color{gray}\text{(continuity)}}\\ \partial_t \mathbf{u} + (\zeta+f)\,\mathbf{u}^{\perp} + \nabla B &= 0 &&\qquad {\small\color{gray}\text{(momentum)}} \end{alignedat} \]

\(\zeta = \partial_x v - \partial_y u\), \(\quad B = g(\eta+b) + \tfrac12|\mathbf{u}|^2\)

• Mass, energy, potential vorticity \(q = (\zeta+f)/\eta\)
• Linear regime:
  • steady geostrophic mode, \(f\mathbf{u}^{\perp} = -g\nabla\eta\)
  • inertia–gravity waves, \(\omega^2 = f^2 + gH|\mathbf{k}|^2\)
• \(f\)-plane, bi-periodic, Earth-sized: no sphere yet

The standard first step for the dynamical cores of global weather and climate models

Williamson et al., J. Comput. Phys. (1992); Vallis, Atmospheric and Oceanic Fluid Dynamics, Cambridge (2017); Zeitlin, Geophysical Fluid Dynamics, Oxford (2018)

(1 min) The standard testbed of horizontal dynamical cores: new schemes for global weather and climate models are first tried on shallow water (Williamson et al. 1992 test set). Vector-invariant form: vorticity and the Bernoulli function appear explicitly; their structure is what we want to keep. The solution carries fast inertia-gravity waves and slow, nearly geostrophic flow. To assess the discretization alone we work on the f-plane, on a bi-periodic, Earth-sized domain, with Voronoi meshes that mimic those of spherical grids. Books: Vallis (2017), Zeitlin (2018).

Discrete space on polygons

Basis function \(L_2(x)\,L_1(y)\) of the bounding box \(\mathcal{B}_\kappa\)

Example of a smooth field (stationary vortex), cell by cell

• \(V_h^p\): degree \(p\) per polygon, discontinuous
• Modal: Legendre polynomials of the bounding box

• Exact integration: quadrature-free or sub-triangles
• Block-diagonal mass matrix, explicit RK4

(1.2 min) Unknowns: polynomials of total degree p on each cell, no continuity between cells. Polygons have no natural nodes (no Gauss-Lobatto points as on quads), so the basis is modal: the unknowns are the coefficients of tensor Legendre polynomials of the cell’s bounding box, which keeps the local matrices well conditioned whatever the cell shape. Integrals are exact: quadrature-free reduction to the boundary (Euler’s identity) for the mass matrices, sub-triangulation for the nonlinear residual. The mass matrix is block-diagonal, so explicit RK4 needs no global solve. In 3D: the relative vorticity of the stationary vortex of the paper projected on P1 and P3; the jumps across faces (red) are what the numerical fluxes act on, and they shrink as p grows.

The weak form

On each polygon \(\kappa\), test with \(\phi\in\mathbb{P}^p(\kappa)\) and integrate by parts:

\[ \begin{aligned} (\partial_t\eta_h,\phi)_\kappa - (\eta_h\mathbf{u}_h,\nabla\phi)_\kappa + \color{#2e7d32}{\langle\,\widehat{\eta\mathbf{u}}\cdot\mathbf{n}_\kappa},\,\phi\,\rangle_{\partial\kappa} &= 0\\[6pt] (\partial_t\mathbf{u}_h,\phi)_\kappa + \big((\zeta_h+f)\,\mathbf{u}_h^{\perp},\phi\big)_\kappa - (B_h,\nabla\phi)_\kappa + \color{#2e7d32}{\langle\,\widehat{B}\,\mathbf{n}_\kappa},\,\phi\,\rangle_{\partial\kappa} &= 0 \end{aligned} \]

• Volume terms: local to \(\kappa\)
• Momentum: only \(\nabla B\) is integrated by parts
• On a face, two traces: \(\avg{a} = \tfrac12(a^-+a^+)\), \(\ \jump{a} = a^- - a^+\)
• Numerical fluxes: one value per face, the only coupling

(1.3 min) Galerkin on each polygon: multiply by a test function of the cell, integrate, and move the derivatives onto the test function. The continuity equation is a conservation law, so its divergence goes to the boundary; the momentum equation is in vector-invariant form, so only the Bernoulli gradient is integrated by parts, and the rotational term stays a volume integral once zeta_h is known. The boundary values are two-valued on each face, one trace from each side, and are replaced by single-valued numerical fluxes (green): the only coupling between cells, and the place where all the stabilization lives. Summing over the cells gives M dw/dt = R(w), with a block-diagonal M: explicit RK4, no global solve. Next: which fluxes.

Fluxes: energy-conjugate, normal penalty

\[ \begin{aligned} \widehat{\eta\mathbf{u}}\cdot\mathbf{n} &= \avg{\mathbf{F}_h}\cdot\mathbf{n} + \frac{\alpha}{g}\,\jump{B_h}\\[6pt] \widehat{B} &= \avg{B_h} + \frac{\alpha}{\bar\eta_F}\,{\color{#2e7d32}\jump{\mathbf{F}_h\cdot\mathbf{n}}} \end{aligned} \]

• \(\mathbf{F}_h = \eta_h\mathbf{u}_h\) and \(B_h\): energy-conjugate pair
  \(\delta E/\delta\eta = B\), \(\ \delta E/\delta\mathbf{u} = \eta\mathbf{u}\)
• Each flux penalized by the jump of the other variable
• Only the normal mass-flux jump is penalized
• Penalty strength \(\alpha \ge 0\): a local wave speed (HLL)

\(\bar\eta_F\): mean depth on the face

\(\alpha = \beta\,\dfrac{|s_L s_R|}{s_R - s_L}\)

Harten–Lax–van Leer (HLL): \(s_L\), \(s_R\) the slowest and the fastest of \(\bar{\mathbf{u}}\cdot\mathbf{n} \mp \sqrt{g\bar\eta}\) on the two sides

Inspired by Ricardo, Lee & Duru, J. Comput. Phys. (2024): collocated vector-invariant DG on the cubed sphere

(1.5 min) Central fluxes plus a penalty. The key choice: each flux is penalized by the jump of its energy-conjugate variable, the depth flux by the jump of the Bernoulli function, the Bernoulli flux by the jump of the normal mass flux. This is the vector-invariant analogue of the entropy-variable dissipation of entropy-stable schemes, and it is what makes the penalty dissipative in the nonlinear energy budget. In the linear regime it reduces to the classical penalties, and with beta = 1 to the upwind (Roe) flux. The momentum equation is penalized only in the normal direction; the coefficient is the local HLL wave speed.

Vorticity: upwind weak curl

• The tangential jump is not penalized: it carries geostrophic balance
• It enters only through the vorticity \(\zeta_h\), a weak curl: \(\zeta = \partial_x v - \partial_y u\) integrated by parts, with one tangential velocity per face, taken upwind

\[ (\zeta_h,\phi)_\kappa = -(\mathbf{u}_h,\nabla^{\perp}\phi)_\kappa + \langle \phi,\, \mathbf{u}_h^{\uparrow}\cdot\mathbf{t}\rangle_{\partial\kappa} \]

• \(\mathbf{u}_h^{\uparrow} = \mathbf{u}_h^-\) or \(\mathbf{u}_h^+\): the upwind side, by the sign of the mean \(\avg{\mathbf{u}_h}\cdot\mathbf{n}\) on the face

• \((\zeta_h+f)\,\mathbf{u}_h^{\perp}\): volume term only

Weak curl with average trace: Ricardo, Lee & Duru, J. Comput. Phys. (2024); upwind trace (\(\gamma = 1\) in the paper): this work

(1.2 min) Penalizing the full velocity jump is the natural first choice, but it damps the tangential jump that a balanced flow carries across faces, and in our runs it needed larger beta and was twice as dissipative. So the tangential jump is left to the vorticity: a weak curl whose face trace is single-valued. The weak curl with the average trace is the one of Ricardo, Lee & Duru (2024) for collocated DG on the cubed sphere; here the trace is upwinded once per face (gamma = 1), from the sign of the mean normal velocity, so it stays polynomial along the face. The rotational term is then a pure volume integral. Figure: the two traces at a point of the face; the mean normal velocity points into kappa^+, so kappa^- is upwind and the weak curl uses u_h^- . t, whatever u_h^+ . t is.

What we prove, on any polygonal mesh and any \(p\)

✓ Mass	locally and globally conserved, every \(\beta\)
✓ Absolute vorticity	\(\int_\Omega(\zeta_h+f) = f\,\lvert\Omega\rvert\) exactly (weak curl)
✓ Linear energy	\(\dfrac{dE_h}{dt} = -\sum_F \alpha_F\left(g\,\lVert\jump{\eta_h'}\rVert^2 + H\,\lVert\jump{\mathbf{u}_h\cdot\mathbf{n}}\rVert^2\right) \le 0\)
✓ Linear balance	\(\eta_h' = \tfrac{f}{g}\psi_h\), \(\ \mathbf{u}_h = \nabla^{\perp}\psi_h\), \(\ \psi_h \in C^0\): steady, every \(\beta\)
✓ Nonlinear energy	penalties dissipate, up to a projection defect
✓ Lake at rest	preserved, with the bottom in \(V_h^p\)

No mesh regularity or convexity needed: stability for every \(\beta \ge 0\) in the linear regime

(1.5 min) Results of Section 3.5, none needs shape regularity or convex cells. Mass and total absolute vorticity are conserved exactly (the element-wise vorticity would not conserve the latter). The linearized scheme has an exact energy identity: the penalties only remove energy, so it is stable for every beta; beta is a matter of accuracy, not stability. Discrete linear geostrophic balance is exact for states built from a continuous streamfunction, for every beta; the full Cartesian penalty would break it. (Caveat if asked: on these Voronoi meshes the C0 subspace of V_h^p is non-trivial from p = 4; for a smooth streamfunction the balance holds up to face jumps that vanish under refinement.) In the nonlinear regime the energy-conjugate pairing makes the penalty terms a sum of squares, up to a defect quadratic in the projection errors, with no Froude-number condition. With topography in the discrete space, the lake at rest is kept. Each row links to a backup slide with the statement and the proof (click it during questions).

Test bed: Voronoi meshes on the plane

• Earth-sized torus, \(H = 10\) km, \(f\)-plane
• \(p = 1,\dots,4\), RK4, \(\beta = 0.5\), up to 65 536 cells

• Same meshes: TRSK (MPAS) and Peixoto (2016)
• Balanced zonal flow, vortex, unstable jet

(0.6 min) Four families of bi-periodic Voronoi meshes, generated with VoronoiMeshes.jl: near-regular hexagons; quasi-uniform centroidal meshes, whose pentagon-heptagon defects (coloured) reproduce the irregularity of spherical grids; locally refined meshes (2x finer core, as in MONAN); and meshes with increasing random distortion at fixed resolution. Reference: the TRSK finite-volume scheme of MPAS and its consistent variant (Peixoto 2016), run on the same meshes.

Convergence: order \(p+1\) on every mesh

Unknowns per cell: TRSK 4 (1 depth, 3 normal velocities) · PolyDG 9, 18, 30, 45 for \(p = 1\)–4 Finest quasi-uniform mesh (65 536 cells): 43 slightly obtuse dual (Delaunay) triangles

PolyDG keeps order \(p+1\) on every mesh; TRSK drops to first order on the quasi-uniform one

(1.2 min) Depth error of the zonal balance after one day against the number of cells. PolyDG converges at p+1 on the three families: observed orders 2.1, 3.0, 4.1, 5.0 on regular hexagons and within 0.1 on the quasi-uniform meshes, p = 4 reaching 1e-13. TRSK (dashed) is second order on regular hexagons, where it beats PolyDG p = 1, but drops to about first order on the quasi-uniform family, and both finite volumes suffer from the slightly obtuse Delaunay triangles of the finest quasi-uniform mesh (43 of 131 072, largest angle 96.8 degrees; the coarser quasi-uniform meshes have none): their error jumps on the last step, PolyDG’s does not. The order of PolyDG does not depend on the cell shape. The stationary vortex gives the same picture (backup slide).

Robust to cell distortion

• 11 600 cells, random displacement of the generators, \(d = 0 \to 0.3\)
• PolyDG error dominated by the projection: insensitive to the cell shape

PolyDG \(p = 1\) depth error, \(d = 0\) and \(d = 0.3\)

Same resolution, distorted cells: PolyDG error bands follow the flow, not the cells

(1 min) Fixed resolution (11 600 cells), cell geometry degraded from d = 0 to 0.3. PolyDG’s total error barely moves (at most 1.2 %, for every p): it is dominated by the projection error, which does not care about the cell shape. TRSK’s error grows 14 times, the consistent variant’s almost doubles. On the undistorted mesh TRSK beats PolyDG p = 1; on the most distorted one it is 4.4 times worse. Right: the p = 1 error keeps the zonal structure of the flow on both meshes.

Local refinement without imprinting

• Zonal balance, 16 384 cells, core 2× finer, day 1
• PolyDG \(p = 1\): error 5× smaller in the refined core

• TRSK: as large in the core as outside, largest on the transition
• Peixoto (2016): dipoles along the transition

(1 min) Depth error maps on the locally refined mesh (PolyDG p = 1 error magnified 5x to share one colour scale with the finite volumes; p = 2 behaves the same, 1000 times smaller). PolyDG behaves as one would like: smaller error in the refined core, no special feature on the transition band. TRSK has errors of the same size inside and outside, sign-changing from cell to cell, and largest on the band. The consistent variant improves the core but puts dipoles along the band.

Barotropically unstable jet

▶ Video, 20 days: www.ime.usp.br/~pedrosp/talks/2026-polimi/figures/video-jet.mp4

(1 min) Planar analogue of Galewsky et al. (2004): two opposite jets in geostrophic balance, perturbed by a small bump; here the northern jet over 20 days, one frame per hour. Top: p = 4 on 4096 cells, 184 320 DOF (eta, u, v), cells visible through the jumps across faces. Bottom: p = 1 on 65 536 cells, 589 824 DOF. Zonal up to day 6, waves by day 8, rolled up by day 10, a turbulent vortex street by day 20, with no grid-scale noise.

Penalty: needed, but not tuned

\(\beta = 0\): blow-up or noise. Any \(\beta \ge 0.1\): clean, results barely depend on \(\beta\)

(0.5 min) The linear scheme is stable for every beta, the nonlinear one needs a penalty: without it the jet diverges within a day at p = 1, 3, 4 and fills with grid-scale noise at p = 2. Any beta from 0.1 removes the noise, energy then decreases monotonically, and from 0.1 to 0.8 the vortices are the same at day 10 and differ only in detail by day 20. We use 0.5 everywhere.

High order pays off

• Day 10, same window, comparable DOF: all unknowns \((\eta, u, v)\)
• \(p = 3\) on 4096 cells forms the vortices that \(p = 1\) smears on 16 384
• \(p = 2\) on 16 384 cells beats \(p = 1\) on 65 536 with half the DOF

At equal cost, high order wins

(1 min) Same window at day 10, low order on the left, higher order on the right, at comparable numbers of DOF (eta, u and v together). p = 3 on 4096 cells (30 DOF per cell, 122 880 DOF) forms the vortices that p = 1 on 16 384 cells (9 DOF per cell, 147 456 DOF) smears out, and p = 2 on 16 384 cells resolves the filaments better than p = 1 on 65 536 cells, with half the DOF. From 16 384 cells on, every p >= 2 gives the same roll-up (paper Fig. 9).

Local refinement in the jet

• 4096 cells, day 14
• North jet crosses the refined core: thinner filaments
• South jet, its mirror image: background cells visible at \(p = 2\)
• No feature marks the transition band

(0.5 min) The initial state is symmetric under a 180-degree rotation, so the southern jet (right, rotated) is what the northern one would be without refinement. The refined core resolves finer scales, and nothing marks the transition band. With 16 384 cells the two trains differ only in detail.

Parallel efficiency

• One ghost layer, one exchange per RK stage
• Block-diagonal mass matrix: no global solve

• ≈ 90 % efficiency on a 48-core node, \(p = 2\)–4
• Beyond a node (TCP): less loss at high \(p\)

(0.7 min) Vulpes, MPI, one rank per core, 65 536 cells (1.2 to 2.9 million unknowns). Within a 48-core node the efficiency is 87-92 % relative to one core, and one core updates about 0.2 million degrees of freedom per second. Across nodes, over 25 Gb/s Ethernet without RDMA, the efficiency at 480 cores is 61-74 % relative to one node, better at higher degree since the work per step grows faster with p than the communication.

Take-home messages

• PolyDG for the rotating shallow water equations on any polygonal grid
• Energy-conjugate, normal penalty + upwind weak curl: proven conservation, linear stability and balance
• Order \(p+1\) on irregular, distorted and refined Voronoi meshes: no grid imprinting
• High order pays off, and stays local in space

Next: curved polygons on spherical Voronoi grids, towards MPAS and MONAN

P. S. Peixoto, F. A. V. B. Alves, I. Fumagalli, P. F. Antonietti, Polytopal discontinuous Galerkin methods for the rotating shallow water equations on polygonal grids (2026)

(0.8 min) Summary in four lines, then the outlook: the spherical geometry, with curved polygonal elements on spherical Voronoi meshes, and the coupling with dynamical cores built on the same grids, such as MPAS and MONAN.

Thanks!

github.com/vulpeslib/vulpes · ppeixoto [at] usp.br · www.ime.usp.br/~pedrosp

Funded by the European Union (ERC Synergy Grant NEMESIS, 101115663). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

Funding: ERC Synergy Grant NEMESIS (101115663), FAPESP (2021/06176-0, 2023/00447-7, 2025/02678-1), CNPq (301553/2026-2). Computing: kami cluster and nemesis server, Dipartimento di Matematica, Politecnico di Milano.

Backup: mass

Proposition (mass). For every \(p \ge 0\), every polygonal mesh and every \(\beta \ge 0\), the scheme is locally conservative and \(\dfrac{d}{dt}\displaystyle\int_\Omega \eta_h = 0\).

Proof. Take \(\phi = 1\) on \(\kappa\) in the continuity equation. The volume term \((\eta_h\mathbf{u}_h, \nabla 1)_\kappa\) vanishes, so \[ \frac{d}{dt}\int_\kappa \eta_h = -\int_{\partial\kappa} \widehat{\eta\mathbf{u}}\cdot\mathbf{n}_\kappa . \] The flux is single-valued and \(\mathbf{n}_\kappa\) changes sign across a face, so the face terms cancel in the sum over the cells of the bi-periodic domain. \(\square\)

• Each cell changes only through the fluxes it shares with its neighbours
• \(\int_\Omega \eta_h\) is linear in the state: it holds for RK4 as well
• Needs single-valued fluxes only: no mesh regularity, any \(\beta\)

Backup: total absolute vorticity

Proposition (total absolute vorticity). With the weak curl, \(\int_\Omega \zeta_h = 0\) at all times, so \(\int_\Omega (\zeta_h + f) = f\,|\Omega|\) is conserved exactly.

Proof. Take \(\phi = 1\) in the weak curl. Since \(\nabla^\perp 1 = 0\), \[ (\zeta_h, 1)_\kappa = \langle 1,\, \mathbf{u}_h^{\uparrow}\cdot\mathbf{t} \rangle_{\partial\kappa} . \] The trace \(\mathbf{u}_h^{\uparrow}\) is single-valued on each face, and the two neighbours run through the face with opposite tangents \(\mathbf{t}\), so the face terms cancel in the sum over the cells. \(\square\)

• The element-wise curl \(\partial_x v_h - \partial_y u_h\) gives \(\int_\Omega \zeta_h = \sum_F \int_F \jump{\mathbf{u}_h\cdot\mathbf{t}} \ne 0\)
• Holds for the average trace and the upwind trace alike
• Potential enstrophy and the other Casimirs are not conserved

Backup: energy, linear regime

Proposition (energy identity). Linearized about rest, with \(E_h = \tfrac12 \int_\Omega \big( g\,(\eta_h')^2 + H\,|\mathbf{u}_h|^2 \big)\), for every \(p\), every polygonal mesh and every \(\beta \ge 0\), \[ \frac{dE_h}{dt} = -\sum_F \alpha_F \Big( g\,\lVert \jump{\eta_h'} \rVert^2_{L^2(F)} + H\,\lVert \jump{\mathbf{u}_h\cdot\mathbf{n}} \rVert^2_{L^2(F)} \Big) \le 0 . \]

Proof. Test the continuity equation with \(g\,\eta_h'\) and the momentum equation with \(H\,\mathbf{u}_h\), both in \(V_h^p\). The Coriolis term gives \(Hf\,\mathbf{u}_h^\perp\cdot\mathbf{u}_h = 0\) pointwise. The volume terms combine into the boundary term \(gH\int_{\partial\kappa} \eta_h'\,\mathbf{u}_h\cdot\mathbf{n}\). On each face, by \(a^-b^- - a^+b^+ = \avg{a}\jump{b} + \jump{a}\avg{b}\), it cancels the central parts of the fluxes, and the penalty parts are left. \(\square\)

• Stable for every \(\beta \ge 0\), and \(\beta = 0\) conserves the energy exactly
• So \(\beta\) is a matter of accuracy and grid-scale noise, not of stability

Backup: energy, linear regime, in the runs

• 3. Plane wave, \(p = 2\): \(\beta = 0\) conserves the energy (RK4 error only), \(\beta > 0\) decreases it monotonically
• 4. Geostrophic balance: the full Cartesian penalty dissipates 1.5–1.9 times more than the normal one

Backup: geostrophic balance, linear regime

Proposition (discrete balance). Let \(\psi_h \in V_h^p\) be continuous, \(\eta_h' = \tfrac{f}{g}\,\psi_h\) and \(\mathbf{u}_h = \nabla^\perp \psi_h\). With the normal penalty, this state is exactly steady for the linear scheme, for every \(\beta \ge 0\).

Proof. In each cell, \(\nabla\cdot\nabla^\perp\psi_h = 0\), and \(f\,\mathbf{u}_h^\perp = -f\,\nabla\psi_h\) cancels \(g\,\nabla\eta_h' = f\,\nabla\psi_h\). On each face the residuals combine \[ \jump{\eta_h'} = \tfrac{f}{g}\,\jump{\psi_h} = 0 \quad\text{and}\quad \jump{\mathbf{u}_h\cdot\mathbf{n}} = -\jump{\partial_{\mathbf{t}}\psi_h} = 0, \] since \(\psi_h\) is continuous and so is its derivative along the face. \(\square\)

• A full velocity penalty adds \(\alpha\jump{\mathbf{u}_h\cdot\mathbf{t}} = \alpha\jump{\partial_{\mathbf{n}}\psi_h} \ne 0\): the balance is lost
• Continuous \(\psi_h\) beyond the constants exist from \(p = 4\) on these meshes
• A smooth \(\psi\) is balanced up to face jumps that vanish under refinement

Backup: energy, nonlinear regime

Proposition (energy budget). With the projection errors \(e_B = (I - \Pi_h) B_h\), \(\mathbf{e}_F = (I - \Pi_h) \mathbf{F}_h\), for every \(p\), mesh and \(\alpha \ge 0\) (no Froude condition): \[ \begin{aligned} \frac{dE_h}{dt} &= -\sum_\kappa \mathcal{R}_\kappa - \sum_F \mathcal{A}_F - \sum_F \mathcal{P}_F,\\ -\mathcal{P}_F &\le \mathcal{D}_F = \int_F \frac{\alpha}{4} \Big( \frac{1}{g}\jump{e_B}^2 + \frac{1}{\bar\eta_F}\jump{\mathbf{e}_F\cdot\mathbf{n}}^2 \Big). \end{aligned} \]

Proof. Test with \(\Pi_h B_h\) and \(\Pi_h \mathbf{F}_h\), the projections of \(\delta E/\delta\eta\) and \(\delta E/\delta\mathbf{u}\). Since \(\jump{\Pi_h B_h} = \jump{B_h} - \jump{e_B}\), completing the square gives \(\jump{B_h}\jump{\Pi_h B_h} = \big(\jump{B_h} - \tfrac12\jump{e_B}\big)^2 - \tfrac14 \jump{e_B}^2\), and the same for \(\mathbf{F}_h\cdot\mathbf{n}\). \(\square\)

• Energy-conjugate pairing: a sum of squares up to \(\mathcal{D}_F\); same-variable penalties are indefinite when \(\avg{\mathbf{u}\cdot\mathbf{t}} \ne 0\)
• \(\mathcal{R}_\kappa\) (rotation) and \(\mathcal{A}_F\) involve projection errors: a full estimate is open

Backup: lake at rest

Remark (lake at rest). With the bottom \(b\) in \(V_h^p\), the state \(\mathbf{u}_h = 0\), \(\eta_h + b = H\) is steady, for every \(\beta \ge 0\).

Proof. At rest \(B_h = g(\eta_h + b) = gH\) is constant and \(\mathbf{F}_h = \eta_h\mathbf{u}_h = 0\). The Bernoulli volume term cancels its central face term (integration by parts of a constant), the mass-flux terms vanish with \(\mathbf{F}_h\), and the penalties vanish since \(\jump{B_h} = 0\) and \(\jump{\mathbf{F}_h} = 0\). \(\square\)

• The energy-conjugate pairing does it: a depth penalty \(\alpha\jump{\eta_h} = -\alpha\jump{b}\) would not
• Positivity of \(\eta_h\) is not guaranteed (deep water here)
• A property of the formulation: topography is not implemented in the solver

Backup: convergence, both steady states

Backup: quadrature-free integration

Euler’s identity + divergence theorem, monomial \(m = x^a y^b\): \[ \int_\kappa m \,d\mathbf{x} = \frac{1}{2+a+b}\sum_{F\subset\partial\kappa} (\mathbf{x}_F\cdot\mathbf{n}_\kappa) \int_F m\, ds \]

• Recursively down to the face endpoints: exact, no sub-triangles
• Bilinear forms (mass matrices): once, at set-up
• Nonlinear residual: sub-triangulation (cheaper for products of fields)

Antonietti, Houston & Pennesi, J. Sci. Comput. (2018); Radley, Houston & Hubbard, Math. Eng. (2024)