Comparing and contrasting size-based particle segregation models

Over the last 12 years, numerous new theoretical continuum models have been formulated to predict particle segregation in the size-based bidisperse granular flows over inclined channels. Despite their presence, to our knowledge, no attempts have been made to compare and contrast the fundamental basis upon which these continuum models have been formulated. In this paper, firstly, we aim to illustrate the difference in these models including the incompatible nomenclature which impedes direct comparison. Secondly, we utilise (i) our robust and efficient in-house particle solver MercuryDPM, and (ii) our accurate micro–macro (discrete to continuum) mapping tool called coarse-graining, to compare several proposed models. Through our investigation involving size-bidisperse mixtures, we find that (i) the proposed total partial stress fraction expressions do not match the results obtained from our simulation, and (ii) the kinetic partial stress fraction dominates over the total partial stress fraction and the contact partial stress fraction. However, the proposed theoretical total stress fraction expressions do capture the kinetic partial stress fraction profile, obtained from simulations, very well, thus possibly highlighting the reason why mixture theory segregation models work for inclined channel flows. However, further investigation is required to strengthen the basis upon which the existing mixture theory segregation models are built upon.


Introduction
Granular materials in nature [10,14] and industry [14] often comprise highly polydisperse particle mixtures. The constituents of these mixtures can vary in size, density, inelasticity, shape, surface roughness, etc. When such polydisperse mixtures are subjected to external forces such as shaking, stirring or shearing [10], these mixtures often segregate, leading to complex pattern formations such as particle segregation. Several factors have been reported to be responsible for segregation or de-mixing in polydisperse mixtures, where individual studies confirm the influences of varying the size [48], density [40], inelasticity [4], shape [30] and surface roughness [43]. However, in rapid free-surface flows over inclined channels, the differences in size and density are the dominant factors [2,5,6] leading to particle segregation.
Kinetic sieving [9] is the dominant mechanism in dense granular flows. Although it is an easy-to-comprehend mechanism, its effects can be mind-boggling. In order to illustrate the idea of kinetic sieving, let us consider a size-bidisperse granular mixture flowing down an inclined channel. As the flow progresses, fluctuations in the local pore space cause particles to fall under gravity into the space/gaps that are created beneath them. As a result, smaller-sized particles fall easily into these gaps leading to their gradual percolation towards the base of the flow (i.e. in the direction of gravity). Simultaneously, force imbalances lever/squeeze particles towards the surface; this process was named as squeeze expulsion by Savage and Lun [31]. The combination of kinetic sieving and squeeze expulsion results in a net migration of large particles upwards and small particles towards the base. As a result, this simple mechanism results in stratified layers, which one terms as particle segregation.
In 2011, Fan and Hill [11] proposed an alternative kineticstress-driven mechanism for segregation. The model was originally derived for vertical chutes; however, it was later extended to include the gravity-driven mechanism [20,21] and then applied to mixture flows over inclined planes and in rotating drums. Their model is very similar to the gravitydriven models and still uses the ideas of kinetic sieving; however, it is driven by gradients in kinetic stress rather than lithostatic pressure as in the case of the gravity-driven mixture theory models. In the kinetic-stress-driven mechanism, all particles are squeezed away from regions of higher fluctuation energy. During this process, smaller particles filter through the matrix of other particles, analogous to the gravity-driven void-filling mechanism, resulting in a net migration of small particles towards regions of lower fluctuation kinetic energy. Previously, Windows-Yule et al. [50] experimentally investigated the competition between gravity-driven, kinetic-stress-driven, and other segregation mechanisms in axially non-uniform drums. Here, we consider a simpler scenario of size-bidisperse mixtures rapidly flowing down an inclined plane.
In rapid free-surface flows, opposing kinetic sieving is diffusive remixing, which causes the random motion of particles as they collide and shear over each other [22]. Based on the relative strength of these competing mechanisms, the mixture strongly or weakly segregates; the relative strength is captured by the segregation Peclet number [18]. However, both experiments [48] and particle simulations [38] have shown that, in rapid chute flows, the effect of diffusive remixing is around 10 % of the strength of segregation.
Apart from kinetic sieving, which is a purely size-based effect, buoyancy effects due to differences in particle density also play a major role in particle segregation [23]. For bidisperse mixtures, varying in size and density, experiments [13] and particle simulations [41] indicate that a balance between the two driving mechanisms is possible, i.e. kinetic sieving and buoyancy effects, which in turn can keep the mixture homogeneously mixed.
Although complete understanding of the dynamics of segregation is beyond the scope of this paper, it becomes increasingly difficult as one allows more particle properties to be varying, i.e. size, density, shape, roughness, etc. In order to carry out the comparison between the existing particle segregation continuum models, we focus on the leading-order effect of bidisperse mixtures varying in size only. As an alternative to experiments, we employ discrete particle simulations (DPMs) [24] from which macroscopic quantities, appearing in the continuum models, are extracted. To perform this accurately, we use an appropriate micromacro mapping technique called coarse-graining [1,16], as this method gives continuum fields that exactly satisfy the continuum equations. The method has previously been extended to include the effects of boundaries or discontinuities [47] and more recently to unsteady polydisperse mixtures [42]. The technique has been successfully applied to investigate monodisperse shallow granular flows [46] and size-bidisperse mixtures [38,45]. However, Weinhart et al. [45] focussed on steady bidisperse flows alone; here, we consider both steady and transient data. For our simulations, we use our in-house open-source code MercuryDPM [36,37,39] which includes all the coarse-graining tools utilised in this paper.
The following section, Sect. 2, gives a brief review of the main theoretical models that have been proposed to describe size-segregation. In this section, special care is taken to express all the continuum segregation models using a uniform notation, thereby facilitating direct comparison between them. A look-up table, see Table 1, is also provided which relates the different notations used in each of these models. On the other hand, Sect. 3 contains a description of the setup of our particle simulations, which we utilised to compare and contrast the theoretical segregation models. More importantly, to perform the discrete-to-continuum (micro-macro) step, Sect. 4 lists the coarse-graining (CG) expressions, which are utilised to obtain the required macroscopic information from the discrete data. This section is provided because in the case of measuring the pressurescaling functions, which are utilised in the segregation models, it is vitally important that the micro-macro step is done in a consistent way; it is very easy with other micromacro methods to produce misleading and even erroneous results for the case of size-segregating flows. Section 5 explains how to apply the coarse-graining formulae from Sect. 4 to microscopic particle data obtained from the particle simulations setup in Sect. 3. The results are presented in Sect. 6, where the different segregation models presented in Sect. 2 are compared and contrasted. Conclusions and future directions are presented in Sect. 7. Moreover, Sects. 3, 4 and 5 are provided for completeness, so that this paper contains a full description of the tools required to validate size-based particle segregation continuum models; however, if you are only interested in the actual validation of the models, Sects. 2, 6 and 7 are written such that they can be treated as a standalone article.

Theoretical models
In this section, we present a brief overview of the existing size-based segregation continuum models, which have been formulated in the past few decades.
Based on the understanding of percolation and diffusion, Bridgwater et al. [3] were the first to formulate a continuum model quantifying particle segregation in a size-bidisperse mixture. Given x, y and z denotes the down-slope, crossslope and depth direction (Fig. 1), their governing equation was in terms of the volume concentration of small particles expressed as a fraction of the solid volume (φ), Fig. 1 A snapshot of a bidisperse mixture, with particle size-ratiô s = 1.3, flowing in a periodic box inclined at 26 • to the horizontal (discrete particle simulation). Colours/shades indicate the base/boundary (yellowish green, F b ), species type-s and type-l (blue, F s and red, F l ). We define the bulk as F s ∪ F l . Here, s and l denote small and large constituents where t is the time, q and D are the segregating velocity and diffusion rates. However, they soon realised that the rate of segregation q, i.e. segregating velocity, was dependent on the shear rate, the particle size-ratio and the normal pressure [3]. A few years later, Savage and Lun [31] used statistical mechanics and information entropy theory to arrive at a segregation model from the first principles. Their model was formulated in terms of number densities and fluxes. Although the model of Savage and Lun [31] considered various functional forms for the shear rate, i.e. different downslope velocity profiles u(z), it certainly had a downside because their model predicts segregation even in the absence of gravity, which is odd given kinetic sieving is a gravity-driven process.
Almost a decade later, from a different perspective, Dolgunin and Ukolov [8] developed a model using an equivalent mass transfer equation, which accounts for the granular mass transfer due to convection, quasi-diffusion and segregation. Their resulting governing equation is very similar to the form of Bridgwater et al. [3], i.e., where u is the down-slope flow velocity. Although, the above model (2) had all the features essential for describing particle segregation, a general framework to derive such models was still missing. In 2005, Gray and Thornton [19] proposed a general continuum framework to model size-based segregation by utilising the principles of mixture theory [29]. Gray and Thornton [19] postulated that in a gravity-driven process of segregation, larger particles bear more of the pressure or normal stress relative to their local concentration in comparison to the smaller particles, which are busy percolating into the voids underneath. This postulate is quantified by subdividing or partitioning the pressure/normal stress among the two constituents (e.g. large and small), where the subdivision is done by the so called pressure-scaling functions [19]. Thereby, the resulting framework governing sizebased particle segregation, including diffusive remixing [18], follows as where ∇ := [∂/∂ x, ∂/∂y, ∂/∂z] T and F (φ) is the function that specifies the segregation flux as a function of the local volume fraction of small particles. The exact details concerning the forms of q and F(φ) depend on the physical assumptions and the choice of the so-called pressure-scaling functions. Thereby, the focus of this work is to compare and contrast these different forms of pressure-scaling functions, using particle simulations and advance micro-macro tools. In the original model of Gray and Thornton, q was taken to be a constant, D = 0, and F(φ) = φ(1 − φ). The framework (3) is itself clearly a simplification, as diffusion only acts in the z-direction and the segregation also occurs in the z-direction, i.e. not the direction defined by gravity. However, the original framework was soon extended to include the effects of diffusive remixing [18], D = 0, and interstitial fluids [35]. More recently, developments by [15,26,41] have introduced revised forms of scaling functions which in turn has led to different forms of F(φ). The origins of these scaling functions are explained, documented and compared in the following section, Sect. 2.1. In 2010, May et al. [27,28] extended the Gray and Thornton model to include non-uniform (exponential) shear profiles (i.e. by modifying q in (3)). Assuming no diffusion, they modified the segregation rate q to depend on the vertical direction (z) and be proportional to the shear rate, γ (z) = du/dz. Based on these assumptions, they proposed the following functional form where s is a dimensionless segregation parameter with q 0 and λ being fitting parameters. However, the extended model could only capture the qualitative features in comparison to their experimental findings. In 2012, Marks et al. [26] significantly extended the particle segregation continuum theories to polydisperse flows, which also allows for varying density and size effects. Besides this, they also attempt to incorporate the shear-rate dependency in a consistent way. By doing so, they include the dependence of percolation velocities on spatially varying shear rate and particle size-ratio, i.e. q in (3) is a function of both shear rate (γ ) and particle size-ratio ( s). Moreover, they were also the first to attempt to quantify the segregation flux, F(φ), in terms of the real particle properties such as the particle size-ratio ( s) itself.
New extensions to the Gray and Thornton [19] continuum framework were added, in the year 2014, by Tunuguntla et al. [41], where they make subtle fundamental changes to the basis upon which Gray and Thornton [19] model is built upon. However, this did not alter the resulting framework presented in (3). Nevertheless, both, Tunuguntla et al. [41] and Gajjar and Gray [15] presented different possible forms for segregation fluxes, F(φ), which intended to quantify segregation in more realistic scenarios. Besides these extensions, the Gray and Thornton framework was further developed to accord for the effects due to density differences [17,26,41], thus enabling us to predict particle segregation in a wider range of applications.
During the same period, 2014 to present, an alternative transport equation approach is adopted by Leuptow and co-workers [12,32,33,51], who use the same general framework (3) to model particle segregation in two-dimensional bounded heaps, circular drums and chute flows. Similar to several extensions of mixture theory segregation models, listed above, they have also extended their models to account for polydisperse size effects [33] and density differences [51]. However, to close the model, they utilise particle simulations to determine the flow kinematics and physical parameters such as the incompressible bulk velocity field and diffusion coefficient. Similar to Marks et al. [26] and Tunuguntla et al. [41], Leuptow and co-workers [12,32,33,51] also consider the percolation velocity to be dependent on the spatially varying shear rates and particle size-ratios. This approach of closing their models with particle simulations does produce relatively good results; however, it is unable to capture flow transitions that lead to stratification patterns [25,49], which the mixture theory models are able to predict. Moreover, it should be noted that they use a binning method to extract their continuum fields, which are required to close their models. The binning method has extra degrees of freedom compared with the coarse-graining method we utilise here, namely how to split the stress between the large and small particles. Making different choices of this split has a large effect on the results, and crucially incorrect splitting of the stress can even change the directions of segregation; see the discussion in Staron and Phillips [34] for more details. For this reason, we utilise our coarse-graining method for bidiperse mixtures, which is summarised in Sect. 4.
In this paper, we focus on size-based mixture theory segregation models and utilise particle simulations as a validation tool alone. The key idea behind mixture theory segregation models is that particle segregation is caused by a gradient in lithostatic pressure caused by gravity, whereas in the case of the Fan-Hill model [11] segregation is caused by a gradient in kinetic stress. Thus, through the use of particle simulations and advance micro-macro tools, we scrutinise and quantify each of this mechanism and, more importantly, see how these two mechanisms play a role in the process of size-based particle segregation. Moreover, as implied from this paper's title, we also compare and contrast the different proposed forms of these models. The ultimate aim is to develop a theoretical model that can accurately predict particle size-segregation, thus eliminating the need for particles simulations entirely. However, the work in this paper is only a stepping stone towards this goal.
As a first step towards our analysis, in the following section, we review the basic background theory upon which these mixture theory models are based on and propose an unified notation to make model comparisons easier, both, for us and for future research.

Mixture theory framework
Mixture theory deals with partial variables that are defined per unit volume of the mixture rather than with the intrinsic variables associated with the material, i.e. the values one would measure experimentally, such as the material density of glass or steel particles.
The basic mixture postulate states that every point in the mixture is simultaneously occupied by all constituents. Hence, at each point in space and time, there exist overlapping fields (displacements, velocities, densities) associated with different constituents.
Since each constituent is assumed to exist everywhere, a volume fraction Φ ν is used to represent the percentage of the local volume occupied by constituent ν. Clearly, for sizebidisperse mixture, where Φ a denotes the fraction of volume corresponding to interstitial pore space filled with a passive fluid, e.g. air. However, for convenience, studies often consider volume fraction of the constituents per unit granular volume rather than per unit mixture volume, e.g. [35]. As the volume fraction of granular constituents per unit mixture is the volume fraction of each constituent per unit granular volume is defined as which also sum to unity, For each individual constituent, conservation laws for mass, momentum, energy and angular momentum can all be obtained, but here for simplicity, we only consider mass and momentum balance for bulk constituents and ignore the interstitial fluid effects. Each bulk 1 constituent satisfies the fundamental laws of balance for mass and momentum [29]. For a flow down a plane inclined at a constant angle θ , these balances are given in terms of partial flow quantities, such as partial stress, σ ν , density, ρ ν , and velocity, u ν = [u ν , v ν , w ν ] T , in the three coordinate directions, corresponding to each constituent indexed ν = s, l. g = (g t , 0, −g n ) T is the gravity vector, with g being the standard acceleration due to free fall; g t = g sin θ and g n = g cos θ . Additionally, the variable β ν represents the interspecies drag force due to resisting motion between the constituents. As these are internal forces residing within the granular mixture, from Newtons' third law the sum of these drags must be zero, i.e., β s + β l = 0. Furthermore, the bulk density ρ, barycentric granular velocity, u and the bulk stress σ are defined as The intrinsic variables defined for each of the constituents also play an integral role in the constitutive theory. These quantities are related to the partial quantities and hence are the key features of mixture theory. The intrinsic density of the constituent ρ ν * , i.e. the mass of the constituent per unit volume of the constituent, is related to the partial density by constituent volume fraction φ ν . The same relation applies to the partial and intrinsic stresses, σ ν and σ ν * , of the constituents. However, in standard mixture theory the partial velocity, u ν , of the constituent is identical to the intrinsic velocity, u ν * , of the constituent, where all the intrinsic quantities are denoted with ' * '. Note, as here we define the density per unit granular volume, φ ν , not per unit mixture volume Φ ν , ρ ν * is still a bulk density (i.e. it is averaged over the grain volume plus the air volume), hence it is the bulk density of glass or steel particles in a nonmixed state. 1 The bulk is defined as F s ∪ F l , see Fig. 1, excluding the interstitial pore-space.

Gravity-driven segregation
Most of the mixture flows involved in industrial and geological applications are shallow in nature, implying that the flow quantities in the down-and cross-slope directions (i.e. xand y-directions) are nearly uniform. Moreover, we assume that the partial densities and momenta become quasi-steady even before the flow segregates, implying that the temporal derivatives ∂ t (ρ ν ) and ∂ t (ρ ν u ν ) vanish after a certain equilibrium time t e . Therefore, we arrive from the momentum equation (9) where the subscript denotes the tensorial/vectorial component. Summing (11) over each particle species ν = s, l for α = z, setting σ zz | z=∞ = 0 implies that the flow is in lithostatic balance As stated earlier, the key idea behind the gravity-driven segregation models [15,17,19,41] is represented by the stress fractions 2 ( f ν ). These stress fractions determine the amount of stress to be distributed among each of the constituents, i.e., as smaller particles percolate downwards through the granular matrix, they carry less of the weight, thus causing the larger particles to proportionately carry more of the weight. In standard mixture theory, the constituent normal stress or pressure is assumed to be linearly related to the bulk normal stress or pressure through the volume fraction, i.e. f ν is assumed to be equal to φ ν . We have, however, a crucial deviation from the standard approach in order to account for the effects of segregation. From the partial stress relation in (10), we have where f ν * is defined as an intrinsic stress fraction and physically f ν * is the fraction of the total stress bore by constituent ν. On the other hand we define f ν as the over/under stress, which indicates the amount of excess stress bore by larger/smaller particles and is a convenient mathematical construct. From (13), it is clear that f ν is simply related to f ν * as f ν = φ ν f ν * . Previously, models have been formulated both in terms of f ν and f ν * and more confusingly both Hill and Tan [21] have been labelled simply f ν . In Table 1, we show the current notation (13) compared with the ones used in previous papers on modelling segregation. For a monodisperse "mixture", the whole mixture weight is supported by the mixture itself, i.e.
This is taken into account by defining the intrinsic stress fraction, f ν * , to take the following functional form, see the work of Tunuguntla et al. [41], where the brackets [ ] denote a functional dependence. The parameter B ν [φ ν ] denotes a material parameter, which is in general a function of the local partial volume fraction of the species type-ν. As we are restricting our attention to just two mixture components (small or large), from this point onwards, without any loss of generality, we will simply define φ s := φ. Thereby, based on the definition of stress fraction (15), Tables 2 and 3 list the intrinsic scaling functions ( f * ν ) used in previous segregation models. In both the tables, variables b, A γ and γ are material or fitting parameters witĥ s being the particle size-ratio d l /d s . Gajjar and Gray [15] determined the values of A γ and γ such that the maximum of the two stress fractions f l gt and f l gg are the same, where f l gt and f l gg denote the stress fractions used in the models of Gray and Thornton [19] and Gajjar and Gray [15]. However, as mentioned by Gajjar and Gray, the values of A γ and γ are dependent upon the actual particle properties. Note: f l gt = φ l f l * gt and f l gg = φ l f l * gg . Besides satisfying the conditions in (14), an important constraint is further imposed by the definition of bulk stress (10), where as is clear from Tables 2 and 3.
Furthermore, all the gravity-driven segregation models [15,17,19,41] consider the interaction drag or interconstituent friction to take the form of Darcy's law, which on neglecting diffusive remixing is where c is an inter-constituent drag coefficient. The interparticle surface interaction force is given by σ ∇ f ν , thus ensuring that segregation is driven by the partial normal stress (≈pressure) gradients. Substituting the expressions for the drag force (17) into the normal momentum balance equation for constituent type-ν (11)-neglecting normal acceleration terms-results in relative percolation velocities between the constituents and the bulk, The above equation illustrates the key idea behind the gravity-driven segregation models, where the relative velocities, between the constituents and the bulk, are dependent on the normal stress or pressure gradients induced by gravity. As ρ * s = ρ * l in our size-bidisperse mixture, on substituting the condition of lithostatic balance (12) in the above expression, (18) is restated in terms of constituents volume fraction as On the other hand, recently, Hill and Tan [21] presented a theory that explicitly combines, both, gravity-and kineticstress-driven mechanisms by modifying the interspecies drag as below where σ con zz and σ kin zz denotes the contact and kinetic stress 3 . Similarly, substituting the above Fan-Hill expression for the drag force (20) into the normal momentum balance equation for constituent type-ν (11) results in the relative percolation velocities as  Gray and Thornton [19] 1 Gajjar and Gray [15] 1 Gajjar and Gray [15] 1 Note that we assume the model from Bridgwater et al. [3] to be a mixture theory model, for which we back compute the pressure scaling for the larger constituent of the mixture Thus, the above equation indicates that the relative percolation velocities responsible for segregation, are dependent on, both, the kinetic stress and the total normal stress gradient. On further simplification, (21) is restated as Moreover, from the above equation, one could also define a dimensionless number as a ratio of its first and second term as The above dimensionless ( V r ) number can also be looked at as a ratio of the kinetic and normal stress gradients, ∂σ kin zz /∂z and ∂σ zz /∂z, thus allowing us to measure the relative strengths of the two segregation mechanisms.
For gravity-driven mechanism, substituting the constituents percolation velocity (19), e.g. for ν = s, in the mass balance equation (9) 1 , gives us back the earlier stated segregation governing equation in terms of the volume fraction of type-s, without accounting for diffusive remixing, Note that the form of segregation flux, F(φ), is dependent on the choice of the pressure scaling functions, f ν , see the works of [15,17,41].
In the following sections, we verify and compare the existing forms of scaling functions, listed in Tables 2 and 3, which subdivide the bulk pressure among the constituents. This is done by utilising information rich discrete particle simulations, which are set up as described in following section.

Simulation setup
Fully three-dimensional (3D) discrete particle simulations (DPMs) are used, as an alternative to experiments, to investigate segregation dynamics in a size-bidisperse mixture flowing over inclined channels. The simulations are set up in our in-house open-source particle solver, MercuryDPM [36,37,39].
To begin with, we consider a cuboidal box inclined at 26 • to the horizontal and is periodic in x-and y-direction. The box has dimensions L × W × H = 30d m × 10d m × 10d m , where d m is the mean particle diameter defined as d m = φ s d s + φ l d l . To create a rough base (bottom), we fill the box with a randomly distributed set of particles with uniform diameter d m and simulate them until a static layer of about 12 particles thickness is produced. Then a slice of particles with centres between z ∈ [9.3, 11]d m is taken and translated 11 mean particle diameters downwards, to form the rough base of the box. To ensure no flowing particles fall through the base, a solid wall is placed underneath this static layer.
Once the rough base is created, the box is inclined and filled with a homogeneously mixed bidisperse mixture of particle diameters d s and d l and equal material densities, i.e. ρ s * = ρ l * , as illustrated in Fig. 1; see Weinhart et al. [46] for more details.
In our DPM simulations, we non-dimensionalise the parameters such that the mean particle diameter d m = 1, the mean particle mass m m = 1, the magnitude of gravity g = 1. This implies that the mean particle density ρ m = m m / V m = 6/π and the mean particle volume V m = π( d m ) 3 /6. The nondimensional quantities are denoted by ' ' . In this paper, we consider three particular bidisperse mixtures withŝ = d l /d s = {1.3, 1.5, 1.7} without any size distribution around its particle size. Hence, the use of term perfectly in the subtitle of the paper.
Finally, we fill the box with the size-bidisperse mixture comprising particles of species type-s and type-l with V box = 30 × 10 × 10 being the non-dimensional volume of the box. The formulae in (25) enforce (i) the dimensionless flow height H to be the same in simulations, when the particle size-ratio is varied and (ii) the ratio of the total volume of species type-s over the total volume of the particles to be φ s , see Appendix 1 for details. Using (25) with homogeneous mixture initial conditions (randomly mixed) and equal particle volume fraction φ s = φ l = 0.5, DPM simulations for different values ofŝ are carried out. For the given box volume and size ratios, the flow thickness is approximately 16 mean particle diameters. Furthermore, a linear spring-dashpot model is used, where the spring stiffness and dissipation for each collision is chosen such that the collision/contact time t c = 0.005 √ d m /g and the coefficient of restitution r c = 0.88 are constant. The microscopic sliding friction coefficient is taken to be 0.5 and no rolling friction is considered. More details about the model can be found in [7,24,46]. Besides the contact model, we use the velocity-Verlet time-stepping algorithm.
Once the particle size-ratio, total number of particles and the contact model parameters are given as an input, the particles are inserted into the box with dimensions V box = 30 × 10 × H where H is defined as (N s + N l )/300. If the inserted particle at any position overlaps with another particle, the insertion is rejected and the insertion domain is enlarged by increasing H to H + 0.01 to ensure that there is enough volume for all the particles, thus leading to a loosely packed mixture initially. Once the simulation starts, the mixture compacts enough, see Fig. 1, giving the particles enough energy to initialise flow. For more details, see Weinhart et al. [46]. Given the particle simulations are setup in the above described manner, we still need to extract continuum fields to compare different forms of stress fractions, f ν , listed in Tables 2 or 3. This is the focus of the following section, i.e. how to perform the micro-macro step accurately?

Micro to Macro: coarse-graining (CG)
Compared to other, simpler methods of averaging such as binning or the method of planes, the coarse-graining method has the following advantages: (i) the resulting macroscopic fields exactly satisfy the equations of continuum mechanics, even near the base of the flow, see Weinhart et al. [47], (ii) the particles are neither assumed to be spherical or rigid, (iii) the resulting fields are even valid for a single particle, as no averaging over an ensemble of particles is required, (iv) the fields are determined at every point in space, not just at the centre of averaging cells as in the case of binning and (v) in a contact between different types of particles i.e. large and small here, the stress-partition is clearly defined. However, the coarse-graining method does assume (i) each particle pair to have a single point of contact, i.e. the particles are convex in shape; (ii) that the contact area can be replaced by a contact point, implying the particles are not too soft; and (iii) that collisions are not instantaneous (i.e., particles cannot be perfectly rigid).
Considering the above advantages and assumptions, in this section, we briefly elucidate the idea behind the coarsegraining technique and, more importantly, the mapping expressions, which will be employed to extract continuum partial densities, velocities, stresses and the interaction force density (interspecies drag force) from the discrete particle simulations setup in Sect. 3. For more information regarding the technique, please see Tunuguntla et al. [42], where they not only derive the coarse-graining expressions systematically but also focus on its application in detail. More importantly, they also present a general mixture-theory-based CG framework that can be easily extended to polydisperse mixtures without any loss of generality.

Nomenclature
Given that we have three different types of constituents: (bulk) type-s, (bulk) type-l and boundary, interstitial porespace of which is filled with a zero-density passive fluid; see Fig. 1. Each particle i ∈ F, where F := F s ∪ F l ∪ F b , will have a radius a i , centre of mass of which is located at r i with mass m i and velocity v i . The total force f i (26), acting on a particle i ∈ F is computed by summing the forces f i j due to interactions with the particles of the same type j ∈ F ν and other type, j ∈ F/F ν , and body forces b i , e.g., gravitational forces (m i g).
For each constituent pair, i and j, we define a contact vector r i j = r i − r j , an overlap δ i j = max(a i + a j − r i j · n i j ,0), where n i j is a unit vector pointing from j to i, n i j = r i j /|r i j |. Furthermore, we define a contact point c i j = r i + (a i − δ i j /2)n i j and a branch vector b i j = r i − c i j , see Fig. 2.
Irrespective of the size of constituent i and j, for simplicity, we place the contact point, c i j , in the centre of the contact area formed by an overlap, δ i j , which for small overlaps has a negligible effect on particle dynamics. In the following sections, we first present the idea of coarse-graining (CG) and then list the CG expressions for the partial and bulk quantities, using the above nomenclature.

Idea behind coarse-graining
To illustrate the idea, we consider the partial microscopic (point) mass density for a system (in a zero-density passive fluid) at point r and time t. From statistical mechanics, it is given as where δ(r) is the Dirac delta function in R 3 . This definition complies with the basic requirement that the integral of the mass density over a volume in space equals the mass of all the particles in this volume.
To extract the partial macroscopic mass density field, ρ ν (r, t), the partial microscopic mass density (27) is convolved with a spatial coarse-graining function ψ(r), e.g. a Heaviside, Gaussian or a class of Lucy polynomials 4 , thus leading to The result is equivalent to replacing the delta function with a spatial coarse-graining function (that is positive semidefinite, integrable, and has finite support), ψ(r), also known as a smoothing function. For simplicity, seen later, we define ψ i = ψ(r − r i (t)).

Coarse-graining expressions: novel micro-macro map
Using the same idea as explained in the previous section, expressions for partial quantities corresponding to constituent type-ν are Density: Velocity: u ν = P ν /ρ ν , Total partial stress: σ ν = σ kin,ν + σ con,ν , Kinetic stress: Contact stress: where in the kinetic stress expression, v i is the fluctuation Furthermore, in the contact stress expression, b i j denotes the branch vector origination from the particle centre to contact point, as illustrated in Fig. 2. Ψ i j denotes a line integral along the branch vector b i j , Ψ i j = 1 0 ψ(r − r i + sb i j )ds, which ensures the distribution of the force (26), between two constituents i and j to the partial stresses to be proportional to the length of the branch vectors. In other words, the stresses are distributed proportionally, based on the fraction of the branch vectors contained within each constituent. Thus, for contacts between a small constituent and a large constituent, the larger-sized constituent receives a larger share of the stress. All the above partial quantities are derived such that both the mass and momentum balance laws are exactly satisfied.
On utilising the above CG expressions, stated in Sect. 4.3, the following section focusses on extracting the continuum fields from, both, the transient and steady particle data of our size-bidisperse mixtures.

Applying coarse-graining to the DPM simulations
In order to compare and contrast the existing mixture theory segregation models, one needs to construct the continuum fields to compute the stress fractions or pressure scalings listed in Tables 2 or 3. Using the same coarse-graining expressions stated in Sect. 4.3, Weinhart et al. [45] performed the micro-macro step on their discrete particle data corresponding to steady state alone. Although their findings did illustrate important aspects of particle segregation, not much was inferred regarding the transient dynamics itself. Thereby, to understand transient segregation dynamics and also see how the stress distribution on each of the mixture constituents varies as a function of, both, space and time, micro-macro step must be performed on unsteady discrete particle data as well. This is the prime focus of this section, how to extract continuum fields from the available unsteady particle data?
In the simulations setup in Sect. 3, the temporal derivatives ∂ t (ρ ν ) and ∂ t (ρ ν u ν ) vanish after a short time interval  Fig. 3 For bidisperse mixture of particle size-ratioŝ = 1.3, evolution of the vertical centre of mass for both large (solid line) and small (dotted line) particles from unsteady to steady state is illustrated. Here,t n−1 tô t n+2 denote a range of points in time about which we would like to temporally average. Note that the points above are just to illustrate, in practice we average around more points than shown above and small particles, are tracked. Therefore, we focus on the time interval before the process of particle segregation has reached a steady state, i.e. whent ∈ [50,2000], where 'ˆ' denotes non-dimensional quantities such ast = t/ √ d m /g. For the purpose of our investigation, particle data is stored at every 5000 simulation time steps 5 .
As in Tunuguntla et al. [42], we use the coarse-graining expressions of Sect. 4.3 to spatially coarse-grain the particle data. For our analysis, we need continuum fields which are a function of, both, time (t) and flow depth (ẑ). To do so, for a given spatial coarse-graining scale ( w = w/d m ), we, first, spatially average the extracted continuum fields in x-and ydirection, thus resulting in averaged quantities,ζ (t,ẑ), as a function of both timet and flow depthẑ = z/d m . However, to construct macroscopic continuum fields in the temporal direction, we further need to averageζ (t,ẑ) temporally over a time interval t − w t ,t + w t , where w t is defined as the temporal averaging scale.
Given a spatial ( w) and temporal ( w t ) averaging scale, temporal averaging of any averaged fieldζ (t,ẑ) can be defined as ẑ)dt, for a given w and w t , (30) wheret denotes the point in time about which we would like to average, and w t determines the width of the averaging time interval, t − w t ,t + w t , see Fig. 3. As seen in (30), the coarse-graining expressions are completely dependent upon, both, spatial and temporal coarse-graining scales. Thereby, in Tunuguntla et al. [42], we focussed on determining optimal spatial and temporal coarse-graining scales, especially for segregating flows over inclined channels. Given these are already determined, we chose the spatial and temporal coarse-graining scale as w = 0.5 and w t = 50. With these at hand, coarse-grained fields are constructed at different timeŝ t n ∈ [100, 3000], see Fig. 3, for particle data corresponding to all the three size ratios.
Given these continuum fields, the following section brings us to the much-awaited section of this paper, where we illustrate and discuss our observations.

Analysis and discussion
With coarse-grained quantities available at the transient stages of particle segregation, we initially begin by looking at the evolution of the local solid volume fraction of, both, type-s and type-l constituents as a function of the flow-depth.

Local mass fractions
From the partial density (28), the partial mass fraction is defined as where ρ ν p is the (constant) material density of constituent type-ν. The bulk mass fraction is defined such that Λ = Λ 1 + Λ 2 .
Utilising these expressions, forŝ = {1.3, 1.5, 1.7}, Fig. 4 illustrates the evolution of local mass fraction for, both, the bulk and the mixture constituents. As seen, flows segregate faster with an increase in the particle size-ratio. Moreover, as the flow segregates, a pure layer of small particles develops at the base. Right above this layer, Λ l shows an oscillating behaviour at t n = [700, 2000, 3000] (Fig. 4), indicating that layers of large particles develop on the small particle bed. Once the flow is fully segregated, at t n = 3000, only a single layer right above the layer of pure small particles remains, see the circle highlighting this in the plots in the rightmost column of Fig. 4a-c. This might be an artefact of using monodisperse constituents, i.e. no size distribution around d s or d l . Apart from these oscillations, Λ l increases steadily towards the free-surface, forming a layer of large particles at the free-surface.

Transient vs. steady state analysis
As a step towards comparing the existing forms of stress fractions or pressure scalings, listed in Tables 2 or 3, we first utilise the computed total bulk and partial normal stresses, where σ zz = σ con zz + σ kin zz and σ ν zz = σ con,ν zz + σ kin,ν zz .
With the coarse-grained normal stress fields at hand, we first construct the total partial stress fractions as and then compute the partial stress fractions corresponding to the contact and kinetic stresses as Note, that from (34) and (32) it follows that f ν = f con,ν + f kin,ν .
The stress fractions determine the amount of normal stresses to be distributed among the small and large constituents. In gravity-driven segregation models, segregation is driven by the pressure gradient which in turn is scaled by the difference between the total partial stress fraction f ν and the partial volume fraction φ ν , i.e. f ν − φ ν , see (19). On the other hand, in the mechanism proposed by Hill and Tan, (22) implies that it is the difference between the contact and kinetic partial stress fractions ( f con,ν − f kin,ν ) and ( f con,ν − φ ν ) which scales the strength of either of the stress gradients, ∂σ kin zz /∂z and ∂σ zz /∂z. If |( f con,ν − f kin,ν )∂σ kin zz /∂z| > |( f con,ν − φ ν )∂σ zz /∂z|, segregation is majorly driven by gradient in kinetic stress. If not, then segregation is driven by the pressure gradients, as seen in (22). This is where our dimensionless number (23) comes into picture. Moreover, on re-examining (19) and (22), we also see that the pressure gradients are scaled by two pre-factors ( f ν − φ ν ) or ( f con,ν − φ ν ), respectively.
To distinguish or determine the profiles of these different pre-factors listed above, we plot the profiles of f ν − φ ν , f con,ν −φ ν , f kin,ν −φ ν and f con,ν − f kin,ν , for different particle size-ratioŝ := {1..3, 1.5, 1.7}, in Figs. 5, 6, 7. The plots are colour-coded-based on the flow depth and the markers • and denote large and small constituents. When closely observed, the plots illustrate the evolving profiles of these pre-factors, thus allowing us to improve our understanding regarding the forces behind segregation in the size-based segregation in inclined channel flows.
The pre-factors or the relative stress fractions of the large and small constituents shown in Figs. 5, 6, 7 are point symmetric, as f s = − f l and φ s = 1 − φ l . In the initial stages of segregation, i.e. att n = 100 to 400 DPM time units, both, the relative total and contact partial stress fraction profiles look very similar, complex and mostly concentrated around a small area between 0.25 < φ ν < 0.75, with values very close to zero except near the free-surface. This indicates that the fraction of the total and contact stress, (σ zz ) and (σ con zz ), bore by the large and small constituents is nearly the same as its local volume fraction, φ ν . Therefore, one can neglect f con,ν − φ ν , such that f con,ν ≈ φ ν . Note, this approximation was also previously used in the work of Hill and Tan [21]. By making this approximation, it gives cleaner data to fit for the kinetic-stress-driven mechanism; as it removes a lot of complex fine structure in the f con,ν profiles, which we believe may be artefacts of the perfectly bidisperse simulations, we considered here. On the other hand, the profiles of the relative kinetic partial stress fractions ( f kin,ν − φ ν ), corresponding to both large and small constituents, are significantly lower/higher than their local volume fractions, φ ν , in sharp contrast to their respective contact or total stress fraction profiles. Although the points near the free-surface (red diamonds and circles) do not exactly lie on an initially formed curve, as seen in the third column of each timestamp in Figs. 5,6,7, the values in the bulk and base of the flow (black, sky-blue, burnt-orange diamonds and circles) align themselves onto a different nonlinear curve.
As the process of segregation progresses,t n = 700 to 2000, the relative total and contact partial stress fractions appear to have unfurled from their complicated initial profiles to a much more structured one. The initially strong oscillations observed, near the free-surface (red circles and diamonds), in the profiles of f ν − φ ν and f con,ν − φ ν , disappear over time and the high concentration of data points observed at intermediate volume fractions, 0.25 < φ < 0.75 initially, also resolve over time. However, the fraction of the total and contact stress, (σ zz ) and (σ con zz ), bore by the large and small constituents is still nearly the same as its local volume fraction. Moreover, when closely observed, the profiles of relative kinetic stress fraction ( f kin,ν −φ ν ) and ( f con,ν − f kin,ν ) stay identical throughout the transient stages of segregation, i.e. the amount of kinetic stress borne by the small and large particles remains the same during the process of segregation (overlooking the points near the free-surface). Strikingly, this is also true compared with the corresponding profile in steady state att n = 3000. Thus, based on the illustrated profiles in Figs. 5, 6, 7, it implies that in a size-bidisperse mixture with given particle size-ratio, the smaller particles support a fraction of the kinetic stress larger than their volume fraction, as f con,ν ≈ φ ν , thereby complementing the findings of Weinhart et al. [45] and Hill and Tan [21].
Moreover, Fig. 8 further compares the steady state profiles of relative kinetic stress fraction ( f kin,ν − φ ν ) and ( f con,ν − f kin,ν ) for increasing particle size-ratio. As the particle size-ratio increases, the smaller-sized constituents support larger fraction of the kinetic stress and, interestingly, we also observe that the relative kinetic stress fractions become more asymmetric with the increase in particle sizeratio. However, more detailed study is required to explain this asymmetry.

Comparison of segregation models
In the previous section, we closely looked at the relative stress fractions, also known as the pre-factors in (19) and (22), computed from the simulations of a size-bidisperse mixture flowing over an inclined channel. In this section, we utilise these coarse-grained profiles and compare them with the existing theoretical forms of stress fractions that are listed in Table 2 and Table 3.
Given this, att n = 3000, we compared the total partial stress fraction profiles corresponding to the large constituents with the scalings listed in Table 3, see Fig. 9. As illustrated, with b = 1, γ = 0.9 and s = {1.3, 1.5, 1.7}, none of the expressions for the stress fractions seem to match the profiles computed from the particle simulations, ( f l sim ). Note that the values of b and γ could be modified such the pro-    Here, diamonds and circles correspond to small (sky-blue) and large (burnt-orange) constituents. Note: We assume particle segregation to have reached steady state, when the vertical centre of masses of, both, the large and small particles do not vary in time posed functional forms are closer to the relative total partial scaling profiles. However, this would not be of much help as the relative total partial scaling profiles obtained from the simulations ( f l sim ), for different bidisperse mixtures, are approximately the same as those of its local volume fraction, φ ν , thus implying negligible relative percolation velocity which in turn implies very weak segregation; see (19), thereby forcing us to ask a very simple question: How do these gravity-driven models even work, when none of the currently proposed scaling forms match in Fig. 9? Even if they do, which of them is relatively good at quantitatively predicting particle segregation? Fig. 9 Compares the relative total partial pressure scalings ( f ν − φ ν ), corresponding to large constituents (ν = l), of different theoretical models with the one obtained from our simulation, Sim. Here, GT , GG, Marks and T T , correspond to the scalings utilised in the gravity-driven segregation models of Gray and Thornton [19], Gajjar and Gray [15], Marks et al. [26] and Tunuguntla et al. [41] models whereas BW corresponds to the scaling back computed from the model of Bridgwater et al. [3]. The results correspond to size-bidisperse simulation with particle size-ratio (a)ŝ

The kinetic-stress model
From the previous sections, we can infer that smaller particles support a fraction of the kinetic stress larger than their local volume fraction. Given that, we postulate the intrinsic partial kinetic stress fractions ( f kin,s * ) to take the following form where B kin,ν , for example when ν = s, takes the form corresponding to the one for the larger constituent listed in Table 3. Basically, we imply that B kin,s := B l as listed in Table 4. By utilising these functional forms, in Fig. 10 we make a comparison between the postulated forms of partial kinetic stress fractions and the one obtained from our particle simulations. When looked upon closely, for a chosen set of fitting parameters, the functional forms proposed by Gray and Thornton and Gajjar and Gray match well with the ones computed from the particle simulations. In fact, with the correctly chosen fitting parameters, the form proposed Gajjar and Gray is in a very close agreement with the relative partial kinetic stress fraction obtained from the simulations. Although the scaling proposed by Marks et al. incorporates the particle size-ratio in its suggested form, it peaks at a different value of φ, see Fig. 10, when compared to the one obtained from simulation.
As a result, one could say that the gravity-driven segregation models work because the suggested functional forms for total partial stress fractions capture the shape of the partial kinetic stress fraction profiles obtained from the simulations. They implicitly support the fact that smaller particles support a fraction of kinetic stress larger than their volume fractions. As seen, three out of the four suggested functional forms are able to capture the segregation dynamics, with Gajjar and Gray [15] 1 Bridgwater et al. [3] 1 + φ 2 φ Fig. 10 For small constituents, this plot illustrates a comparison of relative partial kinetic stress fractions between the ones obtained from simulation and the one computed using the functional forms listed in Table 4 [19] and Gajjar and Gray (GG) [15] appear to closely capture the shape of the profile obtained from our simulation the closest being the form suggested by Gajjar and Gray [15]. However, there is still scope for work to be carried out in this direction, where fundamental changes are needed to be made in the basis upon which the mixture models are constructed. Moreover, on a different note, it also raises the interesting question of why did the model of Tununguntla et al. best determine the zero segregation line [41] if this model is incorrectly capturing the details of the stress distribution.

Summary and conclusions
In this work, we have reviewed the different segregation models that have recently been developed and have also defined a common unifying notation that allows for different models to be easily compared and contrasted. In order to do so, we utilise particle simulations combined with an advance micro-macro tool, called coarse-graining, to analyse the different assumptions made in these models. For the first time, we analyse data corresponding to transient stages of a segregating flow rather than simply using data corresponding to steady state.
Through our investigation, we found the following key findings -The kinetic stress bore by the large and small constituents remains the same during the whole process of segregation, except near the free-surface of the flow where we observe some fluctuations in the initial stages of the flow. -Near the free-surface of the flow, initially the stress (contact and kinetic) profiles are different; but, gradually they relax onto the steady-profile. This process happens from the centre outwards, i.e., the closer your location to the centre the quickly the stress relaxes to the steady-state value. -The contact stress has a much more complicated evolution; however, this seems to be associated with layering effects, that is present in the early stages of segregation. These layers slowly melt as time progress. -We confirm as previously reported by Weinhart et al. [45] and Hill and Tan [21] that, for the given particle size-ratio, the smaller constituents support a larger fraction of the normal kinetic stress. -With rightly chosen set of fitting parameters, the shape of the partial kinetic stress fraction is best captured by the functional form suggested by Gajjar and Gray [15]. -The measured relative stress fractions are asymmetric as observed in the experiments of van der Vaart et al. [44]. However, they considered varying fill fractions, whereas here we consider varying particle size-ratio. Moreover, the asymmetry increases with the increase in particle sizeratio.
The work of this paper could be extended and improved in several ways: -Firstly the strong layer we observed in the earlier stages of segregation, i.e. earlier time data, may be due to the use of a perfectly bidispersed mixture. Thereby, a small size distribution around two distinct mean values should be used. -Also the flows used here are of intermediate thickness; and deeper flows should be used to see the effect of absolute depth.
-None of the partial stress fractions suggested in the literature match the simulation results; however, three of the forms suggested, do closely match the profile of the kinetic-stress fraction obtained from the simulation. So this opens the question such as what is the correct functional form? More interestingly, why did the form of Tunuguntla et al. [41] capture the zero-segregating line so well? -Moreover, the effect of basal roughness and varying fill concentrations on the stress fractions also needs to be studied and included in the theoretical models.
It is clear that their still remains work do be done even for the simple case of size-bidisperse segregation in flows over simple inclined planes. Many different authors have recently made significant contributions to this topic and we believe the solution to fully understanding segregation will require the combination of several of these ideas. Therefore we hope this paper will serve as a useful reference point to understand and compare these distinct models as they often use different and inconsistent notation (relative to each other). taken to be equal for the simulations presented here, i.e.,