Mobility of bidisperse mixtures during bedload transport
Rémi Chassagne, Raphaël Maurin, Julien Chauchat, Philippe Frey
AAPS/123-QED
Mobility of bidisperse mixtures during bedload transport
R´emi Chassagne ∗ and Philippe Frey Univ. Grenoble Alpes, INRAE, IRSTEA,UR ETNA, 38000 Grenoble, France
Rapha¨el Maurin
IMFT, Univ. Toulouse, CNRS - Toulouse, France
Julien Chauchat
Univ. Grenoble Alpes, LEGI, CNRS UMR 5519 - Grenoble, France (Dated: September 3, 2020) a r X i v : . [ phy s i c s . g e o - ph ] S e p bstract The flow of segregated bidisperse assemblies of particles is of major importance for geophysicalflows and bedload transport in particular. In the present paper, the mobility of bidisperse segre-gated particle beds was studied with a coupled fluid discrete element method. Large particles wereinitially placed above small ones and it was observed that, for the same flow conditions, the bedloadtransport rate is higher in the bidisperse configuration than in the monodisperse one. Dependingon the Shields number and on the depth of the interface between small and large particles, differenttransport phenomenologies are observed, ranging from no influence of the small particles to smallparticles reaching the bed surface due to diffusive remixing. In cases where the small particleshardly mix with the overlying large particles and for the range of studied size ratios ( r < µ ( I ) rheology and it is demonstrated that the buried smallparticles are more mobile than larger particles and play the role of a “conveyor belt” for the largeparticles at the surface. Based on rheological arguments, a simple predictive model is proposed forthe additional transport in the bidisperse case. It reproduces quantitatively the DEM results fora large range of Shields numbers and for size ratios smaller than 4. The results of the model areused to identify four different transport regimes of bidisperse mixtures, depending on the mecha-nism responsible for the mobility of the small particles. A phenomenological map is proposed forbidisperse bedload transport and, more generally, for any granular flow on an erodible bed. ∗ [email protected] . INTRODUCTION In mountain rivers, the sediment bed is generally composed of a large range of grainsizes. This polydispersity leads to size segregation, which is largely responsible for ourlimited ability to predict sediment flux [1–3]. When segregating, small particles infiltratethe bed by kinetic sieving, falling down in holes formed by the matrix [4], and large particlesrise to the bed surface [5], resulting in inversely graded beds [6] which can be observed bothin flume experiments and in the field. In 1914, Gilbert [7] was one of the first to observeexperimentally that the introduction of finer sediments leads to an increase of sedimentmobility. This has then been extensively studied due to strong implications for sedimenttransport and fluvial morphology [3, 8, 9]. The mobility of granular assemblies is also akey question in the study of several geophysical flows such as debris flows, pyroclastic flows,snow avalanches and dune behavior. This, together with industrial applications, has led thegranular community to study the influence of the slope [10–12], basal friction [13, 14], totalvolume [15] and polydispersity [16–19] on particle mobility.Size segregation is often identified as the main mechanism responsible for the increasedmobility of a polydisperse bed. In laboratory experiments with natural materials, Bacchi etal. [20] showed that, due to kinetic sieving, small particles smooth the bed roughness andmake the above large particles more mobile. In bedload transport laboratory experimentswith a bidisperse bed, Dudill et al. [3] observed that the finer particles, after having infil-trated the first layers, drastically increased the sediment mobility. With two dimensionaldiscrete element method simulations (DEM), Linares-Guerrero et al. [17] measured the run-out distance of dry bidisperse granular avalanches. They observed an increased mobility ofthe avalanche due to the presence of small particles segregating at the base of the granularflow and acting as a lubrication layer. Similarly, Lai et al. [19] with DEM and laboratoryexperiments of granular collapse with fractal size distributions, observed the formation of abasal small particle layer increasing the total mobility. It seems therefore that size segre-gation, and in particular the formation of a small particle layer below large ones, plays animportant role in the increased mobility process. Despite the few studies presented above,there is still no clear understanding of the physical mechanisms responsible for the increasedmobility.Classically in bedload transport, bed mobility is interpreted in term of transport rate.3he dimensionless transport rate, or Einstein parameter, defined as Q ∗ s = Q s (( ρ p /ρ f − gd ) / , (1)is related to the dimensionless fluid bed shear stress, or Shields number, defined as θ = τ fb ( ρ p − ρ f ) gd , (2)where Q s is the transport rate per unit width, ρ p (resp. ρ f ) is the particle density (resp.fluid density), g is the gravity constant, d is the bed surface particle diameter and τ fb is thefluid bed shear stress. Considering their physical meaning and the link with the transportedgranular layer, the representative diameter for both the Shields and the Einstein numbersshould be taken as the surface layer particle diameter. It is classically chosen as the mediansurface diameter d or d (84% of the sediment is smaller than d ) [21]. However, literaturereview [3, 17, 19, 20] underlines the importance of the depth structure in the mobility ofthe granular bed and, in particular, the influence of buried small particles. Therefore,understanding the impact of the bed depth structure on transport laws is of particularimportance for an accurate description and prediction of turbulent bedload transport.While bedload transport has been mainly studied from the perspective of hydrodynamics,the present analysis illustrates the necessity to consider bedload as a granular phenomenon [2]and to describe the depth behaviour of the granular bed. In this paper, the mobility problemis therefore investigated from a granular perspective in the framework of the µ ( I ) rheology[22–24]. For dense granular flows, the dry inertial number I is the only dimensionlessparameter controlling the system, where I = d ˙ γ (cid:112) P p /ρ p , (3)with ˙ γ the shear rate and P p the granular pressure. The shear to normal granular stressratio µ p therefore depends only on the inertial number as µ p ( I ) = τ p P p = µ + µ − µ I /I + 1 , (4)where τ p is the granular shear stress and µ , µ and I are empirical coefficients fitted ondry experimental data. This rheology has been derived in monodisperse configurations andextended to bidisperse configurations in two dimensions [25] and three dimensions [26]. In arecent work, Maurin et al. [27] studied the rheology of dense granular flows during bedload4ransport using a coupled fluid-DEM model. Despite the presence of water, they showedthat the dry inertial number is still the controlling parameter. They found the µ ( I ) rheologyto be valid in bedload transport over a wider range of inertial numbers and proposed anotherset of parameters than the one proposed by GDR Midi [22] with µ = 0 . µ = 0 .
97 and I = 0 . µ ( I ) rheology framework in section IV and an explanationfor the increased mobility is presented. Based on rheological arguments, a simple predic-tive model for the additional transport is derived and compared with DEM simulations insection V. Finally the results are discussed in section VI. II. NUMERICAL MODEL AND SETUP
Our numerical model is a three dimensional discrete element method (DEM) using theopen source code YADE [28] coupled with a one-dimensional (1-D) turbulent fluid model.It has been derived and validated with particle-scale experiments [29] in [30] and extendedto bi-disperse configurations in [31]. It is briefly presented here but the interested readershould refer to Maurin et al. [30] for more details on the model and its validation. The DEMis a Lagrangian method based on the resolution of contacts. The inter-particle forces aremodelled by a spring-dashpot system [32] of stiffness k n in parallel with a viscous dampercoefficient c n (corresponding to a restitution coefficient of e n = 0 .
5) in the normal direction;and a spring of stiffness k s associated with a slider of friction coefficient µ g = 0 . k n and k s are computed in order to stay within the rigidgrain limit [30, 33]. The particles are additionally submitted to gravity, fluid buoyancy andturbulent drag force [30]. Considering a particle p , the buoyancy force is defined as f pb = − πd p ∇ P f x p , (5)5nd the drag force as f pD = 12 ρ f πd p C D || u f x p − v p || (cid:16) u f x p − v p (cid:17) , (6)where d p denotes the diameter of particle p , u f x p is the mean fluid velocity at the positionof particle p , P f x p is the hydrostatic fluid pressure at the position of particle p and v p isthe velocity of particle p . The drag coefficient takes into account hindrance effects [34] as C D = (0 . . /Re p )(1 − φ ) − . , with φ the packing fraction and Re p = || u f x p − v p || d p /ν f the particle Reynolds number, ν f being the kinematic viscosity.At transport steady state, the total granular phase (of small and large particles) only hasa streamwise component with no main transverse or vertical motion. In such a case, the 3-Dvolume averaged equation for the fluid velocity reduces to a 1-D vertical equation in whichthe fluid velocity is only a function of the wall-normal component, z , and is aligned withthe streamwise direction (see [35]) as ρ f (1 − φ ) ∂u fx ∂t = ∂S xz ∂z + ∂R xz ∂z + ρ f (1 − φ ) g x − n (cid:10) f pf x (cid:11) s , (7)where ρ f is the density of the fluid, S xz is the effective fluid viscous shear stress of a Newto-nian fluid of viscosity ν f . R xz is the turbulent fluid shear stress based on an eddy viscosityconcept R xz = ρ f (1 − φ ) ν t ∂u fx ∂z . (8)The turbulent viscosity ν t follows a mixing length approach that depends on the integral ofthe solid concentration profile to account for the presence of particles [36] ν t = l m | ∂u fx ∂z | , l m ( z ) = κ (cid:90) z φ max − φ ( ζ ) φ max dζ, (9)with κ = 0 .
41 the Von-Karman constant and φ max = 0 .
61 the maximum packing of thegranular medium (random close packing). The term n (cid:10) f pf x (cid:11) s represents the momentumtransfer associated with the interaction forces between fluid and particles. It is computed asthe horizontal solid-phase average of the momentum transmitted by the drag force to eachparticle.The fluid model is classical in sediment transport [30, 35, 37–40] and is only closed usinga mixing length model and a closure for the drag force formulation. The latter are usualin the literature, and it has been shown in [30, 41] that the results obtained in terms of6igure 1: A typical numerical setup. Initially N l layers of large particles ( d l = 6 mm) aredeposited by gravity on N s layers of small particles ( d s = 3 mm). The fluid of depth h w flows by gravity due to the slope angle α and entrains particles.granular behavior are very weakly sensitive to the fluid closure adopted.The numerical setup is presented on figure 1. In the following, subscripts l and s denotequantities for large and small particles respectively. Initially, small particles of diameter d s = 3 mm and large particles of diameter d l = 6 mm are deposited by gravity over a roughfixed bed made of small particles. The size of the 3-D domain is 30 d s × d s in the horizontalplane in order to have converged average values [30] and is periodic in the streamwise andspanwise direction. The number of particles of each class is assimilated to a number oflayers, N s and N l . They represent in terms of particle diameter the height that would beoccupied by the particles if the packing fraction was exactly φ max = 0 .
61, the maximalpacking fraction. Equivalently, at rest, the volume occupied by large particles (resp. smallparticles) is 0 . × d l × d l × N l d l (resp. 0 . × d l × d l × N s d s ). Therefore, specifying N l and N s gives the number of particles in each class. The height of the bed at rest isthus defined by H = N s d s + N l d l . The bed slope is fixed to 10% ( α = 5 . ◦ ), representativeof mountain streams. Since this study mainly focuses on cases where the bed surface iscomposed of only large particles, the Shields number definition is based on the large particlediameter as θ = τ f / (( ρ p − ρ f ) gd l ), where τ f = ρ f gh w sin ( α ) is the fluid bed shear stress, with7 a) (b) (c) Figure 2: Illustration of some considered configurations. (a) N l = 2, (b) N l = 4, (c)monidisperse case h w the water depth. Simulations were performed for Shields numbers ranging from 0 . N l = 1, 2, 3 and 4 which will be compared with a monodisperse largeparticle configuration considered as a reference case (see figure 2). In each case N s variedin order to keep the bed height H constant equal to H = 8 . d l for θ ≤ . H = 10 . d l for 0 . < θ ≤ . H = 16 . d l for larger Shields numbers. This increase in the bedthickness was necessary in order to ensure an erodible bed bottom boundary condition. Theorigin of the vertical axis is set at the top of the particle bed at rest. The interface position,describing the transition between large and small particles, is therefore defined geometricallyas z i = − N l d l .At the beginning of each simulation, the fluid flows by gravity and sets particles intomotion. After approximately 20 seconds, a dynamical equilibrium is achieved between thefluid flow and the transport of sediment. The results are then time-averaged over a 280stime period to ensure converged results. A mixed layer forms at the interface between smalland large particles resulting from an equilibrium between diffusion and size segregation.The present study focuses on the relation between the fluid forcing and sediment transportonce the steady state is achieved. Similarly to the Shields number, the Einstein parameteris defined with the large particle diameter as Q ∗ s = Q s / (cid:0) ( ρ p /ρ f − gd l (cid:1) . , where Q s = (cid:82) z φv px dz is the transport rate per unit width, and v px is the bulk streamwise particle velocity.The horizontal averaged concentration of small (resp. large) particles is defined as φ s (resp.8 l ). By definition, the two concentrations sum to φ the total granular concentration, φ s + φ l = φ. (10) III. ENHANCED MOBILITY DUE TO BIDISPERSITY (a) (b)
Figure 3: (a) Solid transport rate as a function of the Shields number for all simulations.(b) Increased transport rate in percentage compared with monodisperse configuration.In figure 3a is plotted the steady state dimensionless solid transport rate as a function ofthe Shields number. In all configurations, the dimensionless transport rate increases with theShields number. The transport rate is remarkably stronger in all bidisperse configurationswith respect to the monodisperse case, evidencing enhanced particle mobility. Figure 3bshows the bidisperse transport relative to monodisperse configurations, increasing up to 50%.The increase of transport is almost linear with the Shields number and is stronger when thenumber of layers of large particles N l is small. Indeed, for a lower N l , small particles arecloser to the surface (see figure 2) and are more likely to influence transport. This indicatesthat the depth of the interface between large and small particles, z i , plays a role in thetransport efficiency. At low Shields numbers and for N l = 4, almost no increase of transportis observed. In that case, the interface position is too deep to affect the bed mobility, andthe bidisperse bed behaves as if it were monodisperse. Overall, without modification of thefluid forcing, a substantial increase of transport is observed just by changing the particle9ize in the bed depth profile. (a) θ ∼ . N l = 4 (b) θ ∼ . N l = 2 (c) θ ∼ . N l = 1 Figure 4: Transport profiles of each class of particles for different configurations andShields numbers. The transport profile of the large particles in the monodisperse case forthe same Shields number is also plotted for comparison. Note changes in the abscissa scale.To expand the transport description, the local transport rate of each class of particleis defined as q is ( z ) = φ i ( z ) v px ( z ), where φ i ( z ) is the concentration of particle class i = l, s .Figure 4 shows the local transport rate depth profile of each class of particles for differenttypical configurations. The transport rate of large particles in the monodisperse case is alsoplotted in black dashed line for comparison. For θ ∼ . N l = 4 (figure 4a), almostno increase of transport ( ∼ θ ∼ .
55 and N l = 1(figure 4c) the transport of small particles is even stronger and small particles are present upto the bed surface, while they remained buried in the previous configuration (figure 4b). Itis therefore possible to draw two main conclusions. First, the observed increase of transportis a direct consequence of the mobility of the small particles. Second, even the large particletransport is significantly higher than in the monodisperse case.Two types of phenomenology are observed in the results. On the one hand small and10arge particles remain well separated, with small particles buried deep in the bed (figure 4a,b). On the other hand, small and large particles are mixed at the surface (figure 4c). Thewidth of the transition between small and large particles depends on the relative importanceof segregation over diffusion, the ratio of which can be defined as the Peclet number P e [31].If diffusion is strong enough compared to segregation, small buried particles can reach thesurface. To characterise the surface state, the surface diameter is computed as the meanparticle diameter above z = 0 as d surf = (cid:82) + ∞ φ s ( z ) d s + φ l ( z ) d l dz (cid:82) + ∞ φ s ( z ) + φ l ( z ) dz . (11)The non-dimensional surface diameter is set between 0 (only small particles at surface) and1 (only large particles) with the following transformation¯ d surf = d surf − d s d l − d s . (12)Figure 5 shows in scatter plot the value of the surface diameter as a function of the Shieldsnumber and the number of layers of large particles. The domain is clearly separated intotwo parts deliminated by the dashed line. Above the dashed line, the bed surface is onlycomposed of large particles while below it is composed of a mixture of both small andlarge particles. For a given value N l , there exists a transition Shields number θ t ( N l ) whichseparates a monodisperse from a bidisperse bed surface. For θ < θ t , diffusion is weakcompared to segregation, while for θ > θ t it is strong enough to move small particles up tothe bed surface. This therefore indicates that the Peclet number P e depends on the Shieldsnumber. In addition θ t increases with N l . Indeed, when N l increases, the transition depth z i between small and large particles is deeper in the bed and diffusion needs to be even strongerfor the small particles to reach the surface. For N l = 4 the surface is always composed oflarge particles. There is no doubt that increasing again the Shields number will eventuallybring small particles at the surface. Two simulations for N l = 0 . a priori higher11igure 5: Mean surface diameter as a function of the Shields number and the large particlenumber of layers. The dashed line shows the transition between a large particle surfacestate to a mixture surface state.for a mixture surface state. In cases where the small and large particles are well separated(above the dashed line), the increased transport rate cannot be attributed to a fluid effect.Indeed the length over which the fluid shear stress is fully transferred to the granular bedis much smaller than the grain size (see [42], [12]), and it is verified in appendix A thatit is indeed fully transferred to the granular bed below z = 0. The increased transport istherefore necessarily due to a granular process. In the next section, the study focuses onlyon the configurations where small and large particles are well separated and where the bedsurface is composed only of large particles. The granular process responsible for the increaseof mobility is investigated through a mechanical analysis of the granular bed properties. IV. INTERPRETATION AS A GRANULAR PROCESS
The granular stress tensor can be computed from the DEM. Considering a horizontalslice of volume V , the granular stress tensor is calculated as [43, 44] σ pij = − V (cid:88) p ∈ V m p v (cid:48) pi v (cid:48) pj − V (cid:88) c ∈ V f ci b cj , (13)where the sum is performed over the ensemble of particles p and contacts c inside the volume V , v (cid:48) pk = v pk − (cid:104) v pk (cid:105) s is the k component of the spatial velocity fluctuation of particle p , f c isthe interaction force at contact c on particle α by particle β and b c = x β − x α is the branch12ector. Due to the one dimensional structure of the flow, Maurin et al. [27, 41] showed that,in the steady state bedload configuration, σ pzz = T r ( σ p ) / σ pxz . The granular stress can therefore be described by only two scalarparameters which are the granular pressure P p = σ pzz and the shear stress τ p = σ pxz .Figure 6a compares, for θ ∼ .
45, the monodisperse and the bidisperse ( N l = 2) com-ponents of the stress tensor. The pressure and the shear stress exhibit the same behaviorin the monodisperse and bidisperse configurations. For the same forcing, the response ofthe bed in terms of granular stresses is therefore the same whatever the constitution ofthe bed. However, the transport profiles (figure 4b) show that the bidisperse bed is moremobile than the monodisperse one. This means that the dynamical response is dependenton the bed composition. This is analysed within the framework of the µ ( I ) rheology, relat-ing the friction coefficient µ p = τ p /P p to the inertial number I . The diameter to considerin the expression of the inertial number (3) is the local volume-averaged diameter [25, 26] d = φ s d s + φ l d l (which simplifies to d = d l in the monodisperse case). Following GDR Midi[22], the rheology of dense granular flows can be seen as follows. If µ p ≤ µ , where µ is thestatic friction coefficient, no motion is observed and I = 0. If µ p > µ , there exists a one toone correspondence between the friction cofficient µ p and the inertial number I .The friction coefficient is plotted in figure 6b and, as expected from the similarity of thegranular stress profiles (figure 6a), it is the same in the bidisperse and the monodisperseconfiguration. As a consequence, the inertial number profiles should be the same in bothconfigurations and that is indeed the case as observed in figure 6c. The dashed line ( ),defines a depth z such that µ p ( z ) = µ , the theoretical transition between static and densegranular flows. The dashed-dotted line ( ) shows the interface depth z i between smalland large particles.Figure 6d shows the bulk particle velocity for both configurations. For µ < µ or equiv-alently z < z , the inertial number and the velocity are indeed small but not exactly zero.This is due to non-local effects, that the µ ( I ) rheology is not able to capture [45, 46]. Itcorresponds to a quasi-static flow, or creeping regime, in which the velocity is exponentiallydecreasing into the bed ([47], [31]). In order to understand the increased mobility in thebidisperse configuration, the quasi-static regime is assumed to have a negligible impact ontransport and is not considered in this study. For z > z , as the friction coefficient is similar13 a) (b) (c) (d) Figure 6: Comparison of the monodisperse (dotted line) and the bidisperse N l = 2 (fullline) configuration for θ ∼ .
45. (a) Pressure and shear stress profiles, (b) frictioncoefficient profiles, (c) inertial number profiles and (d) velocity profiles. The dotted greenline corresponds to a translation of ∆ v = 0 . √ gd l of the monodisperse velocity profile(blue dotted line). The lower horizontal line at z ( ) separates the quasi-static regimefrom the flowing dense regime. The upper horizontal line at z i ( ) shows the transitionfrom small to large particles in the bidisperse configuration.in both configurations (see figure 6b), the inertial number is also supposed to be the same I b = I m , (14)where subscript b (resp. m ) denotes the bidisperse (resp. monodisperse) configuration. For z < z < z i , the particle diameter in the bidisperse simulation is d b ∼ d s , and d m = d l forthe monodispserse case. Equation 14 becomes d s ˙ γ b (cid:112) P p /ρ p ∼ d l ˙ γ m (cid:112) P p /ρ p . (15)The granular pressure being the same in both configurations (see figure 6a), gives˙ γ b ∼ d l d s ˙ γ m . (16)Integrating equation 16 from z to z ≤ z i , and assuming that the velocities are zero in z ,yields v pb ( z ) ∼ d l d s v pm ( z ) , (17)14nd therefore the velocity is higher in the bidisperse case than in the monodisperse case.This is perfectly observed in figure 6d. It means that for the same granular stress state,small particles are transported more easily than larger particles.For z > z i , the particle diameter is d l in both configurations and equation (15) simplifiesto ˙ γ b ∼ ˙ γ m , (18)and by integration from depth z i to z , v pb ( z ) ∼ v pm ( z ) + ( v pb ( z i ) − v pm ( z i )) ∼ v pm ( z ) + ∆ v, (19)meaning that the particle velocity profile in the bidisperse case is just a translation of thevelocity profile in the monodisperse case. In figure 6d is plotted, in the upper part of thebed, v pm ( z ) + ∆ v , with ∆ v = 0 . √ gd l measured in the DEM simulation. The obtainedcurve is completely superimposed on the velocity profile in the bidisperse configuration. Inboth configurations, the large particles at the top have exactly the same behaviour.The proposed granular analysis explains the observation made previously in figure 4, inwhich a layer of small particles was observed to be transported faster than larger parti-cles at the same depth. Small particles consequentely play the role of a conveyor belt forthe overlying particles and ∆ v represents a slip velocity. It additionally shows that theenhanced mobility is not a roughness effect, due to the reduction of roughness by smallerparticles below the large particle layer. Indeed, if particles do not move at the interface,∆ v = v pb ( z i ) − v pm ( z i ) is zero and no enhanced mobility is observed, as in figure 4a. The fluidorigin for the increased mobility can be discarded because the fluid shear stress is alreadyfully transferred to the granular shear stress below z = 0 (see appendix A). This analysisconfirms that the enhanced mobility originates in the granular rheological properties ofbidisperse beds.This rheological analysis gives a qualitative understanding of the granular bed behaviourin the bidisperse configuration. To be more quantitative, the previous conclusions are usedto predict analytically the additional transport in the bidisperse case.15 . A PREDICTIVE MODEL FOR THE ADDITIONAL TRANSPORT In this section, a simple model is derived, the purpose of which is to predict the additionaltransport observed in the bidisperse case. To obtain a predictive model, the additionaltransport will be expressed as a function of the monodisperse quantities ( φ m , v pm , etc...).The configuration is ideally simplified as a two layer problem in which small and largeparticles are completely separated at the interface depth z i . The mixed layer of small andlarge particles, observed in the bidisperse DEM simulations, is here neglected. Therefore itis assumed that the mixture concentration profiles are identical in the bidisperse and in themonodisperse configuration, ie. φ m ( z ) = φ b ( z ).The transport in the bidisperse case is expressed as, Q b = (cid:90) + ∞−∞ v pb ( z ) φ b ( z ) dz. (20)Below the interface between large and small particles, i.e. z ≤ z i , the previous analysis hasshown that v pb ( z ) = d l /d s v pm ( z ), while for z > z i , v pb ( z ) = v pm ( z ) + ∆ v . Splitting the integralinto two parts, below and above z i , placing the velocity expression into equation 20 andrecalling that φ b ( z ) = φ l ( z ) + φ s ( z ) is assumed to be equal to φ m ( z ), one obtains Q b = (cid:90) z i −∞ d l d s v pm ( z ) φ m ( z ) dz + (cid:90) + ∞ z i ( v pm ( z ) + ∆ v ) φ m ( z ) dz. (21)Distributing the second term and combining it with the first term, it comes Q b = Q m + ( d l d s − (cid:90) z i −∞ v pm ( z ) φ m ( z ) dz + (cid:90) + ∞ z i ∆ vφ m ( z ) dz. (22)where Q m = (cid:82) + ∞−∞ v pm ( z ) φ m ( z ) dz is the monodisperse transport rate. Recalling that ∆ v isindependent of z , the additional transport due to the presence of small particles can thereforebe expressed as∆ Q = ( d l d s − (cid:90) z i −∞ v pm ( z ) φ m ( z ) dz + ∆ v (cid:90) + ∞ z i φ m ( z ) dz = ∆ Q + ∆ Q . (23)The term ∆ Q represents the additional transport below the interface of the small particles,more mobile than larger particles. The term ∆ Q represents the additional transport of thelarge particles at the surface due to the conveyor belt effect. Note that in the monodisperselimit (i.e. d s = d l ), both terms vanish. This is obvious for ∆ Q . For ∆ Q , it is ∆ v = v pb − v pm ,which cancels in the monodisperse limit ( v pb = v pm ). Note that the additional transport in16 a) (b) Figure 7: (a) Dimensionless additional transport measured in the DEM simulations (fullsymbols) and predicted by equation (23) (empty symbols), for different values of theShields number and N l . Only cases for which the surface is exclusively composed of largeparticles are presented for readability. (b) Error between the total transport predicted byequation (23) and the transport computed with the DEM simulations.the bidisperse configuration (equation (23)) is expressed only as a function of monodispersevariables.In order to verify that the model is consistent with the transport mechanisms at play,equation (23) is first tested using DEM monodisperse simulations as inputs. The additionaltransport terms ∆ Q and ∆ Q are computed using the DEM velocity and concentrationprofiles v pm , φ m and estimating the slip velocity ∆ v directly on the DEM simulations. Thepredicted dimensionless additional transport rates are plotted in figure 7. The additionaltransport in the bidisperse case is very well predicted by equation (23) for all values ofShields number and for all numbers of layers of large particles. The small errors obtainedwith equation (23) show that the model contains the significant physical ingredients actingin this transport process.In practice, the concentration and velocity profiles, as well as the slip velocity, are difficultto obtain, and computing the additional transport due to the presence of small particles is17 a) (b) (c) (d) Figure 8: Comparison between idealized (dotted lines) and DEM profiles (full lines) in themonodisperse configuration for θ ∼ .
45. (a) Concentration profiles, (b) granular pressureand shear stress profiles, (c) friction coefficient profiles and (d) velocity profiles.not straightforward. In the following, a method to compute the two additional transportterms is proposed. The particles are assumed to be transported without dilatation of thebed. The concentration is therefore hypothesied constant and equal to φ max = 0 .
61 in thebed with the top of the bed exactly at z = 0 (see figure 8a).To compute the ∆ Q additional small particle transport term, the monodisperse veloc-ity profile for z ≤ z i needs to be estimated. It can be derived using the µ ( I ) rheology(equation 4). The stress state (normal and shear stresses) of the granular bed needs alsoto be computed. Based on the two-phase volume-averaged equations for turbulent bedloadtransport [40, 48] and for the idealized step concentration profile (figure 8a), the granularpressure and shear stress profiles can be expressed as (see appendix A) P p ( z ) = (cid:0) ρ p − ρ f (cid:1) g cos( α ) φ max z, (24) τ p = τ b + (cid:0) ρ f + ( ρ p − ρ f ) φ max (cid:1) z, (25)where τ b = ρ f gh w sin( α ) is the fluid bed shear stress. The friction coefficient can thenbe computed analytically as µ p = τ p /P p with these profiles. Inverting the µ ( I ) rheology(equation 4), replacing the inertial number I by its expression (equation 3) with the large18article diameter and integrating, a velocity profile is obtained v pm ( z ) = , µ p ( z ) < µ , (cid:90) zz (cid:32)(cid:115) P p ( ζ ) ρ p I d l µ p ( ζ ) − µ µ − µ p ( ζ ) (cid:33) dζ, µ ≤ µ p ( z ) < µ , (26)where µ = 0 . µ = 0 .
97 and I = 0 .
69 are the set of parameters proposed by Maurin etal. [27] for bedload transport. The integral can be computed numerically with the analyticalexpression of the granular pressure and of the friction coefficient and without any data fromthe DEM simulations.To verify that this derivation is consistent with the DEM simulations, figure 8 compares,for the monodisperse simulation at θ ∼ .
45, (a) the idealized concentration, (b) the pressureand shear stress, (c) the friction coefficient and (d) the velocity profile with the DEM results.The idealized step concentration profile obviously does not reproduce the dilatation of thebed at the surface. As a result, the pressure and shear stresses correspond with the DEMresults in most part of the bed but differ close to the surface. Similarly discrepancies nearthe bed surface appear for the friction coefficient and the velocity profiles. However, in theexpression of ∆ Q , the velocity and concentration profiles are needed only for z ≤ z i , wherethe idealized concentration and stresses agree very well with the DEM ones. Concerningthe velocity profile (figure 8d), the µ ( I ) rheology can not predict the quasi-static regime asalready mentioned (see inset). The velocity profile is well predicted in the dense regime butthe rheology fails to predict the velocity in the upper part of the bed for µ p ≥ µ , whichcorresponds to a more dilute flow regime. In order to use the predictive model, it is thereforenecessary that µ p ( z i ) < µ , which is the case in all our simulations and should be the casein classical bedload transport configurations. Otherwise, it would mean that small particlesare in the dilute flow regime and would be present at the bed surface, configuration whichhas already been discarded. With the velocity profile (26), it is now possible to computethe first additional transport term ∆ Q without any data from the DEM simulations.To compute the second additional transport term ∆ Q , both the ∆ v slip velocity and the (cid:82) + ∞ z i φ m ( z ) dz term need to be estimated. The second term represents the amount of largeparticles slipping above the small particles. With the idealized concentration profile, it can19e directly computed as (cid:90) + ∞ z i φ m ( z ) dz = φ max N l d l . (27)Lastly, the slip velocity remains to be estimated. By definition, for z ≥ z i , ∆ v = v pb ( z ) − v pm ( z ). It is therefore valid in z = z i , where v pb ( z i ) = d l /d s v pm ( z i ). The slip velocity istherefore finally given by ∆ v = (cid:18) d l d s − (cid:19) v pm ( z i ) , (28)with v pm ( z i ) which can be computed from the velocity profile equation (26) derived previously.All additional transport terms can now be computed and the total additional transport canbe expressed as ∆ Q = ( d l d s − φ max (cid:18)(cid:90) z i v pm ( z ) dz + N l d l v pm ( z i ) (cid:19) , (29)with v pm ( z ) given by equation (26). This additional transport term can be computed withoutany DEM data and uses only the µ ( I ) rheology.Equation (29) is tested and compared with the additional transport rate directly obtainedwith the DEM simulations in figure 9. The model predicts well the additional transportwith a maximum error around 20%, remaining smaller than 10% in most cases. The erroris generally smaller when N l is larger. For each configuration, there is a region where theerror is maximum. The Shields number at which the maximum error is reached seems todepend on the large particle number of layers. These results are discussed and interpretedin the next section. VI. DISCUSSION AND CONCLUSION
This study has shown that the additional transport evidenced in an inversely gradedbidispersed bed is a granular process. In a granular flow, small particles, being more mobilethan larger ones, play the role of a conveyor belt for the overlying large particles. Assumingthat large and small particles are completely separated and are transported without dilata-tion of the bed, a model for the enhanced transport has been derived based on rheologicalarguments. The results have shown that our model contains the significant physical ingredi-ents of the transport process and is able to predict acccurately the additional transport dueto bidispersity in bedload transport. The developed model allows improving upon classical20 a) (b)
Figure 9: (a) Dimensionless additional transport in the bidisperse case obtained with theDEM simulations (full symbols) and computed with equation (29) (empty symbols), fordifferent values of Shields number and N l . Only cases for which the surface is composed oflarge particles only are presented for readability. (b) Error between the total bidispersetransport predicted by equation (29) and the measured transport with the DEMsimulations.transport laws by taking into account not only the classical bed surface state, but the entiremobile granular bed structure.This model can also be used as a tool to interpret the different transport mechanismsobserved in this bidisperse granular flow configuration. The different regimes observed aresummarized in figure 10. The map has been built from the regions of validity of the model,the blue squares showing regions where the error between the model prediction and theDEM is less than 10% while the brown ones show regions where the error is higher. Thiscriterion enables us to define four different regimes of granular flows, corresponding to differ-ent granular depth structure and flowing mechanisms. Regime 1 corresponds to cases wheresmall and large particles are well mixed, with small particles present at the bed surface.In those cases, the additional transport is a combination of granular and fluid processes.Indeed, smaller particles at the surface are more easily entrained by the fluid flow and themixture of small and large particles can affect the flowing properties of the granular mobile21igure 10: Mapping of the four different observed phenomenologies in the bidispersetransport process. Each regime is illustrated with a typical simulation picture where thecreeping flow has been shaded in gray. Results are plotted in colored squares and split intotwo classes : blue (predicted transport error less than 10%) and brown (larger error).layer. Regime 2 corresponds to the domain of validity of the proposed model, where allassumptions are verified. In this regime, the fluid-driven large particles entrain the smallones, which create a so-called conveyor belt effect, due to their higher mobility. The tran-sition depth between small and large particles is here located in the dense granular flowregion. When the transition is located deeper in the bed, near or inside the creeping flowregion, the µ ( I ) rheology is no longer valid and the model predicts erroneously a zero veloc-ity inside the small particle layer (see inset figure 8d). This third regime therefore leads tosmall ( < a) (b) Figure 11: (a) Additional transport rate predicted by equation (29) (orange crosses) andcomputed from the DEM simulations (blue squares) for different size ratios at θ ∼ . N l = 2. (b) Non-dimensional surface diameter in the DEM simulations as defined inequation (12) as a function of the size ratio. Cases r = 2 . r = 4 are illustrated with apicture from the DEM simulation.stress, small particles are more mobile than larger particles and the effect observed for poly-disperse granular collapses [17, 19] or granular avalanches, for example, can be interpretedsimilarly. During the collapse, the small particles segregate and form a basal flowing layer,setting up a conveyor belt effect and increasing the runout distance of the collapse.The results obtained in this study can be put into perspective by considering the depen-dency of the results on the size ratio. Varying the size ratio between r = 1 . r = 4for a given configuration ( θ = 0 . N l = 2), one can evidence that the transport predictedby the model is valid up to r = 2 . r = 3 and r = 4 (see figure 11b). This indicates thatdiffusion remixing increases significantly, and can be related to the onset of inverse segrega-tion as observed in this range of size ratio by Thomas [50]. This link between diffusion and23nverse size segregation challenges our understanding of size segregation and deserves futurework. ACKNOWLEDGEMENTS
This research was funded by the French Agence nationale de la recherche, project ANR-16-CE01-0005 SegSed ’size segregation in sediment transport’. The authors acknowledgethe support of INRAE (formerly Irstea and Cemagref). INRAE, ETNA is member of LabexOsug@2020 (Investissements dAvenir Grant Agreement ANR-10-LABX-0056) and LabexTEC21 (Investissements dAvenir Grant Agreement ANR-11-LABX-0030).We are grateful to M. Church for reviewing and English corrections.
Appendix A: Derivation of the granular stress profiles
The two phase flow equations of bedload transport developed by [35] and [40] are con-sidered. For a unidirectional flow and for steady state condition, they read0 = ∂S xz ∂z + ∂R xz ∂z + ρ f (1 − φ ) g sin( α ) − n (cid:10) f pf x (cid:11) s , (A1)0 = ∂τ p ∂z + ρ p φg sin α + n (cid:10) f pf x (cid:11) s , (A2)0 = ∂P f ∂z + ρ f g cos α, (A3)0 = ∂P p ∂z + ( ρ p − ρ f ) φg cos α, (A4)where S xz and R xz are the viscous and turbulent fluid shear stresses, τ p is the granular shearstress, n (cid:10) f pf x (cid:11) s represents the transfer of momentum from the fluid to the solid phase and P f and P p are the fluid and granular pressure. [30] showed that the viscous fluid shearstress S xz is negligible in the bedload configuration and it will therefore not be taken intoaccount. Considering the following idealized concentration profile φ = φ max = 0 . , if z ≤ , , if z > , (A5)24nd by integration of equation (A4) between an elevation z and 0 where P p (0) is assumedto vanish, the two phase flow model predicts hydrostatic pressure for the granular phase P p ( z ) = − ( ρ p − ρ f ) φ max g cos( α ) z. (A6)Summing equation (A1) and (A2), a mixture momentum balance is obtained0 = ∂R xz ∂z + ∂τ p ∂z + (cid:0) ρ f + ( ρ p − ρ f ) φ (cid:1) g sin( α ) . (A7)In order to understand the partition between the fluid and granular stresses, equation (A7)is integrated between an elevation z and the free water surface h w where both shear stressesare assumed to vanish, leading to R xz ( z ) + τ p ( z ) = (cid:18) ρ f ( h w − z ) + ( ρ p − ρ f ) (cid:90) h w z φ ( ξ ) dξ (cid:19) g sin( α ) . (A8)In the pure fluid phase, where φ = 0 and therefore τ p ( z ) = 0, equation (A8) simplifies to R xz ( z ) = ρ f g sin( α )( h w − z ) , (A9)the classical expression of the turbulent fluid shear stress in a free surface flow. In thegranular bed the fluid shear stress rapidly decreases to zero and only the granular shear stressholds the mixture shear stress. With the idealized concentration profile (A5), equation (A8)simplifies to τ p ( z ) = (cid:2) ρ f ( h w − z ) − ( ρ p − ρ f ) φ max z (cid:3) g sin( α ) , (A10)which can be rewritten as τ p ( z ) = ρ f g sin( α ) h w − (cid:2) ρ p φ max + (1 − φ max ) ρ f (cid:3) gsin ( α ) z. (A11)The expressions of the granular pressure, fluid shear stress and granular shear stress obtainedfor the idealized step concentration are compared with the DEM profiles in figure 12. Theyagree in most parts except in the transition from the compacted granular bed to the purefluid phase that is not modeled by the idealized concentration profile. This step concentra-tion profile corresponds to an idealized situation where the fluid shear stress is completelytransmitted to the granular bed at the discontinuity ( z = 0). Focusing on the granularshear stress, the DEM and analytical profiles correspond almost perfectly as soon as z ≤ a) (b) Figure 12: Monodisperse case at θ ∼ .
45. (a) Granular pressure from DEM simulation(full line) and computed with analytical expression (A6) (dashed line). (b) Fluid andgranular shear stress from DEM simulation (full line) and computed from analyticalexpression A9 (dotted line) and expression (A11) (dashed line) [1] J.C. Bathurst, “Effect of Coarse Surface Layer on Bed-Load Transport,” Journal of HydraulicEngineering , 1192–1205 (2007).[2] P. Frey and M. Church, “Bedload: A granular phenomenon,” Earth Surface Processes andLandforms , 58–69 (2011).[3] A. Dudill, H. Lafaye de Micheaux, P. Frey, and M. Church, “Introducing Finer Grains IntoBedload: The Transition to a New Equilibrium,” Journal of Geophysical Research: EarthSurface , 2602–2619 (2018).[4] G. V. Middleton, “Experimental studies related to problems of flysch sedimentation,” in FlyschSedimentology in North America , Special Paper No. 7 (Business and Economics Science Ltd.,1970) j. lajoie ed., pp. 253–72.[5] S. B. Savage and C. K. K. Lun, “Particle size segregation in inclined chute flow of dry cohe-sionless granular solids,” Journal of Fluid Mechanics , 311–335 (1988).[6] J. M. N. T. Gray, “Particle Segregation in Dense Granular Flows,” Annual Review of FluidMechanics , 407–433 (2018).
7] G. K. Gilbert, “The transportation of d´ebris by running water,” USGS Professional paper 86US Geological Survey: Washington DC (1914), 10.1130/0-8137-2338-8.253.[8] K. M. Hill, J. Gaffney, S. Baumgardner, P. Wilcock, and C. Paola, “Experimental study ofthe effect of grain sizes in a bimodal mixture on bed slope, bed texture, and the transition towashload,” Water Resources Research , 923–941 (2017).[9] A. Dudill, P. Frey, and M. Church, “Infiltration of fine sediment into a coarse mobile bed: Aphenomenological study,” Earth Surface Processes and Landforms , 1171–1185 (2017).[10] A. Mangeney, O. Roche, O. Hungr, N. Mangold, G. Faccanoni, and A. Lucas, “Erosionand mobility in granular collapse over sloping beds,” Journal of Geophysical Research: EarthSurface , F3 (2010).[11] M. Farin, A. Mangeney, and O. Roche, “Fundamental changes of granular flow dynamics,deposition, and erosion processes at high slope angles: Insights from laboratory experiments,”Journal of Geophysical Research: Earth Surface , 504–532 (2014).[12] R. Maurin, J. Chauchat, and P. Frey, “Revisiting slope influence in turbulent bedload trans-port: Consequences for vertical flow structure and transport rate scaling,” Journal of FluidMechanics , 135–156 (2018).[13] C. Chedeville and O. Roche, “Autofluidization of pyroclastic flows propagating on roughsubstrates as shown by laboratory experiments,” Journal of Geophysical Research: Solid Earth , 1764–1776 (2014).[14] A. N. Edwards, S. Viroulet, B. P. Kokelaar, and J. M. N. T. Gray, “Formation of levees,troughs and elevated channels by avalanches on erodible slopes,” Journal of Fluid Mechanics , 278–315 (2017).[15] L. Staron and E. Lajeunesse, “Understanding how volume affects the mobility of dry debrisflows,” Geophysical Research Letters , 12 (2009).[16] J. C. Phillips, A. J. Hogg, R. R. Kerswell, and N. H. Thomas, “Enhanced mobility of granularmixtures of fine and coarse particles,” Earth and Planetary Science Letters , 466–480(2006).[17] E. Linares-Guerrero, C. Goujon, and R. Zenit, “Increased mobility of bidisperse granularavalanches,” Journal of Fluid Mechanics , 475–504 (2007).[18] R. M. Iverson, M. Logan, R. G. LaHusen, and M. Berti, “The perfect debris flow? Aggregatedresults from 28 large-scale experiments,” Journal of Geophysical Research: Earth Surface , , 12,181–12,189(2017).[20] V. Bacchi, A. Recking, N. Eckert, P. Frey, G. Piton, and M. Naaim, “The effects of kineticsorting on sediment mobility on steep slopes,” Earth Surface Processes and Landforms ,1075–1086 (2014).[21] A. Recking, “An analysis of nonlinearity effects on bed load transport prediction,” Journal ofGeophysical Research: Earth Surface , 1264–1281 (2013).[22] GDR MiDi, “On dense granular flows,” The European Physical Journal E , 341–365 (2004).[23] P. Jop, Y. Forterre, and O. Pouliquen, “A constitutive law for dense granular flows,” Nature , 727–730 (2006).[24] Y. Forterre and O. Pouliquen, “Flows of Dense Granular Media,” Annual Review of FluidMechanics , 1–24 (2008).[25] P. G. Rognon, J.N. Roux, M. Naa¨ım, and F. Chevoir, “Dense flows of bidisperse assembliesof disks down an inclined plane,” Physics of Fluids , 058101 (2007).[26] A. Tripathi and D. V. Khakhar, “Rheology of binary granular mixtures in the dense flowregime,” Physics of Fluids , 113302 (2011).[27] R. Maurin, J. Chauchat, and P. Frey, “Dense granular flow rheology in turbulent bedloadtransport,” Journal of Fluid Mechanics , 490–512 (2016).[28] V. Smilauer et al., Yade Documentation 2nd Ed. The Yade Project. (Zenodo, 2015).[29] P. Frey, “Particle velocity and concentration profiles in bedload experiments on a steep slope,”Earth Surface Processes and Landforms , 646–655 (2014).[30] R. Maurin, J. Chauchat, B. Chareyre, and P. Frey, “A minimal coupled fluid-discrete elementmodel for bedload transport,” Physics of Fluids , 113302 (2015).[31] R. Chassagne, R. Maurin, J. Chauchat, J. M. N. T. Gray, and P. Frey, “Discrete and contin-uum modelling of grain size segregation during bedload transport,” Journal of Fluid Mechanics (2020), 10.1017/jfm.2020.274.[32] T. Schwager and T. Poschel, “Coefficient of restitution and linear–dashpot model revisited,”Granular Matter , 465–469 (2007).
33] J.N. Roux and G. Combe, “Quasistatic rheology and the origins of strain,” Comptes RendusPhysique , 131–140 (2002).[34] J. F. Richardson and W. N. Zaki, “The sedimentation of a suspension of uniform spheresunder conditions of viscous flow,” Chemical Engineering Science , 65–73 (1954).[35] T. Revil-Baudard and J. Chauchat, “A two-phase model for sheet flow regime based on densegranular flow rheology,” Journal of Geophysical Research: Oceans , 619–634 (2013).[36] L. Li and M. Sawamoto, “Multi-Phase Model on Sediment Transport in Sheet-Flow RegimeUnder Oscillatory Flow,” Coastal Engineering in Japan , 157–178 (1995).[37] T. G. Drake and J. Calantoni, “Discrete particle model for sheet flow sediment transport inthe nearshore,” Journal of Geophysical Research: Oceans , 19859–19868 (2001).[38] T.-J. Hsu and P. L.-F. Liu, “Toward modeling turbulent suspension of sand in the nearshore,”Journal of Geophysical Research: Oceans , C6 (2004).[39] O. Dur´an, B. Andreotti, and P. Claudin, “Numerical simulation of turbulent sediment trans-port, from bed load to saltation,” Physics of Fluids , 103306 (2012).[40] J. Chauchat, “A comprehensive two-phase flow model for unidirectional sheet-flows,” Journalof Hydraulic Research , 15–28 (2018).[41] R. Maurin, Investigation of Granular Behavior in Bedload Transport Using a Eulerian-Lagragian Model. , PhD thesis, Univ. Grenoble Alpes (2015).[42] M. Ouriemi, P. Aussillous, and E. Guazzelli, “Sediment dynamics. Part 1. Bed-load transportby laminar shearing flows,” Journal of Fluid Mechanics , 295–319 (2009).[43] I. Goldhirsch, “Stress, stress asymmetry and couple stress: From discrete particles to contin-uous fields,” Granular Matter , 239–252 (2010).[44] B. Andreotti, Y. Forterre, and O. Pouliquen, “Granular Media: Between Fluid and Solid,”/core/books/granular-media/E1D234B8D868A9856C9E95B4750470AB (2013).[45] K. Kamrin and G. Koval, “Nonlocal Constitutive Relation for Steady Granular Flow,” PhysicalReview Letters , 178301 (2012).[46] M. Bouzid, M. Trulsson, P. Claudin, E. Cl´ement, and B. Andreotti, “Nonlocal Rheology ofGranular Flows across Yield Conditions,” Physical Review Letters , 238301 (2013).[47] B. Ferdowsi, C. P. Ortiz, M. Houssais, and D. J. Jerolmack, “River-bed armouring as agranular segregation phenomenon,” Nature Communications , 1363 (2017).[48] R. Jackson, The Dynamics of Fluidized Particles (Cambridge University Press, 2000).
49] M. Houssais, C. P. Ortiz, D. J. Durian, and D. J. Jerolmack, “Onset of sediment transportis a continuous transition driven by fluid shear and granular creep,” Nature Communications , 6527 (2015).[50] N. Thomas, “Reverse and intermediate segregation of large beads in dry granular media,”Physical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics , 961–974 (2000)., 961–974 (2000).