Entropy-Stable Discontinuous Spectral-Element Methods for the Spherical Shallow Water Equations in Covariant Form

We introduce discontinuous spectral-element methods of arbitrary order that are well balanced, conservative of mass, and conservative or dissipative of total energy (i.e., a mathematical entropy function) for a covariant flux formulation of the rotating shallow water equations with variable bottom topography on curved manifolds such as the sphere. The proposed methods are based on a skew-symmetric splitting of the tensor divergence in covariant form, which we implement and analyze within a general flux-differencing framework using tensor-product summation-by-parts operators. Such schemes are proven to satisfy semi-discrete mass and energy conservation on general unstructured quadrilateral grids in addition to well balancing for arbitrary continuous bottom topographies, with energy dissipation resulting from a suitable choice of numerical interface flux. Furthermore, the proposed covariant formulation permits an analytical representation of the geometry and associated metric terms while satisfying the aforementioned entropy stability, conservation, and well-balancing properties without the need to approximate the metric terms so as to enforce discrete metric identities. Numerical experiments on cubed-sphere grids are presented in order to verify the schemes’ structure-preservation properties as well as to assess their accuracy and robustness within the context of several standard test cases characteristic of idealized atmospheric flows. Our theoretical and numerical results support the further development of the proposed methodology towards a full dynamical core for numerical weather prediction and climate modelling, as well as broader applications to other hyperbolic and advection-dominated systems of partial differential equations on curved manifolds.

💡 Research Summary

The paper presents a high‑order, entropy‑stable discontinuous spectral‑element (DG) framework for solving the rotating shallow‑water equations (SWE) on curved manifolds such as the sphere. By formulating the SWE in a covariant tensor form, the authors retain the geometric information (metric tensor, Christoffel symbols) analytically, avoiding any discrete metric approximations that would otherwise be required to enforce metric identities.

The continuous system is expressed in terms of covariant momentum components and a flux tensor τᵢⱼ = h vᵢ vⱼ + (g/2) h² Gᵢⱼ. An entropy (total energy) function η = ½ h vᵢ vᵢ + (g/2) h (h + b) is identified, and the associated entropy flux Fⱼ is derived. The authors then rewrite the equations in a split‑form that is skew‑symmetric with respect to the covariant derivative. This split‑form is crucial because it guarantees that the discrete volume terms do not create or destroy entropy; any entropy change originates solely from the numerical interface fluxes.

Spatial discretisation relies on tensor‑product Summation‑by‑Parts (SBP) operators. One‑dimensional SBP operators built on Gauss‑Lobatto or Gauss‑Legendre nodes are tensor‑producted to obtain two‑dimensional operators that mimic integration‑by‑parts on each element. The covariant derivative is approximated by a flux‑differencing formulation that uses the SBP property to move derivatives onto test functions, thereby preserving the discrete analogue of the continuous energy identities.

Two families of two‑point numerical fluxes are introduced. A central flux yields a scheme that is entropy‑conservative (entropy neither grows nor decays). An upwind‑type entropy‑stable flux adds a dissipation term proportional to the jump in entropy variables, guaranteeing entropy decay consistent with the second law of thermodynamics. The choice of flux therefore controls whether the method is strictly entropy‑conservative or entropy‑dissipative.

Source terms arising from the Coriolis force, variable bottom topography, and geometric curvature are split in a manner that preserves the well‑balanced property: the scheme exactly maintains a steady state with constant free‑surface height and zero velocity, even on non‑uniform meshes and with arbitrary smooth bottom topography.

The methodology is implemented on cubed‑sphere grids, a popular quasi‑uniform discretisation of the sphere. Numerical experiments include: (i) an unsteady analytical solution to verify order‑of‑accuracy (demonstrating optimal convergence for polynomial degrees ≥ 4); (ii) the standard isolated‑mountain test (Williamson case 5) and the Rossby‑Haurwitz wave (case 6) to assess well‑balancing and long‑time stability; (iii) the barotropic instability problem of Galewsky et al. to test nonlinear wave growth and transition. In all cases, mass is conserved to machine precision, total energy is either exactly conserved (central flux) or monotonically dissipated (entropy‑stable flux), and the steady‑state balance is maintained to round‑off error. The results also show that the scheme remains stable without any ad‑hoc filtering or artificial viscosity, even on coarse meshes, confirming the robustness imparted by the SBP‑based split‑form and entropy‑stable fluxes.

In summary, the authors extend the SBP‑based split‑form DG methodology to covariant formulations on curved manifolds, delivering a discretisation that simultaneously achieves arbitrary high‑order accuracy, discrete mass and energy conservation (or controlled dissipation), and exact well‑balancing for the spherical shallow‑water equations. The work paves the way for next‑generation dynamical cores in numerical weather prediction and climate modelling, and suggests straightforward extensions to three‑dimensional nonhydrostatic atmospheric models, multi‑physics coupling, and operational forecasting systems.