# Discrete and continuum modelling of grain size segregation during bedload transport

Rémi Chassagne, Raphaël Maurin, Julien Chauchat, J. Gray, Philippe Frey

<figure><img src="image_1.png" /></figure>

## To cite this version:

Rémi Chassagne, Raphaël Maurin, Julien Chauchat, J. Gray, Philippe Frey. Discrete and continuum modelling of grain size segregation during bedload transport. Journal of Fluid Mechanics, 2020, 895, pp.365 - 398. 10.1017/jfm.2020.274. hal-02951912

HAL Id: hal-02951912

https://hal.science/hal-02951912

Submitted on 29 Sep 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Discrete and continuum modelling of grain size segregation during bedload transport

Rémi Chassagne¹·†, Raphaël Maurin³, Julien Chauchat², J. M. N. T. Gray⁴
and Philippe Frey¹

¹Univ. Grenoble Alpes, INRAE, IRSTEA, UR ETNA, 2 rue de la Papeterie-BP 76,
F-38402 St-Martin-d'Héres, France

²Univ. Grenoble Alpes, LEGI, CNRS UMR 5519 - 38000 Grenoble, France

³Institut de Mécanique des Fluides de Toulouse, IMFT, Université de Toulouse, CNRS - 31400 Toulouse, France

⁴Department of Mathematics and Manchester Centre for Nonlinear Dynamics, University of Manchester, Manchester M13 9PL, UK

(Received 30 August 2019; revised 24 March 2020; accepted 1 April 2020)

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.

**Key words:** granular media, sediment transport, multiphase flow

# 1. Introduction

Bedload transport, i.e. the coarser sediment load transported by a flowing fluid in contact with the mobile stream bed by rolling, sliding and/or saltating, has major

†Email address for correspondence: remi.chassagne@inrae.fr
consequences for public safety, water resources and environmental sustainability. In mountains, steep slopes drive an intense transport of a wide range of grain sizes implying size sorting or segregation (see e.g. Gray 2018), which is largely responsible for our limited ability to predict sediment flux and river morphology (Bathurst 2007; Frey & Church 2011; Dudill et al. 2018). Sediment grain interactions can produce vertical, longitudinal or lateral sorting, which lead to very complex and varied morphologies of bed surface and subsurface (such as armouring, bedload sheet, etc.) (Dietrich et al. 1989; Recking et al. 2009; Frey & Church 2011; Bacchi et al. 2014), and can drastically modify the fluvial geomorphological equilibrium (Dudill, Frey & Church 2017; Dudill et al. 2018, 2020).

The present study focuses on vertical size segregation. For large size ratios, the small particles can percolate spontaneously by gravity, without external forcing, into the bed of larger particles (Bridgwater & Ingram 1971; Dudill et al. 2017). For the smaller size ratios, spontaneous percolation is not possible without deformation of the bed, which can be achieved by shearing or vibrating. This dynamic segregation results from the combination of kinetic sieving (Middleton 1970) and squeeze expulsion (Savage & Lun 1988). When sheared, the granular material dilates and creates gaps, as the layers of particles flow past one another. Kinetic sieving is based on the idea that, under the action of gravity, small particles are more likely to fall down into the created gaps than the large particles, because they are more likely to fit into the available space. This process competes against squeeze expulsion, which tends to push all particles upwards with the same probability. The combination of both processes, which for brevity will just be called gravity driven segregation (Gray 2018), results in a net downward flux of small particles and an upward flux of large grains, leading to an inversely graded bed. Gravity driven segregation has been studied experimentally and numerically in multiple configurations, such as dry granular avalanches (Savage & Lun 1988; Dolgunin & Ukolov 1995; Wiederseiner et al. 2011; Fan et al. 2014; Jones et al. 2018), shear cells (Golick & Daniels 2009; May et al. 2010; Fan & Hill 2015; Van der Vaart et al. 2015; Fry et al. 2018), annular rotating drums (Gray & Ancey 2011) and bedload configurations (Ferdowsi et al. 2017; Lafaye de Micheaux, Ducottet & Frey 2018; Frey et al. 2020).

Particle-size segregation has strong implications for the study of both geomorphology and granular media in general. It has been studied in many configurations but no general description has been proposed yet. In order to identify the mechanisms and comment on the literature on gravity driven segregation, a dimensional analysis is performed. This analysis is made in the dry limit case, consistent with the results of Maurin, Chauchat & Frey (2016) showing that the fluid just acts as a forcing mechanism in turbulent bedload transport and is expected to have only secondary effects on the granular behaviour. Therefore, the fluid is not considered as an influencing parameter for this segregation configuration. The different variables of the problem are

$$d_s, d_l, \rho_p, g, \phi_s, \dot{\gamma}_p, P_p, w_s, \qquad (1.1)$$

where $d_s$ (respectively $d_l$) is the diameter of small (respectively large) particles, $\phi_s$ is the local concentration in small particles, $\rho_p$ is the particle density, $g = 9.81 \text{ ms}^{-2}$ is the gravitational acceleration, $\dot{\gamma}_p$ is the particle shear rate, $P_p$ is the granular pressure and $w_s$ is the segregation velocity of the small particles. These eight variables involve only three units (length, time and mass), so the problem depends on five dimensionless groups. Choosing $d_l, \sqrt{d_l/g}$ and $\rho_p d_l^3$ as reference length, time and
mass units, respectively, the following relation is obtained:

$$ \frac{w_s}{\sqrt{gd_l}} = \mathcal{G} \left( r, \phi_s, \sqrt{\frac{d_l}{g}} \dot{\gamma}_p, \frac{P_p}{\rho_p d_l g} \right), \quad (1.2) $$

where $r = d_l/d_s$ is the size ratio between large and small particles and $\mathcal{G}$ is an unknown function. It is assumed that relation (1.2) can be written in a power law form,

$$ \frac{w_s}{\sqrt{gd_l}} \propto \left( \sqrt{\frac{d_l}{g}} \dot{\gamma}_p \right)^{\zeta} \left( \frac{P_p}{\rho_p d_l g} \right)^{\eta} \mathcal{G}_1(r) \mathcal{G}_2(\phi_s), \quad (1.3) $$

where the function $\mathcal{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 $\dot{\gamma}_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 $\eta = -\zeta/2$, which can be rewritten as

$$ \frac{w_s}{\sqrt{gd_l}} \propto I^\zeta \mathcal{G}_1(r) \mathcal{G}_2(\phi_s), \quad (1.4) $$

with $I$ the large particle inertial number defined as

$$ I = \frac{d_l \dot{\gamma}_p}{\sqrt{P_p / \rho_p}}, \quad (1.5) $$

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 $\phi_s$ has also been widely studied, and different forms for $\mathcal{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(\phi_s) = 1 - \phi_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 $\phi_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:

$$ \frac{\partial \phi_s}{\partial t} + \nabla \cdot (\phi_s \mathbf{u}) - \frac{\partial F_s}{\partial z} = \nabla \cdot (D \nabla \phi_s), \qquad (1.6) $$

where $\mathbf{u}$ is the bulk velocity field, $F_s$ the segregation flux and $D$ the diffusive mixing 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
physics relies on the expression of the segregation flux and of the diffusion flux. The former is defined as $F_s = \phi_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:

$$F_s \propto G_1(r) I^\zeta \phi_s (1 - \phi_s). \qquad (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 $\phi_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).

## 2. 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.
## 2.1. 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:

$$m^p \frac{d^2 \mathbf{x}^p}{dt^2} = \mathbf{f}_g^p + \mathbf{f}_c^p + \mathbf{f}_f^p, \qquad (2.1)$$

$$\mathcal{I}^p \frac{d\boldsymbol{\omega}^p}{dt} = \boldsymbol{\mathcal{T}} = \mathbf{x}_c \times \mathbf{f}_c^p, \qquad (2.2)$$

where $m^p$, $\mathbf{x}^p$, $\boldsymbol{\omega}^p$ and $\mathcal{I}^p$ are respectively the mass, position, angular velocity and moment of inertia of particle $p$. The three major forces are: $\mathbf{f}_g^p$ the gravitational force, $\mathbf{f}_c^p$ the inter-particle contact forces and $\mathbf{f}_f^p$ 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 $\mu_p$ in the tangential direction. The normal and tangential contact forces read accordingly:

$$\left. \begin{aligned} F_n &= -k_n \delta_n - c_n \dot{\delta}_n, \\ F_t &= -\min(k_s \delta_t, \mu_p F_n), \end{aligned} \right\} \qquad (2.3)$$

where $\delta_n$ (respectively $\delta_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 $\mu_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 (Maurin et al. 2015). They are defined as

$$\mathbf{f}_b^p = -\frac{\pi d^{p3}}{6} \nabla P_{x^p}^f, \qquad (2.4)$$

$$\mathbf{f}_D^p = \frac{1}{2} \rho_f \frac{\pi d^{p2}}{4} C_D \|\langle \mathbf{u} \rangle_{x^p}^f - \mathbf{v}^p\| (\langle \mathbf{u} \rangle_{x^p}^f - \mathbf{v}^p), \qquad (2.5)$$

where $d^p$ denotes the diameter of particle $p$, $\langle \mathbf{u} \rangle_{x^p}^f$ is the mean fluid velocity at the position of particle $p$, $P_{x^p}^f$ is the hydrostatic fluid pressure at the position of particle $p$ and $\mathbf{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 - \Phi)^{-3.1}$, with $\Phi$ the total solid phase (small and large particle) volume fraction and $Re_p = \|\langle \mathbf{u} \rangle_{x^p}^f - \mathbf{v}^p\|d^p/v^f$ the particle Reynolds number.

## 2.2. 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:

$$ \rho_f(1 - \Phi) \frac{\partial \langle u_x \rangle^f}{\partial t} = \frac{\partial S_{xz}}{\partial z} + \frac{\partial R_{xz}}{\partial z} + \rho_f(1 - \Phi)g_x - n \langle f_{fx}^p \rangle^s, \quad (2.6) $$

where $\rho_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 \langle f_{fx}^p \rangle^s$ represents the momentum transfer associated with the interaction forces between fluid and particles. The viscous shear stress $S_{xz}$ is taken as

$$ S_{xz} = \rho_f(1 - \Phi)\nu^f \frac{\partial \langle u_x \rangle^f}{\partial z}, \quad (2.7) $$

with $\nu^f$ the pure fluid viscosity. The Reynolds shear stress $R_{xz}$ is based on an eddy viscosity concept,

$$ R_{xz} = \rho_f(1 - \Phi)\nu_t \frac{\partial \langle u_x \rangle^f}{\partial z}. \quad (2.8) $$

The turbulent viscosity $\nu_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),

$$ \nu_t = l_m^2 \left| \frac{\partial \langle u_x \rangle^f}{\partial z} \right|, \quad l_m(z) = \kappa \int_0^z \frac{\Phi_{max} - \Phi(\zeta)}{\Phi_{max}} d\zeta, \quad (2.9a,b) $$

with $\kappa = 0.41$ the von Kármán constant and $\Phi_{max} = 0.61$ the maximal packing of the granular medium (random close packing).

The total momentum $n \langle f_{fx}^p \rangle^s$ transmitted by the fluid to the particles is computed as the sum over each class of the horizontal solid-phase average (denoted $\langle \cdot \rangle^s$) of the momentum transmitted by the drag force on each particle

$$ n \langle f_{fx}^p \rangle^s = \frac{\Phi_s}{\pi d_s^3 / 6} \langle f_{fx}^{ps} \rangle^s + \frac{\Phi_l}{\pi d_l^3 / 6} \langle f_{fx}^{pl} \rangle^s, \quad (2.10) $$

and introducing the expression of the drag force (2.5),

$$ n \langle f_{fx}^p \rangle^s = \frac{3}{4} \frac{\Phi_s \rho_f}{d_s} \langle C_D \| \langle \mathbf{u} \rangle_{x^{ps}}^f - \mathbf{v}^{ps} \| (\langle \mathbf{u} \rangle_{x^{ps}}^f - \mathbf{v}^{ps}) \rangle^s + \frac{3}{4} \frac{\Phi_l \rho_f}{d_l} \langle C_D \| \langle \mathbf{u} \rangle_{x^{pl}}^f - \mathbf{v}^{pl} \| (\langle \mathbf{u} \rangle_{x^{pl}}^f - \mathbf{v}^{pl}) \rangle^s, \quad (2.11) $$

where $p_s$ (respectively $p_l$) denotes the ensemble of the small (respectively large) particles, $\Phi_s$ (respectively $\Phi_l$) the volume fraction of small (respectively large) particles and $d_s$ (respectively $d_l$) the diameter of small (respectively large) particles.
The solid-phase volume average $\langle \cdot \rangle^s$ is defined following Jackson (2000) for which a cuboid weighting function $\mathcal{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

$$ \Phi_s(z) = \sum_{p_s} \int_{V_{p_s}} \mathcal{F}(|z - z'|) dz', \qquad (2.12) $$

$$ \Phi_l(z) = \sum_{p_l} \int_{V_{p_l}} \mathcal{F}(|z - z'|) dz', \qquad (2.13) $$

$$ \langle \gamma \rangle^s = \frac{1}{\Phi_s(z)} \sum_{p_s} \int_{V_{p_s}} \gamma(z') \mathcal{F}(|z - z'|) dz' + \frac{1}{\Phi_l(z)} \sum_{p_l} \int_{V_{p_l}} \gamma(z') \mathcal{F}(|z - z'|) dz', \qquad (2.14) $$

where $\gamma$ is a scalar quantity.

The fluid model, which is classical in sediment transport (e.g. Drake & Calantoni 2001; Hsu & Liu 2004; Durán, Andreotti & Claudin 2012; Revil-Baudard & Chauchat 2013; Maurin et al. 2015; Chauchat 2018) is only closed using a mixing length model and a closure for the drag force formulation. The latter are usual in the literature, and it has been shown in Maurin (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.

## 3. 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.

### 3.1. 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 $\rho^p = 2500$ kg m$^{-3}$ and $\rho^f = 1000$ kg m$^{-3}$. The size of the 3-D domain is $30d_l \times 30d_l$ in the horizontal plane in order to have converged average values (Maurin et al. 2015) 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 \times 30d_l \times 30d_l \times N_l d_l$ (respectively $0.61 \times 30d_l \times 30d_l \times 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
<figure><img src="image_3.png" /><figcaption>FIGURE 1. A typical numerical set-up for $N_s = 2$. Initially, $N_s$ layers of small particles ($d_s = 3$ mm) are deposited by gravity above $N_l$ layers of large particles ($d_l = 6$ mm). $N_s$ varies from 0.01 until 2 with $N_l$ varying accordingly in order to always have $N_l d_l + N_s d_s = 10d_l$.</figcaption></figure>

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 $\theta = \tau_f / [(\rho_p - \rho_f)gd_l]$ and simulations are performed at $\theta \simeq 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 $\phi_s$ (respectively $\phi_l$). By definition, the two volume fractions sum to unity,

$$ \phi_s + \phi_l = \frac{\Phi_s}{\Phi_s + \Phi_l} + \frac{\Phi_l}{\Phi_s + \Phi_l} = 1, \qquad (3.1) $$

where $\Phi_s$ (respectively $\Phi_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:

$$ \tilde{z} = \frac{z}{d_l}, \quad \tilde{t} = t \sqrt{g/d_l}, \quad \langle \tilde{v}_x \rangle^p = \frac{\langle v_x \rangle^p}{\sqrt{gd_l}}. \qquad (3.2a-c) $$

In the following, the tildes are dropped for sake of clarity.

## 3.2. 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
<figure><img src="image_4.png" /><figcaption>FIGURE 2. (a) Typical temporal evolution of fine particle concentration profile $\phi_s$ for the case $N_s = 2$, (b) concentration profiles in small particles for different times for the case $N_s = 2$ and (c) time evolution of the maximal value of $\phi_s$ for all simulations.</figcaption></figure>

layer made of large particles. At the beginning, the small particles infiltrate rapidly into the first few layers of large particles. As small particles infiltrate downward, large particles rise to the surface. The DEM simulations exhibit a two-layer structure, with small particles sandwiched between two layers of large particles. While infiltrating,
<figure><img src="image_5.png" /><figcaption>FIGURE 3. (a) Streamwise space-time-averaged particle velocity as a function of the height, (b) time-averaged volume fraction of the granular mixture and (c) evolution of the vertical position of the mass centre of small particles with time.</figcaption></figure>

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 ($\Phi \sim 0.61$). Above $z \sim 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

$$z_c = -a \ln(t) + b. \qquad (3.3)$$

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 $\phi_s$.

### 3.3. Dependence on the inertial number

In figure 4(a), the dimensionless segregation velocity is plotted against the large particle inertial number $I = \dot{\gamma}_p d_l / \sqrt{P_p / \rho_p}$, where $\dot{\gamma}_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,

$$\frac{dz_c(t)}{dt} \propto I^{0.85}(z_c(t)). \qquad (3.4)$$

Interestingly, Fry et al. (2018) obtained a very similar result for dry granular flows at higher inertial numbers ($I \in [10^{-2}, 1]$). They obtained an exponent 0.84, very close to the 0.85 exponent obtained in the present configuration. Our results suggest that the behaviour observed by Fry et al. (2018) is valid in a wider range of inertial numbers, in particular, in the quasi-static regime ($I \in [10^{-5}, 10^{-1}$]).

The scaling with the inertial number (3.4) is in line with the theory of 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 \le 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
<figure><img src="image_6.png" /><figcaption>FIGURE 4. (a) Segregation velocity dependence on the local inertial number, (b) inertial number profile.</figcaption></figure>

exponential profile to the power 0.85, $I^{0.85}$ is also an exponential function of $z$ and can be written as

$$I^{0.85}(z) = I_0 e^{z/c}. \qquad (3.5)$$

Introducing (3.5) into expression (3.4) leads to

$$\frac{dz_c(t)}{dt} = -b_1 e^{z_c(t)/c}, \qquad (3.6)$$

with $b_1$ a positive constant. Equation (3.6) can be integrated

$$\int e^{-z_c/c} dz_c = \int -b_1 dt, \qquad (3.7)$$

leading to

$$ce^{-z_c/c} = b_1 t + \text{Const} \sim b_1 t, \qquad (3.8)$$

for long times. It can be rewritten as

$$z_c(t) = -c \ln(t) + b_2, \qquad (3.9)$$

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$.
<table><caption>TABLE 1. Values of coefficients *a* (slope of the centre of mass) and *c* (exponential decay of the inertial number to the power 0.85) obtained from fitting the curves of figures 3(c) and 4(b) and the error in percentage.</caption><thead><tr><th>$N_s$</th><th>a</th><th>c</th><th>Error (%)</th></tr></thead><tbody><tr><td>0.01</td><td>0.863</td><td>0.873</td><td>1.159</td></tr><tr><td>0.05</td><td>0.850</td><td>0.864</td><td>1.647</td></tr><tr><td>0.25</td><td>0.823</td><td>0.847</td><td>2.916</td></tr><tr><td>0.5</td><td>0.843</td><td>0.862</td><td>2.254</td></tr><tr><td>1</td><td>0.833</td><td>0.857</td><td>2.800</td></tr><tr><td>1.5</td><td>0.820</td><td>0.828</td><td>0.976</td></tr><tr><td>2</td><td>0.849</td><td>0.828</td><td>2.473</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>FIGURE 5. (a) Centre of mass velocity dependence on the inertial number at the bottom of the layer of small particles, (b) small particle concentration profiles at time $t = 40435$ and vertical position of the bottom of the layer (×).</figcaption></figure>

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 $I_b(z_c) = I(z_b)$,

$$ \frac{dz_c(t)}{dt} = -\alpha_0 I_b^{0.85}(z_c(t)), \qquad (3.10) $$

where $\alpha_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
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 = 40435$, 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.

### 3.5. 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:

$$ \frac{dz_c(t)}{dt} = -\alpha_0 \mathcal{G}_1(r) I_b^{0.85}(z_c(t)). \quad (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 $\mathcal{G}_1$ is proposed:

$$ \mathcal{G}_1(r) = 0.45(e^{(r-1)/1.59} - 1), \quad (3.12) $$

empirically fitted from the data points of figure 6(b), and plotted in dashed line in figure 6(b).
<figure><img src="image_8.png" /><figcaption>FIGURE 6. (a) Centre of mass velocity dependence on the inertial number at the bottom of the layer of small particles for different size ratio. Note that for $r = 2.5$, $r = 2.75$ and $r = 3$, simulations were performed with a $15d_l \times 15d_l$ box in order to gain computation time. (b) Best fit coefficient between the segregation velocity and the inertial number to the power 0.85 times $\alpha_0$ as a function of the size ratio minus one (for each curve of figure 6a, one point is obtained). The function $G_1$ is plotted as dashed line.</figcaption></figure>

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:

$$ \frac{\mathrm{d}z_c}{\mathrm{d}t} = -\alpha_0 G_1(r) I_b^{0.85}, \qquad (3.13) $$

where $\alpha_0$ can be computed from figure 5(a) and $G_1(r)$ is given by (3.12).

### 3.6. Dependence on the local concentration $\phi_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).
## 4. Discussion in the framework of continuum modelling

### 4.1. 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,

$$ \frac{\partial \phi_s}{\partial t} - \frac{\partial F_s}{\partial z} = 0. \qquad (4.1) $$

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

$$ F_s = \alpha_0 \mathcal{G}_1(r) I_0^{0.85} \phi_s(1 - \phi_s) = S_{r0} e^{z/c} \phi_s(1 - \phi_s), \qquad (4.2) $$

with $S_{r0} = \alpha_0 \mathcal{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 \times 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,

$$ \phi_s(z, t=0) = \begin{cases} 0, & z < z_i, \\ 1, & z \ge z_i, \end{cases} \qquad (4.3) $$

with $z_i = N_i d_i$. 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 \le t \le t_1$ is

$$ \phi_s(z, t) = \begin{cases} 0, & z < z_1, \\ \phi_{fan}(z, t), & z_1 \le z \le z_2, \\ 0, & z > z_2, \end{cases} \qquad (4.4) $$

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 $\phi_{fan}$ is the solution in the expansion fan given by (4.7):

$$ z_1(t) = z_i - c \ln \left( 1 + \frac{S_{r0}}{c} e^{z_i/c} t \right), \qquad (4.5) $$
<figure><img src="image_9.png" /><figcaption>FIGURE 7. (a) Case $N_s = 1$, $r = 1.5$ space-time evolution of the concentration in fine particles $\phi_s$, and (b) profiles of concentration at time $t = 80870$ with the continuum model without diffusion for several values of $N_s$. Long time expression (4.9) is indicated as dashed line. The same values of $S_{r0}$ and $c$ have been used.</figcaption></figure>

$$ \frac{\mathrm{d}z_2}{\mathrm{d}t} = -S_{r0}e^{z_2/c} (1 - \phi_s(z_2, t)), \qquad (4.6) $$

$$ \phi_{\text{fan}}(z, t) = \frac{1}{2} - \frac{c}{S_{r0}t}e^{-z/c} + \frac{c}{S_{r0}t}e^{-(z+z_i)/(2c)} \sqrt{1 + \left(\frac{S_{r0}t}{2c}\right)^2 e^{(z+z_i)/c}}. \quad (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 \gg 4c/S_{r0}e^{z_i/c} \sim 1000$, it can be shown (see appendix A) that the square root term in (4.7) simplifies to

$$ 1 + \left(\frac{S_{r0}t}{2c}\right)^2 e^{(z+z_i)/c} \sim \left(\frac{S_{r0}t}{2c}\right)^2 e^{(z+z_i)/c}, \qquad (4.8) $$

and the volume fraction profile $\phi_s$ in the rarefaction fan simplifies as

$$ \phi_{\text{fan}}(z, t) = 1 - \frac{c}{S_{r0}t}e^{-z/c}. \qquad (4.9) $$

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

$$ \frac{\mathrm{d}z_2}{\mathrm{d}t} = -\frac{c}{t}, \qquad (4.10) $$

for which the solution is

$$ z_2(t) = -c \ln(t) + \beta, \qquad (4.11) $$

where $\beta$ is an integration constant. Concerning the characteristic curve $z_1$, the term in the logarithm in (4.5) simplifies for long times to

$$ 1 + \frac{S_{r0}}{c} e^{z_i/c} t \sim \frac{S_{r0}}{c} e^{z_i/c} t, \qquad (4.12) $$
and the following expression is obtained for $z_1$:

$$z_1(t) = -c \ln(t) - c \ln\left(\frac{S_{r0}}{c}\right). \qquad (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 \propto I^{0.85}\phi_s(1 - \phi_s)$ implies that the vertical velocity of the small particles $w_s = -F_s/\phi_s$ is

$$w_s(z, t) = -\frac{S_{r0}}{I_0} I^{0.85}(z)(1 - \phi_s(z, t)). \qquad (4.14)$$

Using (4.9), the vertical velocity can be expressed for long times as

$$w_s(z) = -\frac{S_{r0}}{I_0} I^{0.85} \frac{c}{S_{r0}t} e^{-z/c}, \qquad (4.15)$$

and remembering that $I^{0.85} = I_0 e^{z/c}$,

$$w_s(z) = -\frac{c}{t}. \qquad (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,

$$w_s(z, t) = -\frac{S_{r0}}{I_0} I^{0.85}(z)(1 - \phi_s(z, t)) = -c/t. \qquad (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 - \phi_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 $\phi_s = 0$ is considered. It corresponds, in figure 7(b), to the position where all profiles collapse to zero at the bottom. At this position, the segregation velocity is only a function of the inertial number at the bottom of the layer

$$w_s(z, t) = -\frac{S_{r0}}{I_0} I^{0.85}(z_1) = -c/t. \qquad (4.18)$$
<table><caption>TABLE 2. Mean error between the bottom position from the DEM simulations $z_b$ and from the continuum model $z_1$. The maximal error of less than 1% between both positions suggests that they correspond.</caption><tr><td>N<sub>s</sub></td><td>0.01</td><td>0.05</td><td>0.25</td><td>0.5</td><td>1</td><td>1.5</td><td>2</td></tr><tr><td>|z<sub>1</sub> - z<sub>b</sub>| / z<sub>1</sub> (%)</td><td>0.28</td><td>0.22</td><td>0.50</td><td>0.09</td><td>0.04</td><td>0.73</td><td>0.95</td></tr></table>

<figure><img src="image_10.png" /><figcaption>FIGURE 8. Comparison between DEM simulations (---) and segregation model (—) of the profile of concentration $\phi_s$ for the case $N_s = 1$ at time $t = 80870$.</figcaption></figure>

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.

## 4.2. *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
results. The equation to be solved is given by

$$ \frac{\partial \phi_s}{\partial t} - \frac{\partial F_s}{\partial z} = \frac{\partial}{\partial z} \left( D \frac{\partial \phi_s}{\partial z} \right), \qquad (4.19) $$

with a segregation flux of the form $F_s = S_{r0}e^{z/c}\phi_s(1 - \phi_s)$, with $S_{r0} = 3.70 \times 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 $\phi_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:

$$ \tau = t, \quad \xi = z + c \ln(t). \qquad (4.20a,b) $$

Equation (4.19) in the moving frame becomes

$$ \frac{\partial \phi_s}{\partial \tau} + \frac{c}{\tau} \frac{\partial \phi_s}{\partial \xi} - \frac{\partial}{\partial \xi} \left( \frac{S_{r0}}{\tau} e^{\xi/c} \phi_s (1 - \phi_s) \right) = \frac{\partial}{\partial \xi} \left( D(\xi - c \ln(\tau)) \frac{\partial \phi_s}{\partial \xi} \right). \qquad (4.21) $$

For a travelling wave, the solution to (4.21) should not depend on time $\tau$. Assuming $\partial \phi_s / \partial \tau = 0$, this requires that $\tau D(\xi - c \ln(\tau))$ is independent of $\tau$. This implies

$$ \frac{\partial}{\partial \tau} (\tau D(\xi - c \ln(\tau))) = 0, \qquad (4.22) $$

and distributing the derivative with $\tau$, the following differential equation is obtained:

$$ D'(\xi - c \ln(\tau)) - \frac{1}{c} D(\xi - c \ln(\tau)) = 0, \qquad (4.23) $$

for which the solution is

$$ D = D_0 e^{(\xi - c \ln(\tau))/c} = D_0 e^{z/c}. \qquad (4.24) $$

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

$$ c \frac{\partial \phi_s}{\partial \xi} - \frac{\partial}{\partial \xi} \left( S_{r0} e^{\xi/c} \phi_s (1 - \phi_s) \right) = \frac{\partial}{\partial \xi} \left( D_0 e^{\xi/c} \frac{\partial \phi_s}{\partial \xi} \right). \qquad (4.25) $$

By integration, and imposing $\phi_s \to 0$ when $\xi \to \pm\infty$, the following equation is obtained:

$$ \frac{\partial \phi_s}{\partial \xi} - \frac{S_{r0}}{D_0} \phi_s \left[ \phi_s - \left( 1 - \frac{c}{S_{r0}} e^{-\xi/c} \right) \right] = 0, \qquad (4.26) $$
<figure><img src="image_11.png" /><figcaption>FIGURE 9. (a) Travelling wave solution (—) and concentration profiles in the moving frame from DEM simulations at different times for the case $N_s = 1$, $r = 1.5$, $Pe = 3.86$. $t = 40435$ (----), $t = 50543$ (— · —, orange), $t = 60652$ (· · · · ·, green), $t = 70761$ s (-------, red), $t = 80870$ (----·, violet). (b) Value of the Péclet number as a function of the number of layers of small particles.</figcaption></figure>

for which the solution is

$$ \phi_s = \frac{e^{-(S_{r0}/D_0)\xi - (c^2/D_0)e^{-\xi/c}}}{C - \frac{S_{r0}}{D_0} \int_{-\infty}^{\xi} (e^{-(S_{r0}/D_0)\xi' - (c^2/D_0)e^{-\xi'/c}}) d\xi'} \quad (4.27) $$

where $C$ is a constant satisfying the mass conservation $\int_{-\infty}^{+\infty} \phi_s d\xi$. The travelling wave solution has a complex dependency on $\xi$ 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 $\int_{-\infty}^{+\infty} \phi_s d\xi$ 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 $\phi_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 $\phi_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.

## 5. 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 grain-grain forces involved in size segregation as well as collective effects related to particle concentration.

## Acknowledgements

IRSTEA is member of Labex TEC21 (Investissements d'Avenir Grant Agreement ANR-11-LABX-0030) and Labex Osug@2020 (Investissements d'Avenir Grant Agreement ANR-10-LABX-0056). This research was funded by the French 'Agence nationale de la recherche', project ANR-16-CE01-0005 SegSed 'size segregation in sediment transport'. The authors acknowledge the support of Irstea (now INRAE, formerly Cemagref). This research was supported by NERC grants NE/E003206/1 and NE/K003011/1as well as EPSRC grants EP/I019189/1, EP/K00428X/1 and EP/M022447/1. J.M.N.T.G. is a Royal Society Wolfson ResearchMerit Award holder (WM150058) and an EPSRC Established Career Fellow (EP/M022447/1).
Declaration of interests

The authors report no conflict of interest.

Appendix A. Solution of the purely advective segregation model

May et al. (2010) already derived a solution to the model of Thornton et al. (2006) with an exponential dependence of the segregation flux with z. In this appendix, the full solution is presented in the bedload configuration. The segregation equation to solve is

$$ \frac{\partial \phi_s}{\partial t} - \frac{\partial F_s}{\partial z} = 0, \qquad (A 1) $$

with the flux $F_s = S_{r0}e^{z/c}\phi_s(1 - \phi_s)$ and an initial step condition in concentration

$$ \phi_{s0} = \phi_s(z, 0) = \begin{cases} 0, & z < z_i, \\ 1, & z \ge z_i. \end{cases} \qquad (A 2) $$

The boundary conditions ensure that there is no flux at $z = 0$, $H$, which requires that

$$ \phi_s = 0 \text{ or } 1, \quad \text{at } z = 0 \text{ and } H. \qquad (A 3) $$

Substituting the flux into (A 1) implies

$$ \frac{\partial \phi_s}{\partial t} + W(z, \phi_s) \frac{\partial \phi_s}{\partial z} = S(z, \phi_s), \qquad (A 4) $$

where

$$ W(z, \phi_s) = S_{r0}e^{z/c}(2\phi_s - 1), \qquad (A 5) $$

$$ S(z, \phi_s) = \frac{S_{r0}}{c} e^{z/c} \phi_s (1 - \phi_s). \qquad (A 6) $$

Due to dependence of the flux both on $z$ and $\phi_s$, the equation to solve has a nonlinear advective velocity and a source term. In the following, this problem is solved using the method of characteristics.

On a characteristic curve $(Z(s), t(s))$,

$$ \frac{d\phi_s}{ds} = \frac{\partial\phi_s}{\partial t}\frac{dt}{ds} + \frac{\partial\phi_s}{\partial z}\frac{dz}{ds} = S. \qquad (A 7) $$

Comparing coefficients of (A 4) and (A 7)

$$ \frac{d\phi_s}{ds} = S, \qquad (A 8) $$

$$ \frac{dt}{ds} = 1, \qquad (A 9) $$

$$ \frac{dZ}{ds} = W. \qquad (A 10) $$

Eliminating $s$ implies

$$ \frac{d\phi_s}{dt} = S \qquad (A 11) $$
$$ \frac{dZ}{dt} = W. \qquad (A 12) $$

The problem (A 4) is thus equivalent to the following system:

$$ \left. \begin{aligned} \frac{d}{dt} \phi_s(Z(t), t) &= S(Z(t), \phi_s(Z(t), t)), \\ \frac{dZ(t)}{dt} &= W(Z(t), \phi_s(Z(t), t)), \\ Z(0) &= z_0. \end{aligned} \right\} \qquad (A 13) $$

The first equation indicates that $\phi_s$ is not constant along a characteristic curve and the second equation that the characteristic curves are not linear. The method of characteristics consists in solving the value of $\phi_s$ along the characteristics curves. If a characteristic curve passing by each time and space position exists, then the problem can be fully resolved. In order to compute the value of $\phi_s$ along a characteristic curve, the first equation of (A 13) is differentiated with time

$$ \frac{d^2}{dt^2} \phi_s(Z(t), t) = \frac{d}{dt} [S(Z(t), \phi_s(Z(t), t))], \qquad (A 14) $$

distributing the derivative

$$ \frac{d^2}{dt^2} \phi_s(Z(t), t) = \frac{dZ}{dt} \frac{\partial S}{\partial z} + \frac{d\phi_s}{dt} \frac{Z(t), t}{\partial \phi_s}, \qquad (A 15) $$

and introducing the first and second equations of (A 13) it implies

$$ \frac{d^2}{dt^2} \phi_s(Z(t), t) = W \frac{\partial S}{\partial z} + S \frac{\partial S}{\partial \phi_s}. \qquad (A 16) $$

Using (A 5) and (A 6), one can remark that $\partial S/\partial z = S/c$ and $\partial S/\partial \phi_s = -W/c$. Introducing this in (A 16), the second derivation with time of $\phi_s$ along a characteristic curve is identically null,

$$ \frac{d^2}{dt^2} \phi_s(Z(t), t) = 0. \qquad (A 17) $$

Integrating twice with time, a linear dependence on time of $\phi_s$ is obtained:

$$ \phi_s(Z(t), t) = \phi_{s0}(z_0) + S(z_0, \phi_{s0}(z_0))t. \qquad (A 18) $$

The value of $\phi_s$ along the characteristic curve is entirely determined by its initial value and therefore by the initial condition $\phi_{s0}$. It is now possible to compute the position with time of a characteristic curve. Introducing (A 18) in the second equation of (A 13) and integrating, a relation for the characteristic curves is obtained:

$$ e^{-Z(t)/c} - e^{-z_0/c} = -\frac{S_{r0}}{c}(2\phi_{s0}(z_0) - 1)t - \frac{S_{r0}}{c}S(z_0, \phi_{s0}(z_0))t^2, \qquad (A 19) $$

$$ \Leftrightarrow Z(t) = z_0 - c \ln \left( 1 - \frac{1}{c} W(z_0, \phi_{s0}(z_0)) t - \frac{S_{r0}}{c} e^{z_0/c} S(z_0, \phi_{s0}(z_0)) t^2 \right), \qquad (A 20) $$
<figure><img src="image_12.png" /><figcaption>FIGURE 10. Example of solution obtained by the method of characteristics: $z_i = 6$ (a) first time steps; (b) full solution.</figcaption></figure>

and the characteristic curve is entirely determined by the value of the initial solution $\phi_{s0}$ in its origin position $z_0$. Finally the solution to system (A 13) is

$$ \left. \begin{array}{l} \phi_s(Z(t), t) = \phi_{s0}(z_0) + S(z_0, \phi_{s0}(z_0)) t, \\ Z(t) = z_0 - c \ln \left( 1 - \frac{1}{c} W(z_0, \phi_{s0}(z_0)) t - \frac{S_{r0}}{c} e^{z_0/c} S(z_0, \phi_{s0}(z_0)) t^2 \right). \end{array} \right\} \quad (A 21) $$

Equation (A 21) is a general solution for any initial condition. In the present configuration the initial solution (A 2) is a step of concentration with a discontinuity at $z_i$ and represents an initial state where a pure layer of small particles is placed above a bed of only large particles. Figure 10 shows the full solution for the case $z_i = 6$, and the derivation of this solution is presented below. According to the second equation of (A 13), the characteristic curves coming from $z_0 < z_i$ have a negative velocity and are going down. Conversely, those coming from $z_0 > z_i$ are going up. Therefore, there is a rarefaction fan, where no characteristic curves are present.

For $z_0 < z_i$, $\phi_{s0}(z_0) = 0$, the equation of the characteristic curves is $Z(t) = z_0 - c \ln(1 + (S_{r0}/c)e^{z_0/c}t)$. Along these characteristics, $\phi_s(Z(t), t) = 0$. Here, $z_1$ is defined as the characteristic curve at the bottom of the rarefaction zone (coming from $z_i$): $z_1(t) = z_i - c \ln(1 + (S_{r0}/c)e^{z_i/c}t)$. Here, $z_1$ is the lower bound of the rarefaction zone. The time at which $z_1$ meets the bottom of the domain is called $t_1 = (c/S_{r0})(1 - e^{-z_i/c})$ and corresponds to the time when the first small particle reaches the bottom.

For $z_0 > z_i$, $\phi_{s0}(z_0) = 1$, the equation of the characteristic curves is $Z(t) = z_0 - c \ln(1 - (S_{r0}/c)e^{z_0/c}t)$. Along these characteristics, $\phi_s(Z(t), t) = 1$. $z_2$ is defined as the characteristic curve coming from $z_i$: $z_2(t) = z_i - c \ln(1 - (S_{r0}/c)e^{z_i/c}t)$. $z_2$ is the upper bound of the rarefaction zone. The time at which $z_2$ meets the top of the domain is called $t_2 = (c/S_{r0})(e^{-z_i/c} - e^{-H/c})$ and corresponds to the time when the first large particle reaches the top.

Once $z_2$ reaches the top of domain at $t_2$ the boundary condition (A 3) implies that $\phi(H, t_2)$ jumps from 1 to 0. A shock forms and the position of the shock must satisfy the Rankine–Hugoniot relation, which ensures the conservation of momentum through the shock

$$ z_2' = - \frac{F_s(\phi^+) - F_s(\phi^-)}{\phi^+ - \phi^-}. \quad (A 22) $$
Above the shock, $\phi^+ = 0$, and below the shock, $\phi^- = 0$ is the solution in the rarefaction zone

$$z_2' = -S_{r0}e^{z_2/c} (1 - \phi_s(z_2, t)). \quad (A 23)$$

A similar expression applying the Rankine-Hugoniot relation is found for $z_1$ for $t > t_1$. For now, as no solution in the rarefaction has been computed, the shock evolution equation (A 23) cannot be solved. The solution in the rarefaction zone is now computed. The rarefaction zone is delimited by $z_1$ and $z_2$. Due to the discontinuity in the initial solution, no characteristic curves are present in the rarefaction wave ($z_1 \le z \le z_2$). In order to find a solution, it is considered that there is an infinity of characteristic curves coming from the point of discontinuity $z_i$, each of them being associated with a different initial value, noted $\bar{\phi}_{s0}$ and verifying $\phi_{s0}(z_c)^- = 0 \le \bar{\phi}_{s0} \le 1 = \phi_{s0}(z_c)^+$. In other words, for $(z, t)$ in the rarefaction zone, one looks for a characteristic curve coming from $z_i$ associated with an initial value $\bar{\phi}_{s0}$ and verifying system (A 21)

$$\left.
\begin{array}{l}
\phi_s(z, t) = \bar{\phi}_{s0} + S(z_i, \bar{\phi}_{s0}) t, \\
e^{-z/c} - e^{-z_i/c} = -\frac{S_{r0}}{c}(2\bar{\phi}_{s0} - 1)t - \frac{S_{r0}}{c}e^{z_i/c}S(z_i, \bar{\phi}_{s0})t^2.
\end{array}
\right\}
\quad (A 24)$$

To compute the value of $\phi_s$ in the rarefaction zone with the first equation of (A 24), only the value of $\bar{\phi}_{s0}$ is missing. It is computed by an inversion of the second equation of (A 24) (quadratic equation and choosing the relevant solution)

$$\bar{\phi}_{s0} = \frac{1}{2} + \frac{c}{S_{r0}t} e^{-z_i/c} \left( 1 - e^{(z_i-z)/c} \sqrt{1 + \left(\frac{S_{r0}t}{2c}\right)^2 e^{(z+z_i)/c}} \right). \quad (A 25)$$

Inserting (A 25) into the first equation of (A 24) the solution in the rarefaction zone is

$$\phi_s(z, t) = \frac{1}{2} - \frac{c}{S_{r0}t} e^{-z/c} + \frac{c}{S_{r0}t} e^{-(z+z_i)/(2c)} \sqrt{1 + \left(\frac{S_{r0}t}{2c}\right)^2 e^{(z+z_i)/c}}. \quad (A 26)$$

The equation of the shock (A 23) can now be solved. However, due to the complex form of the solution in the rarefaction zone (A 26) the equation of the shock (A 23) is solved numerically. Finally the full solution is presented in figure 10.

A long time asymptotic analysis can be performed on the solution in the rarefaction zone (A 26). Remembering that the lower bound of the rarefaction zone is $z_1$ and focusing on the square root term in (A 26)

$$\left(\frac{S_{r0}t}{2c}\right)^2 e^{(z+z_i)/c} \ge \left(\frac{S_{r0}t}{2c}\right)^2 e^{(z_1+z_i)/c}. \quad (A 27)$$

The expression of $z_1$ is

$$z_1(t) = z_i - c \ln \left( 1 + \frac{S_{r0}}{c} e^{z_i/c} t \right). \quad (A 28)$$

For $t \gg t_0 = c/S_{r0}e^{-z_i/c}$,

$$1 + \frac{S_{r0}}{c} e^{z_i/c} t \sim \frac{S_{r0}}{c} e^{z_i/c} t, \quad (A 29)$$
and the following expression is obtained for $z_1$:

$$z_1(t) = -c \ln \left( \frac{S_{r0}t}{c} \right). \qquad (A 30)$$

Replacing the latter expression of $z_1$ in (A 27),

$$\left(\frac{S_{r0}t}{2c}\right)^2 e^{(z+z_i)/c} \ge \frac{S_{r0}t}{4c} e^{z_i/c}, \qquad (A 31)$$

and for $t \gg 4t_0 = 4c/S_{r0}e^{-z_i/c}$,

$$\left(\frac{S_{r0}t}{2c}\right)^2 e^{(z+z_i)/c} \ge \frac{S_{r0}t}{4c} e^{z_i/c} \gg 1. \qquad (A 32)$$

Finally for $t \gg 4t_0$,

$$1 + \left(\frac{S_{r0}t}{2c}\right)^2 e^{(z+z_i)/c} \sim \left(\frac{S_{r0}t}{2c}\right)^2 e^{(z+z_i)/c}, \qquad (A 33)$$

and the solution in the expansion fan simplifies to

$$\phi_s(z, t) = 1 - \frac{c}{S_{r0}t} e^{-z/c}. \qquad (A 34)$$

In the present configuration, $c \sim 1$, $S_{r0} \sim 1e-8$, $z_i \sim 10$. The simplified expression of $\phi_s$ is therefore valid for dimensionless times larger than $t_0 \sim 1e3$, time achieved very rapidly during simulations.

REFERENCES

BACCHI, V., RECKING, A., ECKERT, N., FREY, P., PITON, G. & NAAIM, M. 2014 The effects of kinetic sorting on sediment mobility on steep slopes. *Earth Surf. Process. Landf.* **39** (8), 1075–1086.

BATHURST, J. C. 2007 Effect of coarse surface layer on bed-load transport. *J. Hydraul. Engng* **133** (11), 1192–1205.

BRIDGWATER, J., FOO, W. S. & STEPHENS, D. J. 1985 Particle mixing and segregation in failure zones – theory and experiment. *Powder Technol.* **41** (2), 147–158.

BRIDGWATER, J. & INGRAM, N. D. 1971 Rate of spontaneous inter-particle percolation. *Trans. Inst. Chem. Engrs* **49** (3), 163–169.

CAMPBELL, C. S. 1997 Self-diffusion in granular shear flows. *J. Fluid Mech.* **348**, 85–101.

CHAUCHAT, J. 2018 A comprehensive two-phase flow model for unidirectional sheet-flows. *J. Hydraul. Res.* **56** (1), 15–28.

DIETRICH, W. E., KIRCHNER, J. W., IKEDA, H. & ISEYA, F. 1989 Sediment supply and the development of the coarse surface layer in gravel-bedded rivers. *Nature* **340** (6230), 215–217.

DOLGUNIN, V. N. & UKOLOV, A. A. 1995 Segregation modeling of particle rapid gravity flow. *Powder Technol.* **83** (2), 95–103.

DRAKE, T. G. & CALANTONI, J. 2001 Discrete particle model for sheet flow sediment transport in the nearshore. *J. Geophys. Res.* **106** (C9), 19859–19868.

DUDILL, A., FREY, P. & CHURCH, M. 2017 Infiltration of fine sediment into a coarse mobile bed: a phenomenological study. *Earth Surf. Process. Landf.* **42** (8), 1171–1185.
DUDILL, A., LAFAYE DE MICHEAUX, H., FREY, P. & CHURCH, M. 2018 Introducing finer grains into bedload: the transition to a new equilibrium. *J. Geophys. Res.* **123** (10), 2602–2619.

DUDILL, A., VENDITTI, J., CHURCH, M. & FREY, P. 2020 Comparing the behaviour of spherical beads and natural grains in bedload mixtures. *Earth Surf. Process. Landf.* **45** (4), 831–840.

DURÁN, O., ANDREOTTI, B. & CLAUDIN, P. 2012 Numerical simulation of turbulent sediment transport, from bed load to saltation. *Phys. Fluids* **24** (10), 103306.

FAN, Y. & HILL, K. M. 2015 Shear-induced segregation of particles by material density. *Phys. Rev. E* **92** (2), 022211.

FAN, Y., SCHLICK, C. P., UMBANHOWAR, P. B., OTTINO, J. M. & LUEPTOW, R. M. 2014 Modelling size segregation of granular materials: The roles of segregation, advection and diffusion. *J. Fluid Mech.* **741**, 252–279.

FERDOWSI, B., ORTIZ, C. P., HOUSSAIS, M. & JEROLMACK, D. J. 2017 River-bed armouring as a granular segregation phenomenon. *Nat. Commun.* **8**, 1363.

FREY, P. 2014 Particle velocity and concentration profiles in bedload experiments on a steep slope. *Earth Surf. Process. Landf.* **39** (5), 646–655.

FREY, P. & CHURCH, M. 2011 Bedload: A granular phenomenon. *Earth Surf. Process. Landf.* **36** (1), 58–69.

FREY, P., LAFAYE DE MICHEAUX, H., BEL, C., MAURIN, R., RORSMAN, K., MARTIN, T. & DUCOTTET, C. 2020 Experiments on grain size segregation in bedload transport on a steep slope. *Adv. Water Resour.* **136**, 103478.

FRY, A. M., UMBANHOWAR, P. B., OTTINO, J. M. & LUEPTOW, R. M. 2018 Effect of pressure on segregation in granular shear flows. *Phys. Rev. E* **97** (6), 062906.

GAJJAR, P. & GRAY, J. M. N. T. 2014 Asymmetric flux models for particle-size segregation in granular avalanches. *J. Fluid Mech.* **757**, 297–329.

GDR MİDI 2004 On dense granular flows. *Eur. Phys. J. E* **14** (4), 341–365.

GOLICK, L. A. & DANIELS, K. E. 2009 Mixing and segregation rates in sheared granular materials. *Phys. Rev. E* **80** (4), 042301.

GRAY, J. M. N. T. 2018 Particle segregation in dense granular flows. *Annu. Rev. Fluid Mech.* **50** (1), 407–433.

GRAY, J. M. N. T. & ANCEY, C. 2011 Multi-component particle-size segregation in shallow granular avalanches. *J. Fluid Mech.* **678**, 535–588.

GRAY, J. M. N. T. & CHUGUNOV, V. A. 2006 Particle-size segregation and diffusive mixing in shallow granular avalanches. *J. Fluid Mech.* **569**, 365–398.

GRAY, J. M. N. T. & THORNTON, A. R. 2005 A theory for particle size segregation in shallow granular free-surface flows. *Proc. R. Soc. Lond. A* **461** (2057), 1447–1473.

GUILLARD, F., FORTERRE, Y. & POULIQUEN, O. 2016 Scaling laws for segregation forces in dense sheared granular flows. *J. Fluid Mech.* **807**, R1.

HERGAULT, V., FREY, P., MÉTIVIER, F., BARAT, C., DUCOTTET, C., BÖHM, T. & ANCEY, C. 2010 Image processing for the study of bedload transport of two-size spherical particles in a supercritical flow. *Exp. Fluids* **49** (5), 1095–1107.

HSU, T.-J. & LIU, P. L.-F. 2004 Toward modeling turbulent suspension of sand in the nearshore. *J. Geophys. Res.* **109**, C06018.

JACKSON, R. 2000 *The Dynamics of Fluidized Particles*. Cambridge University Press.

JONES, R. P., ISNER, A. B., XIAO, H., OTTINO, J. M., UMBANHOWAR, P. B. & LUEPTOW, R. M. 2018 Asymmetric concentration dependence of segregation fluxes in granular flows. *Phys. Rev. Fluids* **3**, 094304.

LAFAYE DE MICHEAUX, H., DUCOTTET, C. & FREY, P. 2018 Multi-model particle filter-based tracking with switching dynamical state to study bedload transport. *Mach. Vision Applics.* **29** (5), 735–747.

LI, L. & SAWAMOTO, M. 1995 Multi-phase model for sediment transport in sheet-flow regime under oscillatory flow. *Coast. Engng J.* **38** (2), 157–178.

MAURIN, R. 2015 Investigation of granular behavior in bedload transport using a Eulerian–Lagragian model. PhD thesis, Univ. Grenoble Alpes.
MAURIN, R., CHAUCHAT, J., CHAREYRE, B. & FREY, P. 2015 A minimal coupled fluid-discrete element model for bedload transport. *Phys. Fluids* **27** (11), 113302.

MAURIN, R., CHAUCHAT, J. & FREY, P. 2016 Dense granular flow rheology in turbulent bedload transport. *J. Fluid Mech.* **804**, 490–512.

MAURIN, R., CHAUCHAT, J. & FREY, P. 2018 Revisiting slope influence in turbulent bedload transport: consequences for vertical flow structure and transport rate scaling. *J. Fluid Mech.* **839**, 135–156.

MAY, L. B. H., GOLICK, L. A., PHILLIPS, K. C., SHEARER, M. & DANIELS, K. E. 2010 Shear-driven size segregation of granular materials: modeling and experiment. *Phys. Rev. E* **81** (5), 051301.

MIDDLETON, G. V. 1970 Experimental studies related to problems of flysch sedimentation. In *Flysch Sedimentology in North America* (ed. J. Lajoie), pp. 253–272. Bussiness and Economics Science Ltd.

RECKING, A., FREY, P., PAQUIER, A. & BELLEUDY, P. 2009 An experimental investigation of mechanisms involved in bed load sheet production and migration. *J. Geophys. Res.* **114**, F03010.

REVIL-BAUDARD, T. & CHAUCHAT, J. 2013 A two-phase model for sheet flow regime based on dense granular flow rheology. *J. Geophys. Res.* **118** (2), 619–634.

RICHARDSON, J. F. & ZAKI, W. N. 1954 The sedimentation of a suspension of uniform spheres under conditions of viscous flow. *Chem. Engng Sci.* **3** (2), 65–73.

ROUX, J. N. & COMBE, G. 2002 Quasistatic rheology and the origins of strain. *C. R. Phys.* **3** (2), 131–140.

SAVAGE, S. B. & LUN, C. K. 1988 Particle size segregation in inclined chute flow of dry cohesionless granular solids. *J. Fluid Mech.* **189**, 311–335.

SCHWAGER, T. & POSCHEL, T. 2007 Coefficient of restitution and linear-dashpot model revisited. *Granul. Matt.* **9** (6), 465–469.

SMILAUER, V. et al. 2015 *Yade Documentation*, 2nd edn. The Yade Project. Zenodo.

THORNTON, A., WEINHART, T., LUDING, S. & BOKHOVE, O. 2012 Modeling of particle size segregation: calibration using discrete particle method. *Intl. J. Mod. Phys. C* **23** (08), 1240014.

THORNTON, A. R., GRAY, J. M. N. T. & HOGG, A. J. 2006 A three-phase mixture theory for particle size segregation in shallow granular free-surface flows. *J. Fluid Mech.* **550**, 1–25.

UTTER, B. & BEHRINGER, R. P. 2004 Self-diffusion in dense granular shear flows. *Phys. Rev. E* **69** (3), 031308.

VAN DER VAART, K., GAJJAR, P., EPELY-CHAUVIN, G., ANDREINI, N., GRAY, J. & ANCEY, C. 2015 Underlying asymmetry within particle size segregation. *Phys. Rev. Lett.* **114**, 238001.

WIEDERSEINER, S., ANDREINI, N., EPELY-CHAUVIN, G., MOSER, G., MONNEREAU, M. L., GRAY, J. M. & ANCEY, C. 2011 Experimental investigation into segregating granular flows down chutes. *Phys. Fluids* **23**, 013301.