Discrete and continuum modelling of grain size segregation during bedload transport

Grain-scale discrete element simulations of bidisperse mixtures during bedload transport are used to understand, and model, bedload transport and particle-size segregation in granular media. For an initial distribution of fine particles on top of a coarse granular bed, this paper investigates the gravity driven percolation/segregation of the fine particles down into the quasi-static part of the bed. The segregation is observed to be driven by the inertial number at the bottom of the fine particle layer, and is independent of the number of fine particles. A novel travelling wave solution for the evolving concentration distribution is constructed using the continuum particle-size segregation model of Thornton, Gray & Hogg (J. Fluid Mech., vol. 550, 2006, pp. 1–25) and Gray & Chugunov (J. Fluid Mech., vol. 569, 2006, pp. 365–398). The observed behaviour is shown to be related to a local equilibrium between the influence of the concentration and of the inertial number. The existence of the exact solution relies on the segregation flux and the diffusion coefficient having the same dependency on the inertial number. This functional dependence allows the continuum model to quantitatively reproduce the discrete simulations. These results significantly improve on our understanding of the size segregation dynamics and represent a step forward in the up-scaling process to polydisperse granular flows in the context of turbulent bedload transport.

mass units, respectively, the following relation is obtained: where r = d l /d s is the size ratio between large and small particles and G is an unknown function. It is assumed that relation (1.2) can be written in a power law form, where the function G 1 should vanish when r = 1 in order to cancel segregation in the monodisperse limit.
In the theory of Savage & Lun (1988), the shear rateγ p is predicted to be an important controlling parameter for the rate of particle-size segregation. This is based on the observation that as particles flow downslope, the shear allows each layer of particles in the flow to move faster than the one beneath and hence the rate of finding gaps to fall into is expected to be proportional to the shear rate. Savage & Lun (1988) showed that this assumption produced good agreement between the theory and bidisperse experimental flows down inclined planes. More recent discrete element method (DEM) simulations of dry bidisperse mixtures in heap flow (Fan et al. 2014) have also found there to be a linear relation between the percolation velocity and the shear rate. However, the annular shear cell experiments of Golick & Daniels (2009) suggest that the segregation rate is also pressure dependent since the segregation rate dramatically slowed when a large pressure was applied on the top plate of their shear cell. This pressure dependence has also been observed in the DEM simulations of Fry et al. (2018), who found a power law relating the segregation velocity and the inertial number I (GDR MiDi 2004). A scaling with the inertial number represents an extension of the Savage & Lun (1988) theory in order to take into account both the effects of the shear rate and pressure on segregation. It corresponds in (1.3) to η = −ζ /2, which can be rewritten as w s √ gd l ∝ I ζ G 1 (r)G 2 (φ s ), (1.4) with I the large particle inertial number defined as suggesting that the segregation velocity is a function of the inertial number, the size ratio and the local concentration. The size ratio between large and small particles r is also a key parameter for size segregation. In the annular shear cell, Golick & Daniels (2009) showed that the segregation velocity is an increasing function of the size ratio with a maximum around r = 2. This effect was also observed in both two-dimensional (2-D) and three-dimensional (3-D) DEM simulations of sheared granular flows (Thornton et al. 2012;Guillard, Forterre & Pouliquen 2016).
The effect of the small particle concentration φ s has also been widely studied, and different forms for G 2 have been proposed (Bridgwater, Foo & Stephens 1985;Savage & Lun 1988;Dolgunin & Ukolov 1995;Gray & Thornton 2005;May et al. 2010;Fan et al. 2014;Gajjar & Gray 2014; Van der Vaart et al. 2015). The small particles have been shown to segregate less easily if they are more concentrated and the downward velocity should vanish for a pure phase of small particles. For this purpose, the simplest form G 2 (φ s ) = 1 − φ s , has been proposed by Dolgunin & Ukolov (1995). In oscillating shear cell experiments, Van der Vaart et al. (2015) measured the time necessary to achieve complete segregation with an initial mixture ranging from one small particle percolating into a bed of large particles, to a single large particle rising up through a small particle bed. They observed that both extreme cases were not symmetric, resulting in an asymmetry of the vertical velocity with the concentration, and Gajjar & Gray (2014) proposed a quadratic dependence of G 2 on φ s . Fan et al. (2014) and Jones et al. (2018) also observed this nonlinearity in dry bidisperse avalanche simulations, but concluded that the form of Dolgunin & Ukolov (1995) is a still good first-order approximation.
The literature review on dry granular flows given above underlines the qualitative understanding of size segregation for simple flows. Bedload transport can be seen as a granular flow, where the coupling with the fluid induces strong gradients in the vertical direction and a complex forcing, which could challenge the classical picture of segregation. Few studies have been made on gravity driven segregation in bedload transport and more remains to be done for a clear understanding of the processes at play. Ferdowsi et al. (2017) experimentally studied size segregation in laminar bedload transport and performed dry granular flow simulations. They studied the formation of armour, i.e. the segregation of large particles to the top, starting from a mixture of small and large particles. They showed that the process seems to be a granular phenomenon and reproduced their experimental results in the framework of continuum segregation modelling, using a generalization of the Gray & Thornton (2005) segregation model. Hergault et al. (2010) and Frey et al. (2020) studied size segregation in turbulent bedload transport, considering a quasi-2-D channel enabling particle tracking. They found that the particles move down as a layer into the bed, and related the segregation velocity to the granular shear rate.
Modelling turbulent bedload transport with segregation phenomenon at the particle scale is computationally demanding and can only be achieved for very small domains. In this context, one of the main goals for segregation in turbulent bedload transport is to be able to do up-scaling to the framework of macroscopic continuum modelling. In this case, it includes the modelling of the fluid phase, the granular phase and the evolution of the different classes of particles with respect to one another. Such a segregation model has been developed by Gray & Thornton (2005), Thornton, Gray & Hogg (2006) and Gray & Chugunov (2006). This three phase model is based on the assumption that the granular overburden pressure is not shared equally between large and small particles (Gray & Thornton 2005). Substituting this into the momentum conservation equation of each constituent and placing the problem within the framework of the mixture theory, an evolution equation for the phase of small particles is found: where u is the bulk velocity field, F s the segregation flux and D the diffusive remixing coefficient of the phase of small particles into the phase of large ones. In this model the fluid has only a passive role and acts through buoyancy to reduce gravity. The Modelling of grain size segregation during bedload transport 895 A30-5 physics relies on the expression of the segregation flux and of the diffusion flux. The former is defined as F s = φ s w s and can therefore directly be linked to the discussion above. Introducing the dimensional analysis (1.4) and considering the simplest form of Dolgunin & Ukolov (1995) for G 2 the following segregation flux is obtained: (1.7) In granular flows, diffusion has been mainly studied in the case of self-diffusion. By analogy to the thermal diffusion of molecules, Campbell (1997) tracked the random motion of particles in numerical simulations of dilute monodisperse sheared flows, and showed that the diffusion coefficient should scale with the shear rate. In dense 2-D granular Couette flow experiments, Utter & Behringer (2004) also observed that the diffusion coefficient was dependent on the shear rate. The diffusion mechanism considered in this paper is the mixing of one class of particles into another. This kind of diffusion is expected to be related to self-diffusion, but to the best of our knowledge, no study has been done in order to determine the physical mechanism controlling this diffusion. While the literature review above underlines the progress in the understanding and modelling of gravity driven segregation, more physically based parameterizations of the continuum segregation model are still lacking, in particular in the complex turbulent bedload transport configuration. In addition, no general description of size segregation, that would be valid in all configurations, has been proposed yet. From a geomorphological point of view, the impact of size segregation on transport rate is not well understood and an effort is still needed in the development of a model able to represent size segregation at the river scale.
In the present contribution, size segregation in turbulent bedload transport is studied numerically considering fluid-DEM simulations. The aim is to understand size segregation in the quasi-static part of the flow and to improve the continuum modelling parameterizations. In particular, the influence of the three main parameters which are the size ratio r, the inertial number I, the small particle concentration φ s , will be investigated through numerical experiments in the bedload configuration.
After presenting the numerical model ( § 2), simulations of turbulent bedload transport will be presented and analysed based on the different expected dependencies from the dimensional analysis ( § 3). The results will then be studied, through a theoretical analysis, in the framework of macroscopic continuum modelling and will highlight the local segregation mechanisms. This analysis shows that diffusion and segregation should have the same dependence on the inertial number and enables quantitatively reproducing the DEM results ( § 4).

Fluid-DEM model for bidisperse system
The numerical model used to simulate size segregation in turbulent bedload transport is a three-dimensional DEM using the open-source code YADE (Smilauer et al. 2015) coupled with a one-dimensional turbulent fluid model. This code has been validated with particle-scale experiments (Frey 2014) in Maurin et al. (2015) for mono-disperse situations. It was used to study bedload rheology (Maurin et al. 2016) and the slope influence (Maurin, Chauchat & Frey 2018). A brief summary of the model formulation as well as a description of its adaptation to bi-disperse situations is now given.

Granular phase
The DEM is a Lagrangian method based on the resolution of contacts between particles. For each spherical particle p, the motion of the particle is obtained from Newton's second law, and a momentum equation and an angular momentum conservation equation are solved: where m p , x p , ω p and I p are respectively the mass, position, angular velocity and moment of inertia of particle p. The three major forces are: f p g the gravitational force, f p c the inter-particle contact forces and f p f the interaction forces with the fluid. The inter-particle contact forces are classically defined as a spring-dashpot system (Schwager & Poschel 2007) composed of a spring of stiffness k n in parallel with a viscous damper of coefficient c n in the normal direction; and a spring of stiffness k s associated with a slider of friction coefficient µ p in the tangential direction. The normal and tangential contact forces read accordingly: where δ n (respectively δ t ) is the overlap between particles in the normal (respectively tangential) direction. For each class of particles, the values of k n and k s are computed in order to stay in the rigid limit of grains (Roux & Combe 2002;Maurin et al. 2015). The normal stiffness, in parallel with the viscous damper, defines a restitution coefficient representative of the loss of energy during collisions, which is fixed to e n = 0.5 for each contact. Finally, considering glass beads, the friction coefficient is fixed to µ p = 0.4 for all particles independently of their diameter. The third type of forces concern the interactions with the fluid which are restricted in this model to the buoyancy force and the drag force . They are defined as where d p denotes the diameter of particle p, u f x p is the mean fluid velocity at the position of particle p, P f x p is the hydrostatic fluid pressure at the position of particle p and v p is the velocity of particle p. The drag coefficient takes into account hindrance effects (Richardson & Zaki 1954) as C D = (0.4 + 24.4/Re p )(1 − Φ) −3.1 , with Φ the total solid phase (small and large particle) volume fraction and Re p = u f x p − v p d p /ν f the particle Reynolds number.

Fluid model
At transport steady state, the total granular phase (of small and large particles) only has a streamwise component with no main transverse or vertical motion. In such a case, it can be shown (Revil-Baudard & Chauchat 2013) that the 3-D volume-averaged equation for the fluid velocity reduces to a one-dimensional (1-D) vertical equation, in which the fluid velocity is only a function of the wall-normal component, z, and is aligned with the streamwise direction. The fluid pressure P f , is in this case given by the hydrostatic fluid pressure, as shown in Revil-Baudard & Chauchat (2013). The proposed fluid model is therefore a one-dimensional turbulent model in the streamwise x-direction with variables only depending on z and is inspired from the Euler-Euler model proposed by Revil-Baudard & Chauchat (2013) and Chauchat (2018). The volume-averaged momentum balance for the fluid phase is as follows: where ρ f is the density of the fluid, S xz is the effective fluid viscous shear stress of a Newtonian fluid, R xz is the turbulent fluid shear stress and n f p f x s represents the momentum transfer associated with the interaction forces between fluid and particles.
The viscous shear stress S xz is taken as with ν f the pure fluid viscosity. The Reynolds shear stress R xz is based on an eddy viscosity concept, The turbulent viscosity ν t follows a mixing length approach that depends on the integral of the solid volume fraction profile to account for the presence of particles (Li & Sawamoto 1995), with κ = 0.41 the von Kármán constant and Φ max = 0.61 the maximal packing of the granular medium (random close packing).
The total momentum n f p f x s transmitted by the fluid to the particles is computed as the sum over each class of the horizontal solid-phase average (denoted . s ) of the momentum transmitted by the drag force on each particle and introducing the expression of the drag force (2.5), where p s (respectively p l ) denotes the ensemble of the small (respectively large) particles, Φ s (respectively Φ l ) the volume fraction of small (respectively large) particles and d s (respectively d l ) the diameter of small (respectively large) particles.

R. Chassagne and others
The solid-phase volume average · s is defined following Jackson (2000) for which a cuboid weighting function F with the same length and width as the 3-D domain is applied. In the vertical direction, in order to capture the strong vertical gradient of the mean flow, the vertical thickness l z of the box is chosen as l z = d s /30. This choice of weighting function has been validated by comparison with mono-disperse turbulent bedload transport experiments in Maurin et al. (2015). Therefore the mean values at position z are computed as (2.14) where γ is a scalar quantity.  (2015) and Maurin et al. (2015) that the results obtained in terms of granular behaviour are very weakly sensitive to the fluid closure adopted. In this paper, the model is used for a bi-disperse situation to study size segregation in bedload sediment transport.

DEM simulations of segregation dynamics
In this section, the infiltration of fine particles into the bed made of larger particles is studied during turbulent bedload transport using the numerical model described in the previous section.

Numerical set-up
The numerical set-up is presented in figure 1. In the following, subscripts l and s denote quantities for large and small particles, respectively. Initially, large particles of diameter d l = 6 mm and fine particles of diameter d s = 4 mm (size ratio r = 1.5) are deposited by gravity over a rough fixed bed made of large particles. The particle and fluid densities are fixed respectively to ρ p = 2500 kg m −3 and ρ f = 1000 kg m −3 . The size of the 3-D domain is 30d l × 30d l in the horizontal plane in order to have converged average values  and is periodic in the streamwise and spanwise directions. The number of particles of each class is assimilated into a number of layers, N l and N s . The number of layers represents in terms of particle diameters the height that would be occupied by the particles if the packing fraction was exactly 0.61. Equivalently, the volume occupied by large particles (respectively small particles) is 0.61 × 30d l × 30d l × N l d l (respectively 0.61 × 30d l × 30d l × N s d s ). Therefore, fixing N l or N s fixes the number of particles of each class. The height of the bed at rest is defined by H = N l d l + N s d s and H is fixed to 10d l in all the simulations. The number of layers of fine particles N s varies from 0.01 (only a few small particles) up to 2 layers (corresponding to figure 1), while N l changes accordingly in order to keep H = 10d l . The bed slope is fixed to 10 % (5.7 • ), representative of mountain streams. The Shields number is defined as the dimensionless fluid bed shear stress θ = τ f /[(ρ p − ρ f )gd l ] and simulations are performed at θ 0.1. For this purpose, the free surface position is fixed and corresponds to a water depth of h = 3.1d l .
At the beginning of each simulation, the fluid flows by gravity and sets particles into motion. A first transient phase takes place, during which fluid and particles are accelerating. During this period, segregation is very fast and at the end of the transient phase, the small particles have already infiltrated into the first layers of the bed. In this paper, the study focuses on the dynamics of segregation once the system is at transport equilibrium and small particles have reached the quasi-static layer.
The horizontal-averaged volume fraction per unit granular volume of small (respectively large) particles is defined as φ s (respectively φ l ). By definition, the two volume fractions sum to unity, where Φ s (respectively Φ l ) is the volume fraction of small (respectively large) particles defined per unit mixture volume as in the previous section. Lastly, z c is defined as the vertical position of the centre of mass of small particles and dz c /dt as its velocity. The results are presented without dimension and the following dimensionless variables are considered:z In the following, the tildes are dropped for sake of clarity.

Results
Simulations have been performed for different numbers of layers of small particles N s . Figure 2(a) shows the temporal evolution of fine particle concentration profiles for the case N s = 2, corresponding to two layers of small particles deposited on top of a 895 A30-10 R. Chassagne and others the thickness of the small particle layer gets slightly larger. Profiles of concentration for different times are presented in figure 2(b). The concentration profiles exhibit a Gaussian-like shape. After the transient phase, neither the maximal value (see figure 2c) nor the width of the profiles evolve in time, suggesting that the small particles infiltrate the bed as a layer having a constant thickness and that this layer is just convected downward by segregation inside the layer of large particles. The vertical time-and volume-averaged streamwise velocity profile of the particle mixture is plotted in figure 3(a) for all the simulations. All the curves are superimposed meaning that the response of the granular medium to the fluid flow forcing is not modified by the number of small particles. The granular forcing is therefore the same for all the simulations independently of the number of small particles in the initial condition. In the quasi-static region (i.e. between z = 3 and 8 approximately), the linearity of the profiles in the semi-log plot indicates that the velocity is exponentially decreasing in the bed. Due to the presence of a fixed layer of particles at the bottom of the domain, the velocity goes to zero there. It is interesting to note that the entire bed is in motion, even if the velocity can be very low at the bottom of the quasi-static region. This is characteristic of a creeping flow. The time mean solid volume fraction of the mixture is plotted in figure 3(b) for all the simulations. The quasi-static region is characterized by the bed at random close packing (Φ ∼ 0.61). Above z ∼ 8, the packing fraction decreases to zero, corresponding to a transition zone between a quasi-static region and a pure fluid phase. Note that due to this decompaction of the bed at the surface, the total height of the bed is slightly larger than 10d l .
Since the small particles infiltrate the bed as a layer, the centre of mass of the small particles, z c , is a representative position of the entire layer. Figure 3(c) shows the temporal evolution of this position for all the simulations, in a semi-logarithmic plot. After the initial transient phase, the curves become linear, meaning that z c is a logarithmic function of time In (3.3), the coefficient a corresponds to the absolute slope of the curve and characterizes the segregation velocity (dz c /dt = −a/t). Whatever the number of small particles in the simulation, figure 3(c) shows that all the curves are parallel to one another, meaning that a is independent of the number of layers of small particles, N s . Therefore the segregation velocity dz c /dt is also independent of the number of layers of small particles. In the following, these simulations will be analysed using the dimensional analysis presented in the introduction (1.4) with the aim to confirm the dependence of the segregation flux on the inertial number I and on the local concentration φ s .

Dependence on the inertial number
In figure 4(a), the dimensionless segregation velocity is plotted against the large particle inertial number I =γ p d l / P p /ρ p , whereγ p is the time mean fluid shear rate and P p is the time mean lithostatic granular pressure. The segregation velocity is higher for larger inertial number. The linearity of the curves shows that the segregation velocity is indeed a power law of the inertial number. For all the simulations a similar exponent is obtained, ranging from 0.81 to 0.88, and a best fit gives a value of 0.85 for the mean exponent, dz c (t) dt ∝ I 0.85 (z c (t)).  Savage & Lun (1988) who suggested relating the segregation velocity to the shear rate. Indeed the inertial number is the dimensionless ratio between the shear rate and the square root of pressure. In the bedload configuration, the pressure increases linearly with depth while the shear rate exponentially decreases in the bed. Therefore, most of the variation of the inertial number is contained in the shear rate. A scaling with the inertial number allows the effect of pressure P p to be taken into account, which is small in this configuration, but the similarity with the results of Fry et al. (2018) is encouraging in the understanding of size segregation in general.
Equation (3.4) allows the temporal evolution of the fine particles centre of mass to be understood. Figure 4(b) shows the inertial number profile for all simulations. The linearity of the curves, for z 8, in the semi-log plot indicates that the inertial number is an exponential function of z in the quasi-static part of the bed. Taking an Modelling of grain size segregation during bedload transport 895 A30-13 leading to for long times. It can be rewritten as where b 2 = −c ln(b 1 /c) is a constant. This analysis shows that the logarithmic descent of the small particle centre of mass is a consequence of the dependence of the segregation velocity on the inertial number. This is confirmed by the comparison between the coefficients a in (3.3) and c in (3.5), that should be equal. Fitting of the elevation of the centre of mass (figure 3c) and of I 0.85 yields coefficients a and c (see table 1). The maximal difference between a and c is approximately 3 % confirming the present analysis and showing that the inertial number is indeed the controlling parameter.
3.4. Bottom controlled segregation While this analysis clearly explains the trends observed, figure 4(a) shows that for a given value of the inertial number, different segregation velocities are obtained depending on the initial number of small particles, N s .   Since the inertial number follows an exponential profile, it varies importantly throughout the layer of small particles. Thus according to (3.4) the segregation velocity should be lower at the bottom than at the top of the layer. This lower segregation velocity at the bottom of the layer implies that all the small particles above cannot move downward faster than the lowest particles in the layer. Defining the position of the bottom of the layer as z b = z c − W/2, with W = 2N s d s /d l the small particle layer thickness, figure 5(a) shows the dependence of the segregation velocity on the inertial number at the bottom of the layer. All the curves collapse on a master curve and a power law relationship is found between the segregation velocity and the inertial number at the bottom of the layer where α 0 is a positive constant independent of the number of layers of small particles. This result shows that the segregation velocity of the layer is indeed completely controlled by the inertial number at the bottom of the layer and it does not contain any dependence on the number of small particle layers, N s . The layer thickness has been chosen to be two times the thickness it would occupy if only small particles were present, W = 2N s d s /d l . This choice is motivated by the fact Modelling of grain size segregation during bedload transport 895 A30-15 that the layer is formed of a mixture of both large and small particles. Figure 5(a) shows that, for the different cases, the width obtained is indeed consistent with the actual thickness. Note that when N s tends to zero, z b tends to z c the centre of mass of the small particles. In the extreme case where only one small particle is present, this position corresponds to its centre. Figure 5(b) shows, at time t = 40 435, the profiles of small particle concentration, where the crosses denote the position of the bottom of the layer. The lower limbs of the concentration profiles are superimposed and the position of the bottom of the layer is identical for all tested values of N s . Whatever the number of small particles, they pile up above the bottom position, where the segregation dynamics is controlled.
It is remarkable that the continuum description (3.10) is applicable even for the cases N s = 0.01 and N s = 0.05, corresponding to very few non-interacting particles and for which a continuum vision seems at first hardly relevant. This can be explained by the bottom controlled dynamics. The isolated particles are placed at the bottom position which is a stable position. Since for cases N s = 0.01 and N s = 0.05, z b and z c are located at the same elevation, the continuum description captures the dynamics of isolated particles.
This analysis highlights the dependence of segregation on the inertial number in bedload transport. Furthermore, the bottom of the small particles layer is shown to be a key position that controls segregation. In particular, it explains why the small particles infiltrate the bed as a layer of constant thickness.

Size ratio influence
The analysis presented above provides a reference position that describes the segregation dynamics. In the following, the effect of the size ratio is investigated. The previous simulations were all performed at the same size ratio, r = 1.5. According to the dimensional analysis, the following relationship is expected: (3.11) A set of simulations in which the diameter of the small particles is varied has been performed. It covers a range of size ratios from r = 1.15 to 3. In all these simulations, the number of small particles layer is fixed to N s = 1. The dependence of the centre of mass velocity on the inertial number at the bottom of the layer of small particles is plotted in figure 6(a). A similar power law relationship is recovered, showing that the inertial number is still the controlling parameter. The best fit provides a 0.82 exponent, which is in the range of values found previously. However, to remain consistent with the previous analysis an exponent of 0.85 is used in what follows. The curves are parallel but not superimposed, which implies that, as expected, the size ratio plays an important role in the segregation dynamics. The dependence with the size ratio is presented in figure 6(b) where the ratio between the segregation velocity and the inertial number to the power 0.85 is plotted as a function of the size ratio minus one. The dependence with the size ratio tends to zero when r is close to unity, cancelling segregation in the monodisperse limit, and increases strongly when the size ratio increases. The following functional form for G 1 is proposed: The maximal efficiency at r = 2 reported by Golick & Daniels (2009), Thornton et al. (2012 and Guillard et al. (2016) is not recovered in this configuration as the segregation velocity still increases for r > 2. Further research is necessary to gain better understanding of the size ratio effect on segregation.
To summarize, the results presented here suggest that the segregation dynamics in bedload transport can be described, to leading order, by the following relationship: (3.13) where α 0 can be computed from figure 5(a) and G 1 (r) is given by (3.12).
3.6. Dependence on the local concentration φ s It is remarkable that the trends observed for the evolution of the centre of mass and for the segregation velocity are independent of the number of small particle layers, N s . Indeed, as the latter is increased from 0.01 to 2, the configuration moves from a few isolated (non-interacting) small particles to two consistent layers of small particles infiltrating collectively inside the coarse bed. Therefore, it is expected that the increase of N s , would increase the local particle concentration and reduce the segregation velocity due to particle hindrance effects.
This independence of the trend observed in the segregation velocity seems in contradiction with the literature (e.g. Savage & Lun 1988;Dolgunin & Ukolov 1995;Van der Vaart et al. 2015;Jones et al. 2018), in which a major influence of the local small particle concentration on the segregation velocity has been evidenced. In order to explain this apparent contradiction and get more insight into the physical processes at play, the results are discussed in the following as a function of local mechanisms within the framework of the continuum model of Gray & Thornton (2005), Thornton et al. (2006) and Gray & Chugunov (2006).
Modelling of grain size segregation during bedload transport 895 A30-17 4. Discussion in the framework of continuum modelling

Local mechanism interpretation
In order to explain the apparent independence of the results on the concentration of small particles, a local analysis in the theoretical framework of continuum modelling is considered. The model of Thornton et al. (2006) (1.6), presented briefly in the introduction, will be used as a baseline. In the present configuration, the velocity field has only a streamwise component (along x) and all quantities only depend on z. As a first step, diffusion is not considered and (1.6) simplifies to a purely hyperbolic equation, Combining the form proposed in the introduction (1.7), the analysis of the DEM simulations (3.13) and the inertial number functional form (3.5), the segregation flux can be written as with S r0 = α 0 G 1 (r)I 0 . The inertial number does not depend on the number of small particles (figure 4b) so S r0 is independent on the initial number of small particles and only depends on the response of the granular mixture to the fluid flow forcing through I 0 and c. For r = 1.5, the corresponding value of S r0 obtained from the DEM simulations is: S r0 = 3.70 × 10 −8 . Due to the exponential decrease of the inertial number in the bed, the form of the segregation flux, F s , is also exponential with z. May et al. (2010) already proposed an analytical solution for an exponential flux, but the flow was shear driven from the bottom. Proceeding similarly, the analytical solution can be derived using the method of characteristics (appendix A) for which the initial condition is simply a step of concentration representing an initial pure layer of small particles lying on a bed made of large particles, with z i = N l d l . The solution is plotted in figure 7(a) for the case N s = 1 and r = 1.5. The initial discontinuity leads to an expansion fan delimited by two characteristics (see appendix A). The first, denoted z 1 , is going down and meets the bottom of the domain at time t 1 , which is not reached during the simulations. The second is going up and is denoted z 2 . When it meets the top of the domain at time t 2 (i.e. when the first large particle reaches the top boundary, which happens very rapidly), a shock is created that travels downward. It can be shown (see appendix A) that the analytical solution for time t 2 t t 1 is where z 1 is the characteristic curve (4.5) delimiting the bottom of the rarefaction zone, z 2 is given by the shock (4.6) delimiting the top of the rarefaction zone and φ fan is the solution in the expansion fan given by (4.7): Despite the complexity of the solution (4.4), some simplifications can be made for long times (smaller than t 1 ). Indeed for t 4c/S r0 e z i /c ∼ 1000, it can be shown (see appendix A) that the square root term in (4.7) simplifies to 1 + S r0 t 2c 2 e (z+z i )/c ∼ S r0 t 2c 2 e (z+z i )/c , (4.8) and the volume fraction profile φ s in the rarefaction fan simplifies as Some simplifications can also be done on the boundary curves of the rarefaction zone z 1 and z 2 . Introducing the long time solution (4.9) in the shock equation of z 2 (4.6), a much more simple equation of z 2 is obtained for long times (4.10) for which the solution is z 2 (t) = −c ln(t) + β, (4.11) where β is an integration constant. Concerning the characteristic curve z 1 , the term in the logarithm in (4.5) simplifies for long times to Modelling of grain size segregation during bedload transport 895 A30-19 and the following expression is obtained for z 1 : (4.13) Therefore, the upper and lower bounds of the small particle layer, z 2 and z 1 , both move down asymptotically as −c ln(t), and the layer thickness z 2 (t) − z 1 (t) is constant. Lastly, the form of the flux considered in the continuum model F s ∝ I 0.85 φ s (1 − φ s ) implies that the vertical velocity of the small particles w s = −F s /φ s is (4.14) Using (4.9), the vertical velocity can be expressed for long times as (4.15) and remembering that I 0.85 (4.16) Considering long time solutions and an imposed forcing by the inertial number, the model exhibits the exact same behaviour observed in the DEM simulations. First, the small particles move down as a layer into the coarse particles, with the lower and upper bounds of the small particle layer showing a logarithmic decrease with time (4.11) and (4.13). This decrease, at a common slope −c related to the forcing, imposes a constant layer thickness with time, as observed in the DEM simulations. The consequence is that the segregation velocity of the small particles is the same throughout the layer (4.16). Considering the formulation of the downward velocity of small particles (4.14), this means that the product of the small particle concentration and the inertial number is constant within the small particle layer, (4.17) Therefore, at a given time, the shape of the small particle concentration profile responds directly to the variation of the inertial number as a function of z. At a point within the layer, the decrease of inertial number observed when going inside the bed is compensated for by an increase of the term 1 − φ s , and so a decrease of the small particle concentration. This competition between the effect of the inertial number and the local concentration results in the shape of solution profiles presented in figure 7(b) for which the expression is given by (4.9). In order to corroborate the explanation of the mechanisms at play, the lower bound of the small particle layer z 1 for which φ s = 0 is considered. It corresponds, in figure 7 This corresponds exactly to the observations in the DEM simulations in § 3.4, where the segregation velocity was shown to collapse as a function of the inertial number at a position considered to be the bottom of the small particle layer. The correspondence between the inertial number at the bottom position z b observed in DEM and the one calculated from the analytical continuum model z 1 (table 2) confirms that the mechanisms described here are indeed the ones at play in the DEM simulations. Finally, it is interesting to remark on figure 7(b), that all the profiles tend to zero at the same position, meaning that the bottom position is the same whatever the quantity of small particles. This confirms the DEM result that the segregation velocity is independent of N s . Indeed, if there are more small particles, they pile up above the others without changing the position of the bottom where the segregation dynamics is controlled.
The continuum model without diffusion enabled understanding of the local segregation mechanisms and gave insights into the underlying physical mechanisms. While this model gives a good agreement at first order, as shown in figure 8, it clearly lacks diffusion in order to reproduce quantitatively the DEM simulations.

Diffusion and upscaling in a continuous framework
The continuum model used in the previous subsection can be expanded to account for diffusion (Gray & Chugunov 2006), in order to quantitatively reproduce the DEM Modelling of grain size segregation during bedload transport 895 A30-21 results. The equation to be solved is given by with a segregation flux of the form F s = S r0 e z/c φ s (1 − φ s ), with S r0 = 3.70 × 10 −8 computed from the DEM simulations, and a diffusion coefficient D(z) which is supposed to depend only on the vertical coordinate z. In this paper, a long time asymptotical solution is presented. At long times, DEM simulations have shown that the shape of the concentration profile φ s is self-similar and is only advected downward like a travelling wave. Indeed, the layer of small particles moves down as a layer of constant thickness as −c ln(t) (see (3.9)) and the shape of the concentration profile does not evolve with time. The following change of variable is proposed to place the problem in the moving frame of the small particles: For a travelling wave, the solution to (4.21) should not depend on time τ . Assuming ∂φ s /∂τ = 0, this requires that τ D(ξ − c ln(τ )) is independent of τ . This implies ∂ ∂τ (τ D(ξ − c ln(τ ))) = 0, (4.22) and distributing the derivative with τ , the following differential equation is obtained: for which the solution is This shows that the only way to obtain a time-independent solution is that the diffusion coefficient has an exponential structure with the same vertical dependence as the advection flux F s . In other words, the diffusion coefficient should also be proportional to the inertial number to the power 0.85 in order to have a concentration profile of constant thickness. Under this assumption, equation (4.21) becomes (4.25) By integration, and imposing φ s → 0 when ξ → ±∞, the following equation is obtained: for which the solution is where C is a constant satisfying the mass conservation +∞ −∞ φ s dξ . The travelling wave solution has a complex dependency on ξ and contains an integral which cannot be computed analytically. A semi-analytical approach is developed where the integral is computed numerically and a numerical optimization is performed on the C coefficient so that +∞ −∞ φ s dξ corresponds to the initial mass of small particles introduced in the model. Figure 9 compares the travelling wave solution (4.27) with concentration profiles from the DEM simulations (case N s = 1, r = 1.5) at different times. The D 0 coefficient was computed in order to minimize the root mean square of the difference between the DEM and the travelling wave solution. The travelling wave solution agrees very well with all DEM profiles. The width and the height of the profile are correctly predicted. However the latter is slightly below the DEM profiles. This may be explained by the uncertainties when computing the value of S r0 from the DEM simulations.
In order to keep the thickness of the small particles layer constant, both the advection and the diffusion coefficients must have the same exponential vertical structure. Defining a Péclet number by the ratio between the segregation intensity and the diffusion coefficient, Pe = S r (z)/D(z) = S r0 /D 0 , the travelling wave solution leads to a constant Péclet number with depth. The relative effect of segregation and diffusion fluxes is constant and independent of z. At each depth, the balance between advection and diffusion has to be the same in order to be consistent with the constant layer thickness observed in the DEM simulations.
The value of the Péclet number depends on the number of layers of small particles N s and is plotted in figure 9(b). The value of Pe increases almost linearly with N s . Since the value of S r0 is the same whatever the number of layers of small particles, this result suggests that the value of the diffusion coefficient decreases with N s . This would indicate that diffusion is more efficient when there are fewer small particles. In addition, the mean value of the local concentration in small particles φ s globally increases with N s (see figure 5b). Therefore, the dependence of the diffusion coefficient with N s could be a global manifestation of a dependence of the diffusion coefficient with the local concentration φ s .
This study has shown that both the segregation flux and diffusion need to be proportional to the inertial number to the power 0.85 in order to represent the dynamics of segregation in bedload transport.

Conclusion
Vertical size segregation in bedload transport has been studied with a discrete and continuum approach. Focusing on the quasi-static region -common to any granular flow on a pile -it has been shown with DEM simulations that small particles infiltrate the coarse quasi-static bed with the same behaviour, irrespective of whether there are just a few isolated particles or two coherent layers of small particles. All the different configurations exhibit the same segregation velocity, related to a power law of the inertial number, generalizing the results observed for moderate inertial numbers in confined granular mixture without fluid (Fry et al. 2018). The dynamics of the infiltration of small particles is totally controlled by the bottom of the small particle layer, as the inertial number decreases exponentially from the top to the bottom of the layer. As a consequence, the small particles segregate down as a layer of constant thickness.
The continuum size-segregation model of Gray & Thornton (2005), Thornton et al. (2006) and Gray & Chugunov (2006) has been shown to reproduce quantitatively the DEM simulations, with a small particle layer of constant thickness segregating at a velocity independent of the thickness of the layer. In the continuum model, this comes from a perfect balance at any elevation between the effect of the inertial number and the local small particle concentration. This analysis demonstrates that the macroscopic segregation behaviour always results from a local equilibrium between the inertial number forcing and the local small particle concentration. Moreover, based on the derivation of an analytical solution with a travelling wave approach, the results show that the diffusion coefficient and the segregation flux in the continuum model should have the same dependency on the inertial number.
This paper represents an original contribution on segregation processes in bedload transport, and more generally in dense granular flows. This is a first step toward upscaling grain size segregation in continuum models for sediment transport. In the future, it would be interesting to use DEM simulations to infer the different graingrain forces involved in size segregation as well as collective effects related to particle concentration.