MULTILAYER SAINT-VENANT EQUATIONS OVER MOVABLE BEDS

We introduce a multilayer model to solve three-dimensional sediment transport by wind-driven shallow water flows. The proposed multilayer model avoids the expensive Navier-Stokes equations and obtains stratified horizontal flow velocities as vertical velocities are relatively small and the flow is still within the shallow water regime. To model the bedload transport we consider an Exner equation for morphological evolution accounting for the velocity field on the bottom layer. As a numerical solver, we apply a kinetic scheme using the finite volume discretization. Preliminary numerical results are presented to demonstrate the performance of the proposed model and to confirm its capability to provide efficient simulations for sediment transport by wind-driven shallow water flows.


Introduction
The main concern of morphodynamics is to determine the evolution of bed levels for hydrodynamic systems such as rivers, estuaries, bays and other nearshore regions where water flows interact with the bed geometry. Example of applications include among others, beach profile changes due to severe wave climates, seabed response to dredging procedures or imposed structures, and harbour siltation. The ability to design numerical methods able to predict the morphodynamic evolution of the coastal seabed has a clear mathematical and engineering relevances [8,9,10,11,12,13]. In practice, morphodynamic problems involve coupling between a hydrodynamic model, which provides a description of the flow field leading to a specification of local sediment transport rates, and an equation for bed level change which expresses the conservative balance of sediment volume and its continual redistribution with time. In the current study, the hydrodynamic model is described by a multi-layer shallow water equations and the sediment transport is modelled by the Exner equation.
Morphodynamic coupling between classical shallow water system and Exner equation has been recently studied in some works [4,5]. Here we would like to extend this coupling strategy to a multilayer shallow water model. The interest of the multilayer approach lies in the possibility to obtain a detailed description of the velocity in the flow while keeping a relative simplicity and a great robustness in the numerical procedure.On one hand and as opposed to what happens for classical shallow water model the multilayer approach allows us to consider quite complex flows as wind driven circulation in a closed basin (lake, estuary...). On the other hand and as opposed to what happens when 3d incompressible hydrostatic free surface Navier-Stokes equations are considered (as it is often the case in oceannographic community) the multilayer model allows us to consider a 2d problem in a fixed domain. Moreover the hyperbolic structure of the system and its link with the classical shallow water model gives some keys to design very robust numerical procedures. Some numerical evidences of the abilities of the model to reproduce physical behavior were obtain in [3].
In this paper, first the governing equations for the morphodynamic problems are formulated. Thereafter, the solution procedure employed to solve the morphodynamic problems is presented and some preliminary results obtained for multilayer models are discussed. Concluding remarks end the paper.

Governing Equations for Sediment Transport Problems
In this section we describe the physical model used for modelling the sediment transport by wind-driven shallow water flows. Here, the Multilayer Saint-Venant equations are briefly recast for the hydrodynamics followed with a short description of the Exner equation for the morphodynamics.

Multilayer Saint-Venant Equations
The governing equations are obtained starting from the three-dimensional hydrostatic incompressible Navier-Stokes equations by considering a vertical P 0 type discretization of the horizontal velocity. This vertical discretization defines some layers in the flow and the equations are vertically integrated on each layer separately. A global mass equation for the whole flow is obtained by adding all the layer mass equations. It is coupled with N momentum equations, one for each of the N layers introduced in the vertical discretization. We refer to [3] for a detailed derivation of the model. We notice that the layers do not refer to physical interface between unmiscible fluids but to a meshless discretization of the flow. Hence the possibility of water exchange between the layers is included in the model. The great interest of the strategy we present is to preserve an accurate description of the velocity profile but to deal with a 2d x − y fluid model and thus to avoid the difficult question of meshing a 3d moving domain for which the free surface may present very sharp profiles -dam break problems, hydraulic jumps. For simplicity we consider in this paper the 1d version of the model. Hence the equations of the multilayer shallow water system with are given by where h(t, x) is the water height of the flow, u α (t, x) is the local water velocity for the αth layer, B(x) is the bottom topography and g the gravitational acceleration. In (1), λ α denotes the relative size of the αth layer such that M α=1 λ α = 1.
The pressure term p α is defined as and F α is the external force acting on the layer α and accounting on the friction and momentum exhange effects. Thus where the three first terms are related to friction effects and the two last ones are related to the momentum exchanges between the layers that are defined through the vertical P 0 discretization of the flow. Term F u is related to advection process whereas term F p takes into account pressure effect. The bed friction forcing term F b is acting only on the lower layer and the wind-driven forcing term F w is acting only on the upper layer. They are given by with δ kα represents the Kronecker delta, ρ w the water density, τ b and τ ω are respectively, the bed shear stress and the shear of the blowing wind defined by the water and wind velocities as where C b is the bed friction coefficient, which may be either constant or estimated as 1 /n b is the Chezy constant, with n b being the Manning roughness coefficient at the bed, ω is the velocity of wind at 10 m above water surface and C w is the coefficient of wind friction defined as [6] where ρ a is the air density. The vertical kinematic eddy viscosity term F µ takes into account the friction between neighbouring layers and is defined as where µ is the kinematic eddy viscosity. The advection term F u is given by where u α+1/2 is the velocity at the interface between neighbouring the αth and (1 + α)th layers, and G α+1/2 is a mass exchange term between the αth and (1 + α)th layers. This mass exchange term can be computed as and the interface velocity is computed by a simple upwinding following the sign of the mass exchange term as Finally the pressure term F p is given by where the term S α+1/2 can be seen as an apparent topography for the (1 + α)th layer since it is defined by Note that for a single layer Saint-Venant problem, the model (1) reduces to where H is still the water depth but u the water velocity of the whole flow.

Exner Equation
To update the bed-load in the multilayer system (1), we used the Exner equation given by where p is the sediment porosity assumed to be constant and the sediment discharge Q can be evaluated by the simple formula introduced in [10] with u 1 is the velocity of the lowest layer, m and A are coefficients usually obtained from experiments taking into account the grain diameter and the kinematic viscosity of the sediment. In practice, the values of the coefficient A are between 0 and 1 depending on the interaction between the sediment transport and the water flow. Another formula frequently used for the sediment discharge Q is given in [11] Q(u 1 , where the grain specific gravity s = ρs ρw , with ρ s is the sediment density. Note that most of existing formulations for sediment transport models are empirical to differing extents and have been derived from experiments and measured data. It should be stressed that the method described in this paper can be applied to other forms of sediment transport fluxes without major conceptual modifications. For instance, the bed-load sediment transport functions proposed in [8,9] can also be handled by the proposed multi-layer model. Notice that the parameters n b , d 50 , ρ s , p and A appeared in above equations are user-defined constants in the sediment transport model. In practice, the selection of these coefficients are problem dependent and their discussion is postponed for section 4 where numerical examples will be presented. In the following we refer as the coupled multilayer system when we consider the solution h, u α , z of the system of the N + 2 equations (1) and (12) and we refer as the coupled single-layer system when we consider the solution h, u, z of the system of 3 equations (11)- (12). We refer to [4,5] for a detailed analysis of the single-layer coupled system when the fluid model is the classical single-layer shallow water system (11). The theoretical analysis of the multilayer coupled system will be performed in a future work. We refer to [3] for an analysis of the multilayer system (1) when the bottom does not evolve in time.

Solution Procedure
In this section we describe the numerical method that is used to solve the coupled multilayer model. Two strategies can be used. The first one is referred as the quasisteady approach. The fluid model is first solved on a fixed bottom until a steady state is reached. Then the bottom is updated by using in the Exner equation (12)  In [3] a numerical strategy is presented for the discretization of the multilayer system (1). This strategy is mostly based on a kinetic interpretation of the model. The kinetic formulations were first introduced for other fluid models [7,1] but is particularly interesting in the multilayer context since it is a way to provide a stable numerical scheme without requiring the computation of the eigenvalues of the (exact or approximated) Jacobian matrix. The advection source term F u can be included in the kinetic interpretation and we use relation (7) to discretize it by using the kinetic fluxes. The real and apparent F p topographic source terms can be handled out by the use of an extension of the hydrostatic reconstruction procedure that was introduced in [2]. We use a simple implicit solver to deal with the friction source terms. In order to derive a numerical strategy for the coupled fluid-bottom problem we have to introduce a kinetic interpretation of the Exner equation (12). Let us first briefly recall the kinetic interpretation of the fluid model. We consider here the simple case of the classical single-layer shallow water system (11) but the ideas are the same for the multilayer case (1) and the details can be found in [3]. We introduce a kinetic velocity ξ and an even probability χ(ξ) with second momentum equal to unity. We define the related kinetic function M f We claim that (h, u) is solution of the shallow water model (11) if and only if M f is solution of the kinetic equation where Q f is a kernel such that its two first integrals in ξ vanish. This property becomes clear when considering the integration of the kinetic equation (16) against 1 and ξ since we recover the conservative part of the shallow water system (11). We now introduce a kinetic interpretation for the Exner equation (12). We define a new density function where δ denotes the Dirac measure (which is a particular case of even probability) and we claim that (h, u, z) is solution of the coupled water flow and bed-load model (11)- (12) if and only if (M f , M b ) is solution of the system of linear equations This kinetic interpretation is very helpful to design a simple and stable finite volume scheme for coupled multilayer system. Let us recall the general form of an explicit finite volume scheme for this system where subscript i and superscript n refer to the space discretization and time discretization respectively, with F i+1/2 the numerical flux whose definition characterizes the chosen finite volume method. In kinetic scheme this numerical flux is first designed at kinetic level for the linear system (18). The linearity makes the choice very natural and we simply apply an upwind scheme following the sign of the kinetic velocity ξ. Since the macroscopic system is obtained from the kinetic equations after a simple integration procedure, we can do the same thing for the numerical scheme and we do recover a numerical scheme for the macroscopic system by considering a simple integration of the kinetic upwind scheme. This strategy is quite powerful since it leads to a stable and simple numerical scheme. Stable means the numerical scheme preserves the positivity of the water height and simple means the kinetic interpretation does not explicitly appears in the numerical scheme since all the computations in the integration process can be perform analytically for reasonable choice of the probability χ, compare [1] among others. For example we obtain the following fluxes for the Exner equation We notice that this numerical flux satisfies the enforced consistancy property

Numerical Results
We present numerical results obtained for mutli-layer model when wind-driven circulations in a closed basin with movable bottom are considered. We focus on this study because it is a quite simple case that is out of the scope of classical monolayer shallow water system. Indeed let us note that the only stationnary state of the classical shallow water system in this configuration is a flow at rest situation where the slope of the free surface is in equilibrium with the wind force. This analytical result is clearly unphysical since we expect that a circulation will take place in the basin : the surface water will be driven by the wind while a reverse current will appear near the bottom.
We present a numerical result that is obtained for such a test case. The initial bottom topography is a parabolic bump. In Figure 4 we present that stationary circulation that takes place in the lake. In Figure 4 we present the initial and final topography of the basin.

Concluding Remarks
In this paper we present a first numerical solution to the morphodynamic problem when the evolution of the flow is modeled by a multilayer approach. The solution is in qualitative agreement with the expected In a near future we would like to investigate some other possible numerical procedure in order to solve this quite complex coupling problem. Moreover we also have to compare the numerical results E. Audusse, F. Benkhaldoun, S. Sari and M. Seaid