# How to make fully discrete entropy stability of three-point well-balanced finite-volume schemes

Christophe Berthon, Manuel Castro Díaz, Ludovic Martaud, Tomas Morales de Luna

► To cite this version:

Christophe Berthon, Manuel Castro Díaz, Ludovic Martaud, Tomas Morales de Luna. How to make fully discrete entropy stability of three-point well-balanced finite-volume schemes. 2025. hal-05378579

HAL Id: hal-05378579

https://hal.science/hal-05378579v1

Preprint submitted on 23 Nov 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
How to make fully discrete entropy stability of three-point
well-balanced finite-volume schemes

Christophe Berthon *

Manuel Castro Díaz †

Ludovic Martaud ‡

Tomas Morales de Luna†

July 24, 2025

Abstract

This paper presents a simple and robust method to enforce discrete entropy stability in first-
order well-balanced finite volume schemes for systems of balance laws, including non-conservative
terms. By leveraging an artificial viscosity technique originally introduced by Tadmor and extended
in recent works, the authors propose an entropy-preserving modification applicable to a wide class
of three-point schemes. The core contribution lies in the construction of a viscosity coefficient that
maintains the well-balanced property while ensuring entropy dissipation at the discrete level. A
theoretical framework is developed to guarantee robustness, well-balancing, and entropy stability
under a suitable CFL condition. The approach is validated through numerical experiments, in-
cluding applications to shallow water systems and two-layer flows, demonstrating accuracy and
stability.

# 1 Introduction

This work is devoted to the numerical approximation of systems of balanced laws, eventually with
non-conservative products, in the following form:

$$
\partial_t w + \partial_x f(w) + A(w) \partial_x w = S(w, x), \quad x \in \mathbb{R}, \quad t > 0. \tag{1}
$$

The unknown state vector $w \in \mathbb{R}^d$, for $d \ge 1$, is assumed to belong to an open convex set $\Omega \subset \mathbb{R}^d$, while $f: \Omega \to \mathbb{R}^d$ is a given flux function, $A: \Omega \to \mathcal{M}_d(\mathbb{R})$ denotes a given smooth matrix and $S: \Omega \times \mathbb{R} \to \mathbb{R}^d$ represents a given source term function. For sake of simplicity, we may enforce the system (1) to be hyperbolic so that $(\nabla f + A)(w)$ is assumed to be diagonalizable in $\mathbb{R}$ for all $w \in \Omega$. From now on, we underline that this last assumption is not necessary in forthcoming numerical developments.

The system (1) is supposed to be endowed with an entropy inequality. Then we assume the existence of an entropy-pair $(\eta, G)$ where $w \mapsto \eta(w)$ is a convex function and $G: \Omega \to \mathbb{R}$ is the entropy flux function satisfying $\nabla G = \nabla \eta \cdot (\nabla f + A)$, so that the following inequality is satisfied:

$$
\partial_t \eta(w) + \partial_x G(w) \leq \nabla \eta(w) \cdot S(w, x), \quad \text{in a weak sense,} \tag{2}
$$

In many particular cases, one may find a general conservative entropy relation as follows:

$$
\partial_t \hat{\eta}(w, x) + \partial_x \hat{G}(w, x) \le 0, \quad (3)
$$

where $(w, x) \mapsto \hat{\eta}(w, x)$ is a convex function and $\hat{G} : \Omega \times \mathbb{R} \to \mathbb{R}$ is a generalized entropy flux according
to the source term. In the sequel, we shall assume that such general entropy-entropy flux pair exists.

<footnote>*Université de Nantes, CNRS UMR 6629, Laboratoire de Mathématiques Jean Leray, 2 rue de la Houssinière, BP 92208, 44322 Nantes, France.</footnote>

<footnote>†Departamento de Análisis Matemático, Estadística e Investigación Operativa, y Matemática Aplicada, Universidad de Málaga, Campus de Teatinos s/n, Málaga 29080, Spain.</footnote>

<footnote>‡Université Rennes, Inria Rennes and IRMAR UMR CNRS 6625, F-35042 Rennes, France.</footnote>
As long as conservative systems are adopted, namely $A(w) = 0$ for all $w \in \Omega$, the entropy inequalities enforce some stability properties. In fact, they impose the system to be hyperbolic and they rule out non-physical shock discontinuous solutions [32,39,51,52,73]. If $\nabla f + A$ does not reformulate as a Jacobian matrix, the role played by the entropy inequalities turns out to be drastically different. Indeed, no longer the entropy inequalities enforce the system to be hyperbolic. As a consequence, the system (1) may simultaneously admit complex eigenvalues and satisfy (3). Moreover, the non-conservative product $A(w)\partial_x w$ has no distributional sense. Then, weak solutions are now defined according to the viscosity approach [8, 71, 72] or, similarly in a sense to be prescribed, by the path-theory [33, 53, 54].

As a consequence, according to the forthcoming numerical derivations, the interests of the entropy inequalities are twofold. On one hand, they ensure the physical relevance of weak solutions for conservative systems. On the other hand, they produce suitable estimations since from (3) we easily deduce for all $\alpha < \beta$ and all $t > 0$,

$$ \int_{\alpha}^{\beta} \hat{\eta}(w(x, t), x) \, dx \leq \int_{\alpha}^{\beta} \hat{\eta}(w(x, 0), x) \, dx - \int_{0}^{t} (\hat{G}(w(\beta, s), \beta) - \hat{G}(w(\alpha, s), \alpha)) \, ds. $$

Since $\hat{\eta}$ is a convex function, in a sense to be prescribed, we then obtain a weighted-$L^2_{\text{loc}}$ estimation of $w$ for all $t > 0$.

In addition, because of the source term $S$, there exists non trivial stationary solutions for the system (1). In the present work, we are concerned by smooth stationary solutions which, according to (1) and (3), satisfy

$$ \partial_x f(w) + A(w)\partial_x w = S(w, x), \quad (4) $$

$$ \partial_x \hat{G}(w, x) = 0. \quad (5) $$

The shallow-water system clearly enters into such a framework. Let us underline that the well-known Euler system with gravity can also be considered in this paper. We notice that the two-layer shallow-water model [20, 21, 23–25] can be here adopted. Since this last model is only conditionally hyperbolic, as emphasized, the entropy inequalities play a major role to get a suitable estimation of the solutions within elliptic regions.

The derivation of numerical schemes to approximate the weak solutions of (1) has been widely discussed in the literature during the last five decades. We do not propose an exhaustive bibliography, but numerous schemes have been designed for conservative hyperbolic systems (for instance, see [1, 13, 14, 39–41, 46, 50, 78, 79] and references therein). Concerning non-conservative systems, the derivation of numerical schemes turns out to be more delicate. We can cite the pioneer work Hou and LeFloch [48] and essential extensions by Parés [63] where the well-known path-consistent scheme is presented. We may also highlight some scheme derivations based on single shock capturing techniques [29, 68, 69] or based on viscous approach [4, 30].

Concerning the numerical approximation of the source term, we refer to the work by Bermudez and Vazquez [3], Greenberg and LeRoux [45] and Gosse [43] where the well-balanced property is introduced. They established the need to accurately approximate the steady solutions in order to avoir some spurious numerical errors. The well-balanced schemes are then derived (for instance, see [2, 7, 16, 18, 22, 27, 28, 36, 37, 49, 60, 62, 64, 80] and the references therein).

Adopting first-order numerical methods, the discrete entropy inequalities can be reached by several homogeneous conservative schemes [12, 13, 31, 41, 46, 50]. Nevertheless, as soon as source term and non-conservation form are considered, up to our knowledge, very few works are concerned with entropy preserving schemes, and the entropy stability is rarely discussed. However, we may cite [10, 35, 36, 58].

In the present work, according to artificial viscosity technique introduced by Tadmor [74–77] and recent extensions proposed in [5], we present a very easy improvement to make entropy preserving any first-order well-balanced scheme. To address such an issue, the remainder of the paper is organized as follows. In the next section, we present an artificial viscosity technique to recover the required discrete entropy inequalities and to preserve the well-balanced property. In fact, according to [5], we give an easy evaluation of the numerical viscosity coefficient to enforce the expected entropy stability. In order to establish such an entropy preservation, we first present, in Section 3, a relevant approximate
Riemann solver to reformulate any 3-point finite volume scheme endowed with artificial viscosity as a Godunov-type scheme. In the two next sections, equipped with such a reformulation we exhibit the evaluation of the numerical viscosity coefficient such that the stationary solutions are preserved and the discrete entropy inequalities are satisfied.

An important aspect that may go unnoticed is the following: the strategy presented in this work, combined, for example, with the one described in the article [28], allows the extension of any consistent first-order solver for the homogeneous problem—without the presence of a source term—into a solver that is well-balanced (using the strategy described in [28]), robust, and also satisfies the discrete entropy inequality. To the best of our knowledge, this is the first time a simple general strategy is presented such that, given any one-dimensional balance law and a consistent first-order numerical flux for the homogeneous problem, one can construct a solver that is well-balanced, robust, and satisfies the discrete entropy inequality.

Finally, the last section is devoted to several numerical simulations to illustrate the interest of the designed numerical procedure.

## 2 Well-balanced artificial viscosity and main result

We consider a standard 3-point finite volume numerical scheme. For the sake of simplicity in the notation, we shall consider constant $\Delta x > 0$ (resp. $\Delta t > 0$) step of space (resp. time). The space mesh $(x_{i+\frac{1}{2}})_{i \in \mathbb{Z}}$ is defined by $x_{i+\frac{1}{2}} = x_{i-\frac{1}{2}} + \Delta x$ and the quantity $w_i^n$ denotes an approximation of the average $\frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} w(x, t^n) dx$. The sequence $(x_i)_{i \in \mathbb{Z}} \subset \mathbb{R}$ is defined by

$$x_i = \frac{1}{2}(x_{i-\frac{1}{2}} + x_{i+\frac{1}{2}}), \quad \forall i \in \mathbb{Z}. \qquad (6)$$

The numerical scheme will be written in the general form

$$\begin{aligned} w_i^{n+1} = w_i^n &- \frac{\Delta t}{\Delta x} (f_{\Delta}(w_i^n, w_{i+1}^n) - f_{\Delta}(w_{i-1}^n, w_i^n)) \\ &- \frac{\Delta t}{2\Delta x} (A_{\Delta}^L(w_i^n, w_{i+1}^n)(w_{i+1}^n - w_i^n) + A_{\Delta}^R(w_{i-1}^n, w_i^n)(w_i^n - w_{i-1}^n)) \\ &+ \frac{\Delta t}{2} (S_{\Delta}^L(w_i^n, x_i, w_{i+1}^n, x_{i+1}) + S_{\Delta}^R(w_{i-1}^n, x_{i-1}, w_i^n, x_i)), \end{aligned} \qquad (7)$$

where $f_{\Delta}: (\mathbb{R}^d)^2 \to \mathbb{R}^d$ stands for the numerical flux satisfying $f_{\Delta}(w, w) = f(w)$, $A_{\Delta}^{L,R}: (\mathbb{R}^d)^2 \to \mathbb{R}^d$ are consistent matrices with $A$ in the sense of $A_{\Delta}^{L,R}(w, w) = A(w)$ and $S_{\Delta}: (\mathbb{R}^{d+1})^2 \to \mathbb{R}^d$ is a consistent approximation of $S$ namely $S_{\Delta}(w, x, w, x) = S(w, x)$. Let us remark that most of the numerical schemes for (1) can be written under the form of (7) (see for instance [13, 26]).

For the sake of completeness, we give the definition of a well-balanced, robust and entropy-satisfying numerical scheme.

**Definition 2.1** (well-balanced, robust, entropy-satisfying numerical scheme). Let us consider a system (1) endowed by a pair of convex entropy-entropy flux $(\hat{\eta}, \hat{G})$ and the corresponding entropy inequality (3) and let us define

$$\hat{w} = (w, x)^{\mathrm{T}} \in \Omega \times \mathbb{R}. \qquad (8)$$

Let us also consider $(x_i)_{i \in \mathbb{Z}} \subset \mathbb{R}$ defined by (6) and let us denote a relation $\mathcal{R}: (\Omega \times \mathbb{R})^2 \to \mathbb{R}^d$ that defines stationary solutions (4) at a discrete level. A three-point finite volume scheme is

(i) well-balanced for the equilibria $\mathcal{R}$ if

$$\mathcal{R}(\hat{w}_i^n, \hat{w}_{i+1}^n) \text{ holds } \forall i \in \mathbb{Z}, \text{ then } w_i^{n+1} = w_i^n, \forall i \in \mathbb{Z}. \qquad (9)$$

(ii) robust if $(w_i^n)_{i \in \mathbb{Z}} \subset \Omega$ then $(w_i^{n+1})_{i \in \mathbb{Z}} \subset \Omega$.
(iii) entropy-satisfying for the convex entropy-entropy flux ($\hat{\eta}$, $\hat{G}$) if there exists $\hat{G}: (\Omega \times \mathbb{R})^2 \to \mathbb{R}$ such that $\hat{G}(\hat{w}, \hat{w}) = \hat{G}(\hat{w})$ and such that $\hat{w}_i^{n+1} = (w_i^{n+1}, x_i)_{i \in \mathbb{Z}}$ verifies

$$ \frac{\hat{\eta}(\hat{w}_i^{n+1}) - \hat{\eta}(\hat{w}_i^n)}{\Delta t} + \frac{\hat{G}(\hat{w}_i^n, \hat{w}_{i+1}^n) - \hat{G}(\hat{w}_{i-1}^n, \hat{w}_i^n)}{\Delta x} \le 0, \quad \forall i \in \mathbb{Z}. \quad (10) $$

The relation $\mathcal{R}$ is a local definition of a discrete stationary solution. It may be explicit or implicit algebraic equations. Remark that we may distinguish, as in [42], two different scenarios whether the scheme exactly preserves the steady states (exactly well-balanced schemes) or just a discrete approximation of them (well-balanced schemes). In this work, we do not assume specific form for the local discrete equilibrium $\mathcal{R}$, but we assume a non-entropy production if $\mathcal{R}$ holds. This assumption writes

$$ \text{if } \mathcal{R}(\hat{w}_i^n, \hat{w}_{i+1}^n) \text{ holds then } \hat{G}(\hat{w}_{i+1}^n) - \hat{G}(\hat{w}_i^n) = 0. \quad (11) $$

The above assumption is a direct discrete counterpart of (5).

According to Definition 2.1-(i), the first-order well-balanced schemes we may find in the literature (for instance, see [2, 7, 16, 18, 22, 27, 28, 36, 37, 49, 60, 62, 64]), in general recast in the generic form (7). Moreover, these known schemes also satisfy the robustness property according to Definition 2.1-(ii). Now, the main objective is to define an improved version of (7) in order to reach the required entropy stability properties given Definition 2.1-(iii). In this regard, we add an artificial viscosity to scheme (7) to get the following viscous scheme:

$$ \begin{aligned} w_i^{n+1} = w_i^n & - \frac{\Delta t}{\Delta x} (f_\Delta(w_i^n, w_{i+1}^n) - f_\Delta(w_{i-1}^n, w_i^n)) \\ & - \frac{\Delta t}{2\Delta x} (A_\Delta^L(w_i^n, w_{i+1}^n)(w_{i+1}^n - w_i^n) + A_\Delta^R(w_{i-1}^n, w_i^n)(w_i^n - w_{i-1}^n)) \\ & + \frac{\Delta t}{2} (S_\Delta^L(\hat{w}_i^n, \hat{w}_{i+1}^n) + S_\Delta^R(\hat{w}_{i-1}^n, \hat{w}_i^n)), \\ & + \frac{\gamma \Delta t}{2\Delta x} (w_{i+1}^n - 2w_i^n + w_{i-1}^n), \end{aligned} \quad (12) $$

where $\gamma \ge 0$ stands for the artificial viscosity coefficient. Our main result given below establishes a suitable choice for $\gamma$ in order to preserve the well-balanced property and to get both robustness and discrete entropy inequality. In this regard, we assume a particular property on the exact solution of the Riemann problem associated with (1). Indeed, for any two constant states $(w_L, w_R) \in \Omega^2$ and considering the initial condition

$$ w(x, t = 0) = \begin{cases} w_L & \text{if } x < 0, \\ w_R & \text{if } x \ge 0, \end{cases} \quad (13) $$

with

$$ S(w(x, t = 0), x) = \begin{cases} S_L & \text{if } x < 0, \\ S_R & \text{if } x \ge 0, \end{cases} \quad (14) $$

we assume the existence of $\lambda^* > 0$ such that the solution $W_{ex}(x/t, w_L, w_R)$ of the Riemann problem made of (1), (13) is such that

$$ W_{ex}(x/t, w_L, w_R) = \begin{cases} w_L, & \forall x < \lambda^*t, \\ w_R, & \forall x > \lambda^*t. \end{cases} \quad (15) $$

Such a condition may appear restrictive but it is naturally satisfied by both one- and two-layer shallow-water model or by Euler with gravity system. In fact, according to [7, 10, 55, 65], the above property is satisfied as soon as the source term can be incorporated within the hyperbolic operator as follows:

$$ \begin{cases} \partial_t w + \partial_x f(w) + A(w)\partial_x w - S(w, X)\partial_x X = 0, \\ \partial_t X = 0. \end{cases} $$

Our main result is now given.
**Theorem 2.1** (Robustness and entropy stability for well-balanced artificial viscosity schemes). Let us consider (1) endowed by a pair of convex entropy-flux entropy $(\hat{\eta}, \hat{G})$ such that the Hessian matrix $\nabla^2\hat{\eta}(w, x)$ is positive definite. Let $(w_i^n)_{i \in \mathbb{Z}}$ be given in $\Omega$ and consider the following quantities

$$ (\tilde{w}_{i+\frac{1}{2}}^n, \tilde{x}_{i+\frac{1}{2}}) := \frac{1}{2}(w_i^n + w_{i+1}^n, x_i + x_{i+1}) \in \Omega \times \mathbb{R}, \quad (16) $$

$$ w_{i+\frac{1}{2}}^{L, \star} := w_i^n - \frac{f_{\Delta}(w_i^n, w_{i+1}^n) - f(w_i^n)}{\lambda_{i+\frac{1}{2}}} - \frac{1}{2\lambda_{i+\frac{1}{2}}} A_{\Delta}^{L}(w_i^n, w_{i+1}^n)(w_{i+1}^n - w_i^n) \\ + \frac{\Delta x}{2\lambda_{i+\frac{1}{2}}} S_{\Delta}^{L}(\hat{w}_i^n, \hat{w}_{i+1}^n) \in \mathbb{R}^d, \quad (17) $$

$$ w_{i+\frac{1}{2}}^{R, \star} := w_{i+1}^n + \frac{f_{\Delta}(w_i^n, w_{i+1}^n) - f(w_{i+1}^n)}{\lambda_{i+\frac{1}{2}}} - \frac{1}{2\lambda_{i+\frac{1}{2}}} A_{\Delta}^{R}(w_i^n, w_{i+1}^n)(w_{i+1}^n - w_i^n) \\ + \frac{\Delta x}{2\lambda_{i+\frac{1}{2}}} S_{\Delta}^{R}(\hat{w}_i^n, \hat{w}_{i+1}^n) \in \mathbb{R}^d, \quad (18) $$

$$ \mathcal{E}_{i+\frac{1}{2}}^{0} := \lambda_{i+\frac{1}{2}} \left( \hat{\eta}(w_{i+\frac{1}{2}}^{L, \star}, x_i) + \hat{\eta}(w_{i+\frac{1}{2}}^{R, \star}, x_{i+1}) - \hat{\eta}(\hat{w}_i^n) - \hat{\eta}(\hat{w}_{i+1}^n) \right) \quad (19) $$

$$ + (\hat{G}(\hat{w}_{i+1}^n) - \hat{G}(\hat{w}_i^n)) \in \mathbb{R}, \quad (20) $$

$$ \mathcal{D}_{i+\frac{1}{2}} := 2\hat{\eta}(\tilde{w}_{i+\frac{1}{2}}^n, \tilde{x}_{i+\frac{1}{2}}) - \hat{\eta}(\hat{w}_i^n) - \hat{\eta}(\hat{w}_{i+1}^n) \in \mathbb{R}, \quad (21) $$

where $\lambda_{i+\frac{1}{2}} > \lambda^*$. Let a well-balanced scheme (7) defined by $f_{\Delta}$, $A_{\Delta}^{L,R}$ and $S_{\Delta}$ and assume

(h1): the matrices $A_{\Delta}^{L,R}$ are such that

$$ \frac{1}{2} (A_{\Delta}^{L}(w_i^n, w_{i+1}^n) + A_{\Delta}^{R}(w_i^n, w_{i+1}^n))(w_{i+1}^n - w_i^n) = \\ \frac{1}{\Delta t} \int_{0}^{\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} A(W_{ex}(x/t; \hat{w}_i^n, \hat{w}_{i+1}^n)) \partial_x W_{ex}(x/t; \hat{w}_i^n, \hat{w}_{i+1}^n) \, dx \, dt, \quad (22) $$

(h2): the discrete source term $S_{\Delta}^{L,R}$ are such that

$$ \frac{\Delta x}{2} (S_{\Delta}^{L}(\hat{w}_i^n, \hat{w}_{i+1}^n) + S_{\Delta}^{R}(\hat{w}_i^n, \hat{w}_{i+1}^n)) = \\ \frac{1}{\Delta t} \int_{0}^{\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} S(W_{ex}(x/t, \hat{w}_i^n, \hat{w}_{i+1}^n), x + x_{i+\frac{1}{2}}) \, dx \, dt, \quad (23) $$

(h3): the states $w_{i+\frac{1}{2}}^{L, \star}$, $w_{i+\frac{1}{2}}^{R, \star}$ are such that $w_{i+\frac{1}{2}}^{L, \star} = w_i^n$ and $w_{i+\frac{1}{2}}^{R, \star} = w_{i+1}^n$ as soon as the local equilibrium relation $R(\hat{w}_i^n, \hat{w}_{i+1}^n)$ holds.

If $\gamma$ in the viscous scheme (12) is defined as follows:

$$ \gamma = \max_{i \in \mathbb{Z}} \left( 0, -\frac{\mathcal{E}_{i+\frac{1}{2}}^{0}}{\mathcal{D}_{i+\frac{1}{2}}} \right), \quad (24) $$

then $\gamma \ge 0$ is bounded and if the time step $\Delta t > 0$ satisfies

$$ (\lambda_{i+\frac{1}{2}} + \gamma) \frac{\Delta t}{\Delta x} \le \frac{1}{2}, \quad \forall i \in \mathbb{Z}, \quad (25) $$

then the viscous scheme (12) is

(i) still well-balanced,
(ii) robust under a suitable CFL condition: i.e. $(w_i^{n+1})_{i \in \mathbb{Z}} \subset \Omega$,

(iii) entropy-preserving: $(w_i^{n+1})_{i \in \mathbb{Z}}$ satisfies a discrete entropy inequality (10).

The proof of (ii) follows the same arguments as in [5], the main ideas being explained in Lemma 3.1. As a consequence, the establishment of the above main result just consists proving the entropy-preservation given by Theorem 2.1-(iii).

Let us underline that the assumptions $(h_1)$ and $(h_2)$ are related to the discrete entropy stability and the assumption $(h_3)$ concerns to the well-balanced property. Their roles will be highlighted in the sequel that deals with the proof of Theorem 2.1. To address such an issue, the following section gives an equivalent reformulation for the viscous scheme (12).

### 3 Godunov-type scheme reformulation

In this section, the numerical scheme with artificial viscosity (12) is reformulated as a Godunov-type scheme. This reformulation is detailed in the following lemma:

**Lemma 3.1 (Reformulation as Simple Riemann Solver).** Let us consider $(w_{i+\frac{1}{2}}^{L,R,*}, \tilde{w}_{i+\frac{1}{2}}^{n}, \lambda_{i+\frac{1}{2}})$ $\in \Omega \times (\mathbb{R}^d)^2 \times \Omega \times \mathbb{R}_*^+$ respectively defined in (17), (18) and let us define

$$ \tilde{w}_{i+\frac{1}{2}}^{\alpha,*} := \left(1 - \frac{\lambda_{i+\frac{1}{2}}}{\lambda_{i+\frac{1}{2}} + \gamma}\right) \tilde{w}_{i+\frac{1}{2}}^{n} + \frac{\lambda_{i+\frac{1}{2}}}{\lambda_{i+\frac{1}{2}} + \gamma} w_{i+\frac{1}{2}}^{\alpha,*} \in \mathbb{R}^d, \quad \forall \alpha \in \{L, R\}. \quad (26) $$

If the CFL restriction (25) holds, then the viscous scheme (12) is equivalent to the following Godunov-type scheme:

$$ w_i^{n+1} = \frac{1}{\Delta x} \int_0^{\frac{\Delta x}{2}} \tilde{w}(x/\Delta t, \hat{w}_i^n, \hat{w}_i^n) \, dx + \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^0 \tilde{w}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n) \, dx, \quad (27) $$

where $\hat{w} = (w, x)^T \in \Omega \times \mathbb{R}$ and the approximated Riemann solver $\tilde{w}(\cdot, \hat{w}_i^n, \hat{w}_{i+1}^n) : \mathbb{R} \times (\Omega \times \mathbb{R})^2 \to \mathbb{R}^d$ satisfies

$$ \tilde{w}(x/t, \hat{w}_i^n, \hat{w}_{i+1}^n) = \begin{cases} w_i^n & \text{if } x \le -(\lambda_{i+\frac{1}{2}} + \gamma)t, \\ \tilde{w}_{i+\frac{1}{2}}^{L,*} & \text{if } -(\lambda_{i+\frac{1}{2}} + \gamma)t < x \le 0, \\ \tilde{w}_{i+\frac{1}{2}}^{R,*} & \text{if } 0 < x \le (\lambda_{i+\frac{1}{2}} + \gamma)t, \\ w_{i+1} & \text{if } (\lambda_{i+\frac{1}{2}} + \gamma)t < x. \end{cases} \quad (28) $$

Moreover, under a suitable CFL condition, one can prove that $\tilde{w}_{i+\frac{1}{2}}^{\alpha,*} \in \Omega$.

*Proof.* First, let us remark that the last part of the lemma is easily proved following [5]. Indeed, $\Omega$ being convex, there exists $r > 0$ such that for any $w \in \mathbb{R}^d$ satisfying $||w - \tilde{w}_{i+1/2}^n|| < r$, we have $w \in \Omega$. Now, in view of (26), it follows that $\tilde{w}_{i+\frac{1}{2}}^{\alpha,*} \in \Omega$ for $\gamma$ large enough.

Let us now then prove the main result of the lemma. Using the above definition of the approximated Riemann solver $\tilde{w}(\cdot, \hat{w}_i^n, \hat{w}_{i+1}^n) : \mathbb{R} \times (\Omega \times \mathbb{R})^2 \to \mathbb{R}^d$ we have

$$ \begin{aligned} \frac{1}{\Delta x} \int_{0}^{\frac{\Delta x}{2}} \tilde{w}(x/\Delta t, \hat{w}_{i-1}^n, \hat{w}_i^n) \, dx &= (\lambda_{i-\frac{1}{2}} + \gamma) \frac{\Delta t}{\Delta x} \tilde{w}_{i-\frac{1}{2}}^{R,*} + \frac{1}{\Delta x} \left( \frac{\Delta x}{2} - (\lambda_{i-\frac{1}{2}} + \gamma) \right) w_i^n, \\ &= \frac{w_i^n}{2} + (\lambda_{i-\frac{1}{2}} + \gamma) \frac{\Delta t}{\Delta x} \left( \tilde{w}_{i-\frac{1}{2}}^{R,*} - w_i^n \right), \\ &= \frac{w_i^n}{2} + (\lambda_{i-\frac{1}{2}} + \gamma) \frac{\Delta t}{\Delta x} \left( \frac{\gamma}{\lambda_{i-\frac{1}{2}}} (\tilde{w}_{i-\frac{1}{2}} - w_i^n) + \frac{\lambda_{i-\frac{1}{2}}}{\lambda_{i-\frac{1}{2}} + \gamma} (w_{i-\frac{1}{2}}^{R,*} - w_i^n) \right), \end{aligned} $$
$$ = \frac{w_i^n}{2} + \frac{\gamma}{2} \frac{\Delta t}{\Delta x} (-w_i^n + w_{i-1}^n) + \frac{\lambda_{i-\frac{1}{2}} \Delta t}{\Delta x} (w_{i-\frac{1}{2}}^{R,\star} - w_i^n). \quad (29) $$

In the same way, we also have

$$ \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{0} \tilde{w}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n) \, dx = \frac{w_i^n}{2} + \frac{\gamma}{2} \frac{\Delta t}{\Delta x} (w_{i+1} - w_i^n) + \frac{\lambda_{i+\frac{1}{2}} \Delta t}{\Delta x} (w_{i+\frac{1}{2}}^{L,\star} - w_i^n), \quad (30) $$

But, according to (17) and (18), $w_{i+\frac{1}{2}}^{L,\star} - w_{i+1}^n$ and $w_{i-\frac{1}{2}}^{R,\star} - w_i^n$ respectively read

$$ \lambda_{i+\frac{1}{2}}(w_{i+\frac{1}{2}}^{L,\star} - w_i^n) = -f_{\Delta}(w_i^n, w_{i+1}^n) + f(w_i^n) - \frac{1}{2}A_{\Delta}^L(w_i^n, w_{i+1}^n)(w_{i+1}^n - w_i^n) + \frac{\Delta x}{2}S_{\Delta}^L(\hat{w}_i^n, \hat{w}_{i+1}^n), $$

$$ \lambda_{i-\frac{1}{2}}(w_{i-\frac{1}{2}}^{R,\star} - w_i^n) = f_{\Delta}(w_{i-1}^n, w_i^n) - f(w_i^n) - \frac{1}{2}A_{\Delta}^R(w_{i-1}^n, w_i^n)(w_i^n - w_{i-1}^n) + \frac{\Delta x}{2}S_{\Delta}^R(\hat{w}_{i-1}^n, \hat{w}_i^n). $$

As a consequence, using the above formulations in the sum of (29) and (30), we deduce

$$ \frac{1}{\Delta x} \int_{0}^{\frac{\Delta x}{2}} \tilde{w}(x/\Delta t, \hat{w}_{i-1}^n, \hat{w}_i^n) \, dx + \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{0} \tilde{w}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n) \, dx = \\ \frac{w_i^n}{2} + \frac{\Delta t}{\Delta x} \left( f_{\Delta}(w_{i-1}^n, w_i^n) - f(w_i^n) - \frac{1}{2}A_{\Delta}^R(w_{i-1}^n, w_i^n)(w_i^n - w_{i-1}^n) \right) \\ + \frac{\Delta t}{2}S_{\Delta}^R(\hat{w}_{i-1}^n, \hat{w}_i^n) + \frac{\gamma}{2}\frac{\Delta t}{\Delta x}(-w_i^n + w_{i-1}^n) \\ + \frac{w_i^n}{2} + \frac{\Delta t}{\Delta x} \left( -f_{\Delta}(w_i^n, w_{i+1}^n) + f(w_i^n) - \frac{1}{2}A_{\Delta}^L(w_i^n, w_{i+1}^n)(w_{i+1}^n - w_i^n) \right) \\ + \frac{\Delta t}{2}S_{\Delta}^L(\hat{w}_{i-1}^n, \hat{w}_{i+1}^n) + \frac{\gamma}{2}\frac{\Delta t}{\Delta x}(w_{i+1}^n - w_i^n). $$

Since a simple reorganization of the right hand-side of the above equality leads to the expected reformulation that concludes the proof. □

The reformulation (27) is not unique and arguing analogous above computations we could also define the following approximated Riemann solver

$$ \tilde{w}(x/t, \hat{w}_i^n, \hat{w}_{i+1}^n) = \begin{cases} w_i^n & \text{if } x \le -(\lambda_{i+\frac{1}{2}} + \gamma)t, \\ \tilde{w}_{i+\frac{1}{2}} & \text{if } -(\lambda_{i+\frac{1}{2}} + \gamma)t \le x < -\lambda_{i+\frac{1}{2}}t, \\ w_{i+\frac{1}{2}}^{L,\star} & \text{if } -\lambda_{i+\frac{1}{2}}t < x \le 0, \\ w_{i+\frac{1}{2}}^{R,\star} & \text{if } 0 < x \le \lambda_{i+\frac{1}{2}}t, \\ \tilde{w}_{i+\frac{1}{2}} & \text{if } \lambda_{i+\frac{1}{2}}t \le x < (\lambda_{i+\frac{1}{2}} + \gamma)t, \\ w_{i+1}^n & \text{if } (\lambda_{i+\frac{1}{2}} + \gamma)t < x. \end{cases} \quad (31) $$

The Godunov-type scheme reformulation (27) will be very useful to establish a discrete entropy stability. In this regard we give in the following lemma a sufficient condition that ensures a discrete entropy inequality (10). This condition is an extension of the work of Harten et al. [46]. For the sake of completeness, we also recall a sufficient condition for the well-balanced property.
**Lemma 3.2** (Entropy-preserving well-balanced Godunov-type scheme). Let us assume the system (1) is endowed by a pair of convex entropy-entropy flux $(\hat{\eta}, \hat{G})$. Let us consider a Godunov-type scheme (27) that approximates the solutions of (1) and defined by an approximated Riemann solver $\tilde{w}(\cdot, \hat{w}_i^n, \hat{w}_{i+1}^n)$ : $\mathbb{R} \times (\Omega \times \mathbb{R})^2 \to \mathbb{R}^d$. Let us define $\tilde{x} : \mathbb{R}^3 \to \mathbb{R}$ such that

$$x_i = \frac{1}{\Delta x} \int_0^{\frac{\Delta x}{2}} \tilde{x}(x/\Delta t, x_{i-1}, x_i) \, dx + \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^0 \tilde{x}(x/\Delta t, x_i, x_{i+1}) \, dx, \quad (32)$$

where $(x_i)_{i \in \mathbb{Z}} \subset \mathbb{R}$ is defined in (6). Assume the CFL condition (25) holds.

(i) If the couple $(\tilde{w}, \tilde{x})$ satisfies the inequality

$$\frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n)) \, dx \le \frac{\hat{\eta}(\hat{w}_i^n) + \hat{\eta}(\hat{w}_{i+1}^n)}{2} - \frac{\Delta t}{\Delta x} (\hat{G}(\hat{w}_{i+1}^n) - \hat{G}(\hat{w}_i^n)), \quad (33)$$

where $\hat{w} = (w, x)^T \in \Omega \times \mathbb{R}$, then the sequence $(w_i^{n+1}, x_i)_{i \in \mathbb{Z}}$ satisfies a discrete entropy inequality (10) in which the entropy numerical flux $\hat{G}: (\mathbb{R}^d \times \mathbb{R})^2 \to \mathbb{R}$ is given by

$$\begin{aligned} \hat{G}(\hat{w}_i^n, \hat{w}_{i+1}^n) = & \frac{\hat{G}(\hat{w}_{i+1}^n) + \hat{G}(\hat{w}_i^n)}{2} - \frac{\Delta x}{4\Delta t} (\hat{\eta}(\hat{w}_{i+1}^n) - \hat{\eta}(\hat{w}_i^n)) \\ & + \frac{1}{2\Delta t} \int_0^{\frac{\Delta x}{2}} \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n)) \, dx \\ & - \frac{1}{2\Delta t} \int_{-\frac{\Delta x}{2}}^0 \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n)) \, dx. \end{aligned} \quad (34)$$

(ii) If the approximate Riemann solver satisfies

$$\tilde{w}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n) = \begin{cases} w_i^n & \text{if } x < 0, \\ w_{i+1}^n & \text{otherwise,} \end{cases} \quad \forall i \in \mathbb{Z}, \quad (35)$$

as soon as the sequence $(\hat{w}_i^n)_{i \in \mathbb{Z}}$ verifies, at each interface of the mesh, the local relation $\mathcal{R}(\hat{w}_i^n, \hat{w}_{i+1}^n)$ then the Godunov-type scheme is well-balanced for the equilibria $\mathcal{R}$.

*Proof.* Concerning (i), considering (27) and (32) we have

$$(w_i^{n+1}, x_i) = \frac{1}{\Delta x} \int_0^{\frac{\Delta x}{2}} (\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_i^n, \hat{w}_i^n) \, dx + \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^0 (\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n) \, dx.$$

Since the entropy $(w, x) \mapsto \hat{\eta}(w, x)$ is a convex function, the Jensen inequality [67] applied to the above formulation leads to

$$\begin{aligned} \hat{\eta}(w_i^{n+1}, x_i) \le & \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^0 \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n)) \, dx + \frac{1}{2\Delta x} \int_0^{\frac{\Delta x}{2}} \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_{i-1}^n, \hat{w}_i^n)) \, dx \\ & + \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_{i-1}^n, \hat{w}_i^n)) \, dx + \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n)) \, dx \\ & - \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^0 \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_{i-1}^n, \hat{w}_i^n)) \, dx - \frac{1}{2\Delta x} \int_0^{\frac{\Delta x}{2}} \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n)) \, dx. \end{aligned}$$
Using the entropy condition (33) in the above inequality, we obtain

$$
\begin{align*}
\eta(w_i^{n+1}, x_i) \le \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{0} \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n)) \, dx &+ \frac{1}{2\Delta x} \int_{0}^{\frac{\Delta x}{2}} \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_{i-1}^n, \hat{w}_i^n)) \, dx \\
&+ \frac{1}{2} \left( \frac{\hat{\eta}(\hat{w}_{i-1}^n) + \hat{\eta}(\hat{w}_i^n)}{2} - \frac{\Delta t}{\Delta x} (\hat{G}(\hat{w}_i^n) - \hat{G}(\hat{w}_{i-1}^n)) \right) \\
&+ \frac{1}{2} \left( \frac{\hat{\eta}(\hat{w}_i^n) + \hat{\eta}(\hat{w}_{i+1}^n)}{2} - \frac{\Delta t}{\Delta x} (\hat{G}(\hat{w}_{i+1}^n) - \hat{G}(\hat{w}_i^n)) \right) \\
&- \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{0} \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_{i-1}^n, \hat{w}_i^n)) \, dx - \frac{1}{2\Delta x} \int_{0}^{\frac{\Delta x}{2}} \hat{\eta}((\tilde{w}, \tilde{x})(x/\Delta t, \hat{w}_{i-1}^n, \hat{w}_{i+1}^n)) \, dx.
\end{align*}
$$

Reorganized the above terms with the definition of the function $\hat{\mathcal{G}}$ given by (34), we deduce the discrete entropy inequality (10).

For the statement (ii), let us assume that the sequence $(\hat{w}_i^n)_{i \in \mathbb{Z}}$ is such that the local equilibrium $\mathcal{R}(\hat{w}_i^n, \hat{w}_{i+1}^n)$ holds at each interface of the mesh. In this case, using the condition (35) in a Godunov-type scheme (27), we have

$$
\begin{align*}
w_i^{n+1} &= \frac{1}{\Delta x} \int_0^{\frac{\Delta x}{2}} \tilde{w}(x/\Delta t, \hat{w}_{i-1}^n, \hat{w}_i^n) \, dx + \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^0 \tilde{w}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n) \, dx, \\
&= \frac{1}{\Delta x} \int_0^{\frac{\Delta x}{2}} w_i^n \, dx + \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^0 w_i^n \, dx, \\
&= w_i^n.
\end{align*}
$$

According to Definition 2.1-(i) this last equality gives the required well-balanced property that con- cludes the proof. $\square$

On the one hand, remark that the relation (35) formulates a sufficient condition for the well-balanced property. This condition is well-known (see for instance [6,9,59,61]) and direct computations could show the non-corrected scheme (7) can also be written as in (28) for the case $\gamma = 0$, so that the assumption $(h_3)$ is equivalent to (35).

On the other hand, the entropy condition (33) needs the definition of $\tilde{x} : (\mathbb{R})^3 \to \mathbb{R}$ and, in this work, we consider the following choice:

**Lemma 3.3.** Let us consider $(x_i)_{i \in \mathbb{Z}} \subset \mathbb{R}$, $\tilde{x}_{i+\frac{1}{2}} \in \mathbb{R}$ defined in (6), (16) and $\tilde{x} : (\mathbb{R})^3 \to \mathbb{R}$ such that

$$
\tilde{x}(x/t, x_i, x_{i+1}) = \begin{cases}
x_i & \text{if } x \le -(\lambda_{i+\frac{1}{2}} + \gamma)t, \\
\tilde{x}_{i+\frac{1}{2}} & \text{if } -(\lambda_{i+\frac{1}{2}} + \gamma)t \le x < -\lambda_{i+\frac{1}{2}}t, \\
x_i & \text{if } -\lambda_{i+\frac{1}{2}}t < x \le 0, \\
x_{i+1} & \text{if } 0 < x \le \lambda_{i+\frac{1}{2}}t, \\
\tilde{x}_{i+\frac{1}{2}} & \text{if } \lambda_{i+\frac{1}{2}}t \le x < (\lambda_{i+\frac{1}{2}} + \gamma)t, \\
x_{i+1} & \text{if } (\lambda_{i+\frac{1}{2}} + \gamma)t < x.
\end{cases} \tag{36}
$$

The above definition for $\tilde{x}(x/t, x_i, x_{i+1})$ satisfies (32).

*Proof.* According to (36), straightforward computations yield to

$$
\frac{1}{\Delta x} \int_{0}^{\frac{\Delta x}{2}} \tilde{x}(x/\Delta t, x_{i-1}, x_i) \, dx + \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{0} \tilde{x}(x/\Delta t, x_i, x_{i+1}) \, dx = x_i + \frac{\gamma \Delta t}{2 \Delta x} (x_{i+1} - 2x_i + x_{i-1}).
$$
But, arguing the definition of $x_i$ given by (6), we have

$$x_{i+1} - 2x_i + x_{i-1} = (\Delta x + x_i - 2x_i + x_i - \Delta x) = 0,$$

that concludes the proof.

Equipped with Lemma 3.2, 3.3, we are now able to show Theorem 2.1. The following section concerns the proof of the preservation of the well-balanced property satisfied by the viscous scheme (12).

# 4 Well-balanced artificial viscosity

In this section, we establish Theorem 2.1-(i) which states that the viscous correction (12) does not spoil the well-balanced property satisfied by the scheme (7). According to Lemma 3.2-(ii), the proof will be completed as soon as (35) is shown for the approximate Riemann solver (28). The proof is now detailed.

*Proof of Theorem 2.1-(i).* Let us assume $(\hat{w}_i^n)_{i \in \mathbb{Z}}$ are such that the local relation $\mathcal{R}(\hat{w}_i^n, \hat{w}_{i+1}^n)$ holds at each interface of the space mesh. As a consequence, according to the assumption $(h_3)$, the equalities $w_{i+\frac{1}{2}}^{L,\star} = w_i^n$, $w_{i+\frac{1}{2}}^{R,\star} = w_{i+1}^n$ hold and $\mathcal{E}_{i+\frac{1}{2}}^0$ defined in (20) reads

$$ \begin{aligned} \mathcal{E}_{i+\frac{1}{2}}^0 &= \lambda_{i+\frac{1}{2}} \left( \hat{\eta}(w_{i+\frac{1}{2}}^{L,\star}, x_i) + \hat{\eta}(w_{i+\frac{1}{2}}^{R,\star}, x_{i+1}) - \hat{\eta}(\hat{w}_i^n) - \hat{\eta}(\hat{w}_{i+1}^n) \right) + (\hat{G}(\hat{w}_{i+1}^n) - \hat{G}(\hat{w}_i^n)), \\ &= \hat{G}(\hat{w}_{i+1}^n) - \hat{G}(\hat{w}_i^n). \end{aligned} $$

Arguing the non-entropy production (11), the above formulation simply reduces to $\mathcal{E}_{i+\frac{1}{2}}^0 = 0$. Hence, the artificial viscosity $\gamma$ defined in (24) reduces to $\gamma = 0$ when the equilibrium relations $\mathcal{R}$ holds. In this case, the approximate Riemann solver (28) simply reads

$$ \tilde{w}(x/t, \hat{w}_i^n, \hat{w}_{i+1}^n) = \begin{cases} w_i^n & \text{if } x \le -\lambda_{i+\frac{1}{2}}t, \\ w_{i+\frac{1}{2}}^{L,\star} & \text{if } -\lambda_{i+\frac{1}{2}}t < x \le 0, \\ w_{i+\frac{1}{2}}^{R,\star} & \text{if } 0 < x \le \lambda_{i+\frac{1}{2}}t, \\ w_{i+1} & \text{if } \lambda_{i+\frac{1}{2}}t < x, \end{cases} $$

where the states $w_{i+\frac{1}{2}}^{L,R,\star}$ are defined in (17), (18). Using once again the assumption $(h_3)$, we deduce the above approximate Riemann solver verifies

$$ \tilde{w}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n) = \begin{cases} w_i^n & \text{if } x < 0, \\ w_{i+1}^n & \text{otherwise,} \end{cases} \quad \forall i \in \mathbb{Z}. $$

According to Lemma 3.2 this last equality achieves the proof.

# 5 Discrete entropy stability

The aim of this section is to show that $\gamma$ given by (24) in the viscous scheme (12) (or in its equivalent reformulation (27), (31)) ensures (33) and thus the discrete entropy inequality (10) holds.
According (31) and (36), the viscous scheme (12) may be equivalently defined by the following couple $(\tilde{w}, \tilde{x}) : \mathbb{R} \times (\Omega \times \mathbb{R})^2 \to \mathbb{R}^{d+1}$:

$$ (\tilde{w}, \tilde{x})(x/t, \hat{w}_i^n, \hat{w}_{i+1}^n) = \begin{cases} (w_i^n, x_i) & \text{if } x \le -(\lambda_{i+\frac{1}{2}} + \gamma)t, \\ (\tilde{w}_{i+\frac{1}{2}}, \tilde{x}_{i+\frac{1}{2}}) & \text{if } -(\lambda_{i+\frac{1}{2}} + \gamma)t \le x < -\lambda_{i+\frac{1}{2}}t, \\ (w_{i+\frac{1}{2}}^{L,\star}, x_i) & \text{if } -\lambda_{i+\frac{1}{2}}t < x \le 0, \\ (w_{i+\frac{1}{2}}^{R,\star}, x_{i+1}) & \text{if } 0 < x \le \lambda_{i+\frac{1}{2}}t, \\ (\tilde{w}_{i+\frac{1}{2}}, \tilde{x}_{i+\frac{1}{2}}) & \text{if } \lambda_{i+\frac{1}{2}}t \le x < (\lambda_{i+\frac{1}{2}} + \gamma)t, \\ (w_{i+1}^n, x_{i+1}) & \text{if } (\lambda_{i+\frac{1}{2}} + \gamma)t < x, \end{cases} $$

where $(\tilde{w}_{i+\frac{1}{2}}, \tilde{x}_{i+\frac{1}{2}}, w_{i+\frac{1}{2}}^{L,R,\star}) \in \Omega \times \mathbb{R} \times (\mathbb{R}^d)^2$ are given by (16), (17), (18). As a consequence, a straightforward computation shows that (33) for the viscous scheme (12) reads

$$ \mathcal{E}_{i+\frac{1}{2}}^{0} + \gamma \mathcal{D}_{i+\frac{1}{2}} \leq 0, $$

where $\mathcal{E}_{i+\frac{1}{2}}^0$ and $\mathcal{D}_{i+\frac{1}{2}}$ are given by (20), (21) respectively. According to the above inequality, the proof of the entropy-preservation for the viscous scheme (12) stated in Theorem 2.1-(iii) will be proved as soon as an upper-bound for the ratio $-\mathcal{E}_{i+\frac{1}{2}}^0/\mathcal{D}_{i+\frac{1}{2}}$ will be established when $\mathcal{E}_{i+\frac{1}{2}}^0 \ge 0$. In this regard, we first give below a lemma that concerns an additional useful inequality.

**Lemma 5.1.** Let $(w_i^n)_{i \in \mathbb{Z}}$ be given in $\Omega$ and $\lambda_{i+\frac{1}{2}} > \lambda^*$ where $\lambda^*$ is such that (15) holds. Assume that the matrices $A_{\Delta}^{L,R}$ that approximate the non-conservative product $A(w)\partial_x w$ satisfy (22) and the discrete source terms $S_{\Delta}^{L,R}$ are such that (23) holds. If $\Delta t > 0$ satisfies $\frac{\lambda_{i+\frac{1}{2}}\Delta t}{\Delta x} \le \frac{1}{2}$ then the following inequality holds

$$ \hat{\eta}(\hat{w}_{i+\frac{1}{2}}^{\text{HLL}}, \tilde{x}_{i+\frac{1}{2}}) \le \hat{\eta}_{i+\frac{1}{2}}^{\text{HLL}}, \quad (37) $$

where $\tilde{x}_{i+\frac{1}{2}} \in \mathbb{R}$ is given by (16) and where we have set

$$ \begin{aligned} \hat{w}_{i+\frac{1}{2}}^{\text{HLL}} := & \frac{1}{2}(w_i^n + w_{i+1}^n) - \frac{1}{2\lambda_{i+\frac{1}{2}}} (f(w_{i+1}^n) - f(w_i^n)) \\ & - \frac{1}{4\lambda_{i+\frac{1}{2}}} \left( A_{\Delta}^{L}(w_i^n, w_{i+1}^n) + A_{\Delta}^{R}(w_i^n, w_{i+1}^n) \right) (w_{i+1}^n - w_i^n) \\ & + \frac{\Delta x}{4\lambda_{i+\frac{1}{2}}} \left( S_{\Delta}^{L}(\hat{w}_i^n, \hat{w}_{i+1}^n) + S_{\Delta}^{R}(\hat{w}_i^n, \hat{w}_{i+1}^n) \right) \in \mathbb{R}^d, \end{aligned} \quad (38) $$

$$ \hat{\eta}_{i+\frac{1}{2}}^{\text{HLL}} := \frac{\hat{\eta}(\hat{w}_i^n) + \hat{\eta}(\hat{w}_{i+1}^n)}{2} - \frac{\hat{G}(\hat{w}_{i+1}^n) - \hat{G}(\hat{w}_i^n)}{2\lambda_{i+\frac{1}{2}}} \in \mathbb{R}. \quad (39) $$

*Proof.* Since the CFL condition $\lambda_{i+\frac{1}{2}}\Delta t/\Delta x \le 1/2$ holds, the exact solution of juxtaposed Riemann's problems $(W_{ex}((x-x_{i+\frac{1}{2}})/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n))_{i \in \mathbb{Z}}$ satisfies

$$ \partial_t W_{ex} + \partial_x f(W_{ex}) + A(W_{ex})\partial_x W_{ex} = S(W_{ex}, x), \quad (40a) $$

$$ \partial_t \hat{\eta}(W_{ex}, x) + \partial_x \hat{G}(W_{ex}, x) \le 0. \quad (40b) $$

As a consequence, applying the operator $\int_0^{\Delta t} \left( \frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_{x_{i+\frac{1}{2}}-\lambda_{i+\frac{1}{2}}\Delta t}^{x_{i+\frac{1}{2}}+\lambda_{i+\frac{1}{2}}\Delta t} \cdot dx \right) dt$ to (40) with $\lambda_{i+\frac{1}{2}}$ large
enough according to (15), we deduce

$$
\begin{aligned}
& \frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} \mathcal{W}_{\text{ex}}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n) \, dx = \\
& \frac{1}{2}(w_i^n + w_{i+1}^n) - \frac{1}{2\lambda_{i+\frac{1}{2}}} (f(w_{i+1}^n) - f(w_i^n)) \\
& \quad - \frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_0^{\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} A(\mathcal{W}_{\text{ex}}(x/t, \hat{w}_i^n, \hat{w}_{i+1}^n)) \partial_x \mathcal{W}_{\text{ex}}(x/t, \hat{w}_i^n, \hat{w}_{i+1}^n) \, dx \, dt \\
& \quad + \frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_0^{\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} S(\mathcal{W}_{\text{ex}}(x/t, \hat{w}_i^n, \hat{w}_{i+1}^n), x + x_{i+\frac{1}{2}}) \, dx \, dt, \\
& \frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} \hat{\eta} \left( \mathcal{W}_{\text{ex}}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n), x + x_{i+\frac{1}{2}} \right) \, dx \le \\
& \quad \frac{1}{2} (\hat{\eta}(\hat{w}_i^n) + \hat{\eta}(\hat{w}_{i+1}^n)) - \frac{1}{2\lambda_{i+\frac{1}{2}}} (\hat{G}(\hat{w}_{i+1}^n) - \hat{G}(\hat{w}_i^n)).
\end{aligned}
$$

Using the assumptions (22), (23), and the definition of $\hat{\eta}_{i+\frac{1}{2}}^{\text{HLL}}$ given by (39), the above formulations write

$$
\begin{aligned}
& \frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} \mathcal{W}_{\text{ex}}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n) \, dx = \\
& \frac{1}{2}(w_i^n + w_{i+1}^n) - \frac{1}{2\lambda_{i+\frac{1}{2}}} (f(w_{i+1}^n) - f(w_i^n)) - \frac{1}{4\lambda_{i+\frac{1}{2}}} \left( A_{\Delta}^{L}(w_i^n, w_{i+1}^n) + A_{\Delta}^{R}(w_i^n, w_{i+1}^n) \right) (w_{i+1}^n - w_i^n) \\
& \quad + \frac{\Delta x}{4\lambda_{i+\frac{1}{2}}} \left( S_{\Delta}^{L}(\hat{w}_i^n, \hat{w}_{i+1}^n) + S_{\Delta}^{R}(\hat{w}_i^n, \hat{w}_{i+1}^n) \right), \\
& \frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} \hat{\eta} \left( \mathcal{W}_{\text{ex}}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n), x + x_{i+\frac{1}{2}} \right) \, dx \le \hat{\eta}_{i+\frac{1}{2}}^{\text{HLL}}.
\end{aligned}
$$

Thanks to the definition of $\hat{w}_{i+\frac{1}{2}}^{\text{HLL}}$ detailed in (38) the above relations simply reduce to

$$
\begin{aligned}
\hat{w}_{i+\frac{1}{2}}^{\text{HLL}} &= \frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} \mathcal{W}_{\text{ex}}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n) \, dx, \\
\frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} \hat{\eta} \left( \mathcal{W}_{\text{ex}}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n), x + x_{i+\frac{1}{2}} \right) \, dx \le \hat{\eta}_{i+\frac{1}{2}}^{\text{HLL}}.
\end{aligned}
$$

But, according to (6), (16), we have

$$
\tilde{x}_{i+\frac{1}{2}} = \frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} (x + x_{i+\frac{1}{2}}) \, dx.
$$

As a consequence, using the Jensen inequality [67], we deduce

$$
\begin{aligned}
\hat{\eta}(\hat{w}_{i+\frac{1}{2}}^{\text{HLL}}, \tilde{x}_{i+\frac{1}{2}}) &= \hat{\eta} \left( \frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} (\mathcal{W}_{\text{ex}}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n), x + x_{i+\frac{1}{2}}) \, dx \right), \\
&\le \frac{1}{2\lambda_{i+\frac{1}{2}}\Delta t} \int_{-\lambda_{i+\frac{1}{2}}\Delta t}^{\lambda_{i+\frac{1}{2}}\Delta t} \hat{\eta}(\mathcal{W}_{\text{ex}}(x/\Delta t, \hat{w}_i^n, \hat{w}_{i+1}^n), x + x_{i+\frac{1}{2}}) \, dx, \\
&\le \hat{\eta}_{i+\frac{1}{2}}^{\text{HLL}},
\end{aligned}
$$

which achieves the proof.

<figure><img src="image_2.png" /></figure>
The state $\hat{w}_{i+\frac{1}{2}}^{\text{HLL}}$ defined in (38) is an extension of the state derived in [46]. This state and the inequality (37) play an essential role in the proof of the entropy-preserving property verified by the viscous scheme (12) (or equivalently (27)). The conclusion of the proof is now given.

*Proof of Theorem 2.1-(iii).* According to Lemma 3.2-(i), it is sufficient to establish the inequality

$$ \mathcal{E}_{i+\frac{1}{2}}^{0} + \gamma \mathcal{D}_{i+\frac{1}{2}} \leq 0, \quad \forall i \in \mathbb{Z}, \qquad (41) $$

with $\mathcal{E}_{i+\frac{1}{2}}^{0}$ and $\mathcal{D}_{i+\frac{1}{2}}$ defined in (20), (21). Since $\mathcal{D}_{i+\frac{1}{2}} \leq 0$, it is enough to show the ratio $-\mathcal{E}_{i+\frac{1}{2}}^{0}/\mathcal{D}_{i+\frac{1}{2}}$ is bounded when $\mathcal{E}_{i+\frac{1}{2}}^{0} \geq 0$. As a consequence, the proof will be completed as soon as an upper-bound for this ratio is be given.

From Lemma 5.1, according to the definition of $\mathcal{E}_{i+\frac{1}{2}}^{0}$, $\mathcal{D}_{i+\frac{1}{2}}$ and $\hat{\eta}_{i+\frac{1}{2}}^{\text{HLL}}$ given by (20), (21), (39) we have

$$ \begin{aligned} -\frac{\mathcal{E}_{i+\frac{1}{2}}^{0}}{\mathcal{D}_{i+\frac{1}{2}}} &= \lambda_{i+\frac{1}{2}} \frac{\hat{\eta}(w_{i+\frac{1}{2}}^{L,*}, x_i) + \hat{\eta}(w_{i+\frac{1}{2}}^{R,*}, x_{i+1}) - 2\hat{\eta}_{i+\frac{1}{2}}^{\text{HLL}}}{\hat{\eta}(\hat{w}_i^n) + \hat{\eta}(\hat{w}_{i+1}^n) - 2\hat{\eta}(\tilde{w}_{i+\frac{1}{2}}, \tilde{x}_{i+\frac{1}{2}})}, \\ &\leq \lambda_{i+\frac{1}{2}} \frac{\hat{\eta}(w_{i+\frac{1}{2}}^{L,*}, x_i) + \hat{\eta}(w_{i+\frac{1}{2}}^{R,*}, x_{i+1}) - 2\hat{\eta}(\hat{w}_{i+\frac{1}{2}}^{\text{HLL}}, \tilde{x}_{i+\frac{1}{2}})}{\hat{\eta}(\hat{w}_i^n) + \hat{\eta}(\hat{w}_{i+1}^n) - 2\hat{\eta}(\tilde{w}_{i+\frac{1}{2}}, \tilde{x}_{i+\frac{1}{2}})}, \end{aligned} $$

where $\hat{w}_{i+\frac{1}{2}}^{\text{HLL}}$ is given by (38). Noticing $(\hat{w}_{i+\frac{1}{2}}^{\text{HLL}}, \tilde{x}_{i+\frac{1}{2}}) = (\frac{1}{2}(w_{i+\frac{1}{2}}^{L,*} + w_{i+\frac{1}{2}}^{R,*}), \frac{1}{2}(x_i + x_{i+1}))$, the above inequality reads

$$ -\frac{\mathcal{E}_{i+\frac{1}{2}}^{0}}{\mathcal{D}_{i+\frac{1}{2}}} \leq \lambda_{i+\frac{1}{2}} \frac{\hat{\eta}(w_{i+\frac{1}{2}}^{L,*}, x_i) + \hat{\eta}(w_{i+\frac{1}{2}}^{R,*}, x_{i+1}) - 2\hat{\eta}(\frac{1}{2}((w_{i+\frac{1}{2}}^{L,*}, x_i) + (w_{i+\frac{1}{2}}^{R,*}, x_{i+1})))}{\hat{\eta}(\hat{w}_i^n) + \hat{\eta}(\hat{w}_{i+1}^n) - 2\hat{\eta}(\tilde{w}_{i+\frac{1}{2}}, \tilde{x}_{i+\frac{1}{2}})}. \qquad (42) $$

However, since $\hat{\eta}$ is a convex function with $\nabla^2 \hat{\eta}$ a definite positive matrix, $\hat{\eta}(\hat{w}_i^n) + \hat{\eta}(\hat{w}_{i+1}^n) - 2\hat{\eta}(\tilde{w}_{i+\frac{1}{2}}, \tilde{x}_{i+\frac{1}{2}}) = 0$ if and only if $\hat{w}_{i+1}^n = \hat{w}_i^n$. Hence, arguing standard, but laborious, Taylor expansion in the neighborhood of $\hat{w}_i^n = \hat{w}_{i+1}^n$ of the right hand side of (42) ensures the expected boundedness as follows:

$$ \frac{\hat{\eta}(w_{i+\frac{1}{2}}^{L,*}, x_i) + \hat{\eta}(w_{i+\frac{1}{2}}^{R,*}, x_{i+1}) - 2\hat{\eta}(\frac{1}{2}((w_{i+\frac{1}{2}}^{L,*}, x_i) + (w_{i+\frac{1}{2}}^{R,*}, x_{i+1})))}{\hat{\eta}(\hat{w}_i^n) + \hat{\eta}(\hat{w}_{i+1}^n) - 2\hat{\eta}(\tilde{w}_{i+\frac{1}{2}}, \tilde{x}_{i+\frac{1}{2}})} \sim \mathcal{O}(1), $$

that concludes the proof. $\square$

# 6 Numerical results

In this section we apply the technique described previously for some particular cases. Specifically, we consider first the shallow water equations as an example of a system of conservation laws with source term. Second, we address the case of the two-layer shallow water system as an example that includes non-conservative products.

## 6.1 Single layer Shallow Water equations

In this section, we consider the single layer shallow water equations given by

$$ \partial_t \begin{pmatrix} h \\ hu \end{pmatrix} + \partial_x \begin{pmatrix} hu \\ hu^2 + gh^2/2 \end{pmatrix} = \begin{pmatrix} 0 \\ -gh \end{pmatrix} \partial_x z, \quad \forall (x, t) \in I \times (0, T], \qquad (43) $$

where $I \subset \mathbb{R}$ is a given interval and $T > 0$ is the final time. This model governs the water height $h \geq 0$ and the horizontal velocity $u \in \mathbb{R}$ of a fluid over a given time independent topography function $z : \mathbb{R} \to \mathbb{R}$. The gravitational is denoted by $g > 0$. The unknown state vector $w = (h, hu)^{\text{T}} \in \mathbb{R}^2$ is assumed to be in the convex set $\Omega = \{(h, hu) \in \mathbb{R}^2 | h \geq 0, hu \in \mathbb{R}\}$. The system (43) is endowed with a given initial data $w_0 : \mathbb{R} \to \Omega$ that prescribes $w(x, t = 0)$.
It is well-known [13,37] that the solutions of (43) satisfy the following conservative entropy inequality

$$ \partial_t (\eta(w) + ghz) + \partial_x (G(w) + ghzu) \le 0, \quad (44) $$

where $\eta(w) = hu^2/2+gh^2/2$, $G(w) = hu^3/2+gh^2u$. However, for a given bottom topography $x \mapsto z(x)$, the entropy $(w, x) \mapsto \eta(w) + ghz(x)$ is not necessarily convex. In order to recover this convexity property, we assume the solutions of (43) are bounded i.e. $w \in L^\infty$ and also $(\partial_x z, \partial_x^2 z) \in (L^+)^2$ then for $C > 0$, we add $\partial_t(gCx^2) = 0$ to $\eta(w) + ghz$. As a consequence, the inequality (44) is also equivalent to

$$ \partial_t \hat{\eta}_C(w, x) + \partial_x \hat{G}(w, x) \le 0, \quad (45) $$

with the following entropy-entropy flux pair functions:

$$ \hat{\eta}_C(w, x) = hu^2/2 + gh^2/2 + ghz(x) + gCx^2, \quad (46a) $$

$$ \hat{G}(w, x) = hu^3/2 + gh^2u + ghuz(x). \quad (46b) $$

Since $(w, \partial_x z, \partial_x^2 z) \in (L^\infty)^3$, a straightforward computation shows the existence of $C > 0$ large enough such that $(w, x) \mapsto \hat{\eta}_C(w, x)$ in (46a) is strictly convex.

On the other hand, the steady states of (43) are given by

$$ hu = \text{uste}, \quad \text{and} \quad u^2/2 + g(h + z) = \text{uste}. \quad (47) $$

Using the above equilibrium formulation, direct computations show that $\hat{G}$ defined in (46b) satisfies the non-entropy production assumption (11). In addition, among equilibria (47), a special attention is often devoted (see for instance [11, 38, 47, 56, 57, 70, 81]) to the lake-at-rest

$$ u = 0, \quad \text{and} \quad h + z = \text{uste}. \quad (48) $$

It is usual to say that a numerical scheme is fully well-balanced if it preserves all the steady states, while it is well-balanced if it preserves lake-at-rest solutions. Let us remark that if the base scheme is either well-balanced or fully well-balanced, the artificial viscosity added as described previously will not change this property.

We will consider different scenarios comparing several well-known numerical schemes. In particular, we consider the schemes and corresponding notations given in Table 1. In what follows, the corresponding viscous versions (12) for the schemes in Table 1 will be denoted with the additional suffix +AV. In practice, the artificial viscosity $\gamma$ defined in (24) is computed according to Algorithm 1. The regularization parameter $\epsilon > 0$ in Algorithm 1 is fixed to $\epsilon = 10^{-14}$ and we denote $\lambda = \max_{i \in \mathbb{Z}} \lambda_{i+\frac{1}{2}}$. An admissible choice for $C > 0$ in (46a) is computed from the exact solution (or a reference solution in a very fine mesh) for each test case.
<table><caption>Table 1: Numerical schemes used for the numerical experiments devote to the single layer shallow water model (43).</caption><tr><td>FullWB:</td><td>Entropic well-balanced scheme defined in [58, Section 6]. This scheme satisfies a fully discrete version of (44) and it is well-balanced for the moving equilibria (47). The well-balanced property may be locally lost depending on the practical implementation of the scheme (see [58, Section 7.3] for the details).</td></tr><tr><td>HR:</td><td>Hydrostatic-reconstruction scheme introduced in [2]. The HR allows, from a given base numerical flux for the case of flat bottom, to define a numerical scheme for shallow water system with topography that is well-balanced. If the base scheme is positive and/or entropy-satisfying, then the resulting scheme is positive preserving and/or it verifies a semi-discrete version of (44). As base scheme we select here the VF-ROE numerical flux [17].</td></tr><tr><td>GenHR:</td><td>Generalized hydrostatic reconstruction introduced in [27]. This technique allows to build a fully well-balanced scheme from a given base scheme for the conservative case. Here VF-ROE numerical flux [17] is selected.</td></tr><tr><td>SubBM:</td><td>Subsonic well-balanced scheme introduced in [15]. This reconstruction technique allows, from a given base scheme for the conservative case, to define a numerical scheme that preserves all subsonic equilibria. If the base scheme is positive and/or entropy-satisfying, then the resulting scheme is positive preserving and/or it satisfies a semi-discrete version of (44). Here we use as base scheme VF-ROE numerical flux [17].</td></tr></table>

<table><tr><td>Algorithm 1 Well-balanced artificial viscosity for the single layer Shallow Water model</td><td></td></tr><tr><td>1: procedure</td><td></td></tr><tr><td>2: Data:</td><td></td></tr><tr><td>3: ε &gt; 0;</td><td>▷ Regularization parameter</td></tr><tr><td>4: (w<sup>n</sup><sub>i</sub>)<sub>i∈Z</sub>, (λ<sub>i+1/2</sub>)<sub>i∈Z</sub>, (η<sub>C</sub>, Ĝ), non-corrected scheme (7);</td><td></td></tr><tr><td>5: Initialization:</td><td></td></tr><tr><td>6: γ &lt; 0;</td><td></td></tr><tr><td>7: Main Loop:</td><td></td></tr><tr><td>8: for i + 1/2 do</td><td>▷ Loop on interface (x<sub>i+1/2</sub>)<sub>i∈Z</sub></td></tr><tr><td>9: E<sup>0</sup><sub>i+1/2</sub> &lt; λ<sub>i+1/2</sub> ( (η<sub>C</sub>(w<sup>L,*</sup><sub>i+1/2</sub>,x<sub>i</sub>)-η<sub>C</sub>(ŵ<sup>n</sup><sub>i</sub>)) / Δx + (η<sub>C</sub>(w<sup>R,*</sup><sub>i+1/2</sub>,x<sub>i+1</sub>)-η<sub>C</sub>(ŵ<sup>n</sup><sub>i+1</sub>)) / Δx ) + (Ĝ(ŵ<sup>n</sup><sub>i+1</sub>)-Ĝ(ŵ<sup>n</sup><sub>i</sub>)) / Δx;</td><td></td></tr><tr><td>10: if E<sup>0</sup><sub>i+1/2</sub> &gt; ε then</td><td></td></tr><tr><td>11: D<sub>i+1/2</sub> &lt; (η<sub>C</sub>(ŵ<sub>i+1/2</sub>,x̂<sub>i+1/2</sub>)-η<sub>C</sub>(ŵ<sup>n</sup><sub>i</sub>)) / Δx + (η<sub>C</sub>(ŵ<sub>i+1/2</sub>,x̂<sub>i+1/2</sub>)-η<sub>C</sub>(ŵ<sup>n</sup><sub>i+1</sub>)) / Δx;</td><td></td></tr><tr><td>12: γ<sub>i+1/2</sub> &lt; -√2 (E<sup>0</sup><sub>i+1/2</sub> D<sub>i+1/2</sub> / √D<sup>4</sup><sub>i+1/2</sub> + max(ε, D<sup>4</sup><sub>i+1/2</sub>))</td><td>(49)</td></tr><tr><td>13: γ &lt; max(γ<sub>i+1/2</sub>, γ);</td><td></td></tr><tr><td>14: end if</td><td></td></tr><tr><td>15: end for</td><td></td></tr><tr><td>16: Result:</td><td></td></tr><tr><td>return γ;</td><td></td></tr><tr><td>17: end procedure</td><td></td></tr></table>
### 6.1.1 Lake-at-rest steady state

Our first numerical experiment is devoted to check the well-balanced property for the lake-at-rest (48). The initial condition and the topography bottom are the following:

$$h_0(x) = \max(2 - z(x), 0), \quad u_0(x) = 0, \quad z(x) = \max(0, 1 - |4x - 2|). \qquad (50)$$

The domain $[0, 1]$ is discretized with 400 cells. The unknown $h$ and $u$ follow periodic boundary conditions on both space boundaries. The constant $C > 0$ in (46a) verifies $C = 8$ and the final time is 1.0. Figure 1 and Table 2 display the compared results.

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Figure 1: Numerical results and time evolution of the viscosities for the equilibrium at rest (50) on a mesh composed of 400 cells. The legend of the schemes is given in Table 1, +AV denotes the viscous version (12).</figcaption></figure>

<table><caption>Numerical errors on $(h, hu)$ for the lake at rest (50).</caption><thead><tr><th></th><th colspan="3">h</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>FullWB+AV</th><td>0.0E-00</td><td>0.0E-00</td><td>0.0E-00</td></tr><tr><th>FullWB</th><td>0.0E-00</td><td>0.0E-00</td><td>0.0E-00</td></tr></tbody></table>

<table><caption>Numerical errors on $(h, hu)$ for the lake at rest (50).</caption><thead><tr><th></th><th colspan="3">hu</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>FullWB+AV</th><td>1.16E-19</td><td>1.86E-19</td><td>2.56E-21</td></tr><tr><th>FullWB</th><td>1.16E-19</td><td>1.86E-19</td><td>2.56E-21</td></tr></tbody></table>

<table><thead><tr><th></th><th colspan="3">h</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>SubBM+AV</th><td>3.01E-16</td><td>3.34E-16</td><td>2.85E-16</td></tr><tr><th>SubBM</th><td>5.41E-18</td><td>6.49E-18</td><td>3.03E-18</td></tr></tbody></table>

<table><caption>Numerical errors on $(h, hu)$ for the lake at rest (50).</caption><thead><tr><th></th><th colspan="3">hu</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>SubBM+AV</th><td>1.83E-15</td><td>2.16E-15</td><td>1.07E-17</td></tr><tr><th>SubBM</th><td>1.24E-17</td><td>1.48E-17</td><td>7.75E-18</td></tr></tbody></table>

<table><thead><tr><th></th><th colspan="3">h</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>GenHR+AV</th><td>7.85E-16</td><td>8.73E-16</td><td>1.08E-15</td></tr><tr><th>GenHR</th><td>1.27E-18</td><td>1.38E-18</td><td>1.30E-18</td></tr></tbody></table>

<table><caption>Numerical errors on $(h, hu)$ for the lake at rest (50).</caption><thead><tr><th></th><th colspan="3">hu</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>GenHR+AV</th><td>9.38E-16</td><td>1.20E-15</td><td>1.32E-17</td></tr><tr><th>GenHR</th><td>7.82E-18</td><td>8.19E-18</td><td>6.69E-18</td></tr></tbody></table>

<table><thead><tr><th></th><th colspan="3">h</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>HR+AV</th><td>3.16E-16</td><td>3.59E-16</td><td>6.50E-16</td></tr><tr><th>HR</th><td>7.20E-19</td><td>7.31E-19</td><td>6.50E-19</td></tr></tbody></table>

<table><caption>Table 2: Numerical errors on $(h, hu)$ committed by the numerical schemes given in Table 1 and their viscous version (12) denoted +AV both for the at rest problem (50).</caption><thead><tr><th></th><th colspan="3">hu</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>HR+AV</th><td>1.10E-15</td><td>1.33E-15</td><td>3.14E-18</td></tr><tr><th>HR</th><td>8.77E-19</td><td>1.06E-18</td><td>1.40E-18</td></tr></tbody></table>

### 6.1.2 Moving equilibria

Our second numerical test concerns the preservation of the well-balanced property for moving equilibria defined by (47). Namely, we consider the following initial condition and topography:

$$ (hu)_0(x) = q_0, \quad u_0^2(x)/2 + g(h_0 + z)(x) = B_0, \quad z(x) = (5/2) \cos^2(4\pi x), \quad (51) $$

where $q_0 = 5/2$, $B_0 = 25/98 + 4g$. Such data involve a subsonic equilibrium that satisfies $u^2 \le gh$.

The space domain $[0, 1]$ is discretized with 400 cells and periodic conditions are imposed on both sides. The constant $C > 0$ in (46a) satisfies $C = 2080$ and the final time is 0.5. Figure 2 shows the compared results and Table 3 reports the numerical errors.
<figure><img src="image_6.png" /><figcaption>Figure 2: Numerical results and time evolution of the viscosities for the moving equilibrium (51) on a mesh composed of 400 cells. The legend of the numerical schemes is given in Table 1, +AV denotes the viscous version (12).</figcaption></figure>
<table><caption>Numerical errors on (h, hu) for the moving steady states (51).</caption></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">h</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>FullWB+AV</th><td>0.0E-00</td><td>0.0E-00</td><td>0.0E-00</td></tr><tr><th>FullWB</th><td>0.0E-00</td><td>0.0E-00</td><td>0.0E-00</td></tr></tbody></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">hu</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>FullWB+AV</th><td>0.0E-00</td><td>0.0E-00</td><td>0.0E-00</td></tr><tr><th>FullWB</th><td>0.0E-00</td><td>0.0E-00</td><td>0.0E-00</td></tr></tbody></table>

<table><caption>Numerical errors on (h, hu) for the moving steady states (51).</caption></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">h</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>SubBM+AV</th><td>3.36E-15</td><td>4.09E-15</td><td>9.72E-15</td></tr><tr><th>SubBM</th><td>2.11E-20</td><td>1.27E-19</td><td>1.62E-38</td></tr></tbody></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">hu</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>SubBM+AV</th><td>8.29E-14</td><td>8.33E-14</td><td>8.43E-14</td></tr><tr><th>SubBM</th><td>5.84E-19</td><td>1.40E-18</td><td>1.97E-36</td></tr></tbody></table>

<table><caption>Numerical errors on (h, hu) for the moving steady states (51).</caption></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">h</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>GenHR+AV</th><td>3.27E-19</td><td>4.66E-19</td><td>9.75E-19</td></tr><tr><th>GenHR</th><td>0.0E-00</td><td>0.0E-00</td><td>0.0E-00</td></tr></tbody></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">hu</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>GenHR+AV</th><td>0.00E-00</td><td>0.00E-00</td><td>0.00E-00</td></tr><tr><th>GenHR</th><td>0.00E-00</td><td>0.00E-00</td><td>0.00E-00</td></tr></tbody></table>

<table><caption>Numerical errors on (h, hu) for the moving steady states (51).</caption></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">h</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>HR+AV</th><td>1.93E-02</td><td>2.30E-02</td><td>5.24E-02</td></tr><tr><th>HR</th><td>1.93E-02</td><td>2.31E-02</td><td>5.24E-02</td></tr></tbody></table>

<table><caption>Table 3: Numerical errors on (h, hu) committed by the schemes detailed in Table 1 and by their viscous version (12) denoted +AV both for the problem (51).</caption><thead><tr><th rowspan="2"></th><th colspan="3">hu</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>HR+AV</th><td>0.56E-00</td><td>0.56E-00</td><td>0.58E-00</td></tr><tr><th>HR</th><td>0.56E-00</td><td>0.56E-00</td><td>0.58E-00</td></tr></tbody></table>

Remark that this equilibrium is preserved by all of the numerical schemes used here with the exception of HR. Since the HR scheme is not fully well-balanced, it does not preserve moving steady states (51). This means that the solution computed presents an error and the equilibrium is perturbed. This results in an artificial viscosity that make in order to satisfy an entropy inequality, as seen in the last picture on the right in Figure 2.

### 6.1.3 Evolution to a steady state

The third case test deals with the Goutal and Maurel problem [44]. The space domain $[0, 25]$ is discretized with 400 cells, the initial condition and the bottom topography are

$$h_0(x) = 2, \quad (hu)_0(x) = 4.42, \quad z(x) = \max(0, 0.2 - 0.05(x - 10)^2). \tag{52}$$

On the left boundary, homogeneous Neumann condition is prescribed for water height with the discharge $hu$ is set to 4.42. On the right boundary, the water height is set to 2 and a homogeneous Neumann condition is considered for the discharge. The initial condition will evolve until a final steady state is reached. The final time is fixed to 500 so that the steady state is attained. The exact solution is computed with the software SWASHES [34]. Figure 3 displays the compared results with $C = \frac{1}{4}$ in (46a). For the sake of clarity, the time evolutions of the ratio $\gamma/\lambda$ is only reported for $t \in [0, 100]$. Table 4 shows the numerical errors.
<figure><img src="image_7.png" /><figcaption>Figure 3: Numerical results and time evolution of the viscosities for the Goutal and Maurel's problem (52) on a mesh composed of 400 cells. The legend of the numerical schemes is given in Table 1, +AV denotes the viscous version (12).</figcaption></figure>
<table><caption>Numerical errors on (h, hu) for Goutal and Maurel's problem (52).</caption></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">h</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>FullWB+AV</th><td>3.86E-02</td><td>1.41E-02</td><td>2.00E-04</td></tr><tr><th>FullWB</th><td>3.86E-02</td><td>1.41E-02</td><td>2.00E-04</td></tr></tbody></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">hu</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>FullWB+AV</th><td>3.77E-02</td><td>2.10E-02</td><td>4.43E-04</td></tr><tr><th>FullWB</th><td>3.77E-02</td><td>2.10E-02</td><td>4.43E-04</td></tr></tbody></table>

<table><caption>Numerical errors on (h, hu) for Goutal and Maurel's problem (52).</caption></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">h</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>SubBM+AV</th><td>3.89E-14</td><td>1.19E-14</td><td>4.07E-17</td></tr><tr><th>SubBM</th><td>3.18E-15</td><td>7.09E-16</td><td>6.50E-19</td></tr></tbody></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">hu</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>SubBM+AV</th><td>1.02E-13</td><td>2.75E-14</td><td>4.93E-15</td></tr><tr><th>SubBM</th><td>2.22E-15</td><td>5.96E-16</td><td>2.85E-16</td></tr></tbody></table>

<table><caption>Numerical errors on (h, hu) for Goutal and Maurel's problem (52).</caption></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">h</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>GenHR+AV</th><td>3.32E-14</td><td>1.01E-14</td><td>2.47E-17</td></tr><tr><th>GenHR</th><td>1.54E-15</td><td>3.93E-16</td><td>6.50E-19</td></tr></tbody></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">hu</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>GenHR+AV</th><td>9.66E-14</td><td>2.46E-14</td><td>2.22E-15</td></tr><tr><th>GenHR</th><td>2.35E-15</td><td>5.28E-16</td><td>7.67E-17</td></tr></tbody></table>

<table><caption>Numerical errors on (h, hu) for Goutal and Maurel's problem (52).</caption></table>

<table><thead><tr><th rowspan="2"></th><th colspan="3">h</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>HR+AV</th><td>1.81E-02</td><td>6.89E-03</td><td>4.75E-05</td></tr><tr><th>HR</th><td>1.81E-02</td><td>6.89E-03</td><td>4.75E-05</td></tr></tbody></table>

<table><caption>Table 4: Numerical errors on (h, hu) committed by the schemes detailed in Table 1 and their viscous version (12) denoted +AV both for the problem (52)</caption><thead><tr><th rowspan="2"></th><th colspan="3">hu</th></tr><tr><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><th>HR+AV</th><td>2.97E-02</td><td>1.67E-02</td><td>2.79E-04</td></tr><tr><th>HR</th><td>2.97E-02</td><td>1.67E-02</td><td>2.79E-04</td></tr></tbody></table>

The numerical results are very satisfying for the GenHR and the SubBM schemes. Once again, a large error is observed for the Hydrostatic Reconstruction (HR) scheme because this scheme is not well-balanced for the moving equilibrium induced by the Goutal and Maurel problem (52). The large numerical error observed for the FullWB scheme is an intrinsic default of the scheme (see [58, Section 7.3] for the details). The artificial viscosity method does not correct this default.

### 6.1.4 Dam-break problem

Our final test case for the shallow water system concerns the break dam problem as detailed in [61]. The spatial domain $[-1, 1]$ is discretized with 400 cells and the initial condition and topography are

$$h_0(x) + z(x) = \begin{cases} 3 & \text{if } x < 0, \\ 1 & \text{otherwise,} \end{cases} \quad u_0(x) = 0, \quad \text{with,} \quad z(x) = \frac{1}{2} \cos^2(\pi x). \quad (53)$$

Homogeneous Neumann boundary conditions are prescribed on both boundaries. The final time is 0.1. A reference solution is computed on a fine grid made of 50 000 cells with the standard HLL scheme [46] coupled to the discrete source term $(-gh\partial_x z)(x_i, t^n) \approx -g \frac{h_i^n+h_{i+1}^n}{2} \frac{z_{i+1}-z_i}{\Delta x}$. Figure 4 displays the results with $C = 17$ in (46a).
<figure><img src="image_8.png" /><figcaption>Figure 4: Numerical results and time evolution of the viscosities for the break dam problem (53) on a mesh composed of 400 cells. The legend of the numerical schemes is given in Table 1, +AV denotes the viscous version (12).</figcaption></figure>

Since the non-corrected FullWB scheme satisfies a discrete entropy inequality, the viscous version coincides with the initial scheme.
## 6.2 Two-layer shallow water equations

Let us now consider the two-layer shallow water system that models two inmiscible fluids in a shallow regime of constant density $\rho_j$ and thickness $h_j$ ($j = 1, 2$). We assume $\rho_1 < \rho_2$ so that index $j = 1$ refers to the upper layer, while $j = 2$ corresponds to the lower layer. We shall denote by $u_j$ the horizontal velocity within each layer and by $q_j = h_j u_j$ the discharge. The fixed bottom topography will be denoted as $z(x)$.

The equations of the system may be written as

$$ egin{cases} 
\partial_t h_1 + \partial_x q_1 = 0, \\
\partial_t q_1 + \partial_x \left( \frac{q_1^2}{h_1} + \frac{1}{2} g h_1^2 \right) = -gh_1 \partial_x h_2 - gh_1 \partial_x z, \\
\partial_t h_2 + \partial_x q_2 = 0, \\
\partial_t q_2 + \partial_x \left( \frac{q_2^2}{h_2} + \frac{1}{2} g h_2^2 \right) = -\frac{\rho_1}{\rho_2} g h_2 \partial_x h_1 - gh_2 \partial_x z.
\end{cases} 
\quad (54) $$

This system can also be written in a compact form

$$ \partial_t w + \partial_x f(w) + A(w) \partial_x w = S(w) \partial_x z, \quad (55) $$

where

$$ w = \begin{pmatrix} h_1 \\ q_1 \\ h_2 \\ q_2 \end{pmatrix}, \quad f(w) = \begin{pmatrix} q_1 \\ \frac{q_1^2}{h_1} + \frac{1}{2} g h_1^2 \\ q_2 \\ \frac{q_2^2}{h_2} + \frac{1}{2} g h_2^2 \end{pmatrix}, $$

$$ A(w) = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & gh_1 & 0 \\ 0 & 0 & 0 & 0 \\ grh_2 & 0 & 0 & 0 \end{pmatrix}, \quad S(w) = \begin{pmatrix} 0 \\ -gh_1 \\ 0 \\ -gh_2 \end{pmatrix}, $$

with $r = \rho_1/\rho_2$.

The usual entropy-pair associated with the above two-layer shallow water system is given by

$$ \hat{\eta}(w, z) = \sum_{j=1}^{2} \rho_j h_j \left( \frac{u_j^2}{2} + g \frac{h_j}{2} + z \right) + g \rho_1 h_1 h_2 \quad (56a) $$

$$ \hat{G}(w, z) = \sum_{j=1}^{2} \rho_j h_j u_j \left( \frac{u_j^2}{2} + gh_j + z \right) + \rho_1 g h_1 h_2 (u_1 + u_2). \quad (56b) $$

We remark that the function $(w, x) \mapsto \eta(w, z(x))$ does not satisfy the required convexity and we may introduce some additional convex correction with respect to $x$ as proposed in (46). However, we take the opportunity to test our numerical method out of the validity domain imposed by Theorem 2.1. As a consequence, we here suggest to adopt the non-convex with respect to $z$ but convex with respect to $w$ entropy function given by (56).

Concerning the steady solutions, we shall focus on lake-at-rest steady states given by

$$ u_1 = u_2 = 0, \quad h_2 + z = Cst, \quad h_1 = Cst. $$

We consider a path-conservative numerical scheme [66], with segments as family of paths combined with a hydrostatic reconstruction [2, 27]. More explicitly, we define

$$ w_{i+1/2}^{\pm} = ((h_1)_{i+1/2}^{\pm}, (q_1)_{i+1/2}^{\pm}, (h_2)_{i+1/2}^{\pm}, (q_2)_{i+1/2}^{\pm})^t, $$
$$ (h_2)_{i+1/2}^{-} = \max(0, (h_2)_i + z_i - z_{i+1/2}), \quad (h_2)_{i+1/2}^{+} = \max(0, (h_2)_{i+1} + z_{i+1} - z_{i+1/2}), \\ (h_1)_{i+1/2}^{-} = (h_1)_i, \quad (h_1)_{i+1/2}^{+} = (h_1)_{i+1}, \quad z_{i+1/2} = \max(z_i, z_{i+1}), $$

and we consider a Roe path-conservative scheme written in a PVM formalism [19]:

$$ w_i^{n+1} = w_i^n - \frac{\Delta t}{\Delta x} \left( f_\Delta(w_{i+1/2}^{-}, w_{i+1/2}^{+}) - f_\Delta(w_{i-1/2}^{-}, w_{i-1/2}^{+}) + \Pi_{i+1/2}^{-} + \Pi_{i-1/2}^{+} \right) \\ - \frac{\Delta t}{2\Delta x} \left( A_\Delta^L(w_{i+1/2}^{-}, w_{i+1/2}^{+}) \cdot (w_{i+1/2}^{+} - w_{i+1/2}^{-}) + A_\Delta^R(w_{i-1/2}^{-}, w_{i-1/2}^{+}) \cdot (w_{i-1/2}^{+} - w_{i-1/2}^{-}) \right), \quad (57) $$

where the numerical flux is given by

$$ f_{\Delta}(w_{i+1/2}^{-}, w_{i+1/2}^{+}) = \frac{1}{2} \left( f(w_{i+1/2}^{-}) + f(w_{i+1/2}^{+}) \right) - \frac{1}{2} |A_{\Delta}| \cdot (w_{i+1/2}^{+} - w_{i+1/2}^{-}), \quad (58) $$

$$ \Pi_{i+1/2}^{-} = \left( 0, 0, 0, \frac{g}{2}((h_2)_i - (h_2)_i^2) \right)^t, \quad \Pi_{i-1/2}^{+} = \left( 0, 0, 0, \frac{g}{2}((h_2)_i - (h_2)_i^2) \right)^t, $$

with $|A_\Delta^n|$ the absolute value of the Roe matrix

$$ A_{\Delta} = \begin{pmatrix} 0 & 1 & 0 & 0 \\ gh_{\Delta,1} - u_{\Delta,1}^2 & 2u_{\Delta,1} & gh_{\Delta,2} & 0 \\ 0 & 0 & 0 & 1 \\ rgh_{\Delta,1} & 0 & gh_{\Delta,2} - u_{\Delta,2}^2 & 2u_{\Delta,2} \end{pmatrix}, \quad (59) $$

$h_{\Delta,j}$ and $u_{\Delta,j}$, $j = 1, 2$ are the Roe averages and are defined as

$$ h_{\Delta,j} = \frac{h_{i,j} + h_{i+1,j}}{2}, \quad u_{\Delta,j} = \frac{\sqrt{h_{i,j}}u_{i,j} + \sqrt{h_{i+1,j}}u_{i+1,j}}{\sqrt{h_{i,j}} + \sqrt{h_{i+1,j}}}, \quad j = 1, 2. $$

Finally,

$$ A_{\Delta}^{L}(w_{i}^{n}, w_{i+1}^{n}) = A_{\Delta}^{R}(w_{i}^{n}, w_{i+1}^{n}) = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & gh_{\Delta,1} & 0 \\ 0 & 0 & 0 & 0 \\ rgh_{\Delta,2} & 0 & 0 & 0 \end{pmatrix}. \quad (60) $$

### 6.2.1 Internal dam-break with small density ratio

Let us first consider a test case with $r = 0.2$. The spatial domain is set to $[-3, 3]$ which is discretized with 500 volumes. The initial condition and bottom are set as

$$ z(x) = 0.5 \exp(-x^2), \quad h_1(x, 0) + h_2(x, 0) + z(x) = 2, \quad h_2(x, 0) + z(x) = \begin{cases} 0.6, & \text{if } x < 0, \\ 1.85, & \text{otherwise,} \end{cases} $$

with zero discharge $q_1(x, 0) = q_2(x, 0) = 0$. Open boundary conditions are considered and the final time is set to $t = 1$.

Figures 5 and 6 show the surface and interface evolution for $t = 0.5$ and $t = 1$ respectively. As we see, the Roe scheme does not satisfy an entropy inequality and it results in a non-physical shock in the interface corresponding to the maximum of the topography. Nevertheless, the AV correction allows to overcome this problem. The corresponding AV correction is shown in Figure 7.

### 6.2.2 Internal dam-break with big density ratio

Let us consider now the same initial condition used in previous example, but now we set the ratio of densities to $r = 0.9$. When the ratio of densities is near 1, the perturbations on the free surface due to the internal dam-break are less noticeable. The internal waves propagate with smaller speed, while the overall dynamics is similar, observing an internal shock that is formed in the interface. This is shown in Figure 8.

In this case, the artifact that can be observed in previous case is now less noticeable and we have to zoom into the interface to observe it. In any case, the AV correction allows to overcome this problem again. The corresponding AV correction is shown in Figure 9.
<figure><img src="image_9.png" /><figcaption>Figure 5: Surface $h_1 + h_2 + z$ and interface $h_2 + z$ evolution of an internal dambreak problem for $r = 0.2$. Comparison for Roe scheme with and without the artificial viscosity method.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 6: Surface $h_1 + h_2 + z$ and interface $h_2 + z$ evolution of an internal dambreak problem for $r = 0.2$. Comparison for Roe scheme with and without the artificial viscosity method.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 7: Time evolution of the maximum ratio $\gamma/\lambda$ for $r = 0.2$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 8: Surface $h_1 + h_2 + z$ and interface $h_2 + z$ evolution of an internal dambreak problem for $r = 0.9$. Comparison for Roe scheme with and without the artificial viscosity method.</figcaption></figure>

# 7 Conclusions

In this work, we have presented a general and efficient framework to ensure discrete entropy stability for first-order well-balanced finite volume schemes applied to systems of balance laws, including those with non-conservative products. Building upon the artificial viscosity approach introduced by Tadmor and further developed in recent literature, we have proposed a minimal correction that can be systematically added to a wide range of existing well-balanced schemes.

The proposed method maintains the well-balanced property and ensures robustness and guarantees entropy stability through a carefully constructed artificial viscosity coefficient. The theoretical analysis confirms that the corrected scheme satisfies discrete entropy inequalities, preserves stationary solutions, and remains robust under a suitable CFL condition.

A reformulation of the modified scheme as a Godunov-type method allowed us to derive sufficient conditions for entropy preservation, well-balancing and robustness. Numerical experiments, particularly on the shallow water equations and the two-layer shallow water model, have demonstrated the effectiveness of the proposed method. The results confirm that entropy stability can be achieved without compromising the accuracy or well-balanced properties of the base schemes.

An important contribution of this work is its flexibility and generality: the proposed correction can be combined, for example, with the strategy described in [28] to extend any consistent first-order solver for the homogeneous problem—i.e., without source terms—into a well-balanced, robust, and
<figure><img src="image_13.png" /><figcaption>Figure 9: Time evolution of the maximum ratio $\gamma/\lambda$ for $r = 0.9$.</figcaption></figure>

entropy-satisfying solver for general one-dimensional balance laws. To the best of our knowledge, this is the first time such a general and constructive procedure is presented, offering a unified framework to develop stable and physically relevant numerical methods for a wide range of balance law systems.

Overall, this approach provides a practical and theoretically sound procedure to enhance existing well-balanced schemes with entropy-preserving capabilities, offering potential for future extensions to higher-order methods.

## Acknowledgments

This work is supported by projects PCI2024-155061-2 and PID2022-137637NB-C21 funded by MCIN/AEI/10.13039/501100011033/ and ERDF A way of making Europe. The authors thank the support of the Centre National de Recherche Scientifique, International Research Project hyPerbolIC models, numerical Analysis and Scientific computation (CNRS, IRP PICASSO).

## References

[1] R. Abgrall. A review of residual distribution schemes for hyperbolic and parabolic problems: the July 2010 state of the art. Commun. Comput. Phys., 11(4):1043–1080, 2012.

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

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

[4] C. Berthon. Schéma nonlinéaire pour l'approximation numérique d'un système hyperbolique non conservatif. Comptes Rendus Mathematique, 335(12):1069–1072, 2002.

[5] C. Berthon, M.J Castro, A. Duran, T. Morales De Luna, and K. Saleh. Artificial viscosity to get both robustness and discrete entropy inequalities. Journal of Scientific Computing, 97(65), 2023.

[6] C. Berthon and C. Chalons. A fully well-balanced, positive and entropy-satisfying Godunov-type method for the shallow-water equations. Mathematics of Computation, 85, 02 2015.

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

[8] C. Berthon, F. Coquel, and P. G. LeFloch. Why many theories of shock waves are necessary: kinetic relations for non-conservative systems. Proc. Roy. Soc. Edinburgh Sect. A, 142(1):1–37, 2012.
[9] C. Berthon, M. M'Baye, M. Le, and D. Seck. A well-defined moving steady states capturing Godunov-type scheme for shallow-water model. *Int. J. Finite Volumes*, 01 2021.

[10] C. Berthon, V. Michel-Dansac, and A. Thomann. Towards a fully well-balanced and entropy-stable scheme for the euler equations with gravity: preserving isentropic steady solutions. *Mathematics of Computation*, 2025.

[11] A. Bihlo and S. MacLachlan. Well-balanced mesh-based and meshless schemes for the shallow-water equations, 2017.

[12] F. Bouchut. Entropy satisfying flux vector splittings and kinetic bgk models. *Numerische Mathematik*, 94(4):623-672, 2003.

[13] F. Bouchut. *Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources*. Frontiers in Mathematics. Birkhäuser Verlag, Basel, 2004.

[14] F. Bouchut and T. Morales de Luna. An entropy satisfying scheme for two-layer shallow water equations with uncoupled treatment. *ESAIM: Mathematical Modelling and Numerical Analysis*, 42:683–698, 07 2008.

[15] F. Bouchut and T. Morales de Luna. A subsonic-well-balanced reconstruction scheme for shallow water flows. *SIAM J. Numerical Analysis*, 48:1733–1758, 01 2010.

[16] S. Bryson, Y. Epshteyn, A. Kurganov, and G. Petrova. Well-balanced positivity preserving central-upwind scheme on triangular grids for the Saint-Venant system. *ESAIM Math. Model. Numer. Anal.*, 45(3):423-446, 2011.

[17] T. Buffard, T. Gallouet, and J.-M. Hérard. Un schéma simple pour les équations de saint-venant. *Comptes Rendus de l'Académie des Sciences - Series I - Mathematics*, 326(3):385-390, 1998.

[18] A. Canestrelli, A. Siviglia, M. Dumbser, and E. F. Toro. Well-balanced high-order centred schemes for non-conservative hyperbolic systems. applications to shallow water equations with fixed and mobile bed. *Advances in Water Resources*, 32(6):834-844, 2009.

[19] M. J. Castro and E. Fernández-Nieto. A class of computationally fast first order finite volume solvers: PVM methods. *SIAM Journal on Scientific Computing*, 34(4):A2173–A2196, January 2012.

[20] M. J. Castro, E. Fernández-Nieto, J. González-Vida, and C. Parés-Madroñal. Numerical treatment of the loss of hyperbolicity of the two-layer shallow-water system. *Journal of Scientific Computing*, 48(1):16–40, 2011.

[21] M. J. Castro, Jörn T. Frings, S. Noelle, C. Parés, and G. Puppo. On the hyperbolicity of two- and three-layer shallow water equations. *Hyperbolic problems, theory, numerics and applications*, 1:337–345, 2012.

[22] M. J. Castro, J. M Gallardo, J. A López-García, and C. Parés. Well-balanced high order extensions of godunov's method for semilinear balance laws. *SIAM Journal on Numerical Analysis*, 46(2):1012–1039, 2008.

[23] M. J. Castro, P. LeFloch, M. L. Muñoz-Ruiz, and C. Parés. Why many theories of shock waves are necessary: Convergence error in formally path-consistent schemes. *Journal of Computational Physics*, 227(17):8107–8129, 2008.

[24] M. J. Castro, J. Macías, and C. Parés. A q-scheme for a class of systems of coupled conservation laws with source term. application to a two-layer 1-d shallow water system. *ESAIM: Mathematical Modelling and Numerical Analysis*, 35(1):107–127, 2001.
[25] M. J. Castro, J. Macías, C. Parés, J. García-Rodríguez, and E. Vázquez-Cendón. A two-layer finite volume model for flows through channels with irregular geometry: computation of maximal exchange solutions: application to the strait of gibraltar. *Communications in Nonlinear Science and Numerical Simulation*, 9(2):241–249, 2004.

[26] M. J. Castro, T. Morales de Luna, and C. Parés. Well-Balanced Schemes and Path-Conservative Numerical Methods. In Rémi Abgrall and Chi-Wang Shu, editor, *Handbook of Numerical Analysis*, volume 18 of *Handbook of Numerical Methods for Hyperbolic Problems Applied and Modern Issues*, pages 131–175. Elsevier, 2017. DOI: 10.1016/bs.hna.2016.10.002.

[27] M. J. Castro, A. Pardo, and C. Parés. Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique. *Mathematical Models and Methods in Applied Sciences*, 5:2055–2113, 12 2007.

[28] M.J. Castro and C. Parés. Well-balanced high-order finite volume methods for systems of balance laws. *J. Sci. Comput.*, 82(2), 2020.

[29] C. Chalons. Path-conservative in-cell discontinuous reconstruction schemes for non conservative hyperbolic systems. *Communications in Mathematical Sciences*, 18(1):1–30, 2020.

[30] C. Chalons and F. Coquel. Navier-stokes equations with several independent pressure laws and explicit predictor-corrector schemes. *Numerische Mathematik*, 101(3):451–478, 2005.

[31] F. Coquel and B. Perthame. Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics. *SIAM Journal on Numerical Analysis*, 35(6):2223–2249, 1998.

[32] C. M. Dafermos. *Hyperbolic conservation laws in continuum physics*, volume 325 of *Grundlehren der Mathematischen Wissenschaften* [Fundamental Principles of Mathematical Sciences]. Springer-Verlag, Berlin, third edition, 2010.

[33] G. Dal Maso, P. G. LeFloch, and F. Murat. Definition and weak stability of nonconservative products. *Journal de mathématiques pures et appliquées*, 74(6):483–548, 1995.

[34] O. Delestre, C. Lucas, P.A. Ksinant, F. Darboux, C. Laguerre, T.N.T Vo, F. James, and S. Cordier. SWASHES: a compilation of Shallow Water Analytic Solutions for Hydraulic and Environmental Studies. *International Journal for Numerical Methods in Fluids*, 72(3):269–300, May 2013.

[35] V. Desveaux, M. Zenk, C. Berthon, and C. Klingenberg. A well-balanced scheme to capture non-explicit steady states in the Euler equations with gravity. *International Journal for Numerical Methods in Fluids*, 81(2):104–127, 2016.

[36] V. Desveaux, M. Zenk, C. Berthon, and C. Klingenberg. Well-balanced schemes to capture non-explicit steady states: Ripa model. *Mathematics of Computation*, 85(300):1571–1602, 2016.

[37] U.S. Fjordholm, S. Mishra, and E. Tadmor. Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography. *Journal of Computational Physics*, 230(14):5587–5609, 2011.

[38] T. Gallouët, J.M. Hérard, and N. Seguin. Some approximate Godunov schemes to compute shallow-water equations with topography. *Computers and Fluids*, 32(4):479–513, 2003.

[39] E. Godlewski and P.-A. Raviart. *Hyperbolic systems of conservation laws*, volume 3/4 of *Mathématiques & Applications (Paris) [Mathematics and Applications]*. Ellipses, Paris, 1991.

[40] E. Godlewski and P.-A. Raviart. *Numerical approximation of hyperbolic systems of conservation laws*, volume 118 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 1996.

[41] S. K. Godunov. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics. *Mat. Sb. (N.S.)*, 47(89):271–306, 1959.
[42] Irene Gómez-Bueno, Manuel Jesús Castro, Carlos Parés, and Giovanni Russo. Collocation methods for high-order well-balanced methods for systems of balance laws. *Mathematics*, 9(15):1799, jul 2021.

[43] L. Gosse. A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms. *Computers & Mathematics with Applications*, 39(9):135–159, 2000.

[44] N. Goutal and F. Maurel. *Proceedings of the 2nd workshop on Dam-Break Wave Simulation*. Department Laboratoire National d'Hydraulique, Groupe Hydraulique Fluviale, 1997.

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

[46] A. Harten, P.D. Lax, and B. Van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. *SIAM review*, 25:35–61, 1983.

[47] J. Hou, Q. Liang, F. Simons, and R. Hinkelmann. A stable 2D unstructured shallow flow model for simulations of wetting and drying over rough terrains. *Computers and Fluids*, 82:132–147, 2013.

[48] T. Y. Hou and P. G. LeFloch. Why nonconservative schemes converge to wrong solutions: error analysis. *Mathematics of computation*, 62(206):497–530, 1994.

[49] R. Käppeli and S. Mishra. Well-balanced schemes for the Euler equations with gravitation. *Journal of Computational Physics*, 259(0):199–219, 2014.

[50] B. Khobalatte and B. Perthame. Maximum principle on the entropy and second-order kinetic schemes. *Math. Comp.*, 62(205):119–131, 1994.

[51] P.D. Lax. Shock waves and entropy. In *Contributions to nonlinear functional analysis (Proc. Sympos., Math. Res. Center, Univ. Wisconsin, Madison, Wis., 1971)*, pages 603–634. Academic Press, New York, 1971.

[52] P.D. Lax. *Hyperbolic systems of conservation laws and the mathematical theory of shock waves*. Society for Industrial and Applied Mathematics, Philadelphia, Pa., 1973. Conference Board of the Mathematical Sciences Regional Conference Series in Applied Mathematics, No. 11.

[53] P. LeFloch. Entropy weak solutions to nonlinear hyperbolic systems under nonconservative form. *Communications in partial differential equations*, 13(6):669–727, 1988.

[54] P. LeFloch and T.-P. Liu. Existence theory for nonlinear hyperbolic systems in nonconservative form. *Forum Math.*, 5:261–280, 1993.

[55] P. LeFloch and M. D. Thanh. Nonlinear hyperbolic systems of balance laws. In *The Riemann Problem in Continuum Physics*, pages 41–85. Springer, 2024.

[56] M. Li, R. Mu, and H. Dong. A well-balanced discontinuous Galerkin method for the shallow water flows on erodible bottomimage 1. *Computers and Mathematics with Applications*, 119:13–20, 2022.

[57] Q. Liang and F. Marche. Numerical resolution of well-balanced shallow water equations with complex source terms. *Advances in Water Resources*, 32(6):873–884, 2009.

[58] L. Martaud and C. Berthon. How to enforce an entropy inequality of (fully) well-balanced godunov-type schemes for the shallow water equations. *ESAIM: Mathematical Modelling and Numerical Analysis*, 59(2):955–997, 2025.

[59] V. Michel-Dansac, C. Berthon, S. Clain, and F. Foucher. A well-balanced scheme for the shallow-water equations with topography. *Computers and Mathematics with Applications*, 72, 09 2015.
[60] V. Michel-Dansac, C. Berthon, S. Clain, and F. Foucher. A well-balanced scheme for the shallow-water equations with topography. *Computers & Mathematics with Applications*, 72(3):568-593, 2016.

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

[62] L. O Müller, C. Parés, and E. F. Toro. Well-balanced high-order numerical schemes for one-dimensional blood flow in vessels with varying mechanical properties. *Journal of Computational Physics*, 2013.

[63] C. Parés. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. *SIAM Journal on Numerical Analysis*, 44(1):300-321, 2006.

[64] C. Parés and M. J. Castro. On the well-balance property of Roe's method for nonconservative hyperbolic systems. Applications to shallow-water systems. *ESAIM: mathematical modelling and numerical analysis*, 38(05):821-852, 2004.

[65] C. Pares and E. Pimentel. The riemann problem for the shallow water equations with discontinuous topography: the wet-dry case. *Journal of Computational Physics*, 378:344-365, 2019.

[66] Carlos Parés. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. *SIAM J. Numer. Anal.*, 44(1):300-321, 2006.

[67] Michael D. Perlman. Jensen's inequality for a convex vector-valued function on an infinite-dimensional space. *Journal of Multivariate Analysis*, 4(1):52-65, 1974.

[68] E. Pimentel-García, M. J. Castro, C. Chalons, T. Morales de Luna, and C. Parés. In-cell discontinuous reconstruction path-conservative methods for non conservative hyperbolic systems-second-order extension. *Journal of Computational Physics*, 459:111152, 2022.

[69] E. Pimentel-García, M. J. Castro, C. Chalons, and C. Parés. High-order in-cell discontinuous reconstruction path-conservative methods for nonconservative hyperbolic systems-dr. mood method. *Numerical Methods for Partial Differential Equations*, 40(6):e23133, 2024.

[70] A. Prieto-Arranz, L. Ramirez, I. Couceiro, I. Colominas, and X. Nogueira. A well-balanced sph-ale scheme for shallow water applications. *Journal of Scientific Computing*, 88, 09 2021.

[71] P.-A. Raviart and L. Sainsaulieu. A nonconservative hyperbolic system modeling spray dynamics. Part I: Solution of the Riemann problem. *Mathematical Models and Methods in Applied Sciences*, 5(03):297-333, 1995.

[72] L Sainsaulieu. Ondes progressives solutions de systemes convectifs-diffusifs et systemes hyperboliques non conservatifs. *Comptes rendus de l'Académie des sciences. Série 1, Mathématique*, 312(7):491-494, 1991.

[73] D. Serre. *Systems of conservation laws. 1*. Cambridge University Press, Cambridge, 1999. Hyperbolicity, entropies, shock waves, Translated from the 1996 French original by I. N. Sneddon.

[74] E. Tadmor. Numerical viscosity and the entropy condition for conservative difference schemes. *Mathematics of Computation*, 43(168):369-381, 1984.

[75] E. Tadmor. The numerical viscosity of entropy stable schemes for systems of conservation laws. *i. Mathematics of Computation*, 49(179):91-103, 1987.

[76] E. Tadmor. Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems. *Acta Numerica*, 12(1):451-512, 2003.
[77] E Tadmor. Entropy stable schemes. *Handbook of Numerical Methods for Hyperbolic Problems: Basic and Fundamental Issues*, edited by R. Abgrall and C.-W. Shu (North-Holland, Elsevier, Amsterdam, 2017), 17:467–493, 2016.

[78] E. F. Toro. *Riemann solvers and numerical methods for fluid dynamics*. Springer-Verlag, Berlin, third edition, 2009. A practical introduction.

[79] E.F. Toro, M. Spruce, and W. Speares. Restoration of the contact surface in the HLL-Riemann solver. *Shock waves*, 4(1):25–34, 1994.

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

[81] Y. Xing, X. Zhang, and C.W. Shu. Positivity-preserving high-order well-balanced discontinuous Galerkin methods for the shallow water equations. *Advances in Water Resources*, 33(12):1476–1493, 2010.