# Arbitrary High Order WENO Finite Volume Scheme with Flux Globalization for Moving Equilibria Preservation

Mirco Ciallella$^{(1)}$; Davide Torlo$^{(2)†}$and Mario Ricchiuto$^{(3)‡}$

(1): École Nationale Supérieure d'Arts et Métiers, Institut de Mécanique et d'Ingénierie I2M, 33400 Talence, France

(2): SISSA mathLab, Mathematics Area, SISSA, via Bonomea 265, Trieste 34136, Italy

(3): INRIA, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, France

January 4, 2024

## Abstract

In the context of preserving stationary states, e.g. lake at rest and moving equilibria, a new formulation of the shallow water system, called Flux Globalization has been introduced by Cheng et al. (2019). This approach consists in including the integral of the source term in the *global flux* and reconstructing the new *global flux* rather than the conservative variables. The resulting scheme is able to preserve a large family of smooth and discontinuous steady state moving equilibria. In this work, we focus on an arbitrary high order WENO Finite Volume (FV) generalization of the global flux approach. The most delicate aspect of the algorithm is the appropriate definition of the *source flux* (integral of the source term) and the quadrature strategy used to match it with the WENO reconstruction of the hyperbolic flux. When this construction is correctly done, one can show that the resulting WENO FV scheme admits exact discrete steady states characterized by constant global fluxes. We also show that, by an appropriate quadrature strategy for the source, we can embed exactly some particular steady states, e.g. the lake at rest for the shallow water equations. It can be shown that an exact approximation of global fluxes leads to a scheme with better convergence properties and improved solutions. The novel method has been tested and validated on classical cases: subcritical, supercritical and transcritical flows.

*Keywords: Flux globalization, WENO, well-balanced, moving equilibria, shallow water.*

## 1 Introduction

The Saint-Venant or shallow water equations model the dynamics of hydrostatic free surface waves under the action of gravity [26]. This system of non-linear partial differential equations is valid under the hypothesis of very large wavelengths, or very shallow depth, and describes the evolution of the water depth and volume flux. They are used in a variety of engineering applications going from river and estuarine hydrodynamics, to urban flood and tsunami risk assessment. The basic set of equations constitutes a system of hyperbolic conservation laws, endowed with all the usual properties [56]: real eigenstructure of the flux Jacobians, Rankine-Hugoniot relations defining weak solutions, and a convex entropy extension to detect the admissible ones. The influence of spatial changes of the height of the bottom topography, of the viscous effects in boundary layers, as well as other physical effects, are accounted for by means of appropriately defined source terms. The addition of these terms 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. Analytical solutions are not always available for such states, but they can often

<footnote>*Corresponding author mirco.ciallella@ensam.eu</footnote>

<footnote>†davide.torlo@sissa.it</footnote>

<footnote>‡mario.ricchiuto@inria.fr</footnote>
be characterized in some implicit form (see e.g. [27]).

The numerical approximation of the shallow water equations is a very active research domain. It is quite impossible to mention all the relevant literature, numerous original methods having been devised in many different settings: finite volumes [6, 30, 47, 13], finite elements [38, 39, 66, 15, 60, 59, 9], residual distribution [55, 54, 53, 52, 4, 5], discontinuous Galerkin [63, 45, 6, 47] and so on. As in many other fields, the final goal is to design robust numerical methods, providing physically relevant solutions in realistic applications with the lowest possible computational effort. In this respect, very high order methods are known for their computational efficiency, especially when dealing with sufficiently smooth solutions. Among these, an excellent compromise between high order of accuracy and robustness for non-smooth solutions or data is provided by WENO finite volume methods [57, 58] which are the object of the present paper.

For a given accuracy order, and on a given mesh, a rewarding strategy to design numerical methods with reduced error is that of the so called property-preserving discretizations. These methods try to mimic additional consistency conditions than just those expressed by the system of equations themselves. A typical example of such derived conditions is the energy/entropy (in)equality. For the shallow water equations, the presence of the source terms makes this issue particularly interesting. As mentioned, the system of balance laws admits a large number of equilibrium solutions, depending on the interaction flux-variation/source. This notion is related to the notion of a well-balanced discretization, usually defined as one capable of reproducing one or more of these equilibria at the discrete level. The most classical version of this notion is related to the discrete conservation of the so called lake at rest state with flat free surface and zero velocity [11, 62, 35]. Since then, this issue is essential to the design of any numerical method for the shallow water equations. The numerical evidence shows that, besides being exact for the lake at rest state, well balanced methods also provide error reduction for other states. This has pushed the community toward the development of a more general notion of well balancedness, hoping to further improve the cost efficiency of the resulting schemes. These generalizations try to account for other equilibria than just the lake at rest state, as well as for more physical effects, and source terms.

There have been several approaches proposed to this end. Among these we mention the fully well balanced generalized Riemann solver proposed by Berthon and Chalons [13, 12], the optimization based corrected reconstruction method by Castro and Parés [18, 33, 34], the steady residual subtraction technique proposed in [10]. In the previously mentioned methods only local information is used to obtain discrete equilibria. Another approach consists of a global flux method first introduced in [32], then applied to design non-linear shock capturing methods for balance laws in [29, 16], and more recently as an alternative way to design fully well balanced methods [21, 19, 22]. The main idea of these techniques is to build a quasi-conservative hyperbolic system starting from a balance law by integrating the source term so that it can be included in a global definition of the flux. Some methods require the explicit knowledge of the exact equilibrium, some others are agnostic of the specific form of the equilibrium, but might need to solve nonlinear equations at each cell and time, or might not preserve the steady states in all configurations, e.g. shocks. In this work we focus on the flux globalization framework, which belong to the latter class. This approach relies on the definition of a single flux (a global flux) whose divergence expresses the effects of both the divergence of the conservative flux, and that of the forcing term. By standard techniques it is possible to devise a scheme which is exactly consistent with a constant global flux. This gives a generalized consistency condition which naturally embeds the effects of the source integral. In its simplest formulation, this approach is a natural way to embed in the finite volume setting the notion of residual distribution, as discussed in [3].

In this work we propose a possible strategy to combine flux globalization with very high order of accu-
racy via a WENO approach. The heart of the method is based on three main elements: a high accurate
quadrature allowing to define cell averages for the global flux; a global flux WENO reconstructions based
on these averages; a standard upwind numerical flux. We can easily show that this strategy guarantees the
exact consistency with constant cell-averages of the global flux. For the Saint-Venant equations with the
<figure><img src="image_1.png" /><figcaption>Figure 1: Shallow Water Equations: definition of the variables.</figcaption></figure>

bathymetric and friction terms, we consider further improvements to guarantee the exact preservation of the lake at rest solution. Although the method is not exactly well balanced with respect to all analytically computable equilibria, we show that for a wide range of solutions the overall gains in accuracy obtained with the proposed method are of several orders of magnitude on a given mesh.

The paper is organized as follows. In Section 2 we introduce the shallow water equations and the steady state equilibria that we want to preserve. In Section 3 we introduce the spatial discretization including the finite volume discretization, the global flux definition, the lake-at-rest well-balanced modification and the high order WENO reconstruction. In Section 4 we describe the temporal discretization, based on a Deferred Correction method. In Section 5 we perform several tests showing the superior capability of the global flux high order method with respect to other high order methods. Finally, in Section 6 we draw some conclusions and we describe some perspective future works.

# 2 Shallow Water Equations

## 2.1 Model

Although the method proposed can be used for more general systems of balance laws, in this paper we will focus on the shallow water equations (SWEs) with bathymetry and the friction source terms. For time independent topographies, the system reads

$$ \frac{\partial \mathbf{u}}{\partial t} + \frac{\partial \mathcal{F}(\mathbf{u})}{\partial x} = \mathcal{S}(\mathbf{u}, x), \quad \text{on} \quad \Omega_T = \Omega \times [0, T] \subset \mathbb{R} \times \mathbb{R}^+, \qquad (1) $$

with conserved variables, flux and source term given by

$$ \mathbf{u} = \begin{bmatrix} h \\ q \end{bmatrix}, \quad \mathcal{F}(\mathbf{u}) = \begin{bmatrix} q \\ \frac{q^2}{h} + g \frac{h^2}{2} \end{bmatrix}, \quad \mathcal{S}(\mathbf{u}, x) = \begin{bmatrix} 0 \\ S(\mathbf{u}, x) \end{bmatrix} = -gh \begin{bmatrix} 0 \\ \frac{\partial b(x)}{\partial x} \end{bmatrix} - gq \begin{bmatrix} 0 \\ \frac{n^2}{h^{7/3}} |q| \end{bmatrix}, \qquad (2) $$

where (cf. Figure 1) $h$ represents the relative water height, $q$ is its discharge (equal to $hu$, where $u$ is the vertically averaged velocity), $g$ is the gravitational acceleration, $b(x)$ is the local bathymetry and $n$ is the Manning friction coefficient. It is also convenient to introduce the free surface water level $\eta := h + b$.

The key idea used in this paper is to recast (1) in the following equivalent pseudo-conservative form:

$$ \frac{\partial \mathbf{u}}{\partial t} + \frac{\partial \mathcal{G}(\mathbf{u}, x)}{\partial x} = 0 \quad \text{such that} \quad \mathcal{G}(\mathbf{u}, x) = \begin{bmatrix} q \\ K \end{bmatrix} = \begin{bmatrix} q \\ \frac{q^2}{h} + g \frac{h^2}{2} + \mathcal{R} \end{bmatrix}, \qquad (3) $$

having set

$$ \mathcal{R}(x, t) := - \int_{x}^{x} S(\mathbf{u}, \xi) \, d\xi = g \int_{x}^{x} \left[ h(\xi, t) \frac{\partial b(\xi)}{\partial \xi} + \frac{n^2}{h^{7/3}(\xi, t)} |q(\xi, t)| q(\xi, t) \right] \, d\xi. \qquad (4) $$
In this case both components of the global flux, q and K, are steady equilibrium variables. Moreover, the form of system (3) allows the use of some classical techniques for conservation laws, being the source term be treated through the definition and the local discrete approximation of $\mathcal{R}$.

## 2.2 Steady state equilibria

As recalled in the introduction, system (2) admits a variety of steady equilibria dependent on the interaction between the variations of $\mathcal{F}$, or of some of its components, and the source. We are interested in well balanced schemes that are capable of preserving in some sense some of these equilibria at the discrete level. Classically, the most studied equilibrium is the lake at rest given by

$$u = 0; \quad \eta(x,t) = h(x,t) + b(x) \equiv \eta_0 \in \mathbb{R}, \quad \forall x \in \Omega, t \in [0, T]. \qquad (5)$$

this can also be considered as a special case of a constant energy state which is not at rest, defined as

$$h(x,t)u(x,t) \equiv q_0; \quad \mathcal{E} = g\eta(x,t) + \frac{u^2}{2} \equiv \mathcal{E}_0 \in \mathbb{R}, \quad \forall x \in \Omega, t \in [0, T]. \qquad (6)$$

There are several works treating specifically these two solutions, and other works deal with other analytical states, see e.g. [53, 48].

A more general standpoint is to consider that steady state solutions are generally characterized by the invariants

$$h(x,t)u(x,t) \equiv q_0 \quad \text{and} \quad K(x,t) = \frac{q_0^2}{h} + g\frac{h^2}{2} + \int_{x_0}^{x} g \left[ h(\xi, t) \frac{\partial b(\xi)}{\partial \xi} + \frac{n^2}{h^{7/3}(\xi, t)} |q(\xi, t)| q(\xi, t) \right] d\xi \equiv K_0. \qquad (7)$$

Obviously, the lake at rest (5) and constant energy (6) are a special cases of (7). For a smooth solution in the frictionless case, for example, we have

$$q \equiv q_0, \qquad (8)$$

$$0 = \partial_x \left( \frac{q^2}{h} + g \frac{h^2}{2} \right) + gh\partial_x b = - \frac{q_0^2}{h^2} \partial_x h + gh\partial_x h + gh\partial_x b = h\partial_x \left( \frac{q_0^2}{2h^2} + g(h+b) \right). \qquad (9)$$

in which the last reduces to the second in (6). Then, the last equation sums up to verify that $\Upsilon := \frac{q_0^2}{2h^2} + g\eta \equiv \Upsilon_0$ is constant in space and time. However, this condition might not hold when bathymetry or water height are discontinuous in some points, e.g. transcritical flows with shock, while (7) is always true.

## 3 Space discretization: Global Flux Finite Volume method

The global flux idea is discussed in some detail in [19] in the context of a piecewise linear approximation. In this section we discuss a possible strategy to construct arbitrarily high order extensions based on WENO reconstructions. To begin with, the computational domain $\Omega$ is discretized into $N_x$ equispaced control volumes $\Omega_i = [x_{i-1/2}, x_{i+1/2}]$ of size $\Delta x$ centered at $x_i = i\Delta x$ with $i = i_\ell, \dots, i_r$. For the control volume $\Omega_i$ we can define the cell average at time t:

$$\bar{U}_i(t) := \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} \mathbf{u}(x, t) \, dx. \qquad (10)$$

The semi-discrete finite volume scheme for the system (3) reads

$$\frac{d\bar{U}_i}{dt} + \frac{1}{\Delta x}(\hat{\mathbf{H}}_{i+1/2} - \hat{\mathbf{H}}_{i-1/2}) = 0, \qquad (11)$$

where $\hat{\mathbf{H}}_{i+1/2}$ is a numerical flux consistent with the global flux $\mathcal{G}$. The global flux differs from the original flux and this makes tricky the development of an upwind scheme based on the solution of Riemann problems
by an approximate solver. In this work we are thus focusing on a relatively simple approach in which upwinding is defined by the homogeneous system. In other words we take

$$ \hat{\mathbf{H}}_{i+1/2} = (L^{-1}\Lambda^{+}L)_{i+1/2}\mathcal{G}_{i+1/2}^{L} + (L^{-1}\Lambda^{-}L)_{i+1/2}\mathcal{G}_{i+1/2}^{R}. \quad (12) $$

Here, $\mathcal{G}_{i+1/2}^{L,R}$ are the discontinuous reconstructed point values of the global flux $\mathcal{G}(\mathbf{U})$ respectively at the left and right side of the cell interface $x_{i+1/2}$. $L$ is the matrix of the left eigenvectors computed from the flux Jacobian of the hyperbolic problem (1) in the Roe averaged state, i.e.,

$$ J(\mathbf{U}^*) = \begin{pmatrix} 0 & 1 \\ -u_*^2 + gh_* & 2u_* \end{pmatrix}, \quad \text{with} \quad \begin{cases} h_* = \frac{h^L+h^R}{2}, \\ u_* = \frac{\sqrt{h^L u^L + \sqrt{h^R u^R}}}{\sqrt{h^L + \sqrt{h^R}}}. \end{cases} \quad (13) $$

$\Lambda^{\pm}$ correspond to the upwinding weights

$$ \Lambda_i^+ = \begin{cases} 1, & \text{if } \lambda_i > 0, \\ 0, & \text{if } \lambda_i < 0, \end{cases} \quad \Lambda_i^- = \begin{cases} 1, & \text{if } \lambda_i < 0, \\ 0, & \text{if } \lambda_i > 0. \end{cases} \quad (14) $$

It is important to notice that the Jacobian of the flux cannot be directly computed as not all the quantities are available at the interface (we reconstruct only the global flux, not the conserved quantities). Hence, one has to recover the value of $h$ from the global flux itself. To do so, in a positive manner, we use the technique proposed in [19] and recalled for completeness in Appendix A.

Other Riemann solvers, e.g. Rusanov, can be used as long as the numerical dissipation term depending on the difference of the states at the interface vanishes when the steady states are reached. For example, this can be achieved using cutoff functions [21].

In order to obtain the reconstructed values at the left and right side of interfaces, we will use a high order WENO reconstruction technique that will be applied directly to cell averages of the global flux. As we will see, this will require introducing WENO polynomials for all the quantities of interests, i.e., $h, q, b, \mathcal{R}, \mathcal{G}$ and $K$. Note in particular that the simple averaging technique for $\mathcal{R}$ used in [19] is not enough to reach more than second order of accuracy. Another delicate point is the definition of the local values of the source term used to obtain $\mathcal{R}$. We deal with this aspect in Section 3.3.

## 3.1 Global flux evaluation and cell quadrature

As a starting point for the WENO procedure, we need cell averages of the global flux, denoted here by

$$ \tilde{\mathcal{G}}_i(\mathbf{u}, x) = \bar{\mathcal{F}}_i(\mathbf{u}) + \begin{bmatrix} 0 \\ \bar{\mathcal{R}}_i \end{bmatrix}. \quad (15) $$

The first step in our discretization is to build a WENO reconstruction of the conserved variables, hereafter denoted by $\tilde{\mathbf{u}}(x)$, which is used as a basis for the quadrature formula

$$ \bar{\mathcal{F}}_i(\mathbf{u}) = \sum_q w_q \mathcal{F}(\tilde{\mathbf{u}}(x_{i,q})), \quad (16) $$

where $(x_{i,q}, w_q)$ are the high order quadrature points and weights used in the cell $[x_{i-1/2}, x_{i+1/2}]$. In particular, we will use a 4-point Gauss-Legendre quadrature. A similar formula is used for the integral source $\mathcal{R}$:

$$ \bar{\mathcal{R}}_i \approx \sum_q w_q \mathcal{R}_{i,q}. \quad (17) $$

To obtain the values at the quadrature points, we suppose a piecewise polynomial reconstruction of the integral source $\mathcal{R}$. Exploiting the relation $\partial_x \mathcal{R} \equiv S$ and we can evaluate the required terms as

$$ \mathcal{R}_{i,q} = \mathcal{R}_{i-1/2}^R - \int_{x_{i-1/2}^R}^{x_{i,q}} \tilde{S}(x) \, dx = \mathcal{R}_{i-1/2}^R - \sum_{\theta} \underbrace{\int_{x_{i-1/2}^R}^{x_{i,q}} \ell_{\theta}(x) \, dx}_{r_{\theta}^q} S(x_{i,\theta}), \quad i \ge i_l, \quad (18) $$
where $\ell_\theta$ are the Lagrangian polynomials associated to the quadrature points, $r_\theta^q$ are computed exactly, and having implicitly assumed that

$$ \tilde{S}(x)|_{[x_{i-1/2}, x_{i+1/2}]} = \tilde{S}_i(x) = \sum_{\theta} \ell_{\theta}(x) S(x_{i,\theta}). \quad (19) $$

We remark that $q$ and $\theta$ are just free indexes and that both $x_{i,q}$ and $x_{i,\theta}$ are the quadrature nodes on the interval $[x_{i-1/2}, x_{i+1/2}]$. Moreover, to maintain the well balancing with respect to lake at rest, a further splitting of the source term has to be done and it will be described in details in Section 3.3, but, for the moment, we maintain this notation to better convey the main idea of the scheme.

The definition of the array of quadrature values of $\mathcal{R}$ requires an appropriate initial value $\mathcal{R}_{i-1/2}^R$. In practice, we impose that $\mathcal{R}$ at the left extrema of the computational domain is equal to 0, i.e.,

$$ \mathcal{R}_{i-1/2}^R = 0. \quad (20) $$

Then, when iterating over the elements to evaluate (18) $\forall i$, we set

$$ \mathcal{R}_{i+1/2}^R = \mathcal{R}_{i+1/2}^L + [\mathcal{R}_{i+1/2}], \quad (21) $$

where we note that $\mathcal{R}_{i+1/2}^L$ is

$$ \mathcal{R}_{i+1/2}^L = \mathcal{R}_{i-1/2}^R - \int_{x_{i-1/2}^R}^{x_{i+1/2}^L} \tilde{S}(x) \, dx = \mathcal{R}_{i-1/2}^R - \Delta x \tilde{S}_i, \quad i > 1, \quad (22) $$

with the average source $\bar{S}_i$ obtained as

$$ \bar{S}_i := \frac{1}{\Delta x} \int_{x_{i-1/2}^R}^{x_{i+1/2}^L} \tilde{S}(x) \, dx = \sum_{\theta} w_{\theta} S(x_{i,\theta}). \quad (23) $$

Notice that we need to introduce the jump of $\mathcal{R}$ to be able to preserve the lake at rest steady state because the bathymetry is discontinuous at the interfaces. Indeed, for continuous reconstruction of the bathymetry this term vanishes.

Finally, we need to provide a precise definition of how the cell average of the source $\bar{S}_i$ are evaluated, as well as of the jumps $[\mathcal{R}_{i+1/2}]$. Concerning the latter, for collocated approximations including the boundary nodes, as those used in [46], this jump could be set to zero. For the non-collocated approximation used here this is not necessarily adequate. For the cell averages $\tilde{S}_i$ a special treatment has to be applied to maintain equilibrium for lake at rest simulations. These aspects are discussed in Section 3.3.

## 3.2 Weighted Essentially Non-Oscillatory (WENO) reconstruction

The WENO reconstruction is used in this work for several quantities:

* for the solution and for the data, in order to be able to define the fluxes and the source values at quadrature points in (16), (18) and (23)

* for the global flux $\mathcal{G}(x) = \mathcal{F}(x) + \mathcal{R}(x)$ in order to obtain the left and right interface values for the upwind scheme in (12).

For each of these quantities, high order approximations are constructed starting from cell averages. We briefly recall the basics of the WENO reconstruction to obtain the polynomial $u_i(x)$ in cell $i$. We consider polynomial reconstructions of order $p$, with $p$ odd. To construct them, we select stencils of $p$ cells around cell $i$:

$$ \{\Omega_{l_x}, \quad l_x = i - r + 1, \dots, i + r - 1\}, \quad (24) $$
where $2r - 1 = p$. On each of these stencils, one constructs a high order polynomial $p^{HO}$ fulfilling the constraints

$$ \frac{1}{\Delta x} \int_{x_{i-j-1/2}}^{x_{i-j+1/2}} p^{HO}(x) dx = u_{i-j}, \quad j = -r+1, \dots, r-1, \qquad (25) $$

and $r$ low order polynomials $p_m(x)$, $m = 0, \dots, r-1$, that fulfill

$$ \frac{1}{\Delta x} \int_{x_{i-r+j+m-1/2}}^{x_{i-r+j+m+1/2}} p_m(x) dx = u_{i-j+m}, \quad j = 1, \dots, r. \qquad (26) $$

The WENO reconstruction aims at combining the low order polynomials in order to obtain the high order reconstruction in case all the low order polynomials are non-oscillatory, while it will prefer the least oscillatory polynomial in case some of these polynomials show oscillations. As an example, the WENO reconstruction of order 5 (WENO5) will use a stencil of length $p = 5$ with $r = 3$ low order reconstructions. Therefore, the involved cells span from $i-2$ to $i+2$.

To achieve high order accuracy optimal linear weights $d_m$ can be defined such that $\sum_m d_m(x)p_m(x) = p^{HO}(x)$ (see e.g. [40, 8]). The linear approximation obtained in this way however suffers from oscillations and Gibbs phenomena. To remove these artifacts, non-linear weights are introduced, defined $\forall m = 0, \dots, r-1$ as

$$ \omega_m = \frac{\alpha_m}{\sum_{k=0}^{r-1} \alpha_k}, \quad \alpha_k = \frac{d_k}{(\beta_k + \epsilon)^2}. \qquad (27) $$

In the last expression $\epsilon$ is a small number controlling the variations detected by the weights, and needed to avoid division by zero ($10^{-6}$ in this paper), while the $\beta_k$ are smoothness indicators defined by

$$ \beta_k = \sum_{l=1}^{r-1} \int_{x_{i-1/2}}^{x_{i+1/2}} \left( \frac{d^l}{dx^l} p_k(x) \right)^2 \Delta x^{2l-1} dx, \quad k = 0, \dots, r-1. \qquad (28) $$

The WENO approximation is finally defined as:

$$ \tilde{u}(x) = \sum_{m=0}^{r-1} \omega_m(x) p_m(x). \qquad (29) $$

The values of the optimal weights $d_m$ and the formulae for computing $\beta_m$ can be found in [40, 8] up to $r = 6$. In the following, we will test the reconstruction with orders $p = 3$ and $p = 5$.

**Proposition 3.1** (Global flux property and steady states). Given a state $\{\bar{U}_i\}_{i=i_l}$ such that $\bar{G}_i = G_0$ for all $i$, then $\{\bar{U}_i\}_{i=i_l}$ is a steady state of scheme (11) with a WENO reconstruction of the global flux $G$.

*Proof.* Since every low order polynomial interpolates the constant values $G_0$, all polynomials are identically the constant $G_0$. In this case, for the upwind flux (12), we trivially get $\bar{H}_{i\pm 1/2} = G_0$ and, hence, $\frac{d\bar{U}_i}{dt} = 0$, i.e., the thesis. $\square$

**Remark 3.2** (ODE integration analogy at equilibrium). The exact solution of the global flux scheme is not (or not necessarily) the analytical solution of the PDE. However, there exists a discrete relation of the steady state equilibrium solution of the global flux method, but it is highly nonlinear due to the high order finite volume formulation and the WENO nonlinear reconstruction. Indeed, at the equilibrium, the PDE (3) becomes an ODE and an analogy to classical ODE solvers may be made considering our discrete solution. In practice, due to WENO reconstruction, the presented global flux method can be recast into a sort of multistep method, where the global flux equations depend on the left and right neighboring cell averages.

The analysis is simpler for other approximation choices, as shown for example in [46] where the same ideas have been developed in the framework of a DG-SEM discretization using Gauss-Lobatto interpolation points. In this case, an explicit equivalence with collocation Runge-Kutta ODE solvers can be shown, and more precisely with the well known LobattoIIIA method. In this case, super-convergence results can be applied to precisely characterize the discrete steady states. We have not been able so far to provide similar explicit analogies for the method proposed here, due to the intricacy of the WENO reconstruction, combined with Gaussian quadrature of the source term.
### 3.3 Source term quadrature for the shallow water equations

For the shallow water equations, we provide here explicit formulas for the evaluation of the source term at quadrature points. We also discuss the issue of the jumps at the interface. We focus mostly on the frictionless case. The treatment of friction is discussed at the end of the paragraph.

We start by defining the WENO reconstructions of $h$, $\eta$ and $b$ in the quadrature points of the cell $\Omega_i$, denoting them by $\tilde{h}_{i,q}$, $\tilde{\eta}_{i,q}$ and $\tilde{b}_{i,q}$. We assume unduly that the same weights are used for $h$, and for the bathymetry $b$, so that the consistency with constants water level is trivially satisfied for the reconstructions $\eta_0 \equiv \tilde{\eta}_{i,q} = \tilde{h}_{i,q} + \tilde{b}_{i,q}$. In practice, the weights are directly computed using $\eta$ to detect this state. Then, we obviously have that, at each quadrature point,

$$ \tilde{h}_{i,q} = \tilde{\eta}_{i,q} - \tilde{b}_{i,q}. \qquad (30) $$

Let us also define a Lagrange interpolation of the bathymetry inside the $i$-th cell and its evaluation at the interfaces:

$$ \tilde{b}_i(x) := \sum_q \ell_q(x) \tilde{b}_{i,q}, \quad \text{and} \quad b_{i+1/2}^L = \tilde{b}_i(x_{i+1/2}), \quad b_{i-1/2}^R = \tilde{b}_i(x_{i-1/2}). \qquad (31) $$

We start by looking at the bathymetry term source which is the most critical as it contains a derivative in space. Following [64] we start by writing this term as

$$ S(\mathbf{u}, x) = -gh(x)\partial_x b(x) = -g\eta(x)\partial_x b(x) + g\partial_x \left(\frac{b^2(x)}{2}\right). \qquad (32) $$

To evaluate the integral of the source we now sample the WENO polynomials of $\eta$ and $b$ in the quadrature points. More precisely, we compute the required values of $\mathcal{R}$ as

$$ \mathcal{R}_{i,q} = \mathcal{R}_{i-1/2}^R - \int_{x_{i-1/2}^R}^{x_{i,q}} S(\mathbf{u}(x), x) \, dx \qquad (33) $$

$$ = \mathcal{R}_{i-1/2}^R + g \int_{x_{i-1/2}^R}^{x_{i,q}} \eta(x) \partial_x b(x) \, dx - g \left( \frac{(b_{i,q})^2}{2} - \frac{(b_{i-1/2}^R)^2}{2} \right) \qquad (34) $$

where in each quadrature point $x_{i,q}$ now we set

$$ \partial_x b(x_{i,q}) = \sum_s \ell'_s(x_{i,q}) b(x_{i,s}). $$

Now, using a Lagrange interpolation of the whole source term in the quadrature points, we obtain

$$ \mathcal{R}_{i,q} = \mathcal{R}_{i-1/2}^R + g \int_{x_{i-1/2}^R}^{x_{i,q}} \sum_{\theta} \ell_{\theta}(x) \tilde{\eta}_{i,\theta} \sum_{s} \ell'_{s}(x_{\theta}) \tilde{b}_{i,s} \, dx - g \left( \frac{(\tilde{b}_{i,q})^2}{2} - \frac{(\tilde{b}_{i-1/2}^R)^2}{2} \right). \qquad (35) $$

Similarly, we can define the cell averages (23) of the source as

$$ \bar{S}_i = \frac{1}{\Delta x} \left[ g \int_{x_{i-1/2}^R}^{x_{i+1/2}^L} \sum_{\theta} \ell_{\theta}(x) \tilde{\eta}_{i,\theta} \sum_{s} \ell'_{s}(x_{\theta}) \tilde{b}_{i,s} \, dx - g \left( \frac{(b_{i+1/2}^L)^2}{2} - \frac{(b_{i-1/2}^R)^2}{2} \right) \right], \qquad (36) $$

so that, the left interface terms, introduced in (22), can be defined as

$$ \mathcal{R}_{i+1/2}^L = \mathcal{R}_{i-1/2}^R - \Delta x \bar{S}_i. \qquad (37) $$
Finally, we define the jump of $\mathcal{R}$ across the interfaces. The objective is to obtain exactness for lake at rest state, i.e., $K$ has to be constant in this equilibrium. For more details see Appendix B and Remark B.1. This leads to

$$ \llbracket \mathcal{R}_{i+1/2} \rrbracket := g \frac{\eta_{i+1/2}^R + \eta_{i+1/2}^L}{2} (b_{i+1/2}^R - b_{i+1/2}^L) - g \left( \frac{(b_{i+1/2}^R)^2}{2} - \frac{(b_{i+1/2}^L)^2}{2} \right). \quad (38) $$

This strategy is similar to classical strategies used in well-balanced path conservative methods [17], where one uses a linear/segment path to connect the left and right states when evaluating the integral. These definitions allow to easily prove the following property.

**Proposition 3.3 (Lake at rest preservation).** The global flux WENO finite volume scheme with quadrature (35) of the bathymetric source, and with definition (38) of the source integral jump at the interface is exactly well balanced for the lake at rest state.

*Proof.* See appendix B. □

**Remark 3.4 (Global Flux with jump of $\mathcal{R}$).** Being the jump of $\mathcal{R}$ included in the definition of the global flux $\mathcal{G}$, it does not affect its preservation when moving equilibria are considered, as Proposition 3.1 is still valid.

To conclude, when adding the friction term we only account for its contributions to the cell integrals. We thus set in the general case

$$ \mathcal{R}_{i,q} = \mathcal{R}_{i-1/2}^R + g \sum_{\theta} \int_{x_{i-1/2}^R}^{x_{i,q}} \ell_{\theta}(x) \, dx \left( \tilde{\eta}_{i,\theta} \sum_s \ell'_s(x_{i,\theta}) \tilde{b}_{i,s} + g \frac{\tilde{q}_{i,\theta} |\tilde{q}_{i,\theta}| n^2}{\tilde{h}_{i,\theta}^{7/3}} \right) - g \left[ \frac{(b_{i,\theta})^2}{2} - \frac{(b_{i-1/2}^R)^2}{2} \right] \quad (39) $$

and

$$ \mathcal{R}_{i+1/2}^L = \mathcal{R}_{i-1/2}^R + g \sum_{\theta} \int_{x_{i-1/2}^R}^{x_{i+1/2}^R} \ell_{\theta}(x) \, dx \left( \tilde{\eta}_{i,\theta} \sum_s \ell'_s(x_{i,\theta}) \tilde{b}_{i,s} + g \frac{\tilde{q}_{i,\theta} |\tilde{q}_{i,\theta}| n^2}{\tilde{h}_{i,\theta}^{7/3}} \right) - g \left[ \frac{(b_{i+1/2}^L)^2}{2} - \frac{(b_{i-1/2}^R)^2}{2} \right]. \quad (40) $$

In Algorithm 1 we summarize the steps of the reconstruction of the source integral.

<caption>**Algorithm 1 Source integral reconstruction**</caption>

```
R_{i,l-1/2} := 0
for i = i_l, ..., i_r do
    Reconstruct the variables h, η and b in each quadrature point θ of the cell, obtaining ̂h_{i,θ}, ̃η_{i,θ} and ̂b_{i,θ}
    using the same WENO weights (computed for η)
    Reconstruct q in the quadrature points obtaining ̃q_{i,θ}
    Define R_{i,q} as in (39)
    Define R_{i+1/2}^L as in (40)
    Define R_{i+1/2}^R as in (38)
    Define R_{i+1/2}^L := R_{i+1/2}^R + R_{i+1/2}
end for
```

# 4 Time discretization

Time integration can be performed with any time discretization method that fits the method of line approach, e.g. RK schemes. In this paper, we employ a Deferred Correction (DeC) method because it is a family of one step methods with arbitrarily high order of accuracy. The original DeC formulation was introduced in [25], then developed and studied in its different forms in [30, 50, 23, 44]. A slightly different
form was presented in [1] for applications to finite element methods. In the reference the DeC is presented as an iterative procedure that involves two operators. The iteration process mimic the Picard-Lindelöf proof at the discrete level with a fixed-point iterative method. Each iteration aims at gaining one order of accuracy, so that the order of accuracy sought can be reached with a finite number of corrections.

<figure><img src="image_2.png" /><figcaption>Figure 2: Time interval divided into sub-time steps</figcaption></figure>

To describe the basic form of the method consider the ODE

$$y'(t) = f(y(t)), \quad y(t_0) = y_0, \tag{41}$$

where $y: \mathbb{R} \to \mathbb{R}^S$ and $f: \mathbb{R}^S \to \mathbb{R}^S$, and in our setting $f$ is the semidiscretized operator from (11). As usual, the temporal domain is discretized in time steps $[t^n, t^{n+1}]$. We also introduce $M$ sub-time steps $\{[t^{n,m-1}, t^{n,m}]\}_{m=1}^M$ where the boundary points coincide with the extrema of the time step, i.e., $t^n = t^{n,0}$ and $t^{n,M} = t^{n+1}$. The DeC method exploits predictions of the solution at each sub-time nodes $t^{n,m}$ denoted as $y^{n,m}$ as explained in Figure 2.

We proceed by first defining a low order approximation of the integral version of the ODE (41), denoted by $\mathcal{L}^1$. The requirement is that the solution of $\mathcal{L}^1(y) = 0$ should be easily obtained. We can for example consider the explicit Euler method. We then consider a high order discretization of the ODE, which may be costly to solve, as e.g. in high order full tableau implicit RK methods. We denote this operator by $\mathcal{L}^2$, and define it as the collocated multi-stage method

$$\mathcal{L}^2(y^{n,0}, \dots, y^{n,M}) := \begin{cases} y^{n,M} - y^{n,0} - \Delta t \sum_{r=0}^M \theta_r^M f(y^{n,r}), \\ \vdots \\ y^{n,1} - y^{n,0} - \Delta t \sum_{r=0}^M \theta_r^1 f(y^{n,r}), \end{cases} \approx \begin{cases} y^{n,M} - y^{n,0} - \int_{t^{n,0}}^{t^{n,M}} f(y(s))ds, \\ \vdots \\ y^{n,1} - y^{n,0} - \int_{t^{n,0}}^{t^{n,1}} f(y(s))ds, \end{cases} \tag{42}$$

where $\theta_r^m$ is the integral of the $r$-th Lagrangian basis function defined on the sub-time nodes over the interval $[t^{n,0}, t^{n,m}]$ [37, 1] normalized by the factor $\Delta t$. This operator is a high order discretization of the integral form of the ODE in each sub-time step. Depending on the chosen sub-time nodes, the order of accuracy varies, for example, with equispaced nodes we can obtain a scheme with order $M+1$, while with Gauss-Lobatto nodes we obtain a $2M$-th accurate scheme, i.e., the Lobatto IIIA schemes [36]. The system of equations $\mathcal{L}^2 = 0$ is often a strongly coupled system of nonlinear algebraic equations. This is certainly the case when $f$ is the nonlinear discretization of a balance law as the shallow water equations.

Concerning the low order operator $\mathcal{L}^1$, it consists of an explicit Euler step at each sub-time step:

$$\mathcal{L}^1(y^{n,0}, \dots, y^{n,M}) := \begin{cases} y^{n,M} - y^{n,0} - \beta^M \Delta t f(y^{n,0}), \\ \vdots \\ y^{n,1} - y^{n,0} - \beta^1 \Delta t f(y^{n,0}), \end{cases} \tag{43}$$

where $\beta^m = (t^{n,m} - t^{n,0})/\Delta t$.

The DeC iterations provide a new value for the array of stage values $\mathbf{y}^{(k)} := (y^{n,0}, \dots, y^{n,M})^{(k)}$, given a prediction $\mathbf{y}^{(k-1)} := (y^{n,0}, \dots, y^{n,M})^{(k-1)}$. The method is first initiated by setting all entries to the last
<figure><img src="image_3.png" /><figcaption>Figure 3: Lake at rest solution: $\eta$ (red) and $b$ (blue).</figcaption></figure>

available solution: $\mathbf{y}^0 = (y(t^n), \dots, y(t^n))$. The new iterates are then computed by evaluating $\forall k = 1, \dots, K$

$$ \mathcal{L}^1(\mathbf{y}^{(k)}) = \mathcal{L}^1(\mathbf{y}^{(k-1)}) - \mathcal{L}^2(\mathbf{y}^{(k-1)}) \quad (44) $$

where $K$ is the final number of iterations. In particular, after $K$ iterations the order of accuracy of the method will be the minimum between $K$ and the accuracy of $\mathcal{L}^2$. As an example, to obtain order 5, we need $K = 5$ and for equispaced sub-time nodes $M = 4$ (for Gauss-Lobatto nodes $M = 3$ suffices). The interested reader, can refer to [37, 2, 61] for more details and recent developments.

# 5 Numerical Simulations

The arbitrary high order well-balanced WENO finite volume scheme based on Flux Globalization has been tested and validated on several test cases to assess convergence properties and performances. For comparison, we have used a standard WENO FV scheme with Rusanov numerical flux [57], without any global flux or other property-preserving features, with the source term computed using a high order Gauss-Legendre quadrature formula starting from the analytical formulation of the bathymetry and its derivative and the WENO reconstruction of $h$. We have also investigated the effect of having the global flux property without including the well-balancedness for lake at rest, by computing the source term with a high order Gauss-Legendre quadrature formula starting from the analytical formulation of the bathymetry and its derivative and the WENO reconstruction of $h$. In this formulation, all the reconstructions are performed with WENO weights obtained by the reconstructed variables themselves, not by the smoothness indicators of the free surface level $\eta$. Still, the global flux is assembled as in (15). We do not compare the presented global flux methods with still-water well-balanced methods, for example the one presented in [65]. Indeed, as observed in [66], such methods do not provide better performances on moving-water equilibria. Both WENO3 and WENO5 schemes are tested. All numerical simulations have been run at CFL = 0.5.

The benchmarks used are quite classical in literature. They allow on one hand to verify first the correct implementation and convergence of the WENO method, even the non well-balanced. On the other, the results prove the dramatic increase in accuracy obtained with the global flux method for cases at rest, moving, with/out friction, and for both continuous and discontinuous bathymetries and solutions.

## 5.1 Lake at rest

Firstly, we consider the lake-at-rest solution characterized by the initial data

$$ h(x, 0) = 1 - b(x), \quad q(x, 0) \equiv 0, \quad (45) $$
<table><caption>Table 1: Lake at rest: errors and estimated order of accuracy (EOA) with WB and non-WB schemes, using WENO3 and WENO5 reconstructions.</caption><thead><tr><th rowspan="3">N<sub>e</sub></th><th colspan="4">Non-WB</th><th colspan="4">WB</th></tr><tr><th colspan="2">h</th><th colspan="2">q</th><th colspan="2">h</th><th colspan="2">q</th></tr><tr><th>L<sub>2</sub> error</th><th>EOA</th><th>L<sub>2</sub> error</th><th>EOA</th><th>L<sub>2</sub> error</th><th>EOA</th><th>L<sub>2</sub> error</th><th>EOA</th></tr></thead><tbody><tr><th colspan="9">GF-WENO3</th></tr><tr><td>25</td><td>1.0384E-4</td><td>-</td><td>4.7943E-5</td><td>-</td><td>9.8858E-14</td><td>-</td><td>1.2228E-15</td><td>-</td></tr><tr><td>50</td><td>1.5496E-5</td><td>2.67</td><td>9.2488E-6</td><td>2.31</td><td>9.8667E-14</td><td>-</td><td>1.4249E-15</td><td>-</td></tr><tr><td>100</td><td>1.2117E-6</td><td>3.62</td><td>3.6777E-7</td><td>4.59</td><td>9.8276E-14</td><td>-</td><td>1.6041E-15</td><td>-</td></tr><tr><td>150</td><td>2.6776E-7</td><td>3.69</td><td>1.5898E-7</td><td>2.05</td><td>1.9644E-13</td><td>-</td><td>3.3908E-15</td><td>-</td></tr><tr><td>200</td><td>9.6323E-8</td><td>3.53</td><td>7.6469E-8</td><td>2.53</td><td>1.9619E-13</td><td>-</td><td>3.6713E-15</td><td>-</td></tr><tr><td>400</td><td>8.2671E-9</td><td>3.53</td><td>6.0441E-9</td><td>3.65</td><td>2.9360E-13</td><td>-</td><td>6.1689E-15</td><td>-</td></tr><tr><td>800</td><td>6.8811E-10</td><td>3.58</td><td>4.7122E-10</td><td>3.67</td><td>5.8655E-13</td><td>-</td><td>1.3035E-14</td><td>-</td></tr><tr><th colspan="9">GF-WENO5</th></tr><tr><td>25</td><td>5.1800E-5</td><td>-</td><td>6.1657E-5</td><td>-</td><td>9.8947E-14</td><td>-</td><td>1.3247E-15</td><td>-</td></tr><tr><td>50</td><td>4.4066E-6</td><td>3.45</td><td>1.5244E-6</td><td>5.18</td><td>9.8661E-14</td><td>-</td><td>1.4060E-15</td><td>-</td></tr><tr><td>100</td><td>6.7998E-7</td><td>2.66</td><td>3.5908E-7</td><td>2.06</td><td>9.8289E-14</td><td>-</td><td>1.5992E-15</td><td>-</td></tr><tr><td>150</td><td>1.5437E-7</td><td>3.63</td><td>8.8535E-8</td><td>3.42</td><td>1.9639E-13</td><td>-</td><td>3.4157E-15</td><td>-</td></tr><tr><td>200</td><td>4.1973E-8</td><td>4.50</td><td>2.3725E-8</td><td>4.55</td><td>1.9611E-13</td><td>-</td><td>3.7034E-15</td><td>-</td></tr><tr><td>400</td><td>1.3952E-9</td><td>4.89</td><td>7.5991E-10</td><td>4.95</td><td>2.9357E-13</td><td>-</td><td>6.2007E-15</td><td>-</td></tr><tr><td>800</td><td>4.3120E-11</td><td>5.01</td><td>2.2633E-11</td><td>5.06</td><td>5.8648E-13</td><td>-</td><td>1.3039E-14</td><td>-</td></tr></tbody></table>

over the computational domain $[0, 25]$ with subcritical inlet/outlet at the two boundaries. The bathymetry employed contains a sinusoidal bump damped at the boundaries (see Figure 3 for better visualization of solution and bathymetry) and it reads

$$b(x) = 0.05 \sin(x - 12.5) \exp(1 - (x - 12.5)^2). \qquad (46)$$

Let us remark that the classical bathymetry provided in [27] is only piecewise polynomial, but globally only $C^0$. Hence, it is not suited to test the accuracy of very high order methods. On the other side, the bathymetry (46) is $C^\infty$ and it has values smaller than machine precision at the boundaries. The gravitational constant is considered to be $g = 1$ and the simulation is run until the final time $T = 1$ with $N_e = \{25, 50, 100, 150, 200, 400, 800\}$ uniform cells using both the non-well-balanced (non-WB) and well-balanced (WB) version of the GF-WENO algorithm to assess the convergence and well-balancing properties. The convergence tests performed with WENO3 and WENO5 reconstructions are listed in Table 1. It can be noticed that the error decay for the non-WB simulations matches the order of the reconstruction for both GF-WENO3 and GF-WENO5, while, for the WB cases, the scheme is able to preserve the exact solution up to machine precision.

## 5.2 Small perturbation of the lake-at-rest solution

For this test case, we analyze the perturbation of the lake-at-rest solution characterized by

$$h(x, 0) = 1 - b(x) + \begin{cases} \alpha\psi(x), & \text{if } 9 < x < 10 \\ 0, & \text{otherwise} \end{cases}, \quad q(x, 0) \equiv 0 \qquad (47)$$

with $\psi(x)$ a perturbation function defined by

$$\psi(x) := \exp \left( 1 - \frac{1}{(1 - r(x))^2} \right), \quad \text{with } r(x) := 4(x - 9.5)^2 \qquad (48)$$

and $\alpha = 10^{-4}$ over the computational domain $[0, 25]$ with subcritical inlet/outlet at the two boundaries. The bathymetry is a rescaling of (46) and it is defined as

$$b(x) = 0.5 \sin(x - 12.5) \exp(1 - (x - 12.5)^2). \qquad (49)$$

A slightly different bathymetry, with respect to Section 5.1, has been chosen in order to introduce more noise in the non-WB simulation and appreciate more the method capabilities. The simulation was run using a
<figure><img src="image_4.png" /><figcaption>Figure 4: Small perturbation of the lake at rest solution computed with the GF-WENO5 WB scheme: $h-h_{eq}$ (red) and rescaled $b$ (blue) with $N_e = 150$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5: Small perturbation of the lake at rest solution computed with the WENO5 scheme: $h-h_{eq}$ (red) and rescaled $b$ (blue) with $N_e = 150$.</figcaption></figure>

mesh with 150 cells. The gravitational constant is considered to be $g = 9.8$ and the simulation is run until the final time $T = 1.5$.

Figure 4 shows the evolution over time of the perturbation added over the lake at rest solution computed with the GF-WENO5 WB scheme. To better present the results, the plots show the relative variable $h-h_{eq}$ where $h_{eq}$ represents the lake at rest solution without perturbation provided in Equation (45). We then compare the results in Figure 4 with those computed using the classical WENO5 approach, pointed out in Figure 5. It should be noticed that the classical approach fails at correctly reproducing the perturbation. Indeed, it generates a discretization error higher than the perturbation itself, eventually spoiling the final result. Contrary to that, the GF-WENO5 WB scheme correctly reproduces the perturbation which splits into two waves traveling at opposite directions and interacting with the bathymetry. The obtained result is fairly accurate, also considering the coarse mesh used for this study. To reach results similar to the GF one with $N_e = 150$, we need around $N_e = 800$ cells for the classical WENO5 method, see Figure 6.

## 5.3 Steady states with smooth bathymetry without friction ($n = 0$)

Here, we test the method for some moving equilibria steady state problems. We run the tests up to convergence towards steady state in different situations. In the subcritical and supercritical tests the bathymetry is smooth and equal to (46), as we want to assess the high order accuracy of the schemes. For the transcritical tests we use a modification of the bathymetry used in [27], to study a very similar discontinuous problem. Depending on the initial and boundary conditions set at the borders of the domain, the flow may be supercritical, subcritical or transcritical. The meshes used are defined by $N_e$ uniform cells, with different $N_e \in [25, 800]$. We consider the following three sets of final time $T$ and, initial and boundary conditions:
<figure><img src="image_6.png" /><figcaption>Figure 6: Small perturbation of the lake at rest solution computed with the WENO5 scheme: $h - h_{eq}$ (red) and rescaled $b$ (blue) with $N_e = 800$.</figcaption></figure>

• Supercritical flow

$$ 
\begin{array}{l}
T = 50, \\
h(x, 0) = 2 - b(x), \quad q(x, 0) \equiv 0, \\
h(0, t) = 2, \quad q(0, t) = 24,
\end{array} 
\qquad (50) $$

• Subcritical flow

$$ 
\begin{array}{l}
T = 200, \\
h(x, 0) = 2 - b(x), \quad q(x, 0) \equiv 0, \\
q(0, t) = 4.42, \quad h(25, t) = 2,
\end{array} 
\qquad (51) $$

• Transcritical flow

$$ b(x) = \begin{cases} 0.2 \exp \left( 1 - \frac{1}{1 - \left( \frac{|x-10|}{5} \right)^2} \right), & \text{if } |x-10| < 5, \\ 0, & \text{else,} \end{cases} \qquad (52) $$

$$ 
\begin{array}{l}
T = 200, \\
h(x, 0) = 0.33 - b(x), \quad q(x, 0) \equiv 0, \\
q(0, t) = 0.18, \quad h(25, t) = 0.33.
\end{array} 
\qquad (53) $$

The gravitational constant is set to $g = 9.812$ for all these tests. For these three cases, we compare the results obtained using a classical WENO finite volume scheme and the new approach based on flux globalization. As already mentioned, for supercritical and subcritical cases, the bathymetry (46) allows us to perform convergence tests for very high order methods, when also the flow is smooth. For the transcritical case with shock, only a qualitative analysis of the test is performed. Hence, when supercritical and subcritical flows are of interest, we can study the convergence properties of the new scheme by finding the exact solution given by the non-linear equations taken from [27]. Both the WB and non-WB versions of the scheme have been run to compare the influence of the formulation on the ability of preserving the balanced steady state solution. Finally, we also run the same test cases with the classical WENO3 and WENO5 schemes.

Convergence curves for supercritical and subcritical flows are depicted in Figures 7 and 8, respectively. All curves, both for WENO3 and WENO5, show the correct third and fifth order accuracy. However, it should be noticed that the GF formulation allows a much better prediction of the solution with errors dropping from 2 to 5 orders of magnitude. This is particularly striking for the well balanced global flux method exactly preserving the lake at rest state. To have a fairer comparison of the methods, one should compare
<figure><img src="image_7.png" /><figcaption>Figure 7: Supercritical flow: convergence tests with WENO3 and WENO5.</figcaption></figure>

the computational time with respect to the obtained error. Anyway, the methods have not been written in order to optimize the computational costs, so the comparison could be improved. We notice that for a fixed mesh we need more or less double the time of the classical FV scheme to run the GF-WB simulation. Let us consider fifth order methods and suppose that the errors of global flux are 10⁴ times smaller than the classical FV ones. Then, the computational time needed for classical FV schemes to obtain the same error of accuracy is 20 times larger than the GF-WB ones.

After the convergence analysis, we focus on a more qualitative analysis of the results computed with the WENO5 and GF-WENO5 methods. To do so we introduce a set of four representative variables that we are going to use to compare the classical WENO5 with the new GF-WENO5 WB. The three main variables that we used to assess the new method are $\eta$, $q$ and $K$ ($K$ is only defined for GF formulation). In particular, it should be noticed that, by construction, the methods based on flux globalization preserve moving equilibria, i.e., $q_x \equiv 0$, $K_x \equiv 0$. This means that the GF-WENO5 will approximate these quantities up to the order of the residual of the time derivative at the end of the simulation. In our case, we are able to preserve the constant $q$ and $K$ up to $\sim 10^{-9}$. The last variable we decided to study is the aforementioned $\Upsilon$, which corresponds to the smooth formulation of the more general global flux introduced herein. Although our approach is developed to preserve other equilibria, this variable allows to get more insights about the capability of the new algorithm since, for smooth flows, the analytical $\Upsilon$ should be constant at equilibrium. Solutions for supercritical and subcritical flows are shown in Figures 9 and 10. For both cases, it is clear that $q$ and $K$ are well-preserved, and $\Upsilon$ is much better predicted with respect to the one computed through the classical WENO5 method. The test case that stands out more among the three situations considered here is the transcritical flow with shock in Figure 11. The WENO5 method introduces spurious oscillations where the shock occurs that are then propagated in the rest of the computational domain making the overall solution spoiled. By using the new GF-WENO5 approach, oscillations are not present and the correct solution is recovered before and after the shock. It can be noticed that, when discontinuous flows are of interest, $\Upsilon$ is not globally constant but features a jump where the shock occurs, see Figure 11.

## 5.4 Perturbation of steady states without friction ($n = 0$)

In this section, we add a perturbation to the tests of Section 5.3 to compare the GF-WENO5 and the classical WENO5 methods for supercritical and subcritical flows. For the subcritical case we use the perturbation
<figure><img src="image_8.png" /><figcaption>Figure 8: Subcritical flow: convergence tests with WENO3 and WENO5.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9: Supercritical flow: relevant variables computed with GF-WENO5 (red continuous line) and WENO5 (black dashed line) schemes with $N_e = 100$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 10: Subcritical flow: relevant variables computed with GF-WENO5 (red continuous line) and WENO5 (black dashed line) schemes with $N_e = 100$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 11: Transcritical flow: relevant variables computed with GF-WENO5 (red continuous line), WENO5 (black dashed line) schemes and $b$ (blue dotted line) with $N_e = 100$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12: Small perturbation of the subcritical solution computed with the WENO5 scheme (black dashed) and GF WENO5 (red continuous): $h - h_{eq}$ with $N_e = 100$.</figcaption></figure>

(48) with $\alpha = 10^{-3}$. In Figure 12 we plot the solution for both methods at different timesteps with $N_e = 100$, while in Figure 13 we use $N_e = 800$. It can be noticed that both methods converge towards the exact solution, but with the GF-WENO5 method, even with a coarse mesh we obtain a very accurate approximation of the perturbation, while the WENO5 method performs poorly. A much more refined mesh is needed for the classical methods to nicely approximate this kind of solutions.

For the supercritical case we use the perturbation (48) with $\alpha = 10^{-4}$. In Figure 14 we plot the solution for both methods at different timesteps with $N_e = 100$, while in Figure 15 we use $N_e = 800$. Similar conclusions can be drawn from these simulations: the GF-WENO5 is accurate in representing the perturbation even for very coarse meshes, while the WENO5 needs more discretization cells.

<figure><img src="image_13.png" /><figcaption>Figure 13: Small perturbation of the subcritical solution computed with the WENO5 scheme (black dashed) and GF WENO5 (red continuous): $h - h_{eq}$ with $N_e = 800$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 14: Small perturbation of the supercritical solution computed with the WENO5 scheme (black dashed) and GF WENO5 (red continuous): $h - h_{eq}$ with $N_e = 100$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 15: Small perturbation of the supercritical solution computed with the WENO5 scheme (black dashed) and GF WENO5 (red continuous): $h - h_{eq}$ with $N_e = 800$.</figcaption></figure>

## 5.5 Steady states with discontinuous bathymetry without friction ($n = 0$)

Here, we change the bathymetry and test the subcritical and supercritical flows with a discontinuous step. The bathymetry is defined by

$$b = \begin{cases} 0.2 & \text{if } 8 < x < 12, \\ 0. & \text{else.} \end{cases} \qquad (54)$$

In Figure 16 we test the subcritical flow with the same initial conditions of Section 5.3 with $T = 500$ and $N_e = 100$. We can observe that the GF-WENO5 performs amazingly without producing any oscillations and obtaining errors of the order of $10^{-10}$, while the WENO5 scheme wildly oscillates after the discontinuities of the bathymetry.

For the supercritical case we have used the initial conditions of Section 5.3 till final time $T = 50$ with $N_e = 100$. In Figure 17 we observe that WENO5 performs better than before, still producing spurious oscillations and not exactly catching the outflow solution. On the other side, GF-WENO5 obtains constant global fluxes with a machine precision accuracy.

## 5.6 Steady states with friction ($n = 0.05$)

In this section we focus on the supercritical and subcritical flows studied in Section 5.3 including the friction term in the source term. As for the previous cases, also when friction with constant Manning coefficient $n$ is present, we can obtain moving equilibria. Again the quantity that are preserved at equilibrium are $q$ and $K$. That is why it is interesting to perform simulations similar to the previous ones comparing standard methods with GF ones. We consider the subcritical case defined in (51) and the supercritical case defined in (50) with the same bathymetry, on which we add the friction term with Manning coefficient $n = 0.05$. In the supercritical case, the friction term implies a slow down of the physical speed from left
<figure><img src="image_16.png" /><figcaption>Figure 16: Subcritical flow: relevant variables computed with GF-WENO5 (red continuous line), WENO5 (black dashed line) schemes and rescaled $b$ (blue dotted line) with $N_e = 100$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 17: Supercritical flow: relevant variables computed with GF-WENO5 (red continuous line) and WENO5 (black dashed line) schemes and rescaled $b$ (blue dotted line) with $N_e = 100$.</figcaption></figure>

to right and a consequent increasing of $\eta$ from left to right. In the subcritical case, conversely, we expect $h$ to decrease from left to right and the speed to increase. The variable $\Upsilon$ is not conserved and, for these tests, there is not another constant variable that can be easily computed analytically, hence, we do not plot it.

We display in Figures 18 and 19 both the solutions computed with the new GF-WENO5 WB scheme and the classical WENO5. It should be noticed that both schemes obtain valid and consistent results with the expected solution. The difference between the schemes is remarkable and the global flux variables clearly highlight it. The WENO5 case, without the global flux, is characterized by strong spurious oscillations around the area where the effect of the bathymetry is stronger. On the other side, the GF-WENO5 results are very precise and are able to preserve the global flux variables up to $\sim 10^{-9}$.

<figure><img src="image_18.png" /><figcaption>Figure 18: Supercritical flow with friction: relevant variables computed with GF-WENO5 (red continuous line) and WENO5 (black dashed line) schemes with $N_e = 100$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 19: Subcritical flow with friction: relevant variables computed with GF-WENO5 (red continuous line) and WENO5 (black dashed line) schemes with $N_e = 100$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 20: Small perturbation of the friction supercritical solution computed with the WENO5 scheme (black dashed) and GF WENO5 (red continuous): $h - h_{eq}$ with $N_e = 100$.</figcaption></figure>

## 5.7 Perturbation of steady states with friction ($n = 0.05$)

Also for these friction tests, we add the perturbation (48) to the final solution obtained at $T = 50$ for the supercritical case with $\alpha = 10^{-4}$ and at $T = 200$ with $\alpha = 10^{-3}$ for the subcritical case. Again, we use $N_e = 100$ cells of the domain. In Figure 20, we notice that the supercritical case with friction is the only case in which the classical WENO scheme is able to reach the steady state (residuals approach machine precision), hence the two solutions appear closer with respect to other cases. Probably, the friction plays a major role in this simulation by introducing some stabilization effects. Nevertheless, the classical WENO scheme exhibits a more oscillatory and unphysical behaviour in between the perturbation waves. In Figure 21, a much better prediction of the perturbation evolution and overall solution is obtained by the GF-WENO5 WB scheme. The classical WENO5 scheme, in this case, is not able to properly achieve the steady state, therefore producing a lot of numerical errors within the computational domain. This spoils the overall accuracy of the solution, that can be reobtained by refining the mesh.

# 6 Conclusion

In this paper, we presented a novel arbitrary high-order well-balanced finite-volume method based on flux globalization for the shallow water equations with source terms. The high order accuracy is obtained through a WENO reconstruction performed on the variables of interest, i.e. free surface level and global fluxes. The flux globalization allows to exactly preserve the constant fluxes in moving water equilibria. The new scheme has been also designed to preserve the “lake-at-rest” solution. This was possible by introducing a particular quadrature procedure for the source flux and considering a jump of the global fluxes at each interface. Several tests have been performed to assess the properties of the scheme. The preservation of the constant flux has
<figure><img src="image_21.png" /><figcaption>Figure 21: Small perturbation of the friction subcritical solution computed with the WENO5 scheme (black dashed) and GF WENO5 (red continuous): $h - h_{eq}$ with $N_e = 100$.</figcaption></figure>

been verified on academic moving water test cases (subcritical, transcritical and supercritical), while the well-balancedness with respect to the lake-at-rest equilibrium has been tested on the classical “lake-at-rest” and “perturbed lake-at-rest” cases. Once the space discretization is well implemented, the introduction of additional source terms is straightforward, e.g. Manning friction, as demonstrated in the simulation section. With the presented scheme, we were able to outperform classical methods in many situations, obtaining much more accurate solutions and useful properties also at the discrete level. Moreover, the high order accuracy of scheme allows to obtain precise solutions also when the equilibria are not reached. Several other extensions are being investigated. These include the DG-SEM approach of [46] as well as a continuous finite element formulation. The WENO work presented here could provide a setting for handling discontinuities in these works. Further improvements of the method might involve a posteriori limiters to deal with wet and dry areas. The benefits of the global flux idea in other settings (more complex equations, multi dimensional problems) is also under investigation.

## Acknowledgements

M. Ricchiuto is a member of the CARDAMOM team at INRIA and University of Bordeaux. During the development of this project, M. Ciallella was funded by an INRIA PhD fellowship and by a postdoctoral fellowship at ENSAM (I2M). D. Torlo has benefitted of INRIA postdoctoral fellowship, and by a postdoctoral fellowship in SISSA.

# A Positive reconstruction of water height at interfaces

*Proof.* From the cell averages a high order WENO reconstruction is performed on the fluxes to have $q_{i+1/2}^{L,R}$ and $K_{i+1/2}^{L,R}$. Equipped with $q_{i+1/2}^{L,R}$, $K_{i+1/2}^{L,R}$, $\mathcal{R}_{i+1/2}^{L,R}$, the point values $h_{i+1/2}^{L,R}$ can be obtained by solving the nonlinear equation coming from the definition of the global variable $K$ in (3):

$$K_{i+1/2}^{L} = \frac{\left(q_{i+1/2}^{L}\right)^{2}}{h_{i+1/2}^{L}} + \frac{g}{2} \left(h_{i+1/2}^{L}\right)^{2} + \mathcal{R}_{i+1/2}^{L}, \quad K_{i+1/2}^{R} = \frac{\left(q_{i+1/2}^{R}\right)^{2}}{h_{i+1/2}^{R}} + \frac{g}{2} \left(h_{i+1/2}^{R}\right)^{2} + \mathcal{R}_{i+1/2}^{R} \quad (55)$$

Let us solve the depressed cubic equation (55) for $h_{i+1/2}^{L}$ (the solution for $h_{i+1/2}^{R}$ is obtained with the same procedure). First of all, it can be noticed that (55) does not have any positive solution unless the determinant is greater than zero, meaning that

$$\left(q_{i+1/2}^{L}\right)^{4} < \frac{8\left(K_{i+1/2}^{L} - \mathcal{R}_{i+1/2}^{L}\right)^{3}}{27g}. \quad (56)$$
If (56) is not satisfied we reconstruct $\eta_{i+1/2}^L$ and then compute $h_{i+1/2}^L$ given that

$$h_{i+1/2}^L = \eta_{i+1/2}^L - b_{i+1/2}^L. \qquad (57)$$

If (56) is satisfied, then we have to deal with two possibilities. First, if $q_{i+1/2}^L = 0$, we obtain the unique positive solution

$$h_{i+1/2}^L = \sqrt{\frac{2(K_{i+1/2}^L - \mathcal{R}_{i+1/2}^L)}{g}},$$

while if $q_{i+1/2}^L \neq 0$, we solve Equation (55) for $h_{i+1/2}^L$ and obtain the following three solutions:

$$h_{i+1/2}^L = 2\sqrt{P} \cos\left(\frac{1}{3}[\Theta + 2\pi k]\right), \quad k = 0, 1, 2, \qquad (58)$$

where

$$P := \frac{2(K_{i+1/2}^L - \mathcal{R}_{i+1/2}^L)}{3g} \quad \text{and} \quad \Theta := \arccos\left(-\frac{(q_{i+1/2}^L)^2}{g P^{3/2}}\right). \qquad (59)$$

It can be shown that one of these roots is negative, whilst the other two roots, corresponding to the subcritical and supercritical cases, are positive. We choose the one closer to the corresponding value of $h_{i+1/2}^L$ given in (57).

## B Proof of proposition 3.3

To prove the result we will prove that, given WENO polynomials of $h$, $q$, and $b$, for the lake at rest state $[h(x) = \eta_0 - b(x) \quad q(x) = 0]^t$ we have $\tilde{G}_{i+1} = \tilde{G}_i \quad \forall i$. Clearly, the first component is the reconstruction of $q \equiv 0$, hence it is zero. Let us focus on $K$, the second component of the global flux. It is then enough to show that $\bar{K}_i$ is constant to have that $\bar{q}_i = 0 \quad \forall i$ at the following timesteps. Let us first compute the values of $K$ at each quadrature point. Simple computations show that

$$K_{i,q} = \mathcal{F}_{i,q} + \mathcal{R}_{i,q} = \mathcal{R}_{i-1/2}^R + g \frac{(\eta_0 - b_{i,q})^2}{2} + g \eta_0 (\tilde{b}_{i,q} - b_{i-1/2}^R) - g \left( \frac{(\tilde{b}_{i,q})^2}{2} - \frac{(b_{i-1/2}^R)^2}{2} \right) \quad (60)$$

$$= \mathcal{R}_{i-1/2}^R + g \frac{\eta_0^2}{2} - g \eta_0 b_{i-1/2}^R + g \frac{(b_{i-1/2}^R)^2}{2}. \qquad (61)$$

This shows that $K$ is constant across the quadrature points, and thus $\bar{K}_i = K_{i,q}$. We now need to show that this constant is the same $\forall i$. This is readily shown by substitution of the relevant quantities and some
simple algebra:

$$ \bar{K}_{i+1} - \bar{K}_i = \mathcal{R}_{i+1/2}^R - \mathcal{R}_{i-1/2}^R - g\eta_0 b_{i+1/2}^R + g\frac{(b_{i+1/2}^R)^2}{2} + g\eta_0 b_{i-1/2}^R - g\frac{(b_{i-1/2}^R)^2}{2} \quad (62a) $$

$$ = \mathcal{R}_{i+1/2}^L - \mathcal{R}_{i-1/2}^R + [\mathcal{R}_{i+1/2}] - g\eta_0 b_{i+1/2}^R + g\frac{(b_{i+1/2}^R)^2}{2} + g\eta_0 b_{i-1/2}^R - g\frac{(b_{i-1/2}^R)^2}{2} \quad (62b) $$

$$ = g\eta_0 \underbrace{(b_{i+1/2}^L - b_{i-1/2}^R) - g \left( \frac{(b_{i+1/2}^L)^2}{2} - \frac{(b_{i-1/2}^R)^2}{2} \right)}_{\mathcal{R}_{i+1/2}^L - \mathcal{R}_{i-1/2}^R} + [\mathcal{R}_{i+1/2}] \quad (62c) $$

$$ -g\eta_0 b_{i+1/2}^R + g\frac{(b_{i+1/2}^R)^2}{2} + g\eta_0 b_{i-1/2}^R - g\frac{(b_{i-1/2}^R)^2}{2} = \quad (62d) $$

$$ = g\eta_0 (b_{i+1/2}^L - b_{i-1/2}^R) - g \left( \frac{(b_{i+1/2}^L)^2}{2} - \frac{(b_{i-1/2}^R)^2}{2} \right) + [\mathcal{R}_{i+1/2}] = 0, \quad (62e) $$

recalling that

$$ [\mathcal{R}_{i+1/2}] = g\eta_0 (b_{i+1/2}^R - b_{i-1/2}^L) - g \left( \frac{(b_{i+1/2}^R)^2}{2} - \frac{(b_{i-1/2}^L)^2}{2} \right), \quad (63) $$

which achieves the proof. □

**Remark B.1** (Definition of the jump of $\mathcal{R}$). Clearly, the definition of $[\mathcal{R}_{i+1/2}]$ is obtained with the goal of having $\bar{K}_{i+1} - \bar{K}_i = 0$ for lake at rest equilibrium, and it can be easily derived from (62e).

## References

[1] R. Abgrall. High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices. *Journal of Scientific Computing*, 73(2-3):461–494, 2017.

[2] R. Abgrall, E. L. Mélédo, P. Öffner, and D. Torlo. Relaxation deferred correction methods and their applications to residual distribution schemes. *arXiv preprint arXiv:2106.05005*, 2021.

[3] R. Abgrall and M. Ricchiuto. *Hyperbolic Balance Laws: Residual Distribution, Local and Global Fluxes*, pages 177–222. Springer Nature Singapore, Singapore, 2022.

[4] L. Arpaia and M. Ricchiuto. r- adaptation for shallow water flows: conservation, well balancedness, efficiency. *Computers & Fluids*, 160:175–203, 2018.

[5] L. Arpaia and M. Ricchiuto. Well-balanced residual distribution for the ALE spherical shallow water equations on moving adaptive meshes. *J. Comput. Phys.*, 405(109173):1–32, 2020.

[6] L. Arpaia, M. Ricchiuto, A. G. Filippini, and R. Pedreros. An efficient covariant frame for the spherical shallow water equations: Well balanced DG approximation and application to tsunami and storm surge. *Ocean Modelling*, 169:101915, 2022.

[7] E. Audusse, F. Bouchut, M.-O. Bristeau, R. Klein, and B. t. Perthame. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. *SIAM Journal on Scientific Computing*, 25(6):2050–2065, 2004.

[8] D. S. Balsara and C.-W. Shu. Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. *Journal of Computational Physics*, 160(2):405–452, 2000.
[9] F. Behzadi and J. C. Newman. An exact source-term balancing scheme on the finite element solution of shallow water equations. *Computer Methods in Applied Mechanics and Engineering*, 359:112662, 2020.

[10] J. P. Berberich, P. Chandrashekar, and C. Klingenberg. High order well-balanced finite volume methods for multi-dimensional systems of hyperbolic balance laws. *Computers & Fluids*, 219:104858, 2021.

[11] A. Bermudez and M. Vazquez. Upwind methods for hyperbolic conservation laws with source terms. *Computers & Fluids*, 23(8):1049–1071, 1994.

[12] C. Berthon, S. Bulteau, F. Foucher, M. M'Baye, and V. Michel-Dansac. A very easy high-order well-balanced reconstruction for hyperbolic systems with source terms. *SIAM Journal on Scientific Computing*, 44(4):A2506–A2535, 2022.

[13] C. Berthon and C. Chalons. A fully well-balanced, positive and entropy-satisfying Godunov-type method for the shallow-water equations. *Mathematics of Computation*, 85(299):1281–1307, 2016.

[14] A. Bollermann, G. Chen, A. Kurganov, and S. Noelle. A well-balanced reconstruction of wet/dry fronts for the shallow water equations. *Journal of Scientific Computing*, 56(2):267–290, 2013.

[15] S. Bunya, S. Yoshimura, and J. J. Westerink. Improvements in mass conservation using alternative boundary implementations for a quasi-bubble finite element shallow water model. *International Journal for Numerical Methods in Fluids*, 51(11):1277–1296, 2006.

[16] V. Caselles, R. Donat, and G. Haro. Flux-gradient and source-term balancing for certain high resolution shock-capturing schemes. *Computers & Fluids*, 38(1):16–36, 2009.

[17] M. Castro, T. Morales de Luna, and C. Parés. Chapter 6 - Well-Balanced Schemes and Path-Conservative Numerical Methods. In R. Abgrall and C.-W. Shu, editors, *Handbook of Numerical Methods for Hyperbolic Problems*, volume 18 of *Handbook of Numerical Analysis*, pages 131–175. Elsevier, 2017.

[18] M. J. Castro and C. Parés. Well-balanced high-order finite volume methods for systems of balance laws. *Journal of Scientific Computing*, 82(2):48, 2020.

[19] 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. *Journal of Scientific Computing*, 80(1):538–554, 2019.

[20] Y. Cheng and A. Kurganov. Moving-water equilibria preserving central-upwind schemes for the shallow water equations. *Communications in Mathematical Sciences*, 14(6):1643–1663, 2016.

[21] 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. *Journal of Computational Physics*, 358:36–52, 2018.

[22] 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. *Journal of Scientific Computing*, 90(1):1–21, 2022.

[23] A. Christlieb, B. Ong, and J.-M. Qiu. Integral deferred correction methods constructed with high order Runge-Kutta integrators. *Mathematics of Computation*, 79(270):761–783, 2010.

[24] M. Ciallella, L. Micalizzi, P. Öffner, and D. Torlo. An arbitrary high order and positivity preserving method for the shallow water equations. *arXiv preprint arXiv:2110.13509*, 2021.

[25] J. W. Daniel, V. Pereyra, and L. L. Schumaker. Iterated deferred corrections for initial value problems. *Acta Cient. Venezolana*, 19:128–135, 1968.
[26] B. De Saint Venant. Theorie du mouvement non-permanent des eaux avec application aux crues des rivers et a l'introduntion des marees dans leur lit. *Academic de Sci. Comptes Redus*, 73(99):148-154, 1871.

[27] O. Delestre, C. Lucas, P.-A. Ksinant, F. Darboux, C. Laguerre, T.-N.-T. Vo, F. James, and S. Cordier. Swashes: a compilation of shallow water analytic solutions for hydraulic and environmental studies. *International Journal for Numerical Methods in Fluids*, 72(3):269-300, 2013.

[28] M. C. Díaz, J. A. López-García, and C. Parés. High order exactly well-balanced numerical methods for shallow water systems. *Journal of Computational Physics*, 246:242-264, 2013.

[29] R. Donat and A. Martínez-Gavara. Hybrid second order schemes for scalar balance laws. *Journal of Scientific Computing*, 48:52-69, 2011.

[30] A. Dutt, L. Greengard, and V. Rokhlin. Spectral Deferred Correction Methods for Ordinary Differential Equations. *BIT Numerical Mathematics*, 40(2):241-266, 2000.

[31] J. M. Gallardo, C. Parés, and M. Castro. On a well-balanced high-order finite volume scheme for shallow water equations with topography and dry areas. *Journal of Computational Physics*, 227(1):574-601, 2007.

[32] L. Gascón and J. Corberán. Construction of second-order TVD schemes for nonhomogeneous hyperbolic conservation laws. *Journal of Computational Physics*, 172(1):261-297, 2001.

[33] I. Gómez-Bueno, M. J. Castro, and C. Parés. High-order well-balanced methods for systems of balance laws: a control-based approach. *Applied Mathematics and Computation*, 394:125820, 2021.

[34] I. Gómez-Bueno, M. J. C. Díaz, C. Parés, and G. Russo. Collocation methods for high-order well-balanced methods for systems of balance laws. *Mathematics*, 9(15):1799, 2021.

[35] J. M. Greenberg and A. Y. Leroux. A well-balanced scheme for the numerical processing of source terms in hyperbolic equations. *SIAM Journal on Numerical Analysis*, 33(1):1-16, 1996.

[36] E. Hairer and G. Wanner. Solving ordinary differential equations. II, Vol. 14 of. Springer Series in Computational Mathematics (Springer Berlin Heidelberg, Berlin, Heidelberg, 1996), 10:978-3, 1996.

[37] M. Han Veiga, P. Öffner, and D. Torlo. DeC and ADER: Similarities, Differences and a Unified Framework. *Journal of Scientific Computing*, 87(1):1-35, 2021.

[38] G. Hauke. A symmetric formulation for computing transient shallow water flows. *Computer Methods in Applied Mechanics and Engineering*, 163(1-4):111-122, 1998.

[39] G. Hauke. A stabilized finite element method for the saint-venant equations with application to irrigation. *International Journal for Numerical Methods in Fluids*, 38(10):963-984, 2002.

[40] G.-S. Jiang and C.-W. Shu. Efficient implementation of weighted ENO schemes. *Journal of Computational Physics*, 126(1):202-228, 1996.

[41] K. Kashiyama, H. Ito, M. Behr, and T. Tezduyar. Three-step explicit finite element computation of shallow water flows on a massively parallel computer. *International Journal for Numerical Methods in Fluids*, 21(10):885-900, 1995.

[42] A. Kurganov. Finite-volume schemes for shallow-water equations. *Acta Numerica*, 27:289-351, 2018.

[43] A. Kurganov and D. Levy. Central-upwind schemes for the saint-venant system. *ESAIM: Mathematical Modelling and Numerical Analysis*, 36(3):397-425, 2002.
[44] Y. Liu, C.-W. Shu, and M. Zhang. Strong stability preserving property of the deferred correction time discretization. *Journal of Computational Mathematics*, pages 633–656, 2008.

[45] Y. Mantri and S. Noelle. Well-balanced discontinuous Galerkin scheme for 2×2 hyperbolic balance law. *Journal of Computational Physics*, 429:110011, 2021.

[46] Y. Mantri, P. Öffner, and M. Ricchiuto. Fully well balanced entropy controlled dgsem for shallow water flows: global flux quadrature and cell entropy correction. *arXiv preprint arXiv:2212.11931*, 2022.

[47] A. Meister and S. Ortleb. On unconditionally positive implicit time integration for the DG scheme applied to shallow water flows. *International Journal for Numerical Methods in Fluids*, 76(2):69–94, 2014.

[48] V. Michel-Dansac, C. Berthon, S. Clain, and F. Foucher. A well-balanced scheme for the shallow-water equations with topography or manning friction. *Journal of Computational Physics*, 335:115–154, 2017.

[49] V. Michel-Dansac, C. Berthon, S. Clain, and F. Foucher. A two-dimensional high-order well-balanced scheme for the shallow water equations with topography and manning friction. *Computers & Fluids*, 230:105152, 2021.

[50] M. L. Minion. Semi-implicit spectral deferred correction methods for ordinary differential equations. *Commun. Math. Sci.*, 1(3):471–500, 2009.

[51] S. Noelle, Y. Xing, and C.-W. Shu. High-order well-balanced finite volume WENO schemes for shallow water equation with moving water. *Journal of Computational Physics*, 226(1):29–58, 2007.

[52] M. Ricchiuto. On the C-property and generalized C-property of residual distribution for the shallow water equations. *Journal of Scientific Computing*, 48(1):304–318, 2011.

[53] M. Ricchiuto. An explicit residual based approach for shallow water flows. *Journal of Computational Physics*, 280:306–344, 2015.

[54] M. Ricchiuto, R. Abgrall, and H. Deconinck. Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes. *Journal of Computational Physics*, 222(1):287–331, 2007.

[55] M. Ricchiuto and A. Bollermann. Stabilized residual distribution for shallow water simulations. *Journal of Computational Physics*, 228(4):1071–1115, 2009.

[56] D. Serre. *Systems of conservation laws I - Hyperbolicity, Entropies, Shock waves*. Cambridge University Press, 1999.

[57] C.-W. Shu. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. In *Advanced numerical approximation of nonlinear hyperbolic equations*, pages 325–432. Springer, 1998.

[58] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes. *Journal of Computational Physics*, 77(2):439–471, 1988.

[59] T. Song, A. Main, G. Scovazzi, and M. Ricchiuto. The shifted boundary method for hyperbolic systems: Embedded domain computations of linear waves and shallow water flows. *Journal of Computational Physics*, 369:45–79, 2018.

[60] S. Takase, K. Kashiyama, S. Tanaka, and T. E. Tezduyar. Space-time supg finite element computation of shallow-water flows with moving shorelines. *Computational Mechanics*, 48(3):293, 2011.

[61] D. Torlo. *Hyperbolic problems: high order methods and model order reduction*. PhD thesis, PhD thesis, University Zurich, 2020.
[62] M. E. Vázquez-Cendón. Improved treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geometry. *Journal of Computational Physics*, 148(2):497-526, 1999.

[63] Y. Xing. Exactly well-balanced discontinuous galerkin methods for the shallow water equations with moving water equilibrium. *Journal of Computational Physics*, 257:536-553, 2014.

[64] Y. Xing and C.-W. Shu. High order finite difference WENO schemes with the exact conservation property for the shallow water equations. *Journal of Computational Physics*, 208(1):206-227, 2005.

[65] Y. Xing and C.-W. Shu. High order well-balanced finite volume weno schemes and discontinuous galerkin methods for a class of hyperbolic systems with source terms. *Journal of Computational Physics*, 214(2):567-598, 2006.

[66] Y. Xing, C.-W. Shu, and S. Noelle. On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations. *Journal of scientific computing*, 48(1):339-349, 2011.

[67] T. Yabe and Y. Ogata. Conservative semi-lagrangian cip technique for the shallow water equations. *Computational Mechanics*, 46(1):125-134, 2010.