# How to enforce an entropy inequality of (fully) well-balanced Godunov-type schemes for the shallow water equations

Ludovic Martaud, Christophe Berthon

► To cite this version:

Ludovic Martaud, Christophe 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, 2025, 59, pp.955 - 997. 10.1051/m2an/2025012. hal-05376635

HAL Id: hal-05376635

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

Submitted on 21 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 ENFORCE AN ENTROPY INEQUALITY OF (FULLY) WELL-BALANCED GODUNOV-TYPE SCHEMES FOR THE SHALLOW WATER EQUATIONS

LUDOVIC MARTAUD<fnref n="1" />,* AND CHRISTOPHE BERTHON<fnref n="2" />

**Abstract.** This work concerns the design of well-balanced entropy-stable numerical schemes for the shallow water equations. The fully discrete entropy inequality is reached by introducing a local entropy condition incorporated in the scheme design. The source term is discretized to preserve both steady states and entropy stability. The method yields explicit schemes which are relevantly illustrated with several test cases.

**Mathematics Subject Classification.** 65M08, 65M12.

Received February 21, 2024. Accepted February 19, 2025.

## 1. INTRODUCTION

The present work is devoted to the numerical approximation of the weak solutions of the shallow water equations with the source term related to topography in one space dimension given by

$$ \partial_t \left( \frac{h}{hu} \right) + \partial_x \left( hu^2 + gh^2/2 \right) = \left( \begin{array}{c} 0 \\ -gh\partial_x z \end{array} \right). \qquad (1) $$

This model governs the water height $h \ge 0$ and the velocity $u \in \mathbb{R}$ of a fluid. In order to deal with dry areas, the water velocity is defined as follows (see [4, 41] for the details):

$$ u = \begin{cases} \frac{(hu)}{h}, & \text{if } h > 0, \\ 0, & \text{otherwise.} \end{cases} $$

The gravitational constant is $g > 0$ and $z : \mathbb{R} \to \mathbb{R}$ is a given time-independent smooth topography function. The unknown state vector $w = (h, hu)^T$ is assumed to be in the convex set $\Omega = \{(h, hu) \in \mathbb{R}^2 | h \ge 0, hu \in \mathbb{R}\}$. The system is endowed with given initial data $w_0 : \mathbb{R} \to \Omega$ at time $t = 0$. As a consequence, we consider the

<footnote>Keywords and phrases. Shallow water equations, finite volume schemes, discrete entropy stability, well-balanced schemes.</footnote>

<footnote n="1">¹ Univ Rennes, Inria Bretagne Atlantique, Rennes, France.</footnote>

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

<footnote>*Corresponding author: ludovic.martaud@inria.fr</footnote>
following Cauchy problem:

$$ \left\{ \begin{array}{ll} \partial_t w + \partial_x f(w) = S(w, z), & x \in \mathbb{R}, t > 0, \\ w(x, t = 0) = w_0(x), & x \in \mathbb{R}, \end{array} \right. \qquad (2) $$

where we have set

$$ f(w) = (hu, hu^2 + gh^2/2)^{\mathrm{T}}, \quad \text{and} \quad S(w, z) = (0, -gh\partial_x z)^{\mathrm{T}}. \qquad (3) $$

For the sake of clarity in the forthcoming notations, we introduce $\hat{w} = (h, hu, z)^{\mathrm{T}}$, which takes its values in the convex set $\hat{\Omega}$ defined by

$$ \hat{\Omega} = \{(h, hu, z) \in \mathbb{R}^3 \mid h \ge 0, hu \in \mathbb{R}, z \in \mathbb{R}\}. $$

In the flat regions (i.e. $\partial_x z = 0$), it is well-known [12] that the homogeneous shallow water system is endowed with the following entropy inequality:

$$ \partial_t \eta(w) + \partial_x G(w) \le 0, \qquad (4) $$

with

$$ \eta(w) = hu^2/2 + gh^2/2 \quad \text{and} \quad G(w) = hu^3/2 + gh^2u. \qquad (5) $$

A smooth generic function $z$ modifies the above inequality with the term $-ghu\partial_x z$ in its right-hand side. As a consequence, the entropy inequality associated to shallow water system (1) now reads

$$ \partial_t \eta(w) + \partial_x G(w) \le -ghu\partial_x z. \qquad (6) $$

According to [12] (for instance, see also [20]) and since $\partial_t z = 0$, the inequality (6) reformulates equivalently in a conservative form given by

$$ \partial_t (\eta(w) + ghz) + \partial_x (G(w) + g(hu)z) \le 0. $$

Let us set

$$ \hat{\eta}(\hat{w}) = \frac{(hu)^2}{2h} + \frac{gh^2}{2} + ghz, \quad \hat{G}(\hat{w}) = \frac{(hu)^3}{2h^2} + g(hu)(h+z), \quad \forall \hat{w} \in \hat{\Omega}, \qquad (7) $$

so that the above entropy inequality now reads

$$ \partial_t \hat{\eta}(\hat{w}) + \partial_x \hat{G}(\hat{w}) \le 0. \qquad (8) $$

Therefore, the presence of the source term $S(w, z)$ in the system (1) leads to the entropy inequality (8) that contains (4) and a contribution of the topography function.

In addition, the presence of the source term $S(w, z)$ involves the existence of non-trivial stationary solutions that satisfy

$$ hu = \text{cst}, \quad \frac{u^2}{2} + g(h+z) = \text{cst}. \qquad (9) $$

Among these steady states, a special attention is paid to the lake at rest (for instance, see [11, 22, 34, 37, 38, 43, 46] for a nonhexaustive list) which is given by

$$ u = 0, \quad h + z = \text{cst}. \qquad (10) $$

From a numerical point of view, the solutions of the shallow water system (1) are approximated on uniform space meshes $(x_{i+\frac{1}{2}})_{i \in \mathbb{Z}}$ in $\mathbb{R}$ of constant size $\Delta x > 0$. Thus, we have $x_{i+\frac{1}{2}} = x_{i-\frac{1}{2}} + \Delta x$ for all $i \in \mathbb{Z}$. Uniform meshes in time $(t^n)_{n \in \mathbb{N}}$ in $[0, +\infty)$ of constant size $\Delta t > 0$ are also considered and they satisfy $t^{n+1} = t^n + \Delta t$
for all $n$ in $\mathbb{N}$. At time $t^0 = 0$, the initial condition $w_0$ and the given smooth function $z$ are discretized by $((w_i^0, z_i))_{i \in \mathbb{Z}}$ in $\hat{\Omega}$ such that

$$ (w_i^0, z_i)^{\mathrm{T}} = \frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} (w_0, z)^{\mathrm{T}}(x) \mathrm{d}x, \quad \forall i \in \mathbb{Z}. \qquad (11) $$

The sequence $((w_i^0, z_i)^{\mathrm{T}})_{i \in \mathbb{Z}}$ defines a piecewise constant approximation of $(w(\cdot, t=0), z)^{\mathrm{T}}$. As a consequence, a numerical approximation of $w(\cdot, t^{n+1})$ is entirely defined by a numerical scheme that gives the updated sequence $(w_i^{n+1})_{i \in \mathbb{Z}}$ from the sequences $(w_i^n)_{i \in \mathbb{Z}}$ and $(z_i)_{i \in \mathbb{Z}}$. However, a suitable updated sequence $(w_i^{n+1})_{i \in \mathbb{Z}}$ has to satisfy some properties. In order to give these properties and for the sake of clarity, the notation $\hat{w}_i^n = (w_i^n, z_i)^{\mathrm{T}}$ is now being considered but we emphasize that $z_i$ is a given quantity.

The sequence $(\hat{w}_i^{n+1})_{i \in \mathbb{Z}}$ has to be well-balanced that means it exactly preserves stationary solutions. On the one hand, the well-balanced property for the lake at rest (10) writes

$$ \text{If } \forall i \in \mathbb{Z}, \quad u_i^n = 0 \quad \text{and} \quad h_i^n + z_i = \text{cst} \quad \text{then} \quad \forall i \in \mathbb{Z}, \quad w_i^{n+1} = w_i^n. \qquad (12) $$

Several schemes satisfying this property have been proposed during the last two decades (for instance, see [11, 22, 34, 37, 38, 43, 46]). On the other hand, the well-balanced property for the moving equilibrium (9) is given by

$$ \text{If } \forall i \in \mathbb{Z}, \quad (hu)_i^n = \text{cst} \quad \text{and} \quad \frac{(u_i^n)^2}{2} + g(h_i^n + z_i) = \text{cst} \quad \text{then} \quad \forall i \in \mathbb{Z}, \quad w_i^{n+1} = w_i^n. \qquad (13) $$

For instance, such a property is satisfied by the schemes described in [15, 25, 35, 41].

In addition to the well-balanced property, the sequence $(\hat{w}_i^{n+1})_{i \in \mathbb{Z}}$ has to satisfy a discrete entropy inequality. Denoting $\hat{G}_{i+\frac{1}{2}}$ a consistent approximation of $\hat{G}$ defined by (7), such an inequality writes

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

In the flat region, i.e. $\partial_x z = 0$, the system (1) is nothing but a conservative hyperbolic system and several schemes verifying the discrete entropy inequality (14) have already been introduced (for a nonhexaustive list, see [23, 24, 28–31, 33]). In the present work, we focus on Godunov-type scheme based on an approximate Riemann solver made of two intermediate constant states (see [29] for details). Considering the well-known Euler equations, the fully discrete entropy inequality (14) can be established for the two intermediate state HLLC scheme [45] or equivalently, for the Suliciu relaxation schemes [6, 12, 17]. Let us underline that Section b.ii in [29] presents an alternative to derive an entropy satisfying Godunov-type scheme with two intermediate constant states.

For a generic function $z$ (i.e. $\partial_x z \neq 0$), the design of a well-balanced scheme that satisfies a discrete entropy inequality (14) turns out to be more challenging. Both properties are satisfied by some Godunov-type schemes [1, 2, 5, 16, 27, 36, 47] or a relaxation well-balanced for the lake at rest scheme [12], but all of these schemes need to solve a set of non-linear equations at each cell of the mesh and for each time iteration. For instance, in [12], we have to solve a third-order polynomial equation to get an entropy preserving and well-balanced lake at rest scheme. Similarly in [7], a fifth-order polynomial equation must be solved to obtain a well-balanced scheme for the moving steady states which satisfies an entropy inequality. Moreover, the fully discrete entropy estimates may contain an error term (for instance, see [8]) in the form $O(\Delta x^2)$ such that

$$ \frac{\hat{\eta}(\hat{w}_i^{n+1}) - \hat{\eta}(\hat{w}_i^n)}{\Delta t} + \frac{\hat{G}_{i+\frac{1}{2}} - \hat{G}_{i-\frac{1}{2}}}{\Delta x} \le O(\Delta x^2). $$

The error term $O(\Delta x^2)$ does not occur in the works of [19, 32] but the proposed well-balanced schemes satisfied a global version of the entropy inequality (14) that writes $\sum_{i \in \mathbb{Z}} \hat{\eta}(\hat{w}_i^N) \Delta x \le \hat{\eta}(\hat{w}_i^0) \Delta x$ where $t^N$ and $t^0$ denote the final and the initial time of the simulation.
In this work, we propose to design numerical schemes to approximate the weak solutions of system (1) which satisfies the discrete entropy inequality (14). The main originality is the introduction of conditions in the scheme design that enforce the discrete entropy inequality (14). As a consequence, the entropy stability is included in the scheme definition. Therefore, it is very easy to establish it whereas in general the proof of the discrete entropy stability is a very difficult task. For instance, the HLLC scheme [45] was defined in 1994 and the proof of its entropy stability was only given in 2004 [12].

In addition, our schemes also have to satisfy the well-balanced property as soon as the topography function $z$ is non-constant. More precisely, we consider schemes written under the following form:

$$w_i^{n+1} = w_i^n - \frac{\Delta t}{\Delta x} \left( \hat{\mathcal{F}}(\hat{w}_i^n, \hat{w}_{i+1}^n) - \hat{\mathcal{F}}(\hat{w}_{i-1}^n, \hat{w}_i^n) \right) + \frac{\Delta t}{2} \left( \hat{S}_{i+\frac{1}{2}}^n + \hat{S}_{i-\frac{1}{2}}^n \right), \quad \forall i \in \mathbb{Z} \quad (15)$$

where $\hat{\mathcal{F}}: (\hat{\Omega})^2 \to \mathbb{R}^2$ stands for the numerical flux function and $\hat{S}_{i+\frac{1}{2}}^n \in \mathbb{R}^2$ denotes a consistent approximation of the source term $(0, -gh\partial_x z)^T$ at the interface $x_{i+\frac{1}{2}}$. To address such an issue, the paper is organized as follows. In Section 2, we adopt Godunov-type schemes to design (15). We give sufficient conditions to obtain the weak consistency, the discrete entropy inequality (14) and the well-balanced property. The first two sufficient conditions are used in Section 3 to define an approximate Riemann solver made of two intermediate states. This approximate Riemann solver is governed by an under-determined system that ensures the consistency and the discrete entropy stability (14). In Section 4, the proposed scheme derivation is exemplified to obtain consistent, entropy satisfying schemes in the case of $z = \text{cst}$. Then Sections 5 and 6 propose other closures in order to define consistent, entropy satisfying, well-balanced schemes for both lake at rest (10) and moving equilibrium (9). In Section 7, numerical tests are carried out to illustrate our numerical schemes.

## 2. GODUNOV-TYPE SCHEME DERIVATION

In this section, we design the finite volume scheme (15) adopting Godunov-type approaches [29]. To address such a derivation, we introduce a given approximate Riemann solver $\tilde{w}(\frac{x}{t}; \hat{w}_L, \hat{w}_R) \in \hat{\Omega}$, with $\hat{w}_L = (w_L, z_L)^T$ and $\hat{w}_R = (w_R, z_R)^T$ in $\hat{\Omega}$, to mimic the behavior of the exact Riemann solution associated to (2) with an initial data given by $w_0(x) = w_L$ if $x < 0$ and $w_0(x) = w_R$ if $x > 0$. This approximate Riemann solver will be enforced to satisfy some consistency conditions given later on. In fact, the sequel of the present paper is devoted to suitably designing $\tilde{w}$.

Next, equipped with an approximate Riemann solver, we are able to derive the Godunov-type scheme as follows:

$$w_i^{n+1} = \frac{1}{\Delta x} \int_0^{\frac{\Delta x}{2}} \tilde{w} \left( \frac{x}{\Delta t}, \hat{w}_{i-1}^n, \hat{w}_i^n \right) dx + \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^0 \tilde{w} \left( \frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n \right) dx, \quad \forall i \in \mathbb{Z}. \quad (16)$$

We emphasize that the time step $\Delta t$ is now assumed to be small enough such that two successive approximate Riemann solvers do not interact; namely $\tilde{w}(\xi, \hat{w}_{i-1}^n, \hat{w}_i^n) = \tilde{w}(-\xi, \hat{w}_i^n, \hat{w}_{i+1}^n)$ with $\xi$ in a neighborhood of $\frac{\Delta x}{2\Delta t}$. This time step restriction is nothing but a CFL-like condition.

The following lemma precisely states the condition to be satisfied by an approximate Riemann solver $\tilde{w}: \mathbb{R} \times (\hat{\Omega})^2 \to \hat{\Omega}$ so that the numerical scheme (15) and the Godunov-type scheme (16) coincide.

**Lemma 2.1** (Reformulation of the numerical scheme (15)). Let us consider $\hat{w} = (w, z)^T$ in $\hat{\Omega}$, a function $\hat{s}: (\hat{\Omega})^2 \to \mathbb{R}$ such that

$$\hat{s}(\hat{w}, \hat{w}) = -gh, \quad \forall \hat{w} \in \hat{\Omega}, \quad (17)$$

and $(\hat{S}_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}}$ a sequence in $\mathbb{R}^2$ that writes

$$\hat{S}_{i+\frac{1}{2}}^n = \begin{pmatrix} 0 \\ \hat{s}(\hat{w}_i^n, \hat{w}_{i+1}^n) \frac{z_{i+1}-z_i}{\Delta x} \end{pmatrix}, \quad \forall i \in \mathbb{Z}, \quad (18)$$
where $(z_i)_{i \in \mathbb{Z}}$ is defined by (11). Let us also consider a Godunov-type scheme (16) defined by an approximate Riemann solver $\tilde{w} : \mathbb{R} \times (\hat{\Omega})^2 \to \Omega$ that satisfies $\tilde{w}(\cdot, \hat{w}, \hat{w}) = w$ and the following integral consistency condition:

$$ \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n\right) dx = \frac{w_i^n + w_{i+1}^n}{2} - \frac{\Delta t}{\Delta x} (f(w_{i+1}^n) - f(w_i^n)) + \Delta t \hat{S}_{i+\frac{1}{2}}^n. \quad (19) $$

Assume $\Delta t$ small enough such that the non-interacting CFL-like condition holds. The Godunov-type scheme (16) then reformulates as the numerical scheme (15) with the numerical flux function $\hat{\mathcal{F}} : (\hat{\Omega})^2 \to \mathbb{R}^2$ given by

$$ \hat{\mathcal{F}}(\hat{w}_i^n, \hat{w}_{i+1}^n) = \frac{f(w_{i+1}^n) + f(w_i^n)}{2} - \frac{\Delta x}{4\Delta t}(w_{i+1}^n - w_i^n) + \frac{1}{2\Delta t} \int_0^{\frac{\Delta x}{2}} \tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n\right) dx \\ - \frac{1}{2\Delta t} \int_{-\frac{\Delta x}{2}}^0 \tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n\right) dx. \quad (20) $$

*Proof.* Since the non-interacting CFL-like condition holds, the Godunov-type scheme (16) can be rewritten as follows:

$$ w_i^{n+1} = \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{0} \tilde{w}\left(\frac{x}{\Delta t}, w_i^n, w_{i+1}^n\right) dx + \frac{1}{2\Delta x} \int_{0}^{\frac{\Delta x}{2}} \tilde{w}\left(\frac{x}{\Delta t}, w_{i-1}^n, w_i^n\right) dx \\ + \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \tilde{w}\left(\frac{x}{\Delta t}, w_{i-1}^n, w_i^n\right) dx + \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \tilde{w}\left(\frac{x}{\Delta t}, w_i^n, w_{i+1}^n\right) dx \\ - \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{0} \tilde{w}\left(\frac{x}{\Delta t}, w_{i-1}^n, w_i^n\right) dx - \frac{1}{2\Delta x} \int_{0}^{\frac{\Delta x}{2}} \tilde{w}\left(\frac{x}{\Delta t}, w_i^n, w_{i+1}^n\right) dx. $$

Using the integral consistency relation (19), the above equality rewrites

$$ \begin{aligned} w_i^{n+1} = & \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{0} \tilde{w}\left(\frac{x}{\Delta t}, w_i^n, w_{i+1}^n\right) dx + \frac{1}{2\Delta x} \int_{0}^{\frac{\Delta x}{2}} \tilde{w}\left(\frac{x}{\Delta t}, w_{i-1}^n, w_i^n\right) dx \\ & + \frac{1}{2} \left( \frac{w_{i-1}^n + w_i^n}{2} - \frac{\Delta t}{\Delta x} (f(w_i^n) - f(w_{i-1}^n)) + \Delta t \hat{S}_{i-\frac{1}{2}} \right) \\ & + \frac{1}{2} \left( \frac{w_i^n + w_{i+1}^n}{2} - \frac{\Delta t}{\Delta x} (f(w_{i+1}^n) - f(w_i^n)) + \Delta t \hat{S}_{i+\frac{1}{2}} \right) \\ & - \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{0} \tilde{w}\left(\frac{x}{\Delta t}, w_{i-1}^n, w_i^n\right) dx - \frac{1}{2\Delta x} \int_{0}^{\frac{\Delta x}{2}} \tilde{w}\left(\frac{x}{\Delta t}, w_i^n, w_{i+1}^n\right) dx. \end{aligned} $$

We then easily recover the finite volume scheme (15) where the numerical flux function $\hat{\mathcal{F}}$ is given by (20) and the source term $(\hat{S}_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}}$ by (17) and (18). The proof is thus completed. $\square$

Now, we state the main conditions to be satisfied by the approximate Riemann solver and the source term approximation in order to get the required entropy stability and the well-balanced property. To address such an issue, we first give the adopted discretization of the steady solutions and the properties to be satisfied by a well-balanced scheme.

**Definition 2.1** (Local equilibrium and well-balanced scheme for the shallow water equations (1)). Consider the shallow water equations (1) endowed with the notation

$$ \hat{B}(\hat{w}) = \frac{u^2}{2} + g(h+z), \quad \forall \hat{w} \in \hat{\Omega}. \quad (21) $$

At time $t^n$, let us also consider $\hat{w}_i^n$ and $\hat{w}_{i+1}^n$ two states in $\hat{\Omega}$ respectively on the left and on the right of the interface $x_{i+\frac{1}{2}}$.
(i) The two states $\hat{w}_i^n$ and $\hat{w}_{i+1}^n$ define a local lake at rest equilibrium at the interface $x_{i+\frac{1}{2}}$ if

$$u_i^n = u_{i+1}^n = 0, \quad h_i^n + z_i = h_{i+1}^n + z_{i+1}. \qquad (22)$$

(ii) The two states $\hat{w}_i^n$ and $\hat{w}_{i+1}^n$ define a local moving equilibrium at the interface $x_{i+\frac{1}{2}}$ if

$$h_i^n u_i^n = h_{i+1}^n u_{i+1}^n, \quad \hat{B}(\hat{w}_i^n) = \hat{B}(\hat{w}_{i+1}^n). \qquad (23)$$

(iii) A numerical scheme is well-balanced for the lake at rest (10) if $w_i^{n+1} = w_i^n$ for all $i \in \mathbb{Z}$ as soon as we have $\forall i \in \mathbb{Z}$

$$u_i^n = u_{i+1}^n = 0 \quad \text{and} \quad h_i^n + z_i = h_{i+1}^n + z_{i+1}.$$

(iv) A numerical scheme is fully well-balanced for the moving equilibrium (9) if $w_i^{n+1} = w_i^n$ for all $i \in \mathbb{Z}$ as soon as we have $\forall i \in \mathbb{Z}$

$$h_i^n u_i^n = h_{i+1}^n u_{i+1}^n \quad \text{and} \quad \hat{B}(\hat{w}_i^n) = \hat{B}(\hat{w}_{i+1}^n).$$

The conditions to be imposed on the approximate Riemann solver $\tilde{w}$ to get both the expected discrete entropy inequality and the well-balanced property are now stated. From now on, we underline that the well-balanced conditions are well-known (for instance, see [7, 10, 40, 41]) and they are here recalled for the sake of completeness. Moreover, although the function $w \to \eta(w)$, given by (5), is a convex function, it is worth noticing that the entropy function $\hat{w} \to \hat{\eta}(\hat{w})$ is not a convex function. As a consequence, the well-known integral entropy consistency introduced by Harten, Lax and van Leer in [29], cannot be applied here. The main originality in the present work is the introduction of an adapted integral entropy consistency based on the convex function $w \to \eta(w)$. In the following statement, this new integral entropy consistency is presented and shown to give the expected discrete entropy inequality (14).

**Lemma 2.2** (Consistent, well-balanced, entropy-stable Godunov-type scheme for the shallow water equations (1)). Consider the functions $\eta, G: \Omega \to \mathbb{R}$ defined by (5). Let us also consider $(\hat{S}_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}}$, which satisfies the consistency definitions (17) and (18), and an approximate Riemann solver $\tilde{w}: \mathbb{R} \times (\hat{\Omega})^2 \to \Omega$ verifying the integral consistency relation (19). Let us denote $\mathcal{F}_{i+\frac{1}{2}}^h$ the numerical flux of the variable $h$ given by the first component of $\hat{\mathcal{F}}(\hat{w}_i^n, \hat{w}_{i+1}^n)$ defined by (20). Assume that the CFL condition of non interaction holds. The following statements hold:

(i) The Godunov-type scheme (16), or equivalently the scheme (15)-(18)-(20), is consistent with the shallow water system (1).

(ii) If the approximate Riemann solver satisfies the inequality

$$\frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \eta \left( \tilde{w} \left( \frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n \right) \right) dx \le \frac{\eta(w_i^n) + \eta(w_{i+1}^n)}{2} - \frac{\Delta t}{\Delta x} (G(w_{i+1}^n) - G(w_i^n)) \\ - g \frac{\Delta t}{\Delta x} (z_{i+1} - z_i) \mathcal{F}_{i+\frac{1}{2}}^h, \quad \forall i \in \mathbb{Z}, \qquad (24)$$

then the Godunov-type scheme (16), or equivalently the scheme (15)-(18)-(20), verifies a discrete entropy inequality (14) where

$$\begin{aligned} \hat{\mathcal{G}}_{i+\frac{1}{2}} = \frac{G(w_{i+1}^n) + G(w_i^n)}{2} &+ g \mathcal{F}_{i+\frac{1}{2}}^h \frac{z_{i+1} + z_i}{2} - \frac{\Delta x}{4\Delta t} (\eta(w_{i+1}^n) - \eta(w_i^n)) \\ &+ \frac{1}{2\Delta t} \int_0^{\frac{\Delta x}{2}} \eta \left( \tilde{w} \left( \frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n \right) \right) dx - \frac{1}{2\Delta t} \int_{-\frac{\Delta x}{2}}^0 \eta \left( \tilde{w} \left( \frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n \right) \right) dx. \end{aligned} \qquad (25)$$
(iii) If the approximate Riemann solver satisfies

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

as soon as the sequence $(\hat{w}_i^n)_{i \in \mathbb{Z}}$ verifies, at each interface of the mesh, a local moving equilibrium (23) (resp. a local lake at rest equilibrium (22)) then the Godunov-type scheme (16), or equivalently the scheme (15)-(18)-(20), is well-balanced for the moving equilibrium (resp. for the lake at rest).

*Proof.* Concerning the first statement (i), according to the definition of $\hat{S}_{i+\frac{1}{2}}^n$ given by (17) and (18) and since the function $z$ is smooth enough, $\hat{S}_{i+\frac{1}{2}}^n$ is immediately consistent with $(0, -gh\partial_x z)^T$. As a consequence, after [29] (see also [21]), an approximate Riemann solver that verifies the consistency integral relation (19) defines a consistent scheme. The statement (i) is then proved.

Concerning the discrete entropy inequality given by (ii), let us consider the convex function $w \mapsto \eta(w)$ defined by (5). Hence, the Jensen inequality applied to the Godunov-type scheme (16) gives

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

Using the inequality (24) in the above estimate, we have

$$ \eta(w_i^{n+1}) \le \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{0} \eta\left(\tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n\right)\right) dx + \frac{1}{2\Delta x} \int_{0}^{\frac{\Delta x}{2}} \eta\left(\tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_{i-1}^n, \hat{w}_i^n\right)\right) dx \\ + \frac{1}{2} \left( \frac{\eta(w_{i-1}^n) + \eta(w_i^n)}{2} - \frac{\Delta t}{\Delta x} (G(w_i^n) - G(w_{i-1}^n)) - \frac{\Delta t}{\Delta x} g(z_i - z_{i-1}) \mathcal{F}_{i-\frac{1}{2}}^h \right) \\ + \frac{1}{2} \left( \frac{\eta(w_i^n) + \eta(w_{i+1}^n)}{2} - \frac{\Delta t}{\Delta x} (G(w_{i+1}^n) - G(w_i^n)) - \frac{\Delta t}{\Delta x} g(z_{i+1} - z_i) \mathcal{F}_{i+\frac{1}{2}}^h \right) \\ - \frac{1}{2\Delta x} \int_{-\frac{\Delta x}{2}}^{0} \eta\left(\tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_{i-1}^n, \hat{w}_i^n\right)\right) dx - \frac{1}{2\Delta x} \int_{0}^{\frac{\Delta x}{2}} \eta\left(\tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n\right)\right) dx. \qquad (27) $$

Now, since there is no contribution of the source term $S(w, z)$ in the first equation of the shallow water system (1), the numerical scheme for the variable $h$ writes

$$ h_i^{n+1} = h_i^n - \frac{\Delta t}{\Delta x} \left( \mathcal{F}_{i+\frac{1}{2}}^h - \mathcal{F}_{i-\frac{1}{2}}^h \right), \quad \forall i \in \mathbb{Z}. \qquad (28) $$

Multiplying the above equation (28) by $gz_i$ and adding to the inequality (27), we then obtain an inequality in the form

$$ \eta(w_i^{n+1}) + gh_i^{n+1}z_i \le \eta(w_i^n) + gh_i^n z_i - \frac{\Delta t}{\Delta x} (\hat{G}_{i+\frac{1}{2}} - \hat{G}_{i-\frac{1}{2}}), \quad \forall i \in \mathbb{Z}, $$

with $\hat{G}_{i+\frac{1}{2}}$ given by (25).
Since we have $\hat{\eta}(\hat{w}) = \eta(w) + ghz$ and since the function $z$ does not depend on the time, this last inequality rewrites in the expected form of (14). Before concluding the proof of the statement (ii), we have to show that the entropy numerical flux $\hat{G}_{i+\frac{1}{2}}$ given by (25) is consistent with the entropy flux function $\hat{G}$ defined by (7).

For a given constant state $w \in \Omega$, we set $\hat{w} = (w, z)^T$. Next, considering the consistency equality $\tilde{w}(\cdot, \hat{w}, \hat{w}) = w$ and the consistency of the water height numerical flux $\mathcal{F}_{i+\frac{1}{2}}^h$, we get

$$ 
\begin{aligned}
\hat{G}(\hat{w}, \hat{w}) &= \frac{G(w) + G(w)}{2} + g\mathcal{F}^h(w, w)z - \frac{\Delta x}{4\Delta t}(\eta(w) - \eta(w)) \\
&\quad + \frac{1}{2\Delta t} \int_0^{\frac{\Delta x}{2}} \eta(\tilde{w}(\frac{x}{\Delta t}, \hat{w}, \hat{w})) dx - \frac{1}{2\Delta t} \int_{-\frac{\Delta x}{2}}^0 \eta(\tilde{w}(\frac{x}{\Delta t}, \hat{w}, \hat{w})) dx, \\
&= G(w) + g(hu)z + \frac{1}{2\Delta t} \int_0^{\frac{\Delta x}{2}} \eta(w) dx - \frac{1}{2\Delta t} \int_{-\frac{\Delta x}{2}}^0 \eta(w) dx, \\
&= G(w) + g(hu)z.
\end{aligned}
$$

Since $\hat{G}(\hat{w}) = G(w) + g(hu)z$, the above equality achieves to prove the statement (ii).

Finally, concerning the statement (iii), let us assume that the sequence $(\hat{w}_i^n)_{i \in \mathbb{Z}}$ is such that a local equilibrium is maintained at each interface of the mesh. In this case, using the condition (26) in a Godunov-type scheme (16), we have

$$ 
\begin{aligned}
w_i^{n+1} &= \frac{1}{\Delta x} \int_0^{\frac{\Delta x}{2}} \tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_{i-1}^n, \hat{w}_i^n\right) dx + \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^0 \tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n\right) 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{aligned}
$$

According to Definition 2.1-(iii) or (iv), this last equality gives the required well-balanced property that concludes the proof. $\square$

### 3. CONSISTENT, ENTROPY-STABLE TWO STATES APPROXIMATE RIEMANN SOLVER FOR THE SHALLOW WATER EQUATIONS

It is now clear that the complete characterization of the Godunov-type scheme (16), or equivalently (15) where the numerical flux function is given by (20) and the source term approximation by (17) and (18), is obtained as soon as a suitable definition of both approximate Riemann solver $\tilde{w}(\xi, \hat{w}_L, \hat{w}_R)$ and source term approximation $\hat{s}(\hat{w}_L, \hat{w}_R)$ are designed. Of course $\tilde{w}$ and $\hat{s}$ must satisfy the conditions stated in Lemma 2.2 in order to get an entropy-stable and well-balanced scheme. At this level, we reformulate these conditions as a system to be solved by the intermediate states when adopting an approximate Riemann solver made of two intermediate states as follows:

$$ \tilde{w}(x/\Delta t, \hat{w}_L, \hat{w}_R) = \begin{cases} w_L, & \text{if } \frac{x}{\Delta t} \le -\lambda, \\ w_L^*, & \text{if } -\lambda < \frac{x}{\Delta t} \le 0, \\ w_R^*, & \text{if } 0 < \frac{x}{\Delta t} \le \lambda, \\ w_R, & \text{if } \lambda < \frac{x}{\Delta t}. \end{cases} \qquad (29) $$

For the sake of simplicity in the forthcoming derivations, we focus on an interface separating $(w_i^n, z_i) = (w_L, z_L) = \hat{w}_L$ and $(w_{i+1}^n, z_{i+1}) = (w_R, z_R) = \hat{w}_R$. Moreover, on this interface we set $\hat{s}_{LR} = \hat{s}(\hat{w}_L, \hat{w}_R)$ to be consistent with $-gh$.
Because of this particular choice of approximate Riemann solver, let us underline that the non-interacting CFL-like condition is satisfied as soon as

$$ \lambda \geq \max_{\alpha \in \{L,R\}} \left| u_{\alpha} \pm \sqrt{gh_{\alpha}} \right| \quad \text{and} \quad \frac{\lambda \Delta t}{\Delta x} \leq \frac{1}{2}. \qquad (30) $$

Now, we give sufficient relations to be satisfied by $w_L^*$ and $w_R^*$ so that the resulting approximate Riemann solver yields a consistent and entropy preserving Godunov-type scheme. At this level, it is essential to notice that the required discrete entropy inequality comes from the estimate (24) which also holds if the equality is imposed. The main idea here is to consider an equality in (24) to get an additional non-linear equation.

In order to exhibit these equations to be satisfied by $w_L^*$ and $w_R^*$, some notations are introduced. Indeed, for the sake of clarity, we set

$$ w^{\text{HLL}} = (h^{\text{HLL}}, (hu)^{\text{HLL}})^{\text{T}} = \frac{w_R + w_L}{2} - \frac{f(w_R) - f(w_L)}{2\lambda} \in \mathbb{R}^2, \qquad (31a) $$

$$ h^{\text{HLL}} \hat{u}^{\text{HLL}} = (hu)^{\text{HLL}} + \frac{\hat{s}_{LR}(z_R - z_L)}{2\lambda} \in \mathbb{R}, \qquad (31b) $$

$$ \hat{w}^{\text{HLL}} = \left( h^{\text{HLL}}, h^{\text{HLL}} \hat{u}^{\text{HLL}}, \frac{z_L + z_R}{2} \right)^{\text{T}} \in \mathbb{R}^3, \qquad (31c) $$

$$ \eta^{\text{HLL}} = \frac{\eta(w_R) + \eta(w_L)}{2} - \frac{G(w_R) - G(w_L)}{2\lambda} \in \mathbb{R}, \qquad (31d) $$

$$ \hat{\eta}^{\text{HLL}} = \frac{\hat{\eta}(\hat{w}_R) + \hat{\eta}(\hat{w}_L)}{2} - \frac{\hat{G}(\hat{w}_R) - \hat{G}(\hat{w}_L)}{2\lambda} \in \mathbb{R}. \qquad (31e) $$

**Lemma 3.1.** Consider $\hat{w}_L$ and $\hat{w}_R$ two states with $h_L > 0$ and $h_R > 0$, $\hat{s}_{LR}$ a consistent discretization of $-gh$ according to (17) and $\tilde{w}(\cdot, \hat{w}_L, \hat{w}_R) : \mathbb{R} \to \mathbb{R}^2$ an approximate Riemann solver in the form (29). Let us also consider the couples $(\eta, G)$ and $(\hat{\eta}, \hat{G})$ given by (5) and (7). Assume that the CFL condition (30) holds. The integral consistency condition (19) completed with the following integral entropy consistency condition:

$$ \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \eta \left( \tilde{w} \left( \frac{x}{\Delta t}, \hat{w}_L, \hat{w}_R \right) \right) dx = \frac{\eta(w_L) + \eta(w_R)}{2} - \frac{\Delta t}{\Delta x} (G(w_R) - G(w_L)) - g \frac{\Delta t}{\Delta x} (z_R - z_L) \mathcal{F}_{LR}^h, \qquad (32) $$

reformulates

$$ \frac{h_L^* + h_R^*}{2} = h^{\text{HLL}}, \qquad (33a) $$

$$ \frac{h_L^* u_L^* + h_R^* u_R^*}{2} = h^{\text{HLL}} \hat{u}^{\text{HLL}}, \qquad (33b) $$

$$ \frac{h_L^* h_R^*}{8h^{\text{HLL}}} (u_R^* - u_L^*)^2 + \frac{g}{8} (h_R^* - h_L^* + z_R - z_L)^2 = \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8} (z_R - z_L)^2. \qquad (33c) $$

Before proving the above statement, we emphasize the quantity $h^{\text{HLL}}\hat{u}^{\text{HLL}}$ given by (31b) depends on the free parameter $\hat{s}_{LR}$. This parameter has to satisfy the consistency condition described in (17) and some admissible choices will be given in the next sections devoted to the well-balanced property. Lemma 3.1 is now established.

*Proof.* First, we show that the system made of (19) and (32) coincides with the system (33). Since $\tilde{w}$ is given by (29), the left-hand side of the equation (19) writes

$$ \begin{aligned} \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \tilde{w} \left( \frac{x}{\Delta t}, \hat{w}_L, \hat{w}_R \right) dx &= \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{-\lambda \Delta t} w_L dx + \frac{1}{\Delta x} \int_{-\lambda \Delta t}^{0} w_L^* dx + \frac{1}{\Delta x} \int_{0}^{\lambda \Delta t} w_R^* dx + \frac{1}{\Delta x} \int_{\lambda \Delta t}^{\frac{\Delta x}{2}} w_R dx, \\ &= \left( -\lambda \frac{\Delta t}{\Delta x} + \frac{1}{2} \right) w_L + \lambda \frac{\Delta t}{\Delta x} w_L^* + \lambda \frac{\Delta t}{\Delta x} w_R^* + \left( \frac{1}{2} - \lambda \frac{\Delta t}{\Delta x} \right) w_R. \end{aligned} $$
As a consequence the equation (19) equivalently rewrites

$$ \left(-\lambda \frac{\Delta t}{\Delta x} + \frac{1}{2}\right) w_L + \lambda \frac{\Delta t}{\Delta x} w_L^* + \lambda \frac{\Delta t}{\Delta x} w_R^* + \left(\frac{1}{2} - \lambda \frac{\Delta t}{\Delta x}\right) w_R = \frac{w_L + w_R}{2} - \lambda \frac{\Delta t}{\Delta x} \frac{f(w_R) - f(w_L)}{\lambda} \\ + \lambda \frac{\Delta t}{\Delta x} \left(0, \frac{\hat{s}_{LR}(z_R - z_L)}{\lambda}\right)^\mathrm{T}. $$

Writing the above equation component by component, we easily obtain the equations (33a) and (33b).

Now, we have to recover the last equation (33c). By expanding the integral of the left-hand side in (32), we get

$$ \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \eta\left(\tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_L, \hat{w}_R\right)\right) dx \\ = \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{-\lambda \Delta t} \eta(w_L) dx + \frac{1}{\Delta x} \int_{-\lambda \Delta t}^{0} \eta(w_L^*) dx + \frac{1}{\Delta x} \int_{0}^{\lambda \Delta t} \eta(w_R^*) dx + \frac{1}{\Delta x} \int_{\lambda \Delta t}^{\frac{\Delta x}{2}} \eta(w_R) dx, \\ = \left(-\lambda \frac{\Delta t}{\Delta x} + \frac{1}{2}\right) \eta(w_L) + \lambda \frac{\Delta t}{\Delta x} \eta(w_L^*) + \lambda \frac{\Delta t}{\Delta x} \eta(w_R^*) + \left(\frac{1}{2} - \lambda \frac{\Delta t}{\Delta x}\right) \eta(w_R). $$

So that, the equation (32) equivalently writes

$$ \left(-\lambda \frac{\Delta t}{\Delta x} + \frac{1}{2}\right) \eta(w_L) + \lambda \frac{\Delta t}{\Delta x} \eta(w_L^*) + \lambda \frac{\Delta t}{\Delta x} \eta(w_R^*) + \left(\frac{1}{2} - \lambda \frac{\Delta t}{\Delta x}\right) \eta(w_R) \\ = \frac{\eta(w_L) + \eta(w_R)}{2} - \lambda \frac{\Delta t}{\Delta x} \frac{G(w_R) - G(w_L)}{\lambda} - g \lambda \frac{\Delta t}{\Delta x} \frac{z_R - z_L}{\lambda} \mathcal{F}_{LR}^h. $$

With $\eta^{\mathrm{HLL}}$ defined by (31d), we immediately reformulate (32) as follows:

$$ \frac{\eta(w_L^*) + \eta(w_R^*)}{2} = \eta^{\mathrm{HLL}} - g \frac{z_R - z_L}{2\lambda} \mathcal{F}_{LR}^h. \qquad (34) $$

We now expand the numerical flux $\mathcal{F}_{LR}^h$ given by (20). Adopting the notation $[X] = X_R - X_L$ for any quantity $X$, since the approximate Riemann solver is given by (29), the water height numerical flux function writes

$$ \mathcal{F}_{LR}^h = \frac{(hu)_R + (hu)_L}{2} - \frac{\lambda}{2}[h] + \frac{\lambda}{2}[h^*]. $$

As a consequence, and using the definition of $h^{\mathrm{HLL}}$ given by (31a), we have

$$ \begin{aligned} -\frac{[z]}{2\lambda} \mathcal{F}_{LR}^h &= -\frac{[z][h^*]}{4} + \frac{[z]}{4} \left( [h] - \frac{(hu)_L + (hu)_R}{\lambda} \right), \\ &= -\frac{[z][h^*]}{4} + \frac{z_R h_R + z_L h_L}{2} - \frac{z_L + z_R}{2} h^{\mathrm{HLL}} - \frac{[h uz]}{2\lambda}. \end{aligned} $$

Now, considering the above relation and $\hat{\eta}^{\mathrm{HLL}}$, given by (31e), expanding $\eta(w_L^*)$ and $\eta(w_R^*)$, we then deduce from (34)

$$ \frac{h_R^*(u_R^*)^2 + h_L^*(u_L^*)^2}{4} + g \frac{(h_L^*)^2 + (h_R^*)^2}{4} + \frac{g}{4}[z][h^*] = \hat{\eta}^{\mathrm{HLL}} - g \frac{z_L + z_R}{2} h^{\mathrm{HLL}}. \qquad (35) $$

Moreover, from the equation (33a) we directly obtain

$$ \frac{g}{8}(h_L^* + h_R^*)^2 = \frac{g}{2}(h^{\mathrm{HLL}})^2, $$
so that, subtracting the above equation to (35), we deduce

$$ \frac{h_R^*(u_R^*)^2 + h_L^*(u_L^*)^2}{4} + \frac{g}{8}[h^*]^2 + \frac{g}{4}[z][h^*] = \hat{\eta}^{\text{HLL}} - \frac{g}{2}(h^{\text{HLL}})^2 - g \frac{z_L + z_R}{2} h^{\text{HLL}}. \quad (36) $$

The next step consists in rewriting $h_L^*(u_L^*)^2+h_R^*(u_R^*)^2$ in order to show that the above equation (36) is equivalent to

$$ \frac{h_R^* h_L^*[u^*]^2}{8h^{\text{HLL}}} + \frac{g}{8}[h^*]^2 + \frac{g}{4}[h^*][z] = \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}). \quad (37) $$

According to the definition of $w^{\text{HLL}}$ given by (31a), a direct computation shows that $h^{\text{HLL}}$ writes

$$ h^{\text{HLL}} = \frac{h_R}{2} \left(1 - \frac{u_R}{\lambda}\right) + \frac{h_L}{2} \left(1 + \frac{u_L}{\lambda}\right). \quad (38) $$

Since the CFL condition (30) holds and since $h_L > 0$ and $h_R > 0$, the above equality leads to $h^{\text{HLL}} > 0$. As a consequence, and using the equations (33a) and (33b), the following computation holds:

$$ \begin{aligned} \frac{h_L^*(u_L^*)^2 + h_R^*(u_R^*)^2}{4} &= \frac{h_L^*(u_L^*)^2 + h_R^*(u_R^*)^2}{4} \frac{h_L^* + h_R^*}{2h^{\text{HLL}}}, \\ &= \frac{(h_L^* u_L^*)^2 + 2h_L^* u_L^* h_R^* u_R^* + (h_R^* u_R^*)^2 + h_R^* h_L^*[u^*]^2}{8h^{\text{HLL}}}, \\ &= \frac{1}{2h^{\text{HLL}}} \left( \frac{h_L^* u_L^* + h_R^* u_R^*}{2} \right)^2 + \frac{h_L^* h_R^*[u^*]^2}{8h^{\text{HLL}}}, \\ &= \frac{(h^{\text{HLL}} \hat{u}^{\text{HLL}})^2}{2h^{\text{HLL}}} + \frac{h_L^* h_R^*[u^*]^2}{8h^{\text{HLL}}}. \end{aligned} \quad (39) $$

Considering the above relation in (36), we eventually deduce that the equation (32) is equivalent to the equation (37). Adding $\frac{g[z]^2}{8}$ on both sides of (37), we deduce the expected equation (33c) and the proof is achieved. $\square$

In the above lemma we have imposed $h_L > 0$ and $h_R > 0$. To deal with the dry-wet transitions, we need a specific treatment [4,39,41]. Indeed, the quantity $\hat{u}^{\text{HLL}}$ given by (31b) is defined up to a multiplication by $h^{\text{HLL}}$. According to (38), since $h^{\text{HLL}}$ is null if and only if $h_L = 0$ and $h_R = 0$, the admitted convention in [4,39,41] requires to define $\hat{u}^{\text{HLL}}$ in $\mathbb{R}$ as follows:

$$ \hat{u}^{\text{HLL}} = \begin{cases} 0, & \text{if } h_L = 0 \text{ and } h_R = 0, \\ \frac{(hu)^{\text{HLL}} + \frac{\hat{s}_{LR}(z_R - z_L)}{2\lambda}}{h^{\text{HLL}}}, & \text{otherwise.} \end{cases} \quad (40) $$

However, from a practical point of view, the dry-wet transitions must be considered but they are not the main goal of this work. As a consequence, and for clarity, with $h_L > 0$ and $h_R > 0$ both definitions (31b) and (40) are equivalent.

Now, we emphasize that intermediate states $w_L^*$ and $w_R^*$, satisfying (19) and (32), because of Lemma 2.2 (i) and (ii), will produce a consistent entropy preserving Godunov-type scheme (16). At this level, we also notice that the system made of (19) and (32) governing $w_L^*$ and $w_R^*$ is under-determined since we have only three equations to define $h_L^*$, $h_R^*$, $u_L^*$, $u_R^*$, and the source term discretization $\hat{s}_{LR}$. The missing relation will be detailed in the next section according to the selected well-balanced property.

To conclude this section, we note a natural compatibility condition to be satisfied in (33c). Indeed, it is clear that a well-posed property of this quadratic equation requires the following inequality:

$$ \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g(z_R - z_L)^2}{8} \geq 0. \quad (41) $$

At this level, this above given inequality cannot be proved in its full generality. However, some properties are now provided.
**Lemma 3.2.** Consider $\hat{w}_L$ and $\hat{w}_R$ two states with $h_L > 0$ and $h_R > 0$. Let us also consider the couples $(\eta, G)$ and $(\hat{\eta}, \hat{G})$ given by (5) and (7), and the quantities $(\hat{w}^{\text{HLL}}, \hat{\eta}^{\text{HLL}})$ and $(w^{\text{HLL}}, \eta^{\text{HLL}})$, defined by (31). Assume that the CFL-like condition (30) holds. The following statements are verified:

(i) For all given quantity $\hat{s}_{LR}$, the following estimate is satisfied:

$$ 
\begin{aligned}
\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}(z_R - z_L)^2 &= \frac{h_R h_L (u_R - u_L)^2}{4(h_L + h_R)} + \frac{g}{8}(h_R - h_L + z_R - z_L)^2 \\
&\quad - \frac{(u_R - u_L)^2 h_R h_L (h_L u_R - h_R u_L)}{4\lambda(h_L + h_R)^2} - \frac{h_L u_L + h_R u_R}{4\lambda} \left( g(z_R - z_L) + \frac{\hat{s}_{LR}(z_R - z_L)}{h^{\text{HLL}}} \right) + \mathcal{O}\left(\frac{1}{\lambda^2}\right).
\end{aligned}
\quad (42) $$

(ii) Assume $z_L = z_R = 0$ and adopt the CFL-like restriction (30), then we have

$$ \eta^{\text{HLL}} - \eta(w^{\text{HLL}}) \geq 0. \quad (43) $$

*Proof.* Concerning (i), with the notation $[X] = X_R - X_L$ for any quantity X and arguing (31) to define the quantities $\hat{u}^{\text{HLL}}, \hat{\eta}^{\text{HLL}}, \hat{\eta}(\hat{w}^{\text{HLL}})$, we obtain

$$ 
\begin{aligned}
\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 &= \hat{\eta}^{\text{HLL}} - \frac{((hu)^{\text{HLL}} + \frac{\hat{s}_{LR}(z_R - z_L)}{2\lambda})^2}{2h^{\text{HLL}}} - \frac{g}{2}(h^{\text{HLL}})^2 - gh^{\text{HLL}} \frac{z_L + z_R}{2} + \frac{g}{8}[z]^2, \\
&= \eta^{\text{HLL}} + g \frac{h_L z_L + h_R z_R}{2} - g \frac{[huz]}{2\lambda} - \eta(w^{\text{HLL}}) - gh^{\text{HLL}} \frac{z_L + z_R}{2} \\
&\quad - \frac{(hu)^{\text{HLL}}}{h^{\text{HLL}}} \frac{\hat{s}_{LR}(z_R - z_L)}{2\lambda} - \frac{1}{8h^{\text{HLL}}} \left( \frac{\hat{s}_{LR}(z_R - z_L)}{\lambda} \right)^2 + \frac{g}{8}[z]^2, \quad (44) \\
&= \eta^{\text{HLL}} - \eta(w^{\text{HLL}}) + \frac{g}{4}[z][h] + \frac{g}{8}[z]^2 \\
&\quad - \frac{g}{4\lambda}[z](h_L u_L + h_R u_R) - \frac{(hu)^{\text{HLL}}}{h^{\text{HLL}}} \frac{\hat{s}_{LR}(z_R - z_L)}{2\lambda} - \frac{1}{8h^{\text{HLL}}} \left( \frac{\hat{s}_{LR}(z_R - z_L)}{\lambda} \right)^2.
\end{aligned}
$$

The right-hand side of the above equality is now expanded with respect to $\frac{1}{\lambda}$ in a neighborhood of 0. Developing $w^{\text{HLL}}$ and $\eta^{\text{HLL}}$, given by (31a) and (31d), we have

$$ \eta^{\text{HLL}} - \eta(w^{\text{HLL}}) = \frac{h_L u_L^2}{4h_L} + \frac{h_R u_R^2}{4h_R} + \frac{g(h_L^2 + h_R^2)}{4} - \frac{[G]}{2\lambda} - \frac{g}{8} \left(h_R + h_L - \frac{[hu]}{\lambda}\right)^2 - \frac{((hu)^{\text{HLL}})^2}{2h^{\text{HLL}}}. \quad (45) $$

Now, it is necessary to expand the quantity $\frac{((hu)^{\text{HLL}})^2}{(2h^{\text{HLL}})}$ in the above equality. Using the definition of $h^{\text{HLL}}$ and $(hu)^{\text{HLL}}$ given by (31a), we have

$$ (hu)^{\text{HLL}} = \frac{h_L u_L + h_R u_R}{2} - \frac{1}{2\lambda} \left[ hu^2 + \frac{gh^2}{2} \right], \quad (46a) $$

$$ \frac{1}{2h^{\text{HLL}}} = \frac{1}{(h_L + h_R) \left( 1 - \frac{[hu]}{\lambda(h_R + h_L)} \right)}. \quad (46b) $$

Considering the square of $(hu)^{\text{HLL}}$ in (46a) and interpreting the quantity $\frac{1}{(1 - \frac{[hu]}{\lambda(h_R + h_L)})}$ as a geometric series, we get

$$ ((hu)^{\text{HLL}})^2 = \left( \frac{h_L u_L + h_R u_R}{2} \right)^2 - \frac{1}{2\lambda} (h_L u_L + h_R u_R) \left[ hu^2 + \frac{gh^2}{2} \right] + \mathcal{O}\left( \frac{1}{\lambda^2} \right), \quad (47a) $$

$$ \frac{1}{1 - \frac{[hu]}{\lambda(h_R + h_L)}} = 1 + \frac{[hu]}{\lambda(h_L + h_R)} + \mathcal{O}\left(\frac{1}{\lambda^2}\right). \quad (47b) $$
As a consequence, using the two above equalities to develop the quantity $\frac{((hu)^{\text{HLL}})^2}{(2h^{\text{HLL}})}$ with respect to $\frac{1}{\lambda}$, from (45) we obtain

$$ 
\begin{aligned}
\eta^{\text{HLL}} - \eta(w^{\text{HLL}}) &= \frac{h_L u_L^2}{4h_L} + \frac{h_R u_R^2}{4h_R} - \frac{(h_L u_L + h_R u_R)^2}{4(h_L + h_R)} + \frac{g}{8}[h]^2 - \frac{[G]}{2\lambda} + \frac{g}{4}(h_L + h_R)\frac{[hu]}{\lambda} \\
&\quad - \frac{(h_L u_L + h_R u_R)^2 [hu]}{4\lambda(h_L + h_R)^2} + \frac{(h_L u_L + h_R u_R)[hu^2 + \frac{g}{2}h^2]}{2\lambda(h_L + h_R)} + \mathcal{O}\left(\frac{1}{\lambda^2}\right), \\
&= \frac{h_R h_L [u]^2}{4(h_L + h_R)} + \frac{g}{8}[h]^2 - \frac{[u]^2 h_R h_L (h_L u_R - h_R u_L)}{4\lambda(h_L + h_R)^2} + \mathcal{O}\left(\frac{1}{\lambda^2}\right).
\end{aligned}
$$

With the two equalities (46), it is also possible to develop $\frac{(hu)^{\text{HLL}}}{h^{\text{HLL}}}$ and $\frac{1}{h^{\text{HLL}}}$ in the equation (44) then using the above estimate, we eventually deduce

$$ 
\begin{aligned}
\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 &= \frac{h_R h_L [u]^2}{4(h_L + h_R)} + \frac{g}{8}[h+z]^2 - \frac{[u]^2 h_R h_L (h_L u_R - h_R u_L)}{4\lambda(h_L + h_R)^2} \\
&\quad - \frac{h_L u_L + h_R u_R}{4\lambda} \left( g[z] + \frac{\hat{s}_{LR}(z_R - z_L)}{h^{\text{HLL}}} \right) + \mathcal{O}\left(\frac{1}{\lambda^2}\right),
\end{aligned}
$$

to achieve the establishment of (i).

Concerning (ii), with $\hat{s}_{LR}|_{z_R=z_L} = 0$, because of (44) we directly get the following equality:

$$ \left( \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}(z_R - z_L)^2 \right) \Big|_{z_L=z_R} = \eta^{\text{HLL}} - \eta(w^{\text{HLL}}). $$

To conclude the proof of this result, we underline that the inequality (43) is just an easy consequence of the discrete entropy inequality satisfied by the well-known one intermediate-state HLL scheme given in [29]. $\square$

It is worth noticing that the statement (i) of the above lemma highlights the role of the parameter $\hat{s}_{LR}$ in the inequality (41). In fact the choice of this parameter has to be motivated on the one hand for the consistency but also, on the other hand, to enforce the inequality (41) in the case $u_L = u_R$ and $h_L + z_L = h_R + z_R$.

Now, this system (33) is adopted to define an approximate Riemann solver (29).

## 4. APPLICATION IN THE FLAT REGIONS

This section concerns the design of an approximate Riemann solver (29) defined by the system (33) but for $z = \text{cst}$. We have to deal with $z_L = z_R$ locally on an interface. As a consequence, the system (33) now reads

$$ \frac{h_L^* + h_R^*}{2} = h^{\text{HLL}}, \qquad (48a) $$

$$ \frac{h_L^* u_L^* + h_R^* u_R^*}{2} = h^{\text{HLL}} u^{\text{HLL}}, \qquad (48b) $$

$$ \frac{h_L^* h_R^*}{8h^{\text{HLL}}} (u_R^* - u_L^*)^2 + \frac{g}{8} (h_R^* - h_L^*)^2 = \eta^{\text{HLL}} - \eta(w^{\text{HLL}}), \qquad (48c) $$

where $w^{\text{HLL}}$ and $\eta^{\text{HLL}}$ are given by (31a) and (31d), and with

$$ u^{\text{HLL}} = \begin{cases} 0, & \text{if } h_L = 0 \text{ and } h_R = 0, \\ \frac{(hu)^{\text{HLL}}}{h^{\text{HLL}}}, & \text{otherwise.} \end{cases} $$
Since the system (48) is under-determined, it is completed with the continuity of $u$ in the intermediate states of the solver (29). Such continuity writes

$$u_L^* = u_R^*, \qquad (49)$$

which has already been proposed in [45] for the HLLC scheme, for instance.

**Lemma 4.1** (Entropy-satisfying approximate Riemann solver for the homogeneous shallow water equations). Let us consider $w_L$ and $w_R$ with $h_L > 0$ and $h_R > 0$. Let us define an approximate Riemann solver $\tilde{w}(\cdot, w_L, w_R): \mathbb{R} \to \mathbb{R}^2$ in the form (29). Assume the CFL-like condition (30) holds. If the intermediate states $w_L^*$ and $w_R^*$ of the approximate Riemann solver are defined by the system (48) and (49) then, there exists two solutions for $w_L^*$ and $w_R^*$, according to the notations (31), given by

$$u_L^* = u_R^* = u^{\text{HLL}}, \qquad (50a)$$

$$h_R^* = h^{\text{HLL}} \pm \sqrt{\frac{2}{g}(\eta^{\text{HLL}} - \eta(w^{\text{HLL}}))}, \qquad (50b)$$

$$h_L^* = h^{\text{HLL}} \mp \sqrt{\frac{2}{g}(\eta^{\text{HLL}} - \eta(w^{\text{HLL}}))}. \qquad (50c)$$

Moreover, a Godunov-type scheme (16) associated with one of the two above approximate Riemann solvers

(i) is consistent with the homogeneous shallow water equations given by (1) with $z = \text{cst}$,

(ii) satisfies a discrete entropy inequality (14)-(25).

In the equalities (50), both symbols $\pm$ and $\mp$ that mean $+$ or $-$ but they are self-dependent. If $\pm$ is positive (resp. negative) then $\mp$ is negative (resp. positive). This notation will be kept in the sequel.

Before we establish the above result, since $z = \text{cst}$, we emphasize that the discrete entropy inequality, given by (14)-(25), reads

$$\frac{1}{\Delta t} \left( \left( \eta(\hat{w}_i^{n+1}) + zh_i^{n+1} \right) - \left( \eta(\hat{w}_i^n) + zh_i^n \right) \right) + \frac{1}{\Delta x} \left( \left( \mathcal{G}_{i+\frac{1}{2}} + z\mathcal{F}_{i+\frac{1}{2}}^h \right) - \left( \mathcal{G}_{i-\frac{1}{2}} + z\mathcal{F}_{i-\frac{1}{2}}^h \right) \right) \le 0,$$

where

$$\mathcal{G}_{i+\frac{1}{2}} = \frac{G(w_{i+1}^n) + G(w_i^n)}{2} - \frac{\Delta x}{4\Delta t} (\eta(w_{i+1}^n) - \eta(w_i^n)) \\ + \frac{1}{2\Delta t} \int_0^{\frac{\Delta x}{2}} \eta\left(\tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n\right)\right) dx - \frac{1}{2\Delta t} \int_{-\frac{\Delta x}{2}}^0 \eta\left(\tilde{w}\left(\frac{x}{\Delta t}, \hat{w}_i^n, \hat{w}_{i+1}^n\right)\right) dx. \qquad (51)$$

Since we have

$$h_i^{n+1} = h_i^n - \frac{\Delta t}{\Delta x} \left( \mathcal{F}_{i+\frac{1}{2}}^h - \mathcal{F}_{i-\frac{1}{2}}^h \right),$$

then we recover the usual discrete entropy inequality

$$\frac{1}{\Delta t} (\eta(\hat{w}_i^{n+1}) - \eta(\hat{w}_i^n)) + \frac{1}{\Delta x} (\mathcal{G}_{i+\frac{1}{2}} - \mathcal{G}_{i-\frac{1}{2}}) \le 0. \qquad (52)$$

*Proof of Lemma 4.1.* First, invoking the notation (31) we solve the system (48) and (49). Since $u_L^* = u_R^* = u^*$, the equation (48b) associated to (48a) gives

$$h^{\text{HLL}} u^{\text{HLL}} = \frac{h_L^* u_L^* + h_R^* u_R^*}{2}, \\ = \frac{h_L^* + h_R^*}{2} u^*,$$
$$= h^{\text{HLL}} u^* \qquad (53)$$

With $h_L > 0$ and $h_R > 0$ then, under the CFL-like condition (30), $h^{\text{HLL}} > 0$ so that (50a) is obtained.

Now, we prove the expressions of $h_L^*$ and $h_R^*$ given by (50b) and (50c). Since the equality $u_L^* = u_R^*$ holds, the equation (48c) reads

$$\frac{g}{8}(h_R^* - h_L^*)^2 = \eta^{\text{HLL}} - \eta(w^{\text{HLL}}).$$

According to Lemma 3.2, the inequality $\eta^{\text{HLL}} - \eta(w^{\text{HLL}}) \ge 0$ is ensured by the CFL-like condition (30). As a consequence, we get

$$h_R^* - h_L^* = \pm \sqrt{\frac{8(\eta^{\text{HLL}} - \eta(w^{\text{HLL}}))}{g}}.$$

The above equation with (48a) immediately gives (50b) and (50c).

Next, we note that $\tilde{w}(\cdot, w, w) = w$ for all $w \in \Omega$. Since the approximate Riemann solver satisfies the integral consistency, according to Lemma 3.1 and relations (33) reformulated here by (48a) and (48b), the resulting Godunov-type scheme (16) is consistent.

Eventually, the discrete entropy stability (ii) is a direct consequence of the equation (48c), of Lemma 2.2 and of Lemma 3.1, which concludes the proof. □

Unfortunately, adopting (50b) and (50c) to define $h_L^*$ and $h_R^*$, it is not possible to satisfy $h_L^* \ge 0$ and $h_R^* \ge 0$ for all $w_L$ and $w_R$ in $\Omega$. To enforce the required positivity of the intermediate water heights, we adopt a conservative cut-off introduced in [4, 9, 41].

**Lemma 4.2** (Robust and entropy-satisfying approximate Riemann solver for the homogeneous shallow water equations). Consider $w_L$ and $w_R$ in $\Omega$. Consider a Godunov-type scheme (16) that approximates the solution of the homogeneous shallow water equations and defined by an approximate Riemann solver $\tilde{w}(\cdot, w_L, w_R) : \mathbb{R} \to \mathbb{R}^2$ in the form (29). Assume the CFL-like condition (30) holds. Let us also consider the quantities $w^{\text{HLL}}$ and $\eta^{\text{HLL}}$ given by (31a) and (31d), and let us denote

$$\tilde{h}_R^* = h^{\text{HLL}} \pm \sqrt{\frac{2}{g}(\eta^{\text{HLL}} - \eta(w^{\text{HLL}}))}, \quad \tilde{h}_L^* = h^{\text{HLL}} \mp \sqrt{\frac{2}{g}(\eta^{\text{HLL}} - \eta(w^{\text{HLL}}))}. \qquad (54)$$

If the intermediate states $w_L^*$ and $w_R^*$ of the approximate Riemann solver write

$$\begin{align}
u_L^* &= u_R^* = u^{\text{HLL}}, \nonumber \\
h_R^* &= \min\left(\max(\tilde{h}_R^*, 0), 2h^{\text{HLL}}\right), \nonumber \\
h_L^* &= \min\left(\max(\tilde{h}_L^*, 0), 2h^{\text{HLL}}\right), \tag{55}
\end{align}$$

then we have

(i) The Godunov-type scheme (16), or equivalently the scheme (15)-(18)-(20), associated to such an approximate Riemann solver is consistent with the homogeneous shallow water equations given by (1) with $z = \text{cst}$ and it preserves the convex set $\Omega$.

(ii) If $h_L^* > 0$ and $h_R^* > 0$ then the Godunov-type scheme (16), or equivalently the scheme (15)-(18)-(20), associated to such an approximate Riemann solver satisfies the discrete entropy inequality (51) and (52).

*Proof.* First, if $\tilde{h}_L^*$ and $\tilde{h}_R^*$ given by (54) are both positive then (55) coincides with (50). In this case, Lemma 4.1 applies and the required consistency and entropy stability are recovered. Moreover, since $h_L^* > 0$ and $h_R^* > 0$, then by the definition of the Godunov-type scheme (16) we immediately get $h_i^{n+1} > 0$ for all $i \in \mathbb{Z}$.
Next, if one of $\tilde{h}_L^*$ or $\tilde{h}_R^*$ is non-positive, say $\tilde{h}_L^* \le 0$, then by (55) we get $h_L^* = 0$ and $h_R^* = 2h^{\text{HLL}}$. As a consequence, the consistency condition (48a) is satisfied and the scheme remains consistent. Once again, since $h_L^* \ge 0$ and $h_R^* \ge 0$, the associated Godunov-type scheme also preserves the positivity of the updated approximate water height $h_i^{n+1} \ge 0$ for all $i \in \mathbb{Z}$. The proof is thus complete. $\square$

We emphasize that the positiveness procedure $\min(\max(\cdot, \cdot), \cdot)$ imposed in (55), ensures the robustness in the regions near the dry areas. This procedure still ensures the consistency relation (48a) when $\tilde{h}_L^* \le 0$ or $\tilde{h}_R^* \le 0$ but it does not necessarily ensure the entropy condition (48c). As a consequence, the discrete entropy stability may be locally lost in the dry-wet transition regions.

## 5. A WELL-BALANCED ENTROPY-STABLE NUMERICAL SCHEME FOR THE LAKE AT REST

In this section, we consider the shallow water equations (1) for a given arbitrary smooth function $z: \mathbb{R} \to \mathbb{R}$. We aim to design a consistent, entropy-satisfying and well-balanced scheme for the lake at rest (10). In this regard, we propose to define an approximate Riemann solver (29) with the under-determined system (33) completed as follows:

$$u_L^* = u_R^*, \qquad (56a)$$

$$\hat{s}_{LR} = \hat{s}_{LR}^{\text{WBAR}} = -gh^{\text{HLL}} - \sqrt{\frac{g}{h^{\text{HLL}}} \frac{h_L u_L + h_R u_R}{h^{\text{HLL}}}} (z_R - z_L), \qquad (56b)$$

where the notation $\hat{s}_{LR}^{\text{WBAR}}$ means Well-Balanced At Rest. This notation is defined in order to distinguish $\hat{s}_{LR}^{\text{WBAR}}$ given by (56b) from a generic discretization $\hat{s}_{LR}$.

According to Lemma 2.2, the under-determined system (33) will be shown to be a sufficient condition for the consistency and the discrete entropy inequality (14) for the couple $(\hat{\eta}, \hat{G})$ defined by (7). The closure (56) is sufficient to ensure the well-balanced property for the lake at rest (10). Before proving this property, we show the existence of the solutions of the system made of (33) and (56).

**Lemma 5.1.** Consider $\hat{w}_L$ and $\hat{w}_R$ two states with $h_L > 0$ and $h_R > 0$ and consider $\tilde{w}(\cdot, \hat{w}_L, \hat{w}_R) : \mathbb{R} \to \mathbb{R}^2$ an approximate Riemann solver in the form (29). Assume the CFL-like condition (30) holds and $\lambda > 0$ large enough. The intermediate states, $w_L^*$ and $w_R^*$, solutions of the system (33) and (56), are given by

$$u_L^* = u_R^* = \hat{u}^{\text{HLL}}, \qquad (57a)$$

$$h_R^* = h^{\text{HLL}} - \frac{(z_R - z_L) \pm \sqrt{8 \frac{(\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL})))}{g} + (z_R - z_L)^2}}{2}, \qquad (57b)$$

$$h_L^* = h^{\text{HLL}} + \frac{(z_R - z_L) \pm \sqrt{8 \frac{(\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL})))}{g} + (z_R - z_L)^2}}{2}, \qquad (57c)$$

where

$$8 \frac{(\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}))}{g} + (z_R - z_L)^2 \ge 0. \qquad (58)$$

A Godunov-type scheme (16) defined by such an approximate Riemann solver

(i) is consistent with the shallow water equations (1),

(ii) satisfies a discrete entropy inequality (14) for the couple $(\hat{\eta}, \hat{G})$ defined by (7), where the numerical entropy flux function $\hat{G}_{i+\frac{1}{2}}$ of this inequality is given by (25),

(iii) is well-balanced for the lake at rest (10).
Proof. First, we solve the system made of (33) and (56). As $u_L^* = u_R^*$, let us denote $u^* = u_L^* = u_R^*$. The equation (33b) associated to (33a) gives

$$ 
\begin{aligned}
h^{\text{HLL}} \hat{u}^{\text{HLL}} &= \frac{h_L^* u_L^* + h_R^* u_R^*}{2}, \\
&= \frac{h_L^* + h_R^*}{2} u^*, \\
&= h^{\text{HLL}} u^*.
\end{aligned} 
\quad (59) $$

Since $h_L > 0$ and $h_R > 0$, we have $h^{\text{HLL}} > 0$ to obtain (57a).

Next, to exhibit $h_R^*$ and $h_L^*$ given by (57b) and (57c), we study the equation (33c). With $u_L^* = u_R^*$ and the notation $[X] = X_R - X_L$ for any quantity X, (33c) rewrites

$$ \frac{g}{8}[h^* + z]^2 = \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2. \quad (60) $$

Since the above equation is quadratic, it is now necessary to show that its right-hand side is positive; namely, we have to prove the estimate (58). Let us first consider the case $z_R - z_L = 0$. Then, since $\hat{s}_{LR}^{\text{WBAR}}$ is consistent and the CFL condition (30) holds, Lemma 3.2-(ii) gives

$$ (\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2) \Big|_{[z]=0} = \eta^{\text{HLL}} - \eta(w^{\text{HLL}}) \ge 0. $$

Next, with $z_R - z_L \neq 0$, Lemma 3.2-(i) gives the following relation:

$$ 
\begin{aligned}
\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 &= \frac{h_R h_L [u]^2}{4(h_L + h_R)} + \frac{g}{8}[h+z]^2 - \frac{[u]^2 h_R h_L (h_L u_R - h_R u_L)}{4\lambda(h_L + h_R)^2} \\
&\quad - \frac{h_L u_L + h_R u_R}{4\lambda} \left( g[z] + \frac{\hat{s}_{LR}(z_R - z_L)}{h^{\text{HLL}}} \right) + \mathcal{O}\left(\frac{1}{\lambda^2}\right).
\end{aligned} 
\quad (61) $$

Therefore, if $[u] \neq 0$ or $[h+z] \neq 0$ then there exists a $\lambda > 0$ large enough such that the inequality $\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g[z]^2}{8} \ge 0$ is satisfied. Now, assume $u_L = u_R = u$ and $[h+z] = 0$. Arguing the definition of $\hat{s}_{LR}^{\text{WBAR}}$, given by (56b), we have

$$ 
\begin{aligned}
\left( g[z] + \frac{\hat{s}_{LR}^{\text{WBAR}}[z]}{h^{\text{HLL}}} \right) \Big|_{\substack{[u]=0, \\ [h+z]=0}} &= -\sqrt{g} \frac{(h_L + h_R)}{(h^{\text{HLL}})^{\frac{5}{2}}} u[z]^2, \\
&= -2^{\frac{5}{2}} \sqrt{g} (h_L + h_R)^{-\frac{3}{2}} u[z]^2 + \mathcal{O}\left(\frac{1}{\lambda}\right).
\end{aligned} 
$$

Considering the above equality in the relation (61), we deduce

$$ 
\begin{aligned}
\left( \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 \right) \Big|_{\substack{[u]=0, \\ [h+z]=0}} &= -\frac{(h_L + h_R)u}{4\lambda} \left( g[z] + \frac{\hat{s}_{LR}^{\text{WBAR}}[z]}{h^{\text{HLL}}} \Big|_{\substack{[u]=0, \\ [h+z]=0}} \right) + \mathcal{O}\left(\frac{1}{\lambda^2}\right), \\
&= \sqrt{\frac{2g}{h_L + h_R}} \frac{u^2[z]^2}{\lambda} + \mathcal{O}\left(\frac{1}{\lambda^2}\right).
\end{aligned} 
\quad (62) $$

According to the above equation, if $u \neq 0$, there exists once again $\lambda > 0$ large enough such that the inequality $\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g[z]^2}{8} \ge 0$ holds for the cases $[h+z] = 0$ and $u_L = u_R \neq 0$. To conclude, we have to consider the last case $[h+z] = 0$ and $u_L = u_R = 0$ which defines a local equilibrium for the lake at rest (22). In this
specific situation, a direct computation using $\hat{w}^{\text{HLL}}$ and $\hat{s}_{LR}^{\text{WBAR}}$ defined by (31c) and (56b) gives

$$ 
\begin{aligned}
& h^{\text{HLL}}|_{\substack{u_L=u_R=0, \\ [h+z]=0}} = \frac{h_L+h_R}{2}, \\
& \hat{s}_{LR}^{\text{WBAR}}|_{\substack{u_L=u_R=0, \\ [h+z]=0}} = -g \frac{h_L+h_R}{2} \frac{z_R-z_L}{\Delta x}, \\
& (h^{\text{HLL}}\hat{u}^{\text{HLL}})|_{\substack{u_L=u_R=0, \\ [h+z]=0}} = -\frac{g}{4\lambda}[h^2] - \frac{g}{4\lambda}(h_L+h_R)[z] = -\frac{g}{4\lambda}(h_L+h_R)[h+z] = 0.
\end{aligned}
\quad (63) $$

Considering the three above relations in the quantity $\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g[z]^2}{8}$, we obtain

$$ (\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2)|_{\substack{u_L=u_R=0, \\ [h+z]=0}} = \frac{h_R h_L [u]^2}{4(h_L + h_R)} + \frac{g}{8}[h+z]^2 = 0. \quad (64) $$

Therefore, there always exists $\lambda > 0$ large enough such that the required inequality (58) is verified.

As a consequence, the quadratic equation (60) is well-posed and we get

$$ [h^* + z] = \pm \sqrt{8 \frac{(\tilde{\eta}^{\text{HLL}} - \tilde{\eta}(w^{\text{HLL}}))}{g} + [z]^2}. $$

Associating the above equation to (33a), we deduce the expected definition of $h_R^*$ and $h_L^*$, given by (57b) and (57c).

Concerning the statements (i) and (ii), $\hat{s}_{LR}^{\text{WBAR}}$ is consistent and a direct computation shows that $\tilde{w}(\cdot, \hat{w}, \hat{w}) = w$. In addition, as the intermediate states are defined by the system (33), Lemma 3.1 ensures that the approximate Riemann solver verifies the consistency integral relation (19) and the inequality (24). Therefore, the consistency of the Godunov-type scheme (i) and the discrete entropy stability (ii) are direct consequences of Lemma 2.2.

To establish the well-balanced property (iii), according to Lemma 2.2, we have to show that $w_L^* = w_L$ and $w_R^* = w_R$ as soon as $\hat{w}_L$ and $\hat{w}_R$ define a local lake at rest equilibrium (22), namely $u_L = u_R = 0$ and $h_L + z_L = h_R + z_R$. With $\hat{s}_{LR}$, $w_L^*$, $w_R^*$ given by (56b) and (57), we now show that $u^* = 0$, $h_L^* = h_L$ and $h_R^* = h_R$.

First, we establish that $\hat{u}^{\text{HLL}} = 0$. According to (31a) and (31b) we have

$$ (h^{\text{HLL}}\hat{u}^{\text{HLL}})|_{\substack{u_L=u_R=0, \\ [h+z]=0}} = -\frac{1}{2\lambda} \left( g \frac{h_R^2}{2} - g \frac{h_L^2}{2} \right) + \frac{(z_R - z_L)}{2\lambda} \hat{s}_{LR}|_{\substack{u_L=u_R=0, \\ [h+z]=0}}. \quad (65) $$

Since

$$ 
\begin{aligned}
\hat{s}_{LR}|_{\substack{u_L=u_R=0, \\ [h+z]=0}} &= -g h^{\text{HLL}}|_{\substack{u_L=u_R=0, \\ [h+z]=0}} \\
&= -\frac{g}{2}(h_L+h_R),
\end{aligned}
$$

we get

$$ (h^{\text{HLL}}\hat{u}^{\text{HLL}})|_{\substack{u_L=u_R=0, \\ [h+z]=0}} = -\frac{g}{4\lambda}(h_R^2 - h_L^2) - \frac{g}{4\lambda}(h_L+h_R)(z_R-z_L). $$

We notice that $z_R - z_L = h_L - h_R$ because of the assumption $[h+z] = 0$, and we immediately obtain $\hat{u}^{\text{HLL}} = 0$ as soon as $u_L = u_R = 0$ and $h_L + z_L = h_R + z_R$ that implies $u_L^* = u_R^* = 0$.
Now, we prove that $h_L^* = h_L$ and $h_R^* = h_R$. When the local lake at rest equilibrium (22) holds, Lemma 5.1 gives in (64) the equality $\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g[z]^2}{8} = 0$. Associating this equality to $h_R^*$ and $h_L^*$ given by (57b) and (57c), we obtain

$$h_R^*|_{\substack{u_L=u_R=0, \\ [h+z]=0}} = \frac{h_L+h_R-[z]}{2} = \frac{h_L+h_R+[h]}{2} = h_R,$$

$$h_L^*|_{\substack{u_L=u_R=0, \\ [h+z]=0}} = \frac{h_L+h_R+[z]}{2} = \frac{h_L+h_R-[h]}{2} = h_L,$$

that concludes the establishment of the well-balanced property (iii). The proof is thus complete. □

Let us remark that the consistent discretization $\hat{s}_{LR}^{\text{WBAR}}$, given by (56b), leads to the following equalities:

$$\hat{u}^{\text{HLL}}|_{[z]=0} = u^{\text{HLL}} \quad \text{and} \quad \left( \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 \right)|_{[z]=0} = \eta^{\text{HLL}} - \eta(w^{\text{HLL}}).$$

As a consequence, as soon as $z_L = z_R$, the intermediate states (57) degenerate toward (50). The approximate Riemann solvers of Lemma 5.1 can be understood as extensions of the solvers presented in Section 4 for flat topography.

The explicit formulations of $w_L^*$ and $w_R^*$ given by (57) completely define an approximate Riemann solver (29). Nevertheless, it is necessary to complete these formulations by a limitation technique that ensures the robustness for the dry-wet transitions.

**Theorem 5.1** (Robust, entropy-satisfying, well-balanced Godunov-type scheme for the lake at rest). Consider $\hat{w}_L$ and $\hat{w}_R$ two states of $\hat{\Omega}$ and $\hat{w}(\cdot, \hat{w}_L, \hat{w}_R) : \mathbb{R} \to \mathbb{R}^2$ an approximate Riemann solver in the form (29). Assume $\lambda > 0$ is such that the CFL-like condition (30) holds and such that the system made of (33) and (56) admits real solutions. Let us also consider the quantities $\hat{w}^{\text{HLL}}$ and $\hat{\eta}^{\text{HLL}}$ defined by (31c) and (31e), and the quantity $(\hat{h}_R^*, \hat{h}_L^*)$ in $\mathbb{R}^2$ such that

$$\tilde{h}_R^* = h^{\text{HLL}} - \frac{(z_R - z_L) \pm \sqrt{8 \frac{(\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}))}{g} + (z_R - z_L)^2}}{2}, \quad (66a)$$

$$\tilde{h}_L^* = h^{\text{HLL}} + \frac{(z_R - z_L) \pm \sqrt{8 \frac{(\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}))}{g} + (z_R - z_L)^2}}{2}. \quad (66b)$$

If $\hat{s}_{LR}$ verifies

$$\hat{s}_{LR} = \begin{cases} 0, & \text{if } h_L = 0 \text{ and } h_R = 0, \\ gh_R^2/2, & \text{if } h_R u_R = 0 \text{ and } h_L = 0 \text{ and } h_R + z_R \le z_L, \\ -gh_L^2/2, & \text{if } h_L u_L = 0 \text{ and } h_R = 0 \text{ and } h_L + z_L \le z_R, \\ -g(h_L + h_R)(z_R - z_L)/2, & \text{if } h_L = 0 \text{ or } h_R = 0, \\ \hat{s}_{LR}^{\text{WBAR}}, & \text{otherwise,} \end{cases} \quad (67)$$
with $\hat{s}_{LR}^{\text{WBAR}}$ defined by (56b) and if the intermediate states, $w_L^*$ and $w_R^*$, write

$$ \begin{pmatrix} h_L^* \\ u_L^* \\ h_R^* \\ u_R^* \end{pmatrix} = \begin{cases} (0,0,0,0)^{\text{T}}, & \text{if } h_L = 0 \text{ and } h_R = 0, \\ (0,0,2h^{\text{HLL}}, \hat{u}^{\text{HLL}})^{\text{T}}, & \text{if } h_L = 0 \text{ and } h_R > 0, \\ (2h^{\text{HLL}}, \hat{u}^{\text{HLL}}, 0, 0)^{\text{T}}, & \text{if } h_L > 0 \text{ and } h_R = 0, \end{cases} \begin{pmatrix} \min(\max(\hat{h}_L^*, 0), 2h^{\text{HLL}}) \\ \hat{u}^{\text{HLL}} \\ \min(\max(\hat{h}_R^*, 0), 2h^{\text{HLL}}) \\ \hat{u}^{\text{HLL}} \end{pmatrix}, \quad \text{otherwise,} \tag{68} $$

then the Godunov-type scheme (16) defined by such an approximate Riemann solver

(i) is consistent with the shallow water equations (1),

(ii) preserves the convex set $\hat{\Omega}$, i.e., if $(\hat{w}_i^n)_{i \in \mathbb{Z}} \subset \hat{\Omega}$ then, $(\hat{w}_i^{n+1})_{i \in \mathbb{Z}} \subset \hat{\Omega}$,

(iii) is robust for the dry-wet transitions,

(iv) is well-balanced for the lake at rest (10).

In addition, if $\tilde{h}_L^* > 0$ and $\tilde{h}_R^* > 0$ then the Godunov-type scheme (16) associated to such an approximate Riemann solver satisfies a discrete entropy inequality (14) for the couple $(\hat{\eta}, \hat{G})$ defined by (7) with a numerical entropy flux function $\tilde{G}_{i+\frac{1}{2}}$ given by (25).

The well-balanced property and the discrete entropy inequality detailed in the above theorem are also given in [12] with a relaxation scheme that needs to solve a cubic equation. The result of Theorem 5.1 overcomes this constraint and its main originality arises from the explicit solution of a quadratic equation. This explicit solving gives two numerical schemes distinguished by the symbol $\pm$ in the definitions (66). Both schemes are entropy-satisfying, well-balanced for the lake at rest and preserves the set $\hat{\Omega}$ thanks to the limitation techniques $\min(\max(\cdot, \cdot), \cdot)$. However, as soon as the limitation is active, the system (33) is no longer necessarily satisfied and consequently, the entropy stability might be locally lost in the dry-wet transitions.

These transitions are also treated with the several cases in the equalities (67) and (68). According to Section 3.1.2.4 of [39], these equalities are robust, preserve the well-balanced property of the scheme defined Lemma 5.1 but it is clear that they are not continuous with it.

We now prove Theorem 5.1.

*Proof.* Concerning (i), the consistency has only to be proved in the wet regions. Therefore, if $\tilde{h}_L^*$ and $\tilde{h}_R^*$ given by (66) are both positive then (68) coincides with (57). In this case, the consistency is shown in Lemma 5.1. Next, if $\tilde{h}_L^* \le 0$ or $\tilde{h}_R^* \le 0$, then the limitation technique $\min(\max(\cdot, \cdot), \cdot)$ is active but, following the arguments as in the proof of Lemma 4.2, the integral consistency relation is shown to be preserved. As a consequence, we deduce the consistency of the numerical scheme.

Now, the preservation of the convex set $\hat{\Omega}$, given by (ii), comes from $h_L^* \ge 0$ and $h_R^* \ge 0$. According to (68), these inequalities are obviously satisfied.

Concerning (iii) and according to Section 3.1.2.4 of [39], the definitions (67) for $\hat{s}_{LR}$ and (68) for the states $w_L^*$ and $w_R^*$ ensure the robustness of the scheme in wet-dry transitions.

The last statement (iv) easily comes from Lemma 5.1-(iii), the source term definition (67) and the intermediate state definition (68).

Finally, concerning the discrete entropy inequality, if $\tilde{h}_L^*$ and $\tilde{h}_R^*$, given by (66), are both positive then (68) coincides with (57). In this case, the discrete entropy inequality is a direct consequence of the approximate Riemann solver definition and Lemma 5.1-(ii), that concludes the proof. $\square$
The schemes designed in Theorem 5.1 are obtained by adopting (33) to enforce both consistency and entropy preservation completed by the source term discretization (56b). An other closure is now presented in order to define an entropy-satisfying well-balanced scheme for the moving equilibrium.

## 6. A FULLY WELL-BALANCED ENTROPY-STABLE NUMERICAL SCHEME FOR THE GENERAL EQUILIBRIUM

In this section, the equations (33) are completed with a suitable source term discretization in order to achieve a well-balanced property for the moving equilibrium (9).

Considering an interface separating $(w_i^n, z_i) = (w_L, z_L) = \hat{w}_L$ and $(w_{i+1}^n, z_{i+1}) = (w_R, z_R) = \hat{w}_R$, the local moving equilibrium (23) now reads

$$h_L u_L = h_R u_R, \quad \text{and} \quad \hat{B}(\hat{w}_L) = \hat{B}(\hat{w}_R), \qquad (69)$$

where $\hat{B}(\hat{w})$ is defined by (21).

For the sake of clarity in the notations, for any quantity X we set

$$X^{\text{eq}} = X|_{[hu]=0}[\hat{B}=0].$$

Since the local moving equilibrium (69) imposes $h_L u_L = h_R u_R$, the notation $q^{\text{eq}} \in \mathbb{R}$ is now defined as follows:

$$q^{\text{eq}} = (h_L u_L)^{\text{eq}} = (h_R u_R)^{\text{eq}}.$$

Finally, $\bar{h}$ denotes the arithmetic mean of $h_L$ and $h_R$ which is written as $\bar{h} = \frac{(h_L+h_R)}{2}$ in $\mathbb{R}^+$.

In addition, we need to introduce some technical quantities, $\zeta_{LR}$, $q_{LR}$, $\beta_{LR}$ and $F_{LR}^2$, which will take particular values when $\hat{w}_L$ and $\hat{w}_R$ define a local moving equilibrium (69). These quantities are given by

$$\zeta_{LR} = [hu]^2 + \bar{h}[\hat{B}(\hat{w})]^2/g, \qquad (70a)$$

$$q_{LR} = \begin{cases} 0, & \text{if } [z] = 0 \text{ and } \zeta_{LR} = 0, \\ \frac{\min(h_L u_L, h_R u_R)(g[z]\Delta t)^2}{(g[z]\Delta t)^2 + \zeta_{LR}}, & \text{otherwise,} \end{cases}, \qquad (70b)$$

$$\beta_{LR} = \begin{cases} 0, & \text{if and } [u] = 0 \text{ and } [h+z] = 0 \text{ and } \zeta_{LR} = 0, \\ \frac{(\bar{h}[u])^2 + g\bar{h}[h+z]^2}{(\bar{h}[u])^2 + g\bar{h}[h+z]^2 + \zeta_{LR}}, & \text{otherwise,} \end{cases}, \qquad (70c)$$

$$F_{LR}^2 = \frac{q_{LR}^2}{gh_L h_R h^{\text{HLL}}}. \qquad (70d)$$

Straightforward computations show the above quantities stay bounded and the following technical lemma gives properties satisfied by $\zeta_{LR}$, $q_{LR}$ and $F_{LR}^2$.

**Lemma 6.1.** Let us consider $\hat{w}_L$ and $\hat{w}_R$ two states with $h_L > 0$ and $h_R > 0$. Let us consider $\lambda > 0$, $\Delta t > 0$ satisfying the CFL-like condition (30). Let us also consider $w^{\text{HLL}}$ given by (31a) and $\zeta_{LR}$, $q_{LR}$, $F_{LR}^2$ defined by (70). The following statements are satisfied:

(i) $\zeta_{LR} = 0$ if and only if $\hat{w}_L$ and $\hat{w}_R$ verify a local equilibrium (69).

(ii) If $[z] \neq 0$ then $q_{LR}^{\text{eq}} = q^{\text{eq}}$.

(iii) If $[z] \neq 0$ then $(F_{LR}^2)^{\text{eq}} = \frac{(q^{\text{eq}})^2}{(gh_L h_R h)}$.
Proof. Concerning (i), since $h_L > 0$ and $h_R > 0$, we have $\bar{h} > 0$. As a consequence, $\zeta_{LR}$ is a positive quantity which vanishes if and only if $[hu] = 0$ and $[\dot{B}(\hat{w})] = 0$.

For the statement (ii) with $[z] \neq 0$, the equality $\zeta_{LR}^{\text{eq}} = 0$ gives

$$q_{LR}^{\text{eq}} = \left( \frac{\min(h_L u_L, h_R u_R)(g[z]\Delta t)^2}{(g[z]\Delta t)^2 + \zeta_{LR}} \right)^{\text{eq}} = (\min(h_L u_L, h_R u_R))^{\text{eq}} = q^{\text{eq}}.$$

Finally, for (iii), writing the definition of $h^{\text{HLL}}$ given by (31a), we obtain

$$(h^{\text{HLL}})^{\text{eq}} = \frac{h_L + h_R}{2} - \frac{[hu]^{\text{eq}}}{2\lambda} = \frac{h_L + h_R}{2} = \bar{h}.$$

Using the above relation and the statement (ii) in the definition of $F_{LR}^2$ given by (70d), we deduce the expected result that concludes the proof. □

Considering the definitions (70), the equations (33) are completed as follows:

$$h_R^* h_L^* (u_R^* - u_L^*)^2 = \frac{q_{LR}^2}{h_R h_L} (h_R^* - h_L^*)^2, \qquad (71a)$$

$$\hat{s}_{LR} = \hat{s}_{LR}^{\text{FWB}} = \hat{s}_{LR}^{\text{WBAR}} + \beta_{LR} \left( \frac{g\bar{h}F_{LR}^2}{4h_Lh_R} \frac{(h_R - h_L)^3}{\Delta x} + \sqrt{g} \frac{h_L u_L + h_R u_R}{(h^{\text{HLL}})^{\frac{3}{2}}} \frac{(z_R - z_L)^2}{\Delta x} \right). \qquad (71b)$$

The notation $\hat{s}_{LR}^{\text{FWB}}$ means *Fully Well-Balanced*. According to (71b), $\hat{s}_{LR}^{\text{FWB}}$ is nothing but a correction of $\hat{s}_{LR}^{\text{WBAR}}$ given by (56b).

The solving of the system made of (33) and (71) is now detailed in the following lemma. From now on, it is worth noticing that the obtained solutions will depend on the usual $\pm$ signs. In fact, $\pm$ signs will be involved twice in the solutions. So, to avoid any possible confusion, we introduce two independent symbols $\pm_1$ and $\pm_2$ such that if $\pm_{1,2}$ is positive (resp. negative) then $\mp_{1,2}$ is negative (resp. positive).

**Lemma 6.2.** Let us consider $\hat{w}_L$ and $\hat{w}_R$ two states with $h_L > 0$ and $h_R > 0$. Let us also consider $\tilde{w}(\cdot, \hat{w}_L, \hat{w}_R) : \mathbb{R} \to \mathbb{R}^2$ an approximate Riemann solver (29). Assume the CFL-like condition (30) holds. Then, there exists $\lambda > 0$ large enough and $\Delta t > 0$ small enough such that the system made of (33) and (71) admits four solutions given by

$$h_L^* = h^{\text{HLL}} - \frac{[h^*]}{2}, \qquad (72a)$$

$$h_R^* = h^{\text{HLL}} + \frac{[h^*]}{2}, \qquad (72b)$$

$$u_L^* = \hat{u}^{\text{HLL}} \mp_2 \frac{q_{LR}[h^*]}{2h^{\text{HLL}}\sqrt{h_L h_R}} \sqrt{\frac{h_R^*}{h_L^*}}, \qquad (72c)$$

$$u_R^* = \hat{u}^{\text{HLL}} \pm_2 \frac{q_{LR}[h^*]}{2h^{\text{HLL}}\sqrt{h_L h_R}} \sqrt{\frac{h_L^*}{h_R^*}}, \qquad (72d)$$

with

$$[h^*] = -\frac{(z_R - z_L)}{1 + F_{LR}^2} \pm_1 \sqrt{\frac{8 \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}})}{g(1 + F_{LR}^2)} + \frac{(z_R - z_L)^2}{(1 + F_{LR}^2)^2}}, \qquad (73)$$

where $\hat{w}^{\text{HLL}}$ and $\hat{\eta}^{\text{HLL}}$ are defined by (31c) and (31e).

Finally, a Godunov-type scheme (16) defined by such approximate Riemann solvers

(i) are consistent with the shallow water equations (1),
(ii) satisfy a discrete entropy inequality (14) for the couple $(\hat{\eta}, \hat{G})$, defined by (7), where the entropy numerical flux function $\hat{G}_{i+\frac{1}{2}}$ is given by (25).

The symbols $\pm_1$ and $\pm_2$, adopted in (72) and (73), may be interpreted as free parameters that will be set in the sequel to guarantee the well-balanced property.

Before proving Lemma 6.2, it is now necessary to give two intermediate results. The first one concerns the behavior of $\hat{s}_{LR}^{\text{FWB}}$ defined by (71b).

**Lemma 6.3 (Property of $\hat{s}_{LR}^{\text{FWB}}$).** Consider $\hat{w}_L$ and $\hat{w}_R$ two states with $h_L > 0$ and $h_R > 0$. Assume that the CFL-like condition (30) holds. In the expression of $h^{\text{HLL}}\hat{u}^{\text{HLL}}$, given by (31b), impose $\hat{s}_{LR} = \hat{s}_{LR}^{\text{FWB}}$, then

$$ (h^{\text{HLL}}\hat{u}^{\text{HLL}})^{\text{eq}} = \left( (hu)^{\text{HLL}} + \frac{\hat{s}_{LR}^{\text{FWB}}(z_R - z_L)}{2\lambda} \right)^{\text{eq}} = q^{\text{eq}}. $$

*Proof.* The proof is led through three exhaustive cases given by

(i) $[z] = 0$,

(ii) $[z] \neq 0$ and $q^{\text{eq}} = 0$,

(iii) $[z] \neq 0$ and $q^{\text{eq}} \neq 0$.

In the case (i), and according to the definition (71b), we immediately have $\hat{s}_{LR}^{\text{FWB}} = 0$. In addition, a local smooth equilibrium defined by $[z] = 0$ necessarily implies $\hat{w}_L = \hat{w}_R$. As a consequence, using the definition of $h^{\text{HLL}}\hat{u}^{\text{HLL}}$ given by (31b), we get

$$ \begin{aligned} (h^{\text{HLL}}\hat{u}^{\text{HLL}})^{\text{eq}}|_{[z]=0} &= ((hu)^{\text{HLL}})^{\text{eq}}, \\ &= \frac{(h_L u_L)^{\text{eq}} + (h_R u_R)^{\text{eq}}}{2} - \frac{[hu^2 + \frac{gh^2}{2}]^{\text{eq}}}{2\lambda}, \\ &= q^{\text{eq}}, \end{aligned} $$

and the result is established in the case (i).

The case (ii) is defined by $[z] \neq 0$ and $q^{\text{eq}} = 0$. A direct computation shows that the local moving equilibrium (69) degenerates towards the local lake at rest equilibrium that is expressed as $u_L = u_R = 0$ and $[h + z] = 0$. But, according to the definition of $\beta_{LR}$ given by (70c), if the local lake at rest equilibrium occurs, then $\beta_{LR} = 0$. Therefore, using the identity (71b), we have

$$ (\hat{s}_{LR}^{\text{FWB}})^{\text{eq}}|_{q^{\text{eq}}=0} = (\hat{s}_{LR}^{\text{WBAR}})^{\text{eq}}|_{q^{\text{eq}}=0} = -g \frac{h_L + h_R}{2} \frac{z_R - z_L}{\Delta x}. $$

As a consequence, arguing (65), we deduce $h^{\text{HLL}}\hat{u}^{\text{HLL}} = 0$. Since $h^{\text{HLL}} > 0$, this concludes the second case.

Let us consider now the case (iii) defined by $[z] \neq 0$ and $q^{\text{eq}} \neq 0$. Since $\hat{w}_L$ and $\hat{w}_R$ define a local moving equilibrium (69) with $q^{\text{eq}} \neq 0$, from (70c) and (69), we obtain

$$ \beta_{LR}^{\text{eq}}|_{q^{\text{eq}} \neq 0} = 1 \quad \text{and} \quad -\frac{(q^{\text{eq}})^2}{2h_L^2 h_R^2} [h^2] + g[h + z] = 0. $$

Arguing the above identities and Lemma 6.1-(iii) to get $(F_{LR}^2)^{\text{eq}} = 2 \frac{(q^{\text{eq}})^2}{(gh_L h_R (h_L + h_R))}$, then we have

$$ \begin{aligned} (\hat{s}_{LR}^{\text{FWB}}(z_R - z_L))^{\text{eq}}|_{q^{\text{eq}} \neq 0} &= -g \frac{h_L + h_R}{2}[z] + \frac{g}{8}[h]^3 \frac{h_L + h_R}{h_L h_R} F_{LR}^2|_{q^{\text{eq}} \neq 0}, \\ &= -g \frac{h_L + h_R}{2}[z] + \frac{[h]^2}{4} \frac{[h^2](q^{\text{eq}})^2}{h_L^2 h_R^2 (h_L + h_R)}, \\ &= -g \frac{h_L + h_R}{2}[z] + \frac{[h]^2}{2} \frac{[h + z]}{h_L + h_R}. \end{aligned} $$
Plugging this last identity into the definition of $h^{\text{HLL}}\hat{u}^{\text{HLL}}$, given by (31b), we deduce

$$ (h^{\text{HLL}}\hat{u}^{\text{HLL}})^{\text{eq}}|_{q^{\text{eq}} \neq 0} = \left(h^{\text{HLL}}u^{\text{HLL}} + \frac{\hat{s}_{LR}^{\text{FWB}}(z_R - z_L)}{2\lambda}\right)^{\text{eq}}|_{q^{\text{eq}} \neq 0}, \\ = q^{\text{eq}} - \frac{1}{2(h_L + h_R)\lambda}\left(-\frac{(q^{\text{eq}})^2[h^2]}{h_R h_L} + \frac{g}{2}(h_L + h_R)^2[h]\right) + \frac{g}{4\lambda}\frac{-(h_L + h_R)^2[z] + [h]^2[h+z]}{h_L + h_R}, \\ = q^{\text{eq}} - \frac{g[h+z]}{4\lambda(h_L + h_R)}(-4h_Lh_R + (h_L + h_R)^2 - [h]^2), \\ = q^{\text{eq}}. $$

The proof is thus complete.

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

The second technical intermediate result to establish Lemma 6.2 is an inequality that defines a necessary and sufficient condition for the existence of $h_L^*$ and $h_R^*$ satisfying the system made of (33) and (71).

**Lemma 6.4.** Consider $\hat{w}_L$ and $\hat{w}_R$ two states with $h_L > 0$ and $h_R > 0$. Assume that the CFL-like condition (30) holds and $\hat{s}_{LR} = \hat{s}_{LR}^{\text{FWB}}$. Let us adopt the notation $h^{\text{HLL}}\hat{u}^{\text{HLL}}$ given by (31b). Consider the quantities $\hat{w}^{\text{HLL}}$, $\hat{\eta}^{\text{HLL}}$ and $F_{LR}^2$ defined by (31c), (31e) and (70d). Then there exists $\lambda > 0$ large enough and $\Delta t > 0$ small enough such that the following inequality is satisfied:

$$ \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 - \frac{g[z]^2}{8} \frac{F_{LR}^2}{1 + F_{LR}^2} \geq 0. \quad (74) $$

*Proof.* The proof is led through three exhaustive cases given by

(i) $[z] = 0$,

(ii) $[z] \neq 0$ and $\zeta_{LR} = 0$,

(iii) $[z] \neq 0$ and $\zeta_{LR} \neq 0$.

First, we consider the case (i) defined by $[z] = 0$. Using the definition (70b), we have $q_{LR} = 0$ and using the equation (71b), we obtain $\hat{s}_{LR}^{\text{FWB}} = 0$. As a consequence, in the case $[z] = 0$, we deduce $F_{LR}^2 = 0$ and $\hat{u}^{\text{HLL}} = u^{\text{HLL}}$. Next, an analogous computation of (44) gives

$$ \left( \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 - \frac{g[z]^2}{8} \frac{F_{LR}^2}{1 + F_{LR}^2} \right) \Big|_{[z]=0} \\ = \eta^{\text{HLL}} - \frac{(h^{\text{HLL}}\hat{u}^{\text{HLL}})^2}{2h^{\text{HLL}}} \Big|_{[z]=0} - \frac{g}{2}(h^{\text{HLL}})^2 + \left( \frac{g}{4}[z] \left( [h] - \frac{h_L u_L + h_R u_R}{\lambda} \right) \right) \Big|_{[z]=0} \\ = \eta^{\text{HLL}} - \eta(w^{\text{HLL}}). $$

Since the CFL condition (30) is satisfied, the inequality $\eta^{\text{HLL}} - \eta(w^{\text{HLL}}) \geq 0$ is ensured (see Lem. 3.2-(ii) for details).

The case (ii) is defined by $[z] \neq 0$ and $\zeta_{LR} = 0$. Since $\zeta_{LR} = 0$, $\hat{w}_L$ and $\hat{w}_R$ define a local moving equilibrium (69) and a direct computation gives $\hat{G}(\hat{w}_R) - \hat{G}(\hat{w}_L) = [\hat{B}]q^{\text{eq}} = 0$. According to Lemma 6.3, the equality $(h^{\text{HLL}}\hat{u}^{\text{HLL}})^{\text{eq}} = q^{\text{eq}}$ is also verified. As a consequence, we have

$$ (\hat{\eta}^{\text{HLL}})^{\text{eq}} - \hat{\eta}(\hat{w}^{\text{HLL}})^{\text{eq}} = \frac{\eta(w_R)^{\text{eq}} + \eta(w_L)^{\text{eq}}}{2} + \frac{g}{2}(h_L z_L + h_R z_R) - \frac{\hat{G}(\hat{w}_R)^{\text{eq}} - \hat{G}(\hat{w}_L)^{\text{eq}}}{2\lambda} \\ - \left( \frac{(h^{\text{HLL}}\hat{u}^{\text{HLL}})^2}{2h^{\text{HLL}}} \right)^{\text{eq}} - \left( \frac{g}{2}(h^{\text{HLL}})^2 \right)^{\text{eq}} - g(h^{\text{HLL}})^{\text{eq}} \frac{z_L + z_R}{2}, $$
$$ 
\begin{align*}
&= \frac{(q^{\text{eq}})^2}{4} \left( \frac{1}{h_L} + \frac{1}{h_R} \right) + \frac{g}{4}(h_L^2 + h_R^2) + \frac{g}{2}(h_L z_L + h_R z_R) \\
&\quad - \frac{(q^{\text{eq}})^2}{h_L + h_R} - \frac{g(h_L + h_R)^2}{8} - \frac{g}{4}(h_L + h_R)(z_L + z_R), \\
&= \frac{(q^{\text{eq}})^2 h_L h_R}{4(h_L + h_R)} \left[ \frac{1}{h} \right]^2 + \frac{g}{8}[h]^2 + \frac{g}{4}[z][h], \\
&= \frac{g}{8} \left( 1 + (F_{LR}^2)^{\text{eq}} \right) [h]^2 + \frac{g}{4}[z][h].
\end{align*}
$$

Using the above equation, we obtain the required estimate (74) as follows:

$$ 
\begin{align*}
\left( \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 - \frac{g[z]^2}{8} \frac{F_{LR}^2}{1 + F_{LR}^2} \right)^{\text{eq}} &= \frac{g}{8} \left( 1 + (F_{LR}^2)^{\text{eq}} \right) [h]^2 + \frac{g}{4}[h][z] + \frac{g}{8} \frac{[z]^2}{1 + (F_{LR}^2)^{\text{eq}}} \\
&= \frac{g}{8} \left( \sqrt{1 + (F_{LR}^2)^{\text{eq}}}[h] + \frac{[z]}{\sqrt{1 + (F_{LR}^2)^{\text{eq}}}} \right)^2 \geq 0.
\end{align*}
$$

Finally, the last case (iii) is such that $[z] \neq 0$ and for couples $\hat{w}_L$ and $\hat{w}_R$ that not define a local moving equilibrium (69). This last case is expressed as $[z] \neq 0$ and $\zeta_{LR} \neq 0$ and Lemma 3.2-(i) gives the following estimate:

$$ 
\begin{align}
\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 &= \frac{h_R h_L [u]^2}{4(h_L + h_R)} + \frac{g}{8}[h+z]^2 - \frac{[u]^2 h_R h_L (h_L u_R - h_R u_L)}{4\lambda(h_L + h_R)^2} \nonumber \\
&\quad - \frac{h_L u_L + h_R u_R}{4\lambda} \left( g[z] + \frac{\hat{s}_{LR}(z_R - z_L)}{h^{\text{HLL}}} \right) + \mathcal{O}\left(\frac{1}{\lambda^2}\right). \tag{75}
\end{align}
$$

Then, as soon as $[u] \neq 0$ or $[h+z] \neq 0$, there exists a $\lambda > 0$ large enough such that the following inequality is satisfied

$$ \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g[z]^2}{8} > 0. \tag{76} $$

For the particular cases $[u] = 0$ and $[h+z] = 0$, the following equations are verified:

$$ \beta_{LR}|_{[u]=0, [h+z]=0} = 0 \quad \text{and} \quad \hat{s}_{LR}^{\text{FWB}}|_{\beta_{LR}=0} = \hat{s}_{LR}^{\text{WBAR}}. $$

As a consequence, the existence of the inequality (76) is linked to the proof of Lemma 5.1, which concerns the estimate (75) with $\hat{s}_{LR} = \hat{s}_{LR}^{\text{WBAR}}$. But, according to the computations (62) and (63), the inequality (76) can be held except for $\hat{w}_L$ and $\hat{w}_R$ defined by a local lake at rest equilibrium (22). Since, the lake at rest coincides with a moving equilibrium with a null velocity, it has already been treated in the case (ii). As a consequence, the inequality (76) can always be ensured for this case (iii). In addition, according to the definitions (70b) and (70d), and since $\zeta_{LR} \neq 0$, then $F_{LR}^2$ is proportional to $\Delta t^4$. Then, if the inequality (76) is verified, the large inequality (74) can be satisfied with $\Delta t > 0$ small enough such that

$$ \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 \geq \frac{g[z]^2}{8} \frac{F_{LR}^2}{1 + F_{LR}^2}. \tag{77} $$

We thus easily deduce the expected inequality (74) for this last case that concludes the proof. $\square$

The inequalities (76), (77) stand for an additional CFL condition to the standard one given by (30). As a consequence, $\lambda$ and $\Delta t$ have to be selected according to the more restrictive of the two. Equipped with both above lemmas, Lemma 6.2 is now established.
Proof of Lemma 6.2. First, we show that the formulations (72) are solutions of the system made of (33) and (71). Plugging (71a) into (33c) and according to the definition of $F_{LR}^2$ given by (70d), we deduce the two following equations satisfied by $h_L^*$ and $h_R^*$:

$$ \frac{h_L^* + h_R^*}{2} = h^{\text{HLL}}, \qquad (78a) $$

$$ \frac{g}{8}(1 + F_{LR}^2)[h^*]^2 + \frac{g}{4}[z][h^*] = \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}). \qquad (78b) $$

Now, we have to show that the above quadratic equation is well-posed. Dividing this quadratic equation (78b) on both sides by $\frac{g}{8(1+F_{LR}^2)}$ then adding $\frac{[z]^2}{(1+F_{LR}^2)^2}$, we obtain

$$ \begin{aligned} \left( [h^*] + \frac{[z]}{1 + F_{LR}^2} \right)^2 &= \frac{8}{g} \frac{\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}})}{1 + F_{LR}^2} + \frac{[z]^2}{(1 + F_{LR}^2)^2}, \\ &= \frac{8}{g(1 + F_{LR}^2)} \left( \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 - \frac{g[z]^2}{8} \frac{F_{LR}^2}{1 + F_{LR}^2} \right). \end{aligned} $$

According to Lemma 6.4, if $\lambda$ (resp. $\Delta t$) is large (resp. small) enough then the right-hand side of the above equation is positive. As a consequence, as soon as $\lambda$ and $\Delta t$ are well-chosen, the previous equation is well-posed and a direct computation leads to $[h^*]$ given by (73). Next, coupling $[h^*] = h_R^* - h_L^*$ to (78a), we deduce the values of $h_L^*$ and $h_R^*$ as given by (72a) and (72b).

Since the quantities $h_R^*$ and $h_L^*$ are known and since they are assumed to be positive, the equation (71a) re-writes

$$ [u^*] = \pm \frac{q_{LR}}{\sqrt{h_L h_R}} \frac{[h^*]}{\sqrt{h_L^* h_R^*}}. \qquad (79) $$

Using the above equation and (78a) in the relation $\frac{(h_L^* u_L^* + h_R^* u_R^*)}{2} = h^{\text{HLL}} \hat{u}^{\text{HLL}}$, we obtain

$$ \begin{aligned} h^{\text{HLL}} \hat{u}^{\text{HLL}} &= \frac{h_L^* u_L^* + h_R^* u_R^*}{2}, \\ &= \frac{h_L^* + h_R^*}{2} \frac{u_L^* + u_R^*}{2} + \frac{[h^*][u^*]}{4}, \\ &= h^{\text{HLL}} \frac{u_L^* + u_R^*}{2} \pm \frac{q_{LR}}{4\sqrt{h_L h_R}} \frac{[h^*]^2}{\sqrt{h_L^* h_R^*}}. \end{aligned} \qquad (80) $$

This equation associated to the equalities (79) and (78a) gives

$$ \begin{aligned} u_R^* &= \hat{u}^{\text{HLL}} \pm \frac{q_{LR}[h^*]}{2\sqrt{h_L h_R h_L^* h_R^*}} \mp \frac{q_{LR}[h^*]^2}{4h^{\text{HLL}}\sqrt{h_L h_R}\sqrt{h_R^* h_L^*}}, \\ &= \hat{u}^{\text{HLL}} \pm \frac{q_{LR}[h^*]}{2\sqrt{h_L h_R h_L^* h_R^*}} \left(1 - \frac{[h^*]}{2h^{\text{HLL}}}\right), \\ &= \hat{u}^{\text{HLL}} \pm \frac{q_{LR}[h^*]}{2\sqrt{h_L h_R h_L^* h_R^*}} \frac{h_L^*}{h^{\text{HLL}}}. \end{aligned} $$

From this last equality, we get the formulation for $u_R^*$ presented in (72d). Since $u_L^*$, given by (72c), can be derived from an analogous computation, this achieves to prove (72).

Concerning the statements (i) and (ii), $\hat{s}_{LR}^{\text{FWB}}$ is consistent and a direct computation shows $\tilde{w}(\cdot, \hat{w}, \hat{w}) = w$. In addition, as the intermediate states are defined by the system (33), Lemma 3.1 ensures that the approximate Riemann solver verifies the integral consistency relation (19) and the entropy inequality (24). Therefore, the consistency of the Godunov-type scheme (i) and the discrete entropy stability (ii) are direct consequences of Lemma 2.2 which concludes the proof. $\square$
As previously underlined, (72) exhibits four solutions which depend on the choices of the symbols $\pm_1$ and $\pm_2$. The selection of one solution must be made according to the well-balanced property.

In order to satisfy this property, $\pm_2$ must be negative (which imposes $\mp_2$ to be positive) and a direct computation shows that the choice of $\pm_1$ must be made according to a condition obtained when $\hat{w}_L$ and $\hat{w}_R$ define a local equilibrium (69). This condition formally writes

$$ \pm_1^{\text{eq}} = \text{sign}((1 + (F_{LR}^2)^{\text{eq}})[h] + [z]). $$

As a consequence, a simple formulation for the symbol $\pm_1$ is given by

$$ \pm_1 = \text{sign}((1 + (F_{LR}^2))[h] + [z]). $$

This definition of $\pm_1$ will be adopted but other choices are possible. This formulation leads to expressions for $(h_L^*, u_L^*)$ and $(h_R^*, u_R^*)$ now detailed in the following result.

**Theorem 6.1** (Robust, entropy-satisfying, well-balanced Godunov-type scheme for all smooth equilibrium). Consider $\hat{w}_L$ and $\hat{w}_R$ two states of $\Omega$ and $\tilde{w}(\cdot, \hat{w}_L, \hat{w}_R) : \mathbb{R} \to \mathbb{R}^2$ an approximate Riemann solver in the form (29). Assume $\lambda > 0$ and $\Delta t > 0$ such that the CFL-like condition (30) holds and such that the system (33), (71) admits real solutions. Let us also consider the quantities $\hat{w}^{\text{HLL}}$, $\hat{\eta}^{\text{HLL}}$, $q_{LR}$ and $F_{LR}^2$ defined by (31c), (31e) (70b), (70d) and the quantity $(\tilde{h}_L^*, \tilde{u}_L^*, \tilde{h}_R^*, \tilde{u}_R^*)$ in $\mathbb{R}^4$ given by

$$ \tilde{h}_L^* = h^{\text{HLL}} - \frac{[h^*]}{2}, \qquad (81a) $$

$$ \tilde{u}_L^* = \hat{u}^{\text{HLL}} + \frac{q_{LR}[h^*]}{2h^{\text{HLL}}\sqrt{h_L h_R}} \sqrt{\frac{h_R^*}{h_L^*}}, \qquad (81b) $$

$$ \tilde{h}_R^* = h^{\text{HLL}} + \frac{[h^*]}{2}, \qquad (81c) $$

$$ \tilde{u}_R^* = \hat{u}^{\text{HLL}} - \frac{q_{LR}[h^*]}{2h^{\text{HLL}}\sqrt{h_L h_R}} \sqrt{\frac{h_L^*}{h_R^*}}, \qquad (81d) $$

with

$$ [h^*] = -\frac{(z_R - z_L)}{1 + F_{LR}^2} + \text{sign}((1 + F_{LR}^2)[h] + [z]) \sqrt{\frac{8 \hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}})}{g(1 + F_{LR}^2)} + \frac{(z_R - z_L)^2}{(1 + F_{LR}^2)^2}}. \qquad (82) $$

In addition, consider $\hat{s}_{LR}$ in the form

$$ \hat{s}_{LR} = \begin{cases} 0, & \text{if } h_L = 0 \text{ and } h_R = 0, \\ gh_R^2/2, & \text{if } h_R u_R = 0 \text{ and } h_L = 0 \text{ and } h_R + z_R \le z_L, \\ -gh_L^2/2, & \text{if } h_L u_L = 0 \text{ and } h_R = 0 \text{ and } h_L + z_L \le z_R, \\ -g(h_L + h_R)(z_R - z_L)/2, & \text{if } h_L = 0 \text{ or } h_R = 0, \\ \hat{s}_{LR}^{\text{FWB}}, & \text{otherwise,} \end{cases} \qquad (83) $$
with $\hat{s}_{LR}^{\text{FWB}}$ defined by (71b). For $\varepsilon > 0$, assume that the intermediate states, $w_L^*$ and $w_R^*$, are given by

$$ \begin{pmatrix} h_L^* \\ u_L^* \\ h_R^* \\ u_R^* \end{pmatrix} = \begin{cases} (0,0,0,0)^{\text{T}}, & \text{if } h_L = 0 \text{ and } h_R = 0, \\ (0,0,2h^{\text{HLL}}, \hat{u}^{\text{HLL}})^{\text{T}}, & \text{if } h_L = 0 \text{ and } h_R > 0, \\ (2h^{\text{HLL}}, \hat{u}^{\text{HLL}}, 0,0)^{\text{T}}, & \text{if } h_L > 0 \text{ and } h_R = 0, \end{cases} \\ \begin{pmatrix} \min(\max(\widetilde{h}_L^*, \varepsilon), 2h^{\text{HLL}} - \varepsilon) \\ \widetilde{u}_L^* \\ \min(\max(\widetilde{h}_R^*, \varepsilon), 2h^{\text{HLL}} - \varepsilon) \\ \widetilde{u}_R^* \end{pmatrix}, & \text{otherwise.} \end{cases} \tag{84} $$

The Godunov-type scheme (16) defined by such an approximate Riemann solver

(i) is consistent with the shallow water equations (1),

(ii) preserves the convex set $\hat{\Omega}$, i.e., if $(\hat{w}_i^n)_{i \in \mathbb{Z}} \subset \hat{\Omega}$ then, $(\hat{w}_i^{n+1})_{i \in \mathbb{Z}} \subset \hat{\Omega}$, and it is robust for the dry-wet transitions,

(iv) is well-balanced for all moving smooth equilibrium (9).

In addition, if $\widetilde{h}_L^* > 0$ and $\widetilde{h}_R^* > 0$ then the Godunov-type scheme (16) associated to such an approximate Riemann solver satisfies a discrete entropy inequality (14) for the couple $(\hat{\eta}, \hat{G})$, defined by (7), with a numerical entropy flux function $\hat{G}_{i+\frac{1}{2}}$ given by (25).

*Proof.* Concerning (i), the consistency only needs to be proved in the wet regions. Therefore, if $\widetilde{h}_L^*$ and $\widetilde{h}_R^*$ given by (81a) and (81c) are both positive then (84) coincides with (81). In this case, the consistency is shown in Lemma 6.2-(i). If $\widetilde{h}_L^* \le 0$ or $\widetilde{h}_R^* \le 0$, then the limitation technique $\min(\max(\cdot, \cdot), \cdot)$ works and an analogous computation to the one done in Lemma 4.2-(i) shows that these procedures preserve the consistency integral relation for all $\varepsilon > 0$ small enough. As a consequence, we deduce the consistency statement (i).

For the statement (ii) related to the preservation of the convex set $\hat{\Omega}$, it is sufficient to prove the inequalities $h_L^* \ge 0$ and $h_R^* \ge 0$ but according to (84), these inequalities are obviously ensured for $\varepsilon > 0$ being small enough. Moreover, the robustness of the scheme in wet-dry transitions is ensured by Section 3.1.2.4 of [39] with $\hat{s}_{LR}$ given by (83) and $w_L^*$ and $w_R^*$ given by (84).

For the well-balanced property (iii), according to Lemma 2.2, it is sufficient to show that if $\hat{w}_L$ and $\hat{w}_R$ define a local moving equilibrium (69) then $w_L^* = w_L$ and $w_R^* = w_R$. Adopting a local moving equilibrium, $\hat{s}_{LR}$, $w_L^*$ and $w_R^*$ are thus given by (71b) and (81), and the proof of Lemma 6.3 yields

$$ (h^{\text{HLL}})^{\text{eq}}(\hat{u}^{\text{HLL}})^{\text{eq}} = \frac{h_L + h_R}{2} (\hat{u}^{\text{HLL}})^{\text{eq}} = q^{\text{eq}}. $$

Next, according to the proof of Lemma 6.4-(ii), the above equation infers

$$ \frac{g}{8} \frac{1}{1 + (F_{LR}^2)^{\text{eq}}} \left(1 + (F_{LR}^2)^{\text{eq}}[h] + [z]\right)^2 = \left(\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}) + \frac{g}{8}[z]^2 - \frac{g[z]^2}{8} \frac{F_{LR}^2}{1 + F_{LR}^2}\right)^{\text{eq}}, \\ = \frac{(\hat{\eta}^{\text{HLL}})^{\text{eq}} - \hat{\eta}(\hat{w}^{\text{HLL}})^{\text{eq}}}{1 + (F_{LR}^2)^{\text{eq}}} + \frac{g}{8} \frac{[z]}{1 + (F_{LR}^2)^{\text{eq}}}. $$
Plugging this last equation into (82), we deduce that $[h^*]^{\text{eq}}$ verifies

$$ 
\begin{aligned}
[h^*]^{\text{eq}} &= -\frac{[z]}{1 + (F_{LR}^2)^{\text{eq}}} + \left(\text{sign}\left((1 + (F_{LR}^2)[h] + [z]\right)\right)^{\text{eq}} \sqrt{\frac{8}{g} \frac{(\hat{\eta}^{\text{HLL}})^{\text{eq}} - \hat{\eta}(\hat{w}^{\text{HLL}})^{\text{eq}}}{1 + (F_{LR}^2)^{\text{eq}}} + \frac{[z]^2}{(1 + (F_{LR}^2)^{\text{eq}})^2}}, \\
&= -\frac{[z]}{1 + (F_{LR}^2)^{\text{eq}}} + \text{sign}\left((1 + (F_{LR}^2)^{\text{eq}})[h] + [z]\right) \frac{|(1 + (F_{LR}^2)^{\text{eq}})[h] + [z]|}{|1 + (F_{LR}^2)^{\text{eq}}|}, \quad (85) \\
&= -\frac{[z] + (1 + (F_{LR}^2)^{\text{eq}})[h] + [z]}{1 + (F_{LR}^2)^{\text{eq}}}, \\
&= [h].
\end{aligned}
$$

Associating the above equality to (81a) and (81c), we eventually have $(h_L^*)^{\text{eq}} = h_L$ and $(h_R^*)^{\text{eq}} = h_R$.

Now, we have to show the equalities

$$ q^{\text{eq}} = (h_R^*)^{\text{eq}}(u_R^*)^{\text{eq}} = (h_L^*)^{\text{eq}}(u_L^*)^{\text{eq}}. $$

Since $(h_L^*)^{\text{eq}} = h_L$ and $(h_R^*)^{\text{eq}} = h_R$, the above equalities are equivalent to $(u_L^*)^{\text{eq}} = \frac{q^{\text{eq}}}{h_L}$ and $(u_R^*)^{\text{eq}} = \frac{q^{\text{eq}}}{h_R}$.

According to Lemma 6.1, we have $q_{LR}^{\text{eq}} = q^{\text{eq}}$ and from (81b) we get

$$ 
\begin{aligned}
(u_L^*)^{\text{eq}} &= (\hat{u}^{\text{HLL}})^{\text{eq}} + \frac{q_{LR}^{\text{eq}}[h^*]^{\text{eq}}}{2(h^{\text{HLL}})^{\text{eq}}\sqrt{h_L h_R}} \sqrt{\frac{(h_R^*)^{\text{eq}}}{(h_L^*)^{\text{eq}}}}, \\
&= \frac{2q^{\text{eq}}}{h_L + h_R} + \frac{q^{\text{eq}}[h]}{(h_L + h_R)h_L}, \\
&= q^{\text{eq}} \frac{2h_L + [h]}{(h_L + h_R)h_L}, \\
&= \frac{q^{\text{eq}}}{h_L} = u_L.
\end{aligned}
$$

Arguing (81d), a similar computation gives $(u_R^*)^{\text{eq}} = u_R$, and we deduce the well-balanced property.

Concerning the discrete entropy inequality, if $h_L^*$ and $h_R^*$ given by (81a) and (81c) are both positive then the formulations (84) coincide with (81). In this case, the discrete entropy inequality is a direct consequence of the Riemann solver definition and this inequality is shown in Lemma 6.2 which completes the proof. $\square$

To conclude this section, let us notice that the dry-wet transitions in the above theorem are addressed with several cases in the formulations (83) and (84). These case distinctions are robust (see [39], Sect. 3.1.2.4) but their formulations are not continuous with the numerical scheme defined in Lemma 6.2.

Moreover, the limitation technique $\min(\max(\cdot, \cdot), \cdot)$ used in (84) is an $\varepsilon$-parametrized version of the procedure defined in Sections 4 and 5. The parameter $\varepsilon > 0$ guarantees the inequalities $h_L^* > 0$ and $h_R^* > 0$ that are essential for solving the system (33)–(71) (see Lem. 6.2).

In addition, we remark that Theorem 6.1 shows the existence of an entropy-satisfying well-balanced numerical scheme for the moving equilibrium. This scheme generalizes the entropy-satisfying well-balanced scheme for the lake at rest established in Theorem 5.1. Indeed, as soon as we have $q_{LR} = 0$ then (84) degenerates toward (68). As a consequence, Theorem 6.1 unifies the numerical schemes introduced in Theorem 5.1 (well-balanced for the lake at rest) and in Lemma 4.2 ($z = \text{cst}$).

## 7. NUMERICAL RESULTS

For all test cases, we fix $g = 9.81$. Section 7.1 concerns the flat regions and we focus on the schemes given by Lemma 4.2. Section 7.2 illustrates the entropy-stable schemes 5.1 that are only well balanced for the lake at rest (10). Finally, the numerical tests of Section 7.3 deal with the fully well-balanced entropy-stable scheme described in Theorem 6.1.
## 7.1. Two states entropy-stable approximate Riemann solvers in the flat regions

In this section, we consider the two schemes described in Lemma 4.1 and their limited versions given by Lemma 4.2. For all interfaces having the states $w_L$, $w_R$ on either on its sides, the numerical artificial viscosity $\lambda > 0$ and the time step $\Delta t > 0$ are taken to be equal to

$$ \lambda = \max_{\alpha \in \{L,R\}} \left| u_{\alpha} \pm \sqrt{gh_{\alpha}} \right|, \quad \frac{\lambda \Delta t}{\Delta x} \le \frac{1}{2}, \quad \forall (L,R). \qquad (86) $$

The results are compared to the Suliciu relaxation scheme and to the solver proposed in Section b.ii of [29]. The two constants $C_1$ and $C_2$ required for the scheme in Section b.ii of [29] are fixed to $C_1 = C_2 = 10^{-7}$.

First, we illustrate the influence of the choice of the intermediate states $h^*$ given by Lemma 4.1. In this regard, we set

$$ \begin{array}{l} \text{EC}_1: \left\{ \begin{array}{l} h_R^* = h^{\text{HLL}} + \sqrt{\frac{2}{g}(\eta^{\text{HLL}} - \eta(w^{\text{HLL}}))}, \\ h_L^* = h^{\text{HLL}} - \sqrt{\frac{2}{g}(\eta^{\text{HLL}} - \eta(w^{\text{HLL}}))}. \end{array} \right. \\ \\ \text{EC}_2: \left\{ \begin{array}{l} h_R^* = h^{\text{HLL}} - \sqrt{\frac{2}{g}(\eta^{\text{HLL}} - \eta(w^{\text{HLL}}))}, \\ h_L^* = h^{\text{HLL}} + \sqrt{\frac{2}{g}(\eta^{\text{HLL}} - \eta(w^{\text{HLL}}))} \end{array} \right. \end{array} \qquad (87) $$

EC_3: random(EC_1, EC_2),

with $w^{\text{HLL}}$, $\eta^{\text{HLL}}$ given by (31a) and (31d) and where random(EC_1, EC_2) denotes a random choice between the two configurations EC_1 and EC_2. The domain $[-1, 1]$ is discretized with 400 cells and we consider the following initial condition:

$$ h_0(x) = \begin{cases} 3, & \text{if } x < 0.5, \\ 1, & \text{otherwise,} \end{cases} \quad u_0(x) = 0. \qquad (88) $$

We impose homogeneous Neumann boundary conditions on both sides. The exact solution consists of a rarefaction wave and a shock wave. The final time is 0.1. Figure 1 shows the results.

The three configurations $(EC_i)_{i \in \{1,2,3\}}$ provide very similar results. As a consequence, we only retain the EC_3 configuration. The second numerical test concerns the role of the limitation technique for the dry areas used in Lemma 4.2.

The domain $[-1, 1]$ is discretized with 400 cells and we consider the following initial condition:

$$ h_0(x) = 0.1, \quad u_0(x) = \begin{cases} 10 & \text{if } x < 0.5, \\ 0 & \text{otherwise.} \end{cases} \qquad (89) $$

We lay down homogeneous Neumann boundary conditions on both sides. The exact solution consists of two strong shock waves near the dry areas. The CFL condition is given by (86). The final time is 0.1. Figure 2 shows the results.

We observe very good agreement with the exact solution. Without the conservative limitation technique used in Lemma 4.2, this problem cannot be carried out with the solvers from Lemma 4.1.

## 7.2. Well-balanced two states entropy-stable approximate Riemann solvers

In this section, we consider the schemes of Theorem 5.1. For these solvers, the numerical artificial viscosity $\lambda > 0$ has to be selected to ensure the existence of solutions of the system made of (33) and (56). We adopt the following selection procedure. Starting from the equation (86), we increase $\lambda$ until the system made of (33) and (56) admits solutions. The time step $\Delta t > 0$ is then determined according to the standard CFL condition
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Numerical results at time 0.1 for the initial condition (88), with the legend ($EC_i$)$_{i \in \{1,2,3\}}$: Lemma 4.1 solvers with the configurations (87), Relax: Suliciu relaxation scheme and HLLSE: solver ([29], Sect. b.ii).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2. Numerical results at time 0.1 for the initial condition (89), with the legend $EC_3$: Lemma's 4.2 solvers, Relax: Suliciu relaxation and HLLSE: solver ([29], Sect. b.ii).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIGURE 3. On the left: numerical results at time 100.0 for the lake at rest problem (90) with the legend WBEC: entropy-satisfying well-balanced solvers for the lake at rest given by Theorem 5.1, HR: hydro-static reconstruction [3] with Rusanov numerical flux [44]. On the right: errors between the exact and the numerical solutions at time 100.0 for the variables $h$, $hu$.</figcaption></figure>

<table><thead><tr><th colspan="4">Errors on (h, hu) for the lake at rest.</th></tr></thead><tbody><tr><th colspan="4">h</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr><tr><th>WBEC</th><td>1.69E-21</td><td>4.79E-21</td><td>2.29E-41</td></tr><tr><th>HR</th><td>6.55E-18</td><td>8.44E-18</td><td>1.25E-17</td></tr></tbody><tbody><tr><th colspan="4">hu</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr><tr><th>WBEC</th><td>9.88E-20</td><td>3.43E-20</td><td>1.14E-22</td></tr><tr><th>HR</th><td>3.17E-17</td><td>2.52E-17</td><td>3.07E-17</td></tr></tbody></table>

$\frac{\lambda\Delta t}{\Delta x} \le \frac{1}{2}$. Unless otherwise stated, to run the following simulations, we randomly select one of two solvers defined in Theorem 5.1 and we compare the results to the standard hydro-static reconstruction [3] coupled to the standard Rusanov numerical flux [44].

The first experiment is devoted to a flow at rest with emerging bottom as introduced in [22]. The spatial domain $[0, 25]$ is discretized with 400 cells. The initial condition and the topography are given by

$$h_0(x) = \max(0.15 - z(x), 0), \quad u_0(x) = 0, \quad \text{with} \quad z(x) = \max(0, 0.2 - 0.05(x - 10)^2). \quad (90)$$

We prescribe periodic boundary conditions. The exact solution is a lake at rest equilibrium (10). The final time is 100. Figure 3 shows the results.

Thanks to the well-balanced property and to the transition toward dry areas, the steady state at rest (90) is preserved up to the machine precision.

The second numerical experiment concerns the three Goutal and Maurel test cases [26]. For these three experiments, the spatial domain $[0, 25]$ is discretized with 400 cells. Using the superscript $GM_1$, $GM_2$ and $GM_3$ to denote each problem, the initial conditions are

$$h_0^{GM_k}(x) = h^{GM_k}, \quad (hu)_0^{GM_k}(x) = q^{GM_k}, \quad \forall k \in \{1, 2, 3\}, \quad (91)$$

where $(h^{GM_k})_{k \in \{1,2,3\}}$ and $(q^{GM_k})_{k \in \{1,2,3\}}$ are given in Table 1. The bottom topography $z$ is given by (90). On the left boundary, the water height satisfies a homogeneous Neumann condition and the discharge $q = hu$ is set to $(q^{GM_k})_{k \in \{1,2,3\}}$. On the right boundary, the water height is set to $h^{GM_k}$ when the flow is subcritical and a homogeneous Neumann boundary condition is prescribed otherwise. The discharge follows a homogeneous Neumann boundary condition.

Such initial and boundary conditions provide a transient state followed by a steady state consisting of a uniform discharge. For $GM_1$ and $GM_2$, this steady state is continuous whereas $GM_3$ involves a stationary shock. The final times are given in Table 1. The exact solutions are computed using the software SWASHES [18] and Figure 4 shows the results.

The free surface $h+z$ given by our schemes may be misplaced: it is particularly obvious for the $GM_2$ problem. For the $GM_3$ problem, the free surface is undervalued before the stationary shock wave and it is sharp after. Despite the fully discrete entropy stability verified by the Theorem 5.1 schemes, the numerical solutions may converge to a non-admissible weak solution. Such arbitrary wrong convergences have already been observed
<table><caption>TABLE 1. Final times, initial values and boundary conditions for the Goutal and Maurel test cases 91.</caption><thead><tr><th colspan="4">Parameters used for the Goutal and Maurel test cases.</th></tr><tr><th></th><th>GM<sub>1</sub></th><th>GM<sub>2</sub></th><th>GM<sub>3</sub></th></tr></thead><tbody><tr><td>Final time</td><td>500</td><td>125</td><td>1000</td></tr><tr><td>Initial height h<sup>GM<sub>k</sub></sup></td><td>2</td><td>0.66</td><td>0.33</td></tr><tr><td>Boundary discharge q<sup>GM<sub>k</sub></sup></td><td>4.42</td><td>1.53</td><td>0.18</td></tr></tbody></table>

[13, 16, 42]. In addition, the schemes of Theorem 5.1 produce spurious oscillations for the variable $hu$. For these three test cases, the standard hydro-static reconstruction [3] coupled to the Rusanov numerical flux [44] is more effective.

The last test case of this section is devoted to the dam break problem as described in [41]. The domain $[-1, 1]$ is discretized with 400 cells. We consider the following initial condition and topography:

$$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). \tag{92}$$

We prescribe homogeneous Neumann boundary conditions on both boundaries. The final time is 0.1. A reference solution is computed using the standard HLL scheme [29] on a fine grid made of 50 000 cells and coupled to the discrete source term $(-gh\partial_x z)(x_i, t^n) \approx -g\frac{h_i^l+h_{i+1}^l}{2}\frac{z_{i+1}-z_i}{\Delta x}$. We also compare the influence of the intermediate states $h^*$ defined in Lemma 5.1. Hence, we set

$$ \text{WBEC}_1: \left\{ \begin{aligned} h_R^* &= h^{\text{HLL}} - \frac{(z_R - z_L) - \sqrt{8 \frac{(\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}))}{g} + (z_R - z_L)^2}}{2}, \\ h_L^* &= h^{\text{HLL}} + \frac{(z_R - z_L) - \sqrt{8 \frac{(\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}))}{g} + (z_R - z_L)^2}}{2}, \end{aligned} \right. \tag{93} \\ \text{WBEC}_2: \left\{ \begin{aligned} h_R^* &= h^{\text{HLL}} - \frac{(z_R - z_L) + \sqrt{8 \frac{(\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}))}{g} + (z_R - z_L)^2}}{2}, \\ h_L^* &= h^{\text{HLL}} + \frac{(z_R - z_L) + \sqrt{8 \frac{(\hat{\eta}^{\text{HLL}} - \hat{\eta}(\hat{w}^{\text{HLL}}))}{g} + (z_R - z_L)^2}}{2}, \end{aligned} \right. \\ \text{WBEC}_3: \text{random}(\text{WBEC}_1, \text{WBEC}_2), \end{aligned} $$

with $w^{\text{HLL}}$, $\hat{\eta}^{\text{HLL}}$ given by (31a)-(31e) and where random $(\text{WBEC}_1, \text{WBEC}_2)$ denotes a random choice between the two configurations $\text{WBEC}_1$, $\text{WBEC}_2$. Figure 5 displays the results.

For all configurations (93), we observe very good agreement to the reference solution. The shock wave is sharper with the schemes detailed in Theorem 5.1.

### 7.3. Fully well-balanced entropy-stable two states approximate Riemann solvers

In this section, we consider the fully well-balanced entropy-stable scheme defined in Theorem 6.1. For this scheme, the numerical artificial viscosity $\lambda > 0$ and the time step $\Delta t > 0$ have to be selected to guarantee the existence of solutions of the system made of (33) and (71). To select such a couple, we initialize $\lambda > 0$ with the value given by the equation (86) and we set $\Delta t > 0$ according to $\frac{\lambda\Delta t}{\Delta x} = \frac{1}{2}$. Then, we increase $\lambda$ and decrease $\Delta t$.
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Numerical results for the Goutal and Maurel problems (91) on a mesh composed of 400 cells. The legend is WBEC: entropy-satisfying well-balanced solvers for the lake at rest given by Theorem 5.1, HR: hydro-static reconstruction [3] with the standard Rusanov numerical flux [44].</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Numerical results at time 0.1 for the dam break problem (92) on a mesh composed of 400 cells with the legend (WBEC<sub>i</sub>)<sub>i∈{1,...,3}</sub>: Theorem 5.1 solvers with the configurations (93), HR: hydro-static reconstruction [3] with the Rusanov numerical flux [44].</figcaption></figure>

until the system composed of (33) and (71) admits solutions. For the sign function required in equation (82) we use the following regularized version:

$$ \mathrm{sign}(r) \approx \frac{r}{|r| + \zeta_{LR}/(g\bar{h})}, \qquad (94) $$
<figure><img src="image_8.png" /><figcaption>FIGURE 6. On the left: numerical results at time 1.0 for the steady state problem (95) with the legend FWBEC: fully well-balanced entropy-satisfying solver given by Theorem 6.1, BM: scheme derived in [14] associated to the Rusanov numerical flux [44], FMT: scheme derived in Section 3.2 of [20]. On the right: errors between the exact and the numerical solutions at time 1.0 for the variables $h, hu$.</figcaption></figure>

<table><thead><tr><th colspan="4">Errors on (h, hu) for the steady state (95).</th></tr></thead><tbody><tr><th colspan="4">h</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr><tr><th>FWBEC</th><td>0.00E-00</td><td>0.00E-00</td><td>0.00E-00</td></tr><tr><th>BM</th><td>2.49E-17</td><td>4.41E-17</td><td>1.27E-17</td></tr><tr><th>FMT</th><td>1.84E-13</td><td>2.03E-13</td><td>3.28E-13</td></tr></tbody><tbody><tr><th colspan="4">hu</th></tr><tr><th></th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr><tr><th>FWBEC</th><td>1.13E-20</td><td>8.74E-20</td><td>7.64E-39</td></tr><tr><th>BM</th><td>8.58E-17</td><td>1.29E-16</td><td>9.19E-17</td></tr><tr><th>FMT</th><td>2.87E-13</td><td>3.75E-13</td><td>4.10E-13</td></tr></tbody></table>

where $\zeta_{LR}$ and $\bar{h}$ are defined by (70). We compare the schemes of Theorem 6.1 to those detailed in Section 3.2 of [20] and in [14] that will be denoted FMT and BM. The scheme BM is associated with the standard Rusanov numerical flux [44]. These two schemes are chosen for the comparative because, like the scheme 6.1, they are well-balanced for equilibria with non-null velocities. In addition, both satisfy a semi-discrete entropy inequality and BM preserves the convex domain $\Omega$.

The first test case examines the fully well-balanced property. We consider the domain $[0, 1]$ discretized with 400 cells. The initial conditions and the bottom topography verify

$$ (hu)_0(x) = q_0, \quad \frac{u_0^2(x)}{2} + g(h_0 + z)(x) = B_0, \quad \text{with} \quad z(x) = \frac{5}{2} \cos^2(4\pi x), \qquad (95) $$

where $q_0 = \frac{5}{2}$, $B_0 = \frac{25}{98} + 4g$. The initial values $(h_0, u_0)(x)$ are computed from the equations (95) using a Newton method. We impose periodic boundary conditions on both sides. The exact solution is a moving steady state (9). The final time is 1.0. The results are reported in Figure 6. We also report the errors between the numerical and the exact solution for several norms.

Thanks to the well-balanced property, the moving equilibrium is preserved up to the machine precision.

Now we repeat the three Goutal and Maurel test cases [26] and we refer to Section 7.2 for the details of the test cases. Figure 7 shows the results and Figure 8 show the errors between the numerical and exact solutions for each problem and for several norms.

Despite our efforts, we were not able to run these three problems with the scheme FMT. The scheme of Theorem 6.1 scheme generates spurious oscillations, in particular for the variable $hu$. Nevertheless, the wrong convergences observed for the $GM_2$ and $GM_3$ problems with the schemes derived in Theorem 5.1 and tested in Section 7.2, do not occur with the fully well-balanced scheme in Figure 7. However, the free surface for the $GM_1$ problem is, once again, slightly misplaced.

Regardless of its fully well-balanced property, the scheme derived in Theorem 6.1 does not reach the exact solution up to the machine precision. This defect could be due to the implementation of the sign function required in equation (82). The choice proposed in the equation (94) is a smooth version of the sign function but other versions are possible and each of them could influence the result. For the three Goutal and Maurel problems the BM scheme is more accurate.

Our next test case concerns the dam break problem as done in Section 7.2. The initial condition and the bottom topography are given by the equations (92). The final time is 0.1 and a reference solution is computed
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Numerical results for the Goutal and Maurel's problems (91) on a mesh composed of 400 cells. The legend is FWBEC: fully well-balanced entropy-satisfying solver given by Theorem 6.1, BM: scheme derived in [14] associated to the Rusanov numerical flux [44].</figcaption></figure>
<table><thead><tr><th colspan="4">Errors on (h, hu) for the GM<sub>1</sub> problem.</th></tr><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>FWBEC</th><td>3.86E-02</td><td>1.41E-02</td><td>2.00E-04</td></tr><tr><th>BM</th><td>1.91E-15</td><td>4.16E-16</td><td>6.11E-17</td></tr></tbody></table><table><thead><tr><th colspan="4">Errors on (h, hu) for the GM<sub>2</sub> problem.</th></tr><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>FWBEC</th><td>9.35E-02</td><td>2.13E-02</td><td>4.46E-03</td></tr><tr><th>BM</th><td>9.06E-03</td><td>2.37E-03</td><td>5.12E-04</td></tr></tbody></table><table><thead><tr><th colspan="4">Errors on (h, hu) for the GM<sub>3</sub> problem.</th></tr><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>FWBEC</th><td>4.28E-02</td><td>5.11E-02</td><td>2.62E-03</td></tr><tr><th>BM</th><td>8.23E-03</td><td>2.20E-02</td><td>4.86E-04</td></tr></tbody></table><table><thead><tr><th colspan="4">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>FWBEC</th><td>2.57E-02</td><td>2.00E-02</td><td>1.59E-03</td></tr><tr><th>BM</th><td>4.21E-03</td><td>1.15E-02</td><td>1.32E-04</td></tr></tbody></table>

<table><thead><tr><th colspan="4">Errors on (h, hu) for the GM<sub>2</sub> problem.</th></tr><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>FWBEC</th><td>9.35E-02</td><td>2.13E-02</td><td>4.46E-03</td></tr><tr><th>BM</th><td>9.06E-03</td><td>2.37E-03</td><td>5.12E-04</td></tr></tbody></table><table><thead><tr><th colspan="4">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>FWBEC</th><td>4.89E-02</td><td>3.53E-02</td><td>1.25E-03</td></tr><tr><th>BM</th><td>7.78E-05</td><td>3.38E-05</td><td>4.33E-08</td></tr></tbody></table>

<table><thead><tr><th colspan="4">Errors on (h, hu) for the GM<sub>3</sub> problem.</th></tr><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>FWBEC</th><td>4.28E-02</td><td>5.11E-02</td><td>2.62E-03</td></tr><tr><th>BM</th><td>8.23E-03</td><td>2.20E-02</td><td>4.86E-04</td></tr></tbody></table><table><thead><tr><th colspan="4">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>FWBEC</th><td>2.57E-02</td><td>2.00E-02</td><td>1.59E-03</td></tr><tr><th>BM</th><td>4.21E-03</td><td>1.15E-02</td><td>1.32E-04</td></tr></tbody></table>

FIGURE 8. Errors between the numerical and the exact solutions for each Goutal and Mau- rel problems (91) and for several norms. The legend is FWBEC: fully well-balanced entropy- satisfying solver given by Theorem 6.1, BM: scheme derived in [14] coupled to the Rusanov numerical flux [44]

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

FIGURE 9. Numerical results at time 0.1 on a mesh composed of 400 cells for the dam break problem (92) with the legend FWBEC: Theorem 6.1 solver, BM: scheme derived in [14] asso- ciated with the Rusanov numerical flux [44], FMT: scheme derived in Section 3.2 of [20].

using the standard HLL scheme [29] on a fine grid made of 50 000 cells 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 9 displays the results.

The scheme Section 3.2 of [20] generates spurious oscillations near the shock wave and we failed to achieve better results. The scheme [14] associated to the Rusanov numerical flux [44] preserves the initial condition and
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Numerical results at time 0.05 on a mesh composed of 400 cells for the problem (96) with the legend FWBEC: Theorem 6.1 solver, BM: scheme derived in [14] associated with the Rusanov numerical flux [44], FMT: scheme derived in Section 3.2 of [20]. At the bottom, numerical results with a mesh refinement for the FWBEC scheme. For the clarity, the display is zoomed in on the shock wave near $x = 1.1$.</figcaption></figure>

captures two steady states at rest. For the scheme of Theorem 6.1, we observe a good agreement to the reference solution.

We now consider a test case in which the fluid is close to the regime $\frac{|u|}{\sqrt{gh}} \approx 1$. The domain $[-\frac{3}{2}, \frac{3}{2}]$ is discretized with 400 cells and we consider the following data:

$$h_0(x) = \sqrt[3]{\frac{(q^{\text{eq}})^2}{g}} \cdot \begin{cases} 4 + \sqrt{14}, & \text{if } x < 0, \\ 4 - \sqrt{14}, & \text{otherwise,} \end{cases} \quad (hu)_0(x) = q^{\text{eq}}, \quad z(x) = 1 - \tanh(10x), \quad (96)$$

with $q^{\text{eq}} = 10$. At $x = 0$, the above initial data involves a mean Froude number equals to one. On both space domain boundary, we enforce Dirichlet boundary conditions according to the initial values given by (96). A reference solution is computed using the standard HLL scheme [29] on a fine grid made of 50000 cells and 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}$. The final time is 0.05 and Figure 10 displays the results.

In the smooth regions, the above numerical results for the schemes BM, FMT, FWBEC are similar but more diffuse for FWBEC. Near the shock wave approximately localized in $x = 1.1$, the FMT scheme is more accurate but it generates some spurious oscillations. These oscillations do not appear with the FWBEC and BM schemes but the numerical solution is very diffusive.
<figure><img src="image_12.png" /><figcaption>FIGURE 11. Numerical results at time 0.05 on a mesh composed of 400 cells for the problem (97) with the legend FWBEC: Theorem 6.1 solver, BM: scheme derived in [14] associated with the Rusanov numerical flux [44], FMT: scheme derived in Section 3.2 of [20]. At the bottom, numerical results with a mesh refinement for the FWBEC scheme. For the clarity, the display is zoomed in on the discontinuity near $x = 0$.</figcaption></figure>

We conclude this section with the numerical simulation of problem (96) endowed with a discontinuous topography function $z$. Namely, we once again consider the domain $[-\frac{3}{2}, \frac{3}{2}]$ discretized with 400 cells, but the data are now as follows:

$$h_0(x) = \sqrt[3]{\frac{(q^{\text{eq}})^2}{g}} \cdot \begin{cases} 4 + \sqrt{14}, & \text{if } x < 0, \\ 4 - \sqrt{14}, & \text{otherwise,} \end{cases} \quad (hu)_0(x) = q^{\text{eq}}, \quad z(x) = \begin{cases} 2, & \text{if } x < 0, \\ 0, & \text{otherwise,} \end{cases} \quad (97)$$

with $q^{\text{eq}} = 10$. On both boundaries, Dirichlet boundary conditions are imposed according to the initial values given in (97). A reference solution is computed with the standard HLL scheme [29] on a fine grid made of 50 000 cells and 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}$. The final time is 0.05. The results are presented in Figure 11.

Except near the discontinuity localized at $x = 0$, the above numerical results are very similar to those in Figure 10. However, the FWBEC scheme derived in Theorem 6.1 commits a large error near $x = 0$. The simulation of this discontinuity is more accurate and relevant with the FMT scheme or the BM scheme.
## 8. CONCLUSION

We have presented three explicit entropy-stable Godunov-type schemes for the shallow water equations. The first one concerns the flat regions, the second scheme is well-balanced for the lake at rest (10) and the third one is well-balanced for all smooth stationary solutions defined by (9).

The discrete entropy inequality is achieved from sufficient conditions used in the design of the scheme. These conditions lead to quadratic equations that are always well-posed under restrictions on the artificial viscosity and on the time step. These restrictions are implicit in the well-balanced schemes.

From a numerical point of view, the scheme devoted to the flat regions provides good results. The well-balanced schemes yield satisfactory results especially in the presence of shock waves. But, they may converge to weak solutions made of non-admissible stationary contact waves. The study of the reasons for this wrong convergence should be conducted.

## ACKNOWLEDGMENTS

The authors thank the ANR project MUFFIN ANR-19-CE46-0004 for the financial support and the France 2030 framework program Centre Henri Lebesgue ANR-11-LABX-0020-01 for its mathematical research programs. The authors also thank Mehdi Badsi for fruitfully discussions.

## DATA AVAILABILITY STATEMENT

The research data associated with this article are included in the article.

## REFERENCES

[1] F. Alcrudo and F. Benkhaldoun, Exact solutions to the Riemann problem of the shallow water equations with a bottom step. *Comput. Fluids* **30** (2001) 643–671.

[2] A. Aleksyuk, M. Malakhov and V. Belikov, The exact Riemann solver for the shallow water equations with a discontinuous bottom. *J. Comput. Phys.* **450** (2022) 110801.

[3] E. Audusse, F. Bouchut, M. Bristeau, R. Klein and B. Perthame, A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. *SIAM J. Sci. Comput.* **25** (2004) 2050–2065.

[4] E. Audusse, C. Chalons and P. Ung, A simple well-balanced and positive numerical scheme for the shallow-water system. *Commun. Math. Sci.* **13** (2015) 1317–1332.

[5] R. Bernetti, V. Titarev and E. Toro, Exact solution of the Riemann problem for the shallow water equations with discontinuous bottom geometry. *J. Comput. Phys.* **227** (2008) 3212–3243.

[6] C. Berthon, Numerical approximations of the 10-moment gaussian closure. *Math. Comput.* **75** (2006) 1809–1832.

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

[8] C. Berthon and V. Michel-Dansac, A simple fully well-balanced and entropy preserving scheme for the shallow-water equations. *Appl. Math. Lett.* **86** (2018) 284–290.

[9] C. Berthon, A. Crestetto and F. Foucher, A well-balanced finite volume scheme for a mixed hyperbolic/parabolic system to model chemotaxis. *J. Sci. Comput.* **67** (2016) 618–643.

[10] 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 Vol.* (2021) https://hal.science/hal-03192954/.

[11] A. Bihlo and S. MacLachlan, Well-balanced mesh-based and meshless schemes for the shallow-water equations. *BIT Numer. Math.* **58** (2018) 579–598.

[12] 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).

[13] F. Bouchut and T.M. de Luna, An entropy satisfying scheme for two-layer shallow water equations with uncoupled treatment. *ESAIM: Math. Model. Numer. Anal.* **42** (2008) 683–698.

[14] F. Bouchut and T.M. de Luna, A subsonic-well-balanced reconstruction scheme for shallow water flows. *SIAM J. Numer. Anal.* **48** (2010) 1733–1758.

[15] M. Castro Díaz, J. López-García and C. Parés, High-order exactly well-balanced numerical methods for shallow water systems. *J. Comput. Phys.* **246** (2013) 242–264.
[16] A. Chinnayya, A. Leroux and N. Seguin, A well-balanced numerical scheme for the approximation of the shallow-water equations with topography: the resonance phenomenon. *Int. J. Finite 1* (2004) 33.

[17] F. Coquel, J.-M. Hérard and K. Saleh, A splitting method for the isentropic Baer-Nunziato two-phase flow model, in ESAIM: Proceedings. Vol. 38. EDP Sciences (2012) 241-256.

[18] O. Delestre, C. Lucas, P. Ksinant, F. Darboux, C. Laguerre, T. Vo, F. James and S. Cordier, SWASHES: a compila- tion of shallow water analytic solutions for hydraulic and environmental studies. *Int. J. Numer. Methods Fluids* **72** (2013) 269-300.

[19] A. Duran, J. Vila and R. Baraille, Semi-implicit staggered mesh scheme for the multi-layer shallow water system. *C. R. Acad. Sci. Paris* **355** (2017) 1298-1306.

[20] U. Fjordholm, S. Mishra and E. Tadmor, Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography. *J. Comput. Phys.* **230** (2011) 5587-5609.

[21] G. Gallice, Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source. *C. R. Math.* **334** (2002) 713-716.

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

[23] T. Gallouët, R. Herbin, J.C. Latché and N. Therme, Consistent internal energy based schemes for the compressible Euler equations, in Numerical Simulation in Physics and Engineering: Trends and Applications: Lecture Notes of the XVIII "Jacques-Louis Lions" Spanish-French School (2021) 119-154.

[24] L. Gastaldo, R. Herbin, J.-C. Latché and N. Therme, A MUSCL-type segregated-explicit staggered scheme for the Euler equations. *Comput. Fluids* **175** (2018) 91-110.

[25] D. George, Augmented Riemann solvers for the shallow water equations over variable topography with steady states and inundation. *J. Comput. Phys.* **227** (2008) 3089-3113.

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

[27] E. Han and G. Warnecke, Exact Riemann solutions to shallow water equations. *Q. Appl. Math.* **72** (2014) 407-453.

[28] A. Harten and P.D. Lax, A random choice finite difference scheme for hyperbolic conservation laws. *SIAM J. Numer. Anal.* **18** (1981) 289-315.

[29] A. Harten, P. Lax and B. Van Leer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. *SIAM Rev.* **25** (1983) 35-61.

[30] R. Herbin, J. Latché and T. Nguyen, Consistent explicit staggered schemes for compressible flows part I: the barotropic Euler equations (2013) https://hal.science/hal-00821069/.

[31] A. Hiltebrand and S. Mishra, Entropy stable shock capturing space-time discontinuous Galerkin schemes for systems of conservation laws. *Numer. Math.* **126** (2014) 103-151.

[32] A. Hiltebrand and S. Mishra, Entropy stability and well-balancedness of space-time DG for the shallow water equations with bottom topography. *Netw. Heterog. Media* **11** (2016) 145-162.

[33] A. Hiltebrand, S. Mishra and C. Parés, Entropy-stable space-time DG schemes for non-conservative hyperbolic systems. *ESAIM: Math. Model. Numer. Anal.* **52** (2018) 995-1022.

[34] 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. *Comput. Fluids* **82** (2013) 132-147.

[35] J. Koellermeier and E. Pimentel-García, Steady states and well-balanced schemes for shallow water moment equa- tions with topography. *Appl. Math. Comput.* **427** (2022) 127166.

[36] P. LeFloch and M. Thanh, A Godunov-type method for the shallow water equations with discontinuous topography in the resonant regime. *J. Comput. Phys.* **230** (2011) 7631-7660.

[37] M. Li, R. Mu and H. Dong, A well-balanced discontinuous Galerkin method for the shallow water flows on erodible bottomimage 1. *Comput. Math. App.* **119** (2022) 13-20.

[38] Q. Liang and F. Marche, Numerical resolution of well-balanced shallow water equations with complex source terms. *Adv. Water Res.* **32** (2009) 873-884.

[39] V. Michel-Dansac, *Development of high-order well-balanced schemes for geophysical flows*. Ph.D. thesis, Université de Nantes, Faculté des sciences et des techniques (2016).

[40] V. Michel-Dansac, C. Berthon, S. Clain and F. Foucher, A well-balanced scheme for the shallow-water equations with topography. *Comput. Math. App.* **72** (2016) 568-593.

[41] V. Michel-Dansac, C. Berthon, S. Clain and F. Foucher, A well-balanced scheme for the shallow-water equations with topography or manning friction. *J. Comput. Phys.* **335** (2017) 115-154.
[42] C. Parés and E. Pimentel, The Riemann problem for the shallow water equations with discontinuous topography: the wet-dry case. *J. Comput. Phys.* **378** (2019) 344–365.

[43] A. Prieto-Arranz, L. Ramirez, I. Couceiro, I. Colominas and X. Nogueira, A well-balanced sph-ale scheme for shallow water applications. *J. Sci. Comput.* **88** (2021) 1–29.

[44] V. Rusanov, The calculation of the interaction of non-stationary shock waves and obstacles. *USSR Comput. Math. Math. Phys.* **1** (1962) 304–320.

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

[46] Y. Xing, X. Zhang and C. Shu, Positivity-preserving high-order well-balanced discontinuous Galerkin methods for the shallow water equations. *Adv. Water Res.* **33** (2010) 1476–1493.

[47] Q. Zhang and W. Sheng, A note on the Riemann problem for shallow water equations with discontinuous topography. *Appl. Math. Lett.* **116** (2021) 107042.

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

Please help to maintain this journal in open access!

This journal is currently published in open access under the **Subscribe to Open model (S2O)**. We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.