Sediment dynamics. Part 1. Bed-load transport by laminar shearing flows

We propose a two-phase model having a Newtonian rheology for the fluid phase and friction for the particle phase to describe bed-load transport in the laminar viscous regime. We have applied this continuum model to sediment transport by viscous shearing flows. The equations are shown to reduce to the momentum equation for the mixture and the Brinkman equation for the fluid velocity. This modelling is able to provide a description of the flow of the mobile granular layer. At some distance from threshold of particle motion where the continuum approach is more realistic as the mobile layer is larger than one particle diameter, there is very little slip between the two phases and the velocities inside the mobile bed have approximately a parabolic profile. When the Poiseuille (or Couette) flow is not significantly perturbed, simple analytical results of the particle flux varying cubically with the Shields number and of the bed-load thickness varying linearly with it can then be obtained. These predictions compares favourably with experimental observations of bed-load transport in pipe flows.

Introduction

When beds constituted of sediment particles are submitted to shearing flows, the particles at the surface of the stream-bed are able to move as soon as hydrodynamic forces acting on them exceed a fraction of their apparent weight. Bed-load refers to the sediment in transport that is carried by intermittent contact with the stream-bed by rolling, sliding, and bouncing. This situation occurs in a wide variety of natural phenomena, such as sediment transport in rivers, and in industrial processes, such as problems due to hydrate formation in pipeline that are encountered in oil production and granular transport in food or pharmaceutical industries. Despite the amount of work devoted to the subject for more than a century, a complete description of the mechanics of bed-load transport is still lacking as the coupling between the granular medium and the fluid flow is poorly understood.

The common way of representing the incipient motion of the particles is to use a dimensionless number called the Shields number, θ, which is constructed as the ratio of the shear stress at the top of the bed to the apparent weight of a single particle and which exceeds a critical value at threshold of motion, θ c . Most of the studies usually attempt to provide this critical Shields number, θ c , as well as to relate the bed-load flow-rate, q p , to the excess Shields number, (θ -θ c ).

The oldest and still widely used model to determine θ c is to write a balance of forces on the grains at the top of the sheared bed, see e.g. [START_REF] White | The equilibrium of grains on the bed of a steam[END_REF] and [START_REF] Vanoni | Sediment Transportation Mechanics: Initiation of Motion[END_REF]. The threshold value is found to be proportional to the tangent of the angle of repose of the

Authors

Method qp q ρ f ∆ρgd 3 [START_REF] Einstein | Formulas for the transportation of bed load[END_REF] empirical e 0.391 1 θ q 2 3 + 36η 2 gd 3 ρ f ∆ρ -q 36η 2 gd 3 ρ f ∆ρ 0.465 [START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF] empirical 8 (c2c3θ -θ c ) 1.5

Einstein (1950) semi-empirical Table 1. Various expressions for qp in turbulent flow, where c2 is a side-wall correction, c3 a bed-form correction, fc the percentage of grains of a given size put into motion, f the correction function obtained experimentally for grain size dispersion, ψ the drag coefficient, ρp the density of the solid, ρ f the density of the fluid (∆ρ = ρpρ f ), η the viscosity of the fluid, d the particle diameter, and tan α the dynamic friction coefficient. Table 2. Various expressions for qp in laminar flow, where Re * = u * dρ f /η is the shear Reynolds number, u * the shear velocity, and N the number of particles in motion per area.

1 -1 √ π R 0.156f 1 θ -2 -0.156f 1 θ -2 e -t 2 dt 27fc " 2 -1 √ π R 0.156f 1 θ -2 -0.156f 1 θ -2 e -t

Authors

grains. Most of the existing experimental data concerns the turbulent regime and present large scatter which is due to the difficulty of defining as well as measuring the quantities of interest such as the shear rate and incipient motion, see [START_REF] Buffington | A systematic analysis of eight decades of incipient motion studies, with special reference to gravel-bedded rivers[END_REF], [START_REF] Dancey | Probability of individual grain movement and threshold condition[END_REF], [START_REF] Paintal | Concept of critical shear stress in loose boundary open channels[END_REF], [START_REF] Vanoni | Sediment Transportation Mechanics: Initiation of Motion[END_REF], and [START_REF] White | The equilibrium of grains on the bed of a steam[END_REF]. Recent experiments in laminar flows have provided reproducible measurements for the threshold of motion and inferred a critical Shields number (θ c ≈ 0.12) independent of the Reynolds number, see [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF]. Note that this threshold is also independent of the bed preparation as it corresponds to a saturated state of the surface microstructure of the bed. It indeed differs from the lower threshold corresponding to the initial motion of single particles and above which the particle flux eventually vanishes due to local trapping and reorganization of the particles near the bed surface, see [START_REF] Charru | Erosion and depositionof particles on a bed sheared by a viscous flow[END_REF].

Several expressions have been proposed for the bed-load flow-rate, q p . Some of them are presented in tables 1 and 2. Again, most of the empirical or semi-empirical laws have been obtained in the turbulent regime. Notably, the same inertial scaling, ∆ρgd 3 /ρ f , for q p found by dimensional analysis is used in these expressions. A viscous scaling, ∆ρgd 3 /η, is generally obtained for laminar flows. Of particular relevance to the present laminar context are the studies of [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF], [START_REF] Charru | Erosion and depositionof particles on a bed sheared by a viscous flow[END_REF], and [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid. Part 1[END_REF] which address different situations. [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF] used the viscous resuspension model of [START_REF] Leighton | Viscous resuspension[END_REF] in which the settled particles are carried entirely in suspension by a competition between shear-induced diffusion and sedimentation. In this situation, the particles interact through hydrodynamic forces and are prevented from touching owing to lubrication forces. Note that one of the first attempts to tackle resuspension in the laminar regime was offered by [START_REF] Bagnold | The flow of cohesionless grains in fluids[END_REF] in part III of his paper. [START_REF] Bagnold | The flow of cohesionless grains in fluids[END_REF] used a viscous rheology for the granular stress which diverges at close packing and examined the case of an infinite flow depth exhibiting a broad concentration profile and that of a limited flow depth having a constant concentration. Conversely, [START_REF] Charru | Erosion and depositionof particles on a bed sheared by a viscous flow[END_REF] and [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid. Part 1[END_REF] dealt with a completely different situation. They developed a phenomenological model for the erosion and deposition of particles rolling and sliding over a fixed bed which only considered a mobile particle monolayer.

In the present paper, we propose a continuum approach describing bed-load transport in the laminar viscous regime. We address a situation in which the bed-load can be considered as a mobile granular medium where the particles mainly interact through contact forces. This continuum description of the particulate or granular phase stems from recent developments made in dense suspension and granular flows, see [START_REF] Nott | Pressure-driven flow of suspensions: simulation and theory[END_REF][START_REF] Nott | Pressure-driven flow of suspensions: simulation and theory[END_REF][START_REF] Gdr Midi | On dense granular flows[END_REF]. We present a two-phase model having a Newtonian rheology for the fluid phase and friction for the particle phase. A description of this model is given in § 2 and is meant to provide a general framework for future threedimensional study. This continuum model is then applied to viscous shearing flows in § 3. The equations are shown to reduce to the momentum equation for the mixture and the Brinkman equation for the fluid velocity. Calculations of bed-load transport by shearing flows in closed two-dimensional channels are performed numerically but also analytically at some distance from threshold where the mobile layer is larger than a grain size. These different calculations allow us to discuss in which conditions the continuum modelling is sensible and realistic. In § 5, the analytical predictions are compared to experimental results obtained for a bed composed of spherical particles in laminar pipe flow, see the description of the experiments in § 4. The closures used in the two-phase model are discussed in § 6 where conclusions are also drawn.

Two-phase model

Many authors have derived two-phase governing equations to describe a system of particles and fluid in an average sense. There is an extensive literature around the process of averaging the equations of motion, for a review see e.g. [START_REF] Jackson | The dynamics of fluidized particles[END_REF]. Each type of averaging is purely a formal process and should lead to essentially the same results if properly done. It generates averaged quantities more numerous than the available equations and therefore there is a closure problem which is a central issue in the modelling. There is no guarantee though that such a closure is possible. Here, we choose to use the standard equations derived by [START_REF] Anderson | A fluid mechanical description of fluidized beds -Equation of motion[END_REF] and revisited by [START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid[END_REF] and propose some closures appropriate to the present problem.

Two-phase equations

The equations of continuity for the fluid and the particle phases are respectively

∂ǫ ∂t + ∂(ǫu f i ) ∂x i = 0, (2.1) ∂φ ∂t + ∂(φu p i ) ∂x i = 0, (2.2)
where u f i is the local mean fluid velocity, u p i the local mean particle velocity, φ the particle volume fraction, and ǫ = 1 -φ the void fraction or fraction of space occupied by the fluid.

Adding equations (2.1) and (2.2), the mixture (particles + fluid) is found to be incompressible

∂U i ∂x i = 0, (2.3)
where U i = φ u p i + ǫ u f i is the volume average velocity. The momentum equations for the fluid and particle phases are respectively

ρ f D f (ǫu f i ) Dt = ρ f ∂(ǫu f i ) ∂t + ∂(ǫu f i u f j ) ∂x j = ∂σ f ij ∂x j -nf i + ǫρ f g i , (2.4) ρ p D p (φu p i ) Dt = ρ p ∂(φu p i ) ∂t + ∂(φu p i u p j ) ∂x j = ∂σ p ij ∂x j + nf i + φρ p g i , (2.5)
where g i is the specific gravity force vector, ρ f the fluid density, ρ p the particle density, n the number density (number of particles per unit volume). The force f i represents the average value of the resultant force exerted by the fluid on a particle. The stress tensors σ f ij and σ p ij may be regarded as effective stress tensors associated with the fluid and particle phases, respectively.

Closures

Following [START_REF] Jackson | The dynamics of fluidized particles[END_REF], the local average force exerted by the fluid on a particle nf i is decomposed into a generalised buoyancy force and a force which gathers all the remaining contributions

nf i = φ ∂σ f ij ∂x j + nf 1 i .
(2.6)

In the case considered here of a viscous flow through a particle bed, the remaining contribution reduces to the dominant viscous drag caused by the relative motion. Extending the well-known empirical law of Darcy, nf 1 i can be written as

nf 1 i = η ǫ 2 K (u f i -u p i ) = η ǫ K (U i -u p i ), (2.7)
where η is the viscosity of the pure fluid. The coefficient of permeability K is empirically related to ǫ and the particle diameter d by the Carman-Kozeny relation

K = ǫ 3 d 2 k (1 -ǫ) 2 .
(2.8)

For k = 180, equation (2.8) is known as the widely-used Carman equation, see e.g. [START_REF] Happel | Low Reynolds number hydrodynamics with special applications to particulate media[END_REF]. We chose to use this relation as it agrees well with the experiments of [START_REF] Goharzadeh | Transition layer thickness at a fluid-porous interface[END_REF] in packed beds and the lattice-Boltzmann simulations of [START_REF] Koch | Particle pressure and marginal stability limits for homogeneous monodisperse gas fluidized bed: kinetic theory and numerical simulations[END_REF]), Hill, Koch & Ladd (2001a), Hill, Koch & Ladd (2001b), [START_REF] Van Der Hoef | Lattice-Boltzmann simulations of low-Reynolds-number flow past mono-and bidisperse arrays of spheres: results for the permeability and drag force[END_REF] for large volume fraction (φ > 0.4) and small particle Reynolds numbers. We assume that the effective stress tensor associated with the fluid phase σ f ij is of Newtonian form and is related to the volume average velocity by the relation

σ f ij = -p f δ ij + τ f ij = -p f δ ij + η e ∂U i ∂x j + ∂U j ∂x i , (2.9)
where η e is the effective viscosity of the mixture. The choice of this Newtonian form for σ f ij can be justified for dilute suspensions of Stokesian particles by exact closure calculations done by [START_REF] Zhang | Momentum and energy equations for disperse two-phase flows and their closure for dilute suspensions[END_REF] and [START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid[END_REF]. In that limit, the familiar Einstein correction to the viscosity of the pure fluid, η e = η(1 + 5φ/2), is recovered. For non dilute suspensions, one needs to rely on empirical relations expressing η e /η as a function of φ, see e.g. [START_REF] Stickel | Fluid mechanics and rheology of dense suspensions[END_REF]. Following [START_REF] Bagnold | The flow of cohesionless grains in fluids[END_REF] [his formula (58)] and [START_REF] Goharzadeh | Transition layer thickness at a fluid-porous interface[END_REF], we choose for simplicity the Einstein viscosity in § 3 and will discuss other empirical laws in § 6.

The stress tensor of the particle phase σ p ij comprises terms coming from particle contact interactions and particle Reynolds stress. Clearly, for a bed of particles, close contacts are dominant. The simplest way to describe this stress coming from direct particle-particle interactions,

σ p ij = -p p δ ij + τ p ij , (2.10)
is to use a Coulomb friction model where: (a) the tangential stress is proportional to the load when the granular shear rate is positive, i.e. is equal to µp p in two-dimensional shearing flows, with a friction coefficient µ which mostly depends upon the particle geometry and which is given by the tangent of the angle of repose, (b) the tangential stress is indeterminate when the granular shear rate is zero. This simple Coulomb model will be used in § 3. In § 6, we will also consider a constitutive law which has been first derived for dry granular flows and has been then found successful for submarine granular flows, see GDR Midi (2004), Jop, Forterre & Pouliquen (2006[START_REF] Cassar | Submarine granular flows down inclined planes[END_REF], and [START_REF] Pouliquen | Flow of dense granular material: towards simple constitutive laws[END_REF]. 3. Calculation of bed-load transport by shearing flows

Formulation of the problem

We consider a flat particle bed of thickness h p submitted to a Poiseuille or a Couette flow in a two dimensional channel of thickness D, see figure 1. The flow is considered stationary and uniform. It is also parallel and the velocities U i , u f i , u p i reduce to their x components denoted U , u f and u p . Moreover, the x-invariance leads to ∂p p /∂x = 0. This also gives ∂p f /∂x = 0 in the case of a Couette flow whereas ∂p f /∂x is a given constant in the case of a Poiseuille flow. With these assumptions, continuity equations (2.1), (2.2), and (2.3) show that U , u f and u p are sole functions of y.

We have the choice of length scale and time scale and since we do not know a priori what would be the correct scalings, we decide to make all the values dimensionless by scaling the length by the channel thickness D and the pressure by the hydrostatic pressure ∆ρgD, and therefore the time by η/∆ρgD where ∆ρ = ρ p -ρ f . Using these scales and the proposed closures (2.7), (2.9), and (2.10), we obtain the dimensionless momentum equations along the x direction

0 = -ǫ ∂ pf ∂ x + ǫ ∂ τ f xy ∂ ȳ - ǫ 2 D 2 K (ū f -ūp ), (3.1) 0 = ∂ τ p xy ∂ ȳ -φ ∂ pf ∂ x + φ ∂ τ f xy ∂ ȳ + ǫ 2 D 2 K (ū f -ūp ), (3.2)
and along the y direction

∂ pf ∂ ȳ = - ρ f ∆ρ , (3.3) ∂ pp ∂ ȳ = -φ. (3.4)
These latter equations show that the pressure of the fluid phase along gravity is simply the hydrostatic pressure and that the pressure of the particle phase is proportional to the apparent weight of the solid phase and increases when penetrating inside the bed. By expressing ūf in terms of the volume average velocity Ū , equation (3.1) becomes the Brinkman equation for the fluid phase, see Brinkman (1947a) and Brinkman (1947b),

- ∂ pf ∂ x = - ∂ τ f xy ∂ ȳ + D 2 K ( Ū -ūp ), (3.5)
where the pressure gradient mostly balances the dominant Darcy term. By summing equation (3.1) and(3.2), we find the mixture (fluid + particles) momentum equation

0 = ∂ τ p xy ∂ ȳ - ∂ pf ∂ x + ∂ τ f xy ∂ ȳ , (3.6)
which describes the exchange between the stress of the fluid phase, τ f xy = (η e /η)(d Ū /dȳ), and that of the solid phase, τ p xy = µp p when dū p /dȳ > 0. The above two equations will be used to solve the present problem.

In the case of a Poiseuille flow, the kinematic boundary conditions are ūf (1) = Ū (1) = 0 and ūf (0) = ūp (0) = Ū (0) = 0.

(3.7)

In the case of a Couette flow, only the first condition is changed to ūf (1) = Ū (1) = Ū0 . At the top of the bed, ȳ = hp , we also have τ p xy ( hp ) = 0 and pp ( hp ) = 0, (3.8) as there is no particle phase for ȳ hp . Using the first condition of formulae (3.8), equation (3.6) can be integrated to give

τ p xy (ȳ) = - ∂ pf ∂ x ( hp -ȳ) + τ f xy ( hp ) -τ f xy (ȳ). (3.9)
At the top of the granular bed, ȳ = hp , only a tangential fluid stress τ f xy ( hp ) is present. Equation (3.9) just means that this tangential fluid stress above the bed should balance the combined fluid + particle stresses inside the bed. The fluid shear stress is equal to τ f xy ( hp ) at ȳ = hp and goes to zero at ȳ = hc below which the granular medium is immobile and behaves as a porous medium. The particle shear stress is zero at ȳ = hp and increases inside the bed until it reaches µp p where the granular medium starts to be sheared. This means that there is a thin layer of granular medium above this depth which is moving but is not sheared. The particle shear stress keeps the value µp p until ȳ = hc inside the bed. Indeed, while the confinement pressure pp increases inside the bed, equation (3.9) shows that the particle stress is bounded. Therefore, the particle stress is eventually unable to match µp p and thus becomes indeterminate. The Coulomb friction model implies then that equation (3.6) must be substituted with

dū p dȳ = 0, (3.10)
for ȳ hc .

The flow can be separated into four regions which are depicted in figure 1. Region (I) corresponds to the pure fluid above the bed ( hp ȳ 1). Region (IV) is the region of the bed where the particles are at rest (0 ȳ hc ). Regions (II) and (III) are those where both the fluid and particles are in motion. Region (III) corresponds to a granular medium in motion with a constant porosity ǫ 0 = 1 -φ 0 ( hc ȳ hpd). The depth-variation of the porosity in this region is very weak and thus considered as negligible. Region (II) is the transition layer between the pure fluid region (I) and the pure Darcy region (III). The thickness of this region has been observed to be of the size of a grain diameter d and to have a porosity which varies linearly with thickness as ǫ(ȳ) = 1 + φ 0 (ȳ -hp )/ d ( hpd ȳ hp ), see [START_REF] Goharzadeh | Transition layer thickness at a fluid-porous interface[END_REF]. The problem is solved both numerically and analytically for the case of a Poiseuille flow in the following subsections and the main results are summarised in appendix A for the case of a Couette flow.

Numerical calculation

As indicated in § 2.2, we use the Einstein viscosity for η e and the Carman-Kozeny relation (2.8) with k = 180 for the coefficient of permeability K. The only unknown parameters (which can be measured independently) are then φ 0 and µ. The system of linear ordinarydifferential-equations (3.5) and (3.9) or (3.10) together with the boundary conditions (3.7) and the integrated pp deduced from (3.4) have been solved using a finite difference method with a non uniform grid which is chosen to be more refined near the interfaces.

The difficulty is to find the different domain of validity of equations (3.9) and (3.10) inside the bed. First, we suppose that equation (3.9) is valid everywhere in the granular medium. The fluid and solid velocity profile can be computed as shown in figure 2 (a) and dū p /dȳ can be found to be negative in some locations. Secondly, equation (3.9) is replaced by equations (3.10) at these locations and the velocity profiles are again calculated as seen in figure 2 (b). This process is iterated until dū p /dȳ 0 everywhere inside the bed, see figure 2 (c).

Figure 3 shows the velocity profile of the fluid and solid phases for three different Shields number. In the pure fluid zone ( hp ȳ 1), the profile is approaching that of a Poiseuille flow (a slightly asymmetric Poiseuille flow having a slip velocity at the top of the bed). A blow-up close to the bed interface shows that for small Shield number, θ = 0.048 (a), the solid phase is not moving at all. Owing to the porosity variation in region (II), the fluid velocity adjusts to the Darcy flow across the transition layer (≈ one particle diameter) as observed experimentally by [START_REF] Goharzadeh | Transition layer thickness at a fluid-porous interface[END_REF]. Note that if one skips region (II), the matching would occur on the scale of the bulk permeability, ≈ √ K, which is much smaller than the particle size. For slightly larger Shields number, θ = 0.094 (b), a very thin layer of solid phase (less than one particle diameter) is put in motion near the interface. There is a slip between the solid and the fluid phases. Note also that a thin layer of granular medium below the top of the bed is moving but is not sheared as in this region the particle stress has not yet reached µp p as explained at the end of § 3.1. For large enough Shields number, θ = 0.64 (c), the mobile layer grows thicker than a grain size and the slip between the fluid and the solid phases eventually becomes negligible as the Darcy term is dominant owing to the weak permeability in region (III).

These last numerical observations reveal the existence a critical Shields number θ c ≈ 0.06 for onset of motion of the solid phase which is a factor two smaller than the ex- perimental value θ c ≈ 0.12, see [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF]. However, the correct order of magnitude is recovered showing that the present continuum modelling having a variable porosity in region (II) is sensible. Note that having a variable porosity in region (II) is important close to incipient motion as skipping it leads to an unrealistic θ c ≈ 0.01 which is significantly smaller. Note also that just above onset of motion less than a monolayer

Region Ū ūp pp (I) ∂ pf ∂ x ȳ2 -1 2 + » - ∂ pf ∂ x hc + µφ0( hp -hc - d 2 ) - (ȳ -1) undefined 0 (II) - η ηe µφ0 (ȳ -hp + d) 3 6 d + η ηe " µφ0 + ∂ pf ∂ x « (ȳ -hc) 2 2 Ū φ0 2 d (ȳ -hp) 2 (III) η ηe " µφ0 + ∂ pf ∂ x « `ȳ -hc ´2 2 Ū φ0( hp -ȳ -d/2) (IV) 0 0 φ0( hp -ȳ -d/2)
Table 3. Volume averaged velocity, particle velocity and pressure.

hp

-hc ηe η (1 -hp) 8 > < > : -1 + 2 6 41 - -µφ0 " η ηe " 2 d2 3(1-hp ) 2 -η ηe µφ0 d 1-hp + η ηe ∂ pf ∂ x ∂ pf ∂ x + µφ0 3 7 5 1/2 9 > = > ; qp η ηe φ0 " µφ0 + ∂ pf ∂ x « » ( hp -hc) 4 -( hp -hc -d) 4 24 d - - η ηe µφ 2 0 d3 120 θ ( hp -hc) d " ∂ pf ∂ x + µφ0 « - µφ0 2
Table 4. Bed-load thickness, particle flow rate, and Shields number.

is mobile. This questions the validity of a continuum modelling close to incipient motion. If one defines the onset of motion as corresponding to a monolayer in motion, the simulation gives θ c ≈ 0.12 in excellent agreement with experiments.

Analytical calculation

For simplicity, we suppose that η e = η e (φ 0 ) and K = K(φ 0 ) are constant in region (II) and identical to their respective values in regions (III), and (IV). However, we choose to keep for now the variation of φ inside region (II) to obtain the pressure of the particle phase through equation (3.4). Indeed, this variation divides the apparent weight of the first monolayer by two which affects the whole granular medium pressure and then the solid friction. Equations (3.4), (3.5), (3.9) or (3.10) are solved to determine pp (ȳ), Ū (ȳ), and ūp (ȳ) in each region by matching Ū and pp at the different interfaces and using the boundary conditions (3.7) in the limit K(φ 0 ) ≪ h 2 , see table 3. Inside the bed, both particle and fluid phases move at the velocity of the mixture as the Darcy drag term is dominant in that limit. In other words, there is very little slip between the two phases, i.e. ūp ≈ ūf ≈ Ū in regions (II), (III), and (IV). Figure 3 shows the good agreement between the numerical and analytical profiles for a large enough Shields number [case (c)]. In that latter case, the assumption K(φ 0 ) ≪ D 2 is verified inside the bed except very close to the interface where the numerical profile of the particle velocity becomes more sharpened as a result of the use of a variable permeability coefficient. Also, at the top of the bed (approximately at half a grain size from the top), the fluid analytical and numerical profiles do not coincide perfectly as a variable effective viscosity is used in the numerical calculation. For smaller Shields numbers [cases (a) and (b)], the numerical simulation shows a slip between the two phases that the analytical calculation cannot capture within its hypotheses. This calculation also yields the bed-load thickness, hphc , the particle and fluid flowrates, qp and qf , and the Shields number, θ, as a function of the bed-load thickness, hphc , see table 4. Interestingly, the bed-load thickness varies linearly with the Shields number. The calculation cannot give the onset of motion that is observed numerically owing to its previously mentioned hypotheses. Nonetheless, considering that the critical Shields number, θ c , is given by the value of θ for hphc = d, i.e. corresponding to a monolayer in motion, yields to

θ c = ∂ pf ∂ x + µ φ 0 2 ≈ µ φ 0 2 , (3.11)
as ∂ pf /∂ x = O( d) at incipient motion. We recover that this threshold value is proportional to the tangent of the angle of repose. Another interesting finding is that it is also proportional to the particle volume fraction in the bulk of the bed. Taking φ 0 = 0.55 and µ = 0.43 as deduced by [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF] for glass spherical particles yields θ c = 0.12 in good agreement with the experiments.

In the limit

Re 2D ≪ (µ s φ 0 /12) (D -h p ) 3 /d 3 Ga with Ga = ρ f ∆ρgd 3 /η 2 and Re 2D = ρ f q f /η or equivalently θ ≪ (µφ 0 /2)(D -h p )/d, one finds θ = 6 d 2 D 2 Re 2D Ga D D -h p 2 , (3.12)
showing that the perturbation induced to the Poiseuille flow by the motion of the granular medium is negligible. In other words, this Shields number corresponds to the case of a Poiseuille flow comprising a flat solid bed of height h p , see [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF].

The bed-load flow-rate given in table 4, q p , can be further simplified to give

q p / ∆ρgd 3 η e = φ 0 θ c 12 θ 2θ c θ 2 θ c 2 + 1 - 1 5 , (3.13)
which has a viscous scaling built on the effective viscosity η e (φ 0 ). The dimensionless particle flux is proportional to the particle volume fraction in the bulk of the bed and to a function of θ c and θ and not simply of the excess Shields number θ -θ c . The different particle flow-rates are plotted versus Shields numbers in figure 4 for h p = 0.5D with φ 0 = 0.55 and µ = 0.43 (yielding to θ c = 0.12). The numerical flow-rate and the analytical approximation given in table 4 are in excellent agreement while the analytical approximation given by equation (3.13) notably deviates for θ 1.5 [in other words, the agreement is good in the range of validity of equations (3.12) and (3.13].

A simpler calculation

For θ ≫ θ c , equation (3.13) leads to a very simple expression

q p / ∆ρgd 3 η e = φ 0 θ c 24 θ θ c 3 , (3.14)
having the particle flux varying cubically with the Shields number. As shown in figure 4, it can be used for 0.5 θ 1.5. If one skips region (II) in which the volume fraction varies from φ 0 to 0 on one particle radius, one easily recovers this simpler expression for the particle flux as velocities inside the moving bed have the parabolic profile of region (III) in table 3, i.e.

U = u p ≈ (µφ 0 ∆ρg + ∂p f /∂x) η e (y -h c ) 2 2 ≈ µφ 0 ∆ρg η e (y -h c ) 2 2 , (3.15)
and, using the results in table 4, the bed-load thickness has a very simple linear variation in θ

h p -h c ≈ θd µφ 0 + ∂p f ∂x /∆ρg ≈ θd µφ 0 = θd 2θ c . (3.16)
Note that within this approximation, the critical Shields number θ c = µφ 0 /2 corresponds to the thickness of the mobile layer being half a particle size since the volume fraction in the top particle layer is φ 0 while, when the transition region (II) is kept, it corresponds to the thickness of the mobile layer being one particle size since the volume fraction is on average φ 0 /2 in this top layer. Note also that the velocity profile within this further simplification is slightly underestimated (dashed line in figure 3) as there is a small shift of the granular pressure due to the removal of region (II).

Clearly, since the thickness of region (II) is one particle radius, the continuum model cannot be expected to be sensible on such small length-scale. Even though it is able to provide a threshold for solid motion in the numerical calculation of § 3.2, it is likely to be only very qualitative just above incipient motion as less than a monolayer of solid phase is mobile. Clearly, the model is expected to be valid at some distance from onset of motion when the mobile layer is larger than a grain size.

Experimental measurements of bed-load transport in pipe flows

Particles and fluids

Three different batches of spheres (polystyrene particles supplied by Maxi-Blast, PMMA particles by Lehmann & Voss & Co., and glass particles by Potters-Ballotini) were used to perform the experiments. The particle size distributions were measured with a digital imaging system. The particle diameter distributions were observed to be approximately Gaussian for all the different batches and therefore well represented by a mean diameter d indicated in table 5 (the error corresponds to one standard deviation). The particle

Experimental apparatus

The apparatus test section was a horizontal glass tube of inner diameter D = 3 cm. The tube had a length L = 1.8 m which is longer than the entry length L e ≈ 60 cm necessary for the laminar flow to fully develop at Re pipe ≈ 300, where Re pipe = 4Q pipe ρ f /πDη is the tube Reynolds number and Q pipe the flow rate. All the experiments were performed for Re pipe < 300. First, the tube was filled with fluid, a low flow rate was imposed, and a mixture of particles and fluid was carefully introduced using a funnel to build an uniform flat bed. Secondly, the pipe flow was driven by gravity using continuous overflow from an overhead tank the elevation of which was varied. At the outlet from the test section, the particles were captured by a mesh while the fluid was run into a thermostated fluid reservoir. From this lower reservoir, the fluid was continuously returned to the overflowing reservoir by a Moineau progressing cavity pump (PCM model MR2.6H24). This arrangement isolated the test section from the pump and insured a constant temperature T = 35 • C across the whole experimental loop. Note that the captured particles were not re-injected into the test section. The flow rate, Q pipe , was determined by measuring a collected volume of fluid at the outlet of the tube in a given time. Note that, even though the flow is controlled by the relative height of two reservoirs, the high pressure losses in the remainder of the flow system result in an imposed flow-rate in the test section despite the changes in the tube section as particles are lost.

Bed profile measurement

For a combination of fluid and particles, the evolution of the bed height was recorded as a function of time for a given flow rate Q pipe and starting with an initial height of the bed, h start p . The bed was illuminated by a laser sheet positioned perpendicularly to its surface and aligned with the tube length in its middle at a distance ≈ 80 cm from the entrance of the tube. The illuminated upper layer of particles intersecting the sheet (of length ≈ 5 cm) was imaged by a digital camera. The images were then analysed (with ImageJ available at http://rsb.info.nih.gov/ij/) to yield the position of the fluidparticle interface. Each image was thresholded to turn this interface into a white curve which was further eroded to a single-pixel-thick curve. After calibration, this provides a precise measurement of the fluid-particle interface. In order to perform the calibration, a grid was inserted into the tube filled with pure fluid. An image of this grid was then recorded under the same optical conditions used in interface-position measurements and the co-ordinates of its points were measured.

Bed profile evolution

We choose to be in the condition of bed-load transport, i.e. above the critical Shields number. Particles at the surface of the bed are then always set in motion by the fluid flow and are observed to roll and slide. The bed-load thickness is of the order of a few particle layers. The temporal evolution of the bed height (averaged over the 5 cm of the measurement zone which is starting at ≈ 80 cm from the entrance of the tube) measured at the middle of the pipe is shown in figure 5. The bed height is always seen to decrease with increasing time as the test section is not fed in with particles. As expected, the bed height decrease is larger when increasing the flow rate and the viscosity for a given batch of particles and when decreasing the particle density for a given fluid. Note that, as the experiments are performed in a closed pipe, the Shields number evolves with the bed height. When the experiment is run for a long enough time (from one day to one week), the bed evolution eventually stops at a final bed height, h end p . This threshold for cessation of motion coincides with that for incipient motion as, by increasing the flow rate by a small amount, particles are set again in motion. Precise measurements of this critical bed height have provided the determination of the critical Shields number θ c = 0.12 ± 0.03 in the same viscous laminar flows, see [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF]. In some cases, the tube can be completely emptied if the Shields number for the empty tube is larger than the critical Shields number.

Comparison with experiments

We now compare the simple analytical predictions of § 3.4, in particular that for the particle flux (3.14), to the present experimental data in pipe flow. Note that these predictions are valid at some distance from threshold (typically θ 0.5) and hold for both Poiseuille and Couette flow (see appendix A).

First, we consider the critical Shields number for having a monolayer of particles in motion, θ c = µφ 0 /2 [equation (3.11)]. It is proportional to the tangent of the angle of repose which has been evaluated to be µ = 0.43 by [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF] for spherical particles and to the particle volume fraction φ 0 in the bulk of the bed. This last quantity has not been measured experimentally but can be estimated to be in the range 0.55 φ 0 0.65, see e.g. [START_REF] Goharzadeh | Transition layer thickness at a fluid-porous interface[END_REF]. This gives a value 0.13 ± 0.01 in good agreement with the experimental measurements 0.12 ± 0.03 obtained in the same experimental conditions of viscous laminar flows by [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF]. We should emphasise that this critical Shields for having a monolayer of particles in motion, θ c = µφ 0 /2, is independent of the specific shear flow used as it is identical for Poiseuille and Couette flows, see § 3 and appendix A. This is also confirmed by experiments as θ c Temporal evolution of the bed height, h p , for (i) batch A in fluid 2 with Q pipe = 6.45 10 -6 m 3 .s -1 (•) and in fluid 4 with Q pipe = 5.21 10 -6 m 3 .s -1 ( ), (ii) for batch B in fluid 1 with Q pipe = 3.37 10 -6 m 3 .s -1 ( ), in fluid 5 with Q pipe = 6.98 10 -6 m 3 .s -1 (△), in fluid 5 with Q pipe = 1.04 10 -5 m 3 .s -1 (⋄), and in fluid 3 with Q pipe = 2.56 10 -6 m 3 .s -1 (▽), and (iii) for batch C in fluid 3 with Q pipe = 2.13 10 -6 m 3 .s -1 (+) and with Q pipe = 5.16 10 -6 m 3 .s -1 (×). ( ), for batch B in fluid 1 with Q pipe = 3.37 10 -6 m 3 .s -1 ( ), for batch B in fluid 5 with Q pipe = 6.98 10 -6 m 3 .s -1 (△), for batch B in fluid 5 with Q pipe = 1.04 10 -5 m 3 .s -1 (⋄), for batch B in fluid 3 with Q pipe = 2.56 10 -6 m 3 .s -1 (▽), for batch C in fluid 3 with Q pipe = 2.13 10 -6 m 3 .s -1 (+), and for batch C with Q pipe = 5.16 10 -6 m 3 .s -1 (×). The solid line corresponds to the numerical integration of equation ( 5.2) for θ start pipe = 0.5, 1.8, and 7.5, respectively. has been found to be ≈ 0.12 in the Couette apparatus of [START_REF] Charru | Erosion and depositionof particles on a bed sheared by a viscous flow[END_REF] and in the Hele-Shaw cell of [START_REF] Loiseleux | Onset of erosion and avalanches for an inclined granular bed sheared by a continuous laminar flow[END_REF].

Secondly, we move to the prediction for the particle flux q p varying cubically with the Shields number [equation (3.14)]. We do not have a direct experimental measurement of the particle flux in the present pipe experiment but have the experimental evolution of the bed height which can be related to the particle flux by the particle flux conservation equation obtained by integrating the continuity equation (2.2) over D φ 0 ∂h p ∂t + ∂q p ∂x = 0.

(5.1)

Using equations (3.12) and (3.14)), the mass conservation (5.1) can be written as a kinematic wave equation in dimensionless form

∂θ ∂t * + c * (θ) ∂θ ∂x * = 0, (5.2)
exhibiting now the Shields number θ instead of h p and q p . The length-scale is again D but the time-scale differs from that (= η/∆ρgD) used in § 3 and has been chosen to be = (6Re 2D /Ga) 1/2 (D/d)(η e /∆ρgd). The dimensionless kinematic velocity

c * (θ) = θ c3/2 4 θ θ c 7/2 , (5.3)
describes the propagation of a perturbation in Shields number which is related to a perturbation in bed height.

To test this latter expression, we need to determine experimentally the kinematic velocity as well as the Shields number. First, we evaluate the velocity of the waves which are triggered at the entrance of the tube as L/t p where t p is the arrival time of the initial perturbation. This time is difficult to estimate as it is the starting time of the decrease in bed height in figure 5. We have instead measured the time for which the initial height has decreased by 0.5 % at the measurement location, which corresponds to t p /α with the numerical factor α = 0.87. Secondly, we deduce the initial height h start p from the time-evolution curves of figure 5. The Shields number is now estimated to be

θ pipe = 3 2 d 2 D 2 πβRe pipe Ga D D -h p 2 , (5.4)
which substitutes to equation (3.12) in the case of a pipe flow with Re pipe = 4ρ f Q pipe /πηD and with the numerical coefficient β = 1.85 in the limit 0.2 1 -hp 0.8 as found by [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF]. The switch to this new formula for the Shields number is due to the switch from a shear rate γ2D = 6(

Q 2D /D 2 )[D/(D -h p )] 2 in a plane flow to a shear rate γpipe = 6β(Q pipe /D 3 )[D/(D -h p )
] 2 in a truncated pipe flow where the coefficient β is found from numerical analysis.

In figure 6, the initial dimensionless kinematic velocity is plotted versus the initial Shields number defined by equation (5.4) with h p = h start p . The data collapse onto a master curve, which suggests that the time-scale found theoretically is appropriate. Equation (5.3) exhibiting the 7/2 power law [with µ = 0.43, φ 0 = 0.55, and η e = η(1 + 5φ 0 /2)] shows good agreement with the data with some departure for θ 1.5 as expected from the range of validity of the particle flow-rate (3.14), 0.5 θ 1.5 (note that θ 0.5 in most of the experiments). We have also tested these data against the law [q p ∝ (θ -θ c ) 3 ] proposed by [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF], see table 2. Importing this law into equation (5.2) while using the same length and time scales gives the new dimensionless kinematic velocity c * (θ) = 2.52 η e η θ 3/2 (θ -θ c ) 2 .

(5.5)

The experimental data are in agreement with this equation for large θ as the same power law 7/2 is recovered. As mentioned in the introduction, the law proposed by [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF] has been adapted from the viscous resuspension model of [START_REF] Leighton | Viscous resuspension[END_REF] which predicted q p ∝ θ 3 . Although this latter relation for the flow rate has the same scaling as relation (3.14), the physical underlying mechanisms are fundamentally different as there is no resuspension of the particles within the present modelling and in the experiments. We also provide comparison with the law [q p ∝ θ(θθ c )] derived by [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid. Part 1[END_REF] which was meant to be applied only to a mobile particle monolayer, see table 2. Again, importing this law into equation (5.2) gives the alternative kinematic velocity

c * (θ) = πc u c e c d 54 φ 0 η e η θ 3/2 θ c 2 (2θ -θ c ) , (5.6)
where, following [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid. Part 1[END_REF] inspired by the experimental measurement of individual particles in the top moving sediment layers by [START_REF] Charru | Erosion and depositionof particles on a bed sheared by a viscous flow[END_REF], the coefficients are chosen to be c u = 3.3, c e = 0.055, c d = 0.067, and θ c = 0.12 and the flow-rate has been adapted to the present two-dimensional geometry by multiplying the flow-rate of [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid. Part 1[END_REF] (expressed in unit of velocity per bed-area) by the volume of a particle. The experimental data are larger than this prediction at least by a factor 100 but they can be satisfactorily fitted by the power law 5/2 suggested by equation (5.6). It is important to mention that different calibrations of the coefficients than those originally . q p /(∆ρgd 3 /η e ) versus θ for the experimental data of [START_REF] Lobkovsky | Erosion of granular bed driven by laminar fluid flow[END_REF] for Q pipe = 500 (•), 750 (△), 1000 ( ), 1500 (⊳), 2000 (⋄), and 2500 (⊲) cm 3 min -1 and equation (3.16) (dashed line) taking the experimental values ∆ρ = 1.474 g/cm 3 , d = 700 µm, and the Einstein viscosity η e = η(1 + 5φ 0 /2) with η = 24.9 cP as well as the adjustable parameter φ 0 /24θ c2 ≈ 0.025 using θ c = 0.12 and φ 0 = 0.55.

given in [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF] or [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid. Part 1[END_REF] would produce lines having the same slope.

More direct comparison is given in figure 7 where the data for the evolution of the bed height of figure 5 are made dimensionless by using the Shields number (5.4) and the dimensionless time. Starting from different initial Shields numbers, the data eventually collapse onto a decay curve. This behavior is very well predicted by the numerical integration of equation (5.2) using a classical Euler scheme and taking the same values for µ, φ 0 and η e . The boundary condition corresponds to zero flux at the entrance as in the experiments. The initial condition for the numerical integration corresponds to a constant height (= constant Shields number) along the integration span to mimic again the experiments.

The present pipe experiment is particularly well suited for detailed experimental investigations of bed height evolution. Beside, the theoretical results in two-dimensional channels can be easily switched to the truncated pipe geometry (pipe flow having a particulate bed). However, it is not amenable to accurate measurement of particle motion inside the moving sediment because of image deformation caused by refraction effects of the glass tube. A recent study by [START_REF] Lobkovsky | Erosion of granular bed driven by laminar fluid flow[END_REF], that we have become aware of during the completion of this manuscript, precisely examine individual particle motion using matched-index technique in a rectangular-pipe at laminar low conditions. [START_REF] Lobkovsky | Erosion of granular bed driven by laminar fluid flow[END_REF] provide the vertical profile of the horizontal granular velocity and of the volume fraction (see their figure 2) and by integrating the product of these quantities over the height obtain the granular flux. Interestingly, the observed volume fraction is approximately constant φ 0 ≈ 0.55 except on the upper particle layer where it goes to zero. This confirms the existence of a transition layer, i.e. the region (II), assumed in the present model. [START_REF] Lobkovsky | Erosion of granular bed driven by laminar fluid flow[END_REF] estimate the Shields number from the flow rate but choose to neglect the side wall effects. While they can be accounted by a constant numerical factor β = 1.85 found by [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF] in the circu-lar pipe geometry, these effects evolve with the fluid gap by an amount of 8 -29% as estimated by [START_REF] Lobkovsky | Erosion of granular bed driven by laminar fluid flow[END_REF] in the rectangular geometry. This leads to a critical Shields number ≈ 0.3 which differs from that ≈ 0.12 found in the experiments of [START_REF] Charru | Erosion and depositionof particles on a bed sheared by a viscous flow[END_REF], [START_REF] Loiseleux | Onset of erosion and avalanches for an inclined granular bed sheared by a continuous laminar flow[END_REF], and [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF]. [START_REF] Lobkovsky | Erosion of granular bed driven by laminar fluid flow[END_REF] have scaled their observed granular flux with an inertial scaling ∆ρgd 3 /ρ f and found the fitted law ∝ (θ -θ c ) n with n = 1.89 ± 0.25 (see their figure 4). We propose an alternative plot against θ in figure 8. The viscous scaling given by equation (3.14) provides a good collapse of the experimental data. Although ranging over half a decade, the agreement with a cubic variation in θ (solid line) seems good as a power law fit gives q p /(∆ρgd 3 /η e ) ∝ θ n with n = 3.9 ± 0.3. However, the coefficient in equation (3.14) φ 0 /24θ c2 ≈ 0.025 using θ c = 0.12 and φ 0 = 0.55 gives an overestimation by a factor 40 (dashed line). This discrepancy is probably due to the neglect of the side wall effects mentioned above (an error of ≈ 30% in the definition of θ yields a factor ≈ 10 in the coefficient). Finally, while these experiments provide flow observation inside the mobile bed, the spatial resolution (4 data points in their figure 2) is too low to yield precise comparison with the parabolic flow profile predicted by equation (3.15) or with the linear variation of the bed-load thickness with the Shields number given by equation (3.16). Clearly, further detailed experiments are needed.

Discussion and conclusion

We now turn to the discussion of the closures used in the two-phase model which have led to the momentum equation for the mixture and the Brinkman equation for the fluid velocity used in the present study.

First, we have assumed that the effective stress tensor associated with the fluid phase was of Newtonian form and have chosen for simplicity to use the Einstein effective viscosity η e = η(1 + 5φ 0 /2) in the present concentrated situation. Figures 6 and7 show that using a time-scale = (6Re 2D /Ga) 1/2 (D/d)(η e /∆ρgd) built on the effective Einstein viscosity provides the good order of magnitude. This may seem surprising as the Einstein viscosity is supposed to be valid for dilute suspensions while, in the present experimental configuration, the particles are in contact and thus the volume fraction is not far from close packing. There are many formulae describing the effective viscosity at high concentrations and most of them diverge at close packing, see e.g. [START_REF] Stickel | Fluid mechanics and rheology of dense suspensions[END_REF]. Clearly, the effective viscosity cannot be described by such formulae as this would increase the time-scale by many orders of magnitude and would not match the experimental data. In fact, as the particles are in contact, higher-order hydrodynamic interactions such as lubrication may be disregarded and the viscosity may just reduce to the Einstein viscosity as suggested to us by [START_REF] Morris | [END_REF]. The issue is yet to be settled. Secondly, we have chosen to describe the stress tensor of the particle phase by a simple Coulomb model. This model neglects the shearing effect of the granular medium. A refined closure would be to consider a constitutive law which has been first derived for dry granular flows and has been then found successful for submarine granular flows, see GDR Midi (2004), Jop, Forterre & Pouliquen (2006[START_REF] Cassar | Submarine granular flows down inclined planes[END_REF], and [START_REF] Pouliquen | Flow of dense granular material: towards simple constitutive laws[END_REF]. The stress tensor of the granular material when sheared at a shear rate ∂u p i /∂x j + ∂u p j /∂x i under a confinement pressure p p is then written as ratio between the time of rearrangement of a particle when it is displaced from its hole and the time taken by the particle to move from one hole to the next. When an interstitial viscous fluid is present, this dimensionless number has been found to be

σ p ij = -p p δ ij + τ p ij = -p p δ ij + µ(I)p p ( ∂u p i ∂x j + ∂u p j ∂x i )/| ∂u p i ∂x j + ∂u p j ∂x i |. ( 6
I = η p p α | ∂u p i ∂x j + ∂u p j ∂x i |, (6.2)
with α = 0.01. The following expression for the friction coefficient has been proposed

µ(I) = µ s + µ 2 -µ s 1 + I0 I , (6.3)
where I 0 , µ 2 , and µ s are constant which depends upon the particle material and shape used. When the granular medium is not sheared, this rheology is not valid and one simply writes that the shear rate is zero as in the simple Coulomb model. Note that µ s is the tangent of the angle of repose and is identical to the friction coefficient µ used in the Coulomb model. We have solved numerically and analytically the Brinkman equation for the fluid phase and the momentum balance equation for the mixture with this new rheology, see figure 9. Near the thresholds of motion, I ≪ I 0 and µ(I) ≈ µ s + (µ 2 -µ s )I/I 0 . This simplified rheology yields the same critical Shields number, θ c = µ s φ 0 /2. For I ≫ I 0 , µ(I) ≈ µ 2 and the same solution as that found in the Coulomb friction case is obtained by substituting µ by µ 2 . The simple expression (3.14) for the bed-load flow-rate becomes

q p / ∆ρgd 3 η e = φ 0 µ s µ 2 2 θ c 24 θ θ c 3 . (6.4)
In the case of glass spherical particles, µ 2 = 0.82, see [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF], and q p is only modified by a factor 1/4. Figure 9 shows that using this sophisticated rheology does not change significantly the velocity profiles. The only difference lies inside the granular medium in motion where the fluid and particles profiles are overestimated when using the Coulomb friction.

In summary, we have proposed a two-phase model to describe bed-load transport in the laminar viscous regime, i.e. the flux of particles in a flat mobile bed submitted to laminar flows. The fluid phase has been assumed to be a Newtonian viscous liquid with Einstein dilute viscosity formula applied to the concentrated situation while the particle phase to have Coulomb solid friction with the shear stress proportional to the pressure. Adopting the more sophisticated granular rheology of [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF] does not change significantly the results as it results in switching the friction coefficient, i.e. the angle of repose, to a larger coefficient.

Region Ū ūp pp (I) µφ0( hp -hc -d/2)(ȳ -1) + Ū0 undefined 0 (II) η ηe µφ0 » - (ȳ -hp + d) 3 6d + (ȳ -hc) 2 2 - Ū φ0 2 d (ȳ -hp) 2 (III) η ηe µφ0 `ȳ -hc ´2 2 Ū φ0 " hp -ȳ - d 2 « (IV) 0 0 φ0 " hp -ȳ - d 2 «
We have applied this continuum model to viscous shearing flows. The relevant equations have been found to be the Brinkman equation for the fluid phase and the momentum balance equation for the mixture. These equations have been solved numerically assuming a constant volume fraction of the solid phase except on a transition region in which the volume fraction goes to zero on one particle radius as observed experimentally by [START_REF] Goharzadeh | Transition layer thickness at a fluid-porous interface[END_REF] and more recently by [START_REF] Lobkovsky | Erosion of granular bed driven by laminar fluid flow[END_REF]. Albeit being questionable and thus only qualitative close to incipient motion, this continuum approach is able to provide an onset of motion for the particle phase and a description of the flow of the mobile granular layer. At some distance from threshold where the continuum approach is more realistic as the mobile layer is larger than one particle diameter, there is very little slip between the two phases and the velocities inside the mobile bed have approximately a parabolic profile. Both the thickness of the mobile layer and the particle flux can be calculated analytically. When the Poiseuille (or Couette) flow is not significantly perturbed, the bed-load thickness varies linearly with the Shields number whereas the particle flux cubically with it.

We have compared these predictions with experimental observations in pipe flow. The cubic law for the particle flux seems satisfactory for describing the velocity of the kinematic waves which are triggered at the entrance of the tube as well as the time-evolution of the bed height in conditions of bed-load transport for pipe flows. This cubic law seems also to compare favourably with the recent experiments of [START_REF] Lobkovsky | Erosion of granular bed driven by laminar fluid flow[END_REF] in a rectangular pipe. to that obtained in the case of a Poiseuille flow. At some distance from this threshold, the mobile layer is thicker than a grain size and the slip between the fluid and the solid phases disappears. This last feature is recovered in the analytical calculation using the same hypotheses as in § 3.3, see tables 7 and 8. Using the same definition as in § 3.3, the critical Shields number for having a particle monolayer in motion is given again by θ c = µφ 0 /2. In the limit 1 -hp ≫ d and showing that the perturbation to the Couette flow due to the motion of the granular medium is small. The bed-load flow-rate given in table 8, can be further simplified to yield q p / ∆ρgd 3 η e = φ 0 θ c 12 1 2

θ θ c 3 + θ θ c - 1 5 . (A 2)
This leads to the simpler expression (3.14) having the particle flux varying cubically with the Shields number at some distance from threshold. This law of practical interest can be used for 0.3 θ 10 as can be seen in figure 11. Again, if one skips region (II) in which the volume fraction varies from φ 0 to 0 on one particle radius and if one assumes θ c = µφ 0 /2 [equation (3.11)], one recovers this simpler expression (3.14) for the particle flux as well as the parabolic profile of the velocities inside the moving bed [equation 3.15] and the linear variation of the bed-load thickness with θ [equation 3.16].

  θθ c ) 1.60

Figure 1 .

 1 Figure 1. Sketch of a particle bed submitted to a Poiseuille (left) or a Couette (right) flow in a two dimensional channel.

Figure 2 .

 2 Figure 2. Numerical velocity profile for the fluid (+) and solid (×) phases at different steps of the computation in the case of particles of batch A in fluid 3 at θ = 0.57.

Figure 3 .

 3 Figure 3. Numerical velocity profiles for the fluid (+) and the particles (×) in the case of particles of batch A in fluid 3 at φ 0 = 0.55 and θ = 0.048 (a), 0.094 (b), 0.62 (c) and analytical velocity profiles given in table 3 (solid line) and obtained by skipping region (II) (dashed line) (left). Blow-up of the profiles for the same conditions (right).

Figure 4 .

 4 Figure 4. Dimensionless numerical solid flux qp (•) versus numerical θ and the various analytical approximations of qp , table 4 (dotted line), equation (3.13) (dashed line), and equation (3.14) showing a cubic variation in θ (solid line), versus analytical θ in the case of particles of batch A in fluid 3.

  Figure5. Temporal evolution of the bed height, h p , for (i) batch A in fluid 2 with Q pipe = 6.45 10 -6 m 3 .s -1 (•) and in fluid 4 with Q pipe = 5.21 10 -6 m 3 .s -1 ( ), (ii) for batch B in fluid 1 with Q pipe = 3.37 10 -6 m 3 .s -1 ( ), in fluid 5 with Q pipe = 6.98 10 -6 m 3 .s -1 (△), in fluid 5 with Q pipe = 1.04 10 -5 m 3 .s -1 (⋄), and in fluid 3 with Q pipe = 2.56 10 -6 m 3 .s -1 (▽), and (iii) for batch C in fluid 3 with Q pipe = 2.13 10 -6 m 3 .s -1 (+) and with Q pipe = 5.16 10 -6 m 3 .s -1 (×).

Figure 6 .Figure 7 .

 67 Figure 6. Dimensionless kinematic velocity, c * = αL * /t * p , as a function of the initial Shields number, θ start pipe [equation (5.4)], for batch A (•), batch B ( ), and Batch C (△) and comparison with the present model [equation (5.3), solid line], the model of Charru & Mouilleron-Arnould (2002) [equation (5.5), dotted line], and that of Charru & Hinch (2006) [equation (5.6), dashed line]. The error bars are only indicated for 3 data points for indication.

  Figure8. q p /(∆ρgd 3 /η e ) versus θ for the experimental data of[START_REF] Lobkovsky | Erosion of granular bed driven by laminar fluid flow[END_REF] for Q pipe = 500 (•), 750 (△), 1000 ( ), 1500 (⊳), 2000 (⋄), and 2500 (⊲) cm 3 min -1 and equation (3.16) (dashed line) taking the experimental values ∆ρ = 1.474 g/cm 3 , d = 700 µm, and the Einstein viscosity η e = η(1 + 5φ 0 /2) with η = 24.9 cP as well as the adjustable parameter φ 0 /24θ c2 ≈ 0.025 using θ c = 0.12 and φ 0 = 0.55.

Figure 9 .

 9 Figure9. (a) Numerical velocity profiles for the fluid (∆) and the solid ( ) with a Coulomb rheology and profiles for the fluid (+) and the particles (×) with the more refined µ(I) rheology in the case of particles of batch A in fluid 3 at φ 0 = 0.55 and q f = 6.2 10 -3 m 2 s -1 . (b) Blow-up for the same conditions.

Figure 10 .

 10 Figure 10. Numerical velocity profile for the fluid (+) and the particles (×) and analytical velocity profiles for the fluid and the particles (-) in the case of particles of batch A in fluid 3 at (a) θ = 0.89 and (b) θ = 4.30 (left). Blow-up of the profiles for the same conditions (right).

Figure 11 .

 11 Figure 11. Same as figure 4 but for Couette flow and with analytical approximations, table 7 (dotted line), equation (A 2) (dashed line), and equation (3.14) showing a cubic variation in θ (solid line)

Table 5 .

 5 Particle characteristics

	Batch Composition d (µm)	ρp (g/cm 3 )
	A	Glass	132 ± 22 2.490 ± 0.003
	B	Polystyrene 538 ± 24 1.051 ± 0.002
	C	PMMA	193 ± 30 1.177 ± 0.002
	Fluid T ( • C) η (cP) ρ f (g/cm 3 )
	1	35	22 ± 2 1.031 ± 0.001
	2	35	30 ± 2 1.036 ± 0.001
	3	35	37 ± 2 1.038 ± 0.001
	4	35	33 ± 2 1.038 ± 0.001
	5	35	17 ± 1 1.029 ± 0.001

Table 6 .

 6 Fluid

characteristics density ρ p (also listed in table 5) was determined using a pycnometer and a fluid of measured density. Experiments were carried out using five different mixtures of distilled water and UCON oil 75H-90000 supplied by Chempoint. The density ρ f , listed in table 6, of these different mixtures was measured using a pycnometer. The viscosity η (also listed in table 6) was measured using a Falling Ball Viscometer from Gilmont Instruments (Barrington, USA). The mixture was maintained at a constant temperature T = 35 • C during all the experiments.

Table 7 .

 7 Same as table 3 but for Couette flow.

Table 8 .

 8 Same as table 4 but for Couette flow.

Acknowledgement

We would like to thank D. Lhuillier for discussions regarding the two-phase equations and closures, Y. Forterre, P. Jop, and O. Pouliquen for discussions regarding the granular rheology, E. J. Hinch for comments and in particular for suggesting the simplification of skipping region (II) in our analytical calculation, F. Charru for critical reading, and P. Cervetti, S. Martinez, and F. Ratouchniak for technical assistance. Funding from Institut Français du Pétrole, Agence Nationale de la Recherche (Project Dunes ANR-07-3 18-3892), and the Indo-French Centre for the Promotion of Advanced Research (IFCPAR Project No. 3404-1) is gratefully acknowledged. This work is part of the thesis of M. Ouriemi on ' Érosion, transport et instabilités d'un lit de particules dans un tube' (Université de Provence, 2007).

Appendix A. Particle bed submitted to a Couette flow

Similar conclusions as those found in § 3 can be drawn in the case of a Couette flow. In the fluid zone, the fluid velocity profile is approaching that of a Couette flow with a slip velocity at the top of the bed, see figure 10. The numerical calculations show the existence of a critical Shields number for onset of granular motion, θ c ≈ 0.06, similar