Global Flux Quadrature for Shallow Water DG-SEM

• Global Flux Quadrature is a strategy for hyperbolic balance laws that embeds the source term into a global flux primitive, ensuring constant flux at steady state.
• It leverages Gauss–Lobatto quadrature within the DG-SEM framework to replace traditional source integrations with collocation-based approximations.
• The method achieves well-balanced, superconvergent discretizations with exact preservation of equilibrium states in shallow water models.

Searching arXiv for the primary paper and closely related works on Global Flux Quadrature. First, retrieving the main paper by arXiv ID (Mantri et al., 2022). Global Flux Quadrature is a discretization strategy for hyperbolic balance laws in which the source term is absorbed into an additional flux primitive, so that the governing equation is rewritten in a quasi-conservative form and the discrete steady states are characterized by constant global flux. In the discontinuous Galerkin spectral element method (DG-SEM) for the shallow water equations, this construction replaces the source integral by a local quadrature based on Gauss–Lobatto nodes, yields a discrete well-balanced property without any a-priori knowledge of the steady equilibrium, does not involve the explicit solution of any local auxiliary problem, and is equivalent at steady state to LobattoIIIA collocation of the steady flux ODE (Mantri et al., 2022).

1. Governing equations and discrete setting

The primary setting is the one-dimensional shallow water system with bathymetry $b(x)$,

$$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$

with

$$U = \begin{bmatrix} h \ hu \end{bmatrix},\qquad F(U) = \begin{bmatrix} hu \ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad S(U,b,\partial_x b) = \begin{bmatrix} 0 \ - g h \,\partial_x b \end{bmatrix}.$$

For the pseudo-1D rotating shallow water model with friction, the source becomes

$$S = -h \begin{bmatrix} 0 \ \partial_x \varphi + c_f u + \omega v \ -\omega u \end{bmatrix},\qquad \varphi = g b.$$

The classical lake-at-rest equilibrium is

$u=0,\qquad h+b=H_0,$

so that the free surface $\zeta=h+b$ is constant (Mantri et al., 2022).

In the DG-SEM formulation, each element is mapped from a reference interval and discretized with nodal Gauss–Lobatto points $\{\xi_i\}_{i=0}^p$ and Lagrange basis functions $\{\phi_i\}$ of degree $p$. The nodal interpolation is

$u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.$

The mass matrix is diagonal,

$$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$0

with Gauss–Lobatto weights $$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$1, and the derivative operator is represented by $$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$2, together with the SBP-scaled form

$$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$3

and the boundary sampling matrix

$$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$4

The standard semi-discrete DG-SEM strong residual form is

$$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$5

or, with explicit numerical flux

$$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$6

equivalently

$$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$7

2. Global flux construction and quadrature formula

The central idea is the global flux approach: construct an additional flux $$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$8 as a primitive of the source so that steady states satisfy

$$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$9

With

$$U = \begin{bmatrix} h \ hu \end{bmatrix},\qquad F(U) = \begin{bmatrix} hu \ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad S(U,b,\partial_x b) = \begin{bmatrix} 0 \ - g h \,\partial_x b \end{bmatrix}.$$0

the equation is rewritten as

$$U = \begin{bmatrix} h \ hu \end{bmatrix},\qquad F(U) = \begin{bmatrix} hu \ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad S(U,b,\partial_x b) = \begin{bmatrix} 0 \ - g h \,\partial_x b \end{bmatrix}.$$1

At steady state, $$U = \begin{bmatrix} h \ hu \end{bmatrix},\qquad F(U) = \begin{bmatrix} hu \ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad S(U,b,\partial_x b) = \begin{bmatrix} 0 \ - g h \,\partial_x b \end{bmatrix}.$$2 is equivalent to

$$U = \begin{bmatrix} h \ hu \end{bmatrix},\qquad F(U) = \begin{bmatrix} hu \ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad S(U,b,\partial_x b) = \begin{bmatrix} 0 \ - g h \,\partial_x b \end{bmatrix}.$$3

so constant global flux is the discrete analogue of equilibrium (Mantri et al., 2022).

In nodal spectral elements, the source primitive is represented through a local cumulative quadrature. If

$$U = \begin{bmatrix} h \ hu \end{bmatrix},\qquad F(U) = \begin{bmatrix} hu \ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad S(U,b,\partial_x b) = \begin{bmatrix} 0 \ - g h \,\partial_x b \end{bmatrix}.$$4

the integration matrix is

$$U = \begin{bmatrix} h \ hu \end{bmatrix},\qquad F(U) = \begin{bmatrix} hu \ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad S(U,b,\partial_x b) = \begin{bmatrix} 0 \ - g h \,\partial_x b \end{bmatrix}.$$5

On an element of size $$U = \begin{bmatrix} h \ hu \end{bmatrix},\qquad F(U) = \begin{bmatrix} hu \ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad S(U,b,\partial_x b) = \begin{bmatrix} 0 \ - g h \,\partial_x b \end{bmatrix}.$$6, the nodal values of the primitive are

$$U = \begin{bmatrix} h \ hu \end{bmatrix},\qquad F(U) = \begin{bmatrix} hu \ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad S(U,b,\partial_x b) = \begin{bmatrix} 0 \ - g h \,\partial_x b \end{bmatrix}.$$7

where $$U = \begin{bmatrix} h \ hu \end{bmatrix},\qquad F(U) = \begin{bmatrix} hu \ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad S(U,b,\partial_x b) = \begin{bmatrix} 0 \ - g h \,\partial_x b \end{bmatrix}.$$8 is a constant vector, taken as the last value of the left neighbor to avoid jumps for continuous data.

The DG-SEM written in terms of the global flux becomes

$$U = \begin{bmatrix} h \ hu \end{bmatrix},\qquad F(U) = \begin{bmatrix} hu \ \dfrac{(hu)^2}{h} + \dfrac{1}{2} g h^2 \end{bmatrix},\qquad S(U,b,\partial_x b) = \begin{bmatrix} 0 \ - g h \,\partial_x b \end{bmatrix}.$$9

Using the explicit numerical flux and the definition of $$S = -h \begin{bmatrix} 0 \ \partial_x \varphi + c_f u + \omega v \ -\omega u \end{bmatrix},\qquad \varphi = g b.$$0, the semi-discrete form is

$$S = -h \begin{bmatrix} 0 \ \partial_x \varphi + c_f u + \omega v \ -\omega u \end{bmatrix},\qquad \varphi = g b.$$1

This yields the global flux quadrature formula: $$S = -h \begin{bmatrix} 0 \ \partial_x \varphi + c_f u + \omega v \ -\omega u \end{bmatrix},\qquad \varphi = g b.$$2 Hence the source contribution reduces to the local matrix–vector product

$$S = -h \begin{bmatrix} 0 \ \partial_x \varphi + c_f u + \omega v \ -\omega u \end{bmatrix},\qquad \varphi = g b.$$3

which requires only the nodal source values and the precomputed integration tableau $$S = -h \begin{bmatrix} 0 \ \partial_x \varphi + c_f u + \omega v \ -\omega u \end{bmatrix},\qquad \varphi = g b.$$4 defined by the Gauss–Lobatto nodes (Mantri et al., 2022).

3. Collocation equivalence, superconvergence, and well-balancedness

At steady state, the method is governed by the flux ODE

$$S = -h \begin{bmatrix} 0 \ \partial_x \varphi + c_f u + \omega v \ -\omega u \end{bmatrix},\qquad \varphi = g b.$$5

with initial condition $$S = -h \begin{bmatrix} 0 \ \partial_x \varphi + c_f u + \omega v \ -\omega u \end{bmatrix},\qquad \varphi = g b.$$6. Provided $$S = -h \begin{bmatrix} 0 \ \partial_x \varphi + c_f u + \omega v \ -\omega u \end{bmatrix},\qquad \varphi = g b.$$7 is invertible and the inverse $$S = -h \begin{bmatrix} 0 \ \partial_x \varphi + c_f u + \omega v \ -\omega u \end{bmatrix},\qquad \varphi = g b.$$8 is uniquely defined and bounded away from singular or critical points, the global flux Gauss–Lobatto DG-SEM admits a nodally continuous discrete steady state obtained by integrating this ODE with the fully implicit LobattoIIIA collocation method having $$S = -h \begin{bmatrix} 0 \ \partial_x \varphi + c_f u + \omega v \ -\omega u \end{bmatrix},\qquad \varphi = g b.$$9 stages at the Gauss–Lobatto nodes. On each element,

$u=0,\qquad h+b=H_0,$0

hence

$u=0,\qquad h+b=H_0,$1

and continuity of $u=0,\qquad h+b=H_0,$2 implies continuity of $u=0,\qquad h+b=H_0,$3 (Mantri et al., 2022).

The collocation system is

$u=0,\qquad h+b=H_0,$4

This is exactly the discrete equation enforced by global flux quadrature at steady state. Under the additional assumption that there exists a flux linearization

$u=0,\qquad h+b=H_0,$5

with $u=0,\qquad h+b=H_0,$6 diagonalizable and eigenvalues bounded away from zero,

$u=0,\qquad h+b=H_0,$7

the discrete steady state inherits the collocation superconvergence:

• at element endpoints, order $u=0,\qquad h+b=H_0,$8 in $u=0,\qquad h+b=H_0,$9;
• at internal Gauss–Lobatto nodes, order $\zeta=h+b$0.

The well-balanced property is formulated directly at the semi-discrete level. The residual is

$\zeta=h+b$1

If

$\zeta=h+b$2

at the Gauss–Lobatto nodes and

$\zeta=h+b$3

then

$\zeta=h+b$4

the interface terms vanish, and the discrete residual is zero. In this sense, any smooth steady state driven by $\zeta=h+b$5 is preserved as a collocation steady state (Mantri et al., 2022).

Exact preservation of lake-at-rest requires a modified nodal source,

$\zeta=h+b$6

with $\zeta=h+b$7 and $\zeta=h+b$8. For $\zeta=h+b$9 constant and $\{\xi_i\}_{i=0}^p$0,

$\{\xi_i\}_{i=0}^p$1

so

$\{\xi_i\}_{i=0}^p$2

and the lake-at-rest state is preserved exactly.

This construction differs from hydrostatic reconstruction and path-conservative formulations. Hydrostatic reconstruction modifies local reconstructions to enforce exact balance at prescribed equilibria and often requires a priori knowledge of the equilibrium form. Path-conservative formulations redesign flux/source discretizations for nonconservative systems relying on prescribed paths or invariants. Global flux quadrature instead embeds the source integral through a collocation-based primitive $\{\xi_i\}_{i=0}^p$3, requires no auxiliary local steady solves, and no explicit reconstruction of the equilibrium (Mantri et al., 2022).

4. Entropy control and cell entropy correction

The method is paired with an entropy framework. For shallow water without bathymetry in the flux, the entropy pair is

$\{\xi_i\}_{i=0}^p$4

with entropy flux

$\{\xi_i\}_{i=0}^p$5

With spatially constant potential, the total entropy and flux are

$\{\xi_i\}_{i=0}^p$6

In frictionless cases,

$\{\xi_i\}_{i=0}^p$7

whereas with friction,

$\{\xi_i\}_{i=0}^p$8

(Mantri et al., 2022).

To control entropy production, the scheme introduces a symmetric positive definite artificial viscosity at the cell level: $\{\xi_i\}_{i=0}^p$9 with

$\{\phi_i\}$0

Dotting by $\{\phi_i\}$1 and summing yields the cell entropy equation

$\{\phi_i\}$2

The correction parameter is chosen as

$\{\phi_i\}$3

Two numerical entropy fluxes are used. The analytical flux-based option is

$\{\phi_i\}$4

for example $\{\phi_i\}$5, with a local entropy balance consistent with analytically balanced steady states. The global-flux-consistent option is

$\{\phi_i\}$6

which aligns the entropy control with the global flux quadrature and preserves discrete collocation steady states exactly under the correction.

The resulting entropy balance is element-wise: $\{\phi_i\}$7 and globally dissipative in the frictional case,

$\{\phi_i\}$8

5. DG-SEM implementation and numerical behavior

Implementation requires only a minimal modification of standard DG-SEM. One chooses the polynomial degree $\{\phi_i\}$9, precomputes Gauss–Lobatto nodes, weights, basis functions, the mass and differentiation matrices,

$p$0

and a surface numerical flux

$p$1

The volume source contribution is then assembled by precomputing

$p$2

and evaluating the nodal source values $p$3, followed by

$p$4

The semi-discrete update is

$p$5

optionally augmented by the local entropy correction. Explicit Runge–Kutta of order $p$6 may be used, with CFL controlled by the eigenvalues $p$7. Relative to standard DG-SEM, the additional operations are multiplication by the dense $p$8 matrix $p$9, one extra derivative matvec, and optional local cell correction integrals, so overall complexity and memory footprint remain comparable to standard DG-SEM (Mantri et al., 2022).

The numerical tests establish several characteristic properties. For steady moving equilibria, the global flux solution is superconvergent, with order $u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.$0 at endpoints and order $u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.$1 at internal nodes. In transient tests, global flux quadrature reduces steady-state error by orders of magnitude versus classical DG-SEM, thereby enabling accurate resolution of small perturbations. For lake-at-rest with nontrivial bathymetry, exact preservation is obtained with the modified source, and perturbations of amplitude down to $u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.$2 are resolved without spurious oscillations, whereas non-well-balanced DG-SEM produces errors larger than the perturbation amplitude on the same mesh. Trans-critical flows also perform well numerically despite the fact that the formal superconvergence proof does not cover critical points (Mantri et al., 2022).

The paper further reports two-dimensional tests, including perturbations of 1D equilibria on 2D domains, stationary vortices with bathymetry, anticyclonic vortex propagation on a $u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.$3-plane, geostrophic adjustment, and an equatorial Kelvin front. The well-balanced formulation suppresses spurious oscillations present in non-well-balanced schemes, is less diffusive in several of these tests, and exhibits robust entropy behavior.

6. Related formulations, extensions, and limitations

Global Flux Quadrature is part of a broader flux globalization program. In finite volume WENO, the same principle is used by reconstructing a global flux $u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.$4 rather than the conservative variables, with a tailored quadrature for the source primitive $u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.$5; exact discrete steady states are then characterized by constant global fluxes (Ciallella et al., 2022). A finite-difference variant computes the source primitive from multi-step ODE weights, so that discrete steady states align with the underlying Adams–Bashforth or Adams–Moulton integrator and the steady-state accuracy is determined solely by the ODE method order (Kazolea et al., 10 Jan 2025). For shallow water moment equations with non-conservative products, flux globalization integrates both source terms and non-conservative products into the divergence term, again reconstructing the global flux to preserve steady states without prior analytical knowledge of them (Ciallella et al., 1 Jul 2025).

The same viewpoint has been extended beyond DG and finite volumes. For nodal continuous finite elements on Cartesian grids, Global Flux quadrature replaces element-wise volume operators by line- or surface-integrated operators, yielding constraint-compatible stabilization and vorticity-preserving schemes for linear acoustics (Barsukow et al., 2024). A closely related multidimensional formulation rewrites the entire spatial operator as a mixed derivative of a single global flux and proves stationarity preservation and steady-state superconvergence for stabilized nodal finite element methods (Barsukow et al., 3 Oct 2025). A positivity-preserving PAMPA formulation uses Gauss–Lobatto global flux quadrature together with a first-order local Lax–Friedrichs blending to preserve still-water equilibria, positivity of water height, and wet–dry fronts in one-dimensional shallow water models (Abgrall et al., 30 Oct 2025).

Within the DG-SEM shallow water setting, the main assumptions and limitations are explicit. The formal superconvergence proof assumes continuous bathymetry $u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.$6 and potential $u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.$7, smooth steady states, and invertibility of $u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.$8 with bounded inverse away from critical points. Exact lake-at-rest preservation requires the modified nodal source; the standard nodal source yields only a high-order approximation. The two-dimensional extension treated in the paper is dimensionally split, so genuinely multidimensional well-balancing remains future work. Dry/wet fronts are not treated. The paper also identifies several extensions: multidimensional tensor DG-SEM and truly multidimensional global flux tensors, application to other balance laws such as Euler with gravity or MHD, use in finite volume and residual distribution methods, and the incorporation of interface jump terms in $u_h(\xi)=\sum_{j=0}^p \phi_j(\xi)u_j.$9 for discontinuous bathymetry (Mantri et al., 2022).

A plausible implication is that Global Flux Quadrature is best understood not as a single algorithm, but as a structural principle: the source term is embedded into a globally reconstructed primitive so that the discrete operator acts on a quantity that is constant at equilibrium. In the DG-SEM shallow water formulation, this principle is realized by the local Gauss–Lobatto quadrature

$$\partial_t U + \partial_x F(U) = S(U,b,\partial_x b),$$00

the steady-state LobattoIIIA equivalence, and the accompanying cell entropy correction, which together provide high-order accuracy, discrete well-balancing, and controlled entropy production (Mantri et al., 2022).