Victor Michel-Dansac
email: victor.michel-dansac@univ-nantes.fr

Christophe Berthon
email: christophe.berthon@univ-nantes.fr

Stéphane Clain
email: clain@math.uminho.pt

Françoise Foucher
email: francoise.foucher@ec-nantes.fr

A well-balanced scheme for the shallow-water equations with topography or Manning friction

Keywords: shallow-water equations, Manning friction, Godunov-type schemes, well-balanced schemes, moving steady states 2000 MSC: 65M08, 65M12

We consider the shallow-water equations with Manning friction or topography, as well as a combination of both these source terms. The main purpose of this work concerns the derivation of a non-negativity preserving and well-balanced scheme that approximates solutions of the system and preserves the associated steady states, including the moving ones. In addition, the scheme has to deal with vanishing water heights and transitions between wet and dry areas. To address such issues, a particular attention is paid to the study of the steady states related to the friction source term. Then, a Godunov-type scheme is obtained by using a relevant average of the source terms in order to enforce the required well-balance property. An implicit treatment of both topography and friction source terms is also exhibited to improve the scheme while dealing with vanishing water heights. A second-order well-balanced MUSCL extension is designed, as well as an extension for the two-dimensional case. Numerical experiments are performed in order to highlight the properties of the scheme.

Introduction

The goal of this paper is to derive a numerical scheme to approximate the solutions of the shallow-water equations with topography and Manning friction. The Manning friction was introduced in [START_REF] Manning | On the flow of water in open channels and pipes[END_REF] (see also [START_REF] Delestre | SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies[END_REF] for an overview of other friction models, and [START_REF] Berthon | An efficient scheme on wet/dry transitions for shallow water equations with friction[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction[END_REF] for related works). The equations of interest consist in the following system:    ∂ t h + ∂ x q = 0,

∂ t q + ∂ x q 2 h + 1 2 gh 2 = -gh∂ x Z - kq|q| h η .
(1.1)

This system is used to model the flow of water in a one-dimensional channel, with a non-flat bottom that applies a friction force on the water. The variables involved in this model are the non-negative water height h(x, t) and the depthaveraged discharge of the water q(x, t). The known quantities are the shape of the channel bottom Z(x), the gravity constant g, the Manning friction coefficient k and a parameter η, equal to 7 3 . We also define the velocity u of the water, such that q = hu.

We define the admissible states space by Ω = W = t (h, q) ∈ R 2 ; h ≥ 0, q ∈ R .

Let us note that the water height may vanish, which accounts for dry areas. By convention, we impose u = 0 as soon as h = 0.

In order to shorten the notations, the system (1.1) is rewritten under a simpler form, as follows:

∂ t W + ∂ x f (W ) = s(W ), W ∈ Ω, (1.2)
where

W = h q , f (W ) =   q q 2 h + 1 2 gh 2   , s(W ) = 0 -gh∂ x Z - kq|q| h η . (1.3)
Omitting the source terms, the homogeneous system deriving from (1.1) is a hyperbolic system. Its characteristic velocities are given by uc and u + c (see [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Leveque | Numerical methods for conservation laws[END_REF] for instance), where c is the sound speed, defined as follows: c = gh.

(1.4)

The main focus of this paper lies in the study of the solutions of (1.2) that satisfy ∂ t W = 0. Such solutions are called steady state solutions. With vanishing partial derivatives in time in (1.1), we get:

   ∂ x q = 0, ∂ x q 2 h + 1 2 gh 2 = -gh∂ x Z - kq|q| h η .
(1.5)

We immediately obtain that, as per the first equation, the discharge q must be uniform. From now on, we denote this uniform discharge as follows:

q = q 0 . (1.6)
Using q 0 within the second equation of (1.5) yields:

∂ x q 2 0 h + 1 2 gh 2 = -gh∂ x Z - kq 0 |q 0 | h η . (1.7)
The solution of the nonlinear ordinary differential equation (1.7) is the water height h for a steady state with uniform discharge q 0 . For a general topography and a nonzero friction contribution, if q 0 = 0, we cannot solve this equation to obtain an analytical expression for h. Thus, we will focus on specific cases; namely, the steady states at rest, and the steady states with topography or friction only.

The most important and most extensively studied steady state of the shallowwater equations is the lake at rest steady state (see for instance the nonexhaustive list [START_REF] Audusse | A simple well-balanced and positive numerical scheme for the shallow-water system[END_REF][START_REF] Berthon | Efficient well-balanced hydrostatic upwind schemes for[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Chertock | Well-balanced positivity preserving central-upwind scheme for the shallow water system with friction terms[END_REF][START_REF] Gallardo | On a wellbalanced high-order finite volume scheme for the shallow water equations with bottom topography and dry areas[END_REF][START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF]). The lake at rest is obtained by assuming that the water is at rest, i.e. q 0 = 0. Thus, this steady state is defined as follows:

q = 0, h + Z = cst .
(1.8)

The importance of this steady state has led to the derivation of many schemes ensuring that the lake at rest is preserved. For instance, in [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF][START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF], the authors introduced the notion of a well-balanced scheme, that is to say a scheme that preserves some or all the steady states. Next, the well-balance approach was extended to nonlinear systems in [START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF]. The techniques from [START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF] were then simplified in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF], to obtain the hydrostatic reconstruction (see also [START_REF] Berthon | Efficient well-balanced hydrostatic upwind schemes for[END_REF][START_REF] Fernández-Nieto | A consistent intermediate wave speed for a well-balanced HLLC solver[END_REF][START_REF] Fjordholm | Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography[END_REF][START_REF] Jin | A steady-state capturing method for hyperbolic systems with geoterms[END_REF][START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF][START_REF] Lukáčová-Medvi ' Dová | Well-balanced finite volume evolution Galerkin methods for the shallow water equations[END_REF][START_REF] Natalini | A well-balanced numerical scheme for a one dimensional quasilinear hyperbolic model of chemotaxis[END_REF]). Now, let us consider a vanishing friction source term (i.e. k = 0) without being at rest (i.e. q 0 = 0). From (1.7), we get the moving steady states, governed by

∂ x q 2 0 h + 1 2 gh 2 = -gh∂ x Z.
(1.9)

The equation (1.9) that defines the moving steady states is known to be harder to satisfy at the discrete level than the lake at rest. Indeed, the lake at rest explicitly yields h with respect the the known quantity Z, while the water height is given by the nonlinear ODE (1.9) in the case of the moving steady states. Such steady states have been studied in the past, for instance in [START_REF] Chinnayya | A well-balanced numerical scheme for the approximation of the shallow-water equations with topography: the resonance phenomenon[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF], where the authors exhibit several properties of the topography steady state solutions. Such steady states have also been studied to provide several properties of numerical schemes, for instance in [START_REF] Goutal | Proceedings of the 2 nd Workshop on Dam-Break Wave Simulation[END_REF], where the authors define some specific moving steady states and use them to test the accuracy of non-well-balanced schemes. Then, in [START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF], the author introduced the notion of fully well-balanced schemes, able to preserve every steady state, and designed such a scheme. This approach was extended in [START_REF] Berthon | A fully well-balanced, positive and entropysatisfying Godunov-type method for the shallow-water equations[END_REF], where the authors derive a fully well-balanced and entropy-stable Godunov-type scheme, based on the resolution of a Bernoullitype equation. Since this equation is difficult to solve, the authors of [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF] have proposed a fully well-balanced Godunov-type scheme that uses a linearization of the Bernoulli-type equation. We mention two additional fully well-balanced numerical schemes based on different techniques: in [START_REF] Castro | Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique[END_REF], the authors generalize the hydrostatic reconstruction mentioned above, while the authors of [START_REF] Castro Díaz | High order exactly well-balanced numerical methods for shallow water systems[END_REF] propose a high-order approach using a relevant reconstruction. High-order techniques were also developed to ensure a very accurate approximation of the moving steady states (see for instance [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF][START_REF] Russo | High order well balanced schemes for systems of balance laws[END_REF][START_REF] Xing | Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium[END_REF][START_REF] Xing | On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations[END_REF]).

Let us now assume a flat topography and a nonzero discharge. The steady states governed by the friction source term are given as follows:

∂ x q 2 0 h + 1 2 gh 2 = - kq 0 |q 0 | h η .
(1.10) Some research has been devoted to Manning friction terms (for instance, see [START_REF] Delestre | SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies[END_REF][START_REF] Berthon | An efficient scheme on wet/dry transitions for shallow water equations with friction[END_REF]), but these works do not focus on steady states and the resulting schemes are not well-balanced.

In the present paper, we exhibit smooth steady solutions of (1.10) and propose a well-balanced scheme, to exactly preserve the lake at rest (1.8), the moving steady states (1.9), and the dominant friction steady states (1.10). This scheme will be derived in order to satisfy some essential properties, namely the exact preservation of all the steady states for the system (1.1) with topography or friction and the preservation of the non-negativity of the water height. In addition, the scheme will be designed to model transitions between wet areas, where h = 0 and dry areas, where h = 0.

The paper is organized as follows. First, in Section 2, we study smooth steady states for the shallow-water equations with friction over a flat bottom. Then, in Section 3, we develop a Godunov-type scheme that is well-balanced for the system (1.1). In this section, we first recall some generalities on Godunovtype schemes, then we derive the well-balanced scheme, according to the friction steady states. Afterwards, in Section 4, we extend this scheme to be wellbalanced for both the topography and the friction source terms. In addition, in this section, we prove the non-negativity preservation. Then, in Section 5, we propose an implicit treatment of the source terms, designed to deal with transitions between wet and dry areas. The scheme is then extended in Section 6 to be second-order accurate. This extension breaks the well-balance property, but we show a convex combination technique that allows the recovery of this essential property. Finally, in Section 7, we propose three kinds numerical experiments, namely well-balance assessments, dam-break simulations, and two dimensional simulations. The paper is concluded in Section 8.

Steady states with friction and a flat bottom

This section is devoted to studying the steady states coming from the friction source term, given by (1.5) with a flat topography (Z = cst), or equivalently by (1.10). Recall that, for steady states and after (1.6), the discharge q = q 0 is uniform throughout the domain. Here, we assume q 0 = 0. Indeed, we assume that the friction contribution vanishes as soon as the discharge is equal to zero. We rewrite (1.10) under the form:

q 2 0 ∂ x 1 h + g 2 ∂ x h 2 = - kq 0 |q 0 | h η . (2.1)
We are looking for smooth solutions h(x) > 0 of the equation (2.1). As a first step, we integrate this equation. Multiplying by h η provides:

- q 2 0 η -1 ∂ x h η-1 + g η + 2 ∂ x h η+2 = -kq 0 |q 0 |. (2.2)
Let us consider x 0 ∈ R as an arbitrary reference point, and introduce the initial condition h(x 0 ) = h 0 . For all x ∈ R, we integrate (2.2) over (x 0 , x), to get

- q 2 0 η -1 h η-1 -h η-1 0 + g η + 2 h η+2 -h η+2 0 + kq 0 |q 0 | (x -x 0 ) = 0. (2.3)
To shorten the notations, we rewrite (2.3) under the following form: ξ(h; x, h 0 , q 0 , x 0 ) = 0, (

where we have set ξ(h; x, h 0 , q 0 , x 0 ) = -

q 2 0 η -1 h η-1 -h η-1 0 + g η + 2
h η+2 -h η+2 0 + kq 0 |q 0 | (xx 0 ) .

(2.5)

We consider x a parameter so that (2.4) is an equation to characterize h := h(x). In order to exhibit such solutions, we first study ξ with a fixed x ∈ R. The derivative of ξ with respect to h is given by

∂ξ ∂h (h; x, h 0 , q 0 , x 0 ) = gh η-2 - q 2 0 g + h 3 .
Let us define a critical water height h c , such that ∂ξ ∂h (h c ; x, h 0 , q 0 , x 0 ) = 0, as follows:

h c = q 2 0 g 1 3
.

(2.6)

We easily obtain that the function h → ξ(h; x, h 0 , q 0 , x 0 ) is strictly increasing on (h c , +∞), while it is strictly decreasing on (0, h c ). As a consequence, it admits a unique minimum in (0, +∞), and this minimum is reached for h = h c . Now, we evaluate the sign of ξ. After straightforward computations, the following limits are proven to be satisfied by ξ:

• lim h→+∞ ξ(h; x, h 0 , q 0 , x 0 ) = +∞,

• lim h→0 + ξ(h; x, h 0 , q 0 , x 0 ) = ξ (x), where we have set

ξ (x) = gh η-1 0 h 3 c η -1 - h 3 0 η + 2 + kq 0 |q 0 |(x -x 0 ).
Note that ξ is finite. Moreover, let us introduce

ξ c (x) := ξ(h c ; x, h 0 , q 0 , x 0 ).
We note that

ξ c (x) = ξ (x) - 3gh η+2 c (η -1)(η + 2) < ξ (x).
Equipped with these new notations and these properties of ξ, we state the following result.

Lemma 1. Assume h > 0 and q 0 = 0. Thus, h c > 0 according to (2.6).

(i) If ξ (x) < 0, then the equation (2.4) admits a unique solution h sub (x).

Moreover, this solution belongs to (h c , +∞).

(ii) If ξ (x) > 0 and ξ c (x) > 0, then there is no solution to the equation (2.4).

(iii) If ξ (x) > 0 and ξ c (x) < 0, then the equation (2.4) admits two solutions, h sup (x) and h sub (x). One solution, h sup (x), stays within (0, h c ), while the second one, h sub (x), is in (h c , +∞). If ξ c (x) = 0, the equation (2.4) admits h c as a unique double solution.

Proof. To establish (i), we use that h → ξ(h; x, h 0 , q 0 , x 0 ) is strictly increasing on (h c , +∞) and that ξ c (x) < ξ (x) < 0, to deduce that the equation (2.4) admits a unique solution, which belongs to (h c , +∞).

Next, we turn to proving (ii). If ξ c (x) > 0, since ξ reaches its unique minimum ξ c (x) for h = h c , it is immediate that there is no solution to the equation (2.4).

Finally, concerning (iii), we assume that ξ (x) > 0 and ξ c (x) < 0. The monotonicity properties of ξ allow the immediate conclusion that the equation (2.4) admits two solutions. One solution stays within (0, h c ) and corresponds to a supercritical state, while the second one is in (h c , +∞) and provides a subcritical state. Moreover, if ξ c (x) = 0, the equation (2.4) admits a double root, h = h c .

The proof is thus achieved.

Remark 1. Arguing the definition of the Froude number Fr = u/c (where c is the sound speed, defined by (1.4)), we obtain for steady states Fr = q 0 / gh 3 . Therefore, it is clear that h > h c (resp. h < h c ) corresponds to Fr < 1 (resp. Fr > 1), i.e. to a subcritical (resp. supercritical) flow. Thus, the solution h of (2.4) that lies in (0, h c ) will henceforth be called the supercritical branch, while the solution of of (2.4) that lies in (h c , +∞) will be called the subcritical branch.

Remark 2. Note that the case ξ (x) = 0 implies that h = 0 is a solution of the equation (2.4). However, this equation has been obtained using the positivity of h: therefore, h = 0 is not a valid solution of (2.4).

Now, we study the solutions h of (2.4) as functions of x. As a first step, we evaluate the monotonicity of h. To this end, we exhibit the derivative of h with respect to x by considering (2.2), which writes

-q 2 0 h η-2 + gh η+1 h (x) = -kq 0 |q 0 |.
Since q 2 0 = gh 3 c from (2.6), we immediately obtain

gh η+1 1 - h 3 c h 3 h (x) = -kq 0 |q 0 |.
As a consequence, the sign of h (x) coincides with the sign of q 0 (h c -h(x)).

Note that both ξ (x) and ξ c (x) are linear in x, and ξ c (x) < ξ (x). Both of these quantities are thus linear functions of x, of slope kq 0 |q 0 |, and are strictly monotonic in x. Therefore, there exists a unique x u such that ξ (x) = 0, and a unique x c = x u such that ξ c (x) = 0. Thus, by Lemma 1, if there exists a solution h(x) of the equation (2.4), then it is either within (0, h c ] or within [h c , +∞), for all x where such a solution exists. Therefore, the sign of h (x) is constant, and h is strictly monotonic. As a consequence, h sub and h sup are bijective on their respective domains, since they are continuous and strictly monotonic. These solutions are displayed on Figure 1. For the sake of simplicity, this figure has been obtained assuming that h 0 = h c (and thus x 0 = x c ) and that q 0 < 0.

Figure 1: Solutions of (2.4) with q 0 = -1 and x 0 = xc = 0, obtained by using Newton's method. The solid line is the subcritical branch (the increasing solution) and the dotted line is the supercritical branch (the decreasing solution).

As a consequence, we have obtained the general form of steady state solutions for the shallow-water equations (1.1) with friction and flat topography, respectively given by the subcritical and supercritical branches. However, the subcritical and supercritical branches that are solution to (2.4) are parametrized by the choice of the initial conditions x 0 and h 0 .

We quickly discuss, in the context of Figure 1 (i.e. q 0 < 0 and x 0 = x c = 0) without loss of generality, the possibility of extending the smooth steady solution for x < x c . Assume that the critical point has been reached, i.e. x = x c = 0, q(x, t) = q 0 and h(x, t) = h c . Since q 0 < 0, we have two distinct steady solutions, the subcritical solution and the supercritical solution (where it exists), for x ≥ x c . However, if we were to compute a smooth steady solution for x < x c , we would still need to have q(x, t) = q 0 , as prescribed by the first steady relation. Lemma 1 states that such a solution to the second steady relation does not exist for x < x c . Therefore, there is no way to extend the smooth steady solution for x < x c .

A well-balanced scheme for a the shallow-water equations with friction

Equipped with the steady states, we now derive a numerical scheme able to exactly capture such essential solutions. We begin by recalling the construction of a Godunov-type scheme that uses a two-state approximate Riemann solver [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] (see also [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF]). The Godunov-type scheme is then adapted to the shallowwater equations with the friction source term. The purpose of this extension is to preserve the steady states presented in the previous section.

Generalities on Godunov-type schemes

First, we introduce the discretization of the space and time domain R × R + . Let ∆x be the space step, assumed to be constant, and ∆t the time step. The space discretization consists in cells (x i-1 2 , x i+ 1 2 ), for all i ∈ Z. The solution W (x, t) of (1.2) is approximated by W n i at some time t n and in cell (x i-1 2 , x i+ 1 2 ). We suppose that this piecewise constant approximation of W (x, t), here denoted W ∆ (x, t n ), is known at time t n for all cells in the space domain. We now evolve W ∆ (x, t n ) in time, which consists in nothing but the solution of the following Riemann problem juxtaposition, defined for all i ∈ Z and all x between two consecutive cells:

     ∂ t W + ∂ x f (W ) = s(W ), W (x, t n ) = W n i if x < x i+ 1 2 , W n i+1 if x > x i+ 1 2 .
(3.1)

From now on, let us emphasize that the consecutive Riemann solutions do not interact as long as t is small enough. Since an exact solution to (3.1) is difficult to determine, we suggest to consider an approximate Riemann solver to obtain an estimate of the unknown solution. We adopt a two-state approximate Riemann solver, illustrated by Figure 2 and defined as follows:

W x t ; W L , W R =          W L if x/t < λ L , W * L if λ L < x/t < 0, W * R if 0 < x/t < λ R , W R if x/t > λ R , (3.2)
where λ L and λ R denote some characteristic wave velocities, and W * L and W * R are the intermediate states, to be detailed later. To ensure that λ L < 0 and λ R > 0, we choose the following expressions of λ L and λ R (see for instance [START_REF] Toro | Riemann solvers and numerical methods dynamics. A practical introduction[END_REF] and references therein):

λ L = min (-|u L | -c L , -|u R | -c R , -ε λ ) , λ R = max (|u L | + c L , |u R | + c R , ε λ ) , (3.3)
with c the sound speed, defined by (1.4), and ε λ a positive real value to be fixed in the numerical applications. The approximate solution of the Riemann problem W ∆ (x, t n + t), depicted on Figure 3, is given for t > 0 and x ∈ (x i-1 2 , x i+ 1 2 ) by:

W L W R λ L λ R W * L W * R x 0 t
∀i ∈ Z, W ∆ (x, t n + t) =          W R, * i-1 2 if x ∈ [x i-1 2 , x i-1 2 + λ R i-1 2 t], W n i if x ∈ [x i-1 2 + λ R i-1 2 t, x i+ 1 2 + λ L i+ 1 2 t], W L, * i+ 1 2 if x ∈ [x i+ 1 2 + λ L i+ 1 2 t, x i+ 1 2 ]. (3.4)
Now, we consider a time step ∆t satisfying the following Courant-Friedrichs-Lewy stability condition :

∆t ≤ ∆x 2Λ , where Λ = max i∈Z -λ L i+ 1 , λ R i+ 1 2 , (3.5)
so as to ensure there are no interactions between the waves from two consecutive Riemann problems (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Toro | Riemann solvers and numerical methods dynamics. A practical introduction[END_REF]). Given such a ∆t, we set t n+1 = t n + ∆t and

x t t n+1 t n x i x i-1 2 x i+ 1 2 W n i W R, * i-1 2 W L, * i+ 1 2 λ R i-1 2 λ L i+ 1 2 W ∆ (x, t n+1 )
Figure 3: The full Godunov-type scheme using an approximate Riemann solver.

focus on the value of W ∆ at time t n+1 . At last, we define W n+1 i by taking the average of W ∆ (x, t n+1 ) over the cell (x i-1 2 , x i+ 1 2 ), as follows:

W n+1 i = 1 ∆x x i+ 1 2 x i-1 2 W ∆ (x, t n+1 )dx. Since W ∆ (x, t n+1
) is given by (3.4), we have:

W n+1 i = 1 ∆x x i-1 2 +λ R i-1 2 ∆t x i-1 2 W R, * i-1 2 dx + 1 ∆x x i+ 1 2 x i+ 1 2 +λ L i+ 1 2 ∆t W L, * i+ 1 2 dx + 1 ∆x x i+ 1 2 +λ L i+ 1 2 ∆t x i-1 2 +λ R i-1 2 ∆t W n i dx.
This relation immediately yields:

W n+1 i = W n i - ∆t ∆x λ L i+ 1 2 W L, * i+ 1 2 -W n i -λ R i-1 2 W R, * i-1 2 -W n i . (3.6)
Note that, in order for W n+1 i to be fully defined, we only need to give explicit values to the intermediate states

W L, * i+ 1 2 and W R, * i+ 1 2
. The remainder of this section is devoted to finding suitable intermediate states, which ensure that the scheme is consistent and preserves the steady states.

A well-balanced approximate Riemann solver for the friction source term

We focus on the derivation of a well-balanced approximate Riemann solver for the shallow-water equations with only the friction source term. First, we recall the system under consideration:

   ∂ t h + ∂ x q = 0, ∂ t q + ∂ x q 2 h + 1 2 gh 2 = -kq|q|h -η . (3.7)
We also recall that the steady states for this model are given by:

     q = q 0 , ∂ x q 2 0 h + 1 2 gh 2 = -kq|q|h -η . (3.8)
Note that the friction source term S f = -kq|q|h -η in the system (3.7) falls under the more generic framework S f = h β ρ(q)∂ x Y , where β = -η, ρ(q) = q|q| and Y = -kx. The computations could be performed in this framework, but they would lead to having to solve a nonlinear equation. However, thanks to the specific form of the friction source term, the scheme can be obtained without having to solve this nonlinear equation. As a consequence, we elect to derive the scheme for the specific case of the friction source term instead of a more generic form. These remarks also apply to the topography source term S t = -gh∂ x Z, which is treated in a later section.

Relevant definitions of the intermediate states

W * L = t (h * L , q * L ) and W * R = t (h * R , q * R )
need to be obtained in order to fully determine the approximate Riemann solver. Therefore, four equations characterizing the four unknowns h * L , h * R , q * L and q * R are required. These four equations will be obtained by considering two essential properties that need to be satisfied by the approximate Riemann solver: consistency with the system (3.7) and well-balance, i.e. preservation of the steady states (3.8).

First, we tackle the issue of the consistency with (3.7), by introducing a necessary consistency condition. From [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF], the average over a cell of the approximate Riemann solver W , defined by (3.2), has to be equal to the average over the same cell of the exact solution of the Riemann problem W R . Therefore, the following equality is imposed:

1 ∆x ∆x/2 -∆x/2 W x ∆t ; W L , W R dx = 1 ∆x ∆x/2 -∆x/2 W R x ∆t ; W L , W R dx. (3.9)
We introduce the following notations for the sake of simplicity:

(λ R -λ L )h HLL = λ R h R -λ L h L -[q], (λ R -λ L )q HLL = λ R q R -λ L q L - q 2 h + 1 2 gh 2 .
Let us underline that t (h HLL , q HLL ) is actually the well-known intermediate state of the HLL approximate Riemann solver introduced in [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]. Also, note that h HLL > 0 for λ L and λ R given by (3.3). We also introduce the notation (-kq|q|h -η ) R that denotes the second component of s(W R ), with s defined by (1.3), to represent the value of the source term of friction for the exact solution of the Riemann problem. Thanks to these notations, (3.9) can be rewritten as follows, after straightforward computations detailed in [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF]:

λ R h * R -λ L h * L = (λ R -λ L )h HLL , (3.10a) λ R q * R -λ L q * L = (λ R -λ L )q HLL + 1 ∆t ∆x/2 -∆x/2 ∆t 0 -kq|q|h -η R x t ; W L , W R dt dx. (3.10b)
Now, we introduce two parameters, q and h -η , to define a consistent approximation of the mean value of the friction source term, as follows:

1 ∆x 1 ∆t ∆x/2 -∆x/2 ∆t 0 -kq|q|h -η R dt dx -kq|q|h -η . (3.11)
The parameters q and h -η are respectively assumed to be consistent with q and h -η , in a sense given later. Then, the equation (3.11) is substituted into (3.10b). As a consequence, the four unknowns are imposed to satisfy the following two relations:

λ R h * R -λ L h * L = (λ R -λ L )h HLL , λ R q * R -λ L q * L = (λ R -λ L )q HLL -kq|q|h -η ∆x. (3.12)
Concerning the well-balance of the approximate Riemann solver, we exhibit a sufficient condition for the preservation of a solution. The scheme (3.6) provides us with such a condition. Indeed, the solution is obviously stationary, i.e.

W n+1 i = W n i for all i ∈ Z, if we have ∀i ∈ Z, W L, * i+ 1 2 = W n i and W R, * i-1 2 = W n i .
Thus, in the context of the approximate Riemann solver, the solution will be stationary if W * L = W L and W * R = W R . As a consequence, the following wellbalance principle will have to be satisfied by the intermediate states:

Principle (WB). The intermediate states W * L and W * R are such that W * L = W L and W * R = W R as soon as W L and W R define a steady state.

Here, the pair (W L , W R ) is said to define a steady state if they satisfy the identity (3.8). We will then look for intermediate states W * L and W * R that satisfy the (WB) principle.

We begin by determining q * L and q * R . After [START_REF] Berthon | A fully well-balanced, positive and entropysatisfying Godunov-type method for the shallow-water equations[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF], we choose q * L = q * R and we set q * := q * L = q * R . This choice helps guarantee that W * L and W * R satisfy the property (WB). The formulas (3.12) are then rewritten using q * as follows:

λ R h * R -λ L h * L = (λ R -λ L )h HLL , (3.13a
)

q * = q HLL - kq|q|h -η ∆x λ R -λ L . (3.13b)
We only need relevant definitions of q and h -η in order to fully determine q * . To that end, we specify the discrete steady states satisfied by the system (3.7). Such discrete steady states are nothing but a discretization of the steady relation at the continuous level (3.8). In order to give an expression of these discrete steady states, we assume that h L = 0 and h R = 0. Now, W L and W R define a steady state if the following relations hold:

   q L = q R = q 0 , (3.14a)
q 2 0 1 h + g 2 h 2 = -kq|q|h -η ∆x. (3.14b)
Moreover, since we are considering smooth steady states for just the friction source term, we can rewrite the relation (2.3) between states W L and W R . Thus, W L and W R satisfy the following algebraic relation:

- q 2 0 η -1 h η-1 + g η + 2 h η+2 = -kq 0 |q 0 | ∆x. (3.15)
We introduce the notation µ 0 = sgn(q 0 ). Therefore, µ 0 represents the direction of the water flow. From (3.15), we obtain the evaluation of q 2 0 as follows:

q 2 0 = g [h η+2 ] η + 2 [h η-1 ] η -1 -kµ 0 ∆x . (3.16)
Now, to recover the behavior of the friction source term when steady states are involved, we assume that the parameter q is equal to q 0 as soon as W L and W R define a steady state. Thus, injecting (3.16) into (3.14b) yields the expression of h -η . We immediately obtain:

h -η = [h 2 ] 2 η + 2 [h η+2 ] - µ 0 k ∆x 1 h + [h 2 ] 2 [h η-1 ] η -1 η + 2 [h η+2 ] . (3.17)
Lemma 2. The expression of h -η given by (3.17) is consistent with h -η .

Proof. With smooth water heights, we fix, in (3.17), h L = h(x) and h R = h(x+O(∆x)). Taylor's formula applied to h R yields h R = h+∆x∂ x h+O(∆x 2 ).

In order to evaluate the Taylor expansions of

[h 2 ], [h η-1 ], [h η+2 ] and [h -1 ],
we now compute a Taylor expansion, for some β ∈ R, of the jump [h β ]:

[h β ] = h β R -h β L = h + ∂ x h∆x + O(∆x 2 ) β -h β = h β 1 + βh -1 ∂ x h∆x + O(∆x 2 ) -h β = βh β-1 ∂ x h∆x + O(∆x 2 ).
Using the above evaluation, we have for the first part of the expression of h -η :

[h 2 ] 2 η + 2 [h η+2 ] = h∂ x h∆x + O(∆x 2 ) h η+1 ∂ x h∆x + O(∆x 2 ) = h -η + O(∆x). (3.18)
Moreover, we have, for the second part of h -η ,

1 h = -h -2 ∂ x h∆x + O(∆x 2 ), (3.19)
[h 2 ] 2 [h η-1 ] η -1 η + 2 [h η+2 ] = (h∂ x h∆x + O(∆x 2 ))(h η-2 ∂ x h∆x + O(∆x 2 )) h η+1 ∂ x h∆x + O(∆x 2 ) = h -2 ∂ x h∆x + O(∆x 2 ). (3.20)
Combining both equations (3.19) and (3.20) immediately yields

- µ 0 k ∆x 1 h + [h 2 ] 2 [h η-1 ] η -1 η + 2 [h η+2 ] = - µ 0 k ∆x O(∆x 2 ) = O(∆x). (3.21)
Using both relations (3.18) and (3.21) gives h -η = h -η + O(∆x), which concludes the proof.

Concerning q, we choose the following average:

   q = 2|q L ||q R | |q L | + |q R | sgn(q L + q R ) if q L = 0 and q R = 0; q = 0 if q L = 0, q R = 0 or k = 0. (3.22)
This average indeed ensures that, if q L = q R , then q = q L = q R . Note that the expression (3.17) of h -η contains µ 0 . This quantity depends on the steady state and has to be determined for non-steady states. To address such an issue, we propose the expression

h -η := [h 2 ] 2 η + 2 [h η+2 ] - µ k ∆x 1 h + [h 2 ] 2 [h η-1 ] η -1 η + 2 [h η+2 ] , (3.23)
where µ is the sign of the quantity q given by (3.22). Since q = q 0 as soon as W L and W R define a steady state, it is clear that h -η , as defined by (3.23), still provides a well-balanced approximate Riemann solver and still satisfies Lemma 2.

Remark 3. With h -η given by (3.23), the average source term -kq|q|h -η , which satisfies (3.11), rewrites as follows:

-

kq|q|h -η = -kq|q| [h 2 ] 2 η + 2 [h η+2 ] + (q) 2 ∆x 1 h + [h 2 ] 2 [h η-1 ] η -1 η + 2 [h η+2 ] . (3.24)
As a consequence, when the Manning coefficient k is almost zero, the second term of this average source term is dominant. From the proof of Lemma 2, we know that this second term is consistent with 0 for a smooth solution. However, this term could become inconsistent within shock waves. To address such an issue, we introduce a cutoff of this inconsistent term. Such a procedure had already been adopted in [START_REF] Berthon | A fully well-balanced, positive and entropysatisfying Godunov-type method for the shallow-water equations[END_REF][START_REF] Berthon | Fully well-balanced, positive and simple approximate Riemann solver for shallow water equations[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF] in the context of the sole topography source term. Therefore, we modify the expression (3.23) of h -η , as follows:

h -η := h -η (h L , h R ) = [h 2 ] 2 η + 2 [h η+2 ] - µ k ∆x [h] c -1 h L h R + h L + h R 2 [h η-1 ] η -1 η + 2 [h η+2 ] , (3.25) where [h] c is a cutoff of [h] = h R -h L , defined as follows: [h] c = h R -h L if |h R -h L | ≤ C∆x, sgn(h R -h L )C∆x otherwise, (3.26)
with C a positive constant that does not depend on ∆x. Equipped with this cutoff (3.26) within the expression (3.25) of h -η , it is clear from the proof of Lemma 2 that the approximate friction source term given by (3.24) becomes consistent, i.e. of order O(∆x), when k tends to 0. Finally, the approximate friction source term reads:

-kq|q|h -η = -kq|q| [h 2 ] 2 η + 2 [h η+2 ] + (q) 2 ∆x [h] c -1 h L h R + h L + h R 2 [h η-1 ] η -1 η + 2 [h η+2 ] .
(3.27) According to Lemma 2, the expression (3.27) is consistent with the friction source term -kq|q|h -η for all k ≥ 0.

Equipped with the expressions of h -η and q, we have fully determined q * . Now, we need to complete the system (3.13) to determine the intermediate water heights h * L and h * R . The following relation is chosen so as to ensure that W * L and W * R satisfy the property (WB):

α f (h * R -h * L ) = -kq|q|h -η ∆x, where α f = -q 2 h L h R + g 2 (h L + h R ). (3.28)
Therefore, h * L and h * R are easily determined by considering the following linear system, defined by (3.13a) and (3.28):

λ R h * R -λ L h * L = (λ R -λ L )h HLL , α f (h * R -h * L ) = -kq|q|h -η ∆x.
(3.29)

The system (3.29) then provides the following formulas:

h * L = h HLL + λ R kq|q|h -η ∆x α f (λ R -λ L ) , (3.30a)
h * R = h HLL + λ L kq|q|h -η ∆x α f (λ R -λ L ) . (3.30b)
Remark 4. Note that, in (3.30), a division by α f occurs. In rare transcritical cases, the term α f could become close to zero, which could lead to a significant error in the numerical solution. However, many simulations have been performed and this problem has not appeared.

Properties satisfied by the approximate Riemann solver

In the previous subsection, we have constructed intermediate states (3.13b) -(3.30) that yield a consistent and well-balanced approximate Riemann solver. However, one can plainly see from (3.30) that the intermediate water heights h * L and h * R can become negative. If that is the case, (3.6) can yield a negative updated water height, which means the updated state would not belong to the admissible space Ω. To recover a physically admissible solution, we apply the procedure from [START_REF] Audusse | A simple well-balanced and positive numerical scheme for the shallow-water system[END_REF][START_REF] Berthon | A Well-Balanced Finite Volume Scheme for a Mixed Hyperbolic/Parabolic System to Model Chemotaxis[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF], which consists in introducing a parameter ε > 0, that controls the positivity of h * L and h * R . We now state the expressions of the intermediate states obtained with this procedure, for given W L and W R . The parameter ε is such that

ε = min (h L , h R , h HLL ) . (3.31)
For the sake of simplicity, we introduce the following notation:

S f := -kq|q|h -η , (3.32)
where q and h -η are respectively defined by (3.22) and (3.25). Then, the intermediate states

W * L = t (h * L , q * L ) and W * R = t (h * R , q * R )
of the approximate Riemann solver (3.2) are given by:

α f = -q 2 h L h R + g 2 (h L + h R ), (3.33a)
q * L = q * R = q * = q HLL + S f ∆x λ R -λ L , (3.33b)
h * L = min max h HLL - λ R S f ∆x α f (λ R -λ L ) , ε , 1 - λ R λ L h HLL + λ R λ L ε , (3.33c)
h * R = min max h HLL - λ L S f ∆x α f (λ R -λ L ) , ε , 1 - λ L λ R h HLL + λ L λ R ε .
(3.33d)

The quantities q * and α f have respectively been defined by (3.13b) and (3.28). Moreover, in (3.33), the positivity correction has been applied, and new expressions (3.33c) and (3.33d) of h * L and h * R have been substituted to (3.30a) and (3.30b). Note that the intermediate states (3.33) are written under the same form as the intermediate states presented in [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF].

The following two results sum up the properties we have obtained thus far, for the approximate Riemann solver and the full scheme. For the sake of conciseness, we do not include their proofs here, since they use classical ingredients that have been presented in [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF].

Lemma 3. Assume ε > 0 such that (3.31) is satisfied. Then, intermediate states W * L = t (h * L , q * L ) and W * R = t (h * R , q * R )
given by (3.33) satisfy the following properties:

(i) consistency: the quantities h * L , h * R , q * L and q * R satisfy the equations (3.13);

(ii) positivity preservation: if

h L > 0, h R > 0 and h HLL > 0, then h * L ≥ ε and h * R ≥ ε; (iii) well-balance: W * L = W L and W * R = W R satisfy the property (WB). Theorem 4. Consider W n i ∈ Ω * for all i ∈ Z,
where Ω * is a restricted admissible states space defined as follows:

Ω * = W = t (h, q) ∈ R 2 ; h > 0, q ∈ R . Assume that the intermediate states W L, * i+ 1 2 and W R, * i+ 1 2 are given, for all i ∈ Z, by W L, * i+ 1 2 = h * L (W n i , W n i+1 ) q * (W n i , W n i+1 ) and W R, * i+ 1 2 = h * R (W n i , W n i+1 ) q * (W n i , W n i+1 )
, where q * , h * L and h * R are given by (3.33b), (3.33c) and (3.33d), respectively. Also, assume ε > 0 given by (3.31). Then the Godunov-type scheme, given by (3.6) under the CFL restriction (3.5), satisfies the following properties:

1. consistency with the shallow-water system (3.7);

positivity preservation

: ∀i ∈ Z, W n+1 i ∈ Ω * ; 3. well-balance: if (W n i ) i∈Z defines a steady state, then ∀i ∈ Z, W n+1 i = W n i .
From the above result, we see that the updated water height never vanishes. This behavior is due to the introduction in (3.33) of the parameter ε. Such a behavior is necessary because the average h -η has been obtained by considering positive water heights. In a later section, devoted to dry/wet transitions, this average will be extended to be able to deal with dry areas. This extension will allow us to consider ε = 0 in the definition (3.33), as well as vanishing water heights.

The case of both topography and friction source terms

In this section, our goal is the derivation of new intermediate states h * L , h * R and q * to approximate the contributions of both topography and friction source terms. For the sake of simplicity in the notations, we write S = S t + S f , where

S t = -gh∂ x Z, S f = -kq|q|h -η ,
which corresponds to the case where both topography and friction source terms are involved. The steady states associated to the full source term S do not admit an algebraic expression. Therefore S is split between S t and S f , since the steady states for the individual source terms of topography and friction can both be rewritten using an algebraic expression. In order to complete the scheme, we need to determine suitable intermediate states in both cases of the topography and the friction. The case of the friction has been treated in the previous section, and the intermediate states are given by (3.33). For the topography, we adopt the approach introduced in [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF], which is now recalled for the sake of completeness.

Approximation of the topography source term

Considering a non-vanishing topography source term but with a vanishing friction, the intermediate states are given by

α t = -(q * ) 2 h L h R + g 2 (h L + h R ), (4.1a)
q * L = q * R = q * = q HLL + S t ∆x λ R -λ L , (4.1b)
h * L = min max h HLL - λ R S t ∆x α t (λ R -λ L ) , ε , 1 - λ R λ L h HLL + λ R λ L ε , (4.1c) h * R = min max h HLL - λ L S t ∆x α t (λ R -λ L ) , ε , 1 - λ L λ R h HLL + λ L λ R ε . (4.1d)
In (4.1), the approximate topography source term S t is defined by:

S t ∆x := S t (h L , h R , Z L , Z R )∆x = -2g[Z] h L h R h L + h R + g 2 [h] 3 c h L + h R , (4.2)
where (3.26). This expression of S t is consistent with the topography source term S t , and it ensures that the scheme is well-balanced (see [START_REF] Berthon | A fully well-balanced, positive and entropysatisfying Godunov-type method for the shallow-water equations[END_REF][START_REF] Berthon | Fully well-balanced, positive and simple approximate Riemann solver for shallow water equations[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF] for more details regarding this expression and how to obtain it). Computations leading to this expression and to the expressions (4.1) have been explained in [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF], to which the reader is referred for more details.

[h] c is the cutoff of [h] = h R -h L defined by

Approximation of the topography and friction source terms

Now, we focus on deriving intermediate states for the shallow-water equations with both topography and friction source terms. We introduce the following discretization of (1.7):

q 2 0 1 h + g 2 h 2 = S t ∆x + S f ∆x, (4.3)
where the approximate source term S has been split into the topography and friction contributions, by writing S = S t + S f . As a consequence, after (3.33b), we define the following intermediate discharge:

q * = q HLL + S t ∆x λ R -λ L + S f ∆x λ R -λ L , (4.4)
where the expressions of S t and S f are given by (4.2) and (3.32), respectively. Concerning the intermediate heights, we first decide to use the following expressions:

h * L = h HLL - λ R S t ∆x α t (λ R -λ L ) - λ R S f ∆x α f (λ R -λ L ) , h * R = h HLL - λ L S t ∆x α t (λ R -λ L ) - λ L S f ∆x α f (λ R -λ L ) , (4.5)
where the quantities α t and α f are respectively defined by (4.1a) and (3.33a).

We immediately see that the intermediate states (4.4) -(4.5) still yield a consistent scheme. Moreover, it is possible to recover the positivity of the water heights. Let ε ≥ 0 that satisfies (3.31), and consider the same procedure as used for the friction source term, which yields the following intermediate states:

q * = q HLL + S t ∆x λ R -λ L + S f ∆x λ R -λ L , h * L = min max h * L , ε , 1 - λ R λ L h HLL + λ R λ L ε , h * R = min max h * R , ε , 1 - λ L λ R h HLL + λ L λ R ε . (4.6)
From Lemma 3, as soon as ε > 0, the intermediate heights h * L and h * R defined by (4.6) are positive.

We note that using the definitions (4.6) and making the friction source term vanish allows the recovery of the intermediate states for topography only. Similarly, if the topography source term vanishes, we recover the intermediate states for friction only. As a consequence, (4.6) yields intermediate states that are well-balanced for the individual source terms of topography or friction.

Let us recall that the steady states relation for the shallow-water system with both topography and friction source terms (1.5) cannot be written under the form of an algebraic relation for all Z. Therefore, we only manage to preserve the steady states up to the chosen discretization (4.3) of (1.7) (see [START_REF] Desveaux | well-balanced scheme for the Euler equation with a gravitational potential[END_REF][START_REF] Käppeli | Well-balanced schemes for the Euler equations with gravitation[END_REF][START_REF] Luo | A well-balanced symplecticity-preserving gaskinetic scheme for hydrodynamic equations under gravitational field[END_REF][START_REF] Xu | A well-balanced gas-kinetic scheme for the shallow-water equations with source terms[END_REF][START_REF] Xu | A well-balanced kinetic scheme for gas dynamic equations under gravitational field[END_REF] where a similar approach is used).

We finally study how the approximate source term averages S t and S f , as well as the terms S t ∆x/α t and S f ∆x/α f , behave when dealing with vanishing water heights. First, we make the following assumption.

Assumption. When the height vanishes, so does the velocity.

From [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF], under this assumption, the quantity S t ∆x satisfies:

S t ∆x := S t (h L , h R , Z L , Z R )∆x =                0 if h L = 0 and h R = 0, -g[Z] h R + h L 2 if h L = 0 or h R = 0, -g[Z] 2h L h R h L + h R + g 2 [h] 3 c h L + h R otherwise, (4.7)
and the quantity S t ∆x/α t satisfies:

S t ∆x α t =                  0 if h L = 0 and h R = 0, -[Z] if h L = 0 or h R = 0, S t ∆x -(q * ) 2 h L h R + g 2 (h L + h R )
otherwise.

We then turn to the approximate friction source term. We have to make the following assumption in order to define it for vanishing water heights.

Assumption. Using S f and α f respectively defined by (3.32) and (3.33a), we impose that S f and S f /α f are zero as soon as h L and/or h R vanishes.

It makes sense to consider such a behavior. Indeed, as previously mentioned, the friction contribution is assumed to vanish as soon as the water height vanishes. In order for both quantities S f and S f /α f to satisfy this requirement, we have to impose that they vanish when h L and/or h R vanishes.

To handle the case where both h L and h R vanish (and thus q L = q R = 0), we have to make sure that, in this case, q * = 0 and h * L = h * R = 0. This requirement is met by taking S t = 0 and S t /α t = 0, as well as S f = 0 and S f /α f = 0, as soon as both h L and h R are zero.

As a consequence, after (3.32) and the above assumption, S f ∆x is given by:

S f ∆x := S f (h L , h R , q L , q R )∆x = 0 if h L = 0 and/or h R = 0, -kq|q|h -η otherwise, (4.8)
where q is given by (3.22) and h -η is defined by (3.25). In addition, the quantity S f ∆x/α is given by:

S f ∆x α =          0 if h L = 0 and/or h R = 0, -kq|q|h -η -(q * ) 2 h L h R + g 2 (h L + h R )
otherwise.

We can now state the following result, concerning the approximate Riemann solver for both source terms of topography and friction. Since its proof uses classical ingredients, it is not presented here; instead, the reader is referred to [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF]. (ii) well-balance: if W L and W R define a steady state, i.e. satisfy (4.3), then

W * L = W L and W * R = W R . Moreover, we have: (iii) with ε > 0, the positivity is preserved: if h L > 0, h R > 0 and h HLL > 0, then h * L ≥ ε and h * R ≥ ε; (iv) with ε = 0, the non-negativity is preserved: if h L ≥ 0, h R ≥ 0 and h HLL ≥ 0, then h * L ≥ 0 and h * R ≥ 0.
This lemma allows us to state the following result, which concerns the full scheme (3.6). Similarly to Lemma 5, we do not prove this result here, and we refer the reader to [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF].

Theorem 6. Consider W n i ∈ Ω for all i ∈ Z. Assume that the intermediate states W L, * i+ 1 2 and W R, * i+ 1 2 are given, for all i ∈ Z, by W L, * i+ 1 2 = h * L (W n i , W n i+1 ) q * (W n i , W n i+1 ) and W R, * i+ 1 2 = h * R (W n i , W n i+1 ) q * (W n i , W n i+1 ) , (4.9)
where q * , h * L and h * R are given by (4.6) with ε = 0, i.e. by

q * = q HLL + S t ∆x λ R -λ L + S f ∆x λ R -λ L , h * L = min h HLL - λ R S t ∆x α t (λ R -λ L ) - λ R S f ∆x α f (λ R -λ L ) + , 1 - λ R λ L h HLL , h * R = min h HLL - λ L S t ∆x α t (λ R -λ L ) - λ L S f ∆x α f (λ R -λ L ) + , 1 - λ L λ R h HLL .
Then the Godunov-type scheme, given by (3.6) under the CFL restriction (3.5), satisfies the following properties:

1. consistency with the shallow-water system (1.1);

non-negativity preservation

: ∀i ∈ Z, W n+1 i ∈ Ω; 3. well-balance: if (W n i ) i∈Z defines a steady state according to (4.3), then ∀i ∈ Z, W n+1 i = W n i .
The same result holds for both topography-only and friction-only steady states, respectively given by (1.9) and (1.10).

Implicitation of the source terms contribution: wet/dry transitions

The scheme (3.6) -(4.9) allows the simulation of wet/dry transitions. However, several spurious oscillations appear when such experiments are performed. The stiffness of the source terms in the vicinity of wet/dry transitions, especially the friction, induces these spurious oscillations. To correct such a failure, we first rewrite this scheme in order to exhibit the numerical flux function and the source terms contribution. Then, we adopt an explicit scheme for the transport part, and an implicit scheme for the source part.

Reformulation of the scheme

In this subsection, we exhibit the numerical flux function and the numerical source terms. To that end, the scheme (3.6), with intermediate states given by (4.6), is rewritten after straightforward computations as follows (see for instance [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]):

W n+1 i = W n i - ∆t ∆x f n i+ 1 2 -f n i-1 2 + ∆t 2 s n i+ 1 2 + s n i-1 2 .
(5.1)

The quantity

f n i+ 1 2 = f (W n i , W n i+1
) is the numerical flux function evaluated at the interface x i+ 1 2 , and the quantity s n i+ 1 2 is the numerical source term at the interface x i+ 1 2 . These quantities are approximations of the flux and the source term, respectively, and are defined by:

f n i+ 1 2 = (f h ) n i+ 1 2 (f q ) n i+ 1 2 and s n i+ 1 2 = 0 (S t ) n i+ 1 2 + (S f ) n i+ 1 2 .
(5.

2)

The quantities

(S t ) n i+ 1 2 and (S f ) n i+ 1 2
are approximations of the topography and the friction source terms, respectively. Adopting extended notations, they are given by:

(S t ) n i+ 1 2 = S t h n i , h n i+1 , Z i , Z i+1 , (5.3a)
(S f ) n i+ 1 2 = S f h n i , h n i+1 , q n i , q n i+1 , (5.3b)
where S t and S f are the approximate source terms already defined by (4.7) and (4.8). The scheme (5.1) then reads:

h n+1 i = h n i - ∆t ∆x (f h ) n i+ 1 2 -(f h ) n i-1 2 ,
(5.4a)

q n+1 i = q n i - ∆t ∆x (f q ) n i+ 1 2 -(f q ) n i-1 2 + ∆t 2 (S t ) n i+ 1 2 + (S t ) n i-1 2 + (S f ) n i+ 1 2 + (S f ) n i-1 2 , (5.4b)
where the approximate fluxes are defined by (5.2) and

f n i+ 1 2 = 1 2 F (W n i ) + F (W n i+1 ) + λ L i+ 1 2 2 W L, * i+ 1 2 -W n i + λ R i+ 1 2 2 W R, * i+ 1 2 -W n i+1 .
(5.5)

For the sake of simplicity, we set the notations:

(S t ) n i = 1 2 (S t ) n i-1 2 + (S t ) n i+ 1 2 and (S f ) n i = 1 2 (S f ) n i-1 2 + (S f ) n i+ 1 2 ,
(5.6) such that (5.4b) rewrites:

q n+1 i = q n i - ∆t ∆x (f q ) n i+ 1 2 -(f q ) n i-1 2 + ∆t(S t ) n i + ∆t(S f ) n i .
(5.7)

Implicitation of the source terms contribution

We now introduce a semi-implicit version of the Godunov-type scheme (5.1). The main idea of this section is to use a splitting method (see for instance [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Toro | Riemann solvers and numerical methods dynamics. A practical introduction[END_REF]) to reduce the impact of the aforementioned instabilities. The splitting strategy we use here is to first consider an explicit treatment of the flux, then an implicit treatment of both source terms. As a consequence, the first step, devoted to the transport part ∂ t W + ∂ x f (W ) = 0, reads as follows:

h n+ 1 3 i = h n i - ∆t ∆x (f h ) n i+ 1 2 -(f h ) n i-1 2 , q n+ 1 3 i = q n i - ∆t ∆x (f q ) n i+ 1 2 -(f q ) n i-1 2 .
(5.8)

During the second step, devoted to the topography source term, we approximately solve the initial value problem

     dh dt = 0, dq dt = S t (W ),
with initial data

   h(0) = h n+ 1 3 i , q(0) = q n+ 1 3 i
, where S t (W ) = -gh∂ t Z. We suggest the following implicit scheme to approximate h

n+ 2 3 i and q n+ 2 3 i : h n+ 2 3 i = h n+ 1 3 i , (5.9a)
q n+ 2 3 i = q n+ 1 3 i + ∆t(S t ) n+ 2 3 i , (5.9b)
where, as per (5.6), (S t )

n+ 2 3 i = (S t ) n+ 2 3 i+ 1 2 + (S t ) n+ 2 3 i-1 2
/2, with, after (5.3a),

(S t ) n+ 2 3 i+ 1 2 = S t h n+ 2 3 i , h n+ 2 3 i+1 , Z i , Z i+1 .
Let us underline that, during the second step, the evaluation of the implicit unknowns is immediate since the topography source term involves known quantities.

The third and last step concerns the friction and consists in solving

     dh dt = 0, dq dt = -kq|q|(h n+1 i ) -η , with initial data    h(0) = h n+ 2 3 i , q(0) = q n+ 2 3 i .
(5.10)

This system can be solved to find an analytic expression of the solution. For t ∈ [0, ∆t], the exact solution of the above system reads as follows: .11) Note that the analytic expression (5.11) guarantees that, for all t ∈ (0, ∆t], the sign of q(t) stays the same as the sign of q(0), and that |q(t)| < |q(0)|. This behavior is consistent with the fact that friction should only slow down the movement of the fluid, rather than changing its direction. Then, evaluating (5.11) at t = ∆t and plugging the initial data yields

   h(t) = h(0), q(t) = h(0) η q(0) h(0) η + k t |q(0)| . ( 5
       h n+1 i = h n+ 2 3 i , (5.12a
)

q n+1 i = h n+1 i η q n+ 2 3 i h n+1 i η + k ∆t q n+ 2 3 i . (5.12b)
Let us note that the updated state t (h n+1 i , q n+1 i ) retains its non-negativity preservation property, since the computations for h n+1 i have not been modified and Theorem 6 applies. Furthermore, the expression (5.12b) is well-posed for h n+1 i = 0 and, when that is the case, forces the discharge q n+1 i to vanish. However, the well-balance property is lost. Indeed, if W n i-1 , W n i and W n i+1 define a steady state, we do not necessarily recover q n+1 i = q n i . As a consequence, we decide to consider an approximation of h n+1 i η , denoted by (h η ) n+1 i , in (5.12b). This approximation is determined in order to ensure the required well-balance of the scheme.

In order to obtain such an expression of (h η ) n+1 i , we momentarily suppose that W n i-1 , W n i and W n i+1 define a steady state. In this case, we need to ensure that W n+1 i = W n i . Since the scheme (5.1) is well-balanced, the equation (5.7) yields q n+1 i = q n i , and can be rewritten as:

(S f ) n i = 1 ∆x (f q ) n i+ 1 2 -(f q ) n i-1 2 -(S t ) n i .
(5.13) Moreover, we have h

n+ 2 3 i = h n i .
Therefore, by the definition (5.3a) of (S t ) n i , we have (S t )

n+ 2 3 i = (S t ) n i . The evaluation of q n+ 2 3 i
is then obtained from (5.9b) as follows:

q n+ 2 3 i = q n i - ∆t ∆x (f q ) n i+ 1 2 -(f q ) n i-1 2 + ∆t(S t ) n i . (5.14)
Since W n i-1 , W n i and W n i+1 define a steady state, the relations (5.13) and (5.14) yield: q

n+ 2 3 i = q n i -∆t(S f ) n i . Thus, using (h η ) n+1 i , (5.12b) reads q n i = (h η ) n+1 i q n i -∆t(S f ) n i (h η ) n+1 i + k ∆t |q n i -∆t(S f ) n i | . ( 5.15)
We are now able to determine the expression of (h η ) n+1 i that ensures the wellbalance of the scheme. We set µ n i = sgn q n i and µ

n+ 2 3 i = sgn q n+ 2 3 i
. With these notations, we get from relation (5.15) the following expression of (h η ) n+1 i :

(h η ) n+1 i = - k(q n i ) 2 µ n+ 2 3 i (S f ) n i + k∆tµ n+ 2 3 i q n i .
(5.16)

Now, recall that the numerical source term (S f ) n i is defined by (5.6). As a consequence, we get the following expression for (S f ) n i :

(S f ) n i = 1 2 -kq n i-1 2 |q n i-1 2 |(h -η ) n i-1 2 -kq n i+ 1 2 |q n i+ 1 2 |(h -η ) n i+ 1 2 ,
where the averages

q n i± 1 2 and (h -η ) n i± 1 2
are given with clear notations by (3.22) and (3.25), respectively. However, recall that the only requirement that has to be fulfilled by the average q is to recover the steady discharge q 0 as soon as a steady state was reached. In the current context, W n i-1 , W n i and W n i+1 define a steady state; hence, we have q n i-1 = q n i = q n i+1 = q 0 . A relevant choice is therefore to take q n i± 1 2 = q n i , which yields the following formula:

(S f ) n i = 1 2 -kq n i |q n i | h -η n+1 i-1 2 -kq n i |q n i | h -η n+1 i+ 1 2 .
(5.17)

In (5.17), we have substituted

(h -η ) n i± 1 2 with (h -η ) n+1 i± 1 2
. This substitution has no effect on the well-balance property, and considering the updated water height makes the scheme contain more implicit terms.

With this simplification in the source term approximation, we get the following expression for (h η ) n+1 i , from (5.16) and (5.17):

(h η ) n+1 i = 2µ n+ 2 3 i µ n i h -η n+1 i-1 2 + h -η n+1 i+ 1 2 + k ∆t µ n+ 2 3 i q n i .
Arguing the formula (3.25), the above expression can then be rewritten as

(h η ) n+1 i = 2kµ n+ 2 3 i ∆x kµ n i ∆x β n+1 i-1 2 + β n+1 i+ 1 2 -γ n+1 i-1 2 + γ n+1 i+ 1 2 + k ∆t µ n+ 2 3 i q n i , (5.18)
where we have set

β n+1 i+ 1 2 = η + 2 2 h n+1 i+1 2 -h n+1 i 2 h n+1 i+1 η+2 -h n+1 i η+2
, and

γ n+1 i+ 1 2 = 1 h n+1 i+1 - 1 h n+1 i + β n+1 i+ 1 2 h n+1 i+1 η-1 -h n+1 i η-1 η -1 .
Computations within the expression of (h η ) n+1 i show that it tends to 0 as soon as h n+1 i-1 , h n+1 i or h n+1 i+1 tends to 0, which is the correct behavior when dealing with wet/dry transitions. We have therefore devised a way to consider the source terms contribution in an implicit way, while still retaining the well-balance property of the scheme. We can thus state the following result. Theorem 7. The scheme (5.8) -(5.9) -(5.12) -(5.18) is positivity-preserving and well-balanced. Indeed, if (W n i ) i∈Z defines a steady state according to (4.3), then ∀i ∈ Z, W n+1 i = W n i . The same result holds for both topography-only and friction-only steady states, respectively given by (1.9) and (1.10).

Proof. We begin this proof by recalling that, from (5.8), (5.9) and (5.12), we have

h n+1 i = h n+ 1 3 i . Therefore, since h n+ 1 3 i
is given by the full scheme (3.6), which is positivity-preserving, the scheme under consideration is positivitypreserving. We now have to show that it is well-balanced. To that end, assume that W n i-1 , W n i and W n i+1 define a steady state, according to (4.3), (1.9) or (1.10). Using a similar chain of arguments as above, we immediately recover that h n+1 i = h n i . To complete the proof, we now have to show that q n+1 i = q n i . The updated discharge q n+1 i is given by (5.12), with q n+ 2 3 i defined by (5.9) and (5.8). Since (h η ) n+1 i is given by (5.18) and has been chosen to ensure q n+1 i = q n i , the proof is concluded.

Second-order MUSCL extension

We devote this section to a second-order extension based on a MUSCL technique (for instance, see [START_REF] Van Leer | Towards the Ultimate Conservative Difference Scheme, V. A Second Order Sequel to Godunov's Method[END_REF][START_REF] Van Leer | On the relation between the upwind-differencing schemes of Godunov, Engquist-Osher and Roe[END_REF][START_REF] Leveque | Numerical methods for conservation laws[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF][START_REF] Toro | Riemann solvers and numerical methods dynamics. A practical introduction[END_REF]), to improve the space accuracy of the scheme. The MUSCL procedure involves a piecewise linear reconstruction, instead of piecewise constant, in the Godunov-type scheme. In addition, we use Heun's method to increase the time accuracy of the scheme.

The MUSCL reconstruction

The variables to be reconstructed are h, q and h + Z. Thus, note that the approximation of Z turns out to be time-dependent. For the sake of simplicity and conciseness in the notations, consider w ∈ {h, q, h + Z}. The reconstruction procedure consists in replacing the constant state w n i with a linear approximation, given in each cell

(x i-1 2 , x i+ 1 2 ) by w n i (x) = w n i + (x -x i ) σ n i ,
where σ n i is the limited slope of the linear reconstruction. A limiter is applied in order to improve the stability of the scheme. Here, we have chosen the classical minmod limiter (the reader is referred for instance to [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] for more details regarding the use of slope limiters and a wider range of limiters). The limited slope σ n i is thus given by

σ n i = minmod w n i+1 -w n i ∆x , w n i -w n i-1 ∆x ,
where the minmod function is defined by

minmod(a, b) =      a if |a| < |b| and ab > 0, b if |a| > |b| and ab > 0, 0 if ab ≤ 0.
Thus, at the inner interfaces of the cell (x i-1 2 , x i+ 1 2 ), the reconstructed variables are given by

w - i := w n i x i - ∆x 2 = w n i - ∆x 2 σ n i , w + i := w n i x i + ∆x 2 = w n i + ∆x 2 σ n i . (6.1)
Since w ∈ {h, q, h+Z}, the reconstructed value of Z at the interfaces is computed from the reconstructed values of h + Z and h at the interfaces. Therefore, the updated states are given by the following three step scheme.

• Transport step:

h n+ 1 3 i = h n i - ∆t ∆x (f h ) n i+ 1 2 -(f h ) n i-1 2 , q n+ 1 3 i = q n i - ∆t ∆x (f q ) n i+ 1 2 -(f q ) n i-1 2 , (6.2)
where

(f h ) n i+ 1 2 := f h (W + i , W - i+1 ) and (f q ) n i+ 1 2 := f q (W + i , W - i+1
) according to (5.5). The water height h

n+ 1 3 i is then reconstructed within the cell (x i-1 2 , x i+ 1
2 ), following (6.1), to get h - i and h + i .

• Topography step:

h n+ 2 3 i = h n+ 1 3 i , q n+ 2 3 i = q n+ 1 3 i + ∆t(S t ) n+ 2 3 i , (S t ) n+ 2 3 i = 1 2 S t (h + i-1 , h - i , Z + i-1 , Z - i ) + S t (h + i , h - i+1 , Z + i , Z - i+1 ) , (6.3)
where S t is defined by (4.7).

• Friction step:

h n+1 i = h n+ 2 3 i , q n+1 i = (h η ) n+1 i η q n+ 2 3 i (h η ) n+1 i η + k ∆t q n+ 2 3 i , (h η ) n+1 i = 2µ n+ 2 3 i µ n i h -η n+1 i-1 2 + h -η n+1 i+ 1 2 + k ∆t µ n+ 2 3 i q n i , (6.4)
where

(h -η ) n+1 i-1 2 := h -η (h + i , h - i+1
), according to (3.25).

Finally, the scheme's time accuracy is improved by the use of the classical Heun's method (see [START_REF] Gottlieb | Strong stability-preserving highorder time discretization methods[END_REF]).

MOOD technique to recover the well-balance property

At this level, the designed MUSCL scheme is not well-balanced. Indeed, if W n i and W n i+1 define a steady state, the reconstruction step (6.1) will provide reconstructed values W + i and W - i+1 that no longer define a steady state. Thus, the scheme will not be able to exactly preserve the initial steady state, but will provide a second-order approximation. To restore this essential property, we use a MOOD-like technique (see [START_REF] Clain | A high-order finite volume method for systems of conservation laws-Multi-dimensional Optimal Order Detection (MOOD)[END_REF] for an overview of such techniques, and [START_REF] Berthon | An entropy preserving MOOD scheme for the Euler equations[END_REF] for more recent applications). We suggest to introduce a convex combination between the reconstructed state and the non-reconstructed one (see [START_REF] Hu | Positivity-preserving method for high-order conservative schemes solving compressible Euler equations[END_REF] for related work). As a consequence, we adopt the following reconstruction:

w - i = (1 -θ n i ) w n i + θ n i w n i - ∆x 2 σ n i = w n i - ∆x 2 σ n i θ n i , w + i = (1 -θ n i ) w n i + θ n i w n i + ∆x 2 σ n i = w n i + ∆x 2 σ n i θ n i , (6.5)
where 0 ≤ θ n i ≤ 1 is the parameter of the convex combination. If θ n i = 1, the full MUSCL scheme is recovered. If θ n i = 0, the states are not reconstructed and the first-order well-balanced scheme is used. Now, we propose a suitable process to define the parameter θ n i . To that end, we first define

∆ψ n i+ 1 2 = (q n i+1 ) 2 h n i+1 - (q n i ) 2 h n i + g 2 (h n i+1 ) 2 -(h n i ) 2 -∆xS t i+ 1 2 -∆xS f i+ 1 2 , where S t i+ 1 2 = S t (h n i , h n i+1 , Z i , Z i+1 ) and S f i+ 1 2 = S f (h n i , h n i+1 , q i , q i+1
), with S t and S f respectively defined by (4.7) and (4.8). Note that, from (4.3), ∆ψ n i+ 1 2 vanishes when W n i and W n i+1 define a steady state. Thus, from this quantity, we define a function to evaluate the deviation with respect to the equilibrium, as follows:

ϕ n i = q n i -q n i-1 ∆ψ n i-1 2 2 + q n i+1 -q n i ∆ψ n i+ 1 2 2 . (6.6) Since ∆ψ n i+ 1 2
vanishes when W n i and W n i+1 define a steady state, this newly introduced quantity ϕ n i vanishes when W n i-1 , W n i and W n i+1 define a steady state. Let M > m > 0. We define the parameter of the convex combination θ n i as follows:

θ n i =          0 if ϕ n i < m ∆x ϕ n i -m∆x M ∆x -m∆x if m ∆x ≤ ϕ n i ≤ M ∆x 1 if ϕ n i > M ∆x. (6.7)
This definition of θ n i enforces the use of the MUSCL scheme if the states are far from defining a steady state, i.e. ϕ n i is large enough; the first-order wellbalanced scheme is used if the equilibrium error ϕ n i is small enough. In addition, the closer the states are to the equilibrium, the more the convex combination will favor the first-order well-balanced scheme.

Numerical experiments

This last section is devoted to numerical tests, in one and two dimensions. We start by recalling the different schemes we shall test:

• the explicit scheme is (3.6) -(4.9);

• the implicit scheme is (5.8) -(5.9) -(5.12) -(5.18);

• the MUSCL scheme is (6.2) -(6.3) -(6.4) -(6.5).

In order to assess the properties of these schemes, we present four sets of numerical experiments. The first set assesses the well-balance of the scheme, by considering steady states at rest and moving steady states with friction, over flat and non-flat topographies. Then, the goal of the second set is to study the influence of the parameters C, m and M on the numerical solution. Afterwards, the focus of the third set is the simulation of wet and dry dam-breaks, over a possibly complicated topography, and with a nonzero Manning coefficient k. Finally, the fourth set of experiments concerns simulations in two dimensions.

In order to provide reference solutions when no exact solution is available, we elect to use the second-order hydrostatic reconstruction scheme (HR) presented in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF], equipped with the HLL flux and an implicit treatment of the friction source term. Note that this scheme preserves the steady states at rest, even with a nonzero friction contribution.

For the numerical experiments involving friction, a suitable value of the Manning coefficient k has to be chosen. For instance, the reader is referred to [START_REF] Chow | McGraw-Hill civil engineering series[END_REF] (page 109), where multiple values of k are given for different types of channel beds. Nevertheless, we deliberately impose stronger Manning coefficients than in reality (up to 10 times). This choice is made to ensure that the friction source term is preponderant compared to the topography source term, in order to enhance the effects of the friction.

An important step in these numerical experiments in the choice of the constant C, introduced in (3.26) to ensure the consistency of the approximate topography source term, and the choice of the constants m and M , introduced in (6.7) to define a steady state detector for the MUSCL scheme. In this article, these constants are chosen heuristically; a deeper study of the stability and of the order of accuracy of the scheme could provide several bounds for these constants. For each set of experiments, we will provide an insight regarding the choice of C, m and M .

The other constants are chosen according to Table 1. Moreover, we will compute errors to make sure the schemes satisfy the required properties. The L 1 , L 2 and L ∞ errors for a bounded function w are computed as follows:

Constant Equation Value

g (1.1) g = 9.81 m.s -2 ε λ (3.3) ε λ = 10 -10 m.s -1 ε (3.33) ε = 0 m
L 1 : 1 N N i=1 |w i -w ex i | ; L 2 : 1 N N i=1 (w i -w ex i ) 2 ; L ∞ : max 1≤i≤N |w i -w ex i |, (7.1
) where w i and w ex i are respectively the approximate and the exact solution at the point x i and at the final physical time t end , and N is the number of discretization cells. Finally, we recall that the CFL condition (3.5) gives the time step ∆t at each iteration, as follows:

∆t ≤ ∆x 2Λ , where Λ = max i∈Z -λ L i+ 1 2 , λ R i+ 1 2 .

Verification of the well-balance

In this first set of experiments, we assess the well-balance of the scheme, i.e. its ability to exactly preserve and capture steady states. Recall that steady states are given by the equation (1.5), which prescribes a uniform discharge over the space domain, denoted q 0 . First, steady states for the friction source term only are studied, that is to say we impose q 0 = 0 and ∂ x Z = 0. Then, we consider steady states for both friction and topography, which are either analytic solutions in specific cases or steady states obtained by approximately solving (1.5).

Steady states for the friction source term

We first focus on the preservation of the friction-only steady states. Thus, we assume q 0 = 0 and a flat topography, i.e. ∂ x Z = 0. The smooth steady states are then given according to (1.10), or equivalently (2.2). The space domain for these experiments is [0, 1], and its discretization is made of 200 cells (i.e. ∆x = 5.10 -3 ).

In Section 2, the water height for a smooth steady state was obtained by considering a zero of the nonlinear function ξ defined by (2.5). We consider the subcritical branch of the steady state obtained by assuming q 0 = -1, x 0 = -∆x and h 0 = h c , with h c defined by (2.6). The Manning coefficient k is chosen equal to 1.

The first experiment concerns the preservation of this steady state. We take q(x) = q 0 and the exact height as initial conditions, and we impose the exact solution at the points -∆x and 1 + ∆x as inhomogeneous Dirichlet boundary conditions. We compute the approximate solution until the final time t end = 1s.

Moreover, the parameters are taken as follows: C = 10 -3 , m = 10 -12 and M = 10 -11 . Such a value for C was chosen to ensur that the term [h] c present in the approximate source terms was small enough. The justification of the choice of the parameters m and M is the same as in the case of the steady states are rest.

The results are presented on Figure 4, and the errors to the steady state are displayed in Table 2. h q From Figure 4 and Table 2, we observe that the friction-only steady state is indeed preserved up to the machine precision by the three schemes, which is not the case when considering non-well-balanced schemes (see for instance [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF] in the context of the shallow-water equations with topography only).

L 1 L 2 L ∞ L 1 L 2 L ∞ explicit 6.00e-
The second experiment is based on the previous steady state, but with a perturbation, as shown on Figure 5. With h ex the exact height, this perturbation is defined by choosing the initial water height as follows:

h(x) =      h ex (x) + 0.2 if x ∈ 3 7 , 4 7
;

h ex (x) otherwise.
The initial discharge is unperturbed, and taken equal to q 0 = -1 throughout the domain. The boundary conditions consist in the unperturbed exact solution. We use the three schemes and 100 discretization cells for the numerical simulation on the domain [0, 1]. Moreover, we take C = 10 -3 , m = 10 -12 and M = 10 -1 . For this perturbed initial data, we have raised the value of the parameter M so that the scheme exactly converges towards the unperturbed steady solution. Indeed, the wellbalanced scheme starts being used as soon as the steady state deviation, defined by (6.6), is lower than M ∆x. Therefore, the value M = 10 -1 that the well-balanced scheme is used soon enough to allow the exact recovery of the unperturbed steady state.

The computations are carried out until the final time t end = 9s. Indeed, such a final time allows the perturbation to be dissipated and a steady state to be reached. In fact, this steady state is the original, unperturbed steady state. The results of the explicit scheme are presented on Figure 5 and an error comparison with the unperturbed steady state is provided in Table 3. h q Figure 5 shows that the perturbation is eventually dissipated and we recover the unperturbed steady state. This assertion is confirmed by the error analysis presented in Table 3, which shows that all three schemes recover the unperturbed steady state up to the machine precision.

L 1 L 2 L ∞ L 1 L 2 L ∞
Remark 5. We have carried out similar experiments with a supercritical branch instead of the subcritical branch, and we obtained the same results, i.e. the steady state solution was exactly recovered. For the sake of conciseness, we do not present them here.

Steady states with both friction and topography

We end the sanity check of the well-balance property with numerical experiments consisting in the preservation of steady states involving both topography and friction (for experiments with the sole topography source term, including the preservation of steady states at rest, the reader is referred to [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF]). Hence, we take k = 0, ∂ x Z = 0, and q 0 = 0. Steady states are therefore given by the full equation (1.5). Recall that this equation cannot be rewritten under an algebraic form. Thus, to find a general steady state solution, we have to numerically solve the equation (1.5).

We by considering two very specific cases, where the height or the free surface are constant, in order to obtain an analytical solution of (1.5).

First, assume h(x) = h 0 , i.e. the water height is uniform and equal to h 0 in the whole space domain. Thus, (1.7) rewrites

gh 0 ∂ x Z + kq 0 |q 0 |h η 0 = 0,
and we immediately obtain a relation characterizing the slope of the topography:

∂ x Z = - kq 0 |q 0 | gh η+1 0 . (7.2)
This specific case is tested numerically by taking h 0 = q 0 = 1, k = 10, Z(0) = 0 and the slope of Z given by (7.2). The space domain is [0, 1] and is discretized using 100 cells. The initial and boundary conditions are the exact solution. The computations are carried out with all three schemes, and we choose t end = 1s.

We take C = +∞, m = 10 -12 and M = 10 -11 , for reasons similar to the ones invoked for the first experiment with friction only. Table 4 shows that the topography and friction steady state with constant height is indeed exactly preserved by the three schemes. To build another exact solution of (1.5), we assume h + Z = H 0 instead of h = h 0 . We therefore have a constant free surface H 0 over the whole space domain [0, 1]. Note that (1.7) rewrites for smooth solutions with positive water height as follows:

h q L 1 L 2 L ∞ L 1 L 2 L ∞ explicit 9.
- q 2 0 h 2 ∂ x h + gh∂ x (h + Z) + kq 0 |q 0 | h η = 0.
Since the free surface h + Z is constant, the above relation rewrites

h η-2 ∂ x h = kµ 0 , (7.3)
where µ 0 = sgn(q 0 ). For x ∈ [0, 1], integrating (7.3) over [0, x] provides an expression of h that depends on the initial condition h 0 := h(0), namely:

h η-1 (x) = h η-1 0 + (η -1) kµ 0 x. (7.4)
Note that (7.4) defines an expression of the height that may become non-positive for some x. Thus, we choose the constants k, q 0 h 0 such that h is positive over the whole domain [0, 1]. In the simulation, we set k = q 0 = h 0 = 1 and H 0 = 1. Equipped with the water height, we define the topography by setting Z(x) = H 0 -h(x). We discretize the space domain with 100 cells and use the three schemes to carry out the simulation until the time t end = 1s. We again take C = +∞, m = 10 -12 and M = 10 -11 . An error comparison of the three schemes is displayed in Table 5, which shows that the steady state is preserved up to the machine precision by the explicit, implicit and MUSCL schemes. Finally, we derive a steady state for the shallow-water equations with topography and friction, without considering a constant height or free surface. Thus, we approximately solve the discretization (4.3) of the full equilibrium relation (1.7). First, we set k = 0.01 and choose [0, 1] to be the space domain. Moreover, the topography is given by

h q L 1 L 2 L ∞ L 1 L 2 L ∞ explicit 4.
Z(x) = 1 2
e cos(4πx) -e -1 e 1 -e -1 . (

We set q(x) = q 0 throughout the domain. The equation (4.3) is then approximately solved using Newton's method, imposing h(0) = 0.3. This procedure allows us to define the water height over the whole domain. This steady state t (h(x), q(x)) is then chosen as the initial and boundary conditions for this experiment. We take 100 discretization cells, and we once again set C = +∞, m = 10 -12 and M = 10 -11 . The numerical simulation of this experiment is done using all three schemes, and runs until a final physical time t end = 1s. The results of the explicit scheme are presented on Figure 6 and the errors to the steady state are displayed in Table 6.

Figure 6 shows that the explicit scheme exactly preserves this topography and friction steady state. Moreover, Table 6 shows that both other schemes also exactly preserve such a steady state.

The very last experiment in the assessment of the well-balance property focuses on a perturbation of the aforementioned steady state. We denote by h ex (x) the water height of the previous topography and friction steady state. h q The initial water height of this last experiment is then given by

L 1 L 2 L ∞ L 1 L 2 L ∞
h(x) =      h ex (x) + 0.05 if x ∈ 2 7 , 3 7 4 7 , 5 7 , h ex (x) otherwise,
and the discharge is defined by

q(x) =      q 0 + 1 2 if x ∈ 2 7 , 3 7 4 7 , 5 7 , q 0 otherwise.
The unperturbed steady state is prescribed as the boundary conditions. For this numerical experiment, the domain [0, 1] is discretized with 100 cells, and the simulation runs until the perturbation has been dissipated and the unperturbed steady state has been recovered. The final physical time we choose for these conditions to be met is t end = 2s. Moreover, since we seek the capture of this perturbed steady state, we set the parameters as C = +∞, m = 10 -13 and M = 0.05. The evolution of the perturbation with the explicit scheme is presented on Figure 7. Then, in Table 7, we present the errors to the original unperturbed steady state when the physical time is elapsed. h q 7 show that the scheme indeed allows to recover the original unperturbed steady state. This experiment emphasizes the ability of the three schemes to exactly capture a steady state, even after a perturbation.

L 1 L 2 L ∞ L 1 L 2 L ∞ explicit 5.

Influence of the parameters C, m and M

This second set of experiments is devoted to the study of the influence of the parameters present in the scheme: C, m and M . Our goal is to emphasize the need for a relevant value of these parameters in numerical simulations.

We first perform three dam-break experiments whose goals are to show the dependence in C of the numerical solution, when considering the topography or the friction source term. Then, the dependence in m and M , related to the order of the numerical scheme, is studied.

Influence of the parameter C for the topography

Recall that the approximate topography source term S t , defined by (4.7), depends on the parameter C. Namely, this parameter is used to ensure the consistency of the approximate source term S t with the actual source term S t , by ensuring that the absolute value of the water height jump [h] is no larger than C∆x. The purpose of the first set of experiments is to highlight the influence of this parameter C.

Shock waves over a flat topography

The first experiment we consider concerns the propagation of shock waves over a flat topography. To that end, we consider the domain [0, 1] with a flat topography (i.e. Z ≡ 0), homogeneous Neumann boundary conditions, and we take the following dam-break initial data: h(0, x) = 1 and q(0, x) = 7.5 if x < 0.5,

-7 if x ≥ 0.5.
Recall that the cutoff involving C had been introduced to ensure the consistency of S t , especially for flat topographies. As a consequence, for this experiment, the value of C should be instrumental in getting the correct shocks waves and an accurate approximation of the intermediate state.

To carry out this experiment, we set a vanishing friction contribution (i.e. k = 0), we use 250 discretization cells, and we take the final time t end = 0.1s. This experiment is performed for C = 10 and C = 1000. The results of the implicit scheme are presented on Figure 8, where we display the exact solution (see for instance [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF][START_REF] Toro | Riemann solvers and numerical methods dynamics. A practical introduction[END_REF] for more details) as well as the numerical approximation. The left panel of Figure 8 shows that, for both values of C, the approximate shock waves are located at a consistent position and seem to have the correct amplitude. However, on the right panel of Figure 8, we note, on the one hand, that the intermediate state obtained with C = 1000 presents spurious oscillations, whose amplitude does not decrease when ∆x decreases. On the other hand, with C = 10, the implicit scheme provides a good approximation of the intermediate state.

Dam-break on an emerging bottom

The second dam-break experiment we perform in order to study the influence of the parameter C presents an emerging bottom. To that end, we modify the flow at rest with emerging bottom experiment from [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF], to add a wave perturbing the water at rest. On the space domain [0, 15], we consider the topography function Z(x) = 0.2 -0.05(x -10) 2 + and the following dambreak initial conditions:

h(0, x) + Z(x) = 0.2 if x < 5, 0.15 if x ≥ 5,
and q(0, x) = 0.

We prescribe homogeneous Neumann boundary conditions, and we consider 4000 discretization cells. We once again take a contribution, i.e. k = 0. The initial condition, as well as the reference solution, computed with the HR scheme at t end = 5s, are displayed on Figure 9. Reference free surface obtained with the HR scheme. On both panels, the gray area is the topography.

In order to study the influence of the parameter C, we consider the implicit scheme, and we carry out the simulation with C = 1 and C = 10. The results are displayed on Figure 10. Indeed, with C = 10, large spurious oscillations appear, a wave is reflected to the left of the domain, and an inconsistent dry area appears on the bump. However, with C = 1, we only get small oscillations near x = 9.75; the amplitude of these oscillations decreases as ∆x decreases. On Figure 10, we once again observe that the consistency is ensured by the cutoff procedure. Indeed, with C = 10, large spurious oscillations appear, a wave is reflected to the left of the domain, and an inconsistent dry area appears on the bump. However, with C = 1, we only get small oscillations near x = 9.75; the amplitude of these oscillations correctly decreases as ∆x decreases.

Influence of the parameter C for the friction

Next, we turn to the study of the influence of the parameter C on the approximate friction source term S f = -kq|q|h -η , defined (3.27). that the purpose of this parameter is to ensure the consistency of S f with the friction source term S f , especially when dealing with discontinuities and a small Manning coefficient k.

In order to highlight the relevance of this parameter, we consider the following numerical experiment. We consider a dam-break problem on a flat topography (i.e. Z ≡ 0) with a very small Manning coefficient (we take k = 10 -9 ). As a consequence, the first term of the expression (3.27) of S f is negligible, and the presence of the cutoff (3.26) makes its second term consistent. The initial conditions of the dam-break problem are taken as follows:

h(0, x) = 20 if x < 0.5, 30 if x ≥ 0.5. and q(0, x) = 250 if x < 0.5, -200 if x ≥ 0.5.
and the space domain is [0, 1], with homogeneous Neumann boundary conditions. Equipped with such initial conditions, we expect two large shock waves to propagate in the domain. This experiment is carried out with the explicit scheme, until the final physical time t end = 0.0175s. Note that we know the exact solution t (h exact , q exact ) to this dam-break problem for k = 0 (see for instance [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]). As a consequence, in order to quantify the error introduced by the inconsistent term and hopefully corrected by the use of C, we consider the following error:

∆q := i∈ 1,N (q n end i ∆x) - 1 0 q exact (t end , x) dx ,
where N is the number of cells and where t n end = t end .

To that end, we compute the errors ∆q, and the associated orders of accuracy, for the following three configurations: 8.

(a) For C = 0, we get a consistent scheme. The values of the error presented in Table 8 are of the same order of magnitude as k, which was to be expected since the approximate solution is compared to the exact solution obtained with k = 0. This is verified by taking, for instance, k = 10 -11 : in this case, ∆q 1.2 × 10 -13 for all N . However, if a steady state solution were considered, the value C = 0 would not allow its exact preservation. (b) For C = 1, since the error is of the order of magnitude of ∆x, we also get a consistent scheme. Indeed, in this case, the formerly inconsistent term of (3.27) becomes a perturbation of size O(∆x). (c) For C = +∞, the scheme is inconsistent: indeed, the error ∆q does decrease with ∆x, but the approximate solution does not converge towards the correct solution. However, this value of C allows steady state solutions to be exactly preserved. This experiment therefore highlights the necessity of the cutoff introduced in (3.27) in order to recover the consistency of the scheme when considering discontinuous solutions and a negligible friction contribution, while ensuring that the scheme is still well-balanced.

Influence of the parameters m and M

The goal of this next experiment is to emphasize the dependence of the numerical solution in m and M . Recall that these parameters were introduced during the MOOD procedure to ensure the well-balance of the scheme. Namely, they were used to define, in (6.7), the parameter θ n i of the convex combination between the non-reconstructed states and the reconstructed ones.

In order to study the influence of these parameters, we consider an experiment devoted to computing the order of accuracy of the scheme. To set up this experiment, we compute a steady solution with both topography and friction contributions. We set k = 0.01 and a topography given by (7.5) over the space domain [0.4, 0.6]. We take the steady discharge q(x) = q 0 = 1 and we approximately solve (4.3) with h(0) = 1 to finally obtain the steady water height h(x) on the whole space domain. This steady state t (h(x), q(x)) is then chosen as the initial and boundary conditions for this experiment.

However, this steady state is exactly preserved by the well-balanced scheme. In order to study the order of accuracy of the scheme, we introduce the following smooth perturbation, of amplitude 10 -6 , of the steady solution on [0.45, 0.55]:

ash p (x) = h(x) + 10 -6 × ω 2 0.05 (x -0.5) ,
where we have set

ω(z) =      2 -|z| 2 4 (1 + 2|z|) if |z| < 2, 0 otherwise.
The initial data of the experiment is perturbed steady state t (h p (x), q(x)). Note that such initial data is smooth (of class C 4 ) on [0.4, 0.6], and it is a steady solution on [0.4, 0.45] ∪ [0.55, 0.6]. As a consequence, on the one hand, a well-balanced scheme will exactly preserve its steady part, and provide a firstorder approximation of the unsteady part. On the other hand, a second-order scheme will provide a second-order approximation of the solution everywhere. Therefore, the advantage of using the aforementioned MOOD procedure lies in ensuring both the exact preservation of the steady part and the second-order approximation of the unsteady part. This experiment is carried out until the final physical time t end = 10 -3 . Such a final time has been chosen to ensure that the perturbation does not reach the boundaries of the domain. We use the MUSCL scheme for several values of m and M , in order to study the order of accuracy for each pair (m, M ) under consideration. To that end, we first compute a reference solution using 20480 cells. Then, to obtain the order of accuracy, we begin by using (7.1) to compute the L 1 , L 2 and L ∞ errors. With e N representing either of these errors for N cells, the order of accuracy ρ is computed as follows: From Table 9, Table 10, Table 11 and Figure 11, we note that the scheme is second-order accurate when using the pairs (a) and (b), and first-order accurate for the pair (c). These observations can be explained as follows: for the first pair, the second-order reconstruction is used almost everywhere, which leads to the scheme being second-order accurate. For the third pair, the secondorder reconstruction is never used, and the scheme is first-order accurate. It is interesting to remark that, for the coarsest grid, the first-order well-balanced scheme is more accurate than the second-order scheme. The second pair leads to a balance between the second-order piecewise linear reconstruction and the first-order piecewise constant values. This balance ensures that the secondorder scheme is used where the solution is unsteady and that the well-balanced scheme is used where the solution is steady, thus leading to a scheme that is more accurate than in the other two cases. This experiment emphasizes that the choice of m and M is of crucial importance in order to ensure a good quality of the numerical approximation.

N L 1 L 2 L ∞ 2.

Dam-break experiments

Now, we focus on several dam-break experiments, whose goal is to provide a validation of the scheme on unsteady cases. We consider a dam located at x = D that instantly breaks at t = 0s, thus releasing the water it was holding. Such a dam-break experiment has initial data

W L = t (h L (x), 0) if x ≤ x D and W R = t (h R (x), 0) if x > x D .
The topography function Z and the Manning coefficient k are different in each experiment.

In addition, the parameters m and M are chosen to ensure the preservation of the situations of water at rest located away from the waves, while still retaining a more accurate approximation within the waves. More details on such water at rest situations are given in the next paragraph. Finally, the parameter C has been chosen heuristically, and it is reported in each experiment.

Wet dam-break

We first carry out a wet dam-break experiment. Here, the topography is given by Z(x) = cos 2 (2πx) /2 and the Manning coefficient is k = 1. The space domain is [0, 1] and the dam is located at x D = 0.5. We prescribe homogeneous Neumann boundary conditions at both boundaries. The initial data is h L (x) + Z(x) = 2 and h R (x) + Z(x) = 1. The results, obtained using 100 discretization cells and t end = 0.05s, are displayed on Figure 12. For this simulation, we choose C = 10, m = 10 -6 and M = 10 -1 . The goal of this experiment is to compare the three schemes to a reference solution, obtained using the HR scheme with 10 4 discretization cells. The three schemes yield a correct approximation of the reference solution. The explicit and the implicit schemes yield very close solutions, while the MUSCL scheme provides a better approximation of the reference solution. addition, we note that the intermediate state between the waves is correctly approximated by the implicit scheme, while the HLL scheme provides a poor approximation, which proves the importance of using a well-balanced scheme for such simulations.

Dry dam-break over a flat topography: explicit scheme disqualification

The goal of this is to show the spurious oscillations induced by the explicit scheme in the presence of the friction, and the improvement obtained by using the semi-implication of the friction in the implicit scheme.

The topography function is Z(x) = 0, and we take k = 5. The initial data is h L (x) = 1.5 and h R (x) = 0, and the space domain is [-1, 1], for a dam located at x D = 0. Hence, the wet/dry front is located within the interval [0, 1]. Figure 14 displays the discharge computed by the explicit and the implicit schemes on this interval. To obtain these results, we used 200 cells, the final physical time t end = 0.03s, and we took C = 1. Finally, we prescribe homogeneous Neumann boundary conditions on both boundaries, and we compute the reference solutions using the HR scheme with 2 × 10 4 cells.

Figure 14: Dry dam-break: instabilities induced by the explicit scheme. The discharge is represented for the explicit scheme (dotted line), the implicit scheme (dashed line) and the reference solution (solid line). Oscillations in the discharge produced by the explicit scheme are clearly visible.

From Figure 14, on the one hand, we note that the discharge obtained with the explicit scheme presents several non-physical oscillations. According to the reference solution, it should stay equal to zero. On the other hand, we remark that no oscillations appear when using the implicit scheme. As a consequence, two boundaries at x = 0m and x = 5m are solid walls. The topography is defined by

Z(x) = 1 2 1 - (x -5/2) 2 1/25 + + 2 1 - (x -4) 2 1/25 + ,
and it indeed consists in two quadratic bumps, a smaller one followed by a larger The dam is located at x D = 0.7m, breaks at t = 0s, and contains an initial water height h L = 6m. The domain x > x D contains no water, i.e. h R = 0. We take a Manning coefficient k equal to 1. Figure 16 shows that the initial water height is significantly larger than the bumps, in order to impact them with a high enough velocity. We choose C = 0.1, m = 10 -5 and M = 10 -1 .

Figure 16: Free for the double bump test case at different times: from left to right and top to bottom, the solution is observed at times 0.38s, 0.74s and 1.70s.

2D simulations

This subsection is devoted to experiments in two space dimensions. In this context of two-dimensional geometries, equations with topography and Manning friction read as follows:

             ∂ t h + ∂ x p + ∂ y q = 0, ∂ t p + ∂ x p 2 h + 1 2 gh 2 + ∂ y pq h = -gh∂ x Z -kpDh -η , ∂ t q + ∂ x pq h + ∂ y q 2 h + 1 2 gh 2 = -gh∂ y Z -kqDh -η , (7.7)
where p is the discharge in the x direction, q is the discharge in the y direction, and D is the Euclidean norm of the discharge vector, i.e. D = p 2 + q 2 . We rewrite (7.7) under the following condensed form:

∂ t W + ∂ x f (W ) + ∂ y g(W ) = s(W ),
where the definitions of W , f , g and s can easily be inferred from (7.7).

The approximation involves a uniform Cartesian grid, made of square cells with sides ∆x = ∆y. For the 2D simulations, we will consider 2D variants of the implicit and MUSCL scheme, by applying the 1D strategy detailed in Section 5 and Section 6 in both directions, horizontal and vertical, of each square cell. This process involves exactly solving the initial value problem (5.10), related to the contribution of the friction only. However, this initial value problem is modified in two dimensions. Indeed, we to solve

             dh dt = 0, dp dt = -k(h n+1 i,j ) -η p p 2 + q 2 , dq dt = -k(h n+1 i,j ) -η q p 2 + q 2 , with initial data          h(0) = h n+ 2 3 i,j , p(0) = p n+ 2 3 i,j , q(0) = q n+ 2 3 i,j . (7.8)
We immediately see that h n+1 i,j = h n+ 2 3 i,j . For the sake of simplicity in the notations, we set τ = -k(h n+1 i,j ) -η < 0, p n = p n+ 2 3 i,j

and q n = q n+ 2 3 i,j . Thus, (7.8) rewrites as the following initial value problem: p = τ p p 2 + q 2 , q = τ q p 2 + q 2 , with initial data p(0) = p n , q(0) = q n .

The exact solution to this system is:

p(t) = p n h n+1 i,j h n+1 i,j + kt p 2 n + q 2 n and q(t) = q n h n+1 i,j h n+1 i,j + kt p 2 n + q 2 n .
Therefore, the friction step of the implicit procedure yields the following updated discharges:

p n+1 i,j = p n+ 2 3 i,j (h n+1 i,j ) η (h n+1 i,j ) η + k∆t (p n+ 2 3 i,j ) 2 + (q n+ 2 3 i,j ) 2 , q n+1 i,j = q n+ 2 3 i,j (h n+1 i,j ) η (h n+1 i,j ) η + k∆t (p n+ 2 3 i,j ) 2 + (q n+ 2 3 i,j ) 2 .
The quantity (h n+1 i,j ) η is then replaced with the average given by (5.18), taken in the x direction for p n+1 i,j and in the y direction for q n+1 i,j . N h p q 900 1.33e-02 -2.54e-02 -3.50e-02 -3600 6.82e-03 0.97 1.36e-02 0.90 1.83e-02 0.94 14400 3.44e-03 0.99 7.02e-03 0.95 9.36e-03 0.96 57600 1.73e-03 0.99 3.57e-03 0.97 4.75e-03 0.98

2D dry dam-break with two bumps

The first 2D experiment is a dry dam-break with a topography presenting two bumps. This experiment is the 2D analogue of the 1D experiment presented in Section 7.3.4. It is heavily inspired from an experiment presented in [START_REF] Berthon | Efficient well-balanced hydrostatic upwind schemes for[END_REF], which did not include the friction source term. The Manning coefficient is k = 0.1, Note that the same color scale for the water height is used in Figure 17 and Figure 18, and the solid gray color represents the topography. Right panel: approximate solution at t = 0.15s, just before the water hits the first bump. Note the shape of the front of the water, due to the nonzero bottom friction.

As in the previous dam-break experiment, we take p(x, y, 0) = 0 and q(x, y, 0) = 0. Moreover, we define the initial water height by: h(x, y, 0) = 10 if x < -5, 5 otherwise.

From such initial data, we present two simulations. The first one is the case without friction, i.e. k = 0, and for the second one we choose a nonzero Manning coefficient k = 2. In both cases, homogeneous Neumann boundary conditions are prescribed on all the boundaries. We choose t end = 7s, and we set C = 0.5, m = 10 -10 and M = 0.5.

The justification of the parameter choice is similar to the one for the previous dam-break experiments, with one difference. Since the current experiment does not involve any water at rest with a non-flat topography, it is relevant to use the MUSCL scheme as much as possible, as evidenced by the relatively low value of m.

The numerical simulation of this experiment is carried out with the implicit and the MUSCL schemes on 202500 cells (450 in each direction). The results are Figure 19 shows that, when no friction is present, vortices appear at both tips of the dam. However, as evidenced by Figure 20, these vortices no longer exist when the friction is nonzero. Both behaviors are exhibited with the implicit and MUSCL schemes. In addition, the tip of the left shock wave has traveled further on Figure 19 than on Figure 20. This behavior is consistent with the fact that the friction tends to slow down the movement of the water. Now, to conclude the analysis of this experiment, we introduce (x T , y T ) and (x B , y B ). They are, respectively, the points where the top and bottom vortices are the deepest when k = 0. In Table 16, we present the height at these two points, obtained with the implicit and MUSCL schemes, and considering k = 0 or k = 2. In addition, we display in this table the approximate position x S of the left shock wave in the four cases under consideration, where it crosses the line y = 0.

Urban topography

The last 2D experiment is a simulation of a city being hit by a wave. We consider the space domain [0, 1000] × [0, 1000]. The topography consists in an upwards slope leading to a flat surface, upon which buildings are placed. Disregarding the buildings, the bottom has the following topography: Z(x, y) =

x/50 if x < 500, 10 otherwise.

The 100 meters high buildings occupy the flat part of the topography, i.e. buildings are only present for x > 500. Figure 21 displays the shapes and the positions of the buildings, for a uniform Cartesian mesh of 10 6 cells (1000 in each direction). The initial conditions are W (x, y, 0) = 0 for all x and y in the space domain. Indeed, the boundary conditions help create the flood and the wave that hits the city. We prescribe homogeneous Neumann boundary conditions for each boundary of the domain, except the left boundary, where a time-dependent boundary condition the x-discharge p is applied, as follows: p(0, y, t) = 15 if t < 350, ∂ x p(0, y, t) = 0 otherwise. (7.9) Such a boundary condition creates water that fills the sloping part of the topography and creates a wave that hits the city. To illustrate the situation, Figure 21 shows the free surface at time t = 300s, when the first buildings are about to be hit by the wave. We consider a nonzero Manning coefficient k = 1 and a final time t end = 850s. The simulation is run using the MUSCL scheme, and we take C = 10 -2 , m = 10 -5 and M = 1 for the same reasons as invoked before.

The left panel of Figure 21 shows the wave created by the Dirichlet boundary condition arriving in the city. Because of the friction, this wave presents a rather steep front. On the left panel of Figure 22, the wave has hit the first buildings located at the south of the city. Note that the space between the first two columns of buildings is still dry. Also note that, as per (7.9), the boundary condition imposed on the x-discharge q x is now a homogeneous Neumann boundary condition, and no more water is injected into the domain. The right panel of Figure 22 displays the wave about to hit the square building located at the middle of the city. As expected, between the southern buildings, the wave is slowed down.

On the left panel of Figure 23, the wave has reflected on the southwestern side of the square building, and it has thus moved faster towards the buildings to the south. On the right panel of Figure 23, the waves reflected from the square building are moving south and north. Moreover, the back of the "S"-shaped building will soon be flooded. The flooding of the back of the "S"-shaped building is happening on Figure 24, with only a small area still dry on the right panel. In addition, on the right panel of Figure 24, the wave has almost hit the small square building on the bottom right of the city.

Figure 25 and Figure 26 display the final phases of the flooding of the city. Note that the southern buildings are mostly uniformly flooded and that the inner courtyard of the square building is still dry. Moreover, the water at the back of the "S"-shaped building is less deep than at other points of the same vertical line.

Conclusion

In this work, we first studied smooth steady state solutions for the shallowwater equations with friction only. Then, we have proposed a Godunov-type scheme for the same equations that is well-balanced. This scheme has then been extended to add the topography contribution as the positivity preservation. Afterwards, we have proposed a semi-implicit scheme to deal with dry/wet interfaces, followed by a well-balanced second-order MUSCL extension. Finally, we have proposed several one-dimensional and two-dimensional numerical experiments to assess the efficiency and robustness of the explicit, semi-implicit and MUSCL scheme.

One issue that remains is the division by α f occurring when computing the intermediate water heights (see Remark 4). Further work could be undertaken in order to try and provide a satisfactory answer to this issue.

Figure 2 :

 2 Figure 2: Structure of the chosen approximate Riemann solver.

Lemma 5 .

 5 The intermediate states (4.6) satisfy the following properties:(i) consistency with the shallow-water equations with topography and friction (1.1);

Figure 4 :

 4 Figure 4: Left panel: initial height for the friction steady state. Right panel: height (solid line) and discharge (dashed line) errors to the steady state after 1s, with the explicit scheme.

Figure 5 :

 5 Figure 5: From left to right: water height for t = 0s, t = 0.06s and t = 9s, with the explicit scheme.

explicit 7 .

 7 42e-15 7.48e-15 8.55e-15 7.06e-15 8.32e-15 1.48e-14 implicit 7.06e-15 7.09e-15 8.10e-15 7.61e-15 8.81e-15 1.51e-14 MUSCL 8.40e-15 8.47e-15 1.04e-14 6.56e-15 7.52e-15 1.73e-14

Figure 6 :

 6 Figure 6: Left panel: initial height for the topography and friction steady state. Right panel: height (solid line) and discharge (dashed line) errors to the steady state with the explicit scheme.

explicit 3 .

 3 54e-16 4.49e-16 1.11e-15 1.66e-15 2.16e-15 3.55e-15 implicit 5.73e-16 8.99e-16 3.00e-15 5.07e-16 6.29e-16 1.55e-15 MUSCL 3.52e-17 5.40e-17 1.67e-16 3.75e-16 4.85e-16 1.33e-15

Figure 7 :

 7 Figure 7: From left to right: water height for t = 0s, t = 0.015s and t = 2s, with the explicit scheme.

Figure 7

 7 Figure 7 and Table7show that the scheme indeed allows to recover the original unperturbed steady state. This experiment emphasizes the ability of the three schemes to exactly capture a steady state, even after a perturbation.

Figure

  Figure Dam-break creating two shock waves over a flat bottom. Left panel: whole domain depicted at t = 0.1s. Right panel: zoom on the intermediate state of the dam-break problem.

Figure 9 :

 9 Figure 9: Dam-break on an emerging bottom. Left panel: Initial free surface. Right panel:Reference free surface obtained with the HR scheme. On both panels, the gray area is the topography.

Figure 10 :

 10 Figure 10: Dam-break on an emerging bottom: zoomed comparison between the HR scheme, the implicit scheme with C = 1, and the implicit scheme with C = 10. The gray area is still the topography.

  (a) C = 0; (b) C = 1; (c) C = +∞.The errors and orders of accuracy are collected in Table

6 )

 6 Equipped with this definition, we compute the L 1 , L 2 and L ∞ discharge errors and the associated orders of accuracy for the following pairs (m, M ): (a) m = 10 -12 and M = 10 -11 ; (b) m = 10 -9 and M = 10 -6 ; (c) m = 10 -5 and M = 10 -2 . The results are collected in

Figure 11 :

 11 Figure 11: Error plots in L 2 -norm for the three (m, M ) pairs.

Figure 12 :

 12 Figure 12: Free surface observed at the final physical time with the three schemes, for the wet dam-break experiment.

Figure 17 :

 17 Figure 17: Left panel: initial condition of the 2D dam-break over a double bump experiment.Note that the same color scale for the water height is used in Figure17and Figure18, and the solid gray color represents the topography. Right panel: approximate solution at t = 0.15s, just before the water hits the first bump. Note the shape of the front of the water, due to the nonzero bottom friction.

Figure 18 :

 18 Figure18: From left to right and top to bottom: approximate solution of the 2D dam-break over a double bump experiment, displayed at times t = 0.3s, t = 0.45s, t = 0.75s, t = 0.9s, t = 1.05s and t = 1.35s.

Figure 19 :

 19 Figure 19: Approximate free surface for the partial dam-break experiment without friction. Left panel: implicit scheme; right panel: MUSCL scheme.

Figure

  Figure Approximate surface for the partial dam-break experiment with nonzero friction. Left panel: implicit scheme; right panel: MUSCL scheme.

Figure 21 :

 21 Figure 21: Wave on an urban topography Left panel: topography of the city. buildings are actually 100 meters high, and are represented in white in this figure. One can see the upwards slope on the left, leading to the city itself. Right panel: free surface at t = 300s. The wave is present to the left of the figure. Note that the same free surface color scale will be used in the next figures.

Figure 22 :

 22 Figure 22: Free surface for the urban topography simulation at t = 355s (left panel) and t = 410s (right panel).

Figure 23 :

 23 Figure 23: Free surface for the urban at t = 465s (left panel) and t = 520s (right panel).

Figure 24 :

 24 Figure 24: Free surface for the urban topography simulation at t = 575s (left panel) and t = 630s (right panel).

Figure 25 :

 25 Figure 25: Free surface for the urban simulation at t = 685s (left and t = 740s (right panel).

Figure 26 :

 26 Figure 26: Free surface for the urban topography simulation at t = 795s (left panel) and t = 850s (right panel).

Table 1 :

 1 Values of the constants within the numerical experiments.

Table 2 :

 2 Height and discharge errors with the three schemes for the friction steady state.

Table 3 :

 3 Height and discharge errors with the three schemes for the perturbed friction steady state.

Table 4 :

 4 Height and discharge errors for the topography and friction steady state with constant height.

  23e-15 4.48e-15 6.66e-15 3.06e-14 3.07e-14 3.60e-14 implicit 1.31e-16 2.00e-16 6.66e-16 1.45e-15 1.83e-15 4.88e-15 MUSCL 8.88e-17 1.57e-16 4.44e-16 4.46e-16 6.81e-16 2.22e-15

Table 5 :

 5 Height and discharge errors for the topography and friction steady state with constant free surface.

Table 6 :

 6 Height and discharge errors with the three schemes for the topography and friction steady state.

  10e-16 7.09e-16 2.44e-15 2.20e-15 2.92e-15 5.33e-15 implicit 4.85e-16 7.97e-16 2.72e-15 9.73e-16 1.15e-15 2.78e-15 MUSCL 6.46e-16 9.40e-16 3.50e-15 1.71e-15 2.23e-15 4.22e-15

Table 7 :

 7 Height and discharge errors with the three schemes for the perturbed topography and friction steady state.

Table 8 :

 8 Dam-break experiment with large shocks: influence of the parameter C.

	N	∆q, C = 0	∆q, C = 1	∆q, C = +∞
	200	2.44e-10 -1.64e-04	-	2.21e-01	-
	400	2.41e-10 -8.03e-05 1.03 2.12e-01 0.06
	800	2.40e-10 -3.98e-05 1.01 2.08e-01 0.03
	1600 2.39e-10 -1.98e-05 1.01 2.06e-01 0.01

Table 9 ,

 9 Table 10 and Table 11, as well as Figure 11.

	N	L 1		L 2		L ∞
	80	3.34e-07	-	3.57e-07	-	4.60e-07	-
	160	8.39e-08 1.99 8.97e-08 1.99 1.16e-07 1.99
	320	2.11e-08 1.99 2.25e-08 1.99 2.90e-08 1.99
	640	5.27e-09 2.00 5.63e-09 2.00 7.26e-09 2.00
	1280 1.31e-09 2.01 1.39e-09 2.01 1.79e-09 2.02

Table 9 :

 9 Discharge errors and order of accuracy for the pair (a).

Table 10 :

 10 Discharge errors and order of accuracy for the pair (b).

	N	L 1		L 2		L ∞
	80	2.00e-07	-	3.33e-07	-	7.74e-07	-
	160	9.98e-08 1.01 1.67e-07 1.00 3.92e-07 0.98
	320	4.96e-08 1.01 8.31e-08 1.01 1.96e-07 1.00
	640	2.45e-08 1.02 4.10e-08 1.02 9.68e-08 1.02
	1280 1.18e-08 1.05 1.99e-08 1.05 4.69e-08 1.05

Table 11 :

 11 Discharge errors and order of accuracy for the pair (c).

Table 12 :

 12 L 1 errors for the friction and topography 2D steady state using the implicit scheme.

	N	h		p		q
	900	4.01e-02	-	7.04e-02	-	7.02e-02	-
	3600	2.27e-02 0.82 3.89e-02 0.86 3.65e-02 0.94
	14400 1.28e-02 0.83 2.00e-02 0.96 1.85e-02 0.98
	57600 7.03e-03 0.86 1.01e-02 0.99 9.22e-03 1.00

Table 13 :

 13 L ∞ errors for the friction and topography 2D steady state using the implicit scheme.

	N	h		p		q
	900	6.89e-04	-	1.43e-03	-	1.04e-03	-
	3600	1.91e-04 1.85 3.90e-04 1.87 2.95e-04 1.81
	14400 5.11e-05 1.90 1.03e-04 1.92 8.05e-05 1.87
	57600 1.33e-05 1.94 2.67e-05 1.95 2.12e-05 1.92

Table 14 :

 14 L 1 errors for the friction and topography 2D steady state using the MUSCL scheme.

	N	h		p		q
	900	2.38e-03	-	4.13e-03	-	2.23e-03	-
	3600	8.05e-04 1.56 1.14e-03 1.86 7.58e-04 1.56
	14400 2.49e-04 1.69 3.04e-04 1.91 2.44e-04 1.63
	57600 7.36e-05 1.76 7.89e-05 1.94 7.25e-05 1.75

Table 15 :

 15 L ∞ errors for the friction and topography 2D steady state using the MUSCL scheme.

Table 16 :

 16 Height at the positions of the two vortices and position of the left shock wave. We have set h T = h(x T , y T , t end ) and h B = h(x B , y B , t end ).

		h T	h B	x S
	implicit, k = 0	2.81	2.81	59.3
	implicit, k = 2	6.08	6.08	53.9
	MUSCL, k = 0	2.33	2.30	59.8
	MUSCL, k = 2	5.86	5.86	52.5

Acknowledgments

C. Berthon, F. Foucher and V. Michel-Dansac would like to thank the ANR-12-IS01-0004-01 GEONUM for financial support. S. Clain acknowledges the Fundação para a Ciência e a Tecnologia for the funding of projects FCT-ANR/MAT-NAN/0122/2012 and UID/MAT/00013/2013.

The authors acknowledge the anonymous referees for very relevant remarks which greatly improved the scheme.

Note that we are dealing with a lake at rest configuration in the regions far from the waves. Indeed, the free surface unaffected by the rarefaction wave or the shock wave remains unperturbed. This means that h(x) + Z(x) = 2 and q(x) = 0 for x before the head of the rarefaction wave, and h(x) + Z(x) = 1 and q(x) = 0 for all x after the shock wave.

Such a property is satisfied by the explicit and implicit schemes. The MUSCL scheme, equipped with the well-balance correction (6.5), also exactly preserves this configuration. In addition, the MUSCL scheme provides a better approximation of the rarefaction wave and the shock wave.

Therefore, such an experiment highlights the interest of a using a wellbalanced scheme for such simulations, even if the whole domain does not involve a steady state. Moreover, it showcases the relevance of the convex combination procedure involved in the MUSCL scheme, to preserve steady states while still ensuring a better approximation of the exact solution.

In order to emphasize even more the importance of using a well-balanced scheme in such situations, we compare the implicit scheme with a non-wellbalanced one. The non-well-balanced scheme we use is the HLL one (see [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]), with a classical splitting treatment of the topography and friction source terms (see [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] for instance). The discharge obtained by both these schemes is compared on Figure 13 with the discharge of the reference solution.

Figure 13: Discharge for the wet dam-break experiment. Results of the implicit scheme and the HLL scheme with 100 discretization cells, compared with a reference solution computed with the MUSCL scheme using 10000 discretization cells.

On Figure 13, we observe that the implicit scheme exactly preserves the lake at rest configurations ahead of the waves. Indeed, the discharge vanishes in those areas. However, the HLL scheme does not preserve this steady state, and the discharge shows perturbations near the boundaries of the domain. In the explicit scheme is not suitable to perform the simulation of a dry area using the same CFL condition (3.5) as the implicit scheme. In the remainder of this section, the implicit scheme will be used for numerical experiments as soon they involve dry areas.

Dry dam-break over a non-flat topography

we consider a dry dam-break. The topography is still given by Z(x) = cos 2 (2πx) /2 and the Manning coefficient is now k = 10. The initial data is h L (x) + Z(x) = 2 and h R (x) = 0. The space domain is [0, 1] and the dam is located at x D = 0.5. Moreover, the boundaries at x = 0 and x = 1 are equipped with homogeneous Neumann boundary conditions. We use 100 discretization cells and display the results on Figure 15 at time t end = 0.05s. We choose the constants C = 7.5, m = 10 -6 and M = 10 -1 . We use the implicit and MUSCL schemes. Indeed, the explicit scheme present instabilities when dealing with dry/wet transitions. We also provide a reference solution for comparison purposes, obtained with the HR scheme using 10 4 discretization cells. The reference solution is once again correctly approximated by the two schemes, and the MUSCL scheme provides a better approximation of the solution than the implicit scheme. We reiterate the same remark concerning the exact preservation of the lake at rest before the head of the rarefaction wave.

Dry dam-break with two bumps

This last dry dam-break experiment concerns a more complicated topography, with two bumps. The space domain is [0, 5] and we choose to use 10 4 discretization cells with the MUSCL scheme to have a relevant simulation. The Such a 2D scheme will be able to exactly preserve the 1D steady states taken in the x direction or the y direction. Moreover, the splitting procedure allows simulating dry/wet interfaces. We use this 2D scheme for the simulations of three relevant situations. First, we focus on assessing the order of accuracy of the implicit and MUSCL schemes. Then, we turn to the simulation of a dry dam-break with two bumps, which is a 2D version of the experiment sented Section 7.3.4. The next experiment is a partial dam-break, and the last experiment represents an urban topography.

Order of accuracy assessment

The purpose of this first experiment is to assess the order of accuracy of the proposed schemes. Namely, the orders of accuracy of the implicit and the MUSCL schemes will be computed. The implicit scheme should be first-order accurate, while the MUSCL scheme should be second-order accurate.

To address this issue, we need a smooth exact solution of the 2D shallowwater equations with topography and friction (7.7). We elect to use a 2D steady state solution, which will not be exactly preserved by the well-balanced scheme, and that will allow computing the order of accuracy.

On the space domain [-0.3, 0.3] × [0.4, 1], we consider the following topography function:

where r = t (x, y). Then, to define a steady solution of (7.7), the water height and discharge are given by: h(x, y, t) = 1, p(x, y, t) = x r 2 , q(x, y, t) = y r 2 .

For the experiment, we take k = 10, and we prescribe the exact solution as both the initial and the boundary conditions. The final physical time is t end = 0.1s, and we set C = +∞, m = 0.05 and M = 1. The parameter C can here be infinite since the exact solution is very smooth and the topography is non-flat, thus making the cutoff procedure irrelevant. In addition, the other two parameters are chosen to ensure that the MUSCL scheme is second-order accurate everywhere.

Recall that the order of accuracy is defined, in the one-dimensional case, by (7.6). In our 2D case, the order of accuracy ρ is computed as follows:

Equipped with this definition, we compute the L 1 and L ∞ errors and the associated orders of accuracy of both schemes. These results are collected in Table 12 and Table 13 for the implicit scheme, and in Table 14 and Table 15 for the MUSCL scheme.

From Table 12, Table 13, Table 14 and Table 15, we conclude that both schemes have the expected order of accuracy. Indeed, the implicit scheme is first-order accurate, while the MUSCL scheme is second-order accurate. and the topography function is given by

The space domain is [0, 5]×[0, 1]. The initial discharge is zero in both directions, i.e. p(x, y, 0) = q(x, y, 0) = 0, and the initial water height is given by h(x, y, 0) = 6 if x < 0.7, 0 otherwise.

In addition, we prescribe wall boundary conditions, that is to say we impose p(0, y, t) = p(5, y, t) = 0 and q(x, 0, t) = q(x, 1, t) = 0 for all x, y and t, and homogeneous Neumann boundary conditions for the other variables.

The simulation runs until a physical time t end = 1.35s with the MUSCL scheme, using C = 1, m = 10 -5 and M = 25. The choice of these constants is motivated by the same reasons as invoked in Section 7.3.4. We take 288000 discretization cells (1200 in the x direction and 240 in the y direction). The results are presented on Figure 17 and Figure 18.

This experiment has been carried out to make sure that the numerical scheme still behaves correctly in two space dimensions and in the presence of dry/ wet transition. We recover a numerical solution involving the friction source term, which can be compared to the numerical solution without friction presented in [START_REF] Berthon | Efficient well-balanced hydrostatic upwind schemes for[END_REF]. In addition, this 2D experiment is similar to the 1D double bump experiment we presented on Figure 16. Indeed, the behavior of the water before it comes into contact with the first bump should be the same in both experiments. As expected, we obtain similar results in 1D and 2D.

Partial dam-break

Next, we focus on a partial dam-break (see for instance [START_REF] Nikolos | An unstructured node-centered finite volume scheme for shallow water flows with wet-dry fronts over complex topography[END_REF][START_REF] Clain | The MOOD method for the non-conservative shallow-water system[END_REF]). This experiment concerns a dam that has partially broken, leaving a corridor where the water flows. We consider the space domain [-100, 100] × [-100, 100]. The dam is located at the middle of the domain, in the x direction, and is 10m wide and 100m high. It breaks in the middle, over a length of 80m. Thus, we take the following topography function: