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

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

A non-negativity preserving and well-balanced scheme that exactly preserves all the smooth steady states of the shallow water system, including the moving ones, is proposed. In addition, the scheme must deal with vanishing water heights and transitions between wet and dry areas. A Godunov-type method is derived by using a relevant average of the source terms within the scheme, in order to enforce the required well-balance property. A second-order well-balanced MUSCL extension is also designed. Numerical experiments are carried out to check the properties of the scheme and assess the ability to exactly preserve all the steady states.

Introduction

During the last two decades, numerous schemes have been derived to preserve exactly (or, at least, accurately) the lake at rest. For instance, we refer the reader to the previous work by Bermudez and Vasquez [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF], and next to Greenberg and LeRoux [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF]. These works introduced the definition and the relevance of the well-balanced procedure. Such approaches were extended by Gosse [START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF] for nonlinear systems, by involving a nonlinear equation to be solved. More recently, in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF], the authors proposed a simplification (by enforcing vanishing velocities) of Gosse's work [START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF], yielding the so-called hydrostatic reconstruction (see [START_REF] Berthon | Efficient well-balanced hydrostatic upwind schemes for shallow-water equations[END_REF][START_REF] Castro | Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique[END_REF][START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction[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 geometrical source terms[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] for related work).

The steady states for the shallow-water equations with nonzero discharge are known to be more difficult to exactly capture than the lake at rest configuration. The critical role played by these specific solutions was illustrated in [START_REF] Goutal | Proceedings of the 2 nd Workshop on Dam-Break Wave Simulation[END_REF], were several benchmarks were exhibited. Next, in [START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF], a pioneer fully well-balanced scheme was designed to deal with these sensitive steady states, where the numerical technique is based on a suitable resolution of the Bernoulli equation. Next, several methods preserving the moving steady states were designed by involving high-order accurate techniques (see [START_REF] Castro Díaz | High order exactly well-balanced numerical methods for shallow water systems[END_REF][START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[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] for high-order and exactly well-balanced schemes, and [START_REF] Russo | High order well balanced schemes for systems of balance laws[END_REF] for a high-order accurate scheme on all steady state configurations).

More recently, 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], the authors have proposed an extension of the work by Gosse [START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF] in order to deal with Godunov-type schemes. Such Godunovtype schemes (see [START_REF] Gallice | Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source[END_REF][START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF]) are based on approximate Riemann solvers, whose intermediate states are obtained by solving a Bernoulli-type equation. This process allows the authors to get a fully well-balanced scheme preserving the entropy stability. Since the resolution of the Bernoulli-type equation has a large computational cost, we adopt in the present paper a linear formulation to deal with a general form of well-balanced states. In addition, a well-balanced HLLC scheme (see [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF]) is designed in [START_REF] Castro Díaz | A HLLC scheme for nonconservative hyperbolic problems. Application to turbidity currents with sediment transport[END_REF] for turbidity currents with sediment transport and in [START_REF] Sánchez-Linares | A HLLC scheme for Ripa model[END_REF] for the Ripa model. Both models are close to the shallowwater equations, and numerical methods can be easily transposed to the shallowwater case. These schemes are similar to the one presented here, since the intermediate states involve the source terms in both cases. Nevertheless, the main difference between the present work and the cited articles is the use of a HLLC-type approximate Riemann solver in [START_REF] Castro Díaz | A HLLC scheme for nonconservative hyperbolic problems. Application to turbidity currents with sediment transport[END_REF][START_REF] Sánchez-Linares | A HLLC scheme for Ripa model[END_REF], instead of the HLL-type solver we develop here. A HLLC-type solver includes more waves than a HLLtype solver, and therefore more unknowns to determine.

In the present paper, we propose a generic approach to provide a wellbalanced scheme. As a consequence, the objectives of the paper are to derive a numerical scheme to approximate the solutions of the shallow-water equations, and that satisfies the following properties:

1. exact preservation of all the smooth steady states for the system with topography;

2. non-negativity preservation for the water height under the usual CFL condition;

3. ability to handle the transitions between wet and dry areas.

The paper is organized as follows. First, we devote Section 2 to the study of smooth steady states for the shallow-water equations with topography. Some comments are also given to determine steady state solutions with a dry/wet transition. Afterwards, Section 3 is dedicated to the construction of a Godunovtype scheme. We then propose a general procedure to obtain a well-balanced scheme for the shallow-water equations with a source term on the discharge equation. At this stage, the definition of the source terms is not specified, and the resulting scheme will depend on the source term discretization. In other words, the well-balance property is obtained according to the PDE governing the steady states. Next, in Section 4, we consider the particular case of the topography source term, of which we propose a suitable discretization in order to exactly preserve the steady states governed by the topography. An extension of the scheme to dry/wet transitions is also designed. A second-order extension with a MUSCL technique is then studied in Section 5, and we conclude the document with numerical experiments in Section 6. Section 7 ends the study by a brief conclusion.

Steady states for the shallow-water equations with topography

The shallow-water model

This paper is devoted to designing a numerical scheme to approximate the solutions of the shallow-water equations with topography. The model of interest is governed by the following system:

   ∂ t h + ∂ x q = 0, ∂ t q + ∂ x q 2 h + 1 2 gh 2 = -gh∂ x Z. (2.1)
Equations (2.1) describe the behavior of water in a one-dimensional channel with a non-flat bottom. The modeled quantities are the water height h(x, t) ≥ 0 and its depth-averaged discharge q(x, t). The depth-averaged velocity u of the water is such that q = hu. The constant g > 0 stands for the gravity, while the function x → Z(x) is the topography. 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.

For the sake of simplicity in the notations, we rewrite (2.1) under the following condensed form:

∂ t W + ∂ x F (W ) = s(W, Z), W ∈ Ω, (2.2)
where

W = h q , F (W ) =   q q 2 h + 1 2 gh 2   , s(W, Z) = 0 -gh∂ x Z .
The homogeneous system deriving from canceling the source terms in (2.1) turns out to be hyperbolic with characteristic velocities 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 by c = gh.

(2.3)

Smooth steady states with positive water heights

In the present paper, we focus on the smooth steady state solutions of (2.2), which thus satisfy ∂ t W = 0. From now on, when discussing steady states, we assume a smooth topography, that is to say continuous and differentiable in space. Such steady states are governed by the following set of partial differential equations:

   ∂ x q = 0, ∂ x q 2 h + 1 2 gh 2 = -gh∂ x Z. (2.4)
The first equation immediately imposes a constant discharge q = q 0 . The second equation of (2.4) then becomes:

∂ x q 2 0 h + 1 2 gh 2 = -gh∂ x Z. (2.5)
This equation is a nonlinear ordinary differential equation, with unknown h. From (2.5), we can extract solutions in specific cases. Some of them have been extensively studied. For instance, by adopting q 0 = 0 and assuming h > 0, we recover the notable lake at rest steady state (see for instance [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 shallow-water equations[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]), defined by q = 0, h + Z = cst .

(

We now turn to the study of the equation (2.5) for q 0 = 0. Since functions are smooth, we rewrite (2.5) under the following algebraic form:

∂ x q 2 0 2h 2 + g (h + Z) = 0. (2.7)
Note that (2.7) is nothing but a statement of Bernoulli's principle. Now, let us consider a fixed x 0 ∈ R. We introduce h(x 0 ) = h 0 and Z(x 0 ) = Z 0 . For all x ∈ R, integration (2.7) over (x 0 , x) provides

q 2 0 2h 2 + g (h + Z) - q 2 0 2h 2 0 -g (h 0 + Z 0 ) = 0. (2.8)
To shorten the notations, we set

ξ(h; Z, h 0 , q 0 , Z 0 ) = q 2 0 2h 2 + g (h + Z) - q 2 0 2h 2 0 -g (h 0 + Z 0 ) . such that (2.8) rewrites ξ(h; Z, h 0 , q 0 , Z 0 ) = 0. (2.9)
Here, Z is effectively a parameter of the function ξ since the topography is given. Now, we study h = h(x), solution of (2.7), or equivalently (2.9), with initial condition h(x 0 ) = h 0 . In order to exhibit the solutions of (2.9), we first differentiate ξ with respect to h:

∂ξ ∂h (h; Z, h 0 , q 0 , Z 0 ) = - q 2 0 h 3 + g.
If h > h c , the function h → ξ(h; Z, h 0 , q 0 , Z 0 ) is strictly increasing, while if h < h c , h → ξ(h; Z, h 0 , q 0 , Z 0 ) is a strictly decreasing function, with h c such that ∂ξ ∂h (h c ) = 0, given by

h c = q 2 0 g 1 3
.

(2.10)

To define solutions h ∈ R * + of (2.9), we evaluate the sign of ξ. After straightforward computations, this function is proven to satisfy the following limits:

• lim h→+∞ ξ(h; Z, h 0 , q 0 , Z 0 ) = +∞,
• lim h→0 + ξ(h; Z, h 0 , q 0 , Z 0 ) = +∞. Moreover, the function ξ is immediately shown to verify the following evaluation for h = h c :

ξ(h c ; Z, h 0 , q 0 , Z 0 ) = q 2 0 2 3 h 2 c - 1 h 2 0 + g (Z -Z 0 -h 0 ) .
Figure 1 displays the function ξ(h; 0.75, 1, √ g, 1) with respect to h. Note that, with this particular choice of parameters, we have h c = 1. Moreover, in this particular case, we have ξ(h c ; 0.75, 1, √ g, 1) < 0, and the function admits two distinct roots. Equipped with these properties of ξ, we have the following statement.

Lemma 1. Assume h > 0 and q 0 = 0. Thus, h c > 0 according to (2.10), and the following properties hold:

(i) If ξ(h c ; Z, h 0 , q 0 , Z 0 ) > 0, then there is no solution to the equation (2.9).

(ii) If ξ(h c ; Z, h 0 , q 0 , Z 0 ) = 0, then the equation (2.9) admits a unique solution. This solution, h = h c , is a double root of the function h → ξ(h; Z, h 0 , q 0 , Z 0 ).

(iii) If ξ(h c ; Z, h 0 , q 0 , Z 0 ) < 0, then the equation (2.9) admits two distinct solutions, h sup ∈ (0, h c ) and h sub ∈ (h c , +∞).

Proof. The proof of this lemma relies on using all the properties of ξ we have obtained above. We recall that ξ tends to infinity as h tends to 0 or infinity. Moreover, the function ξ admits a unique minimum, reached for h = h c . Equipped with these properties, the proofs of (i), (ii) and (iii) are obvious. The proof is thus achieved.

Remark 1. Defining the Froude number as Fr = u/c, we have for steady states Fr = q 0 / gh 3 . Thus, h sub > h c corresponds to the subcritical case (Fr < 1) while h sup < h c is to the supercritical solution (Fr > 1).

Remark 2. Assume that, for all x ∈ R, h(x) = h c . From (2.5), the derivative of h with respect to x writes

∂ x h = h 3 ∂ x Z h 3 c -h 3 .
Therefore, if the solution is subcritical, i.e. h(x) > h c , then the sign of ∂ x h is the opposite of the sign of ∂ x Z, whereas, the sign of ∂ x h is that of ∂ x Z if the solution is supercritical. These results are in accordance with the subcritical and supercritical experiments presented in [START_REF] Goutal | Proceedings of the 2 nd Workshop on Dam-Break Wave Simulation[END_REF] for instance.

As a consequence, we have obtained the general form of steady states for the shallow-water equations (2.1) assuming q 0 = 0 and h > 0.

Smooth steady states with dry areas

Now, we turn to the study of steady states involving dry areas, that is to say areas where the water height is zero. Indeed, it is not straightforward that the condition h = 0 on a subset of the domain implies that q = 0, since the velocity may be unbounded. The goal of this subsection is to obtain a characterization of a steady state where such an area is present.

We begin by defining the kinetic energy in a wet area where h > 0, as follows:

E = 1 2 q 2 h .
Since we inject a bounded quantity of energy at the initial time, the kinetic energy has to be bounded, i.e.

E ∞ < +∞.

The boundedness of the kinetic energy allows us to state the following result.

Proposition 2. As soon as a dry area is involved, smooth steady states must be at rest.

Proof. Since E ∞ < +∞, we have necessarily q = O( √ h) when h tends to 0 + . Thus, we immediately obtain that, if there is some x D ∈ R such that h(x D ) = 0, then q(x D ) = 0. Now, recall from (2.4) that, for a steady state, ∂ x q = 0. Therefore, for all x ∈ R, q(x) = q(x D ) = 0, and we have u = 0 when h > 0: thus, the water is at rest. We conclude that the considered smooth steady state, involving a dry area, must be at rest.

A generic well-balanced scheme

We propose a generic strategy to derive numerical schemes that are able to exactly capture all the steady state solutions. To address such an issue, we first briefly recall some important steps of the construction of a Godunov-type scheme based on a two-state 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]). Then, the scheme is extended to the shallow-water equations with a general source term in order to preserve the steady states.

Godunov-type schemes

We first introduce the discretization of space and time. 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 ), of volume ∆x and centered at x i = x i-1 2 + ∆x 2 , for all i ∈ Z. Now, we assume known a piecewise constant approximation of W (x, t) at time t n , denoted by W ∆ (x, t n ) and defined, for all x ∈ (x i-1 2 , x i+ 1 2 ), by W ∆ (x, t n ) = W n i . In order to evolve this approximation in time, we suggest a Godunov-type finite volume method. The scheme is obtained by setting an approximate Riemann solver, denoted by W ( x t ; W L , W R ), at each interface x i+ 1 2 . This solver is designed in order to mimic the exact solution of a Riemann problem for the initial system (2.1).

From now on, let us consider an approximate Riemann solver made of four constant states, separated by three discontinuities. We define it as follows (see Figure 2):

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.1)
where λ L and λ R denote some characteristic velocities, and W * L and W * R stand for the intermediate states, to be detailed later. To enforce λ L < 0 and λ R > 0, we choose the following expressions of the characteristic velocities (see for instance [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF] and references therein), where c denotes the sound speed, defined by (2.3):

λ L = min (-|u L | -c L , -|u R | -c R , -ε λ ) , λ R = max (|u L | + c L , |u R | + c R , ε λ ) , (3.2)
with ε λ to be fixed in the numerical applications. Remark 3. Note that λ L < 0 and λ R > 0 even for the supercritical case, where both velocities should have the same sign. Because of this choice, both intermediate states W * L and W * R are always considered. This is a very important ingredient in the construction of the well-balanced approximate Riemann solver.

W L W R λ L λ R W * L W * R x 0 t
Equipped with W , we determine the updated states W n+1 i . Indeed, let us introduce the juxtaposition of the approximate Riemann solvers stated at each interface (see Figure 3). Such a juxtaposition is denoted by W ∆ (x, t n + t), for all t ∈ (0, ∆t]. As a consequence, we have

∀x ∈ [x i , x i+1 ) , ∀t ∈ (0, ∆t] , W ∆ (x, t n + t) = W x -x i+ 1 2 t ; W n i , W n i+1 . (3.3
) Moreover, to prevent the waves from interacting, we impose the following CFLlike condition:

∆t ≤ ∆x 2Λ , where Λ = max i∈Z -λ L i+ 1 2 , λ R i+ 1 2 . (3.4) 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 the prescribed approximate Riemann solver.

We now define W n+1 i as the average of W ∆ (x, t n +∆t) 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 + ∆t)dx.
Arguing the definition (3.3) of W ∆ , we get:

W n+1 i = 1 ∆x 0 -∆x/2 W x ∆t ; W n i-1 , W n i dx + 1 ∆x ∆x/2 0 W x ∆t ; W n i , W n i+1 dx. (3.5)
Finally, from the definition (3.1) of the approximate Riemann solver W , we have, after straightforward computations and adopting clear notations, the following expression of W n+1 i : to ensure that the scheme is consistent and preserves all the steady states.

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

A well-balanced approximate Riemann solver for general source terms

We are now interested in deriving generic well-balanced schemes, that is to say schemes that preserve the steady states. Concerning the source term, at this level, we do not specify the definition of s(W ) = t (0, S(W )). Thus, S(W ) may represent the contribution of the topography and/or any other source term on the discharge equation. Note that s may depend on other quantities than W . For instance, in the case of the topography source term, s also depends on the topography function Z. As a consequence, this source term is denoted by s(W, Z) = t (0, S(W, Z)) from now on.

The well-balanced scheme is obtained as soon as relevant definitions of the intermediate states

W * L = t (h * L , q * L ) and W * R = t (h * R , q * R ) are computed. Evalu- ation of W *
L and W * R requires four equations to compute the four unknowns h * L , h * R , q * L and q * R . To address such an issue, we look for intermediate states that enforce both consistency and well-balancedness.

The generic non-conservative shallow-water system is given by

   ∂ t h + ∂ x q = 0, ∂ t q + ∂ x q 2 h + 1 2 gh 2 = S(W, Z). (3.7)
Here, the steady states are governed by

     q = q 0 , ∂ x q 2 0 h + 1 2 gh 2 = S(W, Z). (3.8)
We note that the ordinary differential equation (3.8) cannot be rewritten under an algebraic form for any generic source term S(W, Z). However, in the specific case of the topography source term, we can rewrite (3.8) as an algebraic relation. For instance, the lake at rest is given by (2.6), and the moving steady states are given by (2.7). Such an algebraic relation will allow us to obtain a relevant numerical approximation of the source term.

Remark 4. When no algebraic relation is available for a source term S(W, Z), we can still derive a suitable approximation by splitting the source term into a sum of source terms for which we have algebraic expressions.

From (3.6), we deduce that the solution is stationary, i.e.

W n+1 i = W n i for all i ∈ Z, if we have W * L = W L and W * R = W R .

Consistency

We introduce a necessary condition for the scheme to be consistent with (3.7). We denote by W R x t ; W L , W R the exact solution of the Riemann problem for (3.7). From [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF], the consistency condition states that the average of the approximate Riemann solver W , defined by (3.1), over a cell has to be equal to the average of the exact solution W R over the same cell. As a consequence, we have to impose the following equality:

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)
On the one hand, we rewrite the right-hand side of (3.9) as:

1 ∆x ∆x/2 -∆x/2 W R x ∆t ; W L , W R dx = 1 2 (W L + W R ) - ∆t ∆x (F (W R ) -F (W L )) + 1 ∆x ∆x/2 -∆x/2 ∆t 0 s W R x t ; W L , W R , Z(x) dt dx.
(3.10) On the other hand, a direct computation of the integral in the left-hand side of (3.9) yields:

1 ∆x ∆x/2 -∆x/2 W x ∆t ; W L , W R dx = 1 2 (W L + W R ) + λ L ∆t ∆x (W L -W * L ) -λ R ∆t ∆x (W R -W * R ). (3.11)
Recall that the first component of s(W, Z) is 0. Therefore, combining (3.10) with (3.11), within (3.9), leads us to the following relations:

λ R h * R -λ L h * L = λ R h R -λ L h L -[q], λ R q * R -λ L q * L = λ R q R -λ L q L - q 2 h + 1 2 gh 2 + 1 ∆t ∆x/2 -∆x/2 ∆t 0 S W R x t ; W L , W R , Z(x) dt dx, (3.12)
where

[X] = X R -X L .
For the sake of simplicity, let us now introduce the following notations:

(λ R -λ L )h HLL = λ R h R -λ L h L -[q], (3.13a)
(λ R -λ L )q HLL = λ R q R -λ L q L - q 2 h + 1 2 gh 2 . (3.13b)
Let us emphasize that h HLL and q HLL coincide with the well-known intermediate state introduced by Harten, Lax and van Leer in [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]. We remark that h HLL > 0 for large enough -λ L and λ R . Thanks to these notations, (3.12) can be rewritten as follows:

λ R h * R -λ L h * L = (λ R -λ L )h HLL , (3.14a)
λ R q * R -λ L q * L = (λ R -λ L )q HLL + 1 ∆t ∆x/2 -∆x/2 ∆t 0 S W R x t ; W L , W R , Z(x) dt dx. (3.14b)
After [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], we now state that q * L = q * R , and we denote this value by q * . We will see later on that this choice allows the recovery of the required well-balance property. Using q * , we rewrite (3.14b) under the form

q * = q HLL + 1 λ R -λ L 1 ∆t ∆x/2 -∆x/2 ∆t 0 S W R x t ; W L , W R , Z(x) dt dx.
(3.15)

Well-balancedness parametrization

By adopting a non-explicit source term, given by S(W, Z), we have to introduce two parameters, denoted S and q. These parameters are such that the solution characterized by h * L , h * R and q * corresponds to a well-balanced one according to (3.8), and a consistent one according to (3.14a) and (3.15). In the next section, the source term will be chosen as the topography, and we will impose the additional suitable algebraic relation to govern the steady states and thus to find both additional parameters.

We begin by introducing a first parameter S. It is a consistent approximation of the mean value of the source term S(W, Z), and it implicitly depends on W L and W R , as well as Z L and Z R . The approximation of the source term average is thus given by:

1 ∆x 1 ∆t ∆x/2 -∆x/2 ∆t 0 S(W, Z) dt dx S. (3.16)
We now plug S into (3.15). As a consequence, we impose that the three unknowns h * L , h * R and q * satisfy the following two relations:

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

q * = q HLL + S∆x λ R -λ L . (3.17b)
A relevant definition of S is required in order to fully determine q * . Suitable expressions will be given in the next section for the specific case of varying topography. Before completing the system (3.17) to fully determine the intermediate water heights h * L and h * R , let us assume that h L = 0 and h R = 0. We need to specify discrete steady states associated to (3.8). Indeed, the left and right states, W L and W R , define a steady state if the following relations hold:

   q L = q R = q 0 , (3.18a)
q 2 0 1 h + g 2 h 2 = S∆x. (3.18b)
Such equations are nothing but a discretization of the steady relation (3.8). We rewrite (3.18b) as follows:

α(h R -h L ) = S ∆x, where α = -q 2 0 h L h R + g 2 (h L + h R ).
Now, to complete the determination of the intermediate water heights, we adopt an extension of the above relation, given by:

α(h * R -h * L ) = S ∆x, where α = -q 2 h L h R + g 2 (h L + h R ), (3.19)
where we have introduced a parameter q(W L , W R ) which must be such that q = q 0 as soon as W L and W R define a steady state, i.e. as soon as the relations (3.18) are verified. Therefore, (3.17a) and (3.19) define the following linear system:

λ R h * R -λ L h * L = (λ R -λ L )h HLL , α(h * R -h * L ) = S ∆x, from which we deduce h * L = h HLL - λ R S ∆x α(λ R -λ L ) , (3.20a)
h * R = h HLL - λ L S ∆x α(λ R -λ L )
.

(3.20b)

Properties of the approximate Riemann solver

The goal of this subsection is to assert the properties verified by the approximate Riemann solver. The intermediate states have been chosen to yield a consistent and well-balanced approximate Riemann solver. In addition, it is possible to preserve the positivity of the water height while retaining the consistency and well-balancedness of the scheme. Indeed, let us emphasize that the definitions (3.20) do not ensure positive intermediate water heights. To address such an issue, we follow the procedure proposed in [START_REF] Audusse | A simple well-balanced and positive numerical scheme for the shallow-water system[END_REF] (see also [START_REF] Berthon | A Well-Balanced Finite Volume Scheme for a Mixed Hyperbolic/Parabolic System to Model Chemotaxis[END_REF]) and presented in Figure 4. It consists in enforcing the positivity of h * L and h * R , while still ensuring that they satisfy the consistency relation (3.17a). To that end, we introduce the parameter ε, given by

ε = min (h L , h R , h HLL ) . (3.21) Note that h L > 0, h R > 0 and h HLL > 0: hence ε > 0. If h * L < ε, we take h * L = ε, and h * R is chosen according to (3.17a), which guarantees h * R > 0 (see Figure 4). A similar procedure is applied if h * R < ε. After the correction procedure, we have h * L ≥ ε and h * R ≥ ε. h * L h * R 0 ε ε λ R h * R = λ L h * L + (λ R -λ L )h HLL 2 1
(h * L , h * R ) belongs to domain 1 , then h * L and h * R are not modified. However, if (h * L , h * R ) corresponds to a point within domain 2 , we replace (h * L , h * R ) with (ε, (1 -λ L λ R )h HLL + λ L λ R ε) according to (3.17a).
We now state the expressions of the intermediate states, for given W L , W R , S and q, with the positivity correction. The parameters S and q will be defined in the next section for a specific source term. The intermediate states of the approximate Riemann solver are then given by W

* L = t (h * L , q * L ) and W * R = t (h * R , q * R ), with q * L = q * R = q * = q HLL + S∆x λ R -λ L , (3.22a)
h * L = min max h HLL - λ R S ∆x α(λ R -λ L ) , ε , 1 - λ R λ L h HLL + λ R λ L ε , (3.22b)
h * R = min max h HLL - λ L S ∆x α(λ R -λ L ) , ε , 1 - λ L λ R h HLL + λ L λ R ε , (3.22c)
where α has been defined in (3.19) as follows:

α = -q 2 h L h R + g 2 (h L + h R ), (3.23)
and where the quantities h HLL and q HLL are defined by (3.13).

Lemma 3. Assume h L and h R to be positive. Then, the intermediate states

W * L = t (h * L , q * L ) and W * R = t (h * R , q * R )
given by (3.22) satisfy the following properties:

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

(ii) positivity preservation: h * L ≥ ε and h * R ≥ ε;

(iii) well-balancedness: W * L and W * R satisfy the property (WB). Proof. Since ε is defined by (3.21), we obviously get the required property (ii). Indeed, after (3.22), h * L and h * R stand for the minimum of quantities that are greater than or equal to ε.

Next, let us set

h * L = h HLL - λ R S ∆x α(λ R -λ L ) and h * R = h HLL - λ L S ∆x α(λ R -λ L )
.

We immediately get the following identity:

λ R h * R -λ L h * L = (λ R -λ L )h HLL ,
which means that the heights h * L and h * R satisfy the consistency relation (3.17a). Since (3.17b) is obviously verified by q * , the property (i) is established as soon as h * L and h * R are proven to satisfy (3.17a). Assume that |λ L | and |λ R | are chosen large enough to ensure h HLL > 0 with λ L < 0 < λ R . We have the following three configurations.

• If h * L ≥ ε and h * R ≥ ε, then the relations (3.22) give h * L = h * L and h * R = h * R . • If h * L < ε, then from (3.22) we get h * L = ε and h * R = 1 -λ L λ R h HLL + λ L λ R ε. • Similarly, if h * R < ε, then we have h * R = ε and h * L = 1 -λ R λ L h HLL + λ R λ L ε.
We note that the consistency relation (3.17a) systematically holds and property (i) is proven.

We now have to check that the positivity procedure, involving relations (3.22b) and (3.22c), does not interfere with the well-balance property. In order to prove the well-balancedness, we assume that W L and W R define a steady state, given by (3.18), and we show that, in this case,

W * L = W L and W * R = W R . If h * L = h * L and h * R = h * R
, the property (iii) holds by construction. Now, we complete the proof as the intermediate states involve the water heights h * L and h * R given by (3.22). From the definition (3.22a) of q * and the steady relations (3.18) satisfied by W L and W R , we deduce

q * = λ R q 0 -λ L q 0 λ R -λ L - 1 λ R -λ L q 2 0 h + 1 2 gh 2 + S∆x λ R -λ L = q 0 - 1 λ R -λ L q 2 0 1 h + g 2 h 2 -q 2 0 1 h - g 2 h 2 = q 0 .
We now have to prove that h * L = h L and h * R = h R . First, let us compute S ∆x α at the equilibrium using (3.18) and (3.23). One has:

S ∆x α = q 2 0 1 h + g 2 h 2 -q 2 0 h L h R + g 2 (h L + h R ) = [h].
We then compute h * L at the equilibrium. According to (3.20a), we have:

h * L = λ R h R -λ L h L λ R -λ L - [q] λ R -λ L - λ R S ∆x α(λ R -λ L ) = λ R h R -λ L h L -λ R h R + λ R h L λ R -λ L = h L . Similar computations with (3.20b) lead to h * R = h R . Moreover, from the defini- tion (3.21) of ε, we have h L ≥ ε and h R ≥ ε. Therefore h * L ≥ ε and h * R ≥ ε, and we have h * L = h * L = h L and h * R = h * R = h R .
This proves the (WB) property, that is to say W * L = W L and W * R = W R as soon as W L and W R define a steady state. This concludes the proof of the well-balancedness, and therefore Lemma 3 is established.

Equipped with the suitable intermediate states, we can assert the following result concerning the full scheme (3.6). 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.22a), (3.22b) and (3.22c), respectively. Also, assume that the source term approximation S is consistent with the source term S according to (3.16). Finally, assume that, as soon as (W n i ) i∈Z defines a steady state, S verifies (3.18b) and q = q 0 . Then the Godunov-type scheme, given by (3.6) under the CFL restriction (3.4), satisfies the following properties:

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

positivity preservation

: ∀i ∈ Z, W n+1 i ∈ Ω * ; 3. well-balancedness: if (W n i ) i∈Z defines a steady state, then ∀i ∈ Z, W n+1 i = W n i .
Proof. The consistency is a direct consequence of Lemma 3. We turn to proving that ∀i ∈ Z, h n+1 i > 0 as soon as

h n i > 0. From Lemma 3, we have h L, * i+ 1 2 ≥ ε n i and h R, * i+ 1 2 ≥ ε n i , where ε n i = min(h n i , h n i+1 , h HLL i+ 1 2
) and h HLL i+ 1 2 is given by evaluating (3.13a) between states W n i and W n i+1 . Since the scheme under consideration is given by (3.5), h n+1 i turns out to be the sum of positive quantities.

We finally need to prove the well-balancedness of the scheme. Once again, this property comes from Lemma 3. Indeed, let us consider that (W n i ) i∈Z defines a steady state. Therefore, ∀i ∈ Z, W n i and W n i+1 define a steady state, and Lemma 3 gives

W L, * i+ 1 2 = W n i and W R, * i+ 1 2 = W n i+1 , which in turn yields W n+1 i = W n
i for all i ∈ Z. This concludes the proof of Theorem 4.

Remark 5. Because of the arbitrary small parameter ε > 0, introduced in (3.22) to enforce the positivity of the intermediate water heights, the updated water height never vanishes. In a specific section, devoted to dry/wet transitions, we will present an extension of the scheme to deal with dry areas in the case where the source term is the topography. At this level, we reject vanishing water heights because of the definition of S as well as the term S/α involved within the definitions (3.22) of h * L and h * R . As soon as the full characterization of S is established, the scheme will be extended to allow ε = 0 in the definition (3.22).

Application to the topography source term

In this section, our concern is the derivation of both parameters S and q in order to complete the scheme (3.6). Here, we consider the topography source term, i.e.

S(W, Z) = S t (W, Z) = -gh∂ x Z.

In this specific case, we shall denote these parameters S t and qt .

Determination of the parameters

We begin by determining S t . Let us recall that the steady states are governed by the following PDE:

q 2 0 ∂ x 1 h + g 2 ∂ x h 2 = S t . (4.1)
Since we are only considering the topography source term and smooth steady states, (4.1) can be rewritten under the following algebraic form:

q 2 0 2 1 h 2 + g[h + Z] = 0. (4.2)
Consider W L and W R defining a steady state: the states W L and W R are thus given by (3.18), and they must satisfy the algebraic relation (4.2). As a consequence, W L and W R verify the following two identities:

q 2 0 1 h + g 2 h 2 = S t ∆x, (4.3a
)

q 2 0 2 1 h 2 + g[h + Z] = 0. (4.3b)
The first identity (4.3a) corresponds to the generic steady relation applied to the topography, while the second one (4.3b) is nothing but the algebraic relation specific to the topography source term. We exhibit the expression of S t from the above relations as follows. First, from (4.3b), we extract the expression of q 2 0 :

q 2 0 = 2g[h + Z] h 2 L h 2 R h 2 R -h 2 L ,
which is plugged into (4.3a). We then get the following definition of S t :

S t ∆x = g 2 h 2 -2g[h + Z] h L h R h L + h R ,
which can be rewritten as follows:

S t ∆x = -2g[Z] h L h R h L + h R + g 2 [h] 3 h L + h R . (4.4)
Let us emphasize that such a definition of the approximation of the topography source term can be found in the literature. For instance, the reader is referred to [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] (see also [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF] for related expressions).

An important ingredient in the consistency of the scheme is that the source term approximation S t has to be consistent with the actual source term -gh∂ x Z, assuming positive water heights. For instance, when the topography is flat, i.e.

[Z] = 0, the actual topography source term vanishes. Therefore, in order for S t to be consistent with the actual source term, we need S t = O(∆x) as soon as the topography is flat. However, as underlined 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] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF], S t is no longer consistent with zero when the topography is flat. Indeed, in this case, we have

S t = g 2(h L + h R ) [h] 3 ∆x . ( 4.5)
In order to recover the required consistency, i.e. S t = O(∆x) for a flat topography, we adopt the strategy proposed 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]. We modify S t as follows:

S t ∆x = -2g[Z] h L h R h L + h R + g 2 [h] 3 c h L + h R , (4.6)
where

[h] c is a cutoff of [h] = h R -h L ,
defined as follows, with C a positive constant that does not depend on ∆x:

[h] c = h R -h L if |h R -h L | ≤ C ∆x, sgn(h R -h L )C ∆x otherwise. ( 4.7)
This new expression of S t is consistent with the topography source term S t . Indeed, the cutoff procedure enforces |[h] c | ≤ C∆x, and therefore that S t = O(∆x 2 ) when S t is given by (4.5). In addition, this procedure ensures that the scheme is well-balanced according to smooth steady states. However, it also means that the source term approximation S t does not vanish when the topography is flat, and therefore that the scheme does not reduce to a conservative scheme in that case.

Remark 6. For a smooth water height h, the relation h R -h L = O(∆x) obviously holds. Therefore, for a smooth water height, there exists

K ∈ R * + such that |h R -h L | ≤ K∆x.
As a consequence, for a smooth water height, there exists C such that [h] c = [h], after (4.7). Indeed, taking C < K suffices. In this case, S t is given by (4.4).

Note that the expression of S t given by (4.6) has been obtained by considering W L and W R defining a steady state. Since this expression only depends on the left and right states, it is relevant to extend it to the case where these states do not define a steady state, and actually use this expression for all W L and W R .

To achieve the characterization of the scheme (3.6), we also need to define the parameter q introduced in (3. [START_REF] Clain | A high-order finite volume method for systems of conservation laws-Multi-dimensional Optimal Order Detection (MOOD)[END_REF]. In fact, we can choose any expression for q, as long as it ensures q = q 0 when W L and W R define a steady state. We decide to use the following expression for qt :

qt = q * . (4.8)
From (3.23), we deduce:

α t = -(q * ) 2 h L h R + g 2 (h L + h R ) . (4.9)

Extension to dry/wet transitions

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

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

We now turn to defining a new expression of S t for a vanishing h L or h R . Since the expression (4.6) relied on the assumption that both h L and h R were positive, we cannot use this expression in the present case. In order to determine a new formula for S t , as in the previous subsection, we begin by assuming that W L and W R define a steady state, with vanishing h L or h R . Therefore, from Proposition 2, we have q 0 = 0 as soon as W L and W R define a steady state. Thus, this steady state is a lake at rest steady state, i.e. [h+Z] = 0. In addition, the above assumption ensures that u L = u R = 0. We then use (4.3a) to obtain S t ∆x = g[h 2 ]/2. Now, plugging [h] = -[Z] into this equality, we get the new expression of S t ∆x, to be substituted to (4.6) as soon as h L or h R vanishes:

S t ∆x = -g (Z R -Z L ) h R + h L 2 . (4.10)
Now, recall the definition (4.9) of α t . Note that the quantity α t is ill-defined for h L = 0 or h R = 0. In order to determine a suitable expression of α t for h L = 0 or h R = 0, note that q * = q 0 = 0 as soon as W L and W R define a steady state. Therefore, u L = u R = 0. Thus, as soon as h L = 0 or h R = 0, α t is given by

α t = g 2 (h L + h R ). (4.11)
Finally, 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 enforcing, as soon as both h L and h R are zero, the following expressions: S t ∆x = 0 and S t ∆x α t = 0. (4.12)

We now regroup the three cases (4.6), (4.10) and (4.12), to get the following final expression of S t :

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.13)
with [h] c given by (4.7).

In the same three cases, we define S t ∆x/α t as follows, using (4.9), (4.11) and (4.12):

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.

(4.14)

Note that the expressions (4.13) and (4.14) are not continuous. However, the consistency property holds. Indeed, consider smooth functions h and Z. We fix h L = h(x) and h R = h(x + ∆x). In addition, we take Z L = Z(x) and Z R = Z(x + ∆x). Plugging this data in each of the three cases described by (4.13), i.e. for h L ≥ 0 and h R ≥ 0, we obtain S t = -gh∂ x Z + O(∆x). As a consequence, the consistency property holds. Therefore, the loss of continuity at the discrete level impacts the error term instead of the consistency with the source term at the continuous level.

Equipped with the new expressions (4.13) and (4.14) of S t ∆x and S t ∆x/α t , we can state the following result, concerning the approximate Riemann solver with non-negative water heights. Lemma 5. There exists C ∈ R * + such that the intermediate states (3.22), defined with qt , S t and S t ∆x/α t respectively given by (4.8), (4.13) and (4.14), satisfy the following properties:

(i) consistency with the shallow-water equations with topography (2.1);

(ii) the well-balancedness principle (WB), i.e. preservation of all continuous steady states governed by (4.3). Note that Remark 6 applies in this case.

Moreover, we have:

(iii) with ε > 0 given by (3.21), 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.

Proof. For ε > 0, the proofs of (i), (ii) and (iii) are immediate. They come from Lemma 3 as well as the results obtained in the current section. Now, assume ε = 0. The intermediate states then rewrite:

q * = q HLL + S t ∆x λ R -λ L , (4.15a)
h * L = min h HLL - λ R S t ∆x α t (λ R -λ L ) + , 1 - λ R λ L h HLL , (4.15b)
h * R = min h HLL - λ L S t ∆x α t (λ R -λ L ) + , 1 - λ L λ R h HLL , (4.15c)
with α t given by (4.14). From (4.13) and (4.14), the above expressions are well-defined for all h L ≥ 0 and h R ≥ 0. Moreover, these new intermediate states with ε = 0 can be easily shown to satisfy properties (i), (ii) and (iv).

First, the proof of the consistency stated Lemma 3 is preserved, and thus (i) holds. Then, assertion (iv) is a direct consequence of the fact that h * L and h * R are defined as the minimum of non-negative quantities. Finally, in order to prove (ii) with ε = 0, we assume that W L and W R define a steady state according to (4.3). The goal is now to prove that the (WB) principle holds, i.e. W * L = W L and W * R = W R . If h L > 0, h R > 0 and h HLL > 0, then we know from Lemma 3 that W * L = W L and W * R = W R . Now, we assume that h L = 0, and thus that q L = q R = 0. Recall that, when h L or h R vanishes and W L and W R define a steady state, we have [h + Z] = 0. Then, using (4.13), the definition (4.15a) immediately yields q * = 0 = q 0 . Moreover, from (4.15b) and (4.15c), using (4.14) yields:

h * L = λ R h R λ R -λ L + λ R (Z R -Z L ) λ R -λ L , h * R = min λ R h R λ R -λ L + λ L (Z R -Z L ) λ R -λ L , 1 - λ L λ R λ R h R λ R -λ L . (4.16)
Then, with [h + Z] = 0, the relations (4.16) immediately yield h * L = 0 = h L and h * R = h R . A similar sequence of arguments is used to prove the well-balancedness when h R vanishes, or when both h L and h R are zero.

Thus, the proof is achieved.

This lemma allows us to state the following result, that concerns the full scheme (3.6), providing an extension of Theorem 4 for vanishing h.

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.17)
where q * , h * L and h * R are given by (4.15). Then the Godunov-type scheme, given by (3.6) under the CFL restriction (3.4), satisfies the following properties:

1. consistency with the shallow-water system (2.1); 2. positivity preservation: ∀i ∈ Z, W n+1 i ∈ Ω;

3. well-balancedness: if (W n i ) i∈Z defines a continuous steady state according to (4.3), then ∀i ∈ Z, W n+1 i = W n i . Proof. The same arguments as used in the proof of Theorem 4, while using the results of Lemma 5, yield the proof.

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 for fluid dynamics[END_REF]). The purpose of this extension is to increase the space and time accuracy of the scheme by using a piecewise linear reconstruction instead of piecewise constant. In order to derive this secondorder scheme, we first rewrite the scheme (3.6) in order to exhibit the numerical flux function and the source term contribution. Then, we present the MUSCL technique itself. Finally, we introduce a novel convex combination between the first-order and second-order schemes in order to recover the well-balance property. Note that the scheme can be extended to higher order by using a suitable polynomial reconstruction instead of the linear reconstruction involved in the MUSCL technique.

Rewriting the scheme

The goal of this subsection is to rewrite the scheme (3.6), with intermediate states given by (4.15), under the form proposed above. After straightforward computations (see for instance [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]), the scheme (3.6) can be rewritten as

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 .
The quantities f n i+ 1 2 and s n i+ 1 2 denote numerical approximations of the flux and source term, respectively, at the interface x i+ 1 2 . They 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 .
(5.1)

Here, adopting clear notations, (S t ) n i+ 1 2 is given by:

(S t ) n i+ 1 2 = S t h n i , h n i+1 , Z i , Z i+1 , (5.2)
where S t is defined by (4.13). The scheme then reads as follows:

h n+1 i = h n i - ∆t ∆x (f h ) n i+ 1 2 -(f h ) n i-1 2 , 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 ,
where the approximate fluxes are defined by (5.1) and

f n i+ 1 2 = f (W n i , W n i+1 ) = 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

The MUSCL procedure

Equipped with the numerical flux function and the source term contribution, we can state the MUSCL procedure. Consider w ∈ {h, q, h + Z}. The reconstruction is carried out by 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 slope of the linear reconstruction, defined by

σ n i = minmod w n i+1 -w n i ∆x , w n i -w n i-1 ∆x .
We have applied a limiter to this slope to improve the stability of the scheme.

Here, we have chosen the 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 reconstruction of w n i at the interfaces within the cell (x i-1 2 , x i+ 1 2 ) is then 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 .
(5.

3)

The reconstructed value of Z at the interfaces is finally obtained from the values of h + Z and h. Therefore, the updated states are given by the following scheme:

W n+1 i = W n i - ∆t ∆x f (W + i , W - i+1 ) -f (W + i-1 , W - i ) + ∆t 2 s n i+ 1 2 + s n i-1 2 , (5.4)
where the source term contributions s n i+ 1 2 and s n i-1 2 use the following approximation of the source term instead of (5.2):

(S t ) n i+ 1 2 = S t h + i , h - i+1 , Z + i , Z - i+1 .
Finally, the scheme's time accuracy is improved by Heun's method (see [START_REF] Gottlieb | Strong stability-preserving highorder time discretization methods[END_REF]).

The well-balance property

Because of the reconstruction step (5.3), it is worth noting that the designed MUSCL scheme does not preserve all the steady states. We therefore 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][START_REF] Figueiredo | Second-order finite volume MOOD method for the shallow water with dry/wet interface[END_REF] for more recent applications) to restore this essential property. Consider that W n i and W n i+1 define a steady state. At this level, the MUSCL scheme provides a second-order approximation of this steady state, while the first-order wellbalanced scheme previously derived exactly captures such a state.

Therefore, 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 , (5.5)
where 0 ≤ θ n i ≤ 1 is the parameter of the convex combination. Note that the states are not reconstructed if θ n i = 0, while the full MUSCL scheme is obtained by taking θ n i = 1. The objective is now to propose a suitable process to define the parameter θ n i . In order to have a relevant definition of θ n i , 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
, where S t i+ 1 2 = S t (h n i , h n i+1 , Z i , Z i+1 ), with S t defined by (4.13). Note that this quantity turns out the be a residue that governs steady states, according to (4.3). As a consequence, ∆ψ n i+ 1 2 vanishes when W n i and W n i+1 define a steady state. Next, 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
, which vanishes when W n i-1 , W n i and W n i+1 define a steady state. Equipped with ϕ n i , the parameter of the convex combination θ n i is defined for some M > m > 0 as follows (see Figure 5):

θ 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.
(5.6) Such a definition ensures that the first-order well-balanced scheme is used if the equilibrium error ϕ n i is small enough. Moreover, the MUSCL scheme is used if the states are far from defining a steady state, i.e. ϕ n i is large enough. In addition, the closer the states are to the equilibrium, the more the convex combination will favor the first-order well-balanced scheme.

0 1 m∆x M ∆x θ n i ϕ n i

Numerical experiments

Numerical simulations are carried out to test the scheme devised in the previous sections. We first check that the scheme preserves lake at rest steady states. Then, we assess the well-balancedness of the scheme by studying moving steady states. Afterwards, two experiments designed to validate the scheme are performed. Finally, we present a two-dimensional (2D) experiment to provide an error analysis of the scheme.

We also compare the proposed scheme (3.6) -(4.17) (denoted WBt from now on) with two classical schemes: the HLL scheme (see [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]) and the hydrostatic reconstruction (HR) scheme (see [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]) applied to the HLL flux (see [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]).

Since the HLL scheme is designed for conservative systems, we take the topography contribution into account by using a splitting method. The purpose of carrying out simulations with the HLL scheme is to highlight that the wellbalance is an important property for a scheme to possess. In addition, the choice of the HLL scheme for comparisons is relevant since the construction of our scheme is based on a HLL-like construction.

For a domain discretized with N cells, we compute the L 1 , L 2 and L ∞ errors for a bounded function w using the following expressions:

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 |, (6.1
) where w i and w ex i are respectively the approximate and the exact solution at cell i and at the final physical time t end . We recall that the time step ∆t is given by the CFL condition (3.4):

∆t ≤ ∆x 2Λ , where Λ = max i∈Z -λ L i+ 1 2 , λ R i+ 1 2 .
The constants will be chosen according to Table 1.

h q L 1 L 2 L ∞ L 1 L 2 L ∞
HLL 5.09e-05 4.36e-04 5.04e-03 6.72e-03 9.56e-03 1.56e-02 HR 0 0 0 0 0 0 WBt 1.11e-18 1.11e-17 1.11e-16 0 0 0 h q

L 1 L 2 L ∞ L 1 L 2 L ∞
HLL 2.47e-03 1.30e-02 1.33e-01 6.46e-03 3.48e-02 3.70e-01 HR 0 0 0 0 0 0 WBt 0 0 0 0 0 0 2x-0.5)1 [0.5,1] (x). The initial data are h(x) = (1-Z(x)) + and q(x) = 0. The free surface and topography are displayed Figure 7. We take homogeneous Neumann boundary conditions and we set C = +∞. The results of the HLL, HR and WBt schemes are presented Table 4, at the final physical time t end = 1s and with 200 discretization cells.

As in the previous examples, Table 4 shows that the WBt scheme and the HR scheme allow the preservation of the lake at rest, even with a dry/wet interface. On the contrary, the HLL scheme provides a poor first-order approximation.

h q L 1 L 2 L ∞ L 1 L 2 L ∞
HLL 2.35e-03 1.38e-02 1.32e-01 5.82e-03 3.55e-02 3.74e-01 HR 0 0 0 0 0 0 WBt 0 0 0 0 0 0 Table 4: Height and discharge errors for the steady state at rest with discontinuous topography and a dry/wet transition.

Flow at rest with emerging bottom

This last experiment at rest involves an emerging bottom (see [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF]). The space domain is [0, 25], and the topography is given by Z(x) = 0.2 -0.05(x -10) 2 + . We take h(x) = (0.15 -Z(x)) + and q(x) = 0 as initial data. We present a graph of the free surface and the topography in Figure 8. For this experiment, we set C = +∞ and we use homogeneous Neumann boundary conditions. The simulation is carried out until the physical time t end = 100s, using 200 discretization cells. The results of the three schemes are displayed in Table 5. This last experiment confirms once again the relevance of using a wellbalanced scheme for the simulation of steady states at rest. Indeed, after Table 5, the HLL scheme only provides a first-order approximation of the steady state, while the HR and WBt schemes provide an exact preservation of this lake at rest.

h q L 1 L 2 L ∞ L 1 L 2 L ∞ HLL

Moving steady states

We now assess the scheme's ability to preserve steady states with a nonzero discharge. Three experiments, deriving from Goutal and Maurel's test cases, are carried out. The subcritical flow, the transcritical flow without shock and the transcritical flow with shock, are presented in [START_REF] Goutal | Proceedings of the 2 nd Workshop on Dam-Break Wave Simulation[END_REF], and will be respectively named GM1, GM2 and GM3. The space domain is 0 < x < 25 and the topography is given by Z(x) = 0.2 -0.05(x -10) 2 + . The boundary conditions are given hereafter, in function of two quantities q 0 and h 0 , whose values depend on the experiment studied:

• on the left boundary, the water height satisfies a homogeneous Neumann condition and the discharge is set to some q 0 ;

• on the right boundary, the water height is set to h 0 when the flow is subcritical (and a homogeneous Neumann boundary condition is prescribed otherwise) and the discharge follows a homogeneous Neumann boundary condition.

In addition, the initial conditions are h + Z = h 0 and q = 0 throughout the domain. The values of q 0 and h 0 are:

• GM1: q 0 = 4.42m 3 /s and h 0 = 2m;

• GM2: q 0 = 1.53m 3 /s and h 0 = 0.66m;

• GM3: q 0 = 0.18m 3 /s and h 0 = 0.33m.

We obtain a transient state followed by a steady state, with uniform discharge q 0 . For GM1 and GM2, this steady state is continuous, and it should thus be exactly obtained by the fully well-balanced scheme WBt. However, the steady state in GM3 involves a stationary shock, which should not be exactly captured by the WBt scheme. The final physical time t end and the constant C are chosen with respect to the experiment, as follows:

• GM1: t end = 500s and C = +∞;

• GM2: t end = 125s and C = 2.5;

• GM3: t end = 1000s and C = 1.1.

On the one hand, for GM1 and GM2, note that q = q 0 and that the steady state equation (2.7) is verified. This equation is nothing but a statement of Bernoulli's principle, and can be rewritten:

q 2 0 2h 2 + g (h + Z) = H,
where H is uniform throughout the domain, and is usually called the total head (see [START_REF] Goutal | Proceedings of the 2 nd Workshop on Dam-Break Wave Simulation[END_REF] for instance). As a consequence, to evaluate the well-balancedness of the scheme on GM1 and GM2, we compute the error to the uniform discharge q 0 and the error for the total head.

Since GM3 presents a stationary shock, the discharge is constant but the total head is not. Indeed, it presents a discontinuity where the shock is located. Therefore, only the error to the uniform discharge is computed for this last experiment.

Comparison with other schemes

To assess the relevance of the WBt scheme, we compare it with the HR and the HLL schemes on moving steady states experiments. We display in Figures 9 -10 the results of the WBt scheme for the GM1 and GM2 experiments. Then, we present in the comparison between the WBt scheme and the HR and HLL schemes. These experiments are performed using a mesh of 200 cells. H q

L 1 L 2 L ∞ L 1 L 2 L ∞
HLL 8.24e-03 1.19e-02 7.41e-02 4.31e-03 1.22e-02 5.19e-02 HR 1.32e-02 1.97e-02 7.48e-02 2.37e-03 6.74e-03 2.74e-02 WBt 1.18e-13 1.25e-13 1.53e-13 6.65e-14 6.99e-14 8.26e-14 H q

L 1 L 2 L ∞ L 1 L 2 L ∞
HLL 2.72e-02 3.50e-02 7.45e-02 1.54e-03 6.16e-03 3.70e-02 HR 4.79e-02 6.07e-02 8.12e-02 8.28e-04 3.30e-03 1.82e-02 WBt 1.67e-14 2.13e-14 4.26e-14 1.47e-14 1.58e-14 2.04e-14

Table 7: Total head and discharge errors for the transcritical flow experiment. show that both HLL and HR schemes provide a first-order approximation of the moving steady state configurations GM1 and GM2, while the proposed WBt scheme exactly preserves (i.e. up to the machine precision) such moving steady states. This result is also observed Figures 910. Moreover, the WBt scheme recovers these steady states after a transient state, even if the steady state is not prescribed as initial condition.

Finally, we turn to the GM3 test case. Since it contains a stationary shock, it is not exactly captured by the WBt scheme, which is designed to capture smooth steady states. The results of the WBt scheme are displayed in Figure 11. Comparisons with respect to ∆x and the scheme used are presented Figure 12, and comparisons with the HR and HLL schemes are presented in Table 8. The experiment is first carried out with 1000 discretization cells, and then with 4000 discretization cells.

From Figure 11, we observe, as expected, that the GM3 experiment is not exactly captured by the WBt scheme. Note the presence of a small inconsistent discontinuity on the free surface in the vicinity of the top of the bump. The amplitude of this discontinuity depends heavily on the constant C. Moreover, this amplitude is reduced when ∆x is reduced, which means that the WBt scheme indeed converges towards the required steady state when ∆x tends to 0.

In the left panel of Figure 12, we observe the expected behavior of the discharge error. Although we do not exactly recover the exact solution, the shock becomes narrower as the number of cells increases. The right panel displays the discharge errors with the HLL, HR and WBt scheme, and shows that the WBt scheme is exact everywhere except in the vicinity of the shock. On the contrary, the results of the HLL and HR schemes are exact, up to the machine precision, only in areas where the topography is flat. q, L 1 q, L 2 q, L ∞ HLL 2.99e-04 1.84e-03 3.89e-02 HR 1.54e-04 1.53e-03 4.00e-02 WBt 2.94e-04 3.35e-03 5.39e-02 Table 8 gives the errors to the steady discharge q 0 , and we note that they are of the same order of magnitude.

Note that the WBt scheme can also be compared to other well-balanced schemes that preserve moving steady states. For instance, in [START_REF] Castro | Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique[END_REF][START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF][START_REF] Xing | Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium[END_REF], error tables are provided, to show that the presented schemes indeed exactly preserve the studied moving steady states. However, it is worth noting that there is no evidence that these schemes are able to capture the steady states obtained after a transient state, contrary to the WBt scheme.

Comparison with the MUSCL scheme

A comparison between the first-order WBt scheme and the MUSCL scheme, with or without the well-balance correction, is performed. The MUSCL scheme refers to relations (5.4) -(5.5) without the well-balancedness correction, i.e. m = 0 and M = 0. The hybrid version, denoted θ-MUSCL, is given by (5.4) -(5.5) with the well-balancedness correction.

First, we present the results for the GM1 and GM2 experiments. Recall that these steady states are exactly captured by the WBt scheme. The results obtained with 200 discretization cells are shown in Tables 9 -10. In the numerical experiments, we choose m = 10 -10 and M = 5.10 -1 .

H q L 1 L 2 L ∞ L 1 L 2 L ∞ WBt 1.
18e-13 1.25e-13 1.53e-13 6.65e-14 6.99e-14 8.26e-14 MUSCL 1.32e-03 4.07e-03 3.38e-02 2.36e-04 9.92e-04 8.27e-03 θ-MUSCL 9.32e-14 1.08e-13 1.56e-13 5.51e-14 5.75e-14 8.88e-14 Tables 9 -10 show that the correction of the MUSCL scheme indeed recovers the well-balance property of the first-order scheme. Without this correction, the MUSCL scheme is not exact for these moving steady states, and the accuracy is reduced to a second-order one.

H q L 1 L 2 L ∞ L 1 L 2 L ∞ WBt 1.
Next, we focus on the GM3 experiment. Since this is not a steady state, the WBt scheme is not exact, which makes the θ-MUSCL scheme less efficient than the MUSCL scheme if the well-balance correction is activated too often. Therefore, we choose m = 10 -10 and M = 10 -4 , so as to make sure that the second-order scheme is used within the shock wave, in order to have a better approximation of this shock. The approximate solution is displayed Figure 13, while the errors are presented Table 11. On Figure 13, we observe that the solution given by the θ-MUSCL scheme (right panel) is more accurate than the one given by the WBt scheme (left panel). The shock is sharper and the inconsistent discontinuity almost disappears thanks to the θ-MUSCL scheme. q, L 1 q, L 2 q, L ∞ 13. Both the MUSCL and the θ-MUSCL schemes are more accurate than the WBt scheme on this discontinuous steady state.

Validation experiments

After the numerical study of steady states, we now turn to performing experiments aimed at validating the scheme. We elect to present two experiments, namely the drain on a non-flat bottom and the vacuum occurrence by a double rarefaction wave over a step. These experiments come from [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF].

Drain on a non-flat bottom

The first validation experiment we propose is the drain on a non-flat bottom (see [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF]). The topography is given by Z(x) = 0.2 -0.05(x -10) 2 + , on the space domain [0, 25]. We take the initial data at rest, as follows: h(x) = 0.5 -Z(x) and q(x) = 0.

Concerning the boundary conditions, we assume that the left boundary is a solid wall and that the drain is done by the right boundary, where we impose an outlet condition on a dry bed (see [START_REF] Dubois | Boundary conditions for nonlinear hyperbolic systems of conservation laws[END_REF][START_REF] Buffard | A sequel to a rough Godunov scheme: application to real gases[END_REF][START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF] for more details on this boundary condition). These boundary conditions are given as follows. Let us denote by h L and q L the left boundary conditions, and by h R and q R the right boundary conditions. Let us assume that (W n i ) i∈ 1,N is the vector containing the approximate solution at time t n . Then, the left boundary condition, a solid wall, is taken as follows:

h L = h n i and q L = 0.
Concerning the right boundary condition, the process to obtain an outlet over a dry bed is detailed in [START_REF] Dubois | Boundary conditions for nonlinear hyperbolic systems of conservation laws[END_REF][START_REF] Buffard | A sequel to a rough Godunov scheme: application to real gases[END_REF]. It consists in choosing the following values at the boundary:

h R = min 1 9g u n N + 2 gh n N 2 , h n N and q R = h R 3 u n N + 2 gh n N .
Note that the outlet on a dry bed boundary condition also requires that the flux at the right interface be the physical flux applied to t (h R , q R ). This boundary condition enables the draining of the water through the right boundary.

The simulation is carried out with the θ-MUSCL scheme, using a discretization of 200 cells, and until the final physical time t end = 1000s. In addition, we take C = 1.35, m = 0.5 and M = 10 -10 . The results are presented Figure 14. On Figure 14, we observe that the θ-MUSCL scheme provides results close to the ones from other schemes, given in [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF][START_REF] Berthon | A positive preserving high order VFRoe scheme for shallow-water equations: a class of relaxation schemes[END_REF][START_REF] Xing | Positivity-preserving high order wellbalanced discontinuous Galerkin methods for the shallow water equations[END_REF][START_REF] Berthon | Hydrostatic upwind schemes for shallow-water equations[END_REF] for instance.

Note that this experiment converges to a steady state at rest (i.e. q(x) = 0 over the whole domain), with the free surface equal to 0.2m at the left of the bump, and with a dry state at its right. Table 12 shows the convergence over time of the θ-MUSCL scheme towards this steady state.

h q L 1 L 2 L ∞ L 1 L 2 L ∞ 150s

Vacuum occurrence by a double rarefaction wave over a step

We then turn to the second validation experiment, a vacuum occurrence deriving from a double rarefaction wave over a step, presented in [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF]. We consider the space domain [0, 25], with a topography given by Z

(x) = 1 ( 25 3 , 25
2 ) . The initial data is given as follows:

h(x) = 10 and q(x) =      -350 if x < 50 3 , 350 otherwise.
We prescribe homogeneous Neumann boundary conditions. The mesh consists in 200 discretization cells, and the simulation is carried out with the θ-MUSCL scheme until a final physical time t end = 0.65s. The parameters are set to C = 1, m = 10 -10 and M = 10 4 . The results are displayed Figure 15. We observe from Figure 15 that the θ-MUSCL scheme provides an approximation that is in good accordance with the ones obtained by several other schemes, given in [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF][START_REF] Berthon | A positive preserving high order VFRoe scheme for shallow-water equations: a class of relaxation schemes[END_REF][START_REF] Bollermann | Finite volume evolution Galerkin methods for the shallow water equations with dry beds[END_REF][START_REF] Xing | High-order finite volume WENO schemes for the shallow water equations with dry states[END_REF] for instance.

2D steady vortex

The proposed WBt scheme is exact for one-dimensional (1D) stationary solutions. Therefore, the order of accuracy of the scheme would have to be tested on a smooth unstationary exact solution of the shallow-water equations with topography. To the extent of our knowledge, there is no such solution referenced. We turn to a two-dimensional experiment, by considering the steady vortex (see [START_REF] Clain | The MOOD method for the non-conservative shallow-water system[END_REF]), to assess the order of accuracy of the scheme. Indeed, the scheme is exact for all the steady states in the x and y directions, but the question of the preservation of a fully 2D steady state arises.

Let us state the shallow-water equations with topography in two dimensions: where h(x, y, t) is the water height, u(x, y, t) and v(x, y, t) respectively represent the water velocities in the x and y directions, and Z(x, y) is the smooth topography.

          
For the steady vortex experiment, the space domain is [-3, 3] × [-3, 3]. The topography is defined by Z(x, y) = 0.2e 0.5(1-r 2 ) , where r 2 = x 2 + y 2 . The exact solution is given by h(x, y) = 1 -1 4g e 2(1-r 2 ) -Z(x, y) ; u(x, y) = ye 1-r 2 ; v(x, y) = -xe 1-r 2 .

The exact solution is chosen as the initial condition, and is prescribed as the boundary conditions. This exact solution is displayed on Figure 16. We use a uniform Cartesian mesh of N square cells, whose edges are of length ∆x. We take the final time t end = 1s. The time step ∆t is given by the following CFL-like condition:

∆t ≤ ∆x 4Λ
, where Λ = max

i∈Z -λ L i+ 1 2 , λ R i+ 1 2 .
We take C = +∞, and m = 0 and M = 0 for the θ-MUSCL extension.

In order to compute the order of accuracy, we consider the results from two simulations with smooth data, one carried out on a mesh composed of N discretization cells, and the other one with N > N cells. The errors are then computed according to (6.1). Let e N be the value of error, in any of the three norms, for a mesh with N cells. The order of accuracy p is then defined as follows: In order to have a relevant computation of the order of accuracy, we take N = 4N in (6.2). Thus, the definition of the order of accuracy used Tables 13 -14 -15 -16 is given by: p = ln(e N )ln(e 4N ) ln 2 .

The accuracy assessment of both schemes is presented in Tables 13 -14 -15 -16. Note that only the x-discharge error is presented. Indeed, the y-discharge error is very similar to the x-discharge error, and the orders of accuracy are the same. 13: First-order scheme, height error and order of accuracy. show the proposed WBt scheme does not exactly preserve the steady vortex, but delivers a first-order approximation of the exact solution.

N L 1 L 2 L ∞ 1024
Similarly, Tables 15 -16 show that the θ-MUSCL scheme preserves the steady vortex up to a second-order approximation. As expected, a fully wellbalanced scheme for the 1D problem does not preserve steady states for 2D geometries.

Conclusion

We have derived a scheme for this system that is well-balanced, positivitypreserving and allows dry/wet transitions. We have also obtained a well-balanced second-order MUSCL extension for this scheme. Several numerical experiments have been performed in order to exhibit the relevant properties of this scheme, which is shown to exactly preserve all steady states at rest and the continuous moving steady states.

In addition to preserving all smooth steady states, the proposed scheme is also shown to exactly capture smooth moving steady states obtained after a transient state. On the contrary, there is no evidence that other fully wellbalanced schemes (for instance, those presented in [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF][START_REF] Xing | Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium[END_REF]) are able to exactly capture these steady state solutions after a transient state. Indeed, for these schemes, the steady state solutions are plugged into the scheme, and the transient state is not simulated. The scheme presented in [START_REF] Castro | Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique[END_REF] is able to capture the steady state solutions after a transient state. However, is it not positive in the presence of large discontinuities in the channel bottom. The non-negativity preservation property and the capture of steady states obtained after a transient state are advantages of the WBt scheme over other exactly well-balanced schemes.

The question of the stability of the scheme is not raised in the present paper. The choice (3.2) of the characteristic velocities ensures that λ L < 0 < λ R , which increases the numerical diffusion of the scheme. As a consequence of this increased diffusion, the scheme is more stable. Therefore, the choice of the characteristic velocities allows us to postpone a more precise study of the stability.

Figure 1 :

 1 Figure 1: Sketch of the function ξ(h; 0.75, 1, √ g, 1) with respect to h.

Figure 2 :

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

Figure 4 :

 4 Figure 4: Correction procedure to ensure positive and consistent intermediate water heights. The line represents the consistency equation (3.17a). If the point (h * L , h * R ) belongs to domain 1 , then h * L and h * R are not modified. However, if (h * L , h * R ) corresponds to a point within domain 2 , we replace (h * L , h * R ) with (ε, (1 -λ L

Theorem 6 .

 6 Consider W n i ∈ Ω for all i ∈ Z. Assume that the intermediate states W for all i ∈ Z, by

Figure 5 :

 5 Figure 5: Graph of θ i with respect to ϕ i , according to (5.6).

Figure 7 :

 7 Figure7: Free surface and topography for some dry lake at rest configurations. The gray area represents the discontinuous topography given by Z.

Figure 8 :

 8 Figure 8: Free surface and topography for the flow at rest with emerging bottom. The gray area represents the topography given by Z.

Figure 9 :

 9 Figure 9: Left panel: free surface and topography for the subcritical flow test case. Right panel: errors for the subcritical flow; the solid line is the total head error and the dashed line is the discharge error.

Figure 10 :

 10 Figure 10: Left panel: free surface and topography for the transcritical flow test case. Right panel: errors for the transcritical flow; the solid line is the total head error and the dashed line is the discharge error.

Figure 11 :

 11 Figure 11: Transcritical flow with shock experiment. The topography is the gray area. Left panel: free surface and topography with 1000 discretization cells. Right panel: free surface and topography with 4000 discretization cells.

Figure 12 :

 12 Figure 12: Transcritical flow with shock experiment. Left panel: discharge error in logarithmic scale with the WBt scheme, with respect to the number of cells. Right panel: discharge error in logarithmic scale with 1000 cells, for different schemes.

Figure 13 :

 13 Figure 13: Free surface and topography (gray area) for the transcritical flow with shock experiment with 1000 discretization cells. Left panel: WBt scheme. Right panel: θ-MUSCL scheme.

WBt 2 .

 2 94e-04 3.35e-03 5.39e-02 MUSCL 1.13e-04 1.94e-03 4.78e-02 θ-MUSCL 1.21e-04 1.94e-03 4.76e-02

Figure 14 :

 14 Figure 14: Drain on a non-flat bottom experiment. Left panel: free surface and topography. Right panel: discharge.

Figure 15 :

 15 Figure 15: Vacuum occurrence by a double rarefaction wave over a step experiment. The gray area represents the topography. Left panel: free surface and topography. Right panel: discharge.

Figure 16 :

 16 Figure 16: Exact solution for the steady vortex experiment. Left panel: free surface. Right panel: velocity norm (the flow is clockwise).

  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 the identity (3.8) is satisfied in a sense we shall define in what follows. Such intermediate states will enforce the well-balancedness of our scheme. To that end, we state the following well-balance principle:

	Therefore, we seek W * L and L and W L = W Principle (WB). The intermediate states W * W * R such that W *

* R are such that W * L = W L and W * R = W R

as soon as W L and W R define a continuous steady state.

Table 2 :

 2 Height and discharge errors for the wet steady state at rest with continuous topography Zc.

Table 3 :

 3 Height and discharge errors for the wet steady state at rest with discontinuous topography Z d .

Table 5 :

 5 Height and discharge errors for the flow at rest with emerging bottom.

Table 6 :

 6 Total head and discharge errors for the subcritical flow experiment.

Table 8 :

 8 Discharge errors for the experiment of the transcritical flow with shock for 1000 discretization cells.

Table 9 :

 9 Total head and discharge errors for the subcritical flow experiment.

Table 10 :

 10 Total head and discharge errors for the transcritical flow experiment.

Table 11 :

 11 Discharge errors for the transcritical flow with shock experiment.

Table 11

 11 confirms the conclusions from Figure

Table 12 :

 12 Water height and discharge errors over time for the drain on a non-flat bottom.

		2.42e-03 3.47e-03 5.51e-03 5.39e-04 5.88e-04 2.24e-03
	600s	2.43e-04 3.71e-04 5.87e-04 1.82e-05 1.92e-05 6.55e-05
	2400s	4.76e-05 7.46e-05 1.17e-04 4.44e-07 4.83e-07 1.66e-06
	19200s 3.43e-06 5.40e-06 8.51e-06 9.63e-09 1.04e-08 3.69e-08

    ∂ t h + ∂ x hu + ∂ y hv = 0, ∂ t hu + ∂ x hu 2 + 1 2 gh 2 + ∂ y (huv) = -gh∂ x Z, ∂ t hv + ∂ x (huv) + ∂ y hv 2 + 1 2 gh 2 = -gh∂ y Z,

Table

		6.37e-03	-	1.61e-02	-	1.19e-01	-
	4096	4.31e-03 0.56 1.08e-02 0.58 7.88e-02 0.59
	16384	2.58e-03 0.74 6.43e-03 0.75 4.58e-02 0.78
	65536	1.43e-03 0.85 3.54e-03 0.86 2.45e-02 0.90
	262144 7.52e-04 0.93 1.86e-03 0.93 1.27e-02 0.95

Table 14 :

 14 First-order scheme, x-discharge error and order of accuracy. -05 1.70 1.38e-04 1.98 1.40e-03 1.84 262144 1.99e-05 1.40 3.71e-05 1.90 3.77e-04 1.89

	N	L 1		L 2		L ∞
	1024	3.90e-02	-	6.72e-02	-	3.00e-01	-
	4096	2.42e-02 0.69 4.14e-02 0.70 1.88e-01 0.67
	16384	1.37e-02 0.82 2.35e-02 0.82 1.10e-01 0.77
	65536	7.29e-03 0.91 1.26e-02 0.90 6.05e-02 0.86
	262144 3.77e-03 0.95 6.54e-03 0.95 3.20e-02 0.92
	N	L 1		L 2		L ∞
	1024	2.41e-03	-	7.55e-03	-	5.52e-02	-
	4096	6.65e-04 1.86 2.14e-03 1.82 1.68e-02 1.72
	16384	1.71e-04 1.96 5.44e-04 1.98 5.00e-03 1.75
	65536	5.25e

Table 15 :

 15 Second-order scheme, height error and order of accuracy.

	N	L 1		L 2		L ∞
	1024	1.33e-02	-	2.27e-02	-	1.04e-01	-
	4096	3.55e-03 1.91 6.08e-03 1.90 3.26e-02 1.67
	16384	8.93e-04 1.99 1.60e-03 1.93 1.09e-02 1.58
	65536	2.25e-04 1.99 4.13e-04 1.95 3.25e-03 1.75
	262144 5.93e-05 1.92 1.06e-04 1.96 9.86e-04 1.72

Table 16 :

 16 Second-order scheme, x-discharge error and order of accuracy.

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 would like to thank the FCT-ANR/MAT-NAN/0122/2012 grant for financial support.

Steady states at rest

The first two experiments concern wet steady states at rest governed by (2.6), with a continuous and discontinuous topography, respectively. The third and fourth experiments are a steady states at rest with a transition between wet and dry areas, and with a discontinuous topography and an emerging bottom (see [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF]), respectively.

Wet lake at rest experiments

Two different wet configurations are considered on the space domain [0, 1]. The first one involves a continuous topography, given by Z c (x) = max(0, 0.5 -2|x -0.5|), while we consider a discontinuous topography

for the second one. For Z ∈ {Z c , Z d }, the initial data are h(x) = 1 -Z(x) and q(x) = 0. These configurations are displayed Figure 6. We prescribe homogeneous Neumann boundary conditions, and we set C = +∞. We carry out the simulations with 200 discretization cells, using the three schemes. Errors computed at the physical time t end = 1s are given in Tables 23.

Tables 23show that both the HR and WBt schemes exactly capture these wet configurations of a lake at rest, while the HLL scheme is only first-order accurate.

Lake at rest with a dry/wet transition

We focus on a lake at rest with a transition between a wet area and a dry area. The space domain is [0, 1] and the topography function is Z(x) =


Ricarda Braukmann

Ingrid Dillo

Peter Doorn

A missing link ? Connecting data and dissertations in the Netherlands

>>> Actions

To consistently connect data and dissertations, changes in policy and technical advances are needed.

Linking Data and Dissertations

• Some repositories can locally store the link between the two in the object's metadata • However, dissertations and datasets need to have a PID and the relationship needs to be known • Moreover, a standard for the way in which the relationship is described, stored and shared still has to be developed > Due to lack of access, information and common standards, most connections remain hidden

Dissertation

• PhDs deposit dissertations in the university repositories

• Repository archives thesis and assigns PID • Information on all dissertations is collected and presented by the national portal NARCIS > Good infrastructure available

Changes in Policy

Awareness about the importance of proper data storage and reference needs to be raised • Universities should implement requirements and offer support and training.

• PhD candidates should archive dissertations and data in trustworthy repositories and provide the connections to the repositories.

• Additional support to create a research data management plan at the beginning of a PhD should be provided.

Technical Advances

Agreement on and adaptations of the relevant metadata standards is needed

•

  PhD dissertations and datasets contain innovative research > Important to store and share the work • However: Challenges when it comes to finding and accessing research data and dissertations Problem: Dissertations are often not connected to related datasets(and vice versa)