# Positivity-preserving well-balanced pampa schemes with global flux quadrature for one-dimensional shallow water models

Remi Abgrall, Yongle Liu, Mario Ricchiuto

► To cite this version:

Remi Abgrall, Yongle Liu, Mario Ricchiuto. Positivity-preserving well-balanced pampa schemes with global flux quadrature for one-dimensional shallow water models. 2025. hal-05425468

HAL Id: hal-05425468

https://inria.hal.science/hal-05425468v1

Preprint submitted on 19 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1
2
3
4

POSITIVITY-PRESERVING WELL-BALANCED PAMPA SCHEMES WITH GLOBAL FLUX QUADRATURE FOR ONE-DIMENSIONAL SHALLOW WATER MODELS*

RÉMI ABGRALL†, YONGLE LIU‡, AND MARIO RICCHIUTO§

**Abstract.** We present a novel hydrostatic and non-hydrostatic equilibria preserving Point-Average-Moment PolynomiAl-interpreted (PAMPA) method for solving the one-dimensional hyperbolic balance laws, with applications to the shallow water models including the Saint-Venant system with the Manning friction term and rotating shallow water equations. The idea is based on a global flux quadrature formulation, in which the discretization of the source terms is obtained from the derivative of and additional flux function computed via high order quadrature of the source term. The reformulated system is quasi-conservative with global integral terms computed using Gauss-Lobatto quadrature nodes. The resulting method is capable of preserving a large family of smooth moving equilibria: supercritical and subcritical flows, in a super-convergent manner. We also show that, by an appropriate quadrature strategy for the source, we can exactly preserve the still water states. Moreover, to guarantee the positivity of water depth and eliminate the spurious oscillations near shocks, we blend the high-order PAMPA schemes with the first order local Lax-Friedrichs schemes using the method developed in [2]. The first-order schemes are designed to preserve the still water equilibria and positivity of water height, as well as to deal with wet-dry fronts. Extensive numerical experiments are tested to validate the advantages and robustness of the proposed scheme.

**Key words.** Well-balanced scheme, Point-Average-Moment PolynomiAl-interpreted (PAMPA) method, Shallow water models, Global flux, Positivity-preserving.

MSC codes. 76M10, 76M12, 65M08, 35L40

**1. Introduction.** This paper focus on development of high-order and robust numerical methods for one-dimensional (1-D) shallow water models, which play a fundamental role in various phenomena in natural science and engineering. Such models are given by the following form:

$$
(1.1) \qquad \frac{\partial \mathbf{u}}{\partial t} + \mathbf{f}(\mathbf{u})_x = \mathbf{S}(\mathbf{u}, x), \quad x \in \Omega \subset \mathbb{R}, \quad t \ge 0,
$$

where **u** contains the conserved variables, **f** is the flux function, and **S** denotes the source term that depends on the solution as well as some external data. For example, the influence of spatial changes of the height of the bottom topography, of the Manning friction effect, as well as other physical effects such as Coriolis force, are accounted for by means of appropriately defined this source term.

The addition of the source term results in a system of balance laws which may admit quite a large number of equilibrium solutions dictated by the interaction between different components of the flux variation in space and the forcing terms. Many of such equilibria have some interest in themselves and many physical applications involve small perturbations of such equilibria. The ability of a numerical method to resolve with enhanced accuracy such steady states is a unanimously acclaimed design criterion usually referred to as well-balanced (WB). Even in 1-D case, the challenge of devising WB numerical approximations of shallow water models agnostic of the

<footnote>*Submitted to the editors DATE.</footnote>

<footnote>**Funding:** The work of Y. Liu was supported by UZH Postdoc Grant, 2024 / Verfügung Nr. FK-24-110 and SNSF grant 200020_204917.</footnote>

<footnote>†Institute of Mathematics, University of Zürich, 8057 Zürich, Switzerland</footnote>

<footnote>‡Corresponding author. Institute of Mathematics, University of Zürich, 8057 Zürich, Switzerland (yongle.liu@math.uzh.ch; liuyl2017@mail.sustech.edu.cn).</footnote>

<footnote>§Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, F-33400 Talence, France.</footnote>
form of the steady state is still a very active and open research domain. There is already quite a large literature on the subject, with several different approaches to manage this issue. It is quite impossible to mention all the relevant literatures, we mention a few numerical methods having been proposed in many different settings: finite difference (FD) methods [37, 38, 57, 50, 59, 46], finite volume (FV) methods [7, 14, 18, 25, 44, 36, 33, 30, 23, 11, 22, 12, 21, 24, 13, 60], Discontinuous Galerkin (DG) methods [41, 59, 42, 60, 56, 58, 31, 51], continuous finite element methods [8, 29, 10, 35, 43], residual distribution methods [49, 15, 48, 47, 5, 6], Active Flux methods [9, 4, 40], and so on.

Concerning the preservation of steady states, if the steady-state solutions are explicitly known, one can use a simple and efficient idea that consists in evolving a discrete error with respect to the given equilibrium. This boils down to removing from the discrete equations of any scheme the discrete expression corresponding to the application of the scheme itself to the given solution. A more intricate scenario arises when the relevant equilibrium is not explicitly known, but it can be characterized by a set of (generally nonlinear) algebraic relations that define constant invariants. To put it differently, in this scenario, the equilibrium state can be described through the collection of relationships expressed as $E(u) = E_0$, where $E$ represents a comprehensive set of variables that, ideally, can be leveraged to derive the conserved quantities $u$ and is referred to as equilibrium variable. In this case, one usually has to solve nonlinear equations at each cell and time to proceed the transformation from equilibrium to conserved variables and a nontrivial unisolvent root-finding mechanism should be carefully developed. This idea, known as generalized hydrostatic reconstruction, has been adapted to many discretization approaches; see, e.g., [53, 44, 51, 18, 17, 21, 56, 43].

An alternative WB numerical method, based on the semi-discrete Active Flux approach introduced in [1], was recently proposed in [4] for solving the Saint-Vanent system of shallow water equations with/without friction term. This method is later referred to as Point-Average-Moment PolynomiAl-interpreted (PAMPA) method in [39]. The PAMPA method uses degrees of freedom (DoFs) from quadratic polynomials—all lying on element boundaries—along with the traditional average values within each element. The solution is a global continuous representation of these DoFs (point and average values), which are evolved simultaneously according to two formulations of the same PDEs. The point values can be evolved according to several writings of the same PDEs and can even be nonconservative. For instance, in [4, 40], a primitive formulation with the flux derivative given by equilibrium variable $E$ is used to evolve the point values and the WB property is automatically attained through consistent finite difference approximations of the flux derivative. Whereas, the average values should be updated following the conservative formulation (1.1), due to the well-known Lax–Wendroff theorem. Given the cell boundary point values, the flux is easily computed and the WB evolution of average values only depends on an appropriate approximation of the source term. This not only requires the approximation should be WB but also satisfy exact conservation property or consistent approximation when the source term vanishes. To achieve this, in [4] a cutoff function is used to switch off the inconsistent moving-water equilibrium correction term when the solution is away from steady states. A more elegant way that subtracts a local reference steady-state solution obtained from solving nonlinear equations derived from the relationship $E(u) = E_0$, was used in [40] to ensure the exact conservation property.

In this paper, without assuming either the steady states or local reference steady states a priori, we present a novel flux globalization based positivity-preserving (PP) WB PAMPA scheme for the 1-D shallow water models. This high-order scheme inher-
its the advantages of PAMPA methods in eliminating the need for Riemann solvers and free of any special WB (generalized) hydrostatic reconstructions in solving hyperbolic balance laws [4, 40]. Beyond the positivity-preserving, the low-order schemes are also capable of preserving the still water steady states, which differs from the scheme that we have used in the previous works [4, 40]. Following the idea from [19, 17, 24, 42, 34], we incorporate the source term into the flux and rewrite (1.1) in the following equivalent quasi-conservative form:

$$ (1.2) \qquad \frac{\partial \mathbf{u}}{\partial t} + \partial_x \mathbf{G}(\mathbf{u}, x) = 0, $$

where the global operator (flux) is defined as

$$ (1.3) \qquad \mathbf{G}(\mathbf{u}, x) = \mathbf{f}(\mathbf{u}) - \int_{\hat{x}}^{x} \mathbf{S}(\mathbf{u}, \xi) \, d\xi, $$

which is the natural operator to represent solutions of the steady state ordinary differential equation (ODE)

$$ \partial_x \mathbf{f} = \mathbf{S}. $$

Note that in (1.3) $\hat{x}$ is a given initial integration point, for example the left-end of the domain. We can observe that now any steady-state solution can be naturally expressed by

$$ (1.4) \qquad \mathbf{G} \equiv \text{Const}, $$

thus (1.2) provides a reasonable path to obtain WB schemes capable of preserving a wider class of equilibria. The above condition alone however does not necessarily provide a precise characterization of the discrete equilibria. Here, following [42, 34], we will pursue an approach that allows such characterization in the context of ODE integration methods. To this end a modified version of the PAMPA method will be developed, using the quasi-conservative form (1.2) to update both the cell average values within each element, and point values at the cell boundaries.

We apply the proposed flux globalization based PP WB PAMPA scheme for solving two shallow water models. The first one is the Saint-Venant system with the Manning friction term:

$$ (1.5) \qquad \begin{aligned} h_t + (hu)_x &= 0, \\ (hu)_t + (hu^2 + \frac{g}{2}h^2)_x &= -ghB_x - \frac{gn^2}{h^{\frac{2}{3}}} |hu|hu, \end{aligned} $$

where $h$ is the water depth, $u$ is the velocity, $B$ is bottom topography, $g$ is the constant gravitational acceleration, $n$ is the Manning friction coefficient. The global flux is given by

$$ \mathbf{G} = \begin{pmatrix} hu \\ hu^2 + \frac{g}{2}h^2 + \int_{\hat{x}}^{x} \left[ghB_x + \frac{gn^2}{h^{\frac{2}{3}}} |hu|hu\right] d\xi \end{pmatrix}. $$

The second model is the shallow water equations with the Coriolis forces, which reads as

$$ (1.6) \qquad \begin{aligned} h_t + (hu)_x &= 0, \\ (hu)_t + (hu^2 + \frac{g}{2}h^2)_x &= -ghB_x + fhu, \\ (hv)_t + (huv)_x &= -fhu, \end{aligned} $$
where $u$ and $v$ are the two components of horizontal velocity, $f(x) = f_0 + \beta x$ is the Coriolis parameter with positive constants $f_0$ and $\beta$. The global flux is given by

$$ \mathbf{G} = \left( \begin{array}{c} hu \\ hu^2 + \frac{g}{2}h^2 + \int_{\hat{x}}^{x} [ghB_x - fhv] \, d\xi \\ huv + \int_{\hat{x}}^{x} fhu \, d\xi \end{array} \right). $$

For other 1-D shallow water models, such as the thermal rotating shallow water equations, the shallow water flows in an open channel, and the blood flow models in veins and arteries, we believe that the proposed numerical schemes are extendable, and we will pursue these extensions in future research.

Finally, we point out that the high-order PAMPA scheme is only linear stable under certain CFL (Courant-Friedrichs-Lewy) constraint [1]. When the solution develops discontinuities, the high-order PAMPA scheme will be prone to numerical oscillations. Additionally, the high-order PAMPA scheme cannot guarantee the preservation of water depth positivity, which is crucial due to not only the physical reason but also the numerical reason (code may crash when water depth is negative). We, therefore, require the developed scheme should also be positive-preserving. This can be achieved by considering a convex combination of high-order and low-order fluxes/residuals. We provide PP WB low order schemes based on hydrostatic reconstruction [7]. The optimal combination coefficients are then derived following the idea introduced by bound-preserving PAMPA method developed in [2]. We also impose different conditions on the blending coefficients for the concern of spurious oscillations.

The main findings of this work can be summarized as follows:

* We propose a high-order fully WB PAMPA scheme via the global flux Gauss-Lobatto quadrature. We prove, both theoretically and numerically, that the proposed scheme verifies a nodal super convergence property for moving water steady states. Compared to some existing fully WB flux globalization based finite-volume and discontinuous Galerkin methods, the proposed method avoids solving nonlinear algebraic equations (as required, for example, in [19, 21, 60]). We also show that the high-order WB PAMPA scheme is capable of exactly preserving the still water equilibria.

* We present first-order schemes that preserve simultaneously both still water equilibria and positivity of water height. Positivity is guaranteed by the use of the local Lax-Friedrichs' flavor fluxes and residuals, while still water equilibria are maintained through the hydrostatic reconstructions used in [7, 54, 55]. Numerical tests on perturbed still water steady states confirm these properties and further highlight the superior performance of the high-order schemes in capturing small perturbation dynamics. The proposed PP WB first-order schemes improve upon the purely PP first-order schemes used in our previous works [40, 4].

* We introduce new convex blending methods of the global fluxes and residuals computed by the high-order PAMPA and low-order local Lax-Friedrichs' flavor schemes, following the idea in [2]. By examining the intermediate solution states, we establish a sufficient condition for the PP property and derive the optimal blending coefficients. To reduce oscillations near strong discontinuities, when the solution is away from a local discrete steady state, we incorporate an oscillating-elimination coefficient based on the damping term given in [45], which is then combined with the PP blending coefficient. This approach differs from our previous WB PAMPA methods in [40, 4], where a
posteriori MOOD paradigm is used.

The rest of the paper is organized as follows. In section 2, we introduce the novel flux globalization based PP WB PAMPA scheme. We provide some theoretical analysis of the proposed schemes in section 3, which includes the positivity-preserving and still water equilibria maintaining for the first-order scheme, as well as the positivity and still- and moving- water steady states preserving for the high-order scheme. In section 4, extensive numerical examples are presented to illustrate the high-order accuracy in smooth regions for general solutions, the positivity-preserving and the well-balanced property over non-flat bottom topography, and the robustness in capturing shocks over flat and non-flat bottom topography, of the proposed schemes.

**2. Numerical scheme.** In this section, we describe the flux globalization based PP WB PAMPA scheme for the general quasi-conservative system (1.2). To this end, we consider a tessellation of the 1-D spatial domain $\Omega$ in non overlapping elements $K_{j+\frac{1}{2}} = [x_j, x_{j+1}]$ with uniform size $\Delta x = x_{j+1} - x_j$. For reference, we list other commonly used notations and variables in Table 1.

<table><caption>TABLE 1
Notations and variables to be used in numerical methods.</caption><tr><td>t</td><td>The current time level</td></tr><tr><td>Δt</td><td>The adaptive time step determined by the CFL condition</td></tr><tr><td>Δx</td><td>The spatial grid size</td></tr><tr><td>λ</td><td>The time-space ratio Δt/Δx</td></tr><tr><td>u&#773;<sub>j+1/2</sub></td><td>The cell averages of the numerical solution in K<sub>j+1/2</sub> at time t</td></tr><tr><td>u<sub>j</sub></td><td>The point values of the numerical solution at x<sub>j</sub> and time t</td></tr><tr><td>u&#773;<sub>j+1/2</sub><sup>t+Δt</sup></td><td>The cell averages of the numerical solution in K<sub>j+1/2</sub> at time t + Δt</td></tr><tr><td>u<sub>j</sub><sup>t+Δt</sup></td><td>The point values of the numerical solution at x<sub>j</sub> and time t + Δt</td></tr><tr><td>[·]l,k</td><td>The difference of a function/variable evaluated at two points: (·)l - (·)k</td></tr><tr><td>{·}l,k</td><td>The average of a function/variable evaluated at two points: (·)l + (·)k/2</td></tr></table>

We now assume that the solution given in terms of average values $\bar{u}_{j+\frac{1}{2}}$ and point values $u_{j,j+1}$ for each element $K_{j+\frac{1}{2}}$ are available at the certain time level $t$. On each element $K_{j+\frac{1}{2}}$, we consider a global continuous finite element approximation $u_h(x)^1$ spanned by the quadratic polynomial basis functions:

$$ (2.1) \qquad u_h(x) = \varphi_j(x)u_j + \varphi_{j+1}(x)u_{j+1} + \bar{\varphi}_{j+\frac{1}{2}}(x)\bar{u}_{j+\frac{1}{2}}, \quad x \in K_{j+\frac{1}{2}}, $$

where $u_j \approx u(x_j)$ is the point value of $u$ at the element boundary $x_j$ and $\bar{u}_{j+\frac{1}{2}}$ is the cell average defined as

$$ \bar{u}_{j+\frac{1}{2}} = \frac{1}{|K_{j+\frac{1}{2}}|} \int_{K_{j+\frac{1}{2}}} u_h(x) \, dx. $$

Based on the following three conditions of the finite element approximation space:

<footnote n="1">¹Following the finite element convention, we use the subscript h to indicate a finite element approximation of a variable. For example, $u_h$ denotes a finite element approximation of $u$. Meanwhile, we have also used the letter h for the water depth in this paper. Therefore, $h_h$ would stand for a finite element approximation of the water depth h. Moreover, we have omitted the dependence of $u_h$ and other quantities on t for the sake of brevity from here on.</footnote>
• conservation condition:

$$ \int_{K_{j+\frac{1}{2}}} \bar{\varphi}_{j+\frac{1}{2}}(x) \, dx = |K_{j+\frac{1}{2}}|, \quad \int_{K_{j+\frac{1}{2}}} \varphi_j(x) \, dx = 0, $$

• interpolation condition:

$$ \bar{\varphi}_{j+\frac{1}{2}}(x_j) = 0, \quad \varphi_j(x_i) = \delta_{ij}, $$

• accuracy condition: $\varphi_j(x)$ and $\bar{\varphi}_{j+\frac{1}{2}}(x)$ are quadratic polynomials, we can easily obtain the explicit form of the basis functions ([1, 4]):

$$ \varphi_j = (1 - \xi)(1 - 3\xi), \quad \bar{\varphi}_{j+\frac{1}{2}} = 6\xi(1 - \xi), \quad \varphi_{j+1} = \xi(3\xi - 2), $$

where $\xi = \frac{x-x_j}{\Delta x} \in [0, 1]$.

The next question is to update the DoFs $\bar{u}_{j+\frac{1}{2}}$ and $u_j$ in a PP and WB manner, which can be done by solving the following two semi-discrete forms:

$$ (2.2a) \qquad \frac{d}{dt} \bar{u}_{j+\frac{1}{2}} = - \frac{\hat{G}_{j+1} - \hat{G}_j}{\Delta x} $$

and

$$ (2.2b) \qquad \frac{d}{dt} u_j = -(\hat{\Phi}_{j-\frac{1}{2}} + \hat{\Phi}_{j+\frac{1}{2}}). $$

**2.1. PP WB update for cell averages.** For each element $K_{j+\frac{1}{2}}$, we can define the PP WB blended fluxes $\hat{G}_j$ and $\hat{G}_{j+1}$ in (2.2a) as:

$$ (2.3) \qquad \begin{aligned} \hat{G}_j &= \tilde{G}_j(\bar{u}_{j+\frac{1}{2}}, \bar{u}_{j-\frac{1}{2}}) + \theta_j(G_j - \tilde{G}_j) := \tilde{G}_j + \theta_j \Delta G_j, \\ \hat{G}_{j+1} &= \tilde{G}_{j+1}(\bar{u}_{j+\frac{3}{2}}, \bar{u}_{j-\frac{1}{2}}) + \theta_{j+1}(G_{j+1} - \tilde{G}_{j+1}) := \tilde{G}_{j+1} + \theta_{j+1} \Delta G_{j+1}, \end{aligned} $$

where $G_j$ and $G_{j+1}$ are the high-order fluxes based on the Gauss–Lobatto quadrature rule, $\theta_j, \theta_{j+1} \in [0, 1]$ are the blending coefficients, $\tilde{G}_j = \tilde{G}_j(\bar{u}_{j+\frac{1}{2}}, \bar{u}_{j-\frac{1}{2}})$ and $\tilde{G}_{j+1} = \tilde{G}_{j+1}(\bar{u}_{j+\frac{3}{2}}, \bar{u}_{j-\frac{1}{2}})$ are the first-order local Lax–Friedrichs numerical fluxes, on the left and right interfaces, respectively.

**2.1.1. Low order fluxes.** We begin with constructing the first-order numerical fluxes, which is positivity-preserving, still water state preserving, and wet-dry fronts capturing. We proceed as follows. First we write

$$ (2.4) \qquad \begin{aligned} \tilde{G}_j &= \mathcal{F}^{\text{LLF}}(u_j^+, u_j^-) - \mathbf{R}_{j+\frac{1}{2}} + \int_{x_j}^{x_{j+\frac{1}{2}}} \mathbf{S} \, dx, \\ \tilde{G}_{j+1} &= \mathcal{F}^{\text{LLF}}(u_{j+1}^+, u_{j+1}^-) - \mathbf{R}_{j+\frac{1}{2}} - \int_{x_{j+\frac{1}{2}}}^{x_{j+1}} \mathbf{S} \, dx, \end{aligned} $$
with $\mathcal{F}^{\text{LLF}}$ denoting the local Lax-Friedrich's flux. Then, to compute the left and right states we use the hydrostatic reconstruction idea [7]. In particular, we define:

$$h_j^+ = \bar{h}_{j+\frac{1}{2}} + \bar{B}_{j+\frac{1}{2}} - B_j^+, \quad h_j^- = \bar{h}_{j-\frac{1}{2}} + \bar{B}_{j-\frac{1}{2}} - B_j^-,$$

$$u_j^+ = \frac{(hu)_{j+\frac{1}{2}}}{\bar{h}_{j+\frac{1}{2}}}, \quad u_j^- = \frac{(hu)_{j-\frac{1}{2}}}{\bar{h}_{j-\frac{1}{2}}}$$

$$v_j^+ = \frac{(hv)_{j+\frac{1}{2}}}{\bar{h}_{j+\frac{1}{2}}}, \quad v_j^- = \frac{(hv)_{j-\frac{1}{2}}}{\bar{h}_{j-\frac{1}{2}}}$$

with modified reconstructions of the topography at cell interfaces given by

$$B_j^+ = \min(\bar{h}_{j+\frac{1}{2}} + \bar{B}_{j+\frac{1}{2}}, \max(\bar{B}_{j+\frac{1}{2}}, \bar{B}_{j-\frac{1}{2}})),$$

$$B_j^- = \min(\bar{h}_{j-\frac{1}{2}} + \bar{B}_{j-\frac{1}{2}}, \max(\bar{B}_{j-\frac{1}{2}}, \bar{B}_{j-\frac{1}{2}})).$$

For the Saint-Venant system (1.5), we end with the left and right low order fluxes as:

$$\begin{equation} \begin{split} \tilde{\mathbf{G}}_j = \mathcal{F}^{\text{LLF}}(\mathbf{u}_j^+, \mathbf{u}_j^-) & - g \left( \left\{h\right\}_{j+\frac{1}{2},j^+} \left[\mathbf{B}\right]_{j+\frac{1}{2},j^+} \right) \\ & - \frac{gn^2}{2} \Delta x \left( \left\{ |hu|hu/h^{\frac{7}{3}} \right\}_{j+\frac{1}{2},j^+} \right) - \mathbf{R}_{j+\frac{1}{2}} \end{split} \end{equation}$$

and

$$\begin{equation} \begin{split} \tilde{\mathbf{G}}_{j+1} = \mathcal{F}^{\text{LLF}}(\mathbf{u}_{j+1}^+, \mathbf{u}_{j+1}^-) &+ g \left( \left\{h\right\}_{j+1^-, j+\frac{1}{2}} \left[\mathbf{B}\right]_{j+1^-, j+\frac{1}{2}} \right) \\ &+ \frac{gn^2}{2} \Delta x \left( \left\{ |hu|hu/h^{\frac{7}{3}} \right\}_{j+1^-, j+\frac{1}{2}} \right) - \mathbf{R}_{j+\frac{1}{2}}. \end{split} \end{equation}$$

For the rotating shallow water equations (1.6), we can define the left and right low order fluxes as:

$$\begin{equation} \tilde{\mathbf{G}}_j = \mathcal{F}^{\text{LLF}}(\mathbf{u}_j^+, \mathbf{u}_j^-) - \left( g \left\{h\right\}_{j+\frac{1}{2},j^+} \left[\mathbf{B}\right]_{j+\frac{1}{2},j^+} - \frac{\Delta x}{2} \left\{fhv\right\}_{j+\frac{1}{2},j^+} - \frac{\Delta x}{2} \left\{fhu\right\}_{j+\frac{1}{2},j^+} \right) - \mathbf{R}_{j+\frac{1}{2}}, \end{equation}$$

and

$$\begin{equation} \tilde{\mathbf{G}}_{j+1} = \mathcal{F}^{\text{LLF}}(\mathbf{u}_{j+1}^+, \mathbf{u}_{j+1}^-) + \left( g \left\{h\right\}_{j+1^-, j+\frac{1}{2}} \left[\mathbf{B}\right]_{j+1^-, j+\frac{1}{2}} - \frac{\Delta x}{2} \left\{fhv\right\}_{j+1^-, j+\frac{1}{2}} - \frac{\Delta x}{2} \left\{fhu\right\}_{j+1^-, j+\frac{1}{2}} \right) - \mathbf{R}_{j+\frac{1}{2}}. \end{equation}$$

### 2.1.2. High order fluxes.
Thanks to the global continuous approximation in the PAMPA framework, we do not need to use any Riemann solver and directly compute the high order fluxes at the cell interfaces instead:

$$\begin{align} \mathbf{G}_j &= \mathbf{f}(\mathbf{u}_j) - \mathbf{R}_j, & \mathbf{G}_{j+1} &= \mathbf{f}(\mathbf{u}_{j+1}) - \mathbf{R}_{j+1} \notag \\ \mathbf{R}_j &= \int_{\hat{x}}^{x_j} \mathbf{S}(\mathbf{u}_h, x) \, dx & \mathbf{R}_{j+1} &= \int_{\hat{x}}^{x_{j+1}} \mathbf{S}(\mathbf{u}_h, x) \, dx. \end{align}$$
In order to complete the computations in (2.4) and (2.9), we need to compute the source primitive terms {$\mathbf{R}_j$} and {$\mathbf{R}_{j+\frac{1}{2}}$}. This can be done in several ways. In [42] the authors have used a nodal Lagrange expansion of the source term on the three Gauss-Lobatto points {$x_j, x_{j+1/2}, x_{j+1}$}. This leads to

$$ (2.10) \qquad \begin{aligned} \mathbf{R}_{j+1} &= \mathbf{R}_j + \int_{x_j}^{x_{j+1}} \mathbf{S}_h(\mathbf{u}_h, x) \, dx = \mathbf{R}_j + \sum_{i=1}^3 \omega_i \mathbf{S}(\mathbf{u}_h(x_i), x_i) \\ \mathbf{R}_{j+\frac{1}{2}} &= \mathbf{R}_j + \int_{x_j}^{x_{j+\frac{1}{2}}} \mathbf{S}_h(\mathbf{u}_h, x) \, dx = \mathbf{R}_j + \sum_{i=1}^3 \omega'_i \mathbf{S}(\mathbf{u}_h(x_i), x_i), \quad \forall j \end{aligned} $$

with $\omega_i$ and $\omega'_i$ obtained integrating the Lagrange polynomials on the interval and half-interval respectively. This formulation can be written in compact form as

$$ (2.11) \qquad \left( \begin{array}{c} \mathbf{R}_j \\ \mathbf{R}_{j+\frac{1}{2}} \\ \mathbf{R}_{j+1} \end{array} \right) = \mathbf{R}_j + I_{\text{LIIA}} \left( \begin{array}{c} \mathbf{S}_j \\ \mathbf{S}_{j+\frac{1}{2}} \\ \mathbf{S}_{j+1} \end{array} \right) $$

where, as noted in [42], by construction the matrix

$$ (2.12) \qquad I_{\text{LIIA}} = \begin{bmatrix} 0 & 0 & 0 \\ \omega'_1 & \omega'_2 & \omega'_3 \\ \omega_1 & \omega_2 & \omega_3 \end{bmatrix} = \Delta x \begin{bmatrix} 0 & 0 & 0 \\ 5/24 & 1/3 & -1/24 \\ 1/6 & 2/3 & 1/6 \end{bmatrix} $$

is the Butcher table of the well-known fully implicit Runge-Kutta collocation method LobattoIIIA [32]. For autonomous problems, this method can be shown to have a convergence order $p+2$ for the midpoint, and $2p$ for the end-points (here $p=2$ is the polynomial degree in the finite element approximation space). These estimates have been confirmed on the steady ODE in [42], and lead in our case to a fourth order method. Despite of this, note that the second quadrature rule associated to the second line of (2.11) (or equivalently the second in (2.10)) is not symmetric and only exact for quadratic polynomials.

A different approach is to apply the usual three-points Gauss-Lobatto formula within both the half, and full intervals. In this case, both formulas are exact for cubic polynomials, and we obtain a sub-cell LobattoIII method which can be written as

$$ (2.13) \qquad \left( \begin{array}{c} \mathbf{R}_j \\ \mathbf{R}_{j+\frac{1}{2}} \\ \mathbf{R}_{j+1} \end{array} \right) = \mathbf{R}_j + I_{\text{sc-LIII}} \left( \begin{array}{c} \mathbf{S}(\mathbf{u}_j, x_j) \\ \mathbf{S}(\mathbf{u}_{j+1/4}, x_{j+1/4}) \\ \mathbf{S}(\mathbf{u}_{j+\frac{1}{2}}, x_{j+\frac{1}{2}}) \\ \mathbf{S}(\mathbf{u}_{j+1}, x_{j+1}) \end{array} \right), $$

and

$$ (2.14) \qquad I_{\text{sc-LIII}} = \Delta x \begin{bmatrix} 0 & 0 & 0 & 0 \\ 1/12 & 1/3 & 1/12 & 0 \\ 1/6 & 0 & 2/3 & 1/6 \end{bmatrix}, \quad \mathbf{u}_{j+1/4} = \frac{3}{16} \mathbf{u}_j + \frac{9}{8} \bar{\mathbf{u}}_{j+\frac{1}{2}} - \frac{5}{16} \mathbf{u}_{j+1}, $$

where the last expression is obtained from quadratic interpolation within the cell (2.1). We can show the following result.

**PROPOSITION 2.1** (sc-LIII consistency). *The implicit sub-cell LobattoIII RK method defined by (2.13)-(2.14) is fourth order accurate for smooth enough solutions.*

*Proof.* The result can be easily shown by truncated Taylor series developments. Details are reported in Appendix A for completeness. □
In both cases, the computation of the nodal values of $\mathbf{R}$ can be seen as the integration of the steady-state ODE with the LobattoIIIA or sc-LobattoIII methods. A very important remark is that in practice only differences of values of $\mathbf{R}$ appear in our scheme. This, and the continuity of the approximation, make the method independent on the initial quadrature point. As a consequence, despite the name commonly used in literature, this approach leads to a fully local method. If necessary to prepare the solution, one can initiate the computation of these terms from $\hat{x} = x_0$ (left-end of the domain), and set $\mathbf{R}_0 = 0$. However, as already observed, this initial value plays no role in practice.

The two methods have similar properties. For the shallow water models, the sc-LobattoIII approach has some advantage as it allows to integrate exactly the hydrostatic term. For this reason we will use this approach in the following.

**2.2. PP WB update for point values.** Analogously to the end-point fluxes, the PP blended residuals are defined as

$$ (2.15) \qquad \begin{aligned} \tilde{\mathbf{\Phi}}_{j-\frac{1}{2}} &= \tilde{\mathbf{\Phi}}_{j-\frac{1}{2}} + \theta_{j-\frac{1}{2}}(\tilde{\mathbf{\Phi}}_{j-\frac{1}{2}} - \tilde{\mathbf{\Phi}}_{j-\frac{1}{2}}) := \tilde{\mathbf{\Phi}}_{j-\frac{1}{2}} + \theta_{j-\frac{1}{2}} \Delta \tilde{\mathbf{\Phi}}_{j-\frac{1}{2}}, \\ \tilde{\mathbf{\Phi}}_{j+\frac{1}{2}} &= \tilde{\mathbf{\Phi}}_{j+\frac{1}{2}} + \theta_{j+\frac{1}{2}}(\tilde{\mathbf{\Phi}}_{j+\frac{1}{2}} - \tilde{\mathbf{\Phi}}_{j+\frac{1}{2}}) := \tilde{\mathbf{\Phi}}_{j+\frac{1}{2}} + \theta_{j+\frac{1}{2}} \Delta \tilde{\mathbf{\Phi}}_{j+\frac{1}{2}}. \end{aligned} $$

**2.2.1. Low order residuals.** Mimicking the low order fluxes based on the hydrostatic reconstruction, we construct low order point fluctuations. First, we define the following quantities

$$ \begin{aligned} B_{j+\frac{1}{4}}^{+} &= \min(\bar{h}_{j+\frac{1}{2}} + \bar{B}_{j+\frac{1}{2}}, \max(\bar{B}_{j+\frac{1}{2}}, B_j)), & B_{j+\frac{1}{4}}^{-} &= \min(h_j + B_j, \max(\bar{B}_{j+\frac{1}{2}}, B_j)), \\ h_{j+\frac{1}{4}}^{+} &= \bar{h}_{j+\frac{1}{2}} + \bar{B}_{j+\frac{1}{2}} - B_{j+\frac{1}{4}}^{+}, & h_{j+\frac{1}{4}}^{-} &= h_j + B_j - B_{j+\frac{1}{4}}^{-}, \\ u_{j+\frac{1}{4}}^{+} &= \frac{(hu)_{j+\frac{1}{2}}}{\bar{h}_{j+\frac{1}{2}}}, & u_{j+\frac{1}{4}}^{-} &= u_j, & v_{j+\frac{1}{4}}^{+} &= \frac{(hv)_{j+\frac{1}{2}}}{\bar{h}_{j+\frac{1}{2}}}, & v_{j+\frac{1}{4}}^{-} &= v_j. \end{aligned} $$

Then, the low order residuals are given by

$$ (2.16) \qquad \begin{aligned} \frac{\Delta x}{2} \tilde{\mathbf{\Phi}}_{j-\frac{1}{2}} &= \mathbf{f}(\mathbf{u}_j) - \mathcal{F}^{\text{LLF}}(\mathbf{u}_{j-\frac{1}{4}+}, \mathbf{u}_{j-\frac{1}{4}-}) + \left( g \{h\}_{j,j-\frac{1}{4}+} \llbracket B \rrbracket_{j,j-\frac{1}{4}+} \right) \\ &\quad + \frac{gn^2}{4} \Delta x \left( \{ |hu|hu/h^{\frac{7}{3}} \}_{j,j-\frac{1}{4}+} \right) \end{aligned} $$

and

$$ (2.17) \qquad \begin{aligned} \frac{\Delta x}{2} \tilde{\mathbf{\Phi}}_{j+\frac{1}{2}} &= \mathcal{F}^{\text{LLF}}(\mathbf{u}_{j+\frac{1}{4}+}, \mathbf{u}_{j+\frac{1}{4}-}) - \mathbf{f}(\mathbf{u}_j) + \left( g \{h\}_{j+\frac{1}{4}-,j} \llbracket B \rrbracket_{j+\frac{1}{4}-,j} \right) \\ &\quad + \frac{gn^2}{4} \Delta x \left( \{ |hu|hu/h^{\frac{7}{3}} \}_{j+\frac{1}{4}-,j} \right) \end{aligned} $$

For the rotating shallow water equations (1.6), the low order residuals are given
by

$$ (2.18) \qquad \begin{aligned} \frac{\Delta x}{2} \tilde{\mathbf{\Phi}}_{j-\frac{1}{2}} = \mathbf{f}(\mathbf{u}_j) - \mathcal{F}^{\text{LLF}}(\mathbf{u}_{j-\frac{1}{4}+}, \mathbf{u}_{j-\frac{1}{4}-}) &+ \left( g \left\{ h \right\}_{j,j-\frac{1}{4}+}^{0} \left[ B \right]_{j,j-\frac{1}{4}+}^{0} \right) \\ & - \frac{\Delta x}{4} \left( \begin{array}{c} 0 \\ \left\{ f h v \right\}_{j,j-\frac{1}{4}+} \\ - \left\{ f h u \right\}_{j,j-\frac{1}{4}+} \end{array} \right) \end{aligned} $$

and

$$ (2.19) \qquad \begin{aligned} \frac{\Delta x}{2} \tilde{\mathbf{\Phi}}_{j+\frac{1}{2}} = \mathcal{F}^{\text{LLF}}(\mathbf{u}_{j+\frac{1}{4}}^{+}, \mathbf{u}_{j+\frac{1}{4}}^{-}) - \mathbf{f}(\mathbf{u}_j) &+ \left( g \left\{ h \right\}_{j+\frac{1}{4}-,j}^{0} \left[ B \right]_{j+\frac{1}{4}-,j}^{0} \right) \\ & - \frac{\Delta x}{4} \left( \begin{array}{c} 0 \\ \left\{ f h v \right\}_{j+\frac{1}{4}-,j} \\ - \left\{ f h u \right\}_{j+\frac{1}{4}-,j} \end{array} \right) \end{aligned} $$

2.2.2. **High order residuals.** The high-order residuals in (2.15) are given by

$$ (2.20) \qquad \tilde{\mathbf{\Phi}}_{j-\frac{1}{2}} = \tilde{\mathcal{J}}^{+}(\mathbf{u}_j) \partial_x^{+} \mathbf{G}, \quad \tilde{\mathbf{\Phi}}_{j+\frac{1}{2}} = \tilde{\mathcal{J}}^{-}(\mathbf{u}_j) \partial_x^{-} \mathbf{G}, $$

where $\partial_x^+ \mathbf{G}$ and $\partial_x^- \mathbf{G}$ are the left- and right-biased high-order FD approximations of $\mathbf{G}_x$ at $x = x_j$, respectively. They are given by

$$ (2.21) \qquad \partial_x^+ \mathbf{G} = \frac{1}{\Delta x} (\mathbf{G}_{j-1} - 4\mathbf{G}_{j-\frac{1}{2}} + 3\mathbf{G}_j), \quad \partial_x^- \mathbf{G} = \frac{1}{\Delta x} (-3\mathbf{G}_j + 4\mathbf{G}_{j+\frac{1}{2}} - \mathbf{G}_{j+1}), $$

where $\mathbf{G}_j$ and $\mathbf{G}_{j+1}$ for each element $K_{j+\frac{1}{2}}$ are computed in (2.9) and

$$ \mathbf{G}_{j+\frac{1}{2}} = \mathbf{f}(\mathbf{u}_{j+\frac{1}{2}}) - \mathbf{R}_{j+\frac{1}{2}} $$

with $\mathbf{u}_{j+\frac{1}{2}}$ being computed by the local scaling PP limiter:

$$ (2.22) \qquad \mathbf{u}_{j+\frac{1}{2}} = (1 - \eta)\bar{\mathbf{u}}_{j+\frac{1}{2}} + \eta\hat{\mathbf{u}}_{j+\frac{1}{2}}, \quad \eta \in [0, 1], $$

where $\hat{\mathbf{u}}_{j+\frac{1}{2}} = \frac{3}{2}\bar{\mathbf{u}}_{j+\frac{1}{2}} - \frac{1}{4}(\mathbf{u}_j + \mathbf{u}_{j+1})$ obtained from (2.1) and $\eta$ is given by

$$ \eta = \min \left( 1, \frac{\bar{h}_{j+\frac{1}{2}} - \varepsilon_h}{\bar{h}_{j+\frac{1}{2}} - \hat{h}_{j+\frac{1}{2}}} \right), \quad \varepsilon_h := \min_j \left( 10^{-13}, \bar{h}_{j+\frac{1}{2}} \right). $$

Note that $\hat{\mathbf{u}}_{j+\frac{1}{2}}$ does not guarantee the positivity of water depth and thus the PP limiter is crucial. Similar strategy has also been applied in the computation of $\mathbf{u}_{j+\frac{1}{4}}$ in (2.14).

The sign matrices are defined as

$$ \tilde{\mathcal{J}}^{\pm}(\mathbf{u}_j) = J^{-1}J^{\pm} = \mathcal{R}\Lambda^{\pm}\Lambda^{-1}\mathcal{R}^{-1} = \mathcal{R} \begin{pmatrix} \frac{\lambda_1^{\pm}}{\lambda_1} & & & \\ & \ddots & & \\ & & \frac{\lambda_2^{\pm}}{\lambda_2} & \\ & & & \ddots \end{pmatrix} \mathcal{R}^{-1}, $$
where $\lambda_i$, $i = 1, 2, \dots$ are the eigenvalues of the Jacobian $J = \frac{\partial f}{\partial \mathbf{u}}$ and $\mathcal{R}$ is the corresponding eigenvectors matrix. $\lambda_i^+ = \max(\lambda_i, 0)$ and $\lambda_i^- = \min(\lambda_i, 0)$.

The idea behind defining the high-order residuals in (2.20) are based on the following two conditions:

* linear stable, which is achieved by considering the upwinding flavor and biased FD approximations of the spatial derivatives,

* consistency: if we have $(\partial_x \mathbf{G})_{x_j} = \text{Const}$, we are able to recover

$$ \vec{\mathbf{F}}_{j-\frac{1}{2}} + \vec{\mathbf{F}}_{j+\frac{1}{2}} = (\partial_x \mathbf{G})_{x_j}. $$

With these definitions we can easily show the following property.

**PROPOSITION 2.2 (Compactness).** *The positivity preserving PAMPA method with global flux quadrature is fully local, and has a compact stencil.*

*Proof.* To show the property, it is enough to check that everywhere we only use differences of the values of $\mathbf{R}$. We can start with the numerical fluxes for the update of the cell average. Combining (2.4) and (2.5)-(2.6) (or (2.7)-(2.8)) we can write

$$ \tilde{\mathbf{G}}_{j+1} - \tilde{\mathbf{G}}_j = \mathcal{F}^{\text{LLF}}(\mathbf{u}_{j+1}^+, \mathbf{u}_{j+1}^-) - \int_{x_{j+\frac{1}{2}}}^{x_{j+1}} \mathbf{S} \, dx - \mathcal{F}^{\text{LLF}}(\mathbf{u}_j^+, \mathbf{u}_j^-) - \int_{x_j}^{x_{j+\frac{1}{2}}} \mathbf{S} \, dx $$

where the integrals are replaced by the approximations in (2.5)-(2.6) (or (2.7)-(2.8)). This expression is independent of $\mathbf{R}$. Similarly, we can see that in (2.3)

$$ \Delta \mathbf{G}_j = \mathbf{f}_j - \mathbf{R}_j + \mathbf{R}_{j+\frac{1}{2}} - \mathcal{F}^{\text{LLF}}(\mathbf{u}_j^+, \mathbf{u}_j^-) - \int_{x_j}^{x_{j+\frac{1}{2}}} \mathbf{S} \, dx, $$

$$ \Delta \mathbf{G}_{j+1} = \mathbf{f}_{j+1} - \mathbf{R}_{j+1} + \mathbf{R}_{j+\frac{1}{2}} - \mathcal{F}^{\text{LLF}}(\mathbf{u}_{j+1}^+, \mathbf{u}_{j+1}^-) + \int_{x_{j+\frac{1}{2}}}^{x_{j+1}} \mathbf{S} \, dx, $$

where the integrals are again replaced by the approximations in (2.5)-(2.6) (or (2.7)-(2.8)) and where by definition $\mathbf{R}_{j+\frac{1}{2}} - \mathbf{R}_j$ is obtained from the second in (2.10), or from (2.13)-(2.14), which both only involve computations within the cell (similar to $\mathbf{R}_{j+1} - \mathbf{R}_{j+\frac{1}{2}} = \mathbf{R}_{j+1} - \mathbf{R}_j + \mathbf{R}_j - \mathbf{R}_{j+\frac{1}{2}}$). This shows that the limited update for the cell average is local.

Concerning the nodal updates, we just need to show that the residuals are local. We can recast (2.21) as

$$ \partial_x^+ \mathbf{G} = \frac{4}{\Delta x}(\mathbf{G}_j - \mathbf{G}_{j-\frac{1}{2}}) - \frac{1}{\Delta x}(\mathbf{G}_j - \mathbf{G}_{j-1}). $$

And similarly for $\partial_x^- \mathbf{G}$ (omitted for brevity). For the source part this only involves the differences

$$ \frac{4}{\Delta x}(-\mathbf{R}_j + \mathbf{R}_{j-\frac{1}{2}}) - \frac{1}{\Delta x}(-\mathbf{R}_j + \mathbf{R}_{j-1}) = -\frac{3}{\Delta x}(\mathbf{R}_j - \mathbf{R}_{j-1}) + \frac{4}{\Delta x}(\mathbf{R}_{j-\frac{1}{2}} - \mathbf{R}_{j-1}) $$

which again can be obtained from local evaluations given by (2.11)-(2.12), or (2.13)-(2.14). So the high order fluctuations are also fully local. □

**2.3. Convex combination coefficients.** Here, we derive the convex combination coefficients that used for the blended fluxes and residuals in (2.3) and (2.15) respectively, to close this section. As we mentioned in section 1, the proposed PP
WB PAMPA scheme will be applied for solving the Saint–Venant system (1.5) and the rotating shallow water equations (1.6). Consequently, we only need to ensure the PP property for the water depth h.

We apply the low order schemes introduced in subsection 2.1.1 and consider the forward Euler time stepping for the average value $\bar{h}_{j+\frac{1}{2}}$ from t to t + $\Delta t$:

$$ 
\begin{aligned}
\bar{h}_{j+\frac{1}{2}}^{t+\Delta t} &= \bar{h}_{j+\frac{1}{2}} - \lambda (\hat{\mathbf{G}}_{j+1}^{(1)} - \hat{\mathbf{G}}_{j}^{(1)}) \pm \frac{\lambda}{2}(\alpha_{j+1} + \bar{u}_{j+\frac{1}{2}})h_{j+1}^{-} \pm \frac{\lambda}{2}(\alpha_{j} - \bar{u}_{j+\frac{1}{2}})h_{j}^{+} \\
&= \bar{h}_{j+\frac{1}{2}} - \frac{\lambda}{2}(\alpha_{j+1} + \bar{u}_{j+\frac{1}{2}})h_{j+1}^{-} - \frac{\lambda}{2}(\alpha_{j} - \bar{u}_{j+\frac{1}{2}})h_{j}^{+} \\
&\quad + \lambda\alpha_{j}\left(\frac{1}{2}\left(1 - \frac{\bar{u}_{j+\frac{1}{2}}}{\alpha_{j}}\right)h_{j}^{+} + \frac{\hat{\mathbf{G}}_{j}^{(1)}}{\alpha_{j}}\right) + \lambda\alpha_{j+1}\left(\frac{1}{2}\left(1 + \frac{\bar{u}_{j+\frac{1}{2}}}{\alpha_{j+1}}\right)h_{j+1}^{-} - \frac{\hat{\mathbf{G}}_{j+1}^{(1)}}{\alpha_{j+1}}\right) \\
&= \tilde{h}_{j+\frac{1}{2}}^{C} + \lambda\alpha_{j}\tilde{h}_{j+\frac{1}{2}}^{L} + \lambda\alpha_{j+1}\tilde{h}_{j+\frac{1}{2}}^{R},
\end{aligned}
$$

where

$$ 
\begin{aligned}
\tilde{h}_{j+\frac{1}{2}}^{C} &= \bar{h}_{j+\frac{1}{2}} - \frac{\lambda}{2}(\alpha_{j+1} + \bar{u}_{j+\frac{1}{2}})h_{j+1}^{-} - \frac{\lambda}{2}(\alpha_{j} - \bar{u}_{j+\frac{1}{2}})h_{j}^{+} \\
&\geq \bar{h}_{j+\frac{1}{2}}\left(1 - \lambda\alpha_{j+1}\frac{h_{j+1}^{-}}{\bar{h}_{j+\frac{1}{2}}} - \lambda\alpha_{j}\frac{h_{j}^{+}}{\bar{h}_{j+\frac{1}{2}}}\right) \geq 0.
\end{aligned}
$$

The last inequality due to the fact that $0 \leq \bar{h}_{j+1}^{-} \leq \bar{h}_{j+\frac{1}{2}}$, $0 \leq h_{j}^{+} \leq \bar{h}_{j+\frac{1}{2}}$, and the CFL condition $\lambda\alpha \leq \frac{1}{2}$ with $\alpha = \max_j\{\alpha_j\}$. Moreover,

$$ \tilde{h}_{j+\frac{1}{2}}^{L} = \frac{1}{2}\left(1 - \frac{\bar{u}_{j+\frac{1}{2}}}{\alpha_j}\right)h_j^+ + \frac{\hat{\mathbf{G}}_j^{(1)}}{\alpha_j} = \frac{1}{2}h_j^-\left(1 + \frac{\bar{u}_{j-\frac{1}{2}}}{\alpha_j}\right) + \theta_j \frac{\Delta \mathbf{G}_j^{(1)}}{\alpha_j}, $$

$$ \tilde{h}_{j+\frac{1}{2}}^{R} = \frac{1}{2}\left(1 + \frac{\bar{u}_{j+\frac{1}{2}}}{\alpha_{j+1}}\right)h_{j+1}^{-} - \frac{\hat{\mathbf{G}}_{j+1}^{(1)}}{\alpha_{j+1}} = \frac{1}{2}h_{j+1}^{+}\left(1 - \frac{\bar{u}_{j+\frac{3}{2}}}{\alpha_{j+1}}\right) - \theta_{j+1}\frac{\Delta \mathbf{G}_{j+1}^{(1)}}{\alpha_{j+1}}. $$

Since $\frac{1}{2}h_j^{-}(1 + \frac{\bar{u}_{j-\frac{1}{2}}}{\alpha_j}) \geq 0$ and $\frac{1}{2}h_{j+1}^{+}(1 - \frac{\bar{u}_{j+\frac{3}{2}}}{\alpha_{j+1}}) \geq 0$, by requiring $\tilde{h}_{j+\frac{1}{2}}^{L} \geq 0$ and $\tilde{h}_{j+\frac{1}{2}}^{R} \geq 0$, we therefore obtain the value of blending coefficient, i.e.,

$$ (2.24) \qquad \theta_j = \min \left( 1, \frac{\alpha_j \frac{h_j^+}{2} \left(1 - \frac{\bar{u}_{j+\frac{1}{2}}}{\alpha_j}\right)}{|\Delta \mathbf{G}_j^{(1)}|}, \frac{\alpha_j \frac{h_j^-}{2} \left(1 + \frac{\bar{u}_{j-\frac{1}{2}}}{\alpha_j}\right)}{|\Delta \mathbf{G}_j^{(1)}|} \right). $$

Similarly, we can write down the forward Euler time stepping for the point values as

$$ 
\begin{aligned}
h_j^{t+\Delta t} &= h_j - 2\lambda(\hat{\mathbf{F}}_{j-\frac{1}{2}}^{(1)} + \hat{\mathbf{F}}_{j+\frac{1}{2}}^{(1)}) \pm \lambda(\alpha_{j+\frac{1}{4}} + u_j)h_{j+\frac{1}{4}}^{-} \pm \lambda(\alpha_{j-\frac{1}{4}} - u_j)h_{j-\frac{1}{4}}^{+} \\
&= h_j \left(1 - \lambda(\alpha_{j+\frac{1}{4}} + u_j)\frac{h_{j+\frac{1}{4}}^{-}}{h_j} - \lambda(\alpha_{j-\frac{1}{4}} - u_j)\frac{h_{j-\frac{1}{4}}^{+}}{h_j}\right) \\
&\quad + \lambda(\alpha_{j-\frac{1}{4}} + u_j)h_{j-\frac{1}{4}}^{+} - 2\lambda\hat{\mathbf{F}}_{j-\frac{1}{2}}^{(1)} + \lambda(\alpha_{j+\frac{1}{4}} + u_j)h_{j+\frac{1}{4}}^{-} - 2\lambda\hat{\mathbf{F}}_{j+\frac{1}{2}}^{(1)} \\
&= \tilde{h}_j^C + 2\lambda\alpha_{j-\frac{1}{4}}\tilde{h}_j^L + 2\lambda\alpha_{j+\frac{1}{4}}\tilde{h}_j^R.
\end{aligned}
$$
Since $0 \le h_{j-\frac{1}{4}}^+ \le h_j$, $0 \le h_{j+\frac{1}{4}}^- \le h_j$, and the CFL condition $\lambda\alpha \le \frac{1}{4}$, we immediately obtain $\tilde{h}_j^C \ge 0$. Then, we observe

$$ \tilde{h}_j^L = \frac{1}{2} \left(1 + \frac{u_j}{\alpha_{j-\frac{1}{4}}}\right) h_{j-\frac{1}{4}}^+ - \frac{\hat{\Phi}_{j-\frac{1}{2}}^{(1)}}{\alpha_{j-\frac{1}{4}}} = \frac{h_{j-\frac{1}{4}}^+}{2} \left(1 + \frac{\bar{u}_{j-\frac{1}{2}}}{\alpha_{j-\frac{1}{4}}}\right) - \theta_{j-\frac{1}{2}} \frac{\Delta \hat{\Phi}_{j-\frac{1}{2}}^{(1)}}{\alpha_{j-\frac{1}{4}}}, $$

$$ \tilde{h}_j^R = \frac{h_{j+\frac{1}{4}}^-}{2} \left(1 - \frac{u_j}{\alpha_{j+\frac{1}{4}}}\right) - \frac{\hat{\Phi}_{j+\frac{1}{2}}^{(1)}}{\alpha_{j+\frac{1}{4}}} = \frac{h_{j+\frac{1}{4}}^-}{2} \left(1 - \frac{\bar{u}_{j+\frac{1}{2}}}{\alpha_{j+\frac{1}{4}}}\right) - \theta_{j+\frac{1}{2}} \frac{\Delta \Phi_{j+\frac{1}{2}}^{(1)}}{\alpha_{j+\frac{1}{4}}}. $$

Knowing that $\frac{1}{2}(1 + \frac{\bar{u}_{j-\frac{1}{2}}}{\alpha_{j-\frac{1}{4}}})h_{j-\frac{1}{4}}^+ \ge 0$ and $\frac{1}{2}(1 - \frac{\bar{u}_{j+\frac{1}{2}}}{\alpha_{j+\frac{1}{4}}})h_{j+\frac{1}{4}}^- \ge 0$, we can ask $\tilde{h}_j^L \ge 0$ and $\tilde{h}_j^R \ge 0$. Thus, we obtain (2.25)

$$ \theta_{j-\frac{1}{2}} = \min \left( 1, \alpha_{j-\frac{1}{4}} \frac{\frac{h_{j-\frac{1}{4}}^+}{2} \left( 1 + \frac{\bar{u}_{j-\frac{1}{4}}}{\alpha_{j-\frac{1}{4}}} \right)}{|\Delta \vec{\Phi}_{j-\frac{1}{2}}^{(1)}|} \right), \quad \theta_{j+\frac{1}{2}} = \min \left( 1, \alpha_{j+\frac{1}{4}} \frac{\frac{h_{j+\frac{1}{4}}^-}{2} \left( 1 - \frac{\bar{u}_{j+\frac{1}{4}}}{\alpha_{j+\frac{1}{4}}} \right)}{|\Delta \vec{\Phi}_{j+\frac{1}{2}}^{(1)}|} \right) $$

for any $K_{j+\frac{1}{2}}$.

REMARK 2.3. There are several versions of hydrostatic reconstruction, see, e.g., [16]. We think that many of them can be applied in the proposed monolithic convex limiting framework and the combination coefficients derived in (2.24) and (2.25) will keep the same form.

### 2.3.1. Suppress spurious oscillations.
In order to suppress spurious oscillations as much as possible when approximating discontinuous solutions, we revisit the combination weights determined by the high-order derivative jumps at interfaces. The blending coefficient is designed as (inspired by [3])

$$ (2.26) \qquad \theta_{j+\frac{1}{2}}^{\text{OE}} = \exp \left( - \frac{\alpha_{j+\frac{1}{2}} \Delta t_n \sigma_{j+\frac{1}{2}}(h_h)}{\Delta x} \right), $$

where $\Delta t_n$ is the n-th adaptive time step computed using a suitable CFL condition, $\alpha_{j+\frac{1}{2}}$ is the estimate of the local maximum wave speed around $K_{j+\frac{1}{2}}$. The key of the method is $\sigma_{j+\frac{1}{2}}$, which here is computed based on the solution of water depth, and plays a critical role in balancing accuracy and stability. Ideally, $\sigma_{j+\frac{1}{2}}$ remains very small in smooth regions to preserve high-order accuracy, and becomes sufficiently large near discontinuities to suppress spurious oscillations. Getting inspiration from [45], we define it as

$$ \sigma_{j+\frac{1}{2}}(h_h) = \sum_{m=0}^{2} \sigma_{j+\frac{1}{2}}^{m}(h_h) $$

with

$$ \sigma_{j+\frac{1}{2}}^{m} = \begin{cases} 0, & \text{if } h_h \equiv \langle h_h \rangle_{\Omega}, \\ \frac{\Delta x^m \left| \left| \partial_x^m h_h \right| \right|_j + \left| \left| \partial_x^m h_h \right| \right|_{j+1}}{2 \| h_h - \langle h_h \rangle_{\Omega} \|_{L^{\infty}(\Omega)}} & \text{otherwise,} \end{cases} $$

where $\llbracket \cdot \rrbracket_j$ depicts the jump quantity at $x = x_j$ and $\langle h_h \rangle_{\Omega} = \frac{1}{|\Omega|} \int_{\Omega} h_h dx$ denotes the global average of $h_h(x)$ over the entire computational domain $\Omega$.
As currently formulated, this correction would alter the fully well-balanced property. For this reason, in the numerical implementation we only activate this oscillation eliminating procedure when the numerical solution is far from equilibrium. In particular, we say the solution is near/at a local steady state around cell $K_{j+\frac{1}{2}}$ if $H(\phi_{j+\frac{1}{2}}) \le 10^{-3}$, where

$$ (2.27) \qquad H(\phi_{j+\frac{1}{2}}) = \frac{(C\phi_{j+\frac{1}{2}})^\kappa}{1 + (C\phi_{j+\frac{1}{2}})^\kappa} $$

with the constants $C = 10$ and $\kappa = 20$ being used in all of the numerical experiments, and

$$ (2.28) \qquad \phi_{j+\frac{1}{2}} = \langle \partial \mathbf{G}_x^{(2)} \rangle_{K_{j+\frac{1}{2}}} \cdot \frac{|\Omega|}{\max\{\mathbf{G}_j^{(2)}, \mathbf{G}_{j+\frac{1}{2}}^{(2)}, \mathbf{G}_{j+1}^{(2)}\}} $$

We plot a sketch of this function in Figure 1, where one can clearly see that if $\phi$ is small (near/at a steady state), then the value of $H$ is still close to 0. When $\phi$ increases, the values of $\phi$ rapidly approaches 1.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Sketch of the cut-off function $H(\phi)$.</figcaption></figure>

In the end, the combination coefficients for the fluxes and residuals are taken as the minimal value between (2.26) and (2.24), (2.25), respectively.

REMARK 2.4. In the near dry regions (water height $h$ closes to zero), the potentially harmful division by zero, for example when we need to compute the velocity $u = \frac{hu}{h}$, should be controlled. We achieve this following the idea in [26]: if $h \le 10^{-14}$, we set $u = 0$. Otherwise, if $h > h_0 = 10^{-4}$, the velocity is computed normally as $u = \frac{hu}{h}$. If instead $h < h_0$, then the velocity is computed by a filtered division as

$$ u = hu \frac{h}{h^2 + f(h)\varepsilon}, $$

where $\varepsilon = 5 \times 10^{-9}$ and $f(h)$ is a filter function defined as

$$ f(h) = 2(h/h_0)^3 - 3(h/h_0)^2 + 1. $$

For more details, we refer to [26, Section 3.4].

3. Positivity-preserving and well-balanced properties. In this section, we provide some basic properties of the proposed flux globalization based WB PAMPA schemes.
**PROPOSITION 3.1** (Positivity-preserving of the first-order schemes). In the first-order numerical schemes with the numerical fluxes given by (2.5)-(2.6) or (2.7)-(2.8) and the residuals given by (2.16)-(2.17) or (2.18)-(2.19), if we set

$$ (3.1) \qquad \lambda\alpha \le \frac{1}{4}, \quad \lambda = \frac{\Delta x}{\Delta t}, \quad \alpha = \max_{\forall K_{j+\frac{1}{2}}} (\alpha_j, \alpha_{j+\frac{1}{2}}), $$

the schemes preserve the following property: if $\bar{h}_{j+\frac{1}{2}} \ge 0$ and $h_j \ge 0$, then $\bar{h}_{j+\frac{1}{2}}^{t+\Delta t} \ge 0$, $h_j^{t+\Delta t} \ge 0$, for any j.

*Proof.* If the low order fluxes ((2.5)-(2.6) or (2.7)-(2.8)) and residuals ((2.16)-(2.17) or (2.18)-(2.19)) are taken into account. We have

$$ \begin{aligned} \bar{h}_{j+\frac{1}{2}}^{t+\Delta t} &= \bar{h}_{j+\frac{1}{2}} - \lambda(\tilde{\mathbf{G}}_{j+1}^{(1)} - \tilde{\mathbf{G}}_{j}^{(1)}) \\ &= \bar{h}_{j+\frac{1}{2}} - \frac{\lambda}{2}(\alpha_{j+1} - \bar{u}_{j+\frac{1}{2}})h_{j+1}^{-} - \frac{\lambda}{2}(\alpha_{j} - \bar{u}_{j+\frac{1}{2}})h_{j}^{+} \\ &\quad + \frac{\lambda\alpha_{j}}{2}h_{j}^{-}(1 + \frac{\bar{u}_{j+\frac{1}{2}}}{\alpha_{j}}) + \frac{\lambda\alpha_{j+1}}{2}h_{j+1}^{+}(1 - \frac{\bar{u}_{j+\frac{3}{2}}}{\alpha_{j+1}}). \end{aligned} $$

From (2.23), the non-negativity $h_j^- \ge 0$ and $h_{j+1}^+ \ge 0$, and the fact of $1 + \frac{\bar{u}_{j+\frac{1}{2}}}{\alpha_j} \ge 0$, $1 - \frac{\bar{u}_{j+\frac{3}{2}}}{\alpha_{j+1}} \ge 0$, we have $\bar{h}_{j+\frac{1}{2}}^{t+\Delta t} \ge 0$.

Similarly, for the point values, we have

$$ \begin{aligned} h_j^{t+\Delta t} &= h_j - 2\lambda(\tilde{\mathbf{\Phi}}_{j-\frac{1}{2}}^{(1)} + \tilde{\mathbf{\Phi}}_{j+\frac{1}{2}}^{(1)}) \\ &\stackrel{(2.16),(2.17)}{=} (1 - \lambda\alpha_{j+\frac{1}{4}}\left(1 + \frac{u_j}{\alpha_{j+\frac{1}{4}}}\right)\frac{h_{j+\frac{1}{4}}^-}{h_j} - \lambda\alpha_{j-\frac{1}{4}}\left(1 - \frac{u_j}{\alpha_{j-\frac{1}{4}}}\right)\frac{h_{j-\frac{1}{4}}^+}{h_j})h_j \\ &\quad + \lambda\alpha_{j-\frac{1}{4}}h_{j-\frac{1}{4}}^-\left(1 + \frac{\bar{u}_{j-\frac{1}{2}}}{\alpha_{j-\frac{1}{4}}}\right) + \lambda\alpha_{j+\frac{1}{4}}h_{j+\frac{1}{4}}^+\left(1 - \frac{\bar{u}_{j+\frac{1}{2}}}{\alpha_{j+\frac{1}{4}}}\right) \ge 0, \end{aligned} $$

since $h_{j+\frac{1}{4}}^{\pm} \ge 0$, $-\alpha_{j-\frac{1}{4}} \le \pm u_j, \pm \bar{u}_{j-\frac{1}{2}} \le \alpha_{j-\frac{1}{4}}$, and $-\alpha_{j+\frac{1}{4}} \le \pm u_j, \pm \bar{u}_{j+\frac{1}{2}} \le \alpha_{j+\frac{1}{4}}$. □

**PROPOSITION 3.2** (Preserving still water of the first-order schemes). If the discrete solutions satisfy the still water steady state with $u = v = 0$, $w = h + B = w_0$ (constant) and with arbitrary topography, the first-order schemes can exactly preserve this equilibrium for any time.

*Proof.* At the steady state, we have $h_j + B_j = w_0$, $\bar{h}_{j+\frac{1}{2}} + \bar{B}_{j+\frac{1}{2}} = w_0$, $(hu)_j = 0$, $(hu)_{j+\frac{1}{2}} = 0$, $(hv)_j = 0$, $(\overline{hv})_{j+\frac{1}{2}} = 0$, for any j. It is easy to obtain $\frac{d}{dt}\overline{\mathbf{u}}_{j+\frac{1}{2}}^{(3)} = 0$ and $\frac{d}{dt}\mathbf{u}_j^{(3)} = 0$. We now focus on the first- and second-components.

For the average values, we apply (2.5)-(2.6) and have

$$ \frac{d}{dt} \overline{\mathbf{u}}_{j+\frac{1}{2}}^{(1)} = \frac{\alpha_{j+1}}{2\Delta x} (h_{j+1}^{+} - h_{j+1}^{-}) - \frac{\alpha_j}{2\Delta x} (h_j^{+} - h_j^{-}). $$

Since $h_j^+ = \bar{w}_{j+\frac{1}{2}} - \max(\bar{B}_{j+\frac{1}{2}}, \bar{B}_{j-\frac{1}{2}})$, $h_j^- = \bar{w}_{j-\frac{1}{2}} - \max(\bar{B}_{j+\frac{1}{2}}, \bar{B}_{j-\frac{1}{2}})$, and $\bar{w}_{j+\frac{1}{2}} = \bar{w}_{j-\frac{1}{2}}$ for any j, we obtain $h_j^+ = h_j^-$ for any j and thus $\frac{d}{dt}\overline{\mathbf{u}}_{j+\frac{1}{2}}^{(1)} = 0$. For the second
component, we have

$$
\begin{aligned}
\frac{d}{dt} \bar{\mathbf{u}}_{j+\frac{1}{2}}^{(2)} &= -\frac{g}{2\Delta x} \left( \frac{1}{2}(h_{j+1}^{+})^2 + \frac{1}{2}(h_{j+1}^{-})^2 + (h_{j+1}^{-} + \bar{h}_{j+\frac{1}{2}})(B_{j+1}^{-} - \bar{B}_{j+\frac{1}{2}}) \right. \\
& \qquad \left. - \frac{1}{2}(h_j^+)^2 - \frac{1}{2}(h_j^-)^2 + (h_j^+ + \bar{h}_{j+\frac{1}{2}})(\bar{B}_{j+\frac{1}{2}} - B_j^+) \right) \\
&= -\frac{g}{4\Delta x} \left( (h_{j+1}^+)^2 - (h_{j+1}^-)^2 + (h_j^+)^2 - (h_j^-)^2 \right) = 0,
\end{aligned}
$$

as $h_j^+ = h_j^-$ at the steady state, for all $j$. For the point values, we apply (2.16)–(2.17) and also note that $h_{j \pm \frac{1}{4}}^+ = h_{j \pm \frac{1}{4}}^-$, $h_{j-\frac{1}{4}}^+ = h_{j+\frac{1}{4}}^-$, and obtain

$$
\begin{aligned}
\frac{d}{dt} \mathbf{u}_j^{(1)} &= -\frac{1}{\Delta x} \left( \alpha_{j-\frac{1}{4}}(h_{j-\frac{1}{4}}^{+} - h_{j-\frac{1}{4}}^{-}) + g(h_j + h_{j-\frac{1}{4}}^{+})(B_j - B_{j-\frac{1}{4}}^{+}) \right. \\
& \qquad \left. - \alpha_{j+\frac{1}{4}}(h_{j+\frac{1}{4}}^{+} - h_{j+\frac{1}{4}}^{-}) + g(h_j + h_{j+\frac{1}{4}}^{-})(B_{j+\frac{1}{4}}^{-} - B_j) \right) \\
&= -\frac{1}{\Delta x} \left( g(h_j + h_{j-\frac{1}{4}}^{+})(h_{j-\frac{1}{4}}^{+} - h_j) + g(h_j + h_{j+\frac{1}{4}}^{+})(h_j - h_{j+\frac{1}{4}}^{-}) \right) \\
&= -\frac{g}{\Delta x} \left( (h_{j-\frac{1}{4}}^{+})^2 - (h_{j+\frac{1}{4}}^{-})^2 \right) = 0.
\end{aligned}
$$

For the second component, we have

$$
\begin{aligned}
\frac{d}{dt} \mathbf{u}_j^{(2)} &= -\frac{g}{\Delta x/2} \left( \frac{1}{4}(h_{j-\frac{1}{4}}^{+})^2 + \frac{1}{4}(h_{j-\frac{1}{4}}^{-})^2 + \frac{1}{2}(h_j + h_{j-\frac{1}{4}}^{+})(B_j - B_{j-\frac{1}{4}}^{+}) \right. \\
& \qquad \left. \frac{1}{4}(h_{j+\frac{1}{4}}^{+})^2 + \frac{1}{4}(h_{j+\frac{1}{4}}^{-})^2 + \frac{1}{2}(h_j + h_{j+\frac{1}{4}}^{-})(B_{j+\frac{1}{4}}^{-} - B_j) \right) \\
&= -\frac{g}{\Delta x} \left( -(h_{j-\frac{1}{4}}^{+})^2 + (h_{j+\frac{1}{4}}^{-})^2 + (h_j + h_{j-\frac{1}{4}}^{+})(h_{j-\frac{1}{4}}^{+} - h_j) \right. \\
& \qquad \left. + (h_j + h_{j+\frac{1}{4}}^{-})(h_j - h_{j+\frac{1}{4}}^{-}) \right) = 0.
\end{aligned}
$$

REMARK 3.3. If the discrete solutions satisfy the steady state characterized by $\mathbf{G} = \text{Const}$, namely,

$$ \mathbf{G}_j = \mathbf{G}_{j+1} = \mathbf{G}_{j+\frac{1}{2}}, \quad \forall K_{j+\frac{1}{2}}, $$

we then have the blending coefficients $\theta_j = 1$ and $\theta_{j+\frac{1}{2}} = 1$. This implies that, at a steady state, the scheme reduces to the high-order WB PAMPA method.

Proof. From (2.26), (2.27), and (2.28), it follows directly that $\theta_{j+\frac{1}{2}}^{\text{OE}} = 1$ as the solution is at a discrete steady state. Furthermore, using (2.5)–(2.6) and (2.16)–(2.17), we obtain $\Delta \mathbf{G}_j^{(1)} = 0$, $\Delta \bar{\mathbf{\Phi}}_{j-\frac{1}{2}}^{(1)} = 0$, and $\Delta \bar{\mathbf{\Phi}}_{j+\frac{1}{2}}^{(1)} = 0$. Then, from (2.24) and (2.25), we obtain

$$ \theta_j = 1, \quad \theta_{j-\frac{1}{2}} = 1, \quad \theta_{j+\frac{1}{2}} = 1. $$

PROPOSITION 3.4 (Super-consistency for arbitrary moving equilibria). The WB PAMPA scheme with Gauss-Lobatto flux globalization preserves exactly discrete moving steady states obtained by integrating the nonlinear ODE

$$ (3.2) \qquad \mathbf{f}'(\mathbf{u})(x) = \mathbf{S}(\mathbf{u}(\mathbf{f}), x) $$
with the LobattoIIIA collocation method when using (2.10), or with the sc-LobattoIII method when using (2.13)-(2.14). Provided $\mathbf{f}(\mathbf{u})$ is a one to one mapping and the inverse $\mathbf{u}(\mathbf{f})$ is bounded and uniquely defined, these discrete solutions are formally fourth order accurate.

*Proof.* We write the proof for the case of the LobattoIIIA method. The one for the sc-LobattoIII is almost identical. Steady data obtained integrating (3.2) with the RK-LobattoIIIA method, verify within each element the nonlinear algebraic system

$$ \mathbf{f}_{j+\frac{1}{2}} - \mathbf{f}_j = \sum_{i=1}^{3} \omega_i' \mathbf{S}(\mathbf{u}(\mathbf{f}_{j+\frac{i-1}{2}}), x_{j+\frac{i-1}{2}}) $$

$$ \mathbf{f}_{j+1} - \mathbf{f}_j = \sum_{i=1}^{3} \omega_i \mathbf{S}(\mathbf{u}(\mathbf{f}_{j+\frac{i-1}{2}}), x_{j+\frac{i-1}{2}}) $$

with $\{\omega'_i, \omega_i\}_{i=1,3}$ the LobattoIIIA Butcher coefficients, also used in (2.10) to integrate the source. Combining (2.9) and Remark 3.3 this leads to

$$ \hat{\mathbf{G}}_{j+1} = \mathbf{G}_{j+1} = \mathbf{G}_j = \hat{\mathbf{G}}_j, \quad \forall K_{j+\frac{1}{2}} \Rightarrow \frac{d}{dt} \bar{\mathbf{u}}_{j+\frac{1}{2}} = 0. $$

For the nodal updates, using (2.20)-(2.21) and Remark 3.3, we find easily that

$$ \hat{\bar{\mathbf{F}}}_{j-\frac{1}{2}} = \hat{\bar{\mathbf{F}}}_{j+\frac{1}{2}} = 0 \Rightarrow \frac{d}{dt} \mathbf{u}_j = 0, \quad \forall j. $$

Thus the LobattoIIIA discrete solutions are in the kernel of the flux globalization based WB PAMPA scheme. The invertibility of $\mathbf{f}(\mathbf{u})$ is enough to claim that the nodal accuracy of the LobattoIIIA is preserved, which is order 4 (cf. [32]). The proof for the case of the sc-LobattoIII method follows the same ideas, with the obvious changes in quadrature coefficients and points. □

REMARK 3.5 (Shallow water: invertibility of $\mathbf{f}(\mathbf{u})$ and critical point detection). Let us consider the shallow water flux and set $\mathbf{f}(\mathbf{u}) = (f_1, f_2, f_3)^\top$. For moving-water solutions, so $u \neq 0$, and strictly positive water depth, the invertibility of the flux require detecting critical points where $u = \sqrt{gh}$. It is quite simple to verify that $f_2 > 0$, and

$$ \kappa := g \frac{f_1^4}{f_2^3} = \frac{8\mathrm{Fr}^4}{(2\mathrm{Fr}^2 + 1)^3}, \quad \mathrm{Fr} := \frac{u}{\sqrt{gh}}. $$

This gives an immediate criterion to verify if a flux value corresponds to a critical point where $\kappa = 8/27$.

As discussed in [42, 34] the static case is a particular one which requires a reformulation of the consistency analysis. As already said, for the shallow water equations, the sc-LobattoIII approach has the advantage of always integrating exactly the hydrostatic term. For this reason we now focus on this approach, which is also the one used in the numerical results.

PROPOSITION 3.6 (Preservation of still water states). The WB PAMPA scheme with Gauss-Lobatto global flux quadrature using the sc-LobattoIII approach

1. preserves exactly still water states with $u = v = 0$, $w = h + B = w_0$ constant, and with arbitrary topography;
2. preserves exactly the fourth order accurate discrete states at rest computed with the sc-LobattoIII collocation method applied to

$$f_2'(h) = -ghB_x + fhv$$

with given $v(x) \neq 0$, and for arbitrary topography.

*Proof.* To prove the first part we note that the relation

$$\frac{d}{dt} \bar{u}_{j+\frac{1}{2}} = \frac{d}{dt} u_j = 0$$

is trivial for the first component and last component (appear in (1.6)). We thus focus on the expression of $G$ for the second component for which we can write

$$G_j^{(2)} - G_{j-\frac{1}{2}}^{(2)} = \int_{x_{j-\frac{1}{2}}}^{x_j} \left( g\left(\frac{h^2}{2}\right)_x + ghB_x \right) dx = g \frac{h_j^2 - h_{j-\frac{1}{2}}^2}{2} + \int_{x_{j-\frac{1}{2}}}^{x_j} ghB_x dx$$

and similarly for the $G_{j+\frac{1}{2}}$. For the still-water steady-state solutions satisfying $g(\frac{h^2}{2})_x + ghB_x = 0$. So, once the last integral is evaluated exactly by the Gauss-Lobatto formula, we immediately have $G_j^{(2)} - G_{j-\frac{1}{2}}^{(2)} = 0$. This is achieved by using the sub-cell LobattoIII method in (2.13)-(2.14), which is exactly for the polynomial $hB_x$ of degree 3.

Note that when the LobattoIIIA method is applied, in order to exactly preserve the still-water states, one can use the same modifications of the hydrostatic source proposed in [52, 42, 34] to replace the part of the source terms as

$$-ghB_x = -gwB_x + g\left(\frac{B^2}{2}\right)_x.$$

So we can write with no approximation

$$\int_{x_{j-\frac{1}{2}}}^{x_j} gwB_x dx - \int_{x_{j-\frac{1}{2}}}^{x_j} g\left(\frac{B^2}{2}\right)_x dx = \int_{x_{j-\frac{1}{2}}}^{x_j} gwB_x dx - g\frac{B_j^2 - B_{j-\frac{1}{2}}^2}{2}$$

One easily checks that for $w = h + B = w_0$ constant

$$\begin{aligned} \mathbf{G}_j^{(2)} - \mathbf{G}_{j-\frac{1}{2}}^{(2)} &= g \frac{h_j^2 - h_{j-\frac{1}{2}}^2}{2} + \int_{x_{j-\frac{1}{2}}}^{x_j} gw_0 B_x dx - g \frac{B_j^2 - B_{j-\frac{1}{2}}^2}{2} \\ &= g \frac{h_j^2 - h_{j-\frac{1}{2}}^2}{2} + gw_0 (B_j - B_{j-\frac{1}{2}}) - g \frac{B_j^2 - B_{j-\frac{1}{2}}^2}{2} = 0. \end{aligned}$$

We invoke Remark 3.3 to deduce that $\hat{G}_{j+1}^{(2)} = \hat{G}_j^{(2)} = \hat{G}_{j+\frac{1}{2}}^{(2)} = \text{constant}$, and obtain the result sought.

The second statement of the proposition is a special case of Proposition 3.4, considering only the second component of the flux $f_2 = gh^2/2$ which is always invertible. The proof is identical to that of Proposition 3.4. □

4. **Numerical examples.** We have thoroughly tested the flux globalization based PP WB PAMPA scheme proposed on the Saint-Venant system (1.5) to assess
convergence property, verify the PP and WB properties presented in the previous section, as well as to investigate its application to resolution of Riemann problems and to test its robustness. The benchmarks used are quite classical in literature. They allow to verify the correct implementation and convergence of the proposed method.

The time-dependent ODE systems (2.2a) and (2.2b) are integrated using the three-stage third-order strong stability preserving (SSP) Runge-Kutta method (see, e.g.,[27, 28]). All numerical simulations have been run at CFL=0.2. Unless specified differently, we use the zero-order extrapolation boundary conditions and take the acceleration due to gravity $g = 9.812$. In the following examples, the numerical results obtained using the developed high-order schemes are denoted as “HO-PAMPA” and the corresponding low-order schemes alone are denoted as “LO-PAMPA”.

In Examples 1–7, the PP WB PAMPA scheme is applied to solve the Saint–Venant system with/without Manning friction term (1.5), while in Examples 8–11, the proposed scheme is applied to solve the rotating shallow water equations (1.6).

**Example 1—Accuracy Test.** In the first example, a slight modification from [4] with the Manning coefficient $n = 0.05$, we study a problem with smooth solutions in a finite time frame to check the experimental rate of convergence. The initial conditions and non-flat bottom topography are:

$$h(x, 0) = 0.3 \left[ 1 + e^{-\frac{(x-0.5)^2}{0.05^2}} \right] - 0.2 \cos(6\pi x), \quad u(x, 0) = 0, \quad B(x) = 0.2(1 + \cos(6\pi x)),$$

prescribed in a computational domain $[0, 1]$ and subject to the periodic boundary conditions.

We compute the numerical solution until the final time $t = 0.03$ using the high-order PP WB PAMPA scheme on a sequence of uniform meshes with $\Delta x = 1/64, 1/128, 1/256, 1/512, 1/1024, 1/2048$, and $1/4096$. After that, we measure the discrete $L^1$-errors for the cell boundary point and average values of variables $h$ and $hu$ using the Runge formulae provided in [4]. Equipped with these errors, we compute the corresponding convergence rates and report the results in Table 2 for the point values and cell averages, respectively. From these obtained results, we can clearly see that the expected experimental third order of accuracy is achieved for the proposed schemes.

**Example 2—Still Water Equilibrium and Its Small Perturbation.** In the second example, taken from [4], we examine the still water preservation and positivity-preserving properties of the proposed high-order and first-order schemes. We consider a non-flat bottom topography,

$$B(x) = \begin{cases} 2(\cos(10\pi(x + 0.3)) + 1), & x \in [-0.4, -0.2], \\ 0.5(\cos(10\pi(x - 0.3)) + 1), & x \in [0.2, 0.4], \\ 0, & \text{otherwise,} \end{cases}$$

and the following initial data

$$h(x, 0) = 4.000001 - B(x), \quad u(x, 0) = 0,$$

in the computational domain $[-1, 1]$. The Manning friction coefficient is $n = 0$. We run the simulations until a long final time $t = 10$ using the schemes proposed with 50 uniform cells and report the discrete $L^1$, $L^2$, and $L^\infty$-errors of the average values of $h$ and $hu$ in Table 3. All the calculated errors are within the level of round-off error, clearly validating Proposition 3.2 and Proposition 3.6.
<caption>TABLE 2</caption>

Example 1: $L^1$-errors and experimental convergence rates of the cell boundary point values and average values of variables $h$ and $hu$.

<table><caption>TABLE 3</caption><thead><tr><th colspan="5">errors and rates for u<sub>j</sub></th></tr><tr><th>Δx</th><th>L<sup>1</sup>-error in h</th><th>Rate</th><th>L<sup>1</sup>-error in hu</th><th>Rate</th></tr></thead><tbody><tr><td>1/256</td><td>4.16 × 10<sup>-4</sup></td><td>-</td><td>7.37 × 10<sup>-4</sup></td><td>-</td></tr><tr><td>1/512</td><td>9.69 × 10<sup>-5</sup></td><td>2.10</td><td>1.80 × 10<sup>-4</sup></td><td>2.04</td></tr><tr><td>1/1024</td><td>1.34 × 10<sup>-5</sup></td><td>2.85</td><td>2.59 × 10<sup>-5</sup></td><td>2.79</td></tr><tr><td>1/2048</td><td>1.66 × 10<sup>-6</sup></td><td>3.01</td><td>3.24 × 10<sup>-6</sup></td><td>3.00</td></tr><tr><td>1/4096</td><td>2.02 × 10<sup>-7</sup></td><td>3.04</td><td>3.94 × 10<sup>-7</sup></td><td>3.04</td></tr></tbody></table><table><thead><tr><th colspan="5">errors and rates for u&#773;<sub>j+1/2</sub></th></tr><tr><th>Δx</th><th>L<sup>1</sup>-error in h</th><th>Rate</th><th>L<sup>1</sup>-error in hu</th><th>Rate</th></tr></thead><tbody><tr><td>1/256</td><td>4.72 × 10<sup>-4</sup></td><td>-</td><td>1.01 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>1/512</td><td>8.21 × 10<sup>-5</sup></td><td>2.52</td><td>1.78 × 10<sup>-4</sup></td><td>2.51</td></tr><tr><td>1/1024</td><td>1.14 × 10<sup>-5</sup></td><td>2.85</td><td>2.46 × 10<sup>-5</sup></td><td>2.86</td></tr><tr><td>1/2048</td><td>1.40 × 10<sup>-6</sup></td><td>3.02</td><td>3.01 × 10<sup>-6</sup></td><td>3.03</td></tr><tr><td>1/4096</td><td>1.70 × 10<sup>-7</sup></td><td>3.04</td><td>3.66 × 10<sup>-7</sup></td><td>3.04</td></tr></tbody></table>

Example 2: Errors in $h$ and $hu$ computed by high-order and first-order schemes.

<table><thead><tr><th>Schemes</th><th>Variables</th><th>L<sup>1</sup>-error</th><th>L<sup>2</sup>-error</th><th>L<sup>∞</sup>-error</th></tr></thead><tbody><tr><td rowspan="2">HO-PAMPA</td><td><i>h</i></td><td>2.77 × 10<sup>-14</sup></td><td>2.44 × 10<sup>-14</sup></td><td>3.02 × 10<sup>-14</sup></td></tr><tr><td><i>hu</i></td><td>2.20 × 10<sup>-13</sup></td><td>5.96 × 10<sup>-13</sup></td><td>2.93 × 10<sup>-12</sup></td></tr><tr><td rowspan="2">LO-PAMPA</td><td><i>h</i></td><td>1.29 × 10<sup>-15</sup></td><td>3.51 × 10<sup>-15</sup></td><td>1.02 × 10<sup>-14</sup></td></tr><tr><td><i>hu</i></td><td>1.98 × 10<sup>-14</sup></td><td>1.77 × 10<sup>-14</sup></td><td>2.71 × 10<sup>-14</sup></td></tr></tbody></table>

To further verify the WB property for preserving still-water equilibrium, we add a small Gaussian-shaped perturbation to the stationary water depth and test the ability of the proposed schemes to accurately capture this small perturbation. The small perturbation of the form $10^{-6} \times e^{-200x^2}$ is superimposed onto the initial cell average of water depth. We compute the numerical solution at three different times: $t = 0.02$, 0.04, and 0.06, using 100 uniform cells. A reference solution computed by the first-order scheme on a uniform mesh with 2000 cells is also provided. The differences between the computed and background cell averages of water depth are plotted in Figure 2.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Example 2: Time snapshots of the small perturbation in $h$.</figcaption></figure>
As one can observe, both high-order and low-order PP WB PAMPA schemes can capture the dynamics of the small perturbation, there are no spurious oscillations: the initial perturbation splits into two humps, which are then propagating into opposite directions, respectively. The results generated on coarse and fine meshes are consistent and converge to the same solution profile.

**Example 3—Super-convergent to Moving Water Equilibria without friction.** In the third example, we study the convergence of the numerical solution computed by the proposed PP WB PAMPA scheme towards steady flow over a continuous bottom topography given by

$$ (4.1) \qquad B(x) = \begin{cases} 0.2 - 0.05(x - 10)^2, & \text{if } 8 \le x \le 10, \\ 0, & \text{else.} \end{cases} $$

The goal here is to verify the theoretical predictions on the super convergence on the smooth moving-water equilibria for the scheme. We compute the convergent flow solutions for two classical smooth states: supercritical and subcritical, which depend on the following initial and boundary conditions:

$$ \begin{array}{ll} \text{Case (I):} & \begin{cases} h(x, 0) = 2 - B(x), & hu(x, 0) = 0, \\ h(0, t) = 2, & hu(0, t) = 24; \end{cases} \\ \text{Case (II):} & \begin{cases} h(x, 0) = 2 - B(x), & hu(x, 0) = 0, \\ hu(0, t) = 4.42, & h(25, t) = 2. \end{cases} \end{array} $$

We let the scheme evolve until a given finite time $t = 500$. The numerical solutions computed on a uniform mesh with 100 cells are presented in Figures 3 and 4 for supercritical and subcritical flows, respectively. As one can see, the water depth obtained by the proposed PP WB PAMPA scheme is very close to the corresponding steady states (reported in [4, Example 3]). The global flux $G$ with the first component being the momentum $hu$ have reached to the expected constant states. This indicates that the numerical solutions have converged to the discrete steady state satisfying (1.4).

Next, we qualify the convergence between the discrete water height obtained by this text and the analytical steady-state water height. To this end, we consider the mesh covered by a set of uniform cells with numbers: 50, 100, 200, and 400. We compute the discrete $L^1$- and $L^\infty$-errors between the discrete solutions and the analytical steady states given in [4, Example 3]. The logarithmic plots for errors in water depth versus number of cells are as shown in Figure 5. We can see that the convergence rates are as predicted in Proposition 3.4, i.e., the convergence is of order four for the third-order PP WB PAMPA scheme.

**Example 4—Discrete Moving Water Equilibria and Their Small Perturbations.** In the fourth example, we consider a frictional case with Manning friction parameter $n = 0.05$ and consider the same initial and boundary conditions as in Example 3. We plot the numerical solutions obtained by the proposed PP WB PAMPA scheme at $t = 500$ using 100 uniform grid cells in Figures 6 and 7. One may notice that as in the frictionless case, both momentum $hu$ and global flux $G^{(2)}$ computed by the proposed scheme are very accurate and have reached at the constant steady state. Moreover, with the influence of friction, the water depth and the water surface differ with those obtained in the frictionless cases.

Next, we test the performance of the proposed PP WB PAMPA scheme in the presence of small perturbations of the moving water equilibria. We consider the
<figure><img src="image_4.png" /><figcaption>FIG. 3. Example 3: Convergent solutions for supercritical flow. The deviation parameter $\epsilon = 10^{-12}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Example 3: Convergent solutions for subcritical flow. The deviation parameter $\epsilon = 10^{-12}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. Example 3: logarithmic plots for errors in water depth versus number of cells.</figcaption></figure>

following two sets of the steady-state data:

Case (I): $hu(x, 0) = 24$, $G^{(2)} = 307.624$;

Case (II): $hu(x, 0) = 4.42$, $G^{(2)} = 31.7008$.

We note that the data are given in terms of the global quantities rather than in water depth. In order to start the computation at time $t = 0$, one needs to obtain the values of $\bar{h}_{j+\frac{1}{2}}$ and $h_j$. We, therefore, use Newton's method to numerically solve the
<figure><img src="image_7.png" /><figcaption>FIG. 6. Example 4–Case (I): Convergent solutions for supercritical flow with friction. The deviation parameter $\epsilon = 10^{-12}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Example 4–Case (II): Convergent solutions for subcritical flow with friction. The deviation parameter $\epsilon = 10^{-12}$ and the constant value $G_s^{(2)} = 31.700836562966$.</figcaption></figure>

following nonlinear equations:

$$ (4.2) \qquad \begin{aligned} \mathbf{G}_j^{(2)} &= \frac{((hu)_j)^2}{h_j} + \frac{g}{2}(h_j)^2 - \mathbf{R}_j^{(2)}, \\ \mathbf{G}_{j+\frac{1}{2}}^{(2)} &= \frac{((hu)_{j+\frac{1}{2}})^2}{(h_{j+\frac{1}{2}})^2} + \frac{g}{2}(h_{j+\frac{1}{2}})^2 - \mathbf{R}_{j+\frac{1}{2}}^{(2)}, \end{aligned} \quad \forall j, $$

where $\mathbf{R}_j^{(2)}$ and $\mathbf{R}_{j+\frac{1}{2}}^{(2)}$ are given in (2.13). After obtaining the point values of water depth, it is easy to compute the initial steady-state average value $\bar{h}_{j+\frac{1}{2}}^{eq}$ using Simpson's rule. Using these prepared data, we repeat the simulations from the still water case (Example 2) to verify the exact fully well-balanced property of the proposed HO-PAMPA scheme. We run the simulations using 100 cells until the final time $t = 1000$, and compute the errors between the numerical and discrete prepared steady-state solutions. These errors are within machine precision and are omitted here for saving space.

Finally, we add a small perturbation of the form $10^{-3}e^{-80(x-6)^2}$ to the steady-state water depth. We compute the numerical solutions at $t = 1$ and $t = 1.5$ for the Cases (I) and (II), respectively. The differences between the steady-state and numerical water depths are plotted in Figure 8. The results demonstrate that the proposed PP WB PAMPA scheme accurately captures the propagation of the small perturbation, in agreement with those reported in [4, Example 8].
<figure><img src="image_9.png" /><figcaption>FIG. 8. Example 4: The difference between the computed and the background discrete steady-state cell averages of water depth computed using 100 uniform cells.</figcaption></figure>

**Example 5—Riemann Problems over a Flat Bottom.** In the fifth example, we consider two test cases, taken from [55], to demonstrate the positivity-preserving property of the proposed PP WB PAMPA scheme. The bottom topography is flat ($B(x) = 0$) and the initial data are given by

$$ \text{Test 1: } h(x, 0) = \begin{cases} 10, & x \in [-300, 0], \\ 0, & x \in (0, 300], \end{cases} \qquad hu(x, 0) = 0, $$

$$ \text{Test 2: } h(x, 0) = \begin{cases} 5, & x \in [-200, 0], \\ 10, & x \in (0, 400], \end{cases} \qquad hu(x, 0) = \begin{cases} 0, & x \in [-200, 0], \\ 400, & x \in (0, 400], \end{cases} $$

We compute the numerical solution using the proposed PP WB PAMPA scheme on a uniform mesh with 250 cells. The solutions at three different time snapshots are plotted in Figure 9. No negative water height is generated during the computation, and a good agreement with those obtained in [55].

**Example 6—Riemann Problems over a Non-flat Bottom.** In the sixth example, we consider a dam-break over a non-flat bottom. The initial conditions are given by:

$$ h(x, 0) = \begin{cases} 5 - B(x), & x < 0, \\ 1, & x \ge 0, \end{cases} \qquad u(x, 0) = 0, $$

where the topography function is defined as

$$ B(x) = \begin{cases} 2(\cos(10\pi(x + 0.3)) + 1), & -\frac{2}{5} \le x \le -\frac{1}{5}, \\ \frac{1}{2}(\cos(10\pi(x - 0.3)) + 1), & \frac{1}{5} \le x \le \frac{2}{5}, \\ 0, & \text{otherwise.} \end{cases} $$

The computational domain is $[-1, 1]$, the gravitational constant is $g = 1$, and the numerical solutions are computed on a uniform mesh with 300 cells until the final time $t = 0.3$. We plot the obtained solution in Figure 10, together with a reference solution computed by the LO-PAMPA scheme using 3000 cells. We can see that, the results obtained by the HO-PAMPA and LO-PAMPA are consistent with the
<figure><img src="image_10.png" /><figcaption>FIG. 9. Example 5: Numerical solutions over a flat bottom.</figcaption></figure>

reference solution and the HO-PAMPA scheme outperforms the LO-PAMPA scheme. Moreover, the PP property of the proposed schemes is also verified, as one can clearly observe that the area near $x = -0.3$ is now almost dry and well resolved.

<figure><img src="image_11.png" /><figcaption>FIG. 10. Example 6: Numerical solutions over a non-flat bottom.</figcaption></figure>
**Example 7—Parabolic Bowl.** In the seventh example, taken from [54, 55], we consider the flow with a parabolic bottom topography,

$$B(x) = h_0(x/a)^2, \quad h_0 = 10, \quad a = 3000,$$

prescribed in the domain $[-5000, 5000]$. The analytical water surface for (1.5) without the friction term ($n = 0$) is given by

$$(4.3) \qquad h(x,t) + B(x) = h_0 - \frac{b^2}{4g} \cos(2\omega t) - \frac{b^2}{4g} - \frac{bx}{2a} \sqrt{\frac{8h_0}{g}} \cos(\omega t),$$

where $\omega = \sqrt{2gh_0/a}$ and $b = 5$. The exact location of the wet/dry front takes the form:

$$(4.4) \qquad x_0 = -\frac{b\omega a^2}{2gh_0} \cos(\omega t) \pm a.$$

The initial condition is then given by (4.3)–(4.4) and $hu(x, 0) = 0$. We run the simulations until the final time $t = 6000$ using 250 uniform cells. The averaged water surface level ($\bar{h}_{j+\frac{1}{2}} + \bar{B}_{j+\frac{1}{2}}$, $\forall j$) at different times are shown in Figure 11, together with the exact solution given by (4.3)–(4.4) for a comparison. As one can clearly see, the numerical solutions present a nice agreement with the exact ones and also further confirms the PP property of the proposed PP WB PAMPA schemes.

<figure><img src="image_12.png" /><figcaption>FIG. 11. Example 7: The water surface level in the parabolic bowl problem at different times.</figcaption></figure>

**Example 8—Convergence to Moving-water equilibrium with Coriolis forces.** In the eighth example, taken from [13], we consider the following initial data
and boundary conditions:

$$ \begin{cases} h(x, 0) = 0.33, & hu(x, 0) = 0, & hv(x, 0) = 0, \\ h(25, t) = 0.33, & hu(0, t) = 0.18, & hv(0, t) = 0. \end{cases} $$

The bottom topography is the same one given by (4.1) and the computational domain is also [0, 25]. We set the Coriolis force $f = \frac{2\pi}{50} + 0.01x$ and run the simulations on a uniform mesh with 100 cells until a final time $t = 1000$. The obtained numerical solutions and the discrete values of global flux ($G^{(2)}$ and $G^{(3)}$, where $G^{(1)} = hu$) are reported in Figure 12. From these results, we clearly see that, the discrete steady states are reached. We can also perform the same test as in Example 4 by adding a small perturbation to the steady-state solution and to examine the propagation captured by the proposed PP WB PAMPA scheme. Whereas, for saving space, we omit the results here.

<figure><img src="image_13.png" /><figcaption>FIG. 12. Example 8: Convergent solutions and discrete values of global fluxes for the flow with Coriolis forces. $G^{(2)} = 3.8948194772284$ and $\epsilon = 10^{-2}$.</figcaption></figure>

**Example 9—Geostrophic Equilibrium over a Flat Bottom.** In the ninth example, we consider the following initial conditions, taken from [20]:

$$ u(x, 0) = 0, \quad v(x, 0) = \frac{2g}{f}xe^{-x^2}, $$

where $g = 1$ and $f = 10$. The discrete values for $h(x)^{eq}$ are calculated analytically by solving the equations in (4.2) with $G^{(2)} = 2$ in the domain $[-10, 10]$ covered by 50 cells. In Figure 13, we plot the prepared discrete steady-state water height. We then use these initial setting to run the simulation until the final time $t = 100$ and report
the errors between the numerical solutions and the initial data in Table 4. Again, the WB property of the proposed scheme is verified.

<figure><img src="image_14.png" /><figcaption>FIG. 13. Example 9: Prepared discrete steady-state water height.</figcaption></figure>

<table><caption>TABLE 4
Example 9: Errors in $h$, $hu$, and $hv$ computed by high-order scheme.</caption><thead><tr><th>Variables</th><th>L<sup>1</sup>-error</th><th>L<sup>2</sup>-error</th><th>L<sup>∞</sup>-error</th></tr></thead><tbody><tr><td>h</td><td>0</td><td>0</td><td>0</td></tr><tr><td>hu</td><td>6.63 × 10<sup>-16</sup></td><td>2.49 × 10<sup>-16</sup></td><td>2.29 × 10<sup>-16</sup></td></tr><tr><td>hv</td><td>5.66 × 10<sup>-15</sup></td><td>1.97 × 10<sup>-15</sup></td><td>1.09 × 10<sup>-15</sup></td></tr></tbody></table>

**Example 10—Geostrophic Equilibrium over a Non-flat Bottom.** In the tenth example, we consider a non-flat bottom topography,

$$B(x) = \begin{cases} 0.25(\cos(10\pi(x - 0.8)) + 1), & \text{if } 0.7 \le x \le 0.9, \\ 0, & \text{otherwise.} \end{cases}$$

and the initial conditions again satisfies a geostrophic equilibrium given by

$$G^{(2)} = 2, \quad u(x, 0) = 0, \quad v(x, 0) = 0.05 \sin(2\pi x),$$

over the domain $[0, 1]$. We use $g = 1$ and $f = 10$ here. The initial values of water depth is prepared by solving the nonlinear equations (4.2). We first run the simulation using 20 cells until a final time $t = 20$ and report the errors in Table 5. As one can see again, the proposed PP WB PAMPA scheme is capable of exactly preserving the geostrophic equilibrium within machine accuracy even with a non-flat bottom topography.

To further examine the WB property of the scheme on handling with the geostrophic equilibrium, we study the propagation of a small perturbation. To this end, we add a small perturbation, $10^{-3}e^{-200(x-0.55)^2}$, to the discrete steady-state water depth and compute the difference between the numerical solution and the background steady states. The results at times $t = 0, 0.1$, and $0.2$ are shown in Figure 14. We can clearly see that the proposed PP WB PAMPA scheme well captures a small perturbation.
<table><caption>TABLE 5
Example 10: Errors in $h$, $hu$, and $hv$ computed by high-order scheme.</caption><thead><tr><th>Variables</th><th>L<sup>1</sup>-error</th><th>L<sup>2</sup>-error</th><th>L<sup>∞</sup>-error</th></tr></thead><tbody><tr><td>h</td><td>6.48 × 10<sup>-15</sup></td><td>8.04 × 10<sup>-15</sup></td><td>1.02 × 10<sup>-14</sup></td></tr><tr><td>hu</td><td>1.89 × 10<sup>-16</sup></td><td>2.70 × 10<sup>-16</sup></td><td>7.39 × 10<sup>-16</sup></td></tr><tr><td>hv</td><td>4.62 × 10<sup>-15</sup></td><td>6.53 × 10<sup>-15</sup></td><td>1.58 × 10<sup>-14</sup></td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>FIG. 14. Example 10: Difference between the numerical solution and the background steady states.</figcaption></figure>

**Example 11—Rossby Adjustment in an Open Domain.** In the last example, we numerically investigate the Rossby adjustment problem with the constant Coriolis parameter $f = 1$ and $g = 1$, which was studied previously in [20, 14]. We consider the following initial conditions, which correspond to a jet over a flat $h$:

$$h(x, 0) = 1, \quad u(x, 0) = 0, \quad v(x, 0) = 2N_2(x),$$

where

$$N_2(x) = \frac{(1 + \tanh(2x + 2))(1 - \tanh(2x - 2))}{(1 + \tanh(2))^2}$$

and the bottom is flat $B(x) = 0$. The simulation is run on a uniform mesh with 20000 cells. Time evolution of the water depth along an inertial period $T_f = \frac{2\pi}{f}$ is shown in Figure 15. As one can see, fast inertia gravity waves are emitted as a consequence of the initial momentum disturbance, and two shocks arise in the wave front.

**5. Conclusion.** In this work, we have developed a new formulation of the Point-Average-Moment PolynomiAl-interpreted (PAMPA) method via the flux globalization to construct positivity-preserving and fully well-balanced numerical schemes for one-dimensional shallow water models. The discrete fully well-balanced condition is achieved by incorporating the source terms into the global flux, resulting in a quasi-conservative system that enables a clear characterization of the discrete steady states and provable superconvergence toward general analytical equilibria. In practice, the required integrals are evaluated using direct applications of Gauss-Lobatto quadrature. To ensure positivity of the water depth, capture wet-dry fronts, and suppress spurious oscillations near strong shocks, first-order schemes based on the hydrostatic reconstruction techniques are introduced. Furthermore, a monolithic convex limiting
<figure><img src="image_16.png" /><figcaption>FIG. 15. Example 11: time evolution of $\bar{h}$ at different times.</figcaption></figure>

strategy is employed to determine the blending coefficients for the high- and low-order fluxes and residuals. The proposed schemes are proven to exactly preserve still-water steady states, approximately preserve moving-water and geostrophic equilibria with superconvergent accuracy, and maintain the positivity of the water depth.

Extensive numerical tests for the shallow water equations, with/without Manning friction and Coriolis effect, validate all the theoretical properties on various equilibria, including still-water, moving-water, and geostrophic steady states. Applications to complex problems involving sharp fronts and different wave interactions further demonstrate the accuracy and robustness of the proposed method.

Compared with existing flux globalization based well-balanced methods for one-dimensional shallow water models, the proposed flux globalization based PP WB PAMPA scheme offers several advantages: i) The PP WB PAMPA scheme, using global flux quadrature, is exactly well-balanced for still-water steady states, as proven in Proposition 3.6. The sub-cell LobattoIII method has been used for the source term integral in the half cell. We have shown that for the shallow water equations, it automatically exactly preserves the still-water states. In contrast, the flux globalization based DGSEM methods in [42, 24], where the LobattoIIIA has been used to compute the integral on the half cell, require an additional source term reformulation following the approach of [52], while the flux globalization based central-upwind scheme of [17] does not achieve this property. ii) The PP WB PAMPA scheme is positivity-preserving and capable of handling wet-dry fronts through blending with a positivity- and still-water preserving first-order scheme. This capability is absent in other flux globalization methods, such as those in [42, 17, 24]. iii) A simple oscillation-eliminating parameter is introduced to adaptively adjust the blending coefficients near shock waves, enabling the PP WB PAMPA scheme to robustly capture genuinely dis-
continuous solutions. This improvement is not realized in the DGSEM framework of [42]. iv) The proposed formulation for preserving discrete moving-water steady states is straightforward, fully local, and computationally efficient, whereas the flux globalization based central-upwind schemes of [17, 13, 21] require solving additional nonlinear equations to recover well-balanced reconstructed point values (a generalized hydrostatic reconstruction step). v) We note that, despite the use of integrals to define the global fluxes (see e.g. equation (2.10)), Proposition 2.2 shows that the final numerical method is actually compact. This property enables their blending with the first-order schemes. For the evolution of point values, the computation in (2.21) similarly confirms this compactness. A current limitation of the PP WB PAMPA scheme is its inability to handle discontinuous bottom topography. This could be done with some estimation of the jump in $\mathbf{R}$; see, e.g., [13, 24, 34].

The extensions to multidimensional problems and to more complex systems will be investigated in the near future.

## Appendix A. Consistency of sub-cell LobattoIII collocation method.

We prove here the statement of proposition 2.1. We use explicitly the fact that the main unknowns are approximations of

$$ \mathbf{u}_j \approx \mathbf{u}(x_j), \mathbf{u}_{j+1/2} \approx \mathbf{u}(x_j + \Delta x/2), \mathbf{u}_{j+1} \approx \mathbf{u}(x_j + \Delta x), $$

and we formally replace the above with sampled values of a smooth exact solution in the discrete equations to evaluate the consistency error. Without loss of generality, we consider the problem

$$ \mathbf{f}' = \mathbf{S}(\mathbf{u}(\mathbf{f}), x) $$

which is the steady ODE of the balance law, with classical assumptions that $\mathbf{f}(\mathbf{u})$ is a one-to-one smooth mapping, and $\mathbf{S}$ has bounded derivatives.

Using a truncated Taylor series expansion we can readily show that

$$ \begin{aligned} & \mathbf{f}_{j+1} - \mathbf{f}_j - \frac{\Delta x}{6} \mathbf{S}(\mathbf{u}_j, x_j) - 2 \frac{\Delta x}{3} \mathbf{S}(\mathbf{u}_{j+\frac{1}{2}}, x_{j+\frac{1}{2}}) - \frac{\Delta x}{6} \mathbf{S}(\mathbf{u}_{j+1}, x_{j+1}) \\ &= \mathbf{f}_{j+1} - \mathbf{f}_j - \frac{\Delta x}{6} \mathbf{f}'_j - 2 \frac{\Delta x}{3} \mathbf{f}'_{j+\frac{1}{2}} - \frac{\Delta x}{6} \mathbf{f}'_{j+1} \\ &= \mathbf{f}_{j+1} - \mathbf{f}_j - \mathbf{f}'_j \Delta x - \mathbf{f}''_j \frac{\Delta x^2}{2} - \mathbf{f}'''_j \frac{\Delta x^3}{6} - \mathbf{f}''''_j \frac{\Delta x^4}{24} - \mathbf{f}''''''_j \frac{5 \Delta x^5}{576} = \mathcal{O}(\Delta x^5). \end{aligned} $$

This shows that

$$ \sum_j \{\mathbf{f}_{j+1} - \mathbf{f}_j - \frac{\Delta x}{6} \mathbf{S}(\mathbf{u}_j, x_j) - 2 \frac{\Delta x}{3} \mathbf{S}(\mathbf{u}_{j+\frac{1}{2}}, x_{j+\frac{1}{2}}) - \frac{\Delta x}{6} \mathbf{S}(\mathbf{u}_{j+1}, x_{j+1})\} = \mathcal{O}(\Delta x^4), $$

provided we have at least $\mathbf{u}_{j+\frac{1}{2}} - \mathbf{u}(x_{j+\frac{1}{2}}) = \mathcal{O}(\Delta x^4)$. To show this we first compute (using again a truncated Taylor series expansion)

$$ \mathbf{u}_{j+1/4} = \frac{3}{8}\mathbf{u}_j + \frac{3}{4}\mathbf{u}_{j+1/2} - \frac{1}{8}\mathbf{u}_{j+1} = \mathbf{u}(x_j + \Delta x/4) - \frac{\Delta x^3}{128}\mathbf{u}''' + \mathcal{O}(\Delta x^4). $$
This allows to show that, as long as $\mathbf{S}$ has bounded derivatives, then

$$ 
\begin{aligned}
\mathbf{f}_{j+\frac{1}{2}} - \mathbf{f}_j & - \frac{\Delta x}{12} \mathbf{S}(\mathbf{u}_j, x_j) - \frac{\Delta x}{3} \mathbf{S}(\mathbf{u}_{j+1/4}, x_{j+1/4}) - \frac{\Delta x}{12} \mathbf{S}(\mathbf{u}_{j+\frac{1}{2}}, x_{j+\frac{1}{2}}) \\
&= \mathbf{f}_{j+\frac{1}{2}} - \mathbf{f}_j - \frac{\Delta x}{12} \mathbf{f}'_j - \frac{\Delta x}{3} (\mathbf{f}'_{j+1/4} - \mathbf{S}_u(\mathbf{u}_{j+1/4}, x_{j+1/4}) \frac{\Delta x^3}{128} \mathbf{u}'''') - \frac{\Delta x}{12} \mathbf{f}'_{j+\frac{1}{2}} \\
&= \mathbf{f}_{j+\frac{1}{2}} - \mathbf{f}_j - \frac{\Delta x}{2} \mathbf{f}'_j - \frac{\Delta x^2}{8} \mathbf{f}''_j - \frac{\Delta x^3}{48} \mathbf{f}'''_j - \frac{5\Delta x^4}{1728} \mathbf{f}''''_j \\
&\quad + \mathbf{S}_u(\mathbf{u}_{j+1/4}, x_{j+1/4}) \frac{\Delta x^4}{384} \mathbf{u}''' + \mathcal{O}(\Delta x^5) \\
&= \left( -\frac{5}{1728} \mathbf{f}''''_j + \mathbf{S}_u(\mathbf{u}_{j+1/4}, x_{j+1/4}) \frac{\mathbf{u}'''}{384} \right) \Delta x^4 + \mathcal{O}(\Delta x^5),
\end{aligned}
$$

which achieves the proof.

REFERENCES

[1] R. ABGRALL, A combination of residual distribution and the Active Flux formulations or a new class of schemes that can combine several writings of the same hyperbolic problem: application to the 1d Euler equations, Commun. Appl. Math. Comput., 5 (2023), pp. 370–402.

[2] R. ABGRALL, M. JIAO, Y. LIU, AND K. WU, Bound preserving Point-Average-Moment PolynomiAl-interpreted (PAMPA) scheme: one-dimensional case, Commun. Comput. Phys., (2025). Accepted.

[3] R. ABGRALL, M. JIAO, Y. LIU, AND K. WU, A novel and simple invariant-domain-preserving framework for PAMPA scheme: 1D case, SIAM J. Sci. Comput., (2025). Accepted.

[4] R. ABGRALL AND Y. LIU, A new approach for designing well-balanced schemes for the shallow water equations: a combination of conservative and primitive formulations, SIAM J. Sci. Comput., 46 (2024), pp. A3375–A3400.

[5] L. ARPAIA AND M. RICCHIUTO, r-adaptation for shallow water flows: conservation, well balancedness, efficiency, Commun. Fluids, 160 (2018), pp. 175–203.

[6] L. ARPAIA AND M. RICCHIUTO, Well balanced residual distribution for the ale spherical shallow water equations on moving adaptive meshes, J. Comput. Phys., 405 (2020).

[7] E. AUDUSSE, F. BOUCHUT, M.-O. BRISTEAU, R. KLEIN, AND B. PERTHAME, A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows, SIAM J. Sci. Comput., 25 (2004), pp. 2050–2065.

[8] P. AZERAD, J.-L. GUERMOND, AND B. POPOV, Well-balanced second-order approximation of the shallow water equation with continuous finite elements, SIAM J. Numer. Anal., 55 (2017), pp. 3203–3224.

[9] W. BARSUKOW AND J. BERBERICH, A well-balanced Active Flux method for the shallow water equations with wetting and drying, Commun. Appl. Math. Comput., 6 (2024), pp. 2385–2430.

[10] F. BEHZADI AND J. NEWMAN III, An exact source-term balancing scheme on the finite element solution of shallow water equations, Commun. Method. Appl. M., 359 (2020), p. 112662.

[11] S. BUSTO AND M. DUMBSER, A staggered semi-implicit hybrid finite volume/finite element scheme for the shallow water equations at all Froude numbers, Appl. Numer. Math., 175 (2022), pp. 108–132.

[12] C. CABALLERO-CÁRDENAS, M. J. CASTRO, C. CHALONS, T. MORALES DE LUNA, AND M. MUÑOZ-RUIZ, A semi-implicit fully exactly well-balanced relaxation scheme for the shallow water system, SIAM J. Sci. Comput., 46 (2024), pp. A2503–A2527.

[13] Y. CAO, A. KURGANOV, Y. LIU, AND R. XIN, Flux globalization based well-balanced path-conservative central-upwind schemes for shallow water models, J. Sci. Comput., 92 (2022), p. 69.

[14] M. J. CASTRO, J. A. LÓPEZ, AND C. PARÉS, Finite volume simulation of the geostrophic adjustment in a rotating shallow-water system, SIAM J. Sci. Comput., 31 (2008), pp. 444–477.

[15] W. S. CHANG, M. M. BAZUHAIR, F. ISMAIL, H. CHIZARI, AND M. H. H. ISHAK, Well-balanced energy-stable residual distribution methods for the shallow water equations with varying bottom topography, Comput. Math. Appl., 148 (2023), pp. 188–199.
[16] G. CHEN AND S. NOELLE, *A new hydrostatic reconstruction scheme based on subcell reconstructions*, SIAM J. Numer. Anal., 55 (2017), pp. 758-784.

[17] Y. CHENG, A. CHERTOCK, M. HERTY, A. KURGANOV, AND T. WU, *A new approach for designing moving-water equilibria preserving schemes for the shallow water equations*, J. Sci. Comput., 80 (2019), pp. 538-554.

[18] Y. CHENG AND A. KURGANOV, *Moving-water equilibria preserving central-upwind schemes for the shallow water equations*, Commun. Math. Sci., 14 (2016), pp. 1643-1663.

[19] A. CHERTOCK, S. CUI, A. KURGANOV, Ş. N. ÖZCAN, AND E. TADMOR, *Well-balanced schemes for the Euler equations with gravitation: conservative formulation using global fluxes*, J. Comput. Phys., 358 (2018), pp. 36-52.

[20] A. CHERTOCK, M. DUDZINSKI, A. KURGANOV, AND M. LUKÁČOVÁ-MEDVID'OVÁ, *Well-balanced schemes for the shallow water equations with Coriolis forces*, Numer. Math., 138 (2018), pp. 939-973.

[21] A. CHERTOCK, A. KURGANOV, X. LIU, Y. LIU, AND T. WU, *Well-balancing via flux globalization: Applications to shallow water equations with wet/dry fronts*, J. Sci. Comput., 90 (2022), pp. 1-21.

[22] M. CIALLELLA, L. MICALIZZI, V. MICHEL-DANSAC, P. ÖFFNER, AND D. TORLO, *A high-order, fully well-balanced, unconditionally positivity-preserving finite volume framework for flood simulations*, GEM—International Journal on Geomathematics, 16 (2025), p. 6.

[23] M. CIALLELLA, L. MICALIZZI, P. ÖFFNER, AND D. TORLO, *An arbitrary high order and positivity preserving method for the shallow water equations*, Comput. Fluids, 247 (2022), p. 105630.

[24] M. CIALLELLA, D. TORLO, AND M. RICCHIUTO, *Arbitrary high order WENO finite volume scheme with flux globalization for moving equilibria preservation*, J. Sci. Comput., 96 (2023), p. 53.

[25] M. J. C. DÍAZ, J. A. LÓPEZ-GARCÍA, AND C. PARÉS, *High order exactly well-balanced numerical methods for shallow water systems*, J. Comput. Phys., 246 (2013), pp. 242-264.

[26] M. DUMBSER, O. ZANOTTI, E. GABURRO, AND I. PESHKOV, *A well-balanced discontinuous Galerkin method for the first-order Z4 formulation of the Einstein-Euler system*, J. Comput. Phys., 504 (2024), p. 112875.

[27] S. GOTTLIEB, D. KETCHESON, AND C.-W. SHU, *Strong stability preserving Runge-Kutta and multistep time discretizations*, World Scientific, 2011.

[28] S. GOTTLIEB, C.-W. SHU, AND E. TADMOR, *Strong stability-preserving high-order time discretization methods*, SIAM review, 43 (2001), pp. 89-112.

[29] J.-L. GUERMOND, M. Q. DE LUNA, B. POPOV, C. KEES, AND M. FARTHING, *Well-balanced second-order finite element approximation of the shallow water equations with friction*, SIAM J. Sci. Comput., 40 (2018), pp. A3873-A3901.

[30] W. GUO, Z. CHEN, S. QIAN, G. LI, AND Q. NIU, *A new well-balanced finite volume CWENO scheme for shallow water equations over bottom topography*, Adv. Appl. Math. Mech., 15 (2023), pp. 1515-1539.

[31] A. HAIDAR, F. MARCHE, AND F. VILAR, *A posteriori finite-volume local subcell correction of high-order discontinuous Galerkin schemes for the nonlinear shallow-water equations*, J. Comput. Phys., 452 (2022), p. 110902.

[32] E. HAIRER, G. WANNER, AND S. P. NORSET, *Solving Ordinary Differential Equations I. Nonstiff problems*, Springer, Berlin, Heidelberg, 1993.

[33] J. HWANG AND J. JUNG, *Path-conservative positivity-preserving well-balanced finite volume WENO method for porous shallow water equations*, J. Comput. Phys., 490 (2023), p. 112321.

[34] M. KAZOLEA, C. PAÉS, AND M. RICCHIUTO, *Approximate well-balanced WENO finite difference schemes using a global-flux quadrature method with multi-step ODE integrator weights*, Computers & Fluids, 296 (2025), p. 106646.

[35] P. KNOBLOCH, D. KUZMIN, AND A. JHA, *Well-balanced convex limiting for finite element discretizations of steady convection-diffusion-reaction equations*, J. Comput. Phys., 518 (2024), p. 113305.

[36] A. KURGANOV, Y. LIU, AND V. ZEITLIN, *A well-balanced central-upwind scheme for the thermal rotating shallow water equations*, J. Comput. Phys., 411 (2020), p. 109414.

[37] G. LI, V. CALEFFI, AND Z. QI, *A well-balanced finite difference WENO scheme for shallow water flow model*, Appl. Math. Comput., 265 (2015), pp. 1-16.

[38] P. LI, W. DON, AND Z. GAO, *High order well-balanced finite difference WENO interpolation-based schemes for shallow water equations*, Comput. Fluids, 201 (2020), p. 104476.

[39] Y. LIU, *A well-balanced Point-Average-Moment PolynomiAl-interpreted (PAMPA) method for shallow water equations with horizontal temperature gradients*, SIAM J. Sci. Comput., (2025). Accepted.
[40] Y. LIU AND W. BARSUKOW, *An arbitrarily high-order fully well-balanced hybrid finite element-finite volume method for a one-dimensional blood flow model*, SIAM J. Sci. Comput., 47 (2025), pp. A2041–A2073.

[41] Y. LIU, J. LU, Q. TAO, AND Y. XIA, *An oscillation-free discontinuous Galerkin method for shallow water equations*, J. Sci. Comput., 92 (2022), p. 109.

[42] Y. MANTRI, ÖFFNER, AND M. RICCHIUTO, *Fully well-balanced entropy controlled discontinuous Galerkin spectral element method for shallow water flows: global flux quadrature and cell entropy correction*, J. Comput. Phys., 498 (2024), p. 112673.

[43] L. MICALIZZI, M. RICCHIUTO, AND R. ABGRALL, *Novel well-balanced continuous interior penalty stabilizations*, J. Sci. Comput., 100 (2024), p. 14.

[44] S. NOELLE, Y. XING, AND C.-W. SHU, *High-order well-balanced finite volume WENO schemes for shallow water equation with moving water*, J. Comput. Phys., 226 (2007), pp. 29–58.

[45] M. PENG, Z. SUN, AND K. WU, *OEDG: Oscillation-eliminating discontinuous Galerkin method for hyperbolic conservation laws*, Mathematics of Computation, 94 (2025), pp. 1147–1198.

[46] B. REN, Z. GAO, Y. GU, S. XIE, AND X. ZHANG, *A positivity-preserving and well-balanced high order compact finite difference scheme for shallow water equations*, Commun. Comput. Phys., 35 (2024), pp. 524–552.

[47] M. RICCHIUTO, *On the C-property and generalized C-property of residual distribution for the shallow water equations*, J. Sci. Comput., 48 (2011), pp. 304–318.

[48] M. RICCHIUTO, *An explicit residual based approach for shallow water flows*, J. Comput. Phys., 280 (2015), pp. 306–344.

[49] M. RICCHIUTO AND A. BOLLERMANN, *Stabilized residual distribution for shallow water simulations*, J. Comput. Phys., 228 (2009), pp. 1071–1115.

[50] B. WANG, P. LI, AND Z. GAO, *High-order well-balanced and positivity-preserving finite-difference AWENO scheme with hydrostatic reconstruction for shallow water equations*, Appl. Numer. Math., 181 (2022), pp. 483–502.

[51] Y. XING, *Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium*, J. Comput. Phys., 257 (2014), pp. 536–553.

[52] Y. XING AND C.-W. SHU, *High order finite difference WENO schemes with the exact conservation property for the shallow water equations*, J. Comput. Phys., 208 (2005), pp. 206–227.

[53] Y. XING AND C.-W. SHU, *High-order well-balanced finite difference WENO schemes for a class of hyperbolic systems with source terms*, J. Sci. Comput., 27 (2006), pp. 477–494.

[54] Y. XING AND C.-W. SHU, *High-order finite volume WENO schemes for the shallow water equations with dry states*, Adv. Water Resour., 34 (2011), pp. 1026–1038.

[55] Y. XING, X. ZHANG, AND C.-W. SHU, *Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations*, Adv. Water Resour., 33 (2010), pp. 1476–1493.

[56] Z. XU AND C.-W. SHU, *A high-order well-balanced discontinuous Galerkin method for hyperbolic balance laws based on the Gauss–Lobatto quadrature rules*, J. Sci. Comput., 101 (2024), p. 39.

[57] Z. XU AND C.-W. SHU, *A high-order well-balanced alternative finite difference WENO (AWENO) method with the exact conservation property for systems of hyperbolic balance laws*, Adv. Water Resour., 196 (2025), p. 104898.

[58] R. YANG, Y. YANG, AND Y. XING, *High order sign-preserving and well-balanced exponential Runge–Kutta discontinuous Galerkin methods for the shallow water equations with friction*, J. Comput. Phys., 444 (2021), p. 110543.

[59] Z. ZHANG, J. DUAN, AND H. TANG, *High-order accurate well-balanced energy stable adaptive moving mesh finite difference schemes for the shallow water equations with non-flat bottom topography*, J. Comput. Phys., 492 (2023), p. 112451.

[60] Z. ZHANG, H. TANG, AND K. WU, *High-order accurate structure-preserving finite volume schemes on adaptive moving meshes for shallow water equations: Well-balancedness and positivity*, J. Comput. Phys., 527 (2025), p. 113801.