{'fraction_non_alphanumeric': 0.030775536973274308, 'fraction_numerical': 0.01824889326119036, 'mean_word_length': 4.668308550185873, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bridging the gap between particle-scale forces and continuum modelling of size segregation: application to bedload transport

Gravity-driven size segregation is important in mountain streams where a wide range of grain sizes are transported as bedload. More particularly, vertical size segregation is a multi-scale process that originates in interactions at the scale of particles with important morphological consequences on the reach scale. To address this issue, a volume-averaged multi-phase flow model for immersed bidisperse granular flows was developed based on an interparticle segregation force (Guillard et al. 2016) and a granular Stokesian drag force (Tripathi and Khakhar 2013). An advection-diffusion model was derived from this model yielding parametrisations for the advection and diffusion coefficients based on the interparticle interactions. This approach makes it possible to bridge the gap between grain-scale physics and continuum modelling. Both models were successfully tested against existing Discrete Element Model (DEM) simulations of size segregation in bedload transport (Chassagne et al. 2020). Through a detailed investigation of the granular forces, it is demonstrated that the observed scaling of the advection and diffusion coefficients with the inertial number can be explained by the granular drag force dependency on the viscosity. The drag coefficient was shown to be linearly dependent on the small particle concentration. The scaling relationship of the segregation force with the friction coefficient is confirmed and additional non-trivial dependencies including the inertial number and small particle concentration are identified. Lastly, adding a size ratio dependency in the segregation force perfectly reproduces the DEM results for a large range of small particle concentrations and size-ratios.


Introduction
Bedload transport, the coarser sediment load transported by the water flow in close contact with the mobile river bed, is a major process that shapes the Earth surface with consequences for public safety, water resources, territorial development and fluvial ecology. In mountain streams with steep slopes, large quantities of a wide range of grain sizes are transported leading to grain size sorting more generally named size segregation. Size segregation remains a poorly understood phenomenon (Gray 2018) impairing our ability to model the interplay between sediment transport rates and channel morphological evolution such as armouring (Bathurst 2007), bedload sheets Bacchi et al. 2014), patching  or downstream fining (Paola et al. 1992). The physics of granular media has been advocated to address segregation at the granular scale and understand geomorphological evolution (Frey & Church 2009. Size segregation largely originates from local interparticle interactions but has huge consequences on the particle size repartition both in the downward and streamwise directions over a much larger scale, potentially affecting sediment mobility and the entire channel geomorphological equilibrium (Gilbert & Murphy 1914;Ferguson et al. 2015;Dudill et al. 2017Dudill et al. , 2018. While investigating segregation at the granular scale (usually with discrete methods) is invaluable (Hill & Tan 2014;Ferdowsi et al. 2017;Chassagne et al. 2020), it is also necessary to consider continuum modelling to improve our theoretical understanding and to provide predictions at larger scales. The focus of this paper is therefore to bridge the gap between the granular scale processes and continuum modelling, by determining closures based on local granular mechanisms.
This contribution focuses on vertical size segregation processes due to kinetic sieving and associated squeeze expulsion (Savage & Lun 1988;Gray 2018). The moving particles act as a random fluctuating sieve, in which small particles are more likely to percolate under the action of gravity than larger particles. This downward movement is balanced by upward squeeze expulsion which equally applies on small and large particles resulting in a net downward motion of the small particles. The combination of both processes is called gravity driven segregation (Gray 2018) and is the dominant mechanism in bedload transport. Beyond the few studies made on size segregation in bedload transport (Hergault et al. 2010;Ferdowsi et al. 2017;Frey et al. 2020;Chassagne et al. 2020), these processes have been studied experimentally and numerically in many granular flows such as dry granular avalanches (Savage & Lun 1988;Dolgunin & Ukolov 1995;Wiederseiner et al. 2011;Jones et al. 2018; Thornton et al. 2006;Guillard et al. 2016), shear cells (Golick & Daniels 2009;van der Vaart et al. 2015) or annular rotating drums (Thomas 2000).
While particularly complex segregation phenomena were observed (Thomas 2000), size segregation has been found to be mainly related to the forcing, the size-ratio and the fine particle volume fraction. Savage & Lun (1988) predicted from dimensional analysis that the shear rateγ p should be the controlling parameter for size segregation. Indeed, when a granular medium is sheared, a layer of particles moves relatively faster than the one beneath, allowing particles to find gaps in which to fall by gravity. This theory agrees with experimental bidisperse flow down inclined planes (Savage & Lun 1988). In more recent works, Golick & Daniels (2009) with shear cell experiments, and Fry et al. (2018) with shear cell Discrete Element Model (DEM) simulations, evidenced the effect of granular pressure, observing less efficient segregation when increasing the pressure. Gray (2018) suggested size segregation to depend on the inertial number, classically used to describe granular rheology (GDR MiDi 2004;da Cruz et al. 2005), where d l is the large particle diameter,γ p is the granular shear rate, P p is the granular pressure, and ρ p is the particle density. DEM simulations of dry granular flows (Fry et al. 2018) and turbulent bedload transport (Chassagne et al. 2020) have shown that the segregation velocity indeed scales with the inertial number to a power 0.845 ± 0.05 from quasi-static to dense granular flow regimes. Not surprisingly, a number of studies have also found that the segregation depends on the particle size ratio in the kinetic sieving regime. As kinetic sieving is related to the gaps created by shearing, it appears logical to be related to the size ratio. However, while Chassagne et al. (2020) have found that segregation increases monotonically with the size ratio in quasi-static regimes (I < 10 −3 − 10 −2 ) for size ratio up to 3, Golick & Daniels (2009) and Guillard et al. (2016) found that it experiences a maximum efficiency for a size ratio of two, r = 2, for more dynamic granular regimes (10 −3 < I < 1). Yet, there is still no satisfying theory that explains this difference. Similarly to the hindrance function for the fluid drag force on a particle, size segregation has also been observed to depend on the concentration of fine particles. Indeed, studies indicate that the efficiency of the segregation process is linked to concentration in small (or large) particles in the granular sample (Fan et al. 2014;van der Vaart et al. 2015;Jones et al. 2018).
Size segregation can be analysed both from a particle-scale mechanistic point of view, or a continuum one. On the one hand, considering a large particle in a bath of small particles, size segregation can be seen as a force destabilizing the large particle and leading to a migration with respect to the small particles. A number of authors have adopted this approach and have shown that a particle experiences different kind of forces linked to size segregation (Ding et al. 2011;Tripathi & Khakhar 2013;Guillard et al. 2016;Staron 2018;van der Vaart et al. 2018). The forces can be decomposed into a component that drive the segregation and a resulting resisting component linked to the relative motion of the large particle with respect to the small. These two different forces have been isolated by Guillard et al. (2016) and Tripathi & Khakhar (2013). To assess the segregation forces due to the particle size differences, Guillard et al. (2016) performed 2D DEM simulations of a large disk placed in a bed of small disks in the simple shear flow configuration. Maintaining the large disk at a given position with a virtual spring, they were able to assess the vertical segregation force applied by the small particles to the large one, f seg (see figure 1), without generating a resisting force due to the particle motion. Guillard et al. (2016) found that the vertical segregation force has two contributions: one proportionnal to the pressure gradient ∂P p /∂z arising from the enduring contact between particles; the other proportional to the granular shear stress gradient ∂|τ p |/∂z: where V l = πd 3 l /6 is the volume of the intruder and F and G are empirical functions depending on the friction coefficient µ = |τ p |/P p and on the size ratio r = d l /d s between the intruder and the surrounding small particles. Guillard et al. (2016) studied the dependency on both parameters but only provided a dependency with µ as where µ c is the critical friction coefficient defining the threshold of movement. Tripathi & Khakhar (2013) performed 3D DEM simulations of a settling heavy sphere in a bed of lighter spheres, during a steady dry granular flow on an inclined plan. This density segregation setup generates a relative motion between the heavy sphere and the lighter ones, without generating segregation forces due to size ratio. By analogy with classical hydrodynamics, light particles playing the role of an ambiant fluid, the authors showed that the interaction force could be modeled with a Stokesian form of a solid drag force where v is the settling velocity of the heavier particle, c(Φ) is a drag coefficient depending on the local solid volume fraction Φ, and η p = |τ p |/|γ p | is the viscosity of the granular medium considered as a non-Newtonian fluid. Tripathi & Khakhar (2013) suggested that c(Φ) depends on the local volume fraction Φ but still remains of the same order as the value observed for a Stokes law in Newtonian fluids, i.e. c = 3. These two forces allow one to understand particles migration locally, and to relate the segregation behaviour of particles to the local characteristics of the granular flow.
By contrast, addressing the effect of size segregation processes at the large scale requires a different approach that disregard the particles. Such an approach has been extensively developed in the last few years focusing on a description of segregation as an advectiondiffusion model for the percolation of small particles (Dolgunin et al. 1998;Thornton et al. 2006;Gray & Chugunov 2006;van der Vaart et al. 2015;Ferdowsi et al. 2017;Gray 2018;Cai et al. 2019): where t denotes for time, z for the vertical axis, φ s and φ l are the small and large particle concentration and sum to unity (with φ s + φ l = 1), w s is the advection velocity of segregation and D is the diffusion coefficient. This equation is characterised by the segregation flux φ s w s , and the advective velocity w s , which encompass the physical dependencies of size segregation discussed previously. The advective velocity should therefore have a dependence on the local concentration which is classically taken as proportional to the large particle concentration, w s = φ l S r (Bridgwater et al. 1985;Savage & Lun 1988;Dolgunin & Ukolov 1995;Gray & Thornton 2005;Gajjar & Gray 2014;Fan et al. 2014;Jones et al. 2018). S r is called the advection coefficient and it has been usually taken as an empirical constant for a given application, or determined from semi-empirical analysis. Based on a dimensional analysis and DEM simulations, Chassagne et al. (2020) showed that it should depend on both the inertial number and the size ratio. The diffusion coefficient D models the diffusive remixing of small particles into large particles. Contrary to the advection coefficient, the diffusion coefficient has received less attention in the literature. It has been suggested that it should depend on the volume fraction (Cai et al. 2019) and on the inertial number (Chassagne et al. 2020). Developing a three-phase continuum mixture theory to model a bi-disperse combination of large and small particles with an interstitial passive fluid, Thornton et al. (2006) and Gray & Chugunov (2006) were able to analytically derive the advection-diffusion model (1.5). This derivation represented an important step in the understanding of the physical processes at work in size segregation since the advection and diffusion coefficients of the advection-diffusion equation were linked to the particle-scale interactions. In particular, the derivation is based on the assumption that the size segregation directly takes its origin in the heterogeneous distribution of the granular pressure between small and large particles. However, the form of the interaction forces between large and small particles have been postulated without support from independent physical evidence.
This literature review evidences the absence of direct link between the continuum modelling of grain-size segregation and the local segregation forces experienced by a grain. In this context, the aim of the present paper is to bridge the gap between the granular scale approach and the continuum modelling. Based on the particle-scale forces proposed by Guillard et al. (2016) and Tripathi & Khakhar (2013), a volumeaveraging approach (Jackson 1997(Jackson , 2000 is adopted here to derive a multi-phase continuum model from granular-scale forces. In addition to the novelty of the developed approach, the derivation proposed by Thornton et al. (2006) is used to express the advection-diffusion equation from the new multi-phase flow model, providing improved formulations of the advection and diffusion coefficients that contain the particle-scale granular dependencies. In order to test the proposed models, the bidisperse turbulent bedload transport configurations investigated in Chassagne et al. (2020) are used for comparison. The DEM simulations performed by the authors give a good reference in which granular-scale processes are explicitly resolved. In addition it makes it possible to focus on size segregation by providing an input for the granular rheology. As the study of Chassagne et al. (2020) focused on the quasi-static part of the bed in turbulent bedload transport, the comparison will be mainly performed in this regime. Since the fluid turbulence can be neglected in this regime (Maurin et al. 2016), it will not be taken into account in the derivation of the equations.
The paper is organised as follows: first, the forces acting at the granular scale for a single large intruder in an immersed sheared granular flow are discussed. Then, in section 3, the multi-phase flow model is derived by volume-averaging and the associated advection-diffusion equation is derived in section 4. Finally, both models are compared to the DEM simulations (section 5) and ways to improve the closures are discussed in section 6, including the influence of the size ratio.

A large intruder in a bath of small particles
As a first step, the force balance applied on a single large grain in an immersed granular medium made of smaller particles is presented. This Lagrangian equation of motion for the large intruder is then made dimensionless using classical scalings for granular flows, with the large particle diameter as the length scale. An order of magnitude analysis makes it possible to discriminate the most important forces for bedload transport application.

Force balance on the large intruder
The configuration is sketched in figure 1. The large particle is of diameter d l , of volume V l and of density ρ p in a bed of height h made of small particles. Below this layer the grains are in the quasi-static regime. Applying the Newton's second law, the vertical Lagrangian equation of the intruder can be expressed as In equation 2.1, the large intruder is submitted to five forces (see figure 1): its weight P = −ρ p V l g cos θ, the buoyancy force Π f = −ρ f V l g cos θ, the segregation force f seg , the drag force exerted by the fluid f f d and the drag force exerted by the small particles f p d . In the present configuration, the slope angle is low (tan θ = 0.1) and the streamwise gravity component is negligible, making the contribution of the shear stress gradient to the segregation force negligible (see equation 1.2 from Guillard et al. (2016)). Thus, the pressure gradient contribution is dominant and the segregation force can be simplified as where the function form F(µ) will be called the empirical segregation function.
While segregating at a velocity w l , the large particle is submitted to a fluid drag force. Because the particulate Reynolds number based on the vertical velocity Re p = d l ρ f w l /η f is very small in the bed, fluid inertial effects are negligible at the particle length scale and the vertical fluid drag force f d f may be approximated by the Stokes law (Stokes 1851) where η f is the fluid dynamics viscosity and w f is the vertical velocity of the fluid.
During its segregation motion, the intruder is also submitted to frictional forces from the surrounding small particles. This results in a particle drag force f p d modeled as proposed by Tripathi & Khakhar (2013) where w s is the vertical velocity of small particles and the drag coefficient c is first approximated as a constant equal to 3 (Tripathi & Khakhar 2013).

Dimensionless equation for the large intruder
In order to identify the dominant terms in equation (2.1), it is made dimensionless using classical scalings for granular flows: where k = s, l or f respectively for the surrounding small particles, the large intruder and the fluid. Introducing these variables in (2.1) and taking into account that cos θ ∼ 1, the dimensionless form of the large intruder Lagrangian equation of motion can be written as (2.6) Equation (2.6) contains two dimensionless numbers. The first one is the fluid Stokes number in which W = √ d l g is the characteristic velocity of the large particle. This Stokes number compares the inertia of the large intruder with the viscous friction exerted by the fluid. Similarly, the granular Stokes number compares the inertia of the large intruder with the contact friction exerted by the small particles in the vicinity of the intruder.
Assuming a classical bedload configuration, the fluid flows inside the porous matrix of the granular bed. Only the first layer of particles at the top is in a dense flow regime. Below this layer the grains are in the quasi-static regime. Typical values of the granular viscosity for dense granular flows are very high compared with the water viscosity (typical ranges span from 10 3 at the bed surface to 10 6 at the bed bottom). This results in a fluid Stokes number St f much larger than the granular Stokes number St p whatever the height into the bed. Therefore, the second term in the right hand side of equation (2.6), representing the fluid drag force, can be neglected. In addition, while the large intruder is rising, it only modifies the small particle bed structure locally. Therefore, it is assumed that the averaged vertical velocity of the small particles w s is negligible compared with w l . Focusing on the position of the intruder in the quasi-static part of the bed, it can be deduced from equation (2.6) that the total solid volume fraction is constant, Φ = cste. Therefore, for this configuration, a simple equation for the vertical velocity of the large intruder can be written as This dimensionless equation does not depend on the fluid parameters and should be valid to model the vertical velocity of an intruder segregating in a dry granular flow. In addition, with this equation, it can be shown from an asymptotical analysis that the formulation (1.3) of the empirical segregation function F(µ) does not satisfy equation (2.6) when small particles are not moving (see appendix A for more details). In order to satisfy this constraint, it is proposed to use the following functional form which has the same range of values than the original one but satisfies equation (2.6) for no flow condition. Equation (2.9) allows one to identify the main size segregation mechanisms and shows that the segregation of a large intruder can be seen as a simple relaxation process with characteristic time St p .

Volume averaged multi-phase flow model
As discussed in the previous section, the dynamics of a large intruder in an immersed granular flow made of small particles can be described using interparticle forces published in the literature. In this section, the goal is to upscale this result by volume-averaging the segregation forces over a collection of large particles in order to make the link between this discrete picture and continuum models for size segregation. This is done in the framework of the volume averaged equations from Jackson (1997Jackson ( , 2000 which provides continuum equations for the three phases: large particles, small particles and the interstitial fluid.

3D general governing equations
Following Jackson (1997Jackson ( , 2000, the mass and momentum balance equations for each class are given by where f is the fluid and indices i = l, s denote the large particle phase and the small particle phase respectively (δ = l if i = s and δ = s if i = l). Φ l and Φ s are the volume fractions for the large and small grains and verify Φ s + Φ l = Φ where Φ is the volume fraction of the mixture, i.e. the total solid volume fraction. Consequently, the fluid volume fraction is = 1 − Φ l − Φ s . S k is the stress tensor associated with phase k with k = l, s or f . They can be separated into pressure and shear stress contribution where τ k is the shear stress tensor and p k is the pressure of phase k. It should be noticed that for a solid phase, the static pressure arises from the enduring contacts between the particles. Thanks to the mixture model approaches (Morland 1992), it can be assumed that each particle phase carries the total overburden pressure p m according to their local volume fraction as where m denotes the mixture made of small and large particles. The total overburden pressure p m is computed using the formulation proposed by Johnson & Jackson (1987). For further information, the reader is referred to Chauchat et al. (2017) and Chauchat (2018).
The momentum equations (3.3) and (3.4) contain two terms coming from the momentum exchange between the different phases: n i f f →i and n i f δ→i . The term n i f f →i is the averaged value of the resultant forces exerted by the fluid on the particles of phase i. Jackson (2000) showed that for a collection of immersed particles, this interaction force can be written as ( 3.7) where Φ i ∇p f is the buoyancy force exerted by the fluid phase on the particles and n i f f →i d is the particle averaged viscous drag force between the particles and the fluid phase. The term n i f δ→i is the averaged value of all interacting forces between large and small particle phases. It can be directly expressed in 3D from the local segregation force of Tripathi & Khakhar (2013) and Guillard et al. (2016).
Therefore, the developed model is general and can be applied to 3D configurations. For simplicity and for the purpose of the present study, the model will be only developed for a 1D uniform flow.

Simplified 1D vertical multi-phase flow model
The multi-phase flow model (equations 3.1 to 3.4) is simplified by considering a uniform flow in the streamwise direction. From now, all the variables only depend on the vertical position z. Therefore, the spatially averaged velocity of the phase k can be written as u k = u k (z)e x + w k (z)e z . The mass conservation equations simplify to (3.8) and the momentum balance equations in the vertical direction are: (3.11) In the two last equations, the solid pressures p l and p s are given by equation (3.6). To solve these equations, it is necessary to prescribe closures for the spatially averaged fluid/grain interaction and grain-grain interactions, and for the granular and fluid pressures.
Considering the fluid-grain interaction, both small and large granular phases interact with the fluid phase through Φ i ∇p f and the drag force n i f f →i d . For an assembly of particles, the spatial averaging of the vertical total drag force applied by the fluid gives 3 /18η f is the particle response time and d i is the particle diameter of phase i. The factor (1 − Φ) 3 is a correction proposed by Richardson & Zaki (1954) to take into account hindrance effects. Since the drag is linear, the spatial averaging is simply the drag force applied on one particle (given in (2.3)) multiplied by the number of particles per unit volume (Jackson 2000), where V i is the volume of a single large particle.
The granular phases also interact with each other and the grain-grain interaction closure should be prescribed in the model. For a single large grain in a bath of small particles, it has been shown in section 2 that small particles exert two forces on a large intruder, To extend these forces to a collection of large particles, the interaction force f s→l is spatially averaged. Since this force is linear, it amounts to multiplying f s→l by the number of large particles per unit volume n l = 6Φ l /πd 3 l . Therefore, the total solid interaction force exerted by the small particles on the large ones is given by where t ls = ρ p d 2 l /6cη p is the particle response time for the drag force between small and large particles and F(µ) = 1 − e −70(µ−µc) (see appendix B). According to the Newton's third law, the force exerted by the large particles on the small ones (equation 3.11) is ( 3.15) The solid mixture phase is made of both particle phases and is noted with i = m. Its momentum balance is obtained by summing (3.10) and (3.11). Since the mixture phase does not distinguish between small and large particles, the solid interaction forces should not appear in this equation. Equation (3.15) ensures that these forces vanish when developing the mixture momentum equation.
The proposed volume-averaged multi-phase flow model describes size segregation of a bi-disperse mixture immersed in a fluid. This represents an improvement upon the model of Thornton et al. (2006), which was based on semi-empirical parametrisation of the interparticle forces between small and large particles. The present model provides closures based on forces applied on a single particle, and bridges the gap between granular scale processes and continuum modelling in size segregation. This important result will be used in the following to derive an advection-diffusion model for size segregation.

Derivation of the advection-diffusion model
A classical continuum approach to model size segregation is the advection-diffusion model (Dolgunin & Ukolov 1995;Gray & Thornton 2005). These models can be derived from the multicomponent mixture theory (Thornton et al. 2006; by substituting the percolation velocity of one particle-size into the mass conservation equation. The advection and diffusion coefficients can be modelled using experimental and theoretical closures (Dolgunin et al. 1998;van der Vaart et al. 2015;Ferdowsi et al. 2017;Cai et al. 2019) or can be derived as a simplification from the continuum model of Thornton et al. (2006) and Gray & Chugunov (2006). In the present section, the multi-phase model developed in the previous section (equations 3.9 to 3.11) makes it possible to derive an advection-diffusion model similar to Thornton et al. (2006) and Gray & Chugunov (2006), with advection and diffusion coefficients depending on the segregation and the drag forces (Guillard et al. 2016;Tripathi & Khakhar 2013) determined in independent idealised configurations.
Combining equations (3.11), (3.12), (3.14) and (3.15), the momentum balance of small particles can be written as: The total volume fraction Φ = Φ s +Φ l is assumed to be constant and equal to Φ max = 0.61 since particle velocity fluctuations are small. For a deposited bed, the particle momentum balance in the wall-normal direction reduces to a hydrostatic pressure distribution for both the fluid and the particle phases (Chauchat 2018): Assuming a constant mixture solid phase volume fraction, the pressure gradient can be integrated to give the pressure distributions: Following Thornton et al. (2006), the volume fraction per unit granular volume is introduced as φ i = Φ i /Φ. This notation is more convenient since it ensures φ s + φ l = 1. Using equation (3.6), the momentum equation (4.1) for small particles is rewritten as follows: Using the same scalings as in the Lagrangian equation (2.6) for a single intruder, the equation (4.4) is made dimensionless as follows: As shown in section 2, St f >> St p in the bed and the fluid drag force can be neglected. Furthermore, assuming a quasi-steady state and neglecting inertial terms, equation (4.5) can be rewritten as from which, assuming w m = 0 (Thornton et al. 2006), the flux of small particles can be expressed as (4.7) Equation (4.7) is then substituted in the mass conservation equation (3.8) to obtain the following advection-diffusion equation for the percolation of small particles: Since the pressure gradient is negative, S r is negative which ensures a downward flux for the small particle phase. In equation (4.8), S r is the segregation number or advection coefficient and D is the diffusion coefficient. Here, equation (4.9) provides physical closures, which are directly obtained from volume-averaging of the particle scale segregation forces. The advection coefficient S r is therefore expressed as a product between the segregation term F(µ)∂p m /∂z, which quantifies the ability for the small particles to fall downward, and the solid Stokes number which quantifies the drag force exerted by the other grains counteracting this downward movement. It is interesting to note that the granular Stokes number is present in both coefficients which indicates that it is a key parameter for the advection and the diffusive remixing. Equation (4.9) shows that the mixture pressure is also an important parameter for diffusion.
These results improve upon the original model from Thornton et al. (2006) and Gray & Chugunov (2006) since the experimentally-based closures (Guillard et al. 2016;Tripathi & Khakhar 2013) make it possible to link both the advection and the diffusion coefficients to the local physical parameters of the granular flow. This result not only provides closures for the advection-diffusion model based on local granular forces but also highlights the key local physical mechanisms controlling segregation and diffusion. Since the advection and diffusion coefficients do not contain fluid parameters, it indicates that even though the equation has been set for a bedload configuration, it is also valid for dry granular flows.
In the following, the relevance of the obtained advection-diffusion model with respect to the multi-phase flow model will be tested.

Comparison with existing discrete numerical simulations
In this section, the multi-phase flow model presented in section 3 and the corresponding advection-diffusion model presented in section 4 is tested against the discrete numerical simulations of Chassagne et al. (2020). This study focused on the segregation of small particles initially resting on top of large ones. It provides a comprehensive dataset to evaluate local granular parameters such as the volume fraction and the segregation velocities. In addition, since it is not the purpose of this work to develop a granular rheology, the shear stress and shear rate profiles obtained from the DEM will be used as input parameters for the continuum models. The 3D DEM configuration and the main results from Chassagne et al. (2020) are first summarised (section 5.1). Then, in section 5.2, the multi-phase flow model is compared with the DEM results using default parameters (see sections 3 and 4) for the segregation of the small particles. Finally, in section 5.3, the results predicted by the advectiondiffusion model and the multi-phase flow model are compared to determine the validity of the former with respect to the latter. Chassagne et al. (2020) In this section, the configuration explored and the main results of Chassagne et al. (2020) are briefly presented. The authors investigated grain-size segregation in turbulent bedload transport using a coupled fluid-Discrete Element Model (DEM) originally developed by Maurin et al. (2015). For further details, the interested reader is referred to Chassagne et al. (2020). The 3D bi-periodic DEM set-up consisted in depositing a layer of small particles over large ones, and letting the particles entrained by the fluid flow at a fixed Shields number. The latter is the dimensionless fluid bed shear stress Sh = τ f /[(ρ p − ρ f )gd l ] and was taken equal to 0.1. The bed slope was fixed to 10%, which is representative of mountain streams. The size ratio was taken as r = 1.5 with small particles of diameter d s = 4mm and large particles of diameter d l = 6mm. The amounts of large and small particles were assimilated to a number of layers, N l and N s . The number of layers of a given class represents the height, in terms of particle diameter of this class, occupied by particles if the concentration was equal to the random close packing (Φ max = 0.61). In this way, the bed height at rest was defined as h = N l d l + N s d s and was fixed to h = 10d l with a random close packing volume fraction Φ max = 0.61 (profile in figure 2a). In the study of Chassagne et al. (2020), different simulations have been performed with N s varying from 0.01 (a few isolated particles) to N s = 2. In this section it was decided that the comparison would be made with N s = 1.5. The bulk response of the granular mixture to this fluid forcing is represented by the dimensionless mixture streamwise velocity profile in figure 2a. The inset is a semilog plot of the dimensionless velocity profile and shows that it is exponential. As shown in figure 2b the linearity of the curve in the semilog plot confirms that the shear rate is exponentially decreasing in the quasi-static part of the bed (delimited by the two horizontal black dashed lines). As expected for a uniform flow, the mixture shear stressτ m xz shown in figure 2c is linear with depth. For both quantities, the following fits were proposed and plotted as a red dotted line in figure 2.γ p = γ 0 ez /s0 andτ p xz = a 0z + τ 0 (5.1)

DEM investigation of
The simulations performed by Chassagne et al. (2020) on this configuration were focused on the downward segregation of small particles. It was observed that the layer of small particles percolates rapidly forz > 8.5 (flowing layer) and then slows down below (this limit is marked in figures 2b and 2c). Chassagne et al. (2020) showed that the small particles are advected downward like a travelling wave into the bed made of large particle with a layer of constant thickness. As figure 2e shows, the small particle concentration has a Gaussian-like shape and remains self-similar in time while segregating. The center of mass of small particles,z c , is therefore representative of the dynamics of the entire layer. Chassagne et al. (2020) observed that the small particle layer travels down as a logarithmic function of time:z where a 1 is a constant characterising the segregation velocity (dz c (t)/dt = −a 1 /t). The authors demonstrated that this logarithmic descent of small particles is a consequence of the dependency of the segregation velocity on the inertial number as where S r0 = 0.049. Then, with the help of a travelling wave method, they evidenced that the small particles percolate as a layer and with a self-similar concentration profile because the ratio between the advection coefficient and the diffusion coefficient is constant. The Péclet number reads and is therefore constant withz, so that the diffusion coefficient has to have the same dependency on the inertial number as the segregation coefficient: where D 0 is taken as D 0 = 0.01. This work also demonstrated that the dynamics of the fine particle layer is controlled by its bottom position, which acts as a lower bound for the segregation velocity. In this way, the particles in the layer cannot segregate faster.
The numerical resolution of the 1D multi phase model apriori requires to solve the granular rheology in order to estimate the granular viscosity required to evaluate the granular drag force contribution. The goal of the present study is to focus on grain-size segregation modelling. In addition, the results of Chassagne et al. (2020) were obtained in the quasi-static regime, on which there is still no concensus regarding granular rheology. For these reasons, the granular viscosity is directly determined by the DEM results here. This makes it possible to focus on the effect of the segregation model and to put aside potential discrepancies linked to a non-accurate description of the granular rheology. The granular viscosity is therefore computed from the DEM results using the definition: As shown by the red dashed lines in figures 2b and 2c, the fits presented in equation 5.1 match the DEM results in the region 8.5 >z > 3. Thus, using these expressions in eq (5.7) provides an accurate estimate of the granular viscosity for the particle-particle drag closure. For this reason, the validation of the multi-phase flow model will only be carried out in this part of the bed. Therefore, the initial state consists in placing the small particles in the upper limit of the quasi-static part with the center of massz 0 c = 8.5. The sketch of this configuration is shown in figure 2d. Figure 2e shows the small particle concentration profile of this numerical set-up at the initial state. The initial concentration is taken with a gaussian fit on the DEM initial concentration (figure 2e), and ensures that the mass of particles is the same in the DEM and in the continuum simulations.

Comparison with the multi-phase flow model
The system of partial differential equations (3.8)-(3.11) is solved numerically for the configuration shown in figure 2d with the initial concentration profile of figure 2e. In these equations, the empirical segregation function F(µ) is the one proposed in equation (2.10) and the drag coefficient c is equal to 3 as suggested by Tripathi & Khakhar (2013). Because the fluid is incompressible, there is no equation of state for the fluid pressure. Nevertheless, remembering that + Φ = 1 and defining the volume averaged velocities, w = w f + Φw m , it can be demonstrated that the particle-fluid mixture is incompressible. A PISO (Pressure Implicit with Splitting of Operators) algorithm classically developed for incompressible Navier-Stokes equations is used to solve the pressure-velocity coupling. As the fluid pressure p f is the sum of the hydrostatic pressure and of the excess pore presure p f = p f + ρ f gz, the PISO algorithm consists in solving the momentum balance equations without p f in a predictor step. Then, using the predicted velocity fields, a Poisson equation is solved to find p f . Once the pressure is found, the velocity fields are corrected. This kind of algorithm has already been used to model sediment transport in Chauchat et al. (2017) and Chauchat (2018). Figure 3 shows the results of the spatio-temporal evolution of the small particle concentration for the multi-phase flow model (figure 3a) and for the DEM (figure 3b). First, it can be seen that the dynamics predicted by the multi-phase flow model is similar to the DEM. The position of the bottom of the layer is about the same in both cases. More quantitatively, the center of massz c of the small particle layer as a function of time is compared with the DEM in figure 4a. After a first transient phase (t > 1 × 10 3 ), the center of mass position is linear in the semilog plot, indicating that the logarithmic descent observed in the DEM simulation is well reproduced by the multi-phase flow model. The slope of the curve, representing coefficient a 1 of equation (5.2) is 0.68 in the DEM simulation and 0.49 for the multi-phase flow model, corresponding to an error of 28%. In addition, figure 4b shows that in both models, the bottom of the layer is positioned at the same depth indicating that the multi-phase model reproduces well the bottom controlled behaviour observed by Chassagne et al. (2020) with DEM simulations. However the Gaussian-like profile is not reproduced by the multi-phase flow model and a wider profile is obtained. In figure 3a the maximum concentration max(φ s ) (indicated by in figure 4b) is almost two times smaller than the one predicted by the DEM simulation, while the extent of the small particle layer is much larger. These results indicate that, with the current parametrisation the multi-phase flow model is relevant to qualitatively predict segregation dynamics. However, the error on the segregation velocity and the discrepancies on the concentration profile clearly show that the model needs to be improved.

Evaluation of the advection-diffusion model
In order to determine the ability of the advection-diffusion model to reproduce the same results as the multi-phase flow model, equation (4.8) is solved numerically. The resolution strategy is based on a conservative Godunov schemes where a no flux condition is applied at the bottom and on top of the vertical domain. A vertical discretisation of dz = h/80 is taken and the time step is computed in order to satisfy the CFL condition. The initial solution is the same than in the multi phase flow (figure 2e).
The numerical solution at timet = 60000 is plotted in figure 4. Both the center of mass (figure 4a) and the concentration profile (figure 4b) are almost superimposed with the multi-phase flow model solution, only slightly differing by numerical diffusion. Therefore, both models can be considered as strictly equivalent, meaning that the physical behaviour of the segregation forces added in the momentum equation of the small particles is accurately predicted by the single advection and diffusion coefficients of equation (4.8). Moreover, when deriving the advection-diffusion equation, it was assumed that the mixture volume fraction Φ = cste, that the vertical acceleration of the small particles, the vertical mixture velocity and that the fluid drag were negligible. The strong agreement between the models demonstrate that these assumptions are valid.
This new equation represents an important step in the upscaling since the behaviour of small particles can be predicted without solving the entire multi-phase flow model, providing a speed-up of one thousand for the numerical resolution of the equations. In the light of this result, the bidisperse segregation problem can be simply viewed as a competition between the advection coefficient S r and the diffusion coefficient D of small particles.

Discussion
The multi-phase flow model and the associated advection-diffusion model are able to reproduce qualitatively the DEM simulations and the main properties of segregation in bedload transport obtained by Chassagne et al. (2020) (bottom controlled segregation, logarithmic descent of the small particles). However, the segregation velocity is lower than in the DEM simulation and the shape of the small particle concentration profile is not adequatly reproduced. So far, the inter-particle drag force and the segregation force from Guillard et al. (2016) and Tripathi & Khakhar (2013) have not been modified. Yet, these forces were proposed in very different configurations. Therefore, in section 6.1, a sensitivity analysis to the empirical segregation function is presented. Then, in section 6.2, new formulations of the empirical segregation function F and of the drag coefficient c will be proposed based on the DEM results. As the advection-diffusion equation results are strictly identical with the multi-phase flow equations, the discussion and the associated simulations will only be performed with the advection-diffusion equation.

Investigation of the empirical segregation function F
Numerical data from Guillard et al. (2016) for the empirical segrgeation function F(µ) exhibit a significant scatter and it is possible to show that a constant function could also fit the data (see appendix A). In this section, an analysis of the sensitivity to the empirical segregation function, taken as constant and varying from F = 1 to F = 15, is presented. Figure 5a shows the temporal evolution of the small particle center of mass for the different values of the empirical segregation function F. The linear evolution in the semi-logarithmic plot is conserved with the same slope (coefficient a 1 in equation 5.2), whatever the value of the empirical segregation function, meaning that the segregation velocity dz c /dt is not modified when changing the empirical segregation function. Increasing F only makes the curves shift vertically. Figure 5b shows that the maximum concentration is better predicted when increasing the empirical segregation function. When reaching F = 15 agreement with DEM results is perfect. The previous simulation, where the empirical segregation function F is a function of the friction coefficient µ, is also plotted in these figures. It is interesting to note that simulation with F = 1 is almost superimposed with the one obtained with F(µ). These observations tend to show that the friction coefficient dependency has a small influence on the size segregation configuration investigated. Therefore, taking F as constant is a good approximation, at first order. Chassagne et al. (2020) showed that the advection coefficient S r was a function of the inertial number I (see equation 5.4) and linked the logarithmic descent to the exponential form of I in S r . In the proposed model, the advection coefficient S r = FSt p ∂p m /∂z (equation 4.9). In this coefficient, since the solid pressure gradient is constant in the bed (equation 4.2) and F = cste, there is only one non-constant variable which is the inverse of the granular viscosity 1/η p appearing in the granular Stokes number (equation 2.8). As figure 5a shows that the logarithmic descent is still preserved with this parametrisation, it demonstrates that the viscosity profile is mainly responsible for the logarithmic descent. Therefore, in the proposed advection coefficient S r , the empirical segregation function F controls the strength of the segregation force and can be seen as a forcing parameter. In contrast, the velocity at which the small grains are descending is controlled by the granular viscosity in the granular Stokes number St p . In this way, the segregation problem can simply be seen as the settling of small particles under a force proportional to the empirical segregation function F, into a complex fluid having a variable viscosity.
It should be noted that there is a direct relation between the dimensionless granular viscosity and the inertial number (see appendix C) written whereη p = η p /ρ p d l W is the dimensionless granular viscosity and µ =τ m xz /p m is the friction coefficient. Since the variation of µ and √p m is small compared with the exponential profile ofη p , the inversely proportional relation between the granular viscosity and the inertial number shows that the inertial number dependency found by Chassagne et al. (2020) can be seen as a dependency on the inverse granular viscosity, confirming the important role of the granular viscosity in the segregation dynamics. Figure 5b showed that a better maximum concentration was predicted with F = 15. However, the value of the empirical segregation function F should have no influence on the diffusion as D = φ s /Φp m St p . This influence can be explained by the Péclet number P e = S r /D. Indeed, when F increases, the advection coefficient S r increases as well and makes the Péclet number higher. As a result the diffusive effect becomes small compared with the advection and thus, the small particles stay more concentrated. However, the dynamics of the center of mass shows a better agreement for F = 1 (see figure 5a) than for F = 15, meaning that the diffusion coefficient needs to be improved.
The advection coefficient S r and the diffusion coefficient D have been plotted in figures 6a and 6b with F = 1, for φ s = max(φ s ). Forz > 7, S r is close to the values predicted by the DEM. Under this limit, both curves are exponentially decreasing into the bed but with a different slope leading to discrepancies. This shows that the 1/η p dependency has a fundamental role in the vertical structure of the advection coefficient S r . Yet, another dependency with depth is probably missing in the empirical segregation function or in the drag coefficient c to find again a similar slope to the DEM. Surprisingly, the proposed diffusion coefficient has the same slope as the one predicted by the DEM, which means that it contains the correct depth dependency. However, its value is too high by a factor ten, explaining why the advection-diffusion results are too diffusive.

Missing dependencies in the particle-scale forces
In the present paper, the advection and diffusion coefficients (equation 4.9) have been derived from particle-scale segregation and solid drag forces of Guillard et al. (2016) and Tripathi & Khakhar (2013). However, figure 6 shows that both coefficients should be improved so as to match the DEM results. The segregation and solid drag forces of Guillard et al. (2016) and Tripathi & Khakhar (2013) have been established in idealised granular segregation configurations (e.g. unique intruder, simple forcing, 2D DEM) so that one can expect the two formulations to miss some dependencies when considering more general cases (mixture of small and large particles, 3D modelling or complex forcing). Such dependencies probably lie in the drag coefficient c and the empirical segregation function F contained in the advection and diffusion coefficients. While the drag coefficient c was taken constant in their study, Tripathi & Khakhar (2013) suggested that it could depend on the local concentration of particles. Similarly to an hindrance function in a fluid flow, one indeed expects an increase of the effective solid drag force on a particle with increasing concentration. This dependency of the drag coefficient in the local particle concentration should also impact the diffusion coefficient profile (see figure  6b) and makes it possible to match the DEM. In addition, to correct the slope of the advection coefficient S r (see figure 6a), only the empirical segregation function F should vary with depth.
As detailed in section 5.1, Chassagne et al. (2020) have been able to express the advection and diffusion coefficients dependencies on the inertial number I (see equations 5.4 and 5.6). In the following, both DEM and advection-diffusion coefficients are compared so as to extract the potential missing dependencies of F and c from the DEM coefficients and to propose new formulations of these parameters. Then, it is verified that the results from the advection-diffusion model are consistent when using these proposed closures.
In order to compare the advection and diffusion coefficients to the DEM and to find the missing dependencies, it is first shown that the coefficients of equation (4.9) can be expressed with an inertial number dependency as in the DEM. Indeed, as already shown in equation (6.1), the granular Stokes number can be rewritten as a function of the inertial number I: With this new definition, the advection and diffusion coefficients obtained in equation (4.9) can be rewritten as a function of the inertial number I: √p m ∂p m ∂z and D = φ s √p m I Φ6cµ . (6.3) One can notice that √p m I =γ, which makes the diffusion coefficient directly proportional to the shear rate. Such a dependency for the diffusion coefficient was found in a size bidisperse case by Cai et al. (2019) with DEM simulations. This shows that the particleparticle force based model is able to explain a dependency found in a different experiment. It represents a powerful argument to show the robustness of the proposed model. Yet, the power 0.85 highlighted by the DEM does not seem to appear. In figure 7 the profile of √p m I/µ, the profile of I 0.85 from the DEM and from exponential fitting are plotted. One can observe that values only differ by a factor D 0 /I 0 = 1.13 and that the exponential evolution with depth is the same, which proves that in this case Assuming that the drag coefficient includes the accurate dependencies, the following equality should be obtained between diffusion coefficients of the DEM and the one proposed in equation ( As a consequence of this linear scaling in concentration, the dependency on the small particle volume fraction vanishes in the diffusion coefficient (equation 4.9). This modification will tend to smooth out the small particle concentration profile. In this drag coefficient, when φ s → 0 (i.e. one small particle in a bath of large particles), the drag coefficient vanishes while it should reach a constant value, c = 3 as shown by Tripathi & Khakhar (2013). To ensure a consistent formulation, it is therefore proposed c(φ s ) to read c(φ s ) = 28φ s + 3, (6.7) which tends to 3 when φ s → 0 and to 31 when φ s → 1.
As mentioned in the last section, the empirical segregation function F is expected to depend on depth (see figure 6a). In this way, the advection coefficient should correspond to the DEM and it follows From this equation, the empirical segregation function is easily expressed as Equation (6.9) makes it possible to express the missing dependencies in the empirical segregation function F. One can note that a dependency with the friction coefficient is found as predicted by Guillard et al. (2016). In addition, the empirical segregation function is found to depend on the small particle concentration. The dependency on √p m I −0.15 exhibits an additional more complex mechanism in the segregation force.
With the new formulations of the solid drag coefficient c(φ s ) and the empirical segregation function F(µ,p m , φ s ), it is verified that the DEM results can be reproduced. This is done with simulations accounting for different initial quantities of small particles (N s = 0.5, 1, 1.5, 2). Figure 8 shows the results for the time evolution of the center of mass and the final concentration profile. Both the mass centers and the final concentration profiles are fairly well superimposed with the DEM results. Therefore, the new parametrisation is consistent with the DEM simulations. The concentration profile with the original parametrisations of c and F(µ) (see figure 4b) has also been plotted in figure 8 (N s = 1.5). The shape of the final concentration profile has drastically changed, from a bell-shape to the expected Gaussian-like shape. This is attributed to the small particle concentration dependency in the drag coefficient, which cancels the original concentration dependency in the diffusion coefficient.
Lastly, note that changing the parametrisation of the forces still yields physical solutions which proves that the advection-diffusion model and the corresponding multiphase flow model are physically consistent and robust. 6.3. Influence of the size ratio Guillard et al. (2016) showed that the segregation force depends on the size ratio r = d l /d s and exhibits a maximum for a size ratio of r = 2. Based on DEM simulations, Chassagne et al. (2020) also studied the size-ratio dependency and found the segregation velocity of the percolating small particles to be a monotonic increasing function of r. Best fit of the DEM results suggested the following dependency f (r) = 0.45 e (r−1)/1.59 − 1 for the advection coefficient S r . This dependency is introduced into the empirical segregation function F as follows Using this parametrisation, simulations have been performed for r = 1.25, 1.5, 1.75, 2, 2.25. The results are plotted in figure 9 and compared with DEM simulations. For each case, the center of mass position is in very good agreement with the DEM simulations. For the lower size ratio, the concentration profiles are superimposed with the DEM results while for the higher size ratios (r = 2 and r = 2.25) the maximum concentrations are slightly higher than the DEM. This indicates that the model is not diffusive enough. Indeed, the size ratio dependency has only been introduced in the advection coefficient. As shown in section 6.1, the shape of the concentration profile results from a subtle balance between advection and diffusion through the Péclet number P e. For the highest size ratio, this balance is not perfectly reproduced by the proposed model, which indicates that the diffusion coefficient should also depend on the size ratio. This would imply that the granular Stokes number also depends on the size ratio. It could explain why Guillard et al. (2016) found a maximum segregation force for a size ratio r = 2, while Chassagne et al. (2020) found the advection coefficient S r to increase exponentially with the size ratio. Further research is needed to elucidate this point through a detailed investigation of the granular Stokes number dependency on the size ratio.

Conclusion
In this contribution size segregation in bidisperse systems has been investigated with special focus on bedload transport. The originality of the work presented herein is to propose a new multiphase flow approach, derived from a volume averaging technique, based on the most recent advances on particle-particle forces, namely the segregation force or buoyancy force from Guillard et al. (2016) and the drag force from Tripathi & Khakhar (2013). The proposed multiphase flow model formulation is very general and it can be applied to any immersed granular flow configuration. In a subsequent step, following the same procedure as in Thornton et al. (2006), an advection-diffusion model is derived from the multi-phase flow equations. This derivation makes it possible to identify the dependencies of the advection and diffusion coefficients with the local physical parameters of the flow such as the volume fraction of small particles, the mixture granular pressure and its gradient, the granular Stokes number and the segregation parameter.
Both models have been tested against the Discrete Element Model simulations of Chassagne et al. (2020) for bidisperse turbulent bed load transport. Without any tuning of the forces from Guillard et al. (2016) and Tripathi & Khakhar (2013), both continuum models qualitatively reproduce the main features of size segregation. This demonstrates that the scaling of the advection coefficient with the inertial number observed by Fry et al. (2018) and Chassagne et al. (2020) can be explained thanks to the dependency of the advection coefficient on the granular Stokes number and the underlying presence of the granular viscosity. Using the discrete element simulation results, improved parametrisations for the advection and diffusion coefficients have been proposed. They suggest that the empirical segregation function from Guillard et al. (2016) and the drag coefficient from Tripathi & Khakhar (2013) should incorporate a dependency on the small particle concentration. At last, the influence of the size ratio has been investigated and a dependency of the segregation function on the size ratio has been proposed.
In terms of perspectives for granular flows, the continuum models proposed herein are very general and should as well apply to dense dry granular flows. These models represent a general framework for developing and testing improved parametrisations for the segregation and granular drag forces in different flow configurations. Further work is needed to identify and propose more robust concentration, depth and size ratio dependencies of the empirical dimensionless coefficients appearing in both granular forces. Concerning the upscaling of size segregation processes in sediment transport applications two routes are opened. The first one consists in implementing the proposed multi-phase flow model in a 3D numerical model, such as sedFOAM (Chauchat et al. 2017), for the simulation of size segregation in complex sediment transport application such as riverbed armouring (Frey & Church 2009), scour around an hydraulic structure (Nagel et al. 2020)  The new empirical segregation function has to fill the gap with the advection coefficient found using DEM. Therefore, F is such that: In equation (D 4), the dimensionless pressure gradient is constant and, using the hydrostatic approximation, reads Introducing it in (D 4), it implies that the empirical segregation function must satisfies