# Approximate well-balanced WENO finite difference schemes using a global-flux quadrature method with multi-step ODE integrator weights

Maria Kazolea, Carlos Parés-Madronal, Mario Ricchiuto

► To cite this version:

Maria Kazolea, Carlos Parés-Madronal, Mario Ricchiuto. Approximate well-balanced WENO finite difference schemes using a global-flux quadrature method with multi-step ODE integrator weights. Computers and Fluids, 2025, 296, 10.1016/j.compfluid.2025.106646. hal-04890566

HAL Id: hal-04890566

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

Submitted on 16 Jan 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
# APPROXIMATE WELL-BALANCED WENO FINITE DIFFERENCE SCHEMES USING A GLOBAL-FLUX QUADRATURE METHOD WITH MULTI-STEP ODE INTEGRATOR WEIGHTS

M.Kazolea¹, C. Parés², M. Ricchiuto¹

¹INRIA, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, 200 Avenue de la Vieille Tour, 33405, Talence cedex, France

² University of Málaga, Spain

## ABSTRACT

In this work, high-order discrete well-balanced methods for one-dimensional hyperbolic systems of balance laws are proposed. We aim to construct a method whose discrete steady states correspond to solutions of arbitrary high-order ODE integrators. However, this property is embedded directly into the scheme, eliminating the need to apply the ODE integrator explicitly to solve the local Cauchy problem. To achieve this, we employ a WENO finite difference framework and apply WENO reconstruction to a global flux assembled nodewise as the sum of the physical flux and a source primitive. The novel idea is to compute the source primitive using high-order multi-step ODE methods applied on the finite difference grid. This approach provides a locally well-balanced splitting of the source integral, with weights derived from the ODE integrator. By construction, the discrete solutions of the proposed schemes align with those of the underlying ODE integrator. The proposed methods employ WENO flux reconstructions of varying orders, combined with multi-step ODE methods of up to order 8, achieving steady-state accuracy determined solely by the ODE method's consistency. Numerical experiments using scalar balance laws and shallow water equations confirm that the methods achieve optimal convergence for time-dependent solutions and significant error reduction for steady-state solutions.

**Keywords** approximate well-balanced, global-flux, multi-step methods, hyperbolic balance laws

## 1 Introduction

This paper is devoted to the solution of one-dimensional systems of balance laws of the form

$$U_t + F(U)_x = S(U)H_x + s(U, x), \quad (1)$$

where $U(x, t)$ takes value in $\Omega \subset \mathbb{R}^N$, $F: \Omega \to \mathbb{R}^N$ is the flux function; $S: \Omega \to \mathbb{R}^N$; and $H$ is a known function from $\mathbb{R} \to \mathbb{R}$ (possibly the identity function $H(x) = x$). The last term $s: \Omega \times \mathbb{R} \to \mathbb{R}^N$ accounts for extra point-wise sources. Let us assume initially that $H$ is continuous and $s$ is not singular. The system is supposed to be hyperbolic, i.e. the Jacobian $J(U)$ of the flux function is assumed to have $N$ different real eigenvalues.

Prototypes of scalar and systems of balance laws of the form (1) are:

* Burgers' equations with source term which corresponds to $N = 1$ and

$$F(U) = \frac{U^2}{2}, \quad S(U) = U^p, \quad s(U, x) \equiv 0. \quad (2)$$

* The shallow water system which corresponds to $N = 2$ and

$$U = \begin{bmatrix} h \\ q \end{bmatrix}, \quad F(U) = \begin{bmatrix} \frac{q}{h} \\ \frac{q^2}{2h^2} \end{bmatrix}, \quad S(U) = \begin{bmatrix} 0 \\ gh \end{bmatrix}. \quad (3)$$
The variable $x$ makes reference to the axis of the channel and $t$ is time; $q(x, t)$ and $h(x, t)$ represent the mass-flow and the thickness, respectively; $g$, the acceleration due to gravity; $H(x)$, the depth measured from a fixed level of reference; $q(x, t) = h(x, t)u(x, t)$, with $u$ the depth averaged horizontal velocity. If friction is also accounted for then we have

$$s(U, x) = \begin{bmatrix} 0 \\ -\kappa(U, x)q \end{bmatrix} \qquad (4)$$

with $\kappa(U, x) \ge 0$ a friction coefficient depending usually on hydrodynamic assumptions and on the bottom rugosity. Other effects can also be included, as e.g. wind forcing or variable atmospheric pressure, but will not be considered here.

We consider semi-discrete finite-difference numerical methods of the form

$$\frac{dU_i}{dt} + \frac{1}{\Delta x} (F_{i+1/2} - F_{i-1/2}) + S_i = 0, \qquad (5)$$

Here $U_i(t)$ represents an approximation of $U(x_i, t)$ where $x_0, \dots, x_N$ are the points of a uniform mesh of the computational interval $[a, b]$ and $\Delta x$ is the space step; $F_{i+1/2}$ is the numerical flux; and $S_i$ the numerical source term. We focus here on the development of high-order well-balanced methods for systems of the form (1), i.e. schemes that preserve in some sense all or some of the stationary solutions of the system, i.e. functions $U^*(x)$ that satisfy

$$F(U^*)_x = S(U^*)H_x + s(U^*, x).$$

For the prototypes introduced above, the stationary solutions are as follows:

* For Burgers' equations with source term (1)-(2) the stationary solutions satisfy the ODE

$$U_x = U^{p-1}H_x$$

whose solutions are given by

$$U e^{-H} = C \qquad (6)$$

if $p=2$ and

$$\frac{U^{2-p}}{2-p} - H = C \qquad (7)$$

otherwise. Here, $C$ represents an arbitrary constant.

* For the shallow water system (1)-(3), the stationary solutions satisfy in absence of friction

$$q = C_1, \quad \frac{1}{2} \frac{q^2}{h^2} + gh - gH = C_2, \qquad (8)$$

where $C_i, i = 1, 2$ are arbitrary constants. If $C_1 = 0$ water at rest or lake at rest stationary solutions

$$q = 0, \quad \eta = C_2, \qquad (9)$$

are obtained, where $\eta = h - H$ is the free surface elevation.

In presence of friction, the second in (8) must be replaced by the ODE

$$-C_1^2 h^{-2} h_x + gh(h_x - H_x) + \kappa(h, C_1, x)C_1 = 0 \qquad (10)$$

which can be integrated exactly in some particular cases (cf section §5.3.3).

Well-balanced methods are expected to preserve all the stationary solutions or at least a relevant family of them. They are also expected to correctly handle the waves produced by small perturbations of an equilibrium given by a steady state solution. Different definitions of well-balanced methods can be found in the literature according to the set of stationary solutions that is preserved and to the accuracy with which they are preserved. The following definitions will be adopted in this work:

**Definition 1.** The set of point-values $\{U_i\}_{i=0}^N$ is said to be a discrete stationary solution if it is an equilibrium for the ODE system (5), i.e. if

$$\frac{1}{\Delta x} (F_{i+1/2} - F_{i-1/2}) + S_i = 0.$$

**Definition 2.** The numerical method (5) is said to be exactly well-balanced for a stationary solution $U^*$ of (1) if the set of its point values $\{U^*(x_i)\}$ is a discrete stationary solution.
**Definition 3.** The numerical method (5) is said to be well-balanced for a stationary solution $U^*$ of (1) if, for all $\Delta x$, it is possible to find a discrete stationary solution $\{U_i^*\}$ that approximates $\{U^*(x_i)\}$ with an order of accuracy greater than that of the method (5).

**Definition 4.** The numerical method (5) is said to be fully exactly well-balanced (resp. fully well-balanced) if it is exactly well-balanced (resp. well-balanced) for every stationary solution.

The development of well-balanced numerical methods is a very active front of research. Finite volume methods based on approximate Riemann solvers that preserve exactly a family of stationary solutions were introduced by Roe [32], and then by Vázquez and Bermúdez in [5]: in the later numerical methods that exactly preserve water-at-rest solutions were introduced for the shallow-water equations. This well-balance property was called there C-property. A different strategy based on the use of standard Riemann problems for the homogeneous problem and the so-called hydrostatic reconstruction to discretize the source term was introduced in [2] and generalized in [11] in the framework of path-conservative numerical methods (see also [12]). Fully well-balanced numerical methods based on Riemann solvers were developed by Greenberg and Leroux [24], and then [20], and more recently [6] and subsequent works (e.g. [30] or [7] and references therein).

Strategies based on the computation of local equilibria at every cell and every time step have been followed in different works to obtain fully well-balanced methods. In the context of finite volume methods, such methods were built in [9] for semilinear balance laws by using the so-called well-balanced reconstruction operators: after having computed the local equilibrium, a standard reconstruction operator is applied to the fluctuations of the cell values in the stencil with respect to the local steady state; the reconstruction is then given by the sum of the equilibrium and the reconstructed fluctuations. Since then, this technique has been extended to general systems of balance laws: see [13] and the references therein. In particular, it was applied to the shallow-water equations in [10]. Different applications and extensions of this technique can be found in the literature: in the context of Euler equations with gravity, a similar technique was used in [26] to develop methods that are well-balanced for all the hydrostatic stationary solutions; a strategy based on the use of approximate local equilibria (computed by applying a RK-collocation solver to the ODE system satisfied by the stationary solutions) was introduced in [22] that allows one to use well-balanced reconstructions in problems in which the expression of the stationary solutions is not available neither explicitly nor implicitly; in [19] this technique was applied to obtain well-balanced arbitrary-lagrangian-eulerian finite volume schemes for the Euler equations, etc.

The strategy based on well-balanced reconstructions was extended to finite-difference methods in [13], where WENO reconstructions were applied to flux fluctuations. These methods are related to those based on the subtraction of the equations satisfied by a given stationary solution to be preserved (see [4] for instance).

An alternative strategy to build high-order well-balanced methods is based on the reconstruction in equilibrium variables: see [33], [16] and more recently [35] where a fifth-order fully well-balanced finite-difference method has been introduced.

Another way to construct well balanced approximations is to use the primitive of the evolution operator $G = F - \int (SH_x + s)$ as equilibrium variables. This method is equivalent to recasting (1) in a pseudo-conservation law form, since now the time evolution of the unknown can be formally expressed as the derivative of the global flux $G$. Methods based on the reconstruction of global fluxes have been introduced initially to devise TVD schemes for balance laws [21, 8]. Their use has seen a recent increased interest to produce fully well balanced methods [15, 17] which embed consistency with some integrated form of the steady state ODE. A clear characterization of this consistency has been provided in [28], in which, for finite element based schemes, the authors show that the global flux formulation allows to recover specific continuous collocation discretizations of the steady ODE. This has the advantage of giving a clear definition of the discrete steady state, with rigorous (super)consistency estimates. The formulation in the last reference shows that the iterative computation of the global flux on the mesh is not always necessary. For this reason the approach is rater referred to as global flux quadrature, as it boils down to a non-trivial quadrature of the source terms. Global flux techniques allow to embed the full evolution operator and have been shown to have many similarities to residual based techniques [1, 29]. An approach to generalize this notions to multiple dimensions has also been proposed in [3]. Global flux quadrature methods are not exactly well balanced in general, and fall in the class of methods of Definitions 3 and 4.

Coming back to [31], the family of well-balanced high-order methods for (1) introduced there for the case $s = 0$ write as follows:

$$ \frac{dU_i}{dt} + \frac{1}{\Delta x} (\hat{\mathcal{F}}_{i,i+1/2} - \hat{\mathcal{F}}_{i,i-1/2}) = 0, \qquad (11) $$

where the “numerical fluxes” $\hat{\mathcal{F}}_{i,i\pm 1/2}$ are computed as follows:
1. Look for the solution $U_i^*(x)$ of the Cauchy problem

$$ \begin{cases} F(U_i^*)_x = S(U_i^*)H_x + s(U^*, x), \\ U_i^*(x_i) = U_i. \end{cases} \qquad (12) $$

2. Define

$$ \mathcal{F}_j = F(U_j) - F(U_i^*(x_j)), \quad j = i - 1 - r, \dots, i + s \qquad (13) $$

3. Compute

$$ \hat{\mathcal{F}}_{i,i+1/2} = \mathcal{R}(\mathcal{F}_{i-r}, \dots, \mathcal{F}_{i+s}), \qquad (14) $$

$$ \hat{\mathcal{F}}_{i,i-1/2} = \mathcal{R}(\mathcal{F}_{i-1-r}, \dots, \mathcal{F}_{i-1+s}), \qquad (15) $$

where $\mathcal{R}$ represents a flux-reconstruction operator that gives high-order approximation of the value at the intercells of a function from the values at the points of the stencil $S_i = \{x_{i-r}, \dots, x_{i+s}\}$.

Observe that, (11) can be written in the form (5) with

$$ F_{i+1/2} = \mathcal{R}(F(U_{i-r}), \dots, F(U_{i+s})), \qquad (16) $$

and

$$ S_i = \frac{1}{\Delta x} \left( \mathcal{R}(F(U_i^*(x_{i-r})), \dots, F(U_i^*(x_{i+s}))) - \mathcal{R}(F(U_i^*(x_{i-1-r})), \dots, F(U_i^*(x_{i-1+s}))) \right), \qquad (17) $$

if the reconstruction operator is such that

$$ \mathcal{R}(\mathcal{F}_{i-r}, \dots, \mathcal{F}_{i+s}) = \mathcal{R}(F_{i-r}, \dots, F_{i+s}) - \mathcal{R}(F(U_i^*(x_{i-r})), \dots, F(U_i^*(x_{i+s}))). \qquad (18) $$

These methods are fully exactly well-balanced. Their main drawback is that a local stationary solution has to be computed at every point and every time step. This may be costly or even impossible when the analytic expression in explicit or implicit form of the stationary solutions is not available. The method is for example not applicable in general to flows with friction. In these cases, an alternative is to apply an RK-collocation method to the Cauchy problem (12), like in [22] to obtain approximations $U_{i,j}^* \approx U_i^*(x, j)$ and then define

$$ \mathcal{F}_j = F(U_j) - F(U_{i,j}^*), \quad j = i - 1 - r, \dots, i + s. $$

The numerical methods are then fully well-balanced but not exactly: what is preserved now is not the set of point values of stationary solutions but their approximations using the selected RK collocation method.

The goal of this paper is to obtain well-balanced high-order methods that do not require necessarily the computation of stationary solutions neither exact nor approximate at every mesh point in every time step. We do however require a rigorous characterization of the discrete steady state. To do this, we start from the global flux approach. The PDE system is thus written in the form

$$ U_t + F(U)_x - R_x = 0, $$

where

$$ R = \int_a^x (S(U)H_x + s(U, x)) dx $$

for some $a \in \mathbb{R}$; then, the flux reconstruction operator is applied to reconstruct $F(U) - R$.

In [17], a quadrature formula with points in the cell was considered to compute $R$ and WENO reconstructions were used to approximate the solution at the quadrature points. The key point here is that a quadrature formula

$$ \int_{x_j}^{x_{j+1}} S(U)H_x dx \approx \Delta x \sum_{m=0}^{s} \beta_m [S(U_{j+1-s+m})H_x(x_{j+1-s+m}) + s_{j+1-s+m}] $$

using the approximations of the solution at the left will be considered instead. Moreover, in the spirit of [28], we build some inherent consistency with well known high order ODE integrators. To this end, the coefficients $\beta_m$ of the formula will be those of a multistep ODE solver and it will be shown that the numerical methods are well-balanced in the sense that they preserve the approximations of a stationary solution obtained by applying the multistep method to (12) in an adequate form. This has the neat advantage that, if necessary, with one sweep of the ODE solver one can obtain the discrete steady state to use as well prepared initial data. The methods are extended to the approximations of singular sources, and applied to the scalar Burger equation, and to the shallow water equations with bathymetry, and friction.
The organization of the paper is as follows: in Section 2, the general form of the high-order finite-difference methods based on global-flux quadrature are introduced. Then, the discrete stationary solutions are studied and the well-balanced property is stated and shown. For the sake of completeness, the expression of the WENO reconstruction operators used here is given to conclude the section. Section 3 is devoted to enhancements of the methods: their extension to problems with singular source terms and their modification to satisfy the C-property, i.e. to exactly preserve water-at-rest solution for the shallow-water equations. The implementation details are given in Section 4 and the numerical results are shown in Section 5. There, the numerical methods are applied to the two prototypes equations and also to the shallow-water equations with friction to check their accuracy and well-balance properties. Finally conclusions are drawn.

# 2 Well-balanced scheme using global flux quadrature

## 2.1 Global-flux approach

High-order finite-difference methods for (1) of the form

$$ \frac{dU_i}{dt} + \frac{1}{\Delta x} (\hat{\mathcal{F}}_{i+1/2} - \hat{\mathcal{F}}_{i-1/2}) = 0 \qquad (19) $$

will be considered. The algorithm to compute the generalized numerical fluxes, based on the global-flux approach, is as follows:

1. Compute approximations

$$ R_j \approx \int_{x_0}^{x_j} (S(U)H_x + s(U, x)) \, dx, \quad j = 0, \dots, N, $$

as follows:

$$ R_0 = 0, \quad R_{j+1} = R_j + I_j(U), \quad j = 0, \dots, N-1, \qquad (20) $$

where

$$ \mathcal{I}_j(U) \approx \int_{x_j}^{x_{j+1}} (S(U)H_x + s(U, x)) \, dx. $$

2. Define

$$ \mathcal{F}_j = F(U_j) - R_j, \quad j = i-1-r, \dots, i+s. \qquad (21) $$

3. Use a high-order flux-reconstruction operator to compute

$$ \hat{\mathcal{F}}_{i+1/2} = \mathcal{R}(\mathcal{F}_{i-r}, \dots, \mathcal{F}_{i+s}), \qquad (22) $$

$$ \hat{\mathcal{F}}_{i-1/2} = \mathcal{R}(\mathcal{F}_{i-1-r}, \dots, \mathcal{F}_{i-1+s}). \qquad (23) $$

A quadrature formula that use $U_{j+1-s}, \dots, U_{j+1}$ to approximate the integral of the source term in the interval $[x_j, x_{j+1}]$ will be considered here

$$ \begin{aligned} \mathcal{I}_j(U) := & \Delta x \sum_{m=0}^{s} \beta_m S(U_{j+1-s+m}) H_x(x_{j+1-s+m}) \\ & + \Delta x \sum_{m=0}^{s} \beta_m s(U_{j+1-s+m}, x_{j+1-s+m}) \end{aligned} \qquad (24) $$

whose coefficients $\beta_m$ are obtained from some multi-step ODE integrator of choice, as e.g. those of Adams-Bashforth methods with $s$ steps (of order $s$) if $\beta_s = 0$, or Adams-Moulton with $s$ steps (of order $s+1$) if $\beta_s \neq 0$. In principle, these integrators can be chosen with accuracy greater than that of the spatial reconstruction operator. Other methods as e.g. multi-step and spectral multi-step collocation methods can be used, see e.g. [27, 34].

Observe that, if

$$ \mathcal{R}(\mathcal{F}_{i-r}, \dots, \mathcal{F}_{i+s}) = \mathcal{R}(F_{i-r}, \dots, F_{i+s}) - \mathcal{R}(R_{i-r}, \dots, R_{i+s}), \qquad (25) $$

(19) can be written in the form (5) with (16) and

$$ S_i = \frac{1}{\Delta x} \left( \mathcal{R}(R_{i-r}, \dots, R_{i+s}) - \mathcal{R}(R_{i-1-r}, \dots, R_{i-1+s}) \right). \qquad (26) $$

Note that (13) and (21) provide a similar modification to the scheme, the major difference being of course the way in which the subtracted flux is obtained. Another important difference between the two methods is the following conservation property that is not satisfied by method (19):
**Proposition 1.** Let us assume that (25) holds. If the reconstruction procedure is computed component-by-component and it is exact for the null function, then method (19)-(24) is conservative for the systems of conservation laws included in system (1).

*Proof.* Assume that the lth equation is a conservation law, i.e. the lth components of $S(U)$ and $s(U, x)$ vanish. Then, from (20) and (24) it is clear that the lth component of $R_j$ vanishes for all $j$. Finally, the assumptions on the reconstruction operator together with (26) implies that the lth component of the numerical source term also vanishes. □

In particular, this result implies that, under the assumption of the Proposition, mass-conservation will be satisfied by method (19) for the shallow-water equation, what is not in general the case for (11): see [31].

**Remark 1 (Evaluation of $H_x$).** The derivative of $H$ can be performed analytically or numerically with a sufficiently high accurate approximation. For the shallow water equations, considerations related to a particular steady state will lead to a specific approach which will be discussed in section §3.2.

## 2.2 Discrete stationary solutions

Let us assume that the flux reconstruction operator is such that

$$ \mathcal{R}(\mathcal{F}, \dots, \mathcal{F}) = \mathcal{F}, \quad \forall \mathcal{F}. \qquad (27) $$

Then, a sequence of point values $\{U_j^*\}$ satisfying

$$ F(U_{j+1}^*) - F(U_j^*) = \mathcal{I}_j(U^*), \quad j = 0, \dots, N-1 \qquad (28) $$

is a discrete stationary solution. In effect, if $\{U_j^*\}$ satisfies (28), then one has:

$$ \mathcal{F}_{j+1} = F(U_{j+1}^*) - R_{j+1} = F(U_j^*) - R_j = \mathcal{F}_j, \quad \forall j, $$

so that the numerical fluxes are constant and then

$$ \hat{\mathcal{F}}_{i,i+1/2} - \hat{\mathcal{F}}_{i,i-1/2} = \mathcal{F}_0 - \mathcal{F}_0 = 0, $$

and the right-hand side of (19) vanishes. Observe that (28) is the numerical counterpart of the equality

$$ F(U^*(x_{j+1})) - F(U^*(x_j)) = \int_{x_j}^{x_{j+1}} (S(U^*)H_x + s(U^*, x)) dx, \quad j = 0, \dots, N-1 $$

satisfied by the exact stationary solution of the system.

Discrete stationary solutions satisfying (28) can be computed as follows: given $U_0^*, \dots, U_{s-1}^*$ find $U_j^*$, $j = s, \dots, N$ such that

$$ \begin{aligned} F(U_j^*) = F(U_{j-1}^*) &+ \Delta x \sum_{m=0}^{s} \beta_m S(U_{j-sm}^*) H_x(x_{j-s+m}) \\ &+ \Delta x \sum_{m=0}^{s} \beta_m s(U_{j-sm}^*, x_{j-s+m}), \quad j = s, \dots, N. \end{aligned} \qquad (29) $$

Observe that a nonlinear system has to be solved to compute $U_j^*$ at every step and this happens even if the chosen multistep method is explicit, i.e. if $\beta_s = 0$. In this case, $U_j^*$ can be computed as follows:

1. Compute

$$ F_j^* = F(U_{j-1}^*) + \Delta x \sum_{m=0}^{s-1} \beta_m (S(U_{j-s+m}^*) H_x(x_{j-s+m}) + s(U_{j-s+m}^*, x_{j-s+m})). $$

2. Find $U_j^*$ such that:

$$ F(U_j^*) = F_j^*. $$

Observe that the second step requires the inversion of the flux function. We will assume here that the flux function

$$ U \mapsto F = F(U) $$

is one-to-one, what is true in a neighborhood of every state $U$ that is not sonic, i.e. if the eigenvalues of the Jacobian of the flux $DF(U)$ do not vanish.

The following result holds:
**Theorem 1** (Discrete stationary states). If the flux function is a one-to-one map, then the numerical method (19)-(24) is fully well-balanced with formal accuracy equal to the multi-step ODE solver with weights {$\beta_m$}$_{m=0,s}$.

*Proof.* Let $U^*$ be a stationary solution of (1). Then, $F^* \equiv F(U^*)$ satisfies the Cauchy problem

$$ egin{cases} F_x^* = S(U(F^*))H_x + s(U(F^*), x), \\ F^*(x_0) = F(U^*(x_0)), \end{cases} $$

where $U(F)$ represents the inverse of the flux function. If the selected multistep method is applied to this Cauchy problem starting from $F^*(x_j)$, $j = 0, \dots, s-1$ or from some high-order approximations of them obtained with a one-step method, a sequence of values {$F_i^*$} is obtained that satisfies for $j = s, \dots, N$

$$ F_j^* = F_{j-1}^* + \Delta x \sum_{m=0}^{s} \beta_m [S(U(F_{j-s+m}^*))H_x(x_{j-s+m}) + s(U(F_{j-s+m}^*), x_{j-s+m})]. $$

Then, the set of values

$$ U_j^* = U(F_j^*), \quad j = 0, \dots, N $$

satisfies (29) and, as a consequence, {$U_j^*$} is a discrete stationary solution that approximates $U^*(x_j)$, $j = 0, \dots, N$ with the order of accuracy of the multistep method. $\square$

**Remark 2** (Computation of $U(F)$). The main assumption of the theorem is that $F(U)$ should be a one-to-one mapping. It must be noted that, strictly speaking, for many systems this assumption may be violated in the neighborhood of sonic states. Indeed, in these points the Jacobian of the flux function is singular, and two (or more depending on the system) branches of the function $U(F)$ meet. This aspect is covered for the shallow water equations in many papers, see e.g. [15] section §2 for a discussion in the context of a global flux method. As shown in the reference, the inversion of $F(U)$ is well posed provided that we know on which branch $U$ should lay. We also remark that in practice, unlike the method (11)-(13), the knowledge of $U(F)$ is not strictly necessary here, unless the discrete stationary solutions are used as well-prepared initialization. The theoretical aspects related to sonic points are thus still unclear and under consideration. We will evaluate in practice the behaviour of the scheme in the results section. Finally, the fact that no non-linear system needs to be solved at every point, like e.g. in [22], has some advantages in terms of computational requirements.

Theorem 1 shows that the method proposed exactly preserves discrete stationary solutions obtained with an ODE solver, in a spirit similar to [22, 28]. As we will see in Section 5, the accuracy with which the methods approximate smooth exact stationary solutions can be changed at will to decrease the error almost to machine accuracy. However, some particular cases deserve special attention, either because their simplicity allows an exact preservation, or because they involve inherently discontinuous solutions. This aspect is discussed in sections §3.2 and §3.1.

## 2.3 WENO reconstructions

In the particular case of the WENO reconstruction of order $p = 2k + 1$, two flux reconstructions are computed using the values at the points $x_{i-k}, \dots, x_{i+k}$:

$$ \hat{\mathcal{F}}_{i+1/2}^{L} = \mathcal{R}^{L}(\mathcal{F}_{i-k}, \dots, \mathcal{F}_{i+k}), \qquad (30) $$

$$ \hat{\mathcal{F}}_{i-1/2}^{R} = \mathcal{R}^{R}(\mathcal{F}_{i-k}, \dots, \mathcal{F}_{i+k}). \qquad (31) $$

$\mathcal{R}^L$ and $\mathcal{R}^R$ represent the so-called left and right-biased reconstructions. The right-biased reconstructions can be computed with $\mathcal{R}^L$ by reflecting the arguments through the intercell at which the numerical flux is computed.

To compute the numerical fluxes, first a splitting

$$ \mathcal{F}_j = \mathcal{F}_j^+ + \mathcal{F}_j^-, \quad \forall j, $$

is considered and then they are defined as follows

$$ \hat{\mathcal{F}}_{i+1/2} = \hat{\mathcal{F}}_{i+1/2}^{+} + \hat{\mathcal{F}}_{i+1/2}^{-}, \qquad (32) $$

with

$$ \hat{\mathcal{F}}_{i+1/2}^{+} = \mathcal{R}^{L}(\mathcal{F}_{i-k}^{+}, \dots, \mathcal{F}_{i+k}^{+}), \qquad (33) $$

$$ \hat{\mathcal{F}}_{i+1/2}^{-} = \mathcal{R}^{R}(\mathcal{F}_{i-k+1}^{-}, \dots, \mathcal{F}_{i+k+1}^{-}). \qquad (34) $$
The general upwind biased splitting

$$ \mathcal{F}_{j}^{\pm} = \frac{1}{2} \left( \mathcal{F}_{j} \pm \frac{1}{\Delta x} \mathcal{T} J_{i+1/2} \mathcal{F}_{j} \right), $$

is considered here, where $\mathcal{T}$ is a time-scale matrix and $J_{i+1/2}$ is an approximation of the Jacobian of the flux function $F(U)$ at the intercell (a Roe matrix, for instance). $J_{i+1/2}$ is supposed to have $N$ real different eigenvalues. The choice

$$ \mathcal{T} = \Delta x |J_{i+1/2}|^{-1} $$

leads to the upwind splitting and

$$ \mathcal{T} = \frac{\Delta t}{2}. $$

to a Lax-Wendroff type splitting. The former is used in our implementation. Observe that this choice guarantees (27): if

$$ \mathcal{F}_j = \mathcal{F}, \quad \forall j, $$

then

$$ \hat{\mathcal{F}}_{i+1/2}^{+} = \frac{1}{2} \left( I + \frac{1}{\Delta x} \mathcal{T} J_{i+1/2} \mathcal{F} \right), \quad \hat{\mathcal{F}}_{i+1/2}^{-} = \frac{1}{2} \left( I - \frac{1}{\Delta x} \mathcal{T} J_{i+1/2} \mathcal{F} \right), $$

and then

$$ \hat{\mathcal{F}}_{i+1/2} = \left( \frac{1}{2} \left( I + \frac{1}{\Delta x} \mathcal{T} J_{i+1/2} \mathcal{F} \right) + \frac{1}{2} \left( I - \frac{1}{\Delta x} \mathcal{T} J_{i+1/2} \mathcal{F} \right) \right) \mathcal{F} = \mathcal{F}. $$

On the other hand, it is also clear that $\hat{\mathcal{F}}_{i+1/2} = 0$ if $\mathcal{F}_j = 0$ for all $j$. Moreover, since WENO reconstructions are linear combinations of the flux values, (18) or (25) are true once the coefficients of these linear combination (that depend non-linearly of the data through the smoothness indicators) have been computed. Therefore, the hypothesis of Proposition 1 are satisfied and the well-balanced property holds.

**Remark 3.** The Lax-Friedrichs flux-splitting:

$$ \mathcal{F}_{j}^{\pm} = \frac{1}{2} (\mathcal{F}_{j} \pm \alpha U), $$

where $\alpha$ is the local (WENO-LLF) or global (WENO-LF) maximum of the absolute value of the eigenvalues of $\{J(U_i)\}$, does not satisfy (27) in general. In effect, if

$$ \mathcal{F}_j = \mathcal{F}, \quad \forall j, $$

one would have

$$ \hat{\mathcal{F}}_{i, i+1/2} = \mathcal{F} - \frac{\alpha}{2} (\mathcal{R}^L(U_{i-k}, \dots, U_{i+k}) - \mathcal{R}^R(U_{i-k+1}, \dots, U_{i+k+1})) \neq \mathcal{F}. $$

# 3 Enhancements for singular sources and water at rest

## 3.1 Discontinuous data and singular source terms

Let us assume now the case in which $s = 0$, and the function $H$ has at least one isolated discontinuity. In this case, the definition of weak solutions (and, in particular, of stationary solutions) of (1) becomes more difficult: a solution $U$ is expected to be discontinuous at the discontinuities of $H$ and, in this case, the source term $S(U)H_x$ cannot be defined within the distributional framework. The source term becomes then a non-conservative product that can be defined in infinitely many different forms: see [18]. Like in [13], the following criterion is assumed here to define the weak solutions of (1): a pair $([U^{-}, H^{-}],[U^{+}, H^{+}])$ can be the limits of $U$ and $H$ of an admissible weak solution of (1) to the left and to the right of a discontinuity point of $H$, $x^{*}$, if and only if there exists a solution of the ODE system

$$ \frac{d}{dH} F(U) = S(U) \qquad (35) $$

satisfying

$$ U(H^{\pm}) = U^{\pm}. \qquad (36) $$

This criterion implies for instance that, in the case of Burgers' equation with source term, the relations (6) or (7) have to be satisfied by admissible jumps at the discontinuities of $H$. In the case of the shallow water system, relations (8) must be satisfied at the discontinuities of the bottom function.

Let us suppose for simplicity that $H$ has a unique discontinuity point that is placed at an interface $x_{l+1/2}$. In order to avoid the appearance of oscillations in the computation of discrete stationary solutions with the selected multistep
method, an adaptive technique is used in which the number of steps is adjusted so that $x_{\ell+1/2}$ always stays out of the stencil, except when $U_{\ell+1}^*$ is computed. The notation

$$ \beta_m^r, \quad m = 0, \dots, r, \quad r = 2, \dots, s, $$

will be used in what follows to represent the weights of the r-step method. Using this notation, discrete stationary solutions are computed now as follows: given $U_0^*, \dots, U_{s-1}^*$

* For $j = s, \dots, \ell$ compute $U_j^*$ by

$$ F(U_j^*) = F(U_{j-1}^*) + \Delta x \sum_{m=0}^{s} \beta_m^s S(U_{j-s+m}) H_x(x_{j-s+m}), \quad (37) $$

* For $j = \ell + 1$ compute $U_{\ell+1}^*$ so that $([U_{\ell}^*, H^{-}(x_{\ell})], [U_{\ell+1}^*, H^{+}(x_{\ell+1})])$ is an admissible jump $x_{\ell+1/2}$.

* For $j = \ell + 2, \dots, \ell + s$ compute $U_j^*$ by

$$ F(U_j^*) = F(U_{j-1}^*) + \Delta x \sum_{m=0}^{r} \beta_m^r S(U_{j-r+m}) H_x(x_{j-r+m}), \quad (38) $$

with $r = j - l - 1$.

* For $j = \ell + s + 1, \dots, N$ compute $U_j^*$ by (37).

Once the computation of the discrete stationary solutions has been set, let us introduce the procedure to approximate

$$ \int_{x_j}^{x_{j+1}} S(U) H_x(x) dx $$

so that the numerical method preserves them:

* if $j < \ell$, proceed as usual

$$ \mathcal{I}_j(U) = \Delta x \sum_{m=0}^{s} \beta_m^s S(U_{j+1-s+m}) H_x(x_{j+1-s+m}); \quad (39) $$

* if $j = \ell$ we set

$$ \mathcal{I}_{\ell}(U) = \tilde{S}_{\ell+1/2}[H]_{\ell+1/2} \quad (40) $$

where $\tilde{S}_{i+1/2}$ is a linearization of $S(U)$ such that, across an admissible jump at the discontinuity point of $H$, one has

$$ [F]_{i+1/2} = \tilde{S}_{i+1/2}[H]_{i+1/2}; \quad (41) $$

* if $j = \ell + r$ with $1 \le r \le s$,

$$ \mathcal{I}_j(U) = \Delta x \sum_{m=0}^{r} \beta_m^r S(U_{\ell+1+m}) H_x(x_{\ell+1+m}); \quad (42) $$

* $j > \ell + s$, proceed as usual (39).

Exact formulas for the linearization of $S(U)$ for the prototype equations are provided in appendix B.

An alternative to (40) for computing $\mathcal{I}_{\ell}(U)$ so that (28) is satisfied for $j = \ell$ is given by the so-called Generalized Hydrostatic Reconstruction (GHR): see [11]. To compute the integral in this case, first an intermediate value $H_{\ell+1/2}$ has to be chosen between $H(x_{\ell})$ and $H(x_{\ell+1})$; next two states $U_{\ell+1/2}^{\pm}$ have to be computed such that

($[U_{\ell}, H(x_{\ell})], [U_{\ell+1/2}^{-}, H_{\ell+1/2}]$) and ($[U_{\ell+1/2}^{+}, H_{\ell+1/2}], [U_{\ell+1}, H(x_{\ell+1})]$) are admissible jumps. Then,

$$ \mathcal{I}_{\ell}(U) = F(U_{\ell+1}) - F(U_{\ell+1/2}^{+}) + F(U_{\ell+1/2}^{-}) - F(U_{\ell}). $$

When $([U_{\ell}, H(x_{\ell})], [U_{\ell+1}, H(x_{\ell+1})])$ is an admissible jump then $U_{\ell+1/2}^{-} = U_{\ell+1/2}^{+}$ and (28) is satisfied. In practice, $H_{\ell+1/2}$ can be taken equal to $H(x_{\ell})$ or to $H(x_{\ell+1})$ so that only one intermediate state has to be computed: see [11] for details.
## 3.2 Exact preservation of water-at-rest solutions

As discussed at the end of section §2.2, the proposed method preserves exactly some discrete stationary solutions but not the exact ones. In particular, in the case of the shallow water equations, the numerical methods as described above will not preserve exactly the point values of water-at-rest stationary solutions. To be more precise, they will yield some discrete approximations which satisfy exactly $q = 0$, but verify $\eta = constant$ only up to the order of accuracy of the selected multistep method. In this sense, the C-property stated in [5] is not satisfied. Nevertheless, it is easy to modify the method to recover this property.

To do this, we first compute approximations $\tilde{H}_i \approx H(x_i)$ such that

$$ \tilde{H}_i = \tilde{H}_{i-1} + \Delta x \sum_{m=0}^{s} \beta_m H_x(x_{j-s+m}), $$

i.e. the bottom function is approximated by applying the selected quadrature formula to the expression

$$ H(x_i) = H(x_{i-1}) + \int_{x_{i-1}}^{x_i} H_x(x) dx. $$

We now consider the equality

$$ \int_{x_i}^{x_{i+1}} g h H_x dx = \int_{x_i}^{x_{i+1}} g \eta H_x dx + \frac{g}{2} \llbracket (H^2) \rrbracket_{i+1/2}. \quad (43) $$

Based on this equality and introducing the approximations of $H$, the second component of $\mathcal{I}_j(U)$ is computed as follows

$$ \mathcal{I}_{2,j}(U) = \Delta x \sum_{m=0}^{s} \beta_m g \eta_{j+1-s+m} H_x(x_{j+1-s+m}) + \frac{g}{2} \llbracket (\tilde{H})^2 \rrbracket_{i+1/2}. \quad (44) $$

Let us see that the sequence of point values of a water-at-rest stationary solution $U_j^* = [h_j^*, q_j^*]^T$, where

$$ q_j^* = 0, \quad h_j^* = \eta^* + \tilde{H}_j, \quad \eta^* = constant, $$

is a discrete stationary solution, i.e. that (28) is satisfied. This equality is trivially satisfied for the first components. Let us check the second one: taking into account that the quadrature formula is exact for $H_x^M$, one has

$$ \begin{aligned} \mathcal{I}_{2,j}(U) &= \Delta x \sum_{m=0}^{s} \beta_m g \eta_{j+1-s+m}^* H_x(x_{j+1-s+m}) + \frac{g}{2} \llbracket (\tilde{H})^2 \rrbracket_{i+1/2} \\ &= g \eta^* \llbracket \tilde{H} \rrbracket_{i+1/2} + \frac{g}{2} \llbracket (\tilde{H})^2 \rrbracket_{i+1/2} \\ &= g \eta^* \llbracket \tilde{H} \rrbracket_{i+1/2} + g \tilde{H}_{i+1/2} \llbracket \tilde{H} \rrbracket_{i+1/2} \\ &= g \eta^*_{i+1/2} \llbracket h \rrbracket_{i+1/2} + g \tilde{H}_{i+1/2} \llbracket h \rrbracket_{i+1/2} \\ &= g \bar{h}^*_{i+1/2} \llbracket h \rrbracket_{i+1/2} \\ &= \llbracket g(h^*)^2/2 \rrbracket_{i+1/2}, \end{aligned} $$

where $\bar{w}_{i+1/2}$ represents the arithmetic mean of $w_i$, $w_{i+1}$ for any variable $w$. Therefore, (28) is satisfied as we wanted to prove.

Another possibility is to proceed as in [28] and consider a polynomial interpolation of $H$ of degree $M = s$, denoted by $H^M$. Then, the second component of $\mathcal{I}_j(U)$ is computed as follows

$$ \mathcal{I}_{2,j}(U) = \Delta x \sum_{m=0}^{s} \beta_m g \eta_{j+1-s+m} H_x^M(x_{j+1-s+m}) + \frac{g}{2} \llbracket (H^M)^2 \rrbracket_{i+1/2}. \quad (45) $$

In this case, the sequence of point values of a water-at-rest stationary solution $U_j^* = [h_j^*, q_j^*]^T$, where

$$ q_j^* = 0, \quad h_j^* = \eta^* + H(x_j), \quad \eta^* = constant, $$

is a discrete stationary solution, i.e. that (28) is satisfied: the proof is similar, taking into account that the quadrature formula is exact for $H_x^M$.

Both approaches have different appeals. The first provides in practice lower errors for smooth bathymetries. The second manages better bathymetries with lower regularity. Since for discontinuous bathymetries we use the local modification discussed in the previous section, we have chosen to use the first method in the numerical applications.
# 4 Implementation details

In Section 5 some numerical tests will be shown in which the implementation of the methods introduced in Section 2 is as follows:

* WENO reconstruction operators based on the upwind splitting are used (see 2.3). Like in [31], Jiang-Shu smoothness indicators (see [25]) are used here.

* The coefficients of multi-step Adams-Bashforth method of order $q$ (ABq) or Adams-Moulton of order $q$ (ABq) will be used to compute the integrals of the source term. These coefficients up to order 8 are shown in appendix A.

* TVD-RK time stepping is applied to the semi-discrete in space numerical methods (19) (see [23]).

* CFL stability condition

$$ \Delta t = CFL \frac{\Delta x}{(\max_i\{|\lambda_{j,i+1/2}|, 1 \le l \le N\})} $$

is imposed with $CFL \in (0, 1)$, where $\lambda_{j,i+1/2}$ are the eigenvalues of the Jacobian of the flux function computed at some intermediate states. Unless otherwise stated $CFL = 0.45$.

* Boundary conditions are imposed through the use of ghost nodes.

Taking into account these ingredients, two families of methods are considered:

* WENOpGF-ABq: WENO reconstructions of order $p = 3, 5, 7$ with the Global Flux approach using Adams-Bashforth method of order $q = 4, 6, 8$ ($s = q$ steps) to compute the integrals;

* WENOpGF-AMq: WENO reconstructions of order $p = 3, 5, 7$ with the Global Flux approach using Adams-Moulton method of order $q = 4, 6, 8$ ($s = q - 1$ steps);

In some test cases, the numerical results will be compared with non-well balanced standard WENO implementations that will be called WENOp-nWB. In cetrain cases, the discrete steady-state solution is computed using a sweep of the Adams-Bashforth multi-step solver. This involves finding the roots of a degree three equation which can be done very easily if the super/sub- critical nature of the flow is known (see e.g. [14]). A similar implementation can also been performed for the Adams-Moulton method, which however requires fixed point iterations instead of the direct in the computation of the cubic roots.

Since in many cases the numerical results given by WENOpGF-ABq and WENOpGF-AMq are similar, only those corresponding to the latter will be shown, unless otherwise stated.

# 5 Numerical validation

## 5.1 Burgers' equation

We solve in this section the scalar Burgers' equation

$$ U_t + \left(\frac{U^2}{2}\right)_x = S(U)H_x, \qquad (46) $$

We will consider time dependent, steady smooth, and discontinuous solutions.

### 5.1.1 Manufactured traveling solution

As a first safety check, we consider the manufactured solution

$$ U(x, t) = e^{-(x-x_0-Ct)^2} $$

admitted when setting

$$ S(U) = U - C, \quad H(x, t) = e^{-(x-x_0-Ct)^2} $$

in the Burgers' model. We solve the problem with $C = 1$ and $x_0 = 5$ on the computational domain $[0, 15]$ m. We compute the error of the solution in $L^1$ norm at time $t = 2$ sec for the different WENO reconstruction and multi-step (Adams) methods. Tables 1, 2 and 3 show the errors and measured order of convergence of the different methods. As one might expect for a time dependent solution, the order of accuracy is given by the minimum between the one of the WENO reconstruction, and that of the ODE integrator.
<table><caption>Table 1: Burger's equation. Test 5.1.1: traveling manufactured solution. $L^1$ error convergence for WENO<sub>p</sub>GF-AM4 and WENO<sub>p</sub>GF-AB4, $p = 3, 5, 7$</caption><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AB4</th><th colspan="2">WENO5GF - AB4</th><th colspan="2">WENO7GF - AB4</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>60</td><td>1.393e-01</td><td>-</td><td>1.451e-01</td><td>-</td><td>1.160e-01</td><td>-</td></tr><tr><td>120</td><td>1.515e-02</td><td>3.2</td><td>1.250e-02</td><td>3.5</td><td>1.223e-02</td><td>3.7</td></tr><tr><td>240</td><td>1.293e-03</td><td>3.5</td><td>8.264e-04</td><td>3.9</td><td>8.166e-04</td><td>3.9</td></tr><tr><td>480</td><td>1.212e-04</td><td>3.4</td><td>5.200e-05</td><td>4.0</td><td>5.173e-05</td><td>4.0</td></tr><tr><td>960</td><td>1.306e-05</td><td>3.2</td><td>3.259e-06</td><td>4.0</td><td>3.251e-06</td><td>4.0</td></tr></tbody></table><table><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AM4</th><th colspan="2">WENO5GF - AM4</th><th colspan="2">WENO7GF - AM4</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>60</td><td>3.661e-02</td><td>-</td><td>3.659e-02</td><td>-</td><td>1.273e-01</td><td>-</td></tr><tr><td>120</td><td>5.220e-03</td><td>2.8</td><td>3.103e-03</td><td>3.5</td><td>9.659e-04</td><td>3.9</td></tr><tr><td>240</td><td>7.020e-04</td><td>2.9</td><td>2.531e-04</td><td>3.6</td><td>6.254e-05</td><td>4.0</td></tr><tr><td>480</td><td>9.016e-05</td><td>2.9</td><td>1.816e-05</td><td>3.8</td><td>2.465e-07</td><td>4.0</td></tr><tr><td>960</td><td>1.140e-06</td><td>2.9</td><td>1.195e-06</td><td>3.9</td><td>1.541e-08</td><td>4.0</td></tr></tbody></table>

<table><caption>Table 2: Burger's equation. Test 5.1.1: traveling manufactured solution. $L^1$ error convergence for WENO<sub>p</sub>GF-AM6 and WENO<sub>p</sub>GF-AB6, $p = 3, 5, 7$</caption><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AB6</th><th colspan="2">WENO5GF - AB6</th><th colspan="2">WENO7GF - AB6</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>60</td><td>1.377e-01</td><td>-</td><td>1.204e-01</td><td>-</td><td>1.552e-01</td><td>-</td></tr><tr><td>120</td><td>5.331e-03</td><td>4.6</td><td>4.081e-03</td><td>4.8</td><td>4.306e-03</td><td>5.1</td></tr><tr><td>240</td><td>6.934e-04</td><td>2.9</td><td>6.571e-05</td><td>6.0</td><td>7.551e-05</td><td>5.8</td></tr><tr><td>480</td><td>9.127e-05</td><td>2.9</td><td>9.642e-07</td><td>6.0</td><td>1.217e-06</td><td>5.9</td></tr><tr><td>960</td><td>1.149e-05</td><td>2.9</td><td>1.905e-08</td><td>5.6</td><td>1.905e-08</td><td>6.0</td></tr></tbody></table><table><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AM6</th><th colspan="2">WENO5GF - AM6</th><th colspan="2">WENO7GF - AM6</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>60</td><td>4.419e-02</td><td>-</td><td>1.524e-02</td><td>-</td><td>1.090e-02</td><td>-</td></tr><tr><td>120</td><td>5.840e-03</td><td>2.9</td><td>7.391e-04</td><td>4.3</td><td>2.156e-04</td><td>5.9</td></tr><tr><td>240</td><td>7.320e-04</td><td>2.9</td><td>2.278e-05</td><td>5.0</td><td>3.622e-06</td><td>5.9</td></tr><tr><td>480</td><td>9.193e-05</td><td>2.9</td><td>6.921e-07</td><td>5.0</td><td>5.955e-08</td><td>5.7</td></tr><tr><td>960</td><td>1.150e-05</td><td>2.9</td><td>2.159e-08</td><td>5.0</td><td>7.418e-11</td><td>4.0</td></tr></tbody></table>

## 5.1.2 Approximation of a smooth steady state

We now consider the case

$$S(U) = U^2, \quad H(x) = x,$$

together with its stationary solution

$$U^*(x) = e^x \qquad (47)$$

in the interval $[-1, 1]$. For this case, we expect the ODE-based GF formulation to bring considerable accuracy enhancements. To check this we use the exact point values of the stationary solution as initial data, and then run until a discrete stationary solution is reached. The error in $L^1$-norm of the discrete steady state with respect to the analytical one is then measured. Tables 5-7 show the errors and empirical order of convergence of the different methods. It can be seen that the order of accuracy as well as the error values for the GF case are uniquely determined by the choice of the multi-step method. This allows for error reductions of several orders of magnitude on a given mesh, just by modifying the weights of the quadrature formula.

## 5.1.3 Highly oscillatory solution

As in [31] we consider now (46) with an oscillatory function $H$:

$$S(U) = U^2, \quad H(x) = x + 0.1 \sin(100x), \qquad (48)$$

together with the stationary solution

$$U^*(x) = e^{H(x)}$$
<table><caption>Table 3: Burger's equation. Test 5.1.1: traveling manufactured solution. $L^1$ error convergence for WENO<sub>p</sub>GF-AM8 and WENO<sub>p</sub>GF-AB8, $p = 3, 5, 7$</caption><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AB8</th><th colspan="2">WENO5GF - AB8</th><th colspan="2">WENO7GF - AB8</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>60</td><td>1.2846751</td><td>-</td><td>4.0572682</td><td>-</td><td>5.2918144</td><td>-</td></tr><tr><td>120</td><td>5.750e-03</td><td>7.8</td><td>1.8716079</td><td>1.1</td><td>1.2131538</td><td>2.1</td></tr><tr><td>240</td><td>7.331e-04</td><td>2.9</td><td>1.638e-02</td><td>6.8</td><td>4.032e-02</td><td>4.9</td></tr><tr><td>480</td><td>9.192e-05</td><td>2.9</td><td>9.845e-05</td><td>7.3</td><td>4.848e-06</td><td>13.0</td></tr><tr><td>960</td><td>1.150e-05</td><td>2.9</td><td>5.762e-07</td><td>7.7</td><td>5.588e-07</td><td>3.1</td></tr></tbody><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AM8</th><th colspan="2">WENO5GF - AM8</th><th colspan="2">WENO7GF - AM8</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>60</td><td>4.419e-02</td><td>-</td><td>1.524e-02</td><td>-</td><td>1.090e-02</td><td>-</td></tr><tr><td>120</td><td>5.840e-03</td><td>2.9</td><td>7.391e-04</td><td>4.3</td><td>2.156e-04</td><td>5.9</td></tr><tr><td>240</td><td>7.320e-04</td><td>2.9</td><td>2.278e-05</td><td>5.0</td><td>3.622e-06</td><td>5.9</td></tr><tr><td>480</td><td>9.193e-05</td><td>2.9</td><td>6.921e-07</td><td>5.0</td><td>5.955e-08</td><td>5.7</td></tr><tr><td>960</td><td>1.150e-05</td><td>2.9</td><td>2.159e-08</td><td>5.0</td><td>7.418e-11</td><td>4.0</td></tr></tbody></table>

<table><caption>Table 4: Burgers' equation. Test 5.1.1: traveling manufactured solution. $L^1$ error and convergence rates for WENO<sub>p</sub>-nWB, $p = 3, 5, 7$</caption><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3-nWB</th><th colspan="2">WENO5-nWB</th><th colspan="2">WENO7-nWB</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>60</td><td>5.734e-02</td><td>-</td><td>2.335e-02</td><td>-</td><td>8.761e-03</td><td>-</td></tr><tr><td>120</td><td>8.648e-03</td><td>2.7</td><td>1.634e-03</td><td>3.8</td><td>1.403e-04</td><td>6.0</td></tr><tr><td>240</td><td>1.140e-03</td><td>2.9</td><td>6.376e-05</td><td>4.7</td><td>1.678e-06</td><td>6.3</td></tr><tr><td>480</td><td>1.447e-04</td><td>2.9</td><td>2.062e-06</td><td>4.9</td><td>2.232e-08</td><td>6.2</td></tr><tr><td>960</td><td>1.815e-05</td><td>3.0</td><td>6.546e-08</td><td>5.0</td><td>3.262e-08</td><td>6.0</td></tr></tbody></table>

<table><caption>Table 5: Burger's equation. Test 5.1.2: smooth steady state. Errors in $L^1$ norm and convergence rate for WENO<sub>p</sub>GF-AM4 and WENO<sub>p</sub>GF-AB4, $p = 3, 5, 7$</caption><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AB4</th><th colspan="2">WENO5GF - AB4</th><th colspan="2">WENO7GF - AB4</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>1.380e-03</td><td>-</td><td>1.308e-03</td><td>-</td><td>1.380e-03</td><td>-</td></tr><tr><td>40</td><td>9.628e-05</td><td>3.8</td><td>9.388e-05</td><td>3.8</td><td>9.628e-05</td><td>3.8</td></tr><tr><td>80</td><td>6.359e-06</td><td>3.9</td><td>6.282e-06</td><td>3.9</td><td>6.359e-06</td><td>3.9</td></tr><tr><td>160</td><td>4.086e-07</td><td>4.0</td><td>4.086e-07</td><td>4.0</td><td>4.086e-07</td><td>4.0</td></tr><tr><td>320</td><td>2.590e-08</td><td>4.0</td><td>2.589e-08</td><td>4.0</td><td>2.900e-08</td><td>4.0</td></tr></tbody><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AM4</th><th colspan="2">WENO5GF - AM4</th><th colspan="2">WENO7GF - AM4</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>1.112e-04</td><td>-</td><td>1.115e-04</td><td>-</td><td>1.115e-04</td><td>-</td></tr><tr><td>40</td><td>7.865e-06</td><td>3.9</td><td>7.670e-06</td><td>3.9</td><td>7.669e-06</td><td>3.9</td></tr><tr><td>80</td><td>5.000e-07</td><td>3.9</td><td>5.000e-07</td><td>3.9</td><td>4.939e-07</td><td>3.9</td></tr><tr><td>160</td><td>3.152e-08</td><td>4.0</td><td>3.152e-08</td><td>4.0</td><td>3.155e-08</td><td>4.0</td></tr><tr><td>320</td><td>1.979e-09</td><td>4.0</td><td>1.979e-09</td><td>4.0</td><td>1.978e-09</td><td>4.0</td></tr></tbody></table>

in the interval $[-1, 1]$ (see Figure 1). We take as initial condition the stationary solution and first a 100-cell mesh is considered, so that the period of the oscillations is close to $\Delta x$. In Figure 2 (top-left) the numerical solutions obtained at time $t = 1$ sec with WENO3GF-AMq, $q = 4, 6, 8$ and WENO<sub>p</sub>-nWB, $p = 3, 5$ are compared with the exact stationary solution. As expected the well-balanced methods behave better than the non-WB ones even when their order of accuracy is lower. Nevertheless, the choice $q = 8$ is surprisingly worse than $q = 4$ and $q = 6$. This is a consequence of the large stencil introduced by the eighth order method, which is not suited to approximate the highly oscillatory solution considered. To verify this we consider a slight refinement of the mesh. In top right picture of Figure 2 we thus report the results obtained using 150 points. It can be seen that all methods give improved results. In this case WENO3GF-AM8 is the best one as shown by the zoom in the bottom picture in Figure 2.
<figure><img src="image_2.png" /><figcaption>Figure 1: Burgers' equation. Test 5.1.3: highly oscillatory solution. Graph of the function $H$</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Burgers' equation. Test 5.1.3: highly oscillatory solution. Top: discrete steady states with different schemes for $N = 100$ (left) and $N = 150$ (right). Bottom: zoom of the last oscillations for $N = 150$.</figcaption></figure>
<table><caption>Table 6: Burgers' equation. Test 5.1.2: smooth steady state. Errors in $L^1$ norm and convergence rate for WENO<sub>p</sub>GF-AM6 and WENO<sub>p</sub>GF-AB6, $p = 3, 5, 7$</caption><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AB6</th><th colspan="2">WENO5GF - AB6</th><th colspan="2">WENO7GF - AB6</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>3.917e-05</td><td>-</td><td>3.898e-05</td><td>-</td><td>3.889e-05</td><td>-</td></tr><tr><td>40</td><td>7.909e-07</td><td>5.7</td><td>7.862e-07</td><td>5.7</td><td>7.704e-07</td><td>5.7</td></tr><tr><td>80</td><td>1.370e-08</td><td>5.9</td><td>1.369e-08</td><td>5.9</td><td>1.371e-08</td><td>5.8</td></tr><tr><td>160</td><td>2.256e-10</td><td>5.9</td><td>2.256e-10</td><td>5.9</td><td>2.255e-10</td><td>5.9</td></tr><tr><td>320</td><td>3.737e-12</td><td>5.9</td><td>3.774e-12</td><td>5.9</td><td>3.762e-12</td><td>5.9</td></tr></tbody></table><table><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AM6</th><th colspan="2">WENO5GF - AM6</th><th colspan="2">WENO7GF - AM6</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>2.202e-06</td><td>-</td><td>2.124e-06</td><td>-</td><td>2.088e-06</td><td>-</td></tr><tr><td>40</td><td>3.884e-08</td><td>5.8</td><td>3.879e-08</td><td>5.8</td><td>3.890e-08</td><td>5.8</td></tr><tr><td>80</td><td>6.459e-10</td><td>5.9</td><td>6.455e-10</td><td>5.9</td><td>6.453e-10</td><td>5.9</td></tr><tr><td>160</td><td>1.035e-11</td><td>6.0</td><td>1.035e-11</td><td>6.0</td><td>1.035e-11</td><td>6.0</td></tr><tr><td>320</td><td>4.159e-16</td><td>-</td><td>2.331e-15</td><td>-</td><td>1.993e-15</td><td>-</td></tr></tbody></table>

<table><caption>Table 7: Burgers' equation. Test 5.1.2: smooth steady state. Errors in $L^1$ norm and convergence rate for WENO<sub>p</sub>GF-AM8 and WENO<sub>p</sub>GF-AB8, $p = 3, 5, 7$</caption><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AB8</th><th colspan="2">WENO5GF - AB8</th><th colspan="2">WENO7GF - AB8</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>1.267e-06</td><td>-</td><td>2.071e-05</td><td>-</td><td>1.694e-05</td><td>-</td></tr><tr><td>40</td><td>6.702e-09</td><td>7.6</td><td>5.353e-04</td><td>-</td><td>5.875e-05</td><td>-</td></tr><tr><td>80</td><td>3.056e-11</td><td>7.8</td><td>1.338e-06</td><td>8.6</td><td>4.962e-06</td><td>3.6</td></tr><tr><td>160</td><td>8.677e-14</td><td>8.5</td><td>1.160e-13</td><td>-</td><td>3.425e-08</td><td>7.2</td></tr><tr><td>320</td><td>-</td><td>-</td><td>-</td><td>-</td><td>6.657e-14</td><td>-</td></tr></tbody></table><table><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3GF - AM8</th><th colspan="2">WENO5GF - AM8</th><th colspan="2">WENO7GF - AM8</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>4.792e-08</td><td>-</td><td>4.765e-08</td><td>-</td><td>4.718e-08</td><td>-</td></tr><tr><td>40</td><td>2.317e-10</td><td>7.7</td><td>2.314e-10</td><td>7.7</td><td>2.312e-10</td><td>7.7</td></tr><tr><td>80</td><td>9.666e-13</td><td>7.9</td><td>9.688e-13</td><td>7.9</td><td>9.828e-13</td><td>7.9</td></tr><tr><td>160</td><td>1.346e-16</td><td>-</td><td>1.429e-16</td><td>-</td><td>1.540e-16</td><td>-</td></tr><tr><td>320</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

<table><caption>Table 8: Burgers' equation. Test 5.1.2: smooth steady state. Errors in $L^1$ norm and convergence rates for WENO<sub>p</sub>-nWB, with $p = 3, 5, 7$</caption><thead><tr><th rowspan="2">N</th><th colspan="2">WENO3-nWB</th><th colspan="2">WENO5-nWB</th><th colspan="2">WENO7-nWB</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>9.531e-03</td><td>-</td><td>9.387e-05</td><td>-</td><td>5.073e-07</td><td>-</td></tr><tr><td>40</td><td>1.166e-03</td><td>3.0</td><td>3.468e-06</td><td>4.8</td><td>5.633e-09</td><td>6.5</td></tr><tr><td>80</td><td>1.431e-04</td><td>3.0</td><td>1.112e-07</td><td>5.0</td><td>4.571e-11</td><td>6.9</td></tr><tr><td>160</td><td>1.767e-05</td><td>3.0</td><td>3.503e-09</td><td>5.0</td><td>1.253e-13</td><td>8.5</td></tr><tr><td>320</td><td>2.193e-06</td><td>3.0</td><td>1.094e-10</td><td>5.0</td><td>-</td><td>-</td></tr></tbody></table>

In order to see how the different numerical methods handle perturbations in the next experiment we take initial data

$$U_j = U_j^* + 0.2\chi_{[-0.7, -0.5]}(x_j), \quad j = 1, \dots, 100,$$

where {$U_j^*$} is the discrete stationary solution reached by each method in the previous experiment with 100 points, and $\chi_{[a, b]}$ represents the characteristic function of the interval $[a, b]:$ see Figure 3 (left). The differences between the numerical solutions at time $t = 0.7$ sec, and the underlying discrete stationary solutions are compared with a reference solution on the right in Figure 3. It can be seen that on a coarse mesh the non-WB method exaggerate considerably the evolution of the perturbation, despite of the fact that the underlying steady state is its own discrete stationary state.
<figure><img src="image_4.png" /><figcaption>Figure 3: Burgers' equation. Test 5.1.3: highly oscillatory solution. Perturbed steady state (N = 150). Left: initial solution. Right: evolution of the perturbation of different schemes.</figcaption></figure>

### 5.1.4 Discontinuous data

Let us consider now (46) with a piecewise continuous function $H$:

$$H(x) = \begin{cases} 0.1x & \text{if } x \le 0; \\ 0.5 + x & \text{if } 0 < x \le 0.5; \\ 0.9 + x & \text{otherwise.} \end{cases} \qquad (49)$$

We consider again the interval $[-1, 1]$ and we take as initial condition the stationary solution

$$u(x) = \begin{cases} e^{0.1x} & \text{if } x \le 0; \\ e^{0.5+x} & \text{if } 0 < x \le 0.5; \\ e^{0.9+x} & \text{otherwise.} \end{cases} \qquad (50)$$

We deliberately chose two discontinuities on $H$. The first one, placed at $x = 0$ is aligned with a mesh point while the second one, placed at $x = 0.5$, is positioned at a cell face. Using the stationary solution as the initial condition over the interval $[-1, 1]$ on a mesh of $N = 100$ points, we examine the impact on discontinuity preservation. As illustrated in Figure 4 WENO3GF-AMq, $q = 4, 6, 8$ schemes successfully preserve both discontinuities irrespective of their location while WENO3-nWB smooths out the second discontinuity.

We now focus on the perturbation of the discontinuous solution. For clarity, we simplify the set up by retaining only a single discontinuity in $H$ so that function $H$ reads as:

$$H(x) = \begin{cases} 0.1x & \text{if } x \le 0; \\ 0.9 + x & \text{otherwise;} \end{cases} \qquad (51)$$

We consider again the interval $[-1, 1]$ and we take as initial condition the stationary solution with a small perturbation:

$$\tilde{u}_0(x) = u(x) + 0.3e^{-200(x+1/2)^2},$$

with $u(x) = e^H$, see also Figure 5. The numerical solutions at times $t = 0.3$ and $0.5$ obtained with WENO3GF-AMq, $q = 4, 6$ using a mesh with $N = 100$ and $300$ points are compared with the numerical solutions of the WENO3-nWB using the finer mesh, and with a reference obtained with the WENO5-nWB scheme with upwind fluxes and $5000$ points. As it can be seen, the non well balanced method provides a wrong prediction of the interaction with the discontinuity, while the well balanced approach proposed here with the correction to account for singularities allows to obtain an accurate result even on a very coarse mesh.

## 5.2 Shallow water equations

We now consider the shallow water system (1)-(3) and verify the correct implementation and convergence of our methods using classical benchmarks that involve static and moving equilibria.
<figure><img src="image_5.png" /><figcaption>Figure 4: Burgers’ equation. Test 5.1.3: discontinuous data. Preservation of two discontinuities at $t = 0.2$ sec on a mesh of $N = 100$ points</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Burgers’ equation. Test 5.1.3: discontinuous data. Perturbation of an isolated discontinuity: initial state.</figcaption></figure>

### 5.2.1 Approximation of the lake at rest solution over a discontinuous bottom

First, we focus on the lake at rest solution characterized by the initial data:

$$h(x, 0) = 2 - H(x), \quad q(x, 0) = 0 \tag{52}$$

with the bottom depth given by:

$$H(x) = \begin{cases} 0.05 \sin(x - 12.5) \exp(1 - (x - 12.5)^2), & x \le 14 \\ 0.05 \sin(x - 12.5) \exp(1 - (x - 12.5)^2) - 0.1, & x > 14. \end{cases}$$

in a computational domain of [0, 25] m with subcritical inlet/outlet at the left and right boundary, respectively. The bathymetry consists of a sinusoidal bump with a step of 0.1 m placed at $x = 14$ m. Our goal in this test case is to demonstrate the exact preservation of water at rest solutions when we consider a polynomial approximation of H: see section 3.2. Similar results are obtained when H is approximated using the selected multistep method. To do that, (52) is taken as initial condition and the numerical methods are run until time $t = 2$ sec. The differences between the numerical solutions obtained with WENO3GF-AMq, $q = 4, 6, 8$ and the initial condition are shown in Table 9: as it can be seen, the initial condition is preserved up to machine accuracy.

### 5.2.2 Perturbation of the lake at rest solution over a discontinuous bottom

In this test case, we consider an initial condition which is obtained by adding a small perturbation to the initial solution considered in 5.2.1. The perturbation is of size $\Delta h = 10^{-4}$ and is added in [7.5, 9.5] m. Figure 7, shows the initial
<figure><img src="image_7.png" /><figcaption>Figure 6: Burgers' equation. Test 5.1.3: discontinuous data. Perturbation of an isolated discontinuity: solutions obtained with WENO3GF-AMq, $q = 4, 6$ at time $t = 0.3$ (left) and $t = 0.5$ (right) using $N = 100$ (top) or $N = 300$ points (down) compared with a reference solution and with the numerical results given by WENO3-nWB using 300 points.</figcaption></figure>

<table><caption>Table 9: Shallow water equations. Test 5.2.1: lake at rest solution. Errors in $L^1$ norm for WENO3GF-AMq, $q = 4, 6, 8$ at time $t = 2$ sec.</caption><thead><tr><th>N</th><th>WENO3GF-AM4</th><th>WENO3GF-AM6</th><th>WENO3GF-AM8</th></tr></thead><tbody><tr><td>25</td><td>1.776e-15</td><td>6.661e-15</td><td>2.042e-14</td></tr><tr><td>50</td><td>1.476e-14</td><td>3.586e-14</td><td>3.907e-14</td></tr><tr><td>100</td><td>2.953e-14</td><td>7.382e-14</td><td>8.387e-14</td></tr><tr><td>200</td><td>6.586e-14</td><td>1.691e-13</td><td>2.054e-13</td></tr><tr><td>400</td><td>1.426e-13</td><td>3.784e-13</td><td>4.344e-13</td></tr></tbody></table>

perturbation on the left and the numerical solutions obtained with WENO3GF-AMq, $q = 4, 6, 8$ at time $t = 1$ sec using $N = 100$ points. They are compared with the numerical solution given by WENO3-nWB and a reference solution computed with WENO7GF-AM8. As it can be seen, the GF methods capture correctly the evolution of the initial perturbation while standard WENO3 give errors that are several orders of magnitude bigger than the initial perturbation.

We now change the size of the perturbation, and we place a Riemann problem instead. $\Delta h$ is 1.0 and it is added at the left of $x = 12m$: see Fig. 8 on the left. The numerical solutions for different schemes, after $t = 1$ sec, are plotted on the right. Again, the evolution of the perturbation is better captured by the WB methods as one may expect, especially due to the presence of the discontinuity in $H$.
<figure><img src="image_8.png" /><figcaption>Figure 7: Shallow water equations. Test 5.2.2: small perturbation on a lake at rest. Left: initial condition. Right: differences between the underlying stationary solution and the numerical solutions obtained with WENO3GF-AM$q$, $q = 4, 6, 8$ using $N = 100$ points at time $t = 1$ sec compared with a reference solution and the numerical solutions given by WENO3-nWB. Top: $η$. Down: $q$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Shallow water equations. Test 5.2.2: big perturbation on a lake at rest. Left: initial condition. Right: differences between the underlying stationary solution and the numerical solutions obtained with WENO3GF-AM$q$, $q = 4, 6, 8$ using $N = 100$ points at time $t = 1$ sec compared with a reference solution and the numerical solutions given by WENO3-nWB. Top: $η$. Down: $q$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Shallow water equations. Test 5.2.3: discrete steady state obtained with WENO3GF-AM4 using $N = 100$ points. Left: supercritical case. Right: subcritical case. Top: free surface and bottom. Down: $h$.</figcaption></figure>

### 5.2.3 Approximation of subcritical and supercritical stationary solutions

We consider the shallow water system (1)-(3) with the bottom depth given by:

$$H(x) = 0.05 \sin(x - 12.5) \exp \left(1 - (x - 12.5)^2\right).$$

As initial condition, we take the subcritical stationary solution characterized by

$$q^*(0) = 4.42, \quad h^*(0) = 2,$$

whose values at the points are computed using the relations (8). The numerical methods then run until a discrete steady state is reached up to machine accuracy. The boundary conditions used are:

$$q(0, t) = 4.42, \quad h(25, t) = 2.$$

Figure 9 (right) shows the discrete stationary solution reached by WENO3GF-AM4 using $N = 100$ points.

Table 10 shows the errors in $L^1$ norm and the empirical order of convergence of the different methods using WENO3GF-AMq, $q = 4, 6, 8$. For reference, the errors corresponding to WENO3-nWB are also shown. We should mention that, as observed in the case of Test 5.1.2, the empirical orders and errors for WENOpGF-ABq are independent of p, which is why we are only presenting the table for WENO3. As it can be seen, we confirm that, also for a system of equations, the order of accuracy for stationary solutions is the one of the multi-step method.

<table><caption>Table 10: Shallow water equation. Test 5.2.3: subcritical steady state. Errors in $L^1$ norm and convergence rate for WENO3GF-AMq, $q = 4, 5, 8$ (left) and WENO3-nWB (right).</caption><thead><tr><th rowspan="2">N</th><th colspan="2">4 steps</th><th colspan="2">6 steps</th><th colspan="2">8 steps</th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-weight: bold;">WENO3GF - AM</td></tr><tr><td>25</td><td>2.507e-03</td><td>-</td><td>3.510e-03</td><td>-</td><td>4.861e-03</td><td>-</td></tr><tr><td>50</td><td>1.032e-04</td><td>4.6</td><td>3.700e-04</td><td>3.2</td><td>8.815e-04</td><td>2.5</td></tr><tr><td>100</td><td>1.209e-05</td><td>3.1</td><td>1.098e-05</td><td>5.0</td><td>1.481e-05</td><td>5.9</td></tr><tr><td>200</td><td>8.110e-07</td><td>3.9</td><td>2.298e-07</td><td>5.6</td><td>9.337e-08</td><td>7.3</td></tr><tr><td>400</td><td>5.279e-08</td><td>4.0</td><td>3.817e-09</td><td>5.9</td><td>4.181e-10</td><td>7.8</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">WENO3GF - AB</td></tr><tr><td>25</td><td>1.033e-02</td><td>-</td><td>7.332e-02</td><td>-</td><td>1.555e-01</td><td>-</td></tr><tr><td>50</td><td>1.990e-03</td><td>2.4</td><td>6.356e-03</td><td>3.5</td><td>1.837e-02</td><td>3.0</td></tr><tr><td>100</td><td>1.326e-04</td><td>3.9</td><td>1.583e-03</td><td>2.0</td><td>4.255e-04</td><td>5.4</td></tr><tr><td>200</td><td>1.054e-05</td><td>3.7</td><td>2.940e-05</td><td>5.7</td><td>2.851e-06</td><td>7.2</td></tr><tr><td>400</td><td>6.937e-07</td><td>3.9</td><td>4.185e-09</td><td>6.3</td><td>8.015e-08</td><td>8.0</td></tr></tbody></table>

<table><thead><tr><th colspan="3" style="text-align: center; font-weight: bold;">WENO3-nWB</th></tr><tr><th>N</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>25</td><td>2.157e-01</td><td>-</td></tr><tr><td>50</td><td>1.349e-02</td><td>4.0</td></tr><tr><td>100</td><td>1.840e-03</td><td>2.9</td></tr><tr><td>200</td><td>2.308e-04</td><td>3.0</td></tr><tr><td>400</td><td>2.882e-05</td><td>3.0</td></tr></tbody></table>

Next, we consider as a safety check the same problem with a reversed flow direction. Note that we do not change the ODE weights in the scheme to integrate the source term. Table 11 confirms that, despite their non-symmetric
structure, the GF schemes obtained with both multi-step methods are able to achieve the expected order of accuracy.

The same tests are performed for the supercritical stationary solution characterized by

$$h^*(0) = 2, \quad q^*(0) = 24.$$

The left picture on Figure 9 shows the discrete stationary solution reached by WENO3GF-AM4 using $N = 100$ points. The convergence results are reported on table 12, and show the same behaviour also observed for the subcritical case, while Table 14 reports the same for the reversed supercritical flow confirming once again that the results are independent of the flow direction.

<table><caption>Table 11: Shallow water equation. Test 5.2.3: subcritical reversed steady state. Errors in $L^1$ norm and convergence rate for WENO3GF-AMq, $q = 4, 5, 8$ (left) and WENO3-nWB (right).</caption><thead><tr><th></th><th colspan="2">4 steps</th><th colspan="2">6 steps</th><th colspan="2">8 steps</th></tr><tr><th></th><th colspan="6">WENO3GF - AM</th></tr><tr><th>N</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>50</td><td>1.033e-04</td><td>-</td><td>3.633e-04</td><td>-</td><td>8.636e-04</td><td>-</td></tr><tr><td>100</td><td>1.217e-05</td><td>3.1</td><td>1.098e-05</td><td>5.0</td><td>1.470e-05</td><td>5.9</td></tr><tr><td>200</td><td>8.151e-07</td><td>3.9</td><td>2.289e-07</td><td>5.6</td><td>9.326e-08</td><td>7.3</td></tr><tr><td>400</td><td>5.317e-08</td><td>3.9</td><td>3.809e-09</td><td>5.9</td><td>4.153e-10</td><td>7.8</td></tr><tr><th></th><th colspan="6">WENO3GF - AB</th></tr><tr><td>50</td><td>1.119e-03</td><td>-</td><td>5.626e-03</td><td>-</td><td>1.805e-02</td><td>-</td></tr><tr><td>100</td><td>1.329e-04</td><td>3.8</td><td>2.327e-04</td><td>5.0</td><td>4.228e-04</td><td>5.4</td></tr><tr><td>200</td><td>1.056e-05</td><td>3.7</td><td>4.745e-06</td><td>5.6</td><td>2.849e-06</td><td>7.2</td></tr><tr><td>400</td><td>6.989e-07</td><td>3.9</td><td>8.351e-08</td><td>5.8</td><td>1.300e-08</td><td>7.8</td></tr></tbody></table>

<table><caption>Table 12: Test 5.2.3: Supercritical flow. Errors in $L^1$ norm and convergence rate for WENO3</caption><thead><tr><th colspan="3">WENO3-nWB</th></tr><tr><th>N</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>50</td><td>4.594e-02</td><td>-</td></tr><tr><td>100</td><td>7.163e-03</td><td>2.7</td></tr><tr><td>200</td><td>1.620e-03</td><td>2.1</td></tr><tr><td>400</td><td>2.380e-04</td><td>2.8</td></tr></tbody></table>

<table><thead><tr><th></th><th colspan="2">4 steps</th><th colspan="2">6 steps</th><th colspan="2">8 steps</th></tr><tr><th>N</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><th colspan="7">WENO3GF - AM</th></tr><tr><td>25</td><td>4.270e-04</td><td>-</td><td>6.847e-04</td><td>-</td><td>1.049e-03</td><td>-</td></tr><tr><td>50</td><td>4.181e-05</td><td>3.4</td><td>1.296e-04</td><td>2.4</td><td>3.310e-04</td><td>1.6</td></tr><tr><td>100</td><td>4.113e-06</td><td>3.3</td><td>4.210e-06</td><td>4.9</td><td>5.289e-06</td><td>6.0</td></tr><tr><td>200</td><td>2.749e-07</td><td>3.9</td><td>8.458e-08</td><td>5.6</td><td>3.524e-08</td><td>7.2</td></tr><tr><td>400</td><td>1.778e-08</td><td>4.0</td><td>1.411e-09</td><td>5.9</td><td>1.628e-10</td><td>7.8</td></tr><tr><th colspan="7">WENO3GF - AB</th></tr><tr><td>25</td><td>2.166e-03</td><td>-</td><td>5.047e-03</td><td>-</td><td>1.226e-02</td><td>-</td></tr><tr><td>50</td><td>6.085e-04</td><td>1.8</td><td>2.124e-03</td><td>1.2</td><td>7.079e-03</td><td>0.7</td></tr><tr><td>100</td><td>5.089e-05</td><td>3.6</td><td>8.074e-05</td><td>4.7</td><td>1.394e-04</td><td>5.7</td></tr><tr><td>200</td><td>3.646e-06</td><td>3.8</td><td>1.777e-06</td><td>5.5</td><td>1.092e-06</td><td>7.0</td></tr><tr><td>400</td><td>2.336e-07</td><td>4.0</td><td>3.072e-08</td><td>5.9</td><td>4.928e-09</td><td>7.8</td></tr></tbody></table>

<table><thead><tr><th colspan="3">WENO3-nWB</th></tr><tr><th>N</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>25</td><td>1.678e-01</td><td>-</td></tr><tr><td>50</td><td>2.723e-03</td><td>2.6</td></tr><tr><td>100</td><td>4.562e-04</td><td>2.6</td></tr><tr><td>200</td><td>8.543e-05</td><td>2.4</td></tr><tr><td>400</td><td>1.108e-05</td><td>2.9</td></tr></tbody></table>

### 5.2.4 Perturbations subcritical and supercritical stationary solutions

As in the case of Test 5.2.2, we consider now an initial condition that represents a perturbation of a moving stationary solution. For brevity we consider only the subcritical case of section 5.2.3. The results for the supercritical are very similar.

First, we consider the evolution of a small perturbation $\Delta h = 10^{-4}$ m, which is added to the depth in the subdomain $[7.5, 9.5]$ m. The left picture on figure 10 shows the numerical results obtained with WENO3GF-AMq using a mesh of $N = 100$ points at time $t = 1$ sec. For comparison we report the results of the non-well balanced WENO3 scheme, and a reference solution obtained on 500 points with the WENO7GF-AM8 scheme. We can see that on this coarse mesh the well balanced formulations allow to nicely follow the evolution of such a small perturbation.

For completeness we also consider a perturbation of $\Delta h = 1$ m, which is added for all points to the left of $x = 12$ m. For the smooth bathymetry used in this case all schemes provide a similar evolution of the solution.
<table><caption>Table 13: Shallow water equation. Test 5.2.3: supercritical reversed steady state. Errors in $L^1$ norm and convergence rate for WENO3GF-AM$q$, $q = 4, 5, 8$ (left) and WENO3-nWB (right).</caption><thead><tr><th></th><th colspan="2">4 steps</th><th colspan="2">6 steps</th><th colspan="2">8 steps</th></tr><tr><th>N</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><th colspan="7" style="text-align: center;">WENO3GF - AM</th></tr><tr><td>50</td><td>4.187e-05</td><td>-</td><td>1.303e-04</td><td>-</td><td>3.314e-04</td><td>-</td></tr><tr><td>100</td><td>4.128e-06</td><td>3.4</td><td>4.212e-06</td><td>5.0</td><td>5.292e-06</td><td>6.0</td></tr><tr><td>200</td><td>2.759e-07</td><td>3.9</td><td>8.455e-08</td><td>5.6</td><td>3.524e-08</td><td>7.2</td></tr><tr><td>400</td><td>1.784e-08</td><td>4.0</td><td>1.408e-09</td><td>5.9</td><td>1.618e-10</td><td>7.8</td></tr><tr><th colspan="7" style="text-align: center;">WENO3GF - AB</th></tr><tr><td>50</td><td>6.097e-04</td><td>-</td><td>2.125e-03</td><td>-</td><td>7.092e-03</td><td>-</td></tr><tr><td>100</td><td>5.090e-05</td><td>3.6</td><td>8.095e-05</td><td>4.7</td><td>1.395e-04</td><td>5.6</td></tr><tr><td>200</td><td>3.657e-06</td><td>3.8</td><td>1.776e-06</td><td>5.5</td><td>1.092e-06</td><td>7.0</td></tr><tr><td>400</td><td>2.344e-07</td><td>4.0</td><td>3.070e-08</td><td>5.9</td><td>4.928e-09</td><td>7.8</td></tr></tbody></table>

<table><thead><tr><th colspan="3">WENO3-nWB</th></tr><tr><th>N</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>50</td><td>2.397e-03</td><td>-</td></tr><tr><td>100</td><td>2.996e-04</td><td>3.0</td></tr><tr><td>200</td><td>3.879e-05</td><td>2.9</td></tr><tr><td>400</td><td>4.910e-06</td><td>3.0</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Figure 10: Shallow water equations. Test 5.2.4: perturbation of a subcritical steady state. Left: small perturbation. Right: big perturbation. Differences between the underlying stationary solution and the numerical solutions obtained with WENO3GF-AM$q$, $q = 4, 6, 8$ using $N = 100$ points at time $t = 1$ sec compared with a reference solution and the numerical solutions given by WENO3-nWB. Top row: $\eta$. Bottom row: $q$.</figcaption></figure>

### 5.2.5 Approximation of a transcritical stationary solution

For completeness we also report the results obtained for a transcritical case. To begin with we compute the steady states using the AB ODE solver of different accuracies and perform a grid convergence to the analytical solution obtained solving the exact cubic relation associated to (8). The results are reported on table 14 in which we compare the accuracy of the steady WENO5FG-AB for different accuracies, with those of the standard non well balanced upwind WENO5. Note that we do not modify the flux splitting to account for the presence of the sonic point. The tables show that on the mesh sizes considered also for this case the errors obtained using the global flux quadrature formulation are of orders of magnitude lower than those of the upwind scheme.

<table><caption>Table 14: Shallow water equation. Test 5.2.5: transcritical flow over a smooth bump. Errors in $L^1$ norm and convergence rate for WENO5GF-AB$q$, $q = 4, 6, 8$ (left) and WENO5-nWB (right).</caption><thead><tr><th></th><th colspan="2">4 steps</th><th colspan="2">6 steps</th><th colspan="2">8 steps</th></tr><tr><th>N</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><th colspan="7" style="text-align: center;">WENO5GF - AB</th></tr><tr><td>50</td><td>0.0991049</td><td>-</td><td>0.2174215</td><td>-</td><td>0.2381622</td><td>-</td></tr><tr><td>100</td><td>0.0066493</td><td>3.91</td><td>0.0088290</td><td>4.64</td><td>0.0062804</td><td>5.26</td></tr><tr><td>200</td><td>0.0004984</td><td>3.75</td><td>0.0001740</td><td>5.68</td><td>4.317e-05</td><td>7.21</td></tr><tr><td>400</td><td>4.483e-05</td><td>3.49</td><td>3.375e-06</td><td>5.71</td><td>1.101e-06</td><td>5.31</td></tr></tbody></table>

<table><thead><tr><th colspan="3">WENO5-nWB</th></tr><tr><th>N</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>50</td><td>0.2305202</td><td>-</td></tr><tr><td>100</td><td>0.0913807</td><td>1.33</td></tr><tr><td>200</td><td>0.0194528</td><td>2.24</td></tr><tr><td>400</td><td>0.0004357</td><td>5.50</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>Figure 11: Shallow water equations. Test 5.2.6: perturbation of a transcritical discrete steady state. Left: Discrete steady solution. Right: Differences between the underlying discrete stationary solution and the numerical solutions obtained with WENO5GF-ABq, q = 4, 6, 8 using N = 100 points at time t = 0.7 sec compared with a reference solution and the numerical solutions given by WENO5-nWB. Top row: η. Bottom row: q.</figcaption></figure>

### 5.2.6 Perturbation of a transcritical stationary solution

We also consider the evolution of perturbations of the smooth transcritical flow of the previous section. For each scheme we consider the perturbation of the well prepared data corresponding to its discrete steady state. In the case of the GF-AB schemes this is obtained directly using the ODE solver. We add to this data a small perturbation of $Δh = 10^{-4}$ m. Figure 11 shows on the left the discrete steady solution obtained using the ODE AB and on the right the numerical results obtained with WENO5-ABq,q = 4, 6, 8. We used a mesh of N = 100 points at time t = 0.7 sec. We also report the results of the WENO5-nWB and a reference solution obtained with N = 1000 points and using the WENO7-AM8 scheme. Like before we confirm that the well balanced formulation allow to nicely follow the evolution of the small perturbation.

### 5.2.7 Moving equilibria with bathymetry and friction

**Analytical solutions with $κhq|q|$ friction.** We consider here (4) with the friction closure

$$κ = kh|q|.$$

In this case, assuming that the flow runs left to right, and setting $η = h - H$, the steady ODE (10) can be recast as

$$-C_1^2 h^{-2} h_x + g h \eta_x = -k C_1^2 h$$

which, dividing through by h gives

$$-C_1^2 h^{-3} h_x + g \eta_x = -k C_1^2. \quad (53)$$

This equation can be integrated to obtain the equalities

$$q = C_1, \quad \frac{1}{2} \frac{q^2}{h^2} + gh - gH - kq^2 x = C_2.$$

Given $C_i$, $i = 1, 2$ the point values of a stationary solution can be then computed by solving a cubic equation like in the case without friction.

To obtain a fully analytical expression we proceed in a slightly different manner. We impose a prescribed profile $z(x)$ to the free surface elevation $η$, and then compute $h$ by solving (53). A compatible definition of the depth at rest $H$ is obtained as $H = h - z$. If the boundary condition

$$h(x_0) = h_0, \quad q(x_0) = q_0$$
<figure><img src="image_13.png" /><figcaption>Figure 12: Shallow water equations with friction law $\kappa hq|q|$. Supercritical (left) and subcritical (right) exact solutions. Top: free surface and bottom. Down: $h$.</figcaption></figure>

is imposed, the corresponding stationary solution is given by

$$ 
\begin{aligned}
q^e(x) &= q_0, \\
h^e(x) &= \frac{h_0}{\sqrt{1 + \frac{2h_0^2g}{q_0^2}(z(x_0) - z(x)) - 2kh_0^2(x - x_0)}}, \\
H(x) &= h^e(x) - z^e(x).
\end{aligned} 
\quad (54) $$

We now consider the free surface profiles

$$ z(x) = a - b * \frac{cxe^{\cos(4\pi x)} - 1/e}{e - 1/e} \quad (55) $$

One can show that a supercritical solution is obtained for

$$ (h_0, q_0, a, b, c, k) = (1, 3/2, 5/2, 1/2, 2, 0.3), $$

while a subcritical one is obtained for

$$ (h_0, q_0, a, b, c, k) = (1, 0.3, 5/2, 1/4, 1/2, 1/2). $$

The two solutions are plotted for completeness in Figure 12.

We start by computing the grid convergence of the discrete steady state to the exact one. In Table 15 the $L^1$ errors with respect to the analytical solution for WENO3GF-AMq, $q = 4, 6, 8$ are compared with the corresponding to the standard WENO3 and WENO5 schemes. Similar results are obtained for WENO5GF-AMq or WENOp-ABq, $p = 3, 5, q = 4, 6, 8$. The table shows that also for friction source terms, the accuracy of the underlying ODE solver is retrieved. Compared to the standard WENO3 method, this leads to error reductions which go from 2 to 6 orders of magnitude on the finest resolutions and of 2 or 3 orders of magnitude on the intermediate one. Similar results are obtained for the subcritical stationary solution.

We now add the perturbation $\Delta h = 10^{-4}$ in $[0.1, 0.2]$ to the supercritical stationary solution and run the simulations until the final time $T = 0.08$. A reference solution is computed on $N = 2000$ points with WENO5GF-AM8. Figure 13 shows the numerical solution obtained with WENO3GF-AM4 and WENO3-nWB using $N = 100$ points: as it can be seen, the former captures well the evolution of the perturbation while the latter gives errors of orders of magnitude higher than the initial perturbation on this relatively coarse mesh. For comparison we also report the WENO3 solution on a much refined mesh with $N = 400$ which, of course, follows more accurately the reference. The conclusions are similar for the numerical methods based on WENO5 and/or different ODE solvers as well as for the perturbation of the subcritical stationary solution.
<table><caption>Table 15: Shallow water equations with friction law $\kappa h q |q|$. Supercritical stationary solution: errors in $L^1$ norm for $H$ and convergence rates for WENO3GF-AM$q$, with $q = 4, 6, 8$ and WENO$p$-nWB, $p = 3, 5$.</caption><thead><tr><th colspan="3">4 steps</th><th colspan="2">6 steps</th><th colspan="2">8 steps</th></tr><tr><th>N</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><th colspan="7">WENO3GF - AM</th></tr><tr><td>20</td><td>3.789e-03</td><td>-</td><td>2.960e-03</td><td>-</td><td>2.360e-03</td><td>-</td></tr><tr><td>40</td><td>1.671e-04</td><td>6.3</td><td>9.639e-05</td><td>5.0</td><td>5.734e-05</td><td>5.4</td></tr><tr><td>80</td><td>6.966e-06</td><td>4.6</td><td>1.225e-06</td><td>6.3</td><td>3.060e-07</td><td>7.6</td></tr><tr><td>160</td><td>3.085e-07</td><td>4.5</td><td>1.379e-08</td><td>6.5</td><td>1.184e-09</td><td>8.0</td></tr><tr><td>320</td><td>1.504e-08</td><td>4.4</td><td>1.842e-10</td><td>6.2</td><td>4.554e-12</td><td>8.0</td></tr></tbody></table><table><thead><tr><th colspan="3">WENO3nWB</th><th colspan="3">WENO5nWB</th></tr><tr><th>N</th><th>Error</th><th>Order</th><th>N</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>7.606e-03</td><td>-</td><td>20</td><td>5.706e-03</td><td>-</td></tr><tr><td>40</td><td>1.027e-03</td><td>2.9</td><td>40</td><td>3.913e-04</td><td>3.9</td></tr><tr><td>80</td><td>1.364e-04</td><td>2.9</td><td>80</td><td>2.066e-05</td><td>4.3</td></tr><tr><td>160</td><td>1.743e-05</td><td>3.0</td><td>160</td><td>6.632e-07</td><td>5.0</td></tr><tr><td>320</td><td>2.195e-06</td><td>3.0</td><td>320</td><td>2.037e-08</td><td>5.0</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>Figure 13: Shallow water equations with friction law $\kappa h q |q|$. Small perturbation of the supercritical steady state. Left: initial perturbation and differences between the reference and the stationary solutions. Right: differences between the stationary solution and the numerical solutions obtained with WENO3GF-AM4 using $N = 100$ points at time $t = 0.08$ sec compared with the reference solution and the numerical solutions given by WENO3-nWB with $N = 100, 400$ point. Top: $\eta$. Down: $q$.</figcaption></figure>

**Moving solutions with classical Manning friction**

Last, we consider the shallow water equations with the classical Manning friction closure

$$ \kappa = \frac{k|q|}{h^{\mu}} $$

where $k$ is the Manning friction coefficient and $\mu$ is set to $7/3$. The Manning constant used is set to 0.05.

As in previous cases, moving equilibria can still be achieved even when friction with a constant Manning coefficient is present. In the supercritical case, the friction term causes a reduction in physical speed from left to right, resulting in an increase in the free surface elevation across the same direction.

In this case, we first compute the discrete steady-state solution with a sweep of the Adams-Bashforth multi-step solver (see section 4).

The solution obtained with the AB6 scheme is reported on the left of figure 14. We use the ODE solutions as well prepared initial data. Table 16 presents the errors in $L^1$ norm for the WENO5-AB$q$ with $q = 4, 6, 8$ obtained at time $t = 2$ sec. As it can been seen the discrete steady solution is preserved up to machine accuracy.
<table><caption>Table 16: Shallow water equations with Manning friction law. Supercritical discrete stationary solution. Errors in $L^1$ norm for WENO5GF-AB$\tilde{q}$, $q = 4, 6, 8$ at time $t = 2$ sec.</caption><thead><tr><th>N</th><th>WENO5GF-AB6</th><th>WENO5GF-AB8</th></tr></thead><tbody><tr><td>20</td><td>1.365e-13</td><td>3.912e-10</td></tr><tr><td>40</td><td>2.364e-13</td><td>1.620e-13</td></tr><tr><td>80</td><td>2.033e-13</td><td>3.140e-13</td></tr><tr><td>160</td><td>2.516e-13</td><td>6.079e-13</td></tr><tr><td>320</td><td>1.724e-12</td><td>1.481e-12</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>Figure 14: Perturbation of the discrete steady solutions with Manning friction 0.05, supercritical case. Left: initial discrete steady solution (Weno5-GF-AB6). Right: comparison between reference and the solutions at time $t = 0.7$ obtained with the WENO5-GF-AB6 and WENO5 on $N = 100$ points, and WENO5 on $N = 400$ points.</figcaption></figure>

Next, we perturb the discrete steady solutions of WENO5-GF-AB6 with a small perturbation of $\Delta h = 0.0001$ m placed in the interval [7.5, 9.5] m. Figure 14 (right) compares the scheme's solution, on a supercritical flow, with a reference solution obtained using WENO7GF-AB8, and WENO5 using $N = 100$ points. The figure shows that the error of the WENO5 is of the order of the perturbation resulting to degrade the solution quality. For comparison we also plot the WENO5 solution on a refined mesh with $N = 400$ points.

A similar behavior is observed in the computation of the discrete steady solution for the subcritical case; therefore, we omit the table. In this case, the water depth is expected to decrease from left to right, while the speed increases, see figure 15 on the left. Using the discrete steady solution as well-prepared data for the initialization we add the perturbation and obtain the numerical results at $t = 1$ sec.

# 6 Conclusion

In this work, we proposed high-order well-balanced methods for one-dimensional hyperbolic systems of balance laws. These methods are designed to ensure that the discrete steady states align with solutions of high-order ODE integrators, without explicitly solving the ODEs. By embedding this property directly into the scheme, we avoid applying the ODE integrator to solve the local Cauchy problem. A sweep of the ODE solution can still be used to obtain well prepared initial data if required.

To achieve this result we have combined a WENO finite difference framework, with a global flux quadrature strategy in which the source quadrature employs weights from high-order multi-step ODE methods on the finite difference grid. This approach ensures a well-balanced treatment of the source term and a direct connection between the scheme's discrete solutions and the underlying ODE integrator. In practice, we used WENO reconstructions of orders 3 to 7 with Adams-Moulton and Adams-Bashforth methods of up to order 8.

We performed classical benchmarks for two prototype equations involving moving equilibria to verify the correct implementation and convergence of our methods. For time-dependent solutions, we demonstrated that the order of
<figure><img src="image_16.png" /><figcaption>Figure 15: Perturbation of the discrete steady solutions with Manning friction 0.05, subcritical case. Left: initial discrete steady solution (Weno5-GF-AB6). Right: comparison between reference and the solutions at time $t = 1.0$ obtained with the WENO5-GF-AB6 and WENO5 on $N = 100$ points, and WENO5 on $N = 400$ points.</figcaption></figure>

accuracy is determined by the minimum between the WENO reconstruction order and the ODE integrator order. For smooth steady-state solutions, both the order of accuracy and the error values are uniquely determined by the choice of the multi-step method. This allows for significant error reduction—by several orders of magnitude on a given mesh—simply by adjusting the weights of the quadrature formula.

We also extended our methods to handle singular source terms and to preserve the water-at-rest solution for the shallow water equations. For the latter system, we proposed a strategy to find a particular family of stationary solutions, both with and without friction.

## Acknowledgements

The work of Carlos Parés is partially supported by grants PID2022-137637NB-C21 and PID2022-137637NB-C22 funded by MCIN/AEI/10.13039/501100011033 and by "ERDF A way of making Europe". Maria Kazolea and Mario Ricchiuto are members of the CARDAMOM research team of the Inria center at the University of Bordeaux.

## In memoriam

This paper is dedicated to the memory of Prof. Arturo Hidalgo López (*July 03rd 1966 - †August 26th 2024) of the Universidad Politecnica de Madrid, organizer of HONOM 2019 and active participant in many other editions of HONOM. Our thoughts and wishes go to his wife Lourdes and his sister María Jesús, whom he left behind.

## A Coefficients for the AM and AB schemes

For the sake of completeness, we include here the coefficients of the multi-step methods that are used to compute the integrals of the source term and the discrete stationary solutions. Since we work with schemes of formal accuracy $2k + 1 \ge 3$ we consider ODE solvers whose solutions have accuracy at least $2k + 2 \ge 4$. The coefficients for the relevant methods are provided below up to order $8^1$.

<footnote n="1">¹for higher orders check the tables/codes here: http://www.mymathlib.com/diffeq/adams/</footnote>
**ABq methods** The explicit ABq methods of order q = 4, 6, and 8 use q steps and βq = 0. The remaining coefficients are:

AB8:

$$
\beta_7 = \frac{434241}{120960}
$$

AB6:

$$
\beta_5 = \frac{4277}{1440}
$$

$$
\beta_6 = - \frac{1152169}{120960}
$$

AB4:

$$
\beta_3 = \frac{55}{24}
$$

$$
\beta_4 = - \frac{7923}{1440}
$$

$$
\beta_5 = \frac{2183877}{120960}
$$

$$
\beta_2 = - \frac{59}{24}, \quad \beta_3 = \frac{9982}{1440}
$$

$$
\beta_3 = \frac{9982}{1440},
$$

$$
\beta_4 = - \frac{2664477}{120960}
$$

$$
\beta_1 = \frac{37}{24}
$$

$$
\beta_2 = - \frac{7298}{1440}
$$

$$
\beta_3 = \frac{2102243}{120960}
$$

$$
\beta_0 = - \frac{9}{24}
$$

$$
\beta_1 = \frac{2877}{1440}
$$

$$
\beta_2 = - \frac{1041723}{120960}
$$

$$
\beta_0 = - \frac{475}{1440}
$$

$$
\beta_1 = \frac{295767}{120960}
$$

$$
\beta_0 = - \frac{36799}{120960}
$$

**AMq methods** The AMq methods of order q are implicit and use s = q − 1 steps. The coefficients are the following:

AM8:

$$
\beta_7 = \frac{36799}{120960}
$$

AM6:

$$
\beta_5 = \frac{475}{1440}
$$

$$
\beta_6 = \frac{139849}{120960}
$$

AM4:

$$
\beta_3 = \frac{9}{24}
$$

$$
\beta_4 = \frac{1427}{1440}
$$

$$
\beta_5 = - \frac{121797}{120960}
$$

$$
\beta_2 = \frac{19}{24}, \quad \beta_3 = - \frac{798}{1440}
$$

$$
\beta_3 = - \frac{798}{1440},
$$

$$
\beta_3 = \frac{123133}{120960}
$$

$$
\beta_1 = - \frac{5}{24}
$$

$$
\beta_2 = \frac{482}{1440}
$$

$$
\beta_3 = - \frac{88547}{120960}
$$

$$
\beta_0 = \frac{1}{24}
$$

$$
\beta_1 = - \frac{173}{1440}
$$

$$
\beta_2 = \frac{41499}{120960}
$$

$$
\beta_0 = \frac{27}{1440}
$$

$$
\beta_1 = - \frac{11351}{120960}
$$

$$
\beta_0 = \frac{1375}{120960}
$$

B Discontinuous source linearization

Let us show some examples of linearizations of the source term satisfying (41):
• For the Burgers' equation (1)-(2) we can compute simple linearizations of the source using satisfying (41) from the relations (6), (7) satisfied at an admissible jump:

$$-p = 1: \tilde{S}_{i+1/2} = \bar{U}_{i+1/2} := (U_i + U_{i+1})/2;$$

$$-p = 2: \tilde{S}_{i+1/2} = \frac{\bar{U}_{i+1/2}[U]_{i+1/2}}{\ln(U_{i+1}/U_i)};$$

$$-p \ge 3: \tilde{S}_{i+1/2} = (2-p)\bar{U}_{i+1/2} \frac{[U]_{i+1/2}}{[U^2-p]_{i+1/2}}$$

Note that all the singularities due to the division by jumps can be removed easily for $p \ge 3$. A careful implementation of the case $p = 2$ is required using the $\lim_{\alpha \to 0} \alpha / \ln(1 + \alpha) = 1$.

• For the shallow water equations two relations (8) are satisfied at an admissible jump

$$\begin{aligned} \llbracket q \rrbracket_{i+1/2} &= 0, \\ g \llbracket H \rrbracket_{i+1/2} &= g \llbracket h \rrbracket_{i+1/2} + \llbracket u^2/2 \rrbracket_{i+1/2}. \end{aligned}$$

Therefore one has

$$q_j = \bar{q}_{i+1/2}, \quad u_j = \frac{\bar{q}_{i+1/2}}{h_j}, \quad j = i, i+1.$$

We deduce that across an admissible jump

$$g \llbracket H \rrbracket_{i+1/2} = \left\{ g - \frac{\bar{q}_{i+1/2}^2 \bar{h}_{i+1/2}}{(h_i h_{i+1})^2} \right\} \llbracket h \rrbracket_{i+1/2}.$$

Imposing that

$$\llbracket h u^2 + gh^2/2 \rrbracket_{i+1/2} = g \tilde{h}_{i+1/2} \llbracket H \rrbracket_{i+1/2}$$

which means

$$-\frac{\bar{q}_{i+1/2}^2}{h_i h_{i+1}} \llbracket h \rrbracket_{i+1/2} + g \bar{h}_{i+1/2} \llbracket h \rrbracket_{i+1/2} = \tilde{h}_{i+1/2} \left\{ g - \frac{\bar{q}_{i+1/2}^2 \bar{h}_{i+1/2}}{(h_i h_{i+1})^2} \right\} \llbracket h \rrbracket_{i+1/2},$$

we end up with

$$\tilde{h}_{i+1/2} = \frac{g \bar{h}_{i+1/2} - \frac{\bar{q}_{i+1/2}^2}{h_i h_{i+1}}}{\left\{ g - \frac{\bar{q}_{i+1/2}^2 \bar{h}_{i+1/2}}{(h_i h_{i+1})^2} \right\}} = \bar{h}_{i+1/2} + \frac{\frac{\bar{q}_{i+1/2}^2}{g(h_i h_{i+1})^2} (\bar{h}_{i+1/2}^2 - h_i h_{i+1})}{1 - \frac{\bar{q}_{i+1/2}^2}{g(h_i h_{i+1})^2} \bar{h}_{i+1/2}}.$$

Then,

$$\tilde{S}_{i+1/2} = \begin{bmatrix} 0 \\ -g \tilde{h}_{i+1/2} \end{bmatrix}$$

is defined. If the equality (43) is used to compute the integral of the source term, the linearization would be instead as follows

$$\int_{x_i}^{x_{i+1}} g h H_x \approx g \tilde{\eta}_{i+1/2} \llbracket H \rrbracket_{i+1/2} - \frac{g}{2} \llbracket H^2 \rrbracket_{i+1/2}$$

with

$$\tilde{\eta}_{i+1/2} := \bar{\eta}_{i+1/2} + \frac{\frac{\bar{q}_{i+1/2}^2}{g(h_i h_{i+1})^2} (\bar{h}_{i+1/2}^2 - h_i h_{i+1})}{1 - \frac{\bar{q}_{i+1/2}^2}{g(h_i h_{i+1})^2} \bar{h}_{i+1/2}}$$
References

[1] ABGRALL, R., AND RICCHIUTO, M. *Hyperbolic Balance Laws: Residual Distribution, Local and Global Fluxes*. Springer Nature Singapore, Singapore, 2022, pp. 177–222.

[2] AUDUSSE, E., BOUCHUT, F., BRISTEAU, M.-O., KLEIN, R., AND PERTHAME, B. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. *SIAM Journal on Scientific Computing* 25, 6 (2004), 2050–2065.

[3] BARSUKOW, W., RICCHIUTO, M., AND TORLO, D. Structure preserving nodal continuous finite elements via global flux quadrature, 2024. to appear on Num.Meth. for PDEs.

[4] BERBERICH, J. P., CHANDRASHEKAR, P., AND KLINGENBERG, C. High order well-balanced finite volume methods for multi-dimensional systems of hyperbolic balance laws. *Computers & Fluids* 219 (2021), 104858.

[5] BERMÚDEZ, A., AND VÁZQUEZ-CENDÓN, M. E. Upwind methods for hyperbolic conservation laws with source terms. *Computers & Fluids* 23, 8 (1994), 1049–1071.

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

[7] BERTHON, C., MICHEL-DANSAC, V., AND THOMANN, A. An entropy-stable and fully well-balanced scheme for the Euler equations with gravity. working paper or preprint, June 2024.

[8] CASELLES, V., DONAT, R., AND HARO, G. Flux-gradient and source-term balancing for certain high resolution shock-capturing schemes. *Computers and Fluids* 38, 1 (2009), 16–36.

[9] CASTRO, M., GALLARDO, J., LÓPEZ-GARCÍA, J., AND PARÉS, C. Well-balanced high order extensions of godunov's method for semilinear balance laws. *SIAM Journal on Numerical Analysis* 46 (2008), 1012–1039.

[10] CASTRO, M., LÓPEZ-GARCÍA, J., AND PARÉS, C. High order exactly well-balanced numerical methods for shallow water systems. *Journal of Computational Physics* 246 (2013), 242–264.

[11] CASTRO, M., PARDO MILANÉS, A., AND PARÉS, C. Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique. *Mathematical Models and Methods in Applied Sciences* 17, 12 (2007), 2055–2113.

[12] CASTRO, M. J., MORALES DE LUNA, T., AND PARÉS, C. Chapter 6 - well-balanced schemes and path-conservative numerical methods. In *Handbook of Numerical Methods for Hyperbolic Problems*, R. Abgrall and C.-W. Shu, Eds., vol. 18. Elsevier, 2017, pp. 131–175.

[13] CASTRO, M. J., AND PARÉS, C. Well-balanced high-order finite volume methods for systems of balance laws. *Journal of Scientific Computing*, 82, 42 (2020).

[14] CHEN, G., AND NOELLE, S. A new hydrostatic reconstruction scheme based on subcell reconstructions. *SIAM Journal on Numerical Analysis* 55, 2 (2017), 758–784.

[15] CHENG, Y., CHERTOCK, A., HERTY, M., AND KURGANOV, A. A new approach for designing moving-water equilibria preserving schemes for the shallow water equations. *Journal of Scientific Computing* 80 (2016), 538–554.

[16] CHENG, Y., AND KURGANOV, A. Moving-water equilibria preserving central-upwind schemes for the shallow water equations. *Communications in Mathematical Sciences* 14 (01 2016), 1643–1663.

[17] CIALLELLA, M., TORLO, D., AND RICCHIUTO, M. Arbitrary high order weno finite volume scheme with flux globalization for moving equilibria preservation. *Journal of Scientific Computing* 96 (2023), 53.

[18] DAL MASO, G., LEFLOCH, P., AND MURAT, F. Definition and weak stability of nonconservative products. *Journal de Mathématiques Pures et Appliquées* 74, 6 (1995), 483–548.

[19] GABURRO, E., CASTRO, M. J., AND DUMBSER, M. Well-balanced Arbitrary-Lagrangian-Eulerian finite volume schemes on moving nonconforming meshes for the Euler equations of gas dynamics with gravity. *Monthly Notices of the Royal Astronomical Society* 477, 2 (03 2018), 2251–2275.

[20] GALLOUËT, T., HÉRARD, J.-M., AND SEGUIN, N. Some approximate godunov schemes to compute shallow-water equations with topography. *Computers & Fluids* 32, 4 (2003), 479–513.

[21] GASCÓN, L., AND CORBERÁN, J. Construction of second-order tvd schemes for nonhomogeneous hyperbolic conservation laws. *Journal of Computational Physics* 172, 1 (2001), 261–297.

[22] GÓMEZ-BUENO, I., CASTRO, M., PARÉS, C., AND RUSSO, G. Collocation methods for high-order well-balanced methods for systems of balance laws. *Mathematics* 9, 15 (2021), 1799.
[23] GOTTLIEB, S., AND SHU, C.-W. Total variation diminishing Runge-Kutta schemes. *Mathematics of Computation* 67 (1998), 73–85.

[24] GREENBERG, J. M., AND LEROUX, A. Y. A well-balanced scheme for the numerical processing of source terms in hyperbolic equations. *SIAM Journal on Numerical Analysis* 33, 1 (1996), 1–16.

[25] JIANG, G., AND SHU, C.-W. Efficient implementation of weighted eno schemes. *Journal of Computationas Physcs* 126 (1996), 202–228.

[26] J.P., B., KÄPPELI, R., CHANDRASHEKAR, P., AND KLINGENBERG, C. High order discretely well-balanced methods for arbitrary hydrostatic atmospheres. *Communications in Computational Physics* 30, 3 (2021), 666–708.

[27] LIE, I., AND NØRSETT, S. Superconvergence for multistep collocation. *Mathematics of Computation* 52 (1989), 65–79.

[28] MANTRI, Y., ÖFFNER, P., AND RICCHIUTO, M. Fully well-balanced entropy controlled discontinuous galerkin spectral element method for shallow water flows: Global flux quadrature and cell entropy correction. *Journal of Computational Physics* 498 (2024), 112673.

[29] MICALIZZI, L., RICCHIUTO, M., AND ABGRALL, R. Novel well-balanced continuous interior penalty stabi-lizations. *J.Sci.Comp. 100*, 14 (2024).

[30] MICHEL-DANSAC, V., BERTHON, C., CLAIN, S., AND FOUCHER, F. A well-balanced scheme for the shallow-water equations with topography or manning friction. *Journal of Computational Physics* 335 (2017), 115–154.

[31] PARÉS, C., AND PARÉS-PULIDO, C. Well-balanced high-order finite difference methods for systems of balance laws. *Journal of Computational Physics* 425 (2021), 109880.

[32] ROE, P. Upwind differencing schemes for hyperbolic conservation laws with source terms. In *Nonlinear Hyperbolic Problems* (Berlin, Heidelberg, 1987), C. Carasso, D. Serre, and P.-A. Raviart, Eds., Springer Berlin Heidelberg, pp. 41–51.

[33] RUSSO, G., AND KHE, A. High order well-balanced schemes based on numerical reconstruction of the equilib-rium variables. *Journal of Computational Physics* 246 (2013), 242–264.

[34] SHENG, C.-T., WANG, Z.-Q., AND GUO, B.-Y. A multistep legendre-gauss spectral collocation method for nonlinear volterra integral equations. *SIAM Journal on Numerical Analysis* 52, 4 (2014), 1953–1980.

[35] ZIYAO, X., AND SHU., C.-W. A high-order well-balanced alternative finite difference weno (a-weno) method with the exact conservation property for systems of hyperbolic balance laws. *Preprint* (2024).