Asymptotic derivation and simulations of a non-local Exner model in large viscosity regime

. The present paper deals with the modeling and numerical approximation of bed load transport under the action of water. A new shallow water type model is derived from the stratified two-fluid Navier–Stokes equations. Its novelty lies in the magnitude of a viscosity term that leads to a momentum equation of elliptic type. The full model, sediment and water, verifies a dissipative energy balance for smooth solutions. The numerical resolution of the sediment layer is not trivial since the viscosity introduces a non-local term in the model. Adding a transport threshold makes the resolution even more challenging. A scheme based on a staggered discretization is proposed for the full model, sediment and water.


Introduction
We deal with the derivation and the numerical resolution of a model describing bed load transport. Bed load transport is the material transport occurring close to the surface of the sediment layer at the bottom of a river or of the sea. The main applications pertain to the prediction of the evolution of bed forms (such as dunes) and the evolution of river morphology. In an industrial context, the accumulation of sediment at the bottom of a hydroelectric dam or of a harbor are relevant problems.
For the simulation of rivers, models based on the shallow water equations coupled with the Exner equation are commonly used. These vertically integrated models are suitable for time scales ranging from months to years. The Exner equation is only a mass conservation equation and requires a closure relation for the solid flux. Solid flux formulae such the Grass formula [18] or the Meyer-Peter & Müller formula [32], are widely used but choosing a solid flux formula from the literature is not straightforward. Confrontations with real data seem to show that there is no universal solid flux formula. Most of the formulae only involve the friction at the interface but some of them, as the formula in [36], acknowledge the role of the slope of the sediment layer as a driving force. A transport threshold is frequently involved. It is well known that shallow water-Exner models have difficulty in reproducing dune growth. A phase shift between the maximum shear stress and the maximum of the bed topography is a necessary ingredient to obtain bed form amplification, see [23] and the later works [9,13,26]. This phase shift introduces some non-locality to the model. To obtain this phase shift, a possibility is to refine the water flow model by introducing a boundary layer such as it was proposed in the Triple Deck model [13,26,29] and the multilayer model [4], as well as the bilayer model in [21]. In the present work, we show that the non-local influence can be obtained by improving the description of the sediment layer. In [35] a bilayer shallow water model is proposed, which means that the sediment layer is described as a fluid with its own inertia. This model gave satisfactory results when the density of the sediment layer is close to that of the water layer, but did not perform well for heavier sediment layers. Another improved description was proposed in [12], in which a shallow water-Exner system is derived from a stratified bi-fluid Navier-Stokes equations. It follows that the solid flux depends not only on the water velocity but also on the slope of the free surface. However, none of the models in [12,35] introduce non-local terms. Another strategy for introducing a non-local flux term was proposed in [7,30] based on the minimisation of the energy of the waves.
The approach adopted here is similar to that in [12], but the scaling used for the derivation of the model in the sediment layer is different. In the case of a sediment layer with negligible viscosity, we recover the solid flux established in [12]. Considering a non-negligible viscosity in the sediment layer, a non-local solid flux is obtained. From a mathematical point of view, the model is close to some existing models in the literature, for example the Patlak-Keller-Segel equations [22,31], the Schurtz-Nicolaï model [16,34] and the Stokes-Brinkman model [19]. The works describing these models highlight the difficulty in analyzing them and solving them numerically.
Our main result consists in establishing a new shallow water type model for bed load transport, taking into account the viscosity of the sediment layer. More specifically, the equation for the velocity of the sediment layer reads as a nonlinear elliptic equation, responsible for the non-local character of the model. Even if only a Newtonian viscosity is considered, the model is already very complex to analyze. Considering a Coulomb friction at the substratum, we introduce a transport threshold thanks to a Coulomb friction law at the bottom that is often necessary for applications. A numerical scheme is proposed; it relies on a finite-volume, staggeredgrid discretization and satisfies a dissipative law for the discrete entropy. The scheme for the sediment layer is coupled to an existing scheme for the water layer [8,20], modeled by the shallow water equations, and the coupling strategy does not introduce parasitic entropy. Simulations of water flowing on a sediment dune are performed with and without viscosity. The significant impact of the viscosity is illustrated.
A quick description of the water-sediment system is given in Section 2. The model is derived in Section 3 and a preliminary analysis of the model is made. The numerical scheme is presented in Section 4.2, where 1D numerical results are also shown. In Section 5.1, the coupled water-sediment model is presented, the numerical resolution of which is explained in Section 5.2. The numerical results regarding the coupled system are proposed in Section 5.3.

Bilayer Navier-Stokes equations
We consider a Cartesian coordinate system where x ∈ R , ∈ {1, 2}, is the coordinate in the horizontal plane, ∈ R is the coordinate in the vertical direction and ∈ R + is the time coordinate. In this work, we focus on the modeling of the coupled water-sediment system in the regime of bed load transport. We assume that due to gravity, the flow is well-stratified, i.e. there exists an interface (x, ) that splits the flow in two parts, see Figure 1: below this interface and above a non-erodible substratum (x) from now on referred to as the sediment layer, only the sediment phase is present; above the interface and below the free surface (x, ) from now on referred to as the water layer, only the water phase is present. In addition, we assume that the flow in both layers can be modeled as a continuous medium and more precisely with the incompressible Navier-Stokes equations. Let us specify some notations. In each layer ∈ { , } ( stands for sediment whereas stands for water), the horizontal velocity is denoted by (x, , ) ∈ R , the vertical velocity by (x, , ) ∈ R and the pressure by (x, ). The governing equations in each layer read where the symmetric gradient is used D x = ∇x +(∇x ) 2 . The fluid ∈ { , } is characterized by its density ∈ R * + and its viscosity ∈ R * + fixed. The free surface and the water-sediment interface are respectively governed by the kinematic equations A no-penetration condition is assumed at the water-sediment interface and at the substratum For each surface ∈ { , , }, we define the normal by and a base of tangent vectors by The viscosity tensor is defined for ∈ { , } by and the stress tensor is assumed to be continuous at the free surface and at the water-sediment interface, i.e.
For the sake of simplicity, the atmospheric pressure was set to zero. The friction law at the water-sediment interface reads with the shear reads U ⋆ = U | = − U | = with the velocities U = ( , ) . At the substratum, the friction depends on the fluid, i.e. if the sediment layer vanishes or not. If the sediment layer vanishes, i.e. for any ( , ) such that ( , ) = ( ) we have Otherwise, for any ( , ) such that ( , ) > ( ) we have where the friction force at the substratum is further described in Section 2.2. The friction coefficients > 0, > 0 and the friction force at the substratum can be functions of the horizontal space variable, the surrounding pressure and the shear, i.e. and The system (2.1) must be completed with initial data and (x, 0) = 0 (x).

Introduction of a threshold for the onset of motion
Classical laws of solid sediment transport used in the context of hydraulic engineering have a threshold for incipient motion, and are a power function of the difference between a shear stress and a critical shear stress. The Meyer-Peter and Müller law [32] and van Rijn law [37] are examples of such laws involving a threshold. In the following, we show that a threshold for the onset of motion can also be introduced in the models (3.2) by means of the operator (︀ x, | = , U | = )︀ by defining it with a generalized Coulomb's friction law. The Coulomb's friction law, classically used for contacts between solids, claims that if the relative velocity between the solids does not vanish, the force exerted by each solid on the other one is on a cone (so that the normal and tangential components of the friction force are proportional) and its direction of the horizontal component is opposed to the relative velocity. When the relative velocity vanishes, the friction force is inside the cone. On the other hand, the friction law for a contact between a fluid and a solid is usually proportional to a power of the relative velocity, see for instance the Manning-Strickler law [28]. Since the sediment layer can be considered as a mixture that is neither really a fluid nor a solid, we propose that the friction is a combination of the two classical frictions, i.e.
The parameter (x) ≥ 0 is the Coulomb's friction coefficient, that goes to zero if the layer is almost fluid. The parameter (x) > 0 is the Strickler's coefficient, that goes to zero if the layer is almost solid. The parameter 0 < ≈ 1 depends on the regime of the flow, usually = 1 for a laminar flow and = 2 for a turbulent flow. The shear stress at the surface of the substratum is the resultant of the forces on a column of sediment except the friction at the substratum, i.e.
Even if the model (2.1) to describe the sediment-water system neglects several physical processes such as suspended sediment, water exchanged between the sediment layer and the water layer, as well as complex rheology of the sediment, it is already too complex to be simulated at the scale of a river or of a harbor. To carry out large-scale simulations, a reduced model is needed. Introducing reduced models for the sediment layer is the purpose of the next section.

Derivation of integrated models for the sediment layer
This section focuses on the sediment layer, without considering the water layer. More precisely, the action of water on sediments is solely due to water pressure and water velocity at the sediment-water interface, respectively referred to in this section by = | = and U = U | = . This is a preliminary step before considering the coupled water and sediment layers. Different models are derived for the sediment layer considering different parameters. Then, a numerical scheme is presented and validated with numerical experiments.
To derive reduced models using asymptotical arguments, we introduce the following characteristic values of the system =̃︀ , with = , and respectively the characteristic time and characteristic horizontal and vertical dimensions of the problem, , and respectively the characteristic values of the horizontal and vertical velocities and the pressure, , ∈ { , } the characteristic values of the friction parameters. The dimensionless velocity and the pressure in the water are also introduced, i.e.︀

=̃︀ = and̃︀ = ·
A vertical reference level 0 has been introduced. From now on, we assume that the surfaces are smooth enough, so that the variations of the surfaces are respectively characterized by The flow is therefore characterized by the following dimensionless numbers respectively named the numbers of shallowness, the Froude number, the Reynolds number, the density ratio and what we call the large scale Shields numbers at the substratum and at the sediment-water interface. The shear celerity ⋆ is the characteristic value of the velocity difference at the sediment-water interface The large scale Shields number at the sediment-water interface can be linked to the classical Shields number, , by introducing the characteristic grain diameter . More precisely Θ = 1− . In the following, several vertically-integrated models are formally derived depending on the scaling of the dimensionless numbers, in particular the Reynolds number and the large scale Shields numbers. The models can be summarized by the following general model̃︀̃︀ wherẽ︀ ( , ) correspond to the thickness of the sediment layer (in other words,̃︀ is an approximation of − ) and̃︀ ( , ) correspond to the horizontal velocity in the sediment layer (in other words,̃︀ is an approximation of̃︀ ). The dimensionless friction at the surface of the substratum reads︀̃︀ with the critical shear stress̃︀ = (︁̃︀ +̃︀̃︀ )︁ and the dimensionless shear stress reads︀ The parameters ( , , ) ∈ {0, 1} 3 depend on the scaling, see Propositions 3.1 and 3.2 and activate the following physical processes. The parameter corresponds to the scale of the ratio between the solid velocity and the forcing velocity . For = 0, the order of magnitude of the solid velocity is the same as that of the forcing velocity, while for = 1, the solid velocity is one order smaller than the forcing velocity. The parameter reflects the impact of the gravity on the motion of the sediment layer i.e. when = 0, the gravity term can be neglected whereas when = 1 it cannot. Similarly, the parameter reflects the impact of the viscosity of the sediment layer, i.e. when = 0, the viscosity term can be neglected whereas when = 1 it cannot.
The following models were derived assuming that = 1. The proofs and the models are still valid for lower Froude numbers. However, low Froude regimes are beyond the scope of this work. For details about the low Froude limit, see [24,25,38].

Models with local sediment discharge
Let us now introduce the models with large Reynolds number. We derive three models corresponding to three different asymptotic regimes. In all these models, the solid flux is local, i.e. it only depends of the velocities and the gradient of the free surface at the location where the flux is estimated. The model Proposition 3.1i corresponds to the Exner model with a Grass-type flux [18], where the solid velocity is proportional to the forcing velocity. The model Proposition 3.1ii corresponds to the stratified immiscible bilayer shallow water model. This model was previously used to describe the sediment transport in particular in [39]. In model Proposition 3.1iii the gravity effect remains and the solid velocity is one order smaller than in the other models in Proposition 3.1. The importance of the gravity term was first evidenced in [14,27]. A model very similar to model Proposition 3.1iii was derived in [12] assuming a different scaling, i.e. while in [12], the dimensionless horizontal velocity in the sediment was assumed to be of the order of 2 , we obtain the order of magnitude of︀ as a result of the asymptotic limit, and this order of magnitude is . Then the system (3.2) with the initial conditioñ︀ 5) and, in the case of model Proposition 3.1ii only, the additional initial conditioñ︀ and wherẽ︀̃︀ = 2 ,̃︀̃︀ = , is an approximation of the Navier-Stokes system for the sediment layer with negligible terms in order of (︀ 1+ )︀ .
Proof. First, the mass conservation is obtained classically by integration of the divergence free equation (the first equation of (2.1)). In addition, it yields that = and from the vertical momentum equation (the third equation of (2.1)), it classically yields that 2 = (1). Now, let us consider the main terms of the horizontal momentum equation, of the stress continuity conditions at the surface of the substratum and at the sediment-water interface. In the cases (ii) and (i), we write the following auxiliary problem  Integrating the horizontal momentum equation betweeñ︀ and̃︀ yields︀ This gives the result considering the scaling.
In the case (iii), we write the following auxiliary problem It yields that̃︀ = ( ) and̃︀ does not depend on up to the order 2 , i.e.︀ , integrating the vertical momentum equation leads to the hydrostatic relation (3.6) (up to the order 2 ). Next, we integrate the horizontal momentum equation betweeñ︀ and̃︀; we get Using the orders of magnitude of the Reynolds number = 1 and of the sediment-water interface large-scale Shields number Θ = allows to further simplify equation (3.9). Thus, the result is obtained.
Note that the scaling of the velocity (3.8) might be not satisfied by the initial conditioñ︀ 0 = 0 . Introducing the short timeˆ=̃︀ 2 , the previous auxiliary problem becomeŝ̃︀ thereforẽ︀ becomes of the order of within a characteristic time of the order of ( 2 ).

Models with non-local sediment discharge
The main drawback of the Exner models derived in Proposition 3.1 is that they do not take into account the viscosity in the sediment layer. To derive models involving an operator accounting for the viscosity, the product between the Reynolds number and the Shields number at the surface of the substratum must be of the order of . These models, presented in Proposition 3.2, are viscous versions of the one presented in Proposition 3.1. Then the system (3.2) with the initial condition (3.5) and wherẽ︀̃︀ = | = ,̃︀̃︀ = | = , is an approximation of the Navier-Stokes sediment layer with negligible terms in order of (︀ 1+ )︀ .
Proof. The proof is similar to that of Proposition 3.1. However the auxiliary problem coming from the horizontal momentum equation and the boundary condition reads As in the proof of Proposition 3.1, the initial condition which does not necessary satisfy this scaling vanishes within a characteristic time of the order of (︀ 2 )︀ . Integrating the vertical momentum equation allows to obtain that the pressure is hydrostatic (3.6). The horizontal momentum equation vertically integrated betweeñ︀ and︀ is again equation (3.7). The first two results (i) and (ii) are direct simplifications of (3.7) considering the scaling. For the scaling (iii), the main term of (3.7) reads̃︀̃︀ and it follows that ,0 = 0. To be totally exact, it depends on the boundary condition. However, assuming the velocity at the bound is small (in order of ) or neglecting a boundary layer of thickness in order of , the next term reads This concludes the proof.
In equation (3.2) and if = 1, the term̃︀̃︀̃︀ is of the order of the modeling error, thus can be removed in term of modeling approximation. However, in the perspective of the derivation of the coupled system, it is kept to ensure energy dissipation, see Proposition 5.1.

Governing equations
In the following, we will focus on model (3.2) in the case where ( , , ) = (1, 1, 1), i.e. Proposition 3.2iii. From the physical point of view, this model seems interesting because throughout its derivation, it appears that the solid velocity is much lower than the water velocity, which corresponds to observations, and it takes into account viscous effects in the sediment layer. From the mathematical point of view, this model is much more complex than it seems. Even if this model can naively be considered as simple as the model with inertia ( , , ) = (0, 1, 1), for which numerical strategies are already proposed in the literature [1,6,39], it is not the case. In particular it is well known that the numerical scheme for the case with inertia ( , , ) = (0, 1, 1) does not give satisfactory results in the regime where the inertia terms are negligible, see [5], and a first step is to propose a numerical scheme for the asymptotic model without inertia. Let us first rewrite it in the dimensional framework by multiplying the equation for the conservation of mass by HU and the momentum balance by HU 2 .
with the dimensional friction at the substratum.

Main properties of the model
One can think that in the case of vanishing velocity, i.e. = 0, the minimum function is not needed since the friction is always equal to the shear stress. However, without the minimum function, the solution = 0 (and thus = 0) is always possible. This solution is not consistent with physics and comes from the fact that we neglect the inertia of the sediment layer. The friction formula must be understood such that if for vanishing velocity the shear stress is smaller than the critical shear stress, then this solution is physically valid. It follows that the friction at the substratum can be reformulated as This new formulation (valid only in the regime without inertia) is interesting from a numerical point of view, see Section 4.2. Moreover, the solid velocity reads We recover a power law with a threshold similar to the laws used by hydraulic engineers, see [32,37]. The usual power used in practice is recovered for = 2 3 . However this formula can be misleading since the shear stress is also a function of the velocity.
Let us now highlight some physical properties of model (4.1) to improve the relevance of this model.
Assuming that is bounded, the right hand side vanishes by continuity of . We conclude by considering the initial condition.
Before stating the energy balance, we introduce the potential energy of a column of fluid ℰ of height ℎ placed upon a topography at elevation Proposition 4.2. For smooth enough solutions, the mechanical energy of (4.1) satisfies the following energy balance where the potential energy and the flux of energy in the sediment layer respectively read and the dissipation of energy reads Proof. Let us multiply the continuity equation of (4.1) by ( + ), while the equation on the velocity is multiplied by . Combining the two equations it leads to Then we multiply the continuity equation of (4.1) by but we add the substratum level to the time derivative. We get (︂ ( + ) The mechanical energy of the sediment layer is made only of its potential energy. Note that the estimate of Proposition 4.2 is a dissipation law in the sense that without forcing, i.e. = = 0, the total mechanical energy decreases, i.e. ∫︀ R ℰ dx ≤ 0, since the friction term is dissipative, i.e. · ≥ 0 and (D x ) : (D x ) ≥ 0 because the matrix D x is symmetric.

Asymptotic low-viscosity limit
In this section, the asymptotic limit of the model when the viscosity vanishes is identified. Let us denote with an underline · the solution of the model (4.1) without viscosity = 0. Considering the second equation of (4.1), we observe that the velocity is actually oriented in the same direction (and sense) as the shear stress The norm of the velocity is the unique positive root of a nonlinear polynomial. It reads On the other hand, the continuity equation can be written as an nonlinear advectiondiffusion equation, i.e.

Numerical scheme
In the present section, we propose a numerical strategy to solve (4.1) in one dimension. A staggered grid discretization is used to reduce the stencil of the scheme at the asymptotic limit of low viscosity, see Section 4.2.2. Such a staggered grid discretization was already used for the shallow water system in [20,33] and for the shallow water-Exner system in [10]. Let us consider a Cartesian grid of points )︀ with > 0 the constant space step. The numerical unknown +1 is the approximation of the sediment thickness averaged in a cell ] − 1 /2 , + 1 /2 [ at time +1 = + , where is an adaptative time step defined later on. In addition, + 1 /2 is an approximation of the solid velocity at the interface + 1 /2 and at time . For readability purposes, the following centered discrete operators are used with the number of cells and the number of interfaces of the grid. Let us set the forcing terms + 1 /2 and respectively defined by the values (or an approximation of) . The discrete friction coefficient at the substratum-sediment interface reads , The discrete friction coefficient at the interface depends on the velocities and and on the pressure at the interface , i.e. , Let us focus on the numerical resolution of the non-local model (4.1). We choose the following discretization for the continuity equation with the upwind reconstruction of the sediment layer thickness at the faceŝ +1 The velocity is computed thanks to the following discretization +1 with the friction at the surface of the substratum given by +1 The discrete critical shear stress is defined by , )︁ and the discrete shear stress is defined by +1 The last term in (4.6) is estimated at time +1 because at the limit of low viscosity, this term leads to the diffusion term in (4.3). More details are given in Section 4.2.2.
Upwinding is required to ensure the positivity and the entropy stability of the solution under an hyperbolic CFL condition, see Proposition 4.3. Moreover, the positivity of the thickness is ensured if the time step satisfies the following implicit CFL condition 2 +1 (4.7) The computation of +1 + 1 /2 is not obvious, since it depends on +1 which has yet to be computed. Equations (4.4) and (4.5) form a nonlinear system with + unknowns. This system can be reduced to a system with only unknowns by replacing +1 in (4.5) using the scheme (4.4), which gives if | +1 + 1 /2 | > , + 1 /2 , and +1 + 1 /2 = 0 otherwise. This strategy was already proposed in a simpler case in [16] and a multidimensional version [17], to approximate the solution of the Schurtz-Nicolaï model [34] in plasma physics. The system described by (4.8) is actually nonlinear, because the reconstructionˆ+ 1 + 1 /2 and the shear +1 + 1 /2 depend on +1 + 1 /2 . In practice, a fixed-point method is used. The presence of the threshold makes the system stiff, hence a Newton fixed-point method is used to increase the convergence rate. The velocity is initialized with ,0 + 1 /2 = 0. The convergence criterion used to estimate the convergence of the iterative process is based on the ∞ -norm of the variation of the solution , between two successive iterations. At each iteration, the shear stress , + 1 /2 is compared to the critical shear stress , , the corresponding line in the matrix and in the right-hand side of the system are modified. The diagonal coefficient of the th-line of the matrix is set to 1 while the others coefficients are set to 0 and the th-line of the right-hand side is set to 0. At each iteration, the time step is estimated to satisfy the CFL condition (4.7) at convergence using the relation with , + 1 /2 the velocity at the iteration . We briefly discuss here the implementation of several types of boundary conditions. The computational domain is made of cells. The boundary conditions on the sediment layer thickness are implemented by means of ghost cells numbered 0 and + 1. As regards the boundary conditions on the velocity , the values are imposed at the edges at the bound, i.e. 1 /2 and + 1 /2. The boundary conditions used in practice depend on the test case and will be specified for each test case.
It is not clear that the proposed fixed-point method converges, in particular because of the threshold. However, this work focuses on the modeling of the sediment transport while the numerical results are only illustrations of the behavior of the model. The analysis of the fixed-point method is out of the scope of this work.

Main properties of the scheme
In this section, the energy stability of the scheme is proven. It is the discrete counterpart of the dissipative law of Proposition 4.2.
Proposition 4.3. Assume that the initial condition is non-negative, i.e. 0 ≥ 0, and that the CFL condition (4.7) is satisfied. Then the solution of the scheme (4.4), (4.5) is non-negative, i.e.
≥ 0, and the discrete mechanical energy satisfies the following law with the discrete potential energy ℰ , = ℰ( , ) defined in (4.2), the discrete flux of energy reads

+1
, and the discrete dissipation of energy reads Assume that at the time iteration , the solution is non-negative. Under the CFL condition (4.7), the non-negativity of +1 is a classical property of the upwind scheme.
Let us now focus on the mechanical energy estimate. Equation (4.4) is multiplied by the potential (︀ +1 + +1 )︀ . Using the identity 2 ( − ) = 2 − 2 + | − | 2 (4.10) with = +1 and = , we get where the reconstruction of the potential at the face +1 is used. We also multiply equation . We get An expression for the mechanical work term is obtained by multiplying equation (4.5) by +1 which is then substituted in equation (4.12), i.e.
It remains to rewrite the last terms as a flux and a source term. Expanding this term and using the identity (4.10) allows to write 1 2 We obtain the announced result.

Asymptotic low-viscosity scheme
As presented in Section 4.1.2, the model when the viscosity tends to zero can be written as an advectiondiffusion equation (4.3). At the discrete level, the scheme (4.4) and (4.5) is formally equivalent to the scheme where the sediment layer thickness at the faces is obtained with an upwind reconstruction, i.e. +1 the discrete coefficients are defined by +1 with , + 1 /2 = ( , − ) and the discrete shear stress without velocity and the norm of the velocity is the unique positive root satisfying Even if the scheme (4.13) is not the simplest strategy to approximate the solution of (4.3), one can show that it is consistent using classical arguments. In addition, the stability of the scheme is a corollary of Proposition 4.3 applied with = 0. However the scheme is still implicit and nonlinear, and a fixed-point method is needed for the computation. Here the two strategies are significantly different. In the case of the asymptotic scheme (4.13), the sediment layer thickness is initialized with the value at the previous time step ,0 = , then the shear stress , 0, + 1 /2 is estimated using the approximation of , instead of +1 , then the norm of the velocity ⃒ ⃒ ⃒ +1 Once again this problem is nonlinear and in practice, we simply set Then we compute the coefficients , +1 + 1 /2 and , +1 Finally, the new approximation , +1 is obtained using a linearization of (4.13), i.e. we replace respectively

Numerical validation
In this section, the behavior of the numerical scheme described in Section 4.2 is illustrated. For all the test cases, except when indicated otherwise, the physical parameters are set to = 9.81, = 0.6, = 1, = 0.5, = 0, , + 1 /2 = 1 and , + 1 /2 = 10 −3 . The length of the domain is 1 and it is discretized with a space step set to = 10 −3 . The time step is variable, determined using the CFL condition (4.9) with = 1. Again, except when indicated otherwise, Neumann boundary conditions are imposed, i.e. for any ≥ 0, 0 = 1 and +1 = .

Synthetic forcing
The convergence of the scheme presented in Section 4.2 towards an analytical solution is studied. Though deriving analytical solution is not trivial because of the nonlinearity, one may recover an analytical solution by imposing the adequate forcing. More precisely, fixing = 0, we look for the forcing to recover a given and steady thickness = ( ). The continuity equation implies that the velocity reads = , with a given constant in time and space. Setting ̸ = 0 implies that the shear stress is larger than the threshold value and the second equation of (4.1) gives an expression for the forcing velocity This strategy is applied to assess the convergence of the numerical scheme. We set ( ) = 1 + 0.1 sin (2 ) , = 1. and ( ) = 0.
The initial condition is set to 0 = 1 and Dirichlet boundary conditions are imposed, i.e. for any ≥ 0 we set 0 = 1 − 0.1 sin ( ) and +1 = 1 + 0.1 sin In Figure 2, the errors in 2 -norm are plotted at = 20 for several values of . The time is long enough for the numerical solution to reach the stationary regime. As expected because of the use of an upwind reconstruction in the continuity equation, the convergence order of the scheme is 1.

Asymptotic behavior
As explained in Section 4.2.2, the scheme (4.4), (4.5) converges to the scheme (4.13) as goes to 0. In this section, the convergence of (4.4), (4.5) towards (4.13) with respect to is illustrated numerically. We look for a simulation modeling a sediment layer without water above, so we set = 0, = 0 and coherently, no friction at the interface, i.e.  (4.13) at the time = 2 × 10 −3 . The final time is chosen such that the sediment bump is not entirely flat at even for the lower viscosity values, which allows to quantify the differences between the solutions. To prevent the error due to the time discretization from degrading the comparison between the two schemes, a time step = 10 −5 is imposed. With this value for the time step, the CFL condition (4.7) is always satisfied. The results are shown on Figure 3. As goes to 0, the solution of the scheme (4.4), (4.5) converges towards the solution of the local scheme with order 1.

Influence of the viscosity on the shape of the solution
To discuss on the influence of the viscosity on the shape of the solution, several simulations were performed with various values of the viscosity . In Figure 4a, the initial condition (4.14) is used with = 1 and = 1.1 and the solution is plotted at time = 10 −3 . We remark that the initial discontinuity of the solution is preserved when the viscosity is large enough, in this case > 10 −3 and more generally the viscosity has a large impact on the speed of evolution of the solution. This observation is in agreement with the asymptotic analysis seen in Section 3. In Figure 4b, the initial condition (4.14) is used with = 0 and = 0.5 and the solution is plotted at time = 10 −1 . These simulations are performed to illustrate the ability of the scheme to deal with dry areas. Because of the nonlinearity, the solution is stiff at dry fronts regardless the value of the viscosity.

Non-flat steady state
This section is devoted to the illustration of the necessity of a threshold to obtain a non-flat steady solution. The value of the critical shear stress is set to = 1 and a forcing is introduced through the bottom shape given by ( ) = 0.5 − 0.1 .
The initial condition is described by (4.14) with = 0.1 and = 0.2. Two different simulations are run, one with = 0 and the other with = 0.5. The two simulations are run with a constant time step = 10 −6 , which is small enough for the CFL condition to be satisfied in both cases. This value is also small enough for the numerical scheme not to "miss" the stationary state. The final states are shown on Figure 5. Both simulation reach the non-flat stationary state characterized by the angles of repose which can be determined from the critical shear stress . However, reaching the final state takes longer when the viscosity is larger. The transients of the numerical solutions are not the same either, as observed in Section 4.3.3. Yet the stationary state does not seem to depend on the viscosity. In Figure 6, the shear stress of the steady solution is plotted. With  = 0.5, the shear stress is equal to the critical shear stress (for the points which have moved during the simulation), while for = 0, it is slightly below. Remark that the slope of the steady solution of the solution without viscosity is slightly smaller than the one with velocity. Actually because of the larger velocity in the simulation without viscosity, the numerical scheme misses the exact steady solution and stops slightly after the threshold. This observation is even clearer with a larger time step. The implementation of the threshold raises many theoretical and numerical difficulties which are beyond the scope of the present work. For a discussion of these problems, see for instance [2]. The numerical simulations are only performed to illustrate the relevance of the model (4.1).

Coupled water and sediment system
We are now interested in the modeling and simulation of the whole system, made of a water layer and a sediment layer.

Modeling of the coupled system
Let us introduce the coupled model made of in the water layer and (4.1) in the sediment layer with the forcing now determined by the water dynamics, i.e. = and = ℎ. In order to derive the coupled system, the dimensionless numbers introduced in (3.1) were again considered. Additionally, the Reynolds number in the water layer 2 was introduced. The coupled model is derived from the Navier-Stokes equations with the terms of order of ( ) neglected in the water layer following the derivation of the shallow water equations in the water layer done in [15] and the derivation of the models for the sediment layer from the Navier-Stokes equations done in Proposition 3.2iii (Props. 3.1ii, 3.1iii and 3.2ii can also be considered). The initial condition of the water layer reads The scalings of the Propositions 3.1i and 3.2i can not be considered since the larger friction at the interface Θ is not consistent with the derivation presented in [15].
Remark that the orders of approximations are different in the two layers. However, the fact that the order of approximation in the water layer is does not prevent the approximation from being of the order of 2 in the sediment layer. Indeed, in the second equation of (2.1), is multiplied by , which is of the order of .
where the potential energy, the kinetic energy and the flux of energy in the water layer respectively read The sediment energy and flux as well as the dissipation term are defined in Proposition 4.2.
Proof. Remark that the energy in Proposition 4.2 can be split in two terms. The first term ℰ corresponds to the potential energy of the sediments and the second term ( + ) is the energy due to an external pressure on the sediment layer. Considering the coupled system, hence := ℎ, this term correspond to the part due to the bathymetry in the potential energy in the water layer, i.e.
Let us recall that the shallow water model satisfies the following energy balance This result can be recovered multiplying the continuity equation (first equation) of (5.1) by the potential ℎ and the momentum equation (second equation) of (5.1) by the velocity then summing. Eventually, the right-hand side can be written using the continuity equation as Multiplying (5.2) by the density ratio and summing with the energy balance of the sediment layer given in Proposition 4.2 we get the result.

Numerical scheme for the coupled system
In this section, we present a numerical scheme for the coupled model (4.1) and (5.1). We look for a scheme in the water layer which can be coupled to the scheme (4.4) and (4.5) for the sediment layer. It is important for the friction terms that the two velocities and are discretized at the same location. Since the scheme (4.4) and (4.5) is designed on a staggered grid, we look for a scheme on a staggered grid for the water layer. As for the sediment layer, let the numerical unknown ℎ be the approximation of the water depth ℎ averaged in a cell ] − 1 /2 , + 1 /2 [ at time and + 1 /2 be an approximation of the water velocity at the interface + 1 /2 and at time . The scheme described in [8,20] and used for the simulation of the shallow water-Exner system in [10] meets the design requirements presented at the beginning of the current section. The discrete continuity equation reads with the upwind mass flux This discretization is very close to the discretization used in the sediment layer (4.4), except that the velocity is taken at time and thus the water depth ℎ +1 can be computed explicitly. Once the discrete water depth ℎ +1 is known, the velocity is computed using the scheme ℎ +1 with the centered reconstruction of the water depth at the faces and the approximation of the velocity neglecting the friction ℎ +1 with the upwind momentum flux Neglecting the friction, i.e. considering ⋆ + 1 /2 given by (5.5) as an approximation of the velocity at time +1 , the scheme (5.3) and (5.5) satisfies a discrete counterpart of (5.2).
Lemma 5.2. The scheme (5.3) and (5.5) satisfies the following energy balance where the potential energy ℰ 0, = ℰ (ℎ , 0) is defined in (4.2) and the kinetic energies read The quadratic reconstruction of the water depth at the faces is usedh and the remainder reads Proof. As it was done in [8], we start by estimating the evolution of the kinematic and potential energies separately. Following Proposition 3 of [8], we get ⋆ and Proposition 4 of [8] that can be written under the form We conclude with a straightforward combination.
The remainder ℛ tends to vanish when the time step and the space step go to zero. Unfortunately, Lemma 5.2 is not a proof of stability of the shallow water scheme since the remainder is not signed. However, in practice the numerical scheme seems stable under the CFL condition Let us now focus on the scheme of the coupled model. The following discrete counterpart of Proposition 5.1 can be shown.
Proof. Let us remark that as in the continuous framework we have ℰ 0, = ℰ 0, + ( + ). By summing the estimate of Proposition 4.3 and Lemma 5.2, we will recover a balance of the total energy. In addition, the forcing of the pressure at the interface in Proposition 4.3 can be written using the continuity of the water layer as On the other hand, multiplying (5.4) by +1 + 1 /2 and using the identity (4.10) we get +1 We conclude with a straightforward combination.
Remark that the coupling between the water layer and the sediment layer does not introduce a additional remainder. As for the shallow water scheme (5.3) and (5.5), there is still a remainder not signed. However, we expect the coupled scheme (4.4)-(4.5)-(5.3)-(5.4)-(5.5) to be stable under the CFL conditions (4.7) and (5.6).
The numerical computation of the solution is still not trivial because of the coupling through the time step and the friction term step (5.4). Let us now give some details. Knowing the numerical approximation at time , we compute and a first estimation of the time step. Then we start the iterative strategy to compute the sediment layer presented in Section 4.2. More precisely, we compute an approximation of the water depth ℎ , +1 using (5.3) and of the velocity without friction ⋆, +1 using (5.5). Then using (5.4), we have an explicit formulae of +1 as function of ⋆ , +1 and +1 + 1 /2 , i.e. +1 + 1 /2 = , with , Unfortunately, +1 + 1 /2 as well as +1 + 1 /2 are still unknown. However, the explicit formulae (5.7) is introduced in (4.8) in order to compute , +1 + 1 /2 . It yields , , + 1 /2 1 , Once the velocity , +1 + 1 /2 is known, the sediment thickness , +1 + 1 /2 is computed using (4.4). The iterative process is computed while , +1 , the ℓ ∞ -norm of the residual of the nonlinear scheme (4.5) defined by is larger than a given tolerance > 0. This computation strategy is illustrated in Algorithm 1 where the given numerical parameters are the error tolerance of the iterative process > 0 and the CFL parameter 0 < ≤ 1. Remark that the water depth and velocity are computed in the iterative loop only because of the possible change of time step. In practice, the CFL condition of the water flow (5.6) is more restrictive than the CFL condition coming from the sediment flow (4.9) since the water velocity is higher than the sediment velocity. At the software level, a Boolean can be added to compute the water quantities only if the time step was changed inside the iterative loop. The computational complexity of the numerical scheme for the coupled system is finally of the same order as the sediment scheme alone, even if several explicit calculations are added.

Numerical validation for the coupled system
In what follows, the behavior of the coupled system is illustrated. The length of the domain is = 10. The surface of the substratum is ( ) = 0.5. Unless specified otherwise, the left boundary condition in the water is (ℎ )(0, ) = 0.5 and the right boundary condition is ℎ( , ) = 0.5. In the sediment layer, the boundary , +1 ← (4.9); , +1 = min with the amplitude given in each test case. These boundary conditions and initial sediment shape are those described in [10]. Again, = 1. The friction coefficient at the water-sediment interface is also constant in space. Its value is = 2 × 10 −3 . As in Section 4.3, = 9.81, = 0.6, = 1 and there is no threshold for the onset of motion, i.e. = 0.
To initialize the water layer, we run the simulation with a fixed sediment layer until the water reaches the steady state.

Dune growth test
In this section, we study the influence of the viscosity term on the evolution of a sediment dune under an initially subcritical water flow, i.e. (5.10) with = 1, illustrated in Figure 7a. With the mentioned boundary and initial conditions, the initial water flow is subcritical. The difference with existing models such as the Grass model and that in [12] is shown. Two different studies are performed: one with = 0, and one with = 0.5. A convergence study is performed with = 0. The sediment profiles obtained at = 10 for an increasing number of cells are shown on Figure 8. For the converged solutions, the dune steepens and grows slightly. The sediment accumulates on the downstream side of the bump and a shock seems to appear in the sediment layer at the beginning of the solution. Yet, there is no hydraulic jump in the water layer, the flow remains subcritical everywhere in the domain. The growth of the dune is due to a transient occurring when the sediment layer becomes deformable. Indeed, during the initialization phase for the water, the sediment layer is rigid and the water layer does not transmit energy to it. When the sediment layer becomes erodible, it suddenly receives energy from the water layer. Then, after a short time, the height of the sediment bump decreases due to the diffusion process. This phenomenon is difficult to catch numerically. Large values of are necessary to observe the initial dune growth.
The influence of the viscosity on the time evolution of the dune is assessed. The viscosity is = 0.5. Then, a convergence study is performed with = 0.5. The convergence is harder to obtain than with = 0, in the sense that a larger number of cells is needed to correctly approximate the solution. The converged solutions are shown on Figure 9a and have the following behavior. When the simulation starts, the dune begins to grow. Then, a hydraulic jump appears in the water layer and the dune sharpens, evolving  into a peak under the hydraulic jump and a smaller bump following it. The sediment peak moves downstream along with the hydraulic jump and the phenomenon maintains itself. The sediment velocity is much lower than that of the water; it is the deformation in the sediment layer that moves as fast as the hydraulic jump. A zoom on the hydraulic jump and on the sediment underneath reveals that the sediment profile brutally changes under the shock in the water layer, see Figure 9b. The profiles obtained in the sediment layer are very sharp. But they do not result from a shock in the sediment layer: multiple cells are involved in the peak, see Figure 9b. Figure 10 shows the sediment profiles obtained as the mesh is refined. For the coarsest meshes used, i.e. for ≤ 400, the numerical solution fails to capture the correct behavior: no hydraulic jump is created. The dune grows a bit on the downstream side and then is eroded. Neither the peak height nor the peak velocity are correct   Figure 11. This error is given by the formula The solid flux for the local model is plotted on Figure 12b. Due to the presence of the slope term, the monotonicity of the flux is very different with respect to that of the Grass flux. Yet the local maximum of the flux coincides with the peak of the sediment bump. The additional effect of the viscosity term is illustrated on Figure 12c. The local maximum of the flux is shifted with respect to the peak of the sediment bump, it is located upstream. From the plots of the fluxes, it is not easy to predict the evolution of the dunes. Yet, the fact that they are very different from the plot of the Grass flux -and from each other -provides some insight about why the behaviors subsequently observed are complex.
Numerous authors, among which the authors of [13,26], have argued that such a shift between the shear stress and the peak of the sediment bump is necessary for the dune to grow. In these works, the solid flux directly depends on the shear stress , but in the local and non-local models, the solid flux is , which is why we are interested in the shift of with respect to the sediment bump. While in [13,26] the shift between the maximum solid flux and the sediment bump was achieved via an improved description of the water layer, in the present work, it is obtained thanks to the addition of the viscosity term.

Sensitivity to the boundary condition and initial condition
In this section, we assess the sensitivity of the behavior of the sediment layer with respect to the discharge imposed in the water on the left boundary and the initial condition. The viscosity is set to = 0.5 and the numerical solutions are computed with = 25 600. The sensitivity with respect to the discharge is assessed first. Figure 13a shows the result at time = 18 for two different water discharges at the left boundary condition starting from the same initial condition, i.e. propagate. A small change in the inflow boundary condition can induce a very different behavior, which means that the evolution of the sediment layer is very sensitive to the inflow boundary condition. Then, the behavior of the sediment dune when the steady-state solution in the water is transcritical is illustrated with a viscosity set to = 0.5. The initial shape of the dune is now given by (5.10) with = 1.2. More precisely, the sediment dune is initially higher and narrower, though both dunes have the same mass. When the higher initial shape, the stationary solution in the water presents a hydraulic jump, see Figure 7b. Figure 13b shows the behaviors of the dunes initially described by (5.10) with = 1 (solid line, with the caption "reference geometry") and with = 1.2 (dashed line, with the caption "modified geometry"). At the beginning of the simulations, the two dunes have very different profiles since a peak forms only for the higher initial shape. When the time advances, a peak also forms for the initial condition with = 1 and moreover the two peaks almost overlap after a long enough time.

Sensitivity to the friction and viscosity coefficients
In this section, we investigate the different behaviors that can be obtained with the coupled model as the physical parameters vary. For simplicity, we only study the influences of the friction between the sediment and water layers and of the viscosity in the sediment layer . We perform simulations with a range of coefficients varying from 0 to 10 −1 and varying from 0 to 10 1 . The other parameters are defined as in Section 5.3, the initial condition is given by (5.10) with = 1 and the numerical results were observed until time = 10. All the simulations are performed with = 20 000. We observe several behaviors of the solution that we try to classify. The regions in the ( , ) plane where we observe the different behaviors are illustrated in Figure 14. The different behaviors observed are plotted in Figure 15 with the following plotting style to illustrate the time evolution until time = 10: the lighter the line, the older the solution. Let us now describe each case more precisely.
We start our description by the case without viscosity = 0 and low friction ≪ 1, i.e. the mark in Figure 14. An exemple of this behavior is plotted in Figure 15a. We characterize this behavior by a step down in the sediment layer, following the flow direction. This step is located at the same abscissa as the hydraulic jump in the water layer but exhibits opposite variations, i.e. the sediment depth increases where the water depth decreases. This case does not correspond to the classical Exner regime but to the model proposed in [12] since the main term governing the sediment dynamics in this case is the gravity term. Note that similarly to the shallow water equations, the physical relevance of discontinuous solutions is not clear, since discontinuous solutions are out of the scope of the asymptotic derivation of the sediment layer model, see Proposition 3.2.
When we increase the friction coefficient ≥ 5 × 10 −2 , we observed a new behavior, marked by in Figure 14. An exemple of this behavior is plotted in Figure 15b. We characterize this behavior by the fact that the dune spreads more than it propagates. For a longer simulation, it seems clear that the dune will disappear completely in this case.
Then if we increase the viscosity coefficient ≥ 5, we observe a new behavior, marked by in Figure 14. An exemple of this behavior was plotted in Figure 15c. We characterize this behavior by the fact that the dune propagates with almost no spreading over the duration of the simulation. However, we did not observe any amplification of the amplitude of the dune, hence this regime does not appear to be responsible for the formation of the characteristic sediment transport patterns. Now by decreasing the friction coefficient, we observe that the evolution of the solution is slowed down, until < 5 × 10 −3 , marked by in Figure 14, where the solution is almost unchanged at the final time. This behavior can be explained by the fact that we reduce the energy transferred from the water to the sediments. Therefore, we need to reduce the viscosity to facilitate the deformation of the sediment layer.
By decreasing the viscosity coefficient ≤ 5, we observe a new behavior, marked by in Figure 14. An exemple of this behavior is plotted in Figure 15d. We characterize this behavior by the fact that the dune does not spread but becomes sharper with time. It is obvious that this behavior is a transition and we believe that for a longer simulation, the dune will finally evolve in a similar way to that for smaller viscosity values, i.e. behavior . By decreasing even more the viscosity coefficient ≤ 1, we observe a new behavior, marked by in Figure 14. An example of this behavior is plotted in Figure 15e. We characterize this behavior by the fact that the dune splits into two parts, i.e. a new local maximum is observed. This behavior was already described in the previous simulations, in particular Section 5.3.1. The two parts of the dune have a radically opposite behavior. More precisely the behavior of the upstream part can be compared to the case , i.e. it spreads more than it propagates. The behavior of the downstream part becomes very steep, especially on the upstream side. We explain this formation by the contribution of material by the upstream part.
Finally, when the viscosity coefficient is too small ≤ 10 −1 (except when it vanishes), marked by in Figure 14, the fixed point of the numerical method described in Section 4.2 does not converge anymore. A more robust iterative numerical strategy seems needed in this case, unless the root cause of the problem is actually that the model is ill-posed in this regime. We also observe a region marked by in Figure 14, where the fixed point converged but the solution presents instabilities. An exemple of this behavior is plotted in Figure 15f. We have not been able to determine if these instabilities are related to the model or if they are numerical artifices.
To summarize this description, we can describe more generally that the cases and with large enough friction present a smooth surface of the sediment layer. For small frictions, sharp solutions are observed. Without viscosity, the model recovers the classical behavior that has limited capabilities for the simulation of characteristic sediment transport patterns such as dunes. With viscosity, we observe new behaviors and that are still to be compared to experiments for their physical relevance. There remains a set of parameters and that seems out of reach for our model or our numerical strategy.

Conclusions and perspectives
In this paper, a new bed load transport model is derived through an asymptotic analysis. This new model includes the modeling of the sediment viscosity and the Coulomb friction law at the substratum. From a mathematical point of view, the viscosity is modeled with a non-local operator, which makes the resolution non-trivial, especially with the transport threshold due to the Coulomb friction law. An entropy-dissipative numerical scheme is designed to solve this model and a coupling with the water layer is realized. In particular, the simulation of dune growth and propagation is performed to illustrate the behavior of the model.
The present work raises several questions. From the theoretical point of view, a more thorough analysis of the non-local system could be done. The regularity of the solutions has not yet been investigated, while it is crucial to obtain the positivity of the sediment depth. At the numerical level, two critical points need to be investigated for applications: the stability of the scheme in the low-viscosity, low friction regime; and the efficiency of the scheme which seems in the proposed version too costly to be applied at the scale of a river. The physical relevance of the results is beyond the scope of the present work, the purpose of which was chiefly to investigate the effect of the viscosity term. The simulations of the coupled system could be confronted with experimental data. This would require fitting the parameters , and . Moreover, in this work, the parameters are constant for the simulations, but one could try to determine more complex empirical laws such as a bottom friction term of the Chézy or Manning type for . Several matters are still to be addressed in the modeling. For example, more realistic rheologies, such as the Drucker-Prager rheology, can be considered. Moreover, in the description we have proposed, no mass exchanges occur between the sediment and water layers. The suspension and deposition of grains is also a relevant question. In [3], a layerwise-discretized model of the water for density stratified flows is presented. In [11], a multilayer model for polydisperse sedimentation is proposed. Including mass exchanges between the sediment layer and the water and adopting a layerwise-discretized approach to simulate density variations in the water is a challenging objective.