Emmanuel Audusse
email: audusse@math.univ-paris13.fr

Christophe Chalons
email: chalons@math.jussieu.fr

Olivier Delestre
email: delestre@polytech.unice.fr

Nicole Goutal
email: nicole.goutal@edf.fr

Magali Jodeau
email: magali.jodeau@edf.fr

Jacques Sainte-Marie
email: jacques.sainte-marie@inria.fr

Jan Giesselmann
email: jan.giesselmann@mathematik.uni-stuttgart.de

Georges Sadaka
email: georges.sadaka@u-picardie.fr

Sediment transport modelling : Relaxation schemes for Saint-Venant -Exner and three layer models

are published or not. The documents may come   L'archive ouverte pluridisciplinaire

Sediment transport modelling : A major issue

In rivers, mean sediment discharge may represent several hundred cubic meters of gravels or silt per year. Therefore, the sediments must be taken into account in order to predict the river bed evolutions. For hydroelectricity managers sediment transport modelling is a major issue. Dams stop water and consequently impact sediment transport, as well as sediment deposition may disturb water intakes or bottom gate opening, for example see Fig. 1. Besides, in some rivers like the Loire river, water intakes of nuclear power station may be protected against the deposition of sand. In order to understand sediment transport dynamics and to suggest managing solutions, EDF has been working on sediment transport modelling tools for bedload and suspension sediment transport. This work focuses on the modelling of bedload transport which refers to gravel transport and pushes aside the transport of fine sediments by suspension. Up to now, one very classical approach is to approximate the solid phase equation by a simplified one : the well-known Exner equation [START_REF] Exner | über die wechselwirkung zwischen wasser und geschiebe in flüssen[END_REF]. The Exner equation is obtained by writing a mass conservation on the solid phase in interaction with the fluid. There is no dynamic effect in the solid phase (1) where u is the velocity in the fluid (Grass [START_REF] Grass | Sediments transport by waves and currents[END_REF] formula) or Q s = f (τ b ) where τ b is the boundary shear stress (see for example Meyer-Peter and Muller [START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF], Einstein [START_REF] Einstein | The bed load function for sediment transportation in open channel flows[END_REF] or Engelund and Fredsoe [START_REF] Engelund | A sediment transport model for straight alluvial channels[END_REF] formulas).

The classical fluid model is the shallow water equations. It is coupled with the solid phase by the bottom level evolution. The coupled model stands 1 :

∂H ∂t + ∂Q ∂x = 0, (2)
∂Q ∂t + ∂ ∂x Q 2 H + g 2 H 2 = -gH ∂z b ∂x , (3)
ρ s (1 -p) ∂z b ∂t + ∂Q s ∂x = 0, (4)
where Q = Hu is the water discharge and H the water height.

There are two major ways to solve this three-equation system : the coupled or the uncoupled approaches. The softwares developed at EDF-RD are based on an uncoupled resolution of the system. Firstly, the hydraulic part of the system is solved by using the software MASCARET [START_REF] Goutal | A finite volume solver for 1D shallow-water equations applied to an actual river[END_REF] and then the computed fluid quantities are sent to the software COURLIS [START_REF] Bertier | One dimensional model for reservoir sedimentation management[END_REF] that solves the Exner equation. The fluid part and the solid one are coupled through the time evolution of the bottom level. In the best way, the shallow-water equations and the Exner equation are coupled at each time step without local iterations. EDF actual modelling tools are used in complex 1 For simplicity's sake, equations are written for the one dimension problem in rectangular channels cases, for example the draw-down of a reservoir with a steep slope and deep zones. Engineers are confronted with the limitation of the codes. Fig 2 is an example of instabilities observed in the case of a dam break calculation using an uncoupled scheme (0.3s after dam break). This result is in agreement with [START_REF] Cordier | Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help[END_REF] where the authors show that, with a splitting method, the instabilities cannot be always avoided in supercritical regions.

Thus alternative approaches must be searched. In this note, we investigate two different approaches. First we consider the coupled approach that consists in solving the system (2)-(4) at once with three unknowns which are H, Q and z b . We then have a hyperbolic system to solve and some specific finite volume schemes have been recently developed, see [START_REF] Benkhaldoun | Solution of the sediment transport equations using a finite volume method based on sign matrix[END_REF][START_REF] Castro Daz | Sediment transport models in shallow water equations and numerical approach by high order finite volume methods[END_REF]. For stiff cases this approach seems to be more robust than the uncoupled one [START_REF] Cordier | Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help[END_REF]. Here we propose a relaxation approach [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] of the model that allows us to deal with a large class of sediment fluxes and we construct the related relaxation solver. Second we go one step further and we derive a new three layers model that is able to propose a more accurate modelling of the phenomenon than the classical Exner approach. The obtained model is quite similar to some models studied in [START_REF] Spinewine | Two-layer flow behaviour and the effects of granular dilatancy in dam-break induced sheet-flow[END_REF]. We also propose a relaxation approach for this new model.

In Section 2, we use a relaxation approach to develop a stable coupled numerical solver for the classical SW-Exner model ( 2)-( 4) and we propose some numerical test cases in Section 3. Then in Sections 4 and 5 we derive the new three layers model and we propose a related uncoupled relaxation model.

A relaxation solver for the Saint-Venant -Exner model

In this section, we consider the classical SW-Exner model ( 2)-( 4) and we propose a relaxation approach for this model. Then we construct the related numerical solver, see also [START_REF] Martin | Modélisation des écoulements à surface libre : étude du ruissellement des eaux de pluie[END_REF]. This work has two main motivations : first, to construct a stable coupled solver for the numerical simulation of model ( 2)-(4) ; second, to introduce the main ideas that will be used to derive the relaxation model associated to the three layers model that is derived in the following sections.

The hyperbolic nature of the SW-Exner model ( 2)-( 4) strongly depends on the formula that is chosen for the sediment flux Q s [START_REF] Cordier | Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help[END_REF]. Moreover the computation of the eigenvalues (that are needed if we want to compute an approximate Riemann solver) of the system is also related to this formula. Here we overcome this difficulty by considering a relaxation model where we relax the sediment flux. We also relax the hydrostatic pressure as it is usually done for shallow water system. This method was first introduced in [START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF]. It consists in the introduction of auxiliary variables and leads to a larger but more suitable system. The method is well-described in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] where a large number of classical hyperbolic solver are interpreted as relaxation schemes. In our case we finally obtain the following model, see also [START_REF] Martin | Modélisation des écoulements à surface libre : étude du ruissellement des eaux de pluie[END_REF] ∂H ∂t + ∂Hu ∂x = 0 (5)

∂Hu ∂t + ∂ ∂x Hu 2 + Π = -gH ∂z b ∂x (6)
∂Π ∂t + u ∂Π ∂x + a 2 H ∂u ∂x = 1 λ (Π - gH 2 2 ) ( 7
)
∂z b ∂t + ∂ Qs ∂x = 0 (8) ∂ Qs ∂t + b 2 H 2 -u 2 ∂z b ∂x + 2u ∂ Qs ∂x = 1 λ ( Qs -Q s ) ( 9
)
where λ is a small parameter and with Π and Qs the relaxed quantities for the fluid pressure and for the sediment flux. The parameters a and b have to be chosen sufficiently large to ensure the stability of the model. This system obviously tends (at least formally) to the original SW-Exner model ( 2)-( 4) when λ tends to zero. The main advantage of the relaxation system (5)-( 9) is that all the eigenvalues are easy to compute whatever the sediment flux Q s is and they are ordered whatever the choice of a and b is, provided a = b. Moreover all the fields are linearly degenerated and then the solution of the Riemann problem is (quite) easy to compute. We can write system (5)-( 9) on the quasi-linear form

∂ t V f + A(V f )∂ x V f = S z + S λ with V f =       H u Π z b Qs       , A(V f ) =        u H 0 0 0 0 u 1/H g 0 0 a 2 /H u 0 0 0 0 0 0 1 0 0 0 b 2 H 2 -u 2 2u       
In the following, we suppose that the parameter a is chosen as it is usually done for the classical Saint-Venant system and that the parameter b is large enough to ensure the stability of the relaxation system (this criterium will be precised later). It follows that a < b and the eigenvalues of matrix A are ordered

u - b H < u - a H < u < u + a H < u + b H . (10)
The associated right eigenvectors are (with δ = ±1)

r u+δ b H =                      1 δ b H 2 a 2 H 2 - a 2 -b 2 gH 3 - a 2 -b 2 gH 3 u + δ b H                      , r u+δ a H =                 1 δ a H 2 a 2 H 2 0 0                 , r u =               1 0 0 0 0               . ( 11
)
We can now compute the solution to the Riemann problem where the relaxation system (5)-( 9) is associated to the initial conditions

V (0, x) = V l if x < 0 V r if x ≥ 0
where we suppose that the water heights H l (x) and H r (x) are bounded below by a positive constant. The solution contains five waves, each associated to one eigenvalue [START_REF] Exner | über die wechselwirkung zwischen wasser und geschiebe in flüssen[END_REF] of the matrix A. These five waves separate six constant states that we denote from left to right : V l , Vl , V * l , V * r , Vr and V r . It follows immediately from the particular form of the right eigenvectors [START_REF] Gerbeau | Derivation of Viscous Saint-Venant System for Laminar Shallow Water; Numerical Validation[END_REF] that Z and Q s are continuous through the three internal waves. Then Zl = Z * l = Z * r = Zr and we will denote this common value by Z * (resp. Q * s ). For u and Π, they are continuous through the intermediate u-wave and then u * l = u * r will be denoted by u * (resp. Π * ). Finally some computations lead to

Z * = u r + δ b Hr Z r -u l -δ b H l Z l u r + δ b Hr -u l -δ b H l - 1 u r + δ b Hr -u l -δ b H l Qr -Ql , and
Q * = u r + δ b Hr Qr -u l -δ b H l Ql u r + δ b Hr -u l -δ b H l - u r + δ b Hr u l -δ b H l u r + δ b Hr -u l -δ b H l (Z r -Z l ) .
Then we can compute the intermediate water heights

1 Hl = 1 H 2 l - 2g a 2 -b 2 (Z l -Z * ) 1 2 , 1 Hr = 1 H 2 r - 2g a 2 -b 2 (Z r -Z * ) 1 2
, and the intermediate velocities and pressure

ūl = u l + b 1 Hl - 1 H l , ūr = u r -b 1 Hr - 1 H r ,
and

Πl = Π l -a 2 1 Hl - 1 H l , Πr = Π r -a 2 1 Hr - 1 H r .
We finally compute the internal quantities

Π * = Πl + Πr 2 - a 2 (ū r -ūl ), u * = ūl + ūr 2 - 1 2a ( Πr -Πl ), and 1
H * l = 1 Hl - 1 a 2 (Π * -Πl ), 1 H * r = 1 Hr - 1 a 2 (Π * -Πr ).
The relaxation model can be extended to the case where one of the water height H l or H r vanishes by considering the method proposed by Bouchut [START_REF] Bouchut | Gravity driven shallow water models for arbitrary topography[END_REF]. We now build the numerical scheme as follows. Starting from a given sequence V n i computed at time t n , we first compute the time step by considering the CFL condition

∆t n ≤ ∆x max i (|u n i | + b/H n i )
.

Let us now denote by V i+1/2 (t, x) the exact solution of the Riemann problem with initial conditions

V (0, x) = V i if x < 0 V i+1 if x ≥ 0
We compute the numerical solution V n+1 i at the next time step t n+1 by considering the mean value of this exact solution on the cell C i at time t n+1 , i.e.

V n+1 i = 1 x i -x i-1/2 xi x i-1/2 V i-1/2 (∆t n , x)dx + 1 x i+1/2 -x i x i+1/2 xi V i+1/2 (∆t n , x)dx.
Notice that for the conservative quantities H and Z we can also consider a more classical formula by using fluxes at interfaces

H n+1 i = H n i - ∆t n |C i | F H (V i+1/2 (∆t n , x i+1/2 )) -F H (V i-1/2 (∆t n , x i-1/2 )) .

Numerical test cases

Analytical solution

We first consider a very simple test case for which it is possible to compute an analytical solution [START_REF] Berthon | An analytical solution of shallow water system coupled to exner equation[END_REF]. For the Exner law, we consider the Grass formula (1) with A g = 0.005 and m = 3. The initial data are the following :

q(0, x) = q 0 , h(0, x) = q 0 (αx + β) 1/3 , q 2 0 2h(0, x) 2 + g(z b (0, x) + h(0, x)) = C.
It follows from easy computations, see [START_REF] Berthon | An analytical solution of shallow water system coupled to exner equation[END_REF], that this solution is a stationary solution for the fluid part

q(t, x) = q(0, x), h(t, x) = h(0, x), ∀t > 0,
and that the bottom topography is given by the relation

z b (t, x) = z b (0, x) -A g αt, ∀t > 0.
In Figure 3, we compare the numerical solution with the analytical one for two different number of grid points. We can see that the numerical solution converges to the analytical one. Nevertheless the scheme is quite diffusive since for a small number of point the numerical solution is quite far from the analytical one.

Dam break on movable bottom

We now consider the very classical case of a dam break over a flat bottom but we allow the bottom to evolve in time with respect to the Exner law (4). As initial conditions, the dam is located at x = 5m, with H f = 2m upstream the dam and H f = 0.125m downstream. The topography is initially flat and the flow is at rest u f = 0m/s. For the Exner law, we still consider the Grass formula (1) with the same parameters and we let evolve in time until t = 1s, for J = 4000 cells in space. The results obtained with a splitting approach where the fluid equations and the Exner law are solved separately, presented in Figure 2 in the introduction, exhibit the presence of numerical instabilities (with J = 400 cells in space). The results obtained with our coupled approach are shown in Figure 4. It shows that for this kind of stiff cases the coupled approach is much more stable. It confirms the conclusion of [START_REF] Cordier | Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help[END_REF].

"Steady" flow over a movable bump

The last numerical test case is concerned with the extension to moveable bottom of the so-called "steady flow over a bump" experiment. Here the flow will not reach a steady state since the bump will be deformed by the flow. As initial conditions, we consider a subcritical steady state over a bump for the shallow water system as in [START_REF] Cordier | Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help[END_REF], i.e.

       Q f (x, t = 0) = 0.5 z b (x, t = 0) = 0.1 + 0.1e -(x-5) 2 u 2 f 2 + g(H f + z b ) = 6.386
.

We still consider the Grass formula with the same parameters and we let evolve in time until t = 30s, for J = 256 and J = 4000 cells in space. The results are presented in Figure 5. With J = 256 cells, we notice that the scheme is very diffusive. For J = 4000 cells, we recover the result obtained in [START_REF] Cordier | Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help[END_REF]. A second order accurate extension of the scheme has to be developed to enhance the precision of the results. We then consider the same initial conditions but with a slightly larger A g = 0.007, thus the ground is more erodible. The simulation is ended at t = 5s. It was noticed in [START_REF] Cordier | Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help[END_REF] that the splitting approach was not adapted for this test case, i.e. they observed some oscillations. Here for J = 4000 cells, we do not observe any oscillation, see Figure 6, and the solution is relevant. This last numerical test confirms the interest of the coupled approach since it is stable for a large class of test cases.

A two phases Navier Stokes system

The erosion phenomenon is a fluid/structure interaction problem. It can be formulated in the following way: the displacement of a pure fluid phase that interacts with a movable bed made of a mixture between the pure fluid phase and the static bed components, see Fig 7 . The density of the mixture fluid/solid varies with respect to the concentration of the static bed components. In the following superscripts b, s and f deal respectively with the static bed, mobile bed and pure water layer. We also denote by H the depths of each layer, u their mean horizontal velocities, c the sediment concentrations, ρ the mean density, µ the dynamic viscosity and g the gravity acceleration. The free surface is written η(x, t), the bottom topography is z b (x, t) and z s (x, t) stands for the interface between the fluid and the mobile bed.

The fluid part

We begin by considering the two-dimensional hydrostatic Navier-Stokes system [START_REF] Lions | Incompressible models[END_REF] describing a free surface gravitational flow moving over a bottom topography. We denote with x and z the horizontal and vertical directions, respectively. The system has the form

∂ρ f ∂t + ∂ρ f u f ∂x + ∂ρ f w f ∂z = 0, ( 12
)
∂ρ f u f ∂t + ∂ρ f (u f ) 2 ∂x + ∂ρ f u f w f ∂z + ∂p f ∂x = ∂Σ f xx ∂x + ∂Σ f xz ∂z , (13)
∂p f ∂z = -ρ f g + ∂Σ f zx ∂x + ∂Σ f zz ∂z , (14)
and we consider solutions of the equations for, see Fig.

7

t > t 0 , x ∈ R, z s (x, t) ≤ z ≤ η(x, t).
The chosen form of the viscosity tensor is

Σ f xx = 2µ f ∂u f ∂x , Σ f xz = µ f ∂u f ∂z + ∂w f ∂x , Σ f zz = 2µ f ∂w f ∂z , Σ f zx = µ f ∂u f ∂z + ∂w f ∂x .
The fluid density ρ f (x, t) is assumed to be constant or to depend on the spatial and temporal distribution of the concentration of the sediment C f (x, t) in the fluid layer , namely

ρ f = ρ(C f ), (15)
and C f is governed by a transport-diffusion equation

∂ρ f C f ∂t + ∂ρ f u f C f ∂x + ∂ρ f w f C f ∂z = µ f C ∂ 2 C f ∂x 2 + µ f C ∂ 2 C f ∂z 2 + S, (16)
where µ f C is the sediment diffusivity and S the source terms (chemical reactions...).

In the previous system we may consider the Boussinesq assumption :

ρ f = ρ f 0 + ǫρ f 1 (C f ), ( 17
)
where ǫ is a small parameter. In that case, at the leading order, the density depends on the concentration only in the gravity term.

The mobile bed

We also consider that the mobile bed can be modelled by the two-dimensional hydrostatic Navier-Stokes system but with a varying density ρ s

∂ρ s ∂t + ∂ρ s u s ∂x + ∂ρ s w s ∂z = 0, (18)
∂ρ s u s ∂t + ∂ρ s (u s ) 2 ∂x + ∂ρ s u s w s ∂z + ∂p s ∂x = ∂Σ s xx ∂x + ∂Σ s xz ∂z , (19)
∂p s ∂z = -ρ s g + ∂Σ s zx ∂x + ∂Σ s zz ∂z , (20)
and we consider solutions of the equations for, see Fig. 7 t > t 0 , x ∈ R, z b (x, t) ≤ z ≤ z s (x, t).

The viscosity tensor Σ s has the same form as in the fluid part. The density ρ s (x, t) is assumed to be constant or to depend on the spatial and temporal distribution of the concentration of the sediment C s (x, t) in the mobile bed, see relations ( 15) and ( 16).

The boundary conditions

The systems ( 12)-( 14) and ( 18)-( 20) are completed with boundary conditions. The outward unit normal vector to the free surface n η and the upward unit normal vectors n f s and n sb respectively to the interfaces z s and z b are given by

n η = 1 1 + ∂η ∂x 2 -∂η ∂x 1 , n j = 1 1 + ∂zj ∂x 2 - ∂zj ∂x 1 , j = f s, sb.
We then denote Σ j T the total stress tensor, which has the form:

Σ j T = -p j I d + Σ j xx Σ j xz Σ j zx Σ j zz , j = f, s.

Free surface conditions

At the free surface we have the kinematic boundary condition

∂η ∂t + u f η ∂η ∂x -w f η = 0, (21)
where the subscript η indicates the value of the considered quantity at the free surface.

We denote τ w the surface wind stress, for which the following expression [START_REF] Imberger | Physical processes in lakes and oceans[END_REF] is considered:

τ w = C D ρ a |V w |V w ,
where V w is the wind velocity, ρ a the air density and C D the wind drag coefficient. Assuming negligible the air viscosity, the continuity of stresses at the free boundary imposes

Σ f T n η = -p a n η + τ w t η , (22)
where p a = p a (x, t) is a given function corresponding to the atmospheric pressure and t η is the unit vector orthogonal to n η . Relation ( 22) is equivalent to

n η • Σ f T n η = -p a , t η • Σ f T n η = τ w .
In the following, we will assume p a = 0.

Boundary conditions at the interface z s

At the interface z = z s (x, t), we impose the continuity of the normal velocity

u f f s ∂z s ∂x -w f f s = u s f s ∂z s ∂x -w s f s (23)
and we can define a quantity e f s that describes the exchange between the fluid and the mobile bed, see Fig. 7 e f s = -

∂z s ∂t -u f f s ∂z s ∂x + w f f s = - ∂z s ∂t -u s f s ∂z s ∂x + w s f s (24)
Moreover the continuity of stresses at the moving boundary z s also imposes the continuity of the pressure and the continuity of tangential stresses

t f s • Σ f T n f s = t f s • Σ s T n f s = τ f s , (25)
where τ f s models inter-phase drag forces which can be expressed as

τ f s = D(H f + H s )(u f f s -u s f s ),
where

D = D(|u f f s -u s f s |; σ) is a drag function depending in general on |u f f s -u s f s |
, and a set of physical parameters σ (e.g. specific densities, particle diameter), see [START_REF] Pelanti | A riemann solver for single-phase and two-phase shallow flow models based on relaxation. relations with roe and vfroe solvers[END_REF] for a more complete description.

Boundary conditions at the bottom interface z b

At the interface z = z b (x, t), we have to impose the continuity of the normal velocity (which vanishes since the bottom is supposed to be at rest) and we can define a quantity e f s that describes the exchange between the static and mobile beds, see Fig. 7 e sb = -

∂z b ∂t . ( 26
)
For the stresses at the bottom, we consider a wall law under the form

t sb • Σ s T n sb = τ b = κ(u s sb , H)u s sb • t sb , ( 27
)
where t sb is a unit vector satisfying t sb • n sb = 0. Due to thermo-mechanical considerations, in the sequel we will suppose κ(u s sb , H) ≥ 0, and κ(u s sb , H) will be often simply denoted by κ. If κ(u s sb , H) is constant then we recover a classical Navier friction condition.

A two-phase shallow granular flow model

Derivation

Let us introduce two small parameters ε f and ε s defined by

ε f = h f λ , ε s = h s λ , ( 28
)
where λ is a characteristic dimension along x and h f ,h s are characteristic dimensions along the vertical axis in the fluid and mobile bed layers respectively. We consider in this section that the shallow water assumption is valid for the fluid and for the mobile bed. This means the models ( 12)-( 14) and ( 18)-( 20) completed with boundary conditions ( 21)-( 27) can be approximated up to O((ε f ) 2 , (ε s ) 2 ) terms by the coupled system (29)-( 33)

Pure Water :

∂H f ∂t + ∂ ∂x (H f u f ) = ρ s ρ f max(0, e f s ) + min(0, e f s ), (29)
∂(H f u f ) ∂t + ∂ ∂x H f (u f ) 2 + g 2 (H f ) 2 = -gH f ∂(z b + H s ) ∂x + τ w ρ f - τ f s ρ f + ρ s ρ f u s max(0, e f s ) + u f min(0, e f s ), (30)
Mobile Bed Layer :

∂H s ∂t + ∂ ∂x (H s u s ) = - ρ f ρ s min(0, e f s ) -max(0, e f s ) + ρ b ρ s max(0, e b ) + min(0, e b ), (31)
∂(H s u s ) ∂t + ∂ ∂x H s (u s ) 2 + ρ f ρ s gH s H f + g 2 (H s ) 2 = g ρ f ρ s H f ∂z s ∂x -g H s + ρ f ρ s H f ∂z b ∂x + τ f s ρ s - τ b ρ s -u s max(0, e f s ) - ρ f ρ s u f min(0, e f s ) + u s min(0, e b ), (32)
Static Bed Layer :

∂z b ∂t = -max(0, e sb ) - ρ s ρ b min(0, e sb ). (33)
This shallow water type model is derived from the original Navier-Stokes equations ( 12)-( 14) by using the shallow water assumption (28) and an integration along the vertical direction. This process is very similar to the one used in [START_REF] Gerbeau | Derivation of Viscous Saint-Venant System for Laminar Shallow Water; Numerical Validation[END_REF] to obtain the classical shallow water model. Let us give some details for the fluid layer in the constant density case. First we integrate the mass equation ( 12) on the fluid layer to obtain

∂ ∂t η zs ρ f dz + ∂ ∂x η zs ρ f u f (t, x, z)dz -ρ f η ∂η ∂t + ρ f f s ∂z s ∂t + -ρ f η u f η ∂η ∂x + ρ f f s u f f s ∂z s ∂x + ρ f η w f η -ρ f f s w f f s = 0. ( 34
)
Then we use boundary conditions ( 21) and (24) to write

∂ ∂t η zs ρ f dz + ∂ ∂x η zs ρ f u f (t, x, z)dz + ρ f f s e f s = 0.
Next we introduce the layer depth H f = η -z s and the layer velocity

ūf = 1 ρ f H f η zs ρ f u f (t, x, z)dz, to obtain ∂H f ∂t + ∂ ∂x (H f ūf ) + ρ f s
ρ f e f s = 0. We finally consider an upwind estimation of the interface density ρ f s to derive equation (29). The process is similar for the momentum equation and for the mobile bed layer. Note that comparison between equations (26) and (33) is quite surprising. It should be noted that equation (26) has to be seen as an interface condition between static and mobile bed layers whereas equation ( 33) is obtained after an integration process over the static bed layer.

In the previous derivation, we assumed that the densities are constant in each layer. The derivation can be done for the general case where the densities depend on the concentration of the sediment in each layer. It is also possible to consider an intermediate version by taking into account the Boussinesq assumption [START_REF] Martin | Modélisation des écoulements à surface libre : étude du ruissellement des eaux de pluie[END_REF].

It remains to give an expression for the exchange terms e f s and e sb . In the constant densities case, the total mass conservation imposes a relation between these two parameters

e f s = - ρ b -ρ s ρ s -ρ f e sb . (35)
Possible expressions for e sb are discussed in [START_REF] Spinewine | Two-layer flow behaviour and the effects of granular dilatancy in dam-break induced sheet-flow[END_REF]. Here we choose a relation where e sb is proportional to the jump of the shear stress, and we consider a Mohr-Coulomb bed strength

e sb = -1 ρ b |u s | (ρ s -ρ f )H s tan φ + τ b , τ b = C b ρ s |u s | 2 sgn(u s ).labeleq : esb (36)

Steady states

For a special choice of constitutive functions for e f s , e sb , τ f s , τ b , τ w and constant densities ρ b > ρ s > ρ f > 0 of the layers we determine steady states of the problem. Let us first give the missing consitutive relation for the interface shear stress :

τ f s = C s ρ f ρ s (u f -u s ) 2 sgn(u f -u s ). (37)
We also consider there is no wind and so τ w = 0. To investigate steady states of the system(29)-(33) with the above constitutive laws, we introduce the following expansions, for a small parameter ε:

H i (x, t) = H i 1 (x) + εH i 1 (x, t) + o(ε), u i (x, t) = 0 + εu i 1 (x, t) + o(ε), for i ∈ {f, s} (38)
Then we obtain the following Lemma, whose proof is immediate by inserting the expansions (38), into the equations ( 29)-(33)

Lemma 5.1. Provided (38) is a solution of (29)-(33) with consitutive laws given by (35)-(37), it satisfies the following conditions

H s 0 = 0, ∂ ∂x (H f 0 + z b ) = 0, H s 1 = 0.
This shows that we can only have steady states when H s 0 . Thus, the (expected) "lake at rest" condition

∂ ∂x (H f 0 + H s 0 + z b ) = 0
is simplified accordingly. The latter condition H s 1 = 0 shows that steady states are incompatible even with small (order ε) perturbations of H s . This indicates some sort of instability of the steady states with respect to this variable.

Linear Stability

We also tried to determine how the fluxes between the phases in (29)-(33) have to be chosen, such that a linearly stable system is obtained. To this end we started with a simplified model where there is no mass exchange between the phases, the bottom is flat and the pure water and the mobile bed have constant prescribed densities ρ f , ρ s , respectively. The only exchange between the phases, we considered, was a strong friction term. As we were not able to determine in which cases the system is linearly stable, we will just explain the problems we encountered, when trying classical linear stability analysis.

The simplified system with quadratic friction reads

∂ ∂t H f + ∂ ∂x (H f u f ) = 0 (39) ∂ ∂t (H f u f ) + ∂ ∂x (H f (u f ) 2 + g 2 (H f ) 2 ) + gH f ∂ ∂x H s = Cλ 2 |u f -u s |(u f -u s ) (40) ∂ ∂t H s + ∂ ∂x (H s u s ) = 0 (41) ∂ ∂t (H s u s ) + ∂ ∂x (H s (u s ) 2 + g 2 (H s ) 2 ) + rgH s ∂ ∂x H f = -rCλ 2 |u f -u s |(u f -u s ) -Kλ 2 |u s |u s , (42)
where Kλ 2 , Cλ 2 are friction coefficients and r := ρ f /ρ s . To determine the behavior for strong friction and long times we change the time variable and expanded the velocities

s = λ -1 t, u i = 0 + λ -1 v i + o(λ -1 ) for i ∈ {f, s}. (43)
We insert (43) in ( 39)-( 42) and only consider the leading order equations, which are of order λ -1 for the mass conservation equations and of order λ 0 for the momentum balances. We obtain

∂ ∂s H f + ∂ ∂x (H f v f ) = 0 (44) ∂ ∂x ( g 2 (H f ) 2 ) + gH f ∂ ∂x H s = C|v s -v f |(v s -v f ) (45)
∂ ∂s H s + ∂ ∂x (H s v s ) = 0 (46) ∂ ∂x ( g 2 (H s ) 2 ) + rgH s ∂ ∂x H f = -rC|v s -v f |(v s -v f ) -K|v s |v s . (47)
The equations ( 44)-(47) can be rewritten as a problem for H s , H f only, by solving (45),(47) for v s , v s and inserting the result into (44),(46). Thereby we obtain

∂ ∂s H f + ∂ ∂x H f -sgn(A) |A| K -sgn(H f ∂ ∂x (H f + H s )) g C H f ∂ ∂x (H f + H s ) = 0 (48) ∂ ∂s H s + ∂ ∂x -H s sgn(A) |A| K = 0, (49)
where

A := gr(H f + H s ) ∂ ∂x H f + g(rH f + H s ) ∂ ∂x H s .
To find out whether the system (48),(49) is linear stable or not one investigates the behavior of perturbations of solution H s 0 , H f 0 to (48),(49) as follows. For small values of ε we define

H i = H i 0 + εH i 1 + o(ε) for i ∈ {f, s}, (50)
insert (50) into (48),(49) and study the arising system of equations for H f 1 , H s 1 . It has the form

∂ ∂s H f 1 H s 1 + ∂ ∂x A H f 1 H s 1 + D ∂ ∂x H f 1 H s 1 = 0 (51)
where A, D : [0, ∞) × R → R 2×2 will be given in the sequel. The classical strategy is to exploit the properties of A and D to determine whether or not perturbations are damped over time. This seems to be impossible here because of the involved structure of A and D. In particular D is singular for

A 0 := gr(H f 0 + H s 0 ) ∂ ∂x H f 0 + g(rH f 0 + H s 0 ) ∂ ∂x H s 0 = 0 or H f ∂ ∂x (H f 0 + H s 0 ) = 0.
The latter condition indicates that we cannot consider lineralizations around constant solutions, in this setting. Finally let us give the component-wise definitions of A and D such that (51) is valid,

A 11 = -sgn(A 0 ) |A 0 | K -sgn(H f 0 ∂ ∂x (H f 0 + H s 0 )) g C H f 0 ∂ ∂x (H f 0 + H s 0 ) -sgn(A 0 )H f 0 gr 4|A 0 |K ∂ ∂x (H f 0 + H s 0 ) -sgn(H f 0 ∂ ∂x (H f 0 + H s 0 )) g 4CH f 0 ∂ ∂x (H f 0 + H s 0 ) ∂ ∂x (H f 0 + H s 0 ) A 12 = -sgn(A 0 ) g 4|A 0 |K H f 0 ∂ ∂x (rH f 0 + H s 0 ) A 21 = -sgn(A 0 ) gr 4|A 0 |K H s 0 ∂ ∂x (H f 0 + H s 0 ) A 22 = -sgn(A 0 ) |A 0 | K -sgn(A 0 ) g 4|A 0 |K H s 0 ∂ ∂x (rH f 0 + H s 0 ) D 11 = -sgn(A 0 ) gr 4|A 0 |K H f 0 (H f 0 + H s 0 ) -sgn(H f 0 ∂ ∂x (H f 0 + H s 0 )) g 4CH f 0 ∂ ∂x (H f 0 + H s 0 ) H f 0 D 12 = -sgn(A 0 ) g 4|A 0 |K H f 0 (rH f 0 + H s 0 ) -sgn(H f 0 ∂ ∂x (H f 0 + H s 0 )) g 4CH f 0 ∂ ∂x (H f 0 + H s 0 ) H f 0 D 21 = -sgn(A 0 ) gr 4|A 0 |K H s 0 (H f 0 + H s 0 ) D 22 = -sgn(A 0 ) g 4|A 0 |K H s 0 (rH f 0 + H s 0 ).

A relaxation approach

The system (29)-( 33) is quite complex : the eigenvalues are impossible to compute analytically and the field are fully non linear. It follows that it is not possible to use approximate Riemann solvers. As in the first part of this work, a possible way to cure this problem is to use a relaxation approach. So we introduce three auxiliary quantities, namely the fluid pressure, the "fluid bottom" and the related "fluid bottom" flux. For simplicity we consider the model with constant densities and we give the details of the method for the upper fluid layer. We first consider the "fluid bottom"

z s (t, x) = z b (t, x) + H s (t, x).
Using equations (31) and (33), we obtain (52)

∂z s ∂t + ∂ ∂x (H s u s ) = - ρ f ρ s min(0,
Then we can write the following relaxation model for the fluid (we do not consider the zero order right hand side and so we omit the algebraic exchange terms)

∂ t H f + ∂ x H f ūf = 0 ( 53
)
∂ t H f ūf + ∂ x [H f ū2 f + π f ] = -gH f ∂ x z f ( 54
)
∂ t π f + ūf ∂ x π f + a 2 f H f ∂ x ūf = 1 λ f [g H 2 f 2 -π f ] ( 55
)
∂ t z f + ∂ x Q f = 1 λ f [z s -z f ] ( 56
)
∂ t Q f + b 2 H 2 f -ū2 f ∂ x z f = 1 λ f [H s u s -Q f ], (57)
where Π f (resp. z f and Q f ) denotes the auxiliary quantity that allows us to relax the fluid pressure (resp. the "fluid bottom" and "fluid bottom" flux). Here the parameters a f and b f have to be chosen sufficiently large to ensure the stability of the model. Starting from equations (31)-(32) and using equations ( 29

∂ t z s + ∂ x Q s = 1 λ s [z b + ρ f ρ s H f -z s ] (61)
∂ t Q s + b 2 H 2 s -ū2 s ∂ x z s = 1 λ s [ ρ f ρ s H f u f -Q s ]. ( 62
)
Many remarks have now to be done. First, the relaxed model ( 53)-(61) obviously tends (at least formally) to the original model ( 29)-(33) when the parameters λ f and λ s tend to zero. Second, we note that the two systems (53)-( 57) and ( 58)-(62) are only coupled through the right hand side relaxation terms that appear in equations ( 56)-( 57) and ( 60)-(62). Since we see in the first part of this work that these relaxation terms are not explicitely computed (they are taken into account at once when we update the auxiliary quantities at the new time step), it follows that we obtained two decoupled systems ! Third, we remark that systems (53)-( 57) and ( 58)-( 62) are very similar to system (5)-( 9) that we obtained in the first part of this work. It follows that the relaxation solver that has been presented previously can be applied without any big change to this much more complex case.

Conclusion

In this note, we consider the morphodynamic coupling in rivers. We first propose a coupled finite volume numerical scheme for the classical Saint-Venant -Exner model that is based on a relaxation approach. We compute some numerical test cases to show the stability of this new solver even for stiff cases. Then we derived a new three layers model that we hope to be more relevant in the physical sense than the classical approach and we show that the relaxation framework can be extended to this new model. Further work has to be done to improve the properties of the relaxation solver and to extend it to the three-layer model.

Figure 1 .

 1 Figure 1. Example of silt deposition pattern in reservoirs (water flows from the back to the front).

  z b is the bed elevation, Q s the bed load, p the porosity of the gravel bed. The bed load may be expressed by empirical formulae of the form Q s = A g (u)|u| m-1 u

Figure 2 .

 2 Figure 2. Dam break over a moveable bottom with a splitting approach -Free surface (top) and Bottom topography (bottom)

Figure 3 .

 3 Figure 3. Comparison between numerical and analytical solutions for the Saint-Venant -Exner model ; 100 (left) and 4000 (right) grid points

  Bottom topography at t=1 s -J=4000 z b +H f at t=1 s

Figure 4 .

 4 Figure 4. Dam break over a moveable bottom -Free surface (top) and bottom topography (bottom)

Figure 5 .

 5 Figure 5. Flow over a moveable bump, coarse (left) and fine (right) meshes

Figure 6 .

 6 Figure 6. Flow over a moveable bump, a more erodible bottom (A g = 0.007)

Figure 7 .

 7 Figure 7. Description of the domain.

  e f s ) -max(0, e f s ) + ρ b ρ s max(0, e b ) + min(0, e b ) -max(0, e sb ) -ρ s ρ b min(0, e sb ).

  ) and (33), we obtain the same kind of model for the mobile bed layer where all the quantities are denoted by subscript s∂ t H s + ∂ x H s u s = 0 (58) ∂ t H s u s + ∂ x [H s u 2 s + π s ] = -gH s

						∂z s ∂x			(59)
	∂ t π s + ūs ∂ x π s +	a 2 s H s	∂ x u s =	1 λ s	[g	H 2 s 2	-g	ρ f ρ s	H s H f -π s ]	(60)

Acknoledgments

The authors thank B. Spinewine for discussions about the three-layer model and C. Berthon for discussion about the relaxation framework. The authors also thank M.O. Bristeau for her comments that greatly improve the manuscript. O. Delestre, J. Sainte-Marie, G. Sadaka and J. Giesselmann thank EDF for financial support during the CEMRACS. E. Audusse and C. Chalons thank CNRS and INSMI for PEPS financial support. The authors finally greatly thank the CIRM for its warm welcome during the CEMRACS.