Coupling rheology and segregation in granular flows

Abstract During the last fifteen years there has been a paradigm shift in the continuum modelling of granular materials; most notably with the development of rheological models, such as the $\mu (I)$-rheology (where $\mu$ is the friction and I is the inertial number), but also with significant advances in theories for particle segregation. This paper details theoretical and numerical frameworks (based on OpenFOAM) which unify these currently disconnected endeavours. Coupling the segregation with the flow, and vice versa, is not only vital for a complete theory of granular materials, but is also beneficial for developing numerical methods to handle evolving free surfaces. This general approach is based on the partially regularized incompressible $\mu (I)$-rheology, which is coupled to the gravity-driven segregation theory of Gray & Ancey (J. Fluid Mech., vol. 678, 2011, pp. 353–588). These advection–diffusion–segregation equations describe the evolving concentrations of the constituents, which then couple back to the variable viscosity in the incompressible Navier–Stokes equations. A novel feature of this approach is that any number of differently sized phases may be included, which may have disparate frictional properties. Further inclusion of an excess air phase, which segregates away from the granular material, then allows the complex evolution of the free surface to be captured simultaneously. Three primary coupling mechanisms are identified: (i) advection of the particle concentrations by the bulk velocity, (ii) feedback of the particle-size and/or frictional properties on the bulk flow field and (iii) influence of the shear rate, pressure, gravity, particle size and particle-size ratio on the locally evolving segregation and diffusion rates. The numerical method is extensively tested in one-way coupled computations, before the fully coupled model is compared with the discrete element method simulations of Tripathi & Khakhar (Phys. Fluids, vol. 23, 2011, 113302) and used to compute the petal-like segregation pattern that spontaneously develops in a square rotating drum.


Introduction
Despite nearly all natural and man-made granular materials being composed of grains of varying size, shape and frictional properties, the majority of continuum flow modelling has largely been restricted to perfectly monodisperse aggregates. The purpose of this work is therefore to extend the current granular flow models by introducing multiple phases, with different properties, and to model inter-phase segregation. Coupling the flow rheology to the local constituent concentrations is important because the mobility of a granular flow is strongly affected by the local frictional properties of the grains. In turn, the bulk flow controls the strength and direction of the segregation as well as the advection of the granular phases.
Striking examples of segregation induced feedback on the bulk flow are found during levee formation (Iverson & Vallance 2001;Johnson et al. 2012;Kokelaar et al. 2014) and fingering instabilities (Pouliquen, Delour & Savage 1997;Pouliquen & Vallance 1999;Woodhouse et al. 2012;Baker, Johnson & Gray 2016b), which commonly occur during the run-out of pyroclastic density currents, debris flows and snow avalanches. Many other examples of segregation-flow coupling occur in industrial settings (Williams 1968;Gray & Hutter 1997;Makse et al. 1997;Hill et al. 1999;Ottino & Khakhar 2000;Zuriguel et al. 2006). Storage silo filling and emptying, stirring mixers and rotating tumblers all have the common features of cyclic deformation and an ambition of generating well-mixed material. However, experiments consistently suggest that these processes have a tendency to promote local segregation, which can feedback on the bulk flow velocities. Considering the inherent destructive potential of geophysical phenomena and the implications of poor efficiency in industrial mixing, a continuum theory which captures the important physics of flow and of segregation simultaneously is therefore highly desirable.
To date, the leading approaches for solving coupled flow and segregation have come from either discrete particle simulations (Tripathi & Khakhar 2011;Thornton et al. 2012) or from depth-averaged equations (Woodhouse et al. 2012;Baker et al. 2016b;Viroulet et al. 2018). Particle simulations, using the discrete element method (DEM), provide important rheological information as evolving velocities, stresses and constituent concentrations can be directly computed given only minimal approximations. Such results can then be used to motivate models for the bulk flow (GDR MiDi 2004;Jop, Forterre & Pouliquen 2006;Singh et al. 2015) and also to form connections between flow and segregation processes (Hill & Fan 2008;Staron & Phillips 2015). Unfortunately, the discrete particle approach is naturally limited by computational expense as many flows of interest include such a large number of particles that direct DEM calculations are unfeasible. Recently efforts have been made to overcome this limitation with the development of hybrid schemes (e.g. Yue et al. 2018;Xiao et al. 2019) which couple discrete particle dynamics to continuum solvers, but these approaches naturally invoke additional complexity and new assumptions are required in order to map properly and consistently between the somewhat disparate approaches.
Depth-averaged models, which reduce the full three-dimensional flow to two dimensions by integrating though the depth and assuming shallowness, lead to efficient numerical codes which are widely used in geophysical modelling (see e.g. Grigorian, Eglit & Iakimov 1967;Savage & Hutter 1989;Iverson 1997;Gray, Wieland & Hutter 1999;Pouliquen & Forterre 2002;Sheridan et al. 2005;Mangeney et al. 2007;Christen, Kowalski & Bartelt 2010;Gray & Edwards 2014;Rauter & Tuković 2018;Rocha, Johnson & Gray 2019). However, depth-averaged approaches are limited to geometries in which there is a clear dimension that remains shallow throughout the dynamics. This approximation holds well for thin flows on inclined planes and for flows over certain gradually varying terrain, but breaks down in many flows of practical interest, such as those in hoppers, silos and rotating drums.
Historical attempts to construct three-dimensional continuum models for monodisperse granular materials focused on quasi-static deformations and lead to elasto-plastic formulations of models such as the Drucker-Prager yield condition (Lubliner 2008) and critical state soil mechanics (Schofield & Wroth 1968). Despite successes in modelling the point of failure of materials under load, calculations of the subsequent time-dependent flow proved to be problematic, because the results are grid-size dependent. Schaeffer (1987) showed that this was because the underlying equations are mathematically ill posed, i.e. in the small wavelength limit the growth rate of linear instabilities becomes unbounded in certain directions.
Despite the Mohr-Coulomb/Drucker-Prager plasticity theory being designed for the flow of monodisperse grains, the grain diameter d does not appear in the constitutive model. It can be incorporated by making the friction μ a function of the non-dimensional inertial number, which is defined as whereγ is the shear rate, p is the pressure and ρ * is the intrinsic grain density (Savage 1984;Ancey, Coussot & Evesque 1999;GDR MiDi 2004). Jop et al. (2006) generalized the scalar μ(I)-rheology to tensorial form. The resultant incompressible μ(I)-rheology leads to a significantly better posed system of equations (Barker et al. 2015). For the μ(I) curve suggested by Jop, Forterre & Pouliquen (2005), the equations are well posed for a large range of intermediate values of I and are only ill posed for very low or relatively high inertial numbers.  derived a new functional form for the μ(I) relation, which is known as the partially regularized μ(I)-rheology. This ensures well posedness for 0 < I < I max , where I max is a very large value, and leads to stable and reliable numerical schemes. It also provides a better fit to experimental (Holyoake & McElwaine 2012; and DEM data (Kamrin & Koval 2012) than the original μ(I) curve, but also introduces a creep state (i.e. μ = 0 when I = 0) so the granular material no longer has a yield stress. It is possible to formulate well-posed models with a yield stress by introducing bulk compressibility Schaeffer et al. 2019) or non-locality (Henann & Kamrin 2013). However, in this paper the partially regularized μ(I)-rheology is chosen for the bulk flow, both for simplicity and because it is most readily compatible with existing numerical methods and particle segregation models.
Initially well-mixed granular materials have a strong propensity of ordering spatially when they undergo flow. Chief among these effects is that of particle-size segregation, made famous through the moniker 'the Brazil nut effect' (Rosato et al. 1987), whereby particles move relative to the bulk flow based on their size compared with their neighbours. The resultant vertical distribution, in which larger particles are often concentrated at the surface of a flow, can also be observed in many geophysical mass flows, forming so-called inversely graded deposits (e.g. Middleton 1970;Festa et al. 2015). The origin of this effect was explained through statistical entropic arguments by Savage & Lun (1988) who proposed a means of 'kinetic sieving' (Middleton 1970) in which smaller grains are more likely to fall (by gravity) into voids that are created as layers of particles are sheared over one another. Force imbalances then drive particles out of the denser layer, which is known as 'squeeze expulsion'. The combination of kinetic sieving and squeeze expulsion produces a net upward motion of large particles as the smaller grains percolate downwards.
These concepts formed the basis of the theory of Gray & Thornton (2005) who focused on this form of gravity-driven segregation in granular free-surface flows. The theory was later extended by Gray & Chugunov (2006), in order to account for diffusive mixing, and has been successfully applied to a range of gravity-driven flows (Gray 2018). However, Fan & Hill (2011) found that the direction of segregation was not always aligned with the vector of gravitational acceleration. Instead gradients in kinetic stress were found to drive and orient segregation in a range of geometries (Hill & Tan 2014). These findings have since inspired many investigations into the micromechanical origin of size segregation (Staron & Phillips 2015;Guillard, Forterre & Pouliquen 2016;van der Vaart et al. 2018), but a unified and compelling theory is still lacking.
In order to accommodate different models for size segregation and different flow rheologies, this paper first introduces a very general framework for multi-component flows in § 2. In particular, the multicomponent segregation theory of  is generalized to allow sub-mixtures to segregate in different directions and with differing diffusion rates. In § 3 the three primary coupling mechanisms are discussed in detail. Section 4 documents the general numerical method, which is then extensively tested against the one-way coupled simulations in § 5. Two-way fully coupled simulations are then presented for flow down an inclined plane, in § 6, and in § 7 simulations are performed for a square rotating drum. The new experimental segregation law of Trewhela, Ancey & Gray (2021) is tested against the steady-state DEM solutions of Tripathi & Khakhar (2011) in § 6.3 and then used in § 7 for the rotating drum simulations, which are able to spontaneously generate petal-like patterns that have previously been seen in the experiments of Hill et al. (1999), Ottino & Khakhar (2000) and Mounty (2007).

The partially regularized μ(I)-rheology for the bulk flow
The granular material is assumed to be composed of a mixture of particles that may differ in size, shape and surface properties, but have the same intrinsic particle density ρ * . If the solids volume fraction Φ is constant, which is a reasonable first approximation (GDR MiDi 2004;Tripathi & Khakhar 2011;Thornton et al. 2012), then the bulk density ρ = Φρ * is constant and uniform throughout the material. Mass balance then implies that the bulk velocity field u is incompressible where ∇ is the gradient and · is the dot product. The momentum balance is where p is the pressure, η is the viscosity, D = (∇u + (∇u) T )/2 is the strain-rate tensor and g is the gravitational acceleration. Assuming alignment of the shear-stress and strain-rate tensors the μ(I)-rheology (Jop et al. 2006) implies that the granular viscosity is where the second invariant of the strain-rate tensor is defined as and the inertial number, defined in (1.1), in this notation becomes The meaning of the particle size d in a polydisperse mixture will be clarified in § 3.2. Note that this paper is restricted to two-dimensional deformations with an isotropic Drucker-Prager yield surface. However, as shown by Rauter, Barker & Fellin (2020), this framework can be extended to include three-dimensional deformations through further modification of the granular viscosity i.e. dependence on det(D). The viscosity (2.3) is a highly nonlinear function of the inertial-number-dependent friction μ = μ(I), pressure p and the second invariant of the strain rate D . Barker et al. (2015) examined the linear instability of the system, to show that the growth rate becomes unbounded in the high wavenumber limit, and hence the incompressible μ(I)-rheology is mathematically ill posed, when the inequality is satisfied, where μ = ∂μ/∂I. Ill posedness of this type is not only unphysical, but results in two-dimensional time-dependent numerical computations that do not converge with mesh refinement (see e.g. Barker et al. 2015;Martin et al. 2017). If the friction is not inertial number dependent (μ = const.) the ill-posedness condition (2.6) is satisfied for all inertial numbers and the system of equations is always ill posed (Schaeffer 1987). The equations are also ill posed if the friction μ is a decreasing function of I, since all the terms in (2.6) are strictly positive. The original form of the μ(I)-curve proposed by Jop et al. (2005) is a monotonically increasing function of I starting at μ s at I = 0 and asymptoting to μ d at large I, where I 0 is a material specific constant. The inertial number dependence in (2.7) gives the rheology considerably better properties than the original, constant friction coefficient, The range of well posedness was extended by  to 0 ≤ I ≤ I max , where I max is a large maximal value, by changing the shape of the μ(I)-curve. This paper uses the μ(I)-curve proposed by     friction law, which were measured for 143 μm glass beads. The value I 1 I N 1 is set by the lower bound for well posedness in Jop et al.'s (2006) friction law using the parameters above. Unless stated otherwise, the remaining parameters are the values chosen in the numerical simulations. Note that the air viscosity is higher than the physical value of η a * = 1.81 × 10 −5 kg (ms) −1 to prevent the convective Courant number limiting the time-step size.  Jop et al. (2006) (red line) and the partially regularized law of  (blue line). The Jop et al. (2006) curve has a finite yield stress μ s (red dot) and asymptotes to μ d at large inertial number (dashed line). For the parameters summarized in table 1, it is well posed in the range [I N 1 , I N 2 ] = [0.00397, 0.28016] (red shaded region). A necessary condition for well posedness is that the friction μ is zero at I = 0 (blue dot).  curve therefore introduces a creep state for I ∈ [0, I 1 ] to the left of the green dot (see inset) and becomes linear at large inertial numbers. The value of I 1 = 0.004 is chosen to be very slightly larger than I N 1 . The resulting partially regularized law is well posed for I ∈ [0, 16.9918].
where α and μ ∞ are new material constants and is chosen to ensure continuity between the two branches at I = I 1 . As shown in figure 1 this curve stays close to (2.7) in the well-posed region of parameter space, but passes though μ = 0 at I = 0 and is asymptotically linear in I at large inertial numbers. For the parameters given in table 1, the matching occurs at I 1 = 0.004 (which is very slightly larger than I N 1 ) and the maximum well-posed inertial number is I max = 16.9918.
The partially regularized μ(I)-rheology not only ensures well posedness for I < I max , but it also provides better fitting to experimental and DEM results. For instance, relative to (2.7) the new μ(I)-curve (2.8) predicts higher viscosities for large values of I, as seen in the chute flow experiments of Holyoake & McElwaine (2012) and . For low values of I, the partially regularized μ(I)-rheology predicts very slow creeping flow, since μ → 0 as I → 0. This behaviour is seen, to a certain extent, in DEM simulations (Kamrin & Koval 2012;Singh et al. 2015) and has been postulated by Jerolmack & Daniels (2019) to play an important role in soil creep. The lack of a yield stress may, however, be viewed as a disadvantage of the theory. It is important to note that by allowing some bulk compressibility, it is possible to formulate granular rheologies that are always well posed mathematically Heyman et al. 2017;Goddard & Lee 2018;Schaeffer et al. 2019) and support a yield stress.

Generalized polydisperse segregation theory
The granular material is assumed to be composed of a finite number of grain-size classes, or species ν, which have different sizes d ν , but all have the same intrinsic density ρ ν * = ρ * . Note that the inclusion of density differences between the particles implies that the bulk velocity field is compressible, which significantly complicates the theory (Tripathi & Khakhar 2013;Gilberg & Steiner 2020) and is therefore neglected. Even for a bidisperse mixture of particles of the same density, the grains can pack slightly denser in a mixed state than in a segregated one (Golick & Daniels 2009). However, the DEM simulations (Tripathi & Khakhar 2011;Thornton et al. 2012) suggest these packing effects are small, and for simplicity, and compatibility with the incompressible μ(I)-rheology, these solids volume fraction changes are neglected. Each grain-size class is therefore assumed to occupy a volume fraction φ ν ∈ [0, 1] per unit granular volume, and the sum over all grain sizes therefore equals unity ∀ν φ ν = 1.
(2.10) Many models to describe particle segregation have been proposed (see e.g. Bridgwater, Foo & Stephens 1985;Savage & Lun 1988;Dolgunin & Ukolov 1995;Khakhar, Orpe & Hajra 2003;Gray & Thornton 2005;Gray & Chugunov 2006;Fan & Hill 2011;Schlick et al. 2015) and these all have the general form of an advection-segregation-diffusion equation where F ν is the segregation flux and D ν is the diffusive flux. Provided that these fluxes are independent, this formulation is compatible with the bulk incompressibility provided ∀ν F ν = 0, and The form of the segregation flux is motivated by early bidisperse models (Bridgwater et al. 1985;Dolgunin & Ukolov 1995;Gray & Thornton 2005). These all had the property that the segregation shut off when the volume fraction of either species reached zero. This is satisfied if the segregation flux for species ν and λ is proportional to φ ν φ λ . In polydisperse systems,  proposed that the segregation flux for species ν was simply the sum of the bidisperse segregation fluxes with all the remaining constituents λ. This paper proposes a significant generalization of this concept, by allowing the local direction of segregation to be different for each bidisperse sub-mixture, so that the segregation flux takes the general polydisperse form where f νλ is the segregation velocity magnitude and e νλ is the unit vector in the direction of segregation, for species ν relative to species λ. This segregation flux function satisfies the summation constraint (2.12a,b) provided f νλ = f λν and e νλ = −e λν . (2.14a,b) In contrast to the theory of  the segregation velocity magnitude is the same for species ν with species λ and species λ with species ν, and it is instead the direction of segregation that now points in the opposite sense. This approach has the property that individual sub-mixtures may segregate in different directions, which allows the theory to be applied to polydisperse problems where gravity-driven segregation (e.g. Gray 2018) competes against segregation driven by gradients in kinetic stress (Fan & Hill 2011). This would require the constituent vector momentum balance to be solved in order to determine the resultant magnitude and direction of segregation (Hill & Tan 2014;Tunuguntla, Weinhart & Thornton 2017). In this paper the inter-particle segregation is always assumed to align with gravity. However, the direction of segregation for the particles and air can be chosen to be different. This proves to be advantageous in the numerical method that will be developed to solve the coupled system of equations in § 4. It is also very useful in the numerical method to allow the rate of diffusion between the various sub-mixtures to be different. By direct analogy with the Maxwell-Stefan equations (Maxwell 1867) for multi-component gas diffusion, the diffusive flux vector is therefore assumed to take the form where D νλ is the diffusion coefficient of species ν with species λ. Equation (2.15) satisfies the summation constraint (2.12a,b), provided D νλ = D λν , and reduces to the usual Fickian diffusion for the case of bidisperse mixtures (see e.g. Gray & Chugunov 2006). For a mixture of n distinct species, it is necessary to solve n − 1 separate equations of the form (2.11) together with the summation constraint (2.10) for the n concentrations φ ν , assuming that the bulk velocity field u is given.
In the absence of diffusion, concentration shocks form naturally in the system (see e.g. Gray & Thornton 2005;Thornton, Gray & Hogg 2006;. The jumps in concentration across such boundaries can be determined using jump conditions that are derived from the conservation law (2.11) (see e.g. Chadwick 1976). These jump conditions are also useful when formulating boundary conditions with diffusion. The most general form of the jump condition for species ν is where n is the normal to the shock, v n is the normal speed of the shock and the jump bracket [[ ]] is the difference of the enclosed quantity on the forward and rearward sides of the shock. In particular, if the flow is moving parallel to a solid stationary wall, then the jump condition reduces to the one-sided boundary condition (2.17) This implies that the segregation and diffusive fluxes balance and that there is no mass lost or gained through the wall.
2.3. Reduction to the bidisperse case For the case of a mixture of large and small particles, which will be referred to by the constituent letters ν = s, l respectively, the summation constraint (2.10) becomes φ s + φ l = 1. (2.18) Assuming that the gravitational acceleration vector g points downwards and that the segregation aligns with this direction, the concentration equation (2.11) for small particles reduces to ∂φ s ∂t where f sl is the segregation velocity magnitude and D sl is the diffusivity of the small and large particles. The functional dependence of these quantities on the shear rate, pressure, gravity, particle size and the particle-size ratio, will be discussed in detail in § 3.3.

Coupling the bulk flow with the segregation
One of the key advances of this paper is to develop a coupled framework that solves for the bulk velocity field u, the pressure p and the particle concentrations φ ν at the same time. This framework allows us to explore some of the intimate couplings between the segregation and the bulk flow. A variety of couplings are envisaged, that may act singly or all at once, to generate very complex behaviour. The models fall into two classes: (i) one-way coupled and (ii) two-way coupled, and both forms of coupling are investigated in this paper.

Advection by the bulk flow field
Many important practical segregation problems involve a time-dependent spatially evolving bulk flow that cannot easily be prescribed or determined from DEM simulations. Since the particle concentrations are advected by the bulk velocity u, the most basic one-way coupling involves the solution of the mass (2.1) and momentum (2.2) balances to determine this velocity field. This enables the segregation equation (2.11) to be solved within a physically relevant flow field, provided the segregation velocity magnitudes and diffusivities are prescribed. Computations of this nature may give a good indication of where differently sized particles are transported, in a flow field that does not experience strong frictional feedback from the evolving species concentrations. This simplification implicitly assumes that an essentially monodisperse flow field provides a reasonable approximation for the dynamics of a much more complex polydisperse mixture of particles, and that there is no feedback of this local flow field on the segregation and diffusion rates. This simple coupling is investigated in § 5 for a time-dependent spatially evolving flow down an inclined plane. Importantly, this simple one-way coupling also enables the accuracy of the numerical method to be tested against known exact travelling wave and steady-state solutions for the bulk flow field and the particle concentrations. In general, the particle concentrations are always transported by the bulk flow field, so this mechanism is also active in models with more complex couplings, which will be investigated in § § 6 and 7.

Segregation induced frictional feedback on the bulk flow
Each distinct granular phase may have differing particle size, shapes or surface properties, that lead to different macroscopic friction and/or rheological parameters. In this next stage of coupling these rheological differences are built into the model, so that the evolving particle concentrations feedback on the bulk flow through the evolving macroscopic friction of the mixture. There are two basic ways to introduce this coupling.
A key finding of the μ(I)-rheology (GDR MiDi 2004) was that the inertial number (2.5) is a function of the particle size d. This is clearly defined in a monodisperse mixture, but an important generalization is needed for polydisperse systems. Based on DEM simulations of bidisperse two-dimensional assemblies of disks, Rognon et al. (2007) proposed an inertial number in which the particle size d was replaced by the local volume fraction weighted average particle sized. The same law was also proposed by Tripathi & Khakhar (2011) and shown to agree with three-dimensional DEM simulations of spheres. Generalizing this concept to polydisperse systems, implies that the average particle sizē (3.1) evolves as the local concentrations φ ν of each particle species change. As a result, given the same local shear rate 2 D , pressure p and intrinsic grain density ρ * , the new inertial number will be larger for a mixture composed of larger particles than one made of smaller grains. As well as differences in size, the particles may also differ in shape and/or surface properties. A prime example of this are segregation induced fingering instabilities, which develop with large angular (resistive) particles and finer spherical particles (Pouliquen et al. 1997;Pouliquen & Vallance 1999;Woodhouse et al. 2012;Baker et al. 2016b). The effect of particle shape and surface properties can certainly be modelled in monodisperse flows by changing the assumed macroscopic frictional parameters (see e.g. Pouliquen & Forterre 2002;Forterre 2006;Edwards et al. 2019;Rocha et al. 2019). Furthermore, the results of Baker et al.'s (2016b) granular fingering model suggest that a good approach is to assume that each phase satisfies a monodisperse friction law μ ν = μ ν (I) of the form (2.8) and then compute the effective friction by the weighted sum of these laws, i.e. μ = ∀ν φ ν μ ν . (3.3) On the other hand, it is also possible to assume that there is a single μ(I)-curve, given by (2.8), but that the parameters in it evolve as the mixture composition changes, i.e.
where μ ν s , μ ν d , μ ν ∞ and I ν 0 are the frictional parameters for a pure phase of constituent ν. There is clearly potential for a great deal of complexity here that needs to be explored. However, to the best of our knowledge there are no DEM studies that measure the effective frictional properties of mixtures of particles of different sizes, shapes and surface properties that could further guide the model formulation. Segregation mobility feedback on the bulk flow will be investigated further in § 6.

Feedback of the bulk flow on the segregation rate and diffusivity
The shear rateγ = 2 D , the pressure p, gravity g and the particle properties also enter the equations more subtly through the functional dependence of the segregation velocity magnitude f νλ and diffusivity D νλ in the fluxes (2.13) and (2.15). Even in bidisperse granular mixtures very little is known about their precise functional dependencies. However, dimensional analysis is very helpful in constraining the allowable forms.
Consider a bidisperse mixture of large and small grains of sizes d l and d s , respectively, which have the same intrinsic density ρ * . The small particles occupy a volume fraction φ s = 1 − φ l per unit granular volume and the total solids volume fraction is Φ. The system is subject to a bulk shear stress τ , a pressure p and gravity g, which results in a shear rateγ . Even though these variables are spatially varying, they are considered here as inputs to the system, whereas the segregation velocity magnitude f sl and the diffusivity D sl are outputs. Since there are nine variables, with three primary dimensions (mass, length and time), dimensional analysis implies that there are six independent non-dimensional quantities whered is the volume fraction weighted average grain size defined in (3.1), P is the non-dimensional pressure and R is the grain-size ratio. For a monodisperse system in the absence of gravity, only the first three quantities are relevant and GDR MiDi (2004) made a strong case for the friction μ being purely a function of the inertial number I. This led to the development of the incompressible μ(I)-rheology (GDR MiDi 2004;Jop et al. 2006;, which is used in this paper. Using the monodisperse scalings, it follows that in the absence of gravity the self-diffusion of grains should scale as where F is an arbitrary function of the friction, the inertial number and the solids volume fraction, and with no dependence on P, R and φ s . In both the incompressible and compressible μ(I)-rheologies (GDR MiDi 2004;da Cruz et al. 2005;Jop et al. 2006;Forterre & Pouliquen 2008) the friction μ and the solids volume fraction Φ are rigidly bound to the inertial number I, so it is not necessary to retain their dependence in (3.6). However, in the latest well-posed compressible theories (e.g. Heyman et al. 2017;Schaeffer et al. 2019) the μ = μ(I) and Φ = Φ(I) laws only hold at steady state, and so the general form of the diffusivity (3.6) applies. Utter & Behringer (2004) showed experimentally that the self-diffusion coefficient scaled with the shear rate and the particle size squared. This suggests that the simplest model for the diffusion of the grains in a polydisperse system is where A = 0.108 is a universal constant (Utter & Behringer 2004) andd is now interpreted as the average, locally evolving, particle size defined in (3.1). Some evidence for this is provided by the experiments of Trewhela et al. (2021) which show that a single small intruder in a matrix of large grains performs larger random walks than a single large intruder in a matrix of fine grains. In general, however, the diffusivity could be multiplied by an arbitrary function of the other non-dimensional quantities in (3.5a-f ). Gravity-driven percolation (kinetic sieving) and squeeze expulsion (Middleton 1970;Bridgwater et al. 1985;Savage & Lun 1988;Gray & Thornton 2005;Gray 2018) combine to create the dominant mechanism for segregation in dense sheared granular flows. Assuming that the segregation is independent of the diffusion, dimensional analysis suggests that the segregation velocity magnitude in a bidisperse mixture of large and small particles should scale as where G is an arbitrary function. It has long been known that the segregation velocity magnitude f sl is strongly dependent on the strain rate and the particle-size ratio (see e.g. Bridgwater et al. 1985;Savage & Lun 1988). Gray & Thornton (2005) also suggested that there should be a dependence on gravity. Evidence for this is provided by the fact that granular segregation experiments, with a density matched interstitial fluid, do not segregate (Vallance & Savage 2000;Thornton et al. 2006), i.e. when gravity is effectively reduced, so is the rate of segregation. Inclusion of the gravitational acceleration suggests that the segregation velocity magnitude should also be pressure dependent, since g only appears in the non-dimensional pressure P. This is supported by the experiments of Golick & Daniels (2009), who observed a dramatic slowing in the segregation rate when they applied a normal force on their ring shear cell. This pressure-dependent suppression of segregation has been investigated further in the DEM simulations of Fry et al. (2018), who suggested that the segregation velocity magnitude should scale with the reciprocal of the square root of the pressure. When this is combined with the shear-rate dependence this implies that f sl is linear in the inertial number. In this paper, the segregation velocity magnitude is based on the refractive index matched shear box experiments of Trewhela et al. (2021). They measured the trajectories of (i) a single large and (ii) a single small intruder for a wide range of shear-rateṡ γ ∈ [0.26, 2.3] and size ratios R ∈ [1.17, 4.17]. Trewhela et al. (2021) made four key observations (a-d below) that allowed them to collapse all their data. (a) Both the large and small intruder data showed a linear dependence of f sl on the shear rateγ . (b) Large intruders have a linear dependence on the size ratio that shuts off when R equals unity, i.e. linear in (R − 1), while (c) small intruders have the same linear dependence at small size ratios, but develop a quadratic dependence on (R − 1) at larger size ratios. Finally, (d) both large and small intruders do not move linearly through the depth of the cell, but describe approximately quadratic curves as they rise up, or percolate down, through it. Since the pressure is linear with depth, this suggests a 1/(C + P) dependence, where the non-dimensional constant C is introduced to prevent a singularity when the pressure is equal to zero. Trewhela et al. (2021) therefore suggested that the segregation velocity A = 0.108, B = 0.3744, C = 0.2712, E = 2.0957, TABLE 2. Non-dimensional constants A, B, C and E in the diffusion (3.7) and segregation laws (3.9) of Trewhela et al. (2021).
magnitude has the form where B, C and E are universal constants. This expression encapsulates the key processes of gravity, shear and pressure, which drive the dominant mechanism for gravity-driven segregation of particles of different sizes and size ratios in shear flows. Moreover, as a consequence of thed 2 dependence, (3.9) automatically gives rise to asymmetric flux functions (Gajjar & Gray 2014;van der Vaart et al. 2015), whose asymmetry is size-ratio dependent (Trewhela et al. 2021). The function (3.9) not only collapses all the single intruder experiments of Trewhela et al. (2021), but it also quantitatively matches the time and spatial evolution of van der Vaart et al. 's (2015) shear box experiments, with a 50 : 50 mix of 4 mm and 8 mm glass beads, using the same values of B, C and E and the generalized diffusion law (3.7). The values of all the non-dimensional parameters are given in table 2. Note that, since the segregation velocity magnitude (3.9) is pressure dependent, but the generalized diffusivity (3.7) is not, Trewhela et al.'s (2021) theory also exhibits the segregation suppression with increased pressure, observed by Golick & Daniels (2009) and Fry et al. (2018). The formula (3.9) cannot be pushed too far, because, for size ratios greater than five, spontaneous percolation is known to occur for low small particle concentrations (Cooke, Bridgwater & Scott 1978), while isolated large intruders may exhibit intermediate or reverse segregation (Thomas 2000;Thomas & D'Ortona 2018).

Numerical method
In order to solve the coupled system of equations the mass and momentum equations (2.1) and (2.2) are written in conservative form where is now the mixture density and ⊗ is the dyadic product. This paper focusses on solving fully coupled bidisperse segregation problems with an evolving free surface using a multiphase approach based on the segregation theory of § 2.2. The method assumes that there are three coexisting phases; large particles, small particles and excess air, which occupy volume fractions ϕ l , ϕ s and ϕ a per unit mixture volume, respectively. In this representation the granular phases are implicitly assumed to retain some air between the grains, so that the overall solids volume fraction in a purely granular state is still Φ as before. Assuming that there is no diffusion of the excess air phase with respect to the particles (i.e. D al = D as = 0) the three conservation laws (2.11) for large particles, small particles and excess air are respectively, where the concentration of grains is defined as When ϕ a = 0, both the large and small particle segregation equations, (4.3) and (4.4), reduce to the bidisperse segregation equation (2.19), and (4.5) is trivially satisfied. As will be demonstrated in § 5, this approach provides a simple and effective way of segregating the large and small particles from one another, while simultaneously expelling unwanted air bubbles and sharpening the free-surface interface. The excess air is assumed to segregate from the grains with constant segregation velocity magnitude f ag along the direction e. The excess air segregation velocity magnitude has no physical significance and the approach should be thought of as a convenient numerical interface sharpening method. The rate is chosen to expel the excess air quickly enough to prevent bubble trapping. For the inclined plane simulations in § § 5 and 6, the direction e is chosen to be the upwards pointing normal to the plane in order to avoid air being segregated through the advancing front. This is not a concern in the rotating drum simulations in § 7 and the direction e is therefore chosen to point in the opposite direction to gravity g.
The system of (4.1)-(4.5) is solved numerically with OpenFOAM assuming that the density and viscosity are given by the local volume fraction weighted averaged values The intrinsic density of the air a * is equal to a constant and the intrinsic densities of the large and small particles are both constant and equal to one another, i.e. l * = s * = Φρ * a * , where the solids volume fraction Φ accounts for the interstitial air that is always present in the granular matrix. The intrinsic viscosity of the air η a * is also assumed to be constant, while the intrinsic viscosity of the grains is calculated from the viscosity (2.3) of the μ(I)-rheology, with the friction μ and inertial number I calculated using the couplings discussed in § 3.2. The parameters used in the simulations in § § 5 and 6 are summarized in tables 1 and 3.
Equations (4.1) and (4.2) are of the form of the incompressible Navier-Stokes equations and the pressure-velocity coupling is solved by the PISO algorithm (Issa 1986). The MUlti-dimensional Limiter for Explicit Solution (MULES) algorithm (Weller 2006), is used to solve the concentration equations (4.4) and (4.5). The first two terms in (4.4) and (4.5) are the same as those in classic multi-phase flow problems, and the inclusion of segregation actually simplifies the problem, as it provides a physical mechanism to counteract the inherent and unwanted numerical diffusion. The numerical treatment of the segregation flux can yield phase fractions outside the interval [0, 1]. Limiting of the respective fluxes (to avoid this discrepancy) is accomplished with the MULES algorithm. The diffusive flux in polydisperse flows is numerically unproblematic and is treated in a similar way to the convective flux, but without a limiter. The coupling of phase fractions with the bulk flow equations for the velocity and pressure is achieved with iterative coupling (Picard iteration) through the respective calculation of local viscosity and density in (4.7a,b).
Numerical diffusion leads to a smearing of the free-surface interface, which has to be suppressed by the numerical scheme. These issues are not limited to the present problem but appear in similar form in many multi-phase problems (e.g. Marschall et al. 2012). In OpenFOAM, this effect is normally corrected with an artificial flux, that compresses the interface (Rusche 2002;Weller 2008). For a general multi-phase mixture the interface sharpening equation for phase fraction ϕ ν is where u νλ is the relative velocity between phases ν and λ. This relative velocity is specifically constructed to be similar in magnitude to the bulk velocity and directed towards regions of higher concentration of phase ν, i.e.
The parameter c νλ is usually chosen to be of order 1 and regulates the amount of counter-gradient transport between phases ν and λ. The counter-gradient flux sharpens the interface, but can lead to results that are outside the range [0, 1] and the MULES algorithm is used again to keep all cell values within this interval.
For the case of a mixture of air and grains, (4.8) and (4.9) reduce to which has the same φ a φ g structure to the air concentration equation (4.5). The key difference, is that (4.5) allows the user to choose the direction e and magnitude f ag of the air segregation, rather than being constrained to the counter-gradient direction. Since many problems of practical interest involve dense granular free-surface flows, with a region of air above them, choosing the direction to segregate the air is not difficult, and completely avoids the unfortunate tendency of interface sharpening methods to create bubbles of air within the body of grains that may remain permanently stuck. The magnitude of the air segregation velocity magnitude may also be chosen to parameterize the typical time scales over which excess air is physically expelled from the body of grains. The polydisperse segregation theory, developed in § 2.2, provides a promising general method of interface sharpening that can be applied to a wide range of problems.
Time stepping is conducted in the ordinary time marching manner. However, special consideration is required due to the spatially varying and high viscosity. In OpenFOAM, each velocity component is solved individually and coupling is achieved explicitly (in a numerically segregated approach). The explicit terms introduce a strict Courant-Friedrichs-Lewy (CFL) criterion which incorporates the local viscosity (Moukalled, Mangani & Darwish 2016). The CFL number is defined as (4.11) and should be limited to a value that is characteristic for the time integration scheme (e.g. 1 for forward Euler). In most multi-phase flows the first term (convection) dominates and the second term (viscosity or diffusion) is neglected. In granular flows with stationary zones, the opposite is the case, since the granular viscosity tends towards infinity in the limit D → 0. To avoid infinitely small time steps, the granular viscosity is therefore limited to a reasonably high value (see e.g. Lagrée, Staron & Popinet 2011;Staron, Lagrée & Popinet 2012), i.e. η = min(η max , η), (4.12) so that η max is the maximum viscosity when the pressure is large and/or the strain rate is small. This is a purely numerical regularization rather than a physically motivated one (see e.g. . The viscous part is still the dominating contribution in the CFL number and granular flow simulations require much smaller time steps than comparable simulations with low-viscosity liquids. Note that computations can be sped up considerably by giving the air phase an artificially high viscosity. This reduces inertial effects in the air, whilst still resulting in a negligible influence of the air on the grains.
The general multi-component segregation-diffusion equations have been implemented into a custom solver based on the OpenFOAM solver multiphaseInterFoam, which makes extensive use of the MULES algorithm provided in the OpenFOAM library. The original solver implements a system of multiple immiscible phases. The system requires an additional diffusion term and replaces the counter gradient transport term with the segregation fluxes. The granular rheology is implemented in a separate library, making use of the respective OpenFOAM programming interface. A similar interface has been created to allow for different expressions for segregation and diffusion coefficients.

Segregation in an uncoupled bulk flow down an inclined plane
The various couplings and feedbacks between segregation and the bulk flow, discussed in § 3, are now explored in more detail. In order to test the numerical method against known steady-state and travelling wave solutions, § 5 examines the one-way coupled model, in which the segregation velocity magnitudes and diffusivities are prescribed, and the bulk flow field is computed with a monodisperse model (as described in § 3.1). The parameters for the bulk flow are summarized in table 1 and are based on the monodisperse glass bead experiments of . The segregation velocity magnitudes and diffusivities are given in table 3 and are chosen to rapidly segregate the air from the grains to produce a sharp free surface, whilst simultaneously allowing a diffuse inversely graded steady-state segregation profile to develop (see e.g. Wiederseiner et al. 2011).

Inflow conditions and boundary conditions
A rectangular Cartesian coordinate system is defined with the x-axis pointing down the slope, which is inclined at ζ = 24 • to the horizontal, and the z-axis being the upward pointing slope normal. The unit vectors in each of these directions are e x and e z , respectively. Numerical simulations are performed on a rectangular grid in the region 0 ≤ x ≤ L x , 0 ≤ z ≤ L z , where L x and L z define the box size. In order to represent an initially empty domain, a Newtonian air phase ν = a is used, which initially fills the box and is stationary, so that ϕ a = 1 and u = 0 everywhere at time t = 0. Granular material, composed of a bidisperse mixture of large ν = l and small ν = s grains, is then injected at the left boundary using Dirichlet conditions on the velocity and on the constituent volume fractions where h is the height of the interface between air and grains at the inflow, and u a and u g = u s = u l are prescribed air and grain velocities. This corresponds to a 50 : 50 mix by volume of large and small grains, with an air phase above. Along the solid base of the chute (z = 0) the no slip and no penetration condition u = 0 is enforced, as well as the no normal flux condition (2.17) for all of the phases. At the outlet wall at x = L x a free outlet condition is applied. This means that there is free outflow (i.e. zero gradient) unless the velocity vector points into the domain (inflow). If inflow is predicted, then the condition switches to Dirichlet and (ϕ a , ϕ s , ϕ l ) = (1, 0, 0) i.e. there is only air inflow and not granular inflow. A similar free-outflow condition applies for the concentration on the top boundary, z = L z . Here the normal velocity has zero gradient, but the pressure is prescribed to be a small constant . Simulations have been performed with p(L z ) = 10 −3 N m −2 and 10 −6 N m −2 and are insensitive to this change.
5.2. Steady uniform bulk flow velocity As this becomes an effectively monodisperse problem for the bulk flow u and pressure p, fully developed steady uniform flow should correspond to the Bagnold flow solution (see e.g. Silbert et al. 2001;GDR MiDi 2004;Gray & Edwards 2014;Barker et al. 2015). Assuming a flow of thickness h, the exact solution to the μ(I)-rheology implies that the pressure is lithostatic the downslope velocity is given by the Bagnold profile 4) and the inertial number I is equal to the constant For the partially regularized form of the friction law (2.8) used in this paper, it follows, that for μ ∞ > 0 and I > I N 1 , the inertial number is equal to The granular inflow velocity is therefore set to u g = u Bagnold e x . The velocity in the air phase above is set to the Newtonian flow solution u a = u p (z) e x , where the Poiseuille profile is This implicitly assumes no slip at the lower free-surface interface with the moving grains, i.e. u p (h) = u Bagnold (h).

Comparison between the different methods of interface tracking
For this simple case, it is instructive to compare the alternative interface sharpening techniques that were discussed in § 4. As shown in figure 2(a), when there is no interface sharpening, numerical diffusion leads to a very wide diffuse layer between the air and the grains, rather than a sharp free surface. In addition, a large vortex of dilute granular material is thrown into the air at the front and a thin layer of air is trapped next to the basal solid wall. This trapping of air next to the boundary is a serious problem, because it prevents direct contact of the grains with the lower boundary and consequently affects the effective friction experienced by the grains as they flow downslope. In reality, any air that is trapped adjacent to the lower wall is free to percolate up through the pore space between the particles and escapes. This unphysical air trapping is also observed in the simulation with active counter-gradient transport as shown in figure 2(b). Although the free surface is much sharper than before, there is a tendency for the trapped air to form bubbles. This effect is especially strong in high viscosity flows because the bubbles become stuck and are unable to escape. The results, both with and without interface sharpening, are also found to be sensitive to the numerical mesh and time step used in the calculation. Figure 2(c) shows the new method of tracking the interface using (4.4) and (4.5) assuming that trapped air is segregated upwards, i.e. e = e z . The segregation velocity magnitude and diffusion coefficients (see table 3) are chosen to give diffuse segregation inside the granular mixture, but also to generate a sharp interface between the granular phases and the air above. It is clear from figure 2(c) that with this method there is no trapped air next to the basal boundary, the free-surface interface is sharp and there is no vortex shedding at the flow front. Moreover, the results are grid converged. The new method of treating the free surface is therefore very promising, and provides a simple way of parameterizing the physics that is actually taking place.

Numerical simulation of the bulk flow and the segregation
Armed with this improved and reliable method of interface capture, the full transient evolution of the travelling front can be explored. Figure 3 shows the results of a calculation performed in a long aspect ratio domain with dimensions (L x , L z ) = (0.62, 6.2 × 10 −3 ) m i.e. 100 : 1. As the front progresses into the domain, there is dynamic evolution of both the front shape and the distribution of the granular phases. In particular, a steadily travelling front forms with a well-defined shape (Pouliquen 1999a;Gray & Ancey 2009;Saingier, Deboeuf & Lagrée 2016). Behind the advancing front, the initially evenly mixed concentration of large and small grains is swept downstream from the inflow and is gradually eroded by a growing layer of large particles at the surface and a growing layer of fines adjacent to the base of the flow. By 20 cm downstream the homogeneously mixed region completely disappears and further downstream there is a thin layer with high concentrations of large grains at the surface and a thicker layer with high concentrations of fine grains at the base. This is known as an inversely graded particle-size distribution. The difference in thickness is due to the large particles being concentrated in the faster moving region of the flow, so a much thinner layer can transport the same mass flux as the thick, slow moving layer beneath, which contains high concentrations of fines. An immediate consequence of the large particles being segregated into the faster moving near surface layers is that they are preferentially transported to the flow front, as shown in figure 3(b-d). As large grains reach the front, they are over-run, but can rise back towards the surface again by particle segregation, to form a recirculating frontal cell of large particles that grows in size with increasing time (Pierson 1986;Pouliquen et al. 1997;Iverson & Vallance 2001;Gray & Kokelaar 2010b, a;Johnson et al. 2012;Woodhouse et al. 2012;Baker et al. 2016b;Denissen et al. 2019). The large particle rich flow front and the inversely graded body of the flow are connected by what is known as a breaking size segregation wave (Thornton & Gray 2008;Johnson et al. 2012;Gajjar et al. 2016). This travels steadily downslope, but at a slower speed than the front. It is this wave that segregates the large slow moving particles, close to the base of the flow, up into faster moving regions allowing them to be recirculated, and conversely, allows fast moving small grains to percolate down into slower moving basal layers. The breaking wave shown here includes the effects of diffusion as well as segregation for the first time. Eventually the flow front and the breaking size segregation wave propagate out of the domain, to leave the approximately steady uniform flow shown in figure 3(e). For comparison, Gray & Thornton's (2005) concentration shock solution (see appendix A) is also plotted in figure 3(e) using the Bagnold velocity profile (5.4). For an inflow small particle concentration ϕ s 0 = 0.5 this accurately predicts the position of the centre of the final steady-state height of the inversely graded layer, with the large particles occupying a thinner faster moving region than the fines. However, the solution neglects diffusion in both the downslope and slope normal directions, and only resolves the segregation flux in the slope normal direction, so it does not capture the precise point at which the solution reaches steady state. Figure 4(a) shows excellent agreement between the computed two-dimensional steady uniform flow solution for the downslope velocity u and the Bagnold velocity profile (5.4). The only slight difference occurs near the free surface, where the weight of the column of air above produces the largest relative change in the pressure within the granular material. With the μ(I)-rheology, this changes the balances in the inertial number and hence the computed velocity profile. For steady uniform flows, Gray & Chugunov (2006) derived an exact solution for the small particle concentration, assuming that the segregation and diffusion rates were constant. This solution takes the form

Comparison with steady uniform solutions for the bulk flow and the segregation
where A GC is a constant and Pe is the Péclet number for segregation. Note that in this solution the z-coordinate has been non-dimensionalized using the scaling z = hẑ, where h is the slope normal flow depth. In terms of the dimensional segregation and diffusion rates, given in table 3, the Péclet number is defined as where the factor cos ζ arises from the fact that the segregation is inclined at an angle ζ to the slope normal z-axis, i.e. e z · g/|g| = − cos ζ . The constant A GC alters the position of the transition between large and small particles in the solution. If the depth-averaged concentration is equal to   Gray & Chugunov (2006). The parameters are summarized in tables 1 and 3.

Comparison of the frontal shape with depth-averaged solutions
The basal friction law of Pouliquen (1999b) predates the full tensorial μ(I)-rheology and was designed to model the frictional source term in the shallow avalanche equations of Savage & Hutter (1989) on chutes with rough bases. The fully developed numerical front solution, shown in figure 3, is indeed very shallow, so it is appropriate to compare it with solutions of these reduced equations. The depth-averaged theory provides a very simple means of predicting the shape of a steadily travelling granular flow front (Pouliquen 1999a;Gray & Ancey 2009;Saingier et al. 2016). In a frame ξ = x − u F t moving with the front speed u F the steady-state depth-averaged mass and momentum balances are where h is the avalanche thickness, and the depth-averaged velocityū, the depth-average of the velocity squared u 2 and the shape factor χ are defined as respectively. Many theories assume that the shape factor χ = 1, which corresponds to plug flow, and which dramatically simplifies the characteristic structure of this hyperbolic system of equations. For the Bagnold velocity profile (5.4), the shape factor χ = 5/4. Saingier et al. (2016) showed that with Pouliquen & Forterre's (2002) effective basal friction law this led to the formation of a thin precursor layer ahead of the main front that extended to infinity, which is unphysical. The depth-averaged mass balance equation (5.12) can be integrated directly, subject to the condition that the thickness is zero at the flow front, to show that for non-trivial solutions the depth-averaged velocity is equal to the front speed, i.e. u = u F , (5.15) everywhere in the flow. Far upstream the flow is steady and uniform. The front speed can therefore be determined by integrating the Bagnold solution (5.4) through the avalanche depth to show that where h ∞ andū ∞ are the steady uniform thickness and downslope velocity far upstream. Expanding (5.13), dividing through by hg cos ζ and using (5.15) yields an ordinary differential equation (ODE) for the flow thickness Since, the depth-averaged velocity is the same as the front velocity (5.15) everywhere in the flow, (5.16) and (5.19) can be equated to determine the inertial number at a general position ξ . Substituting this expression into the high-I branch of the full μ(I) curve (2.8) gives the regularized depth-averaged basal friction The significance of this expression is made clear by taking the limit as h → 0. Unlike for the previous expression for μ, in which μ ∞ = 0, the friction now tends to infinity for vanishingly thin layers. This means that the ODE (5.17) naturally predicts an infinite slope and therefore the front always pins to the boundary and this system is guaranteed to preclude infinite precursor layers. The front shape predicted by this newly derived regularized depth-averaged formulation is compared with the full two-dimensional numerics in figure 5. In order to guarantee that the full solution does indeed correspond to a steady travelling front, the simulation is continued from t = 4 s in a moving frame. This change is applied simply by shifting all the velocities and the boundary conditions by the depth-averaged velocity (5.16) i.e. u new = u(t = 4 s) −ūe x everywhere. The following analysis applies to the long-time solution in this moving frame, which is found to be numerically invariant of time after another ∼5 s of simulation. Upstream of the front (for low values of ξ ) the flow is almost uniform, so the Bagnold solution, which has a shape factor χ = 5/4, is observed as expected. However, closer to the flow front the assumption of uniformity breaks down and the two solutions differ. As shown in figure 5, the front computed with the multi-phase approach lies between the depth-averaged solution with χ = 5/4 and that with χ = 1, which corresponds to pure plug flow, where u no longer depends on z. This comparison therefore highlights the expected discrepancies between full two-dimensional theories and depth-averaged equivalents when the dynamics varies in a non-shallow manner.

The two-dimensional internal flow fields in the moving frame
Given that the two-dimensional transient flow front has developed into a steady travelling state, the detailed flow fields inside the granular material are of particular interest. These are plotted in figure 6. Figure 6(a) shows the downstream velocity, shifted back to the laboratory frame by addingūe x , which is monotonically increasing in z for all x in a similar manner to the Bagnold velocity profile. Only at the tip of the front is the vertical velocity non-zero (figure 6b) and there is a downwards motion. As these two velocity components define a steady travelling front, the streamlines which result from them coincide with the particle paths. However, these trajectories, which are plotted in figure 6(c), only correspond to the paths of monodisperse particles. The large and small particle trajectories, which couple to these flow fields, but not vice versa, are not steady in this frame, or any frame of reference as the large particle recirculation region at the head is forever growing in size. Just like the similarity to the Bagnold velocity solution, the pressure field in figure 6(d) is close to the lithostatic profile (5.3) except that the flow thickness is not constant. Similarly, the inertial number (figure 6e) takes its steady uniform value upstream, but gets larger as the front is approached, as predicted by (5.20). It should be noted that any potential issues of ill posedness at high inertial numbers, close to the very tip of the flow, are suppressed by the maximum viscosity cutoff (4.12) in the numerical method.

Segregation mobility feedback on the bulk flow
The one-way coupled simulations in § 5 demonstrate the effectiveness of the numerical method developed in § 4, and also show qualitatively how large and small particles are advected, segregated and diffused within the bulk flow field. To produce quantitative results, it is necessary to couple the evolving particle-size distribution to the bulk flow dynamics, as discussed in § 3.2. There are essentially two ways of producing frictional feedback; namely (i) indirectly through the evolving average local grain size, which changes the inertial number and hence the friction, and (ii) directly through the modification of the frictional parameters associated with each of the species. Both couplings are investigated in this section, and the results of the inertial number coupling are compared directly with the DEM simulations of Tripathi & Khakhar (2011).

Steady uniform flow down an inclined plane with segregation mobility feedback
Consider once again a steady uniform flow down an inclined plane, but this time incorporating feedback of the steady-state concentration distribution. If the segregation and diffusion rates are constant, then the volume fractions ϕ ν = ϕ ν (z) can be solved for with the polydisperse theory in § 2.2, completely independently of the bulk flow. These concentrations will therefore be assumed to be known in what follows. The normal component of the momentum balance then implies that the pressure is lithostatic (5.3). The only difference to the classical Bagnold solution (Silbert et al. 2001;GDR MiDi 2004;Gray & Edwards 2014) is that, with the volume fraction weighted friction (3.3), the downslope momentum balance reduces to ∀ν ϕ ν μ ν (I) = tan ζ, (6.1) where μ ν is the friction law for constituent ν. For the purposes of illustration, let us assume that each phase satisfies the classical μ(I) friction law, which is of the form where I 0 is assumed to be the same for all the phases. Substituting (6.2) into (6.1) and solving for the inertial number, it follows that whereμ s andμ d are now the volume fraction weighted averages that are depth dependent Importantly, (6.3) shows that, if there are frictional differences between the particles, then the inertial number is dependent on the normal coordinate z rather than being equal to the constant I ζ defined in (5.5). Using the definition of the generalized inertial number for polydisperse systems (3.2) and assuming steady uniform flow, it follows that the ODE for the velocity profile is whered is the local average particle size, which is also depth dependent (6.6) This averaged particle-size dependence is important, because even if the particles have the same shape and the same effective frictional properties, the velocity profile will no longer be the classical Bagnold solution (5.4), but will depend on the local changes in particle size. Figure 7 shows a specific example of the qualitative types of solution that are generated for a bidisperse mixture of large and small particles. The solutions assume Gray & Chugunov's (2006) exact solution for the concentration profile (5.8)-(5.11) using the same constant segregation velocity magnitude f sl , constant diffusivity D sl , flow depth h as in table 3, as well as the same slope angle ζ = 24 • . The only difference is that the depth-averaged concentration ϕ s is chosen to be equal to 50 % in order to produce flowing velocity for a bidisperse mixture of large and small particles (black lines) on a slope inclined at ζ = 24 • to the horizontal. The solutions assume a small particle concentration profile given by Gray & Chugunov's (2006) exact solution in (5.8)-(5.11), with ϕ s = 0.5 and using the parameters in table 3. Here, all bulk flow parameters are identical to those in table 1 except that the large particles have μ l s = 1.2μ s and μ ν ∞ = 0 for both phases. The dashed lines indicate uniform concentration solutions with red corresponding to pure large, blue corresponding to pure small particles and green being the solution for a mixture with ϕ s = 0.5 everywhere.
layers of large and small particles that are the same depth. For consistency with the assumed friction law (6.2), μ ν ∞ = 0 for both the large and small particles. All the other parameters are the same for both species, and identical to those given in table 1, except that μ l s = 1.2μ s . This small change is sufficient to make the inertial number (6.3) depth dependent, as shown in figure 7(a). The increase in μ l s for the large particles decreases the inertial number in the near surface regions, where the large particles are located. Integrating the ODE (6.5) through the flow depth, subject to the no slip condition at the base, gives the velocity profile in figure 7(b). The solution lies between the velocity profiles for pure large and for pure small particles, and closely follows the small particle velocity profile in the lower part of the flow, where the small particles are concentrated. In the upper part of the flow it rapidly transitions onto a curve that is parallel to that of the pure large particles, but they attain a much higher speed than if there were no small particles in the flow. Or indeed, if the particles were evenly mixed throughout the column with ϕ s = 1/2 everywhere. The small particles therefore provide an important lubricating mechanism that can significantly increase flow speeds and the overall run-out (Kokelaar et al. 2014).

Formation of a large rich bulbous flow front on an inclined plane
Given the steady solution in § 6.1, it is also interesting to consider the transient behaviour of a granular flow front when the large particles are more frictional than the fines. Analogously to the DEM study of Denissen et al. (2019), the solution detailed in figure 7 is used as the boundary condition at the inlet wall x = 0, so that material entering the domain is already stratified and well developed. All parameters are the same as those in § 6.1. As shown in figure 8(a), the two-dimensional transient dynamics generates a bulbous head of large particles in front of an approximately uniform thickness upstream flow. This bulging of the surface differs from the monotonically decreasing free-surface shape, observed when there is no feedback of the segregation on the bulk, as shown in § 5 and figures 3-6. The fundamental cause of this effect is that pure regions of large particles are much less FIGURE 8. Contour plots of (a) the concentration of small particles and (b) the base 10 logarithm of the inertial number at t = 5.2 s for a flow in which the large particles are more frictional than the fines. Here, as in figure 7, the parameters for each species are identical to those in table 1 except that μ l s * = 1.2μ s and μ ν ∞ = 0 for both species. The inflow concentration is assumed to be a steady uniform solution (5.8) of the segregation equations assuming the parameters in table 3 and with ϕ s = 0.5. A movie of the full dynamics is available in the online supplementary movie 2. mobile than the inversely graded flows behind, which are lubricated by the fine particles at the base. The preferential transport of large particles to the front, where they recirculate and accumulate (by a combination of the bulk flow field and particle segregation) causes the front to grow in size and become increasingly resistive. This causes it to bulge upwards until it (i) stops and blocks the flow, (ii) permanently deposits some of the large grains on the substrate and flows over them (Gray & Ancey 2009), (iii) pushes some of the large particles to the side to form static levees (Pierson 1986;Pouliquen et al. 1997;Pouliquen & Vallance 1999;Iverson & Vallance 2001;Woodhouse et al. 2012;Kokelaar et al. 2014;Baker, Barker & Gray 2016a) or (iv) becomes sufficiently thick that a flow of large particles can form that moves slightly faster than the thinner upstream inversely graded layer behind, to accommodate the continued supply of large particles to the front (Denissen et al. 2019).
This problem therefore has a very strong two-way coupling between the bulk flow and the segregation. As shown in figure 8(b), the inertial number in the flow front provides a clear demonstration of this coupling. Upstream of the head, where the flow is uniform, I approximately matches the two-layer solution from figure 7(a) and close to the flow head the fields are reminiscent of the monodisperse case detailed in figure 6(e). A diffuse breaking size segregation wave (Thornton & Gray 2008;Johnson et al. 2012;Gajjar et al. 2016) allows the two regions to connect to one another. It is located at x 450 mm and is clearly evident in both the small particle concentration distribution as well as in the inertial number distribution. This is therefore the first fully coupled breaking size segregation wave to be computed. Tripathi & Khakhar (2011) To provide a quantitative comparison for the steady-state behaviour, the theory is now compared with the bidisperse DEM simulations of Tripathi & Khakhar (2011) Behringer's (2004) diffusivity to bidisperse systems (rather than prescribed rates). The results shown in Tripathi & Khakhar's (2011) figure 9 correspond to flow down a plane inclined at an angle ζ = 25 • , in which the large particle diameter is one and a half times the small grain diameter, i.e. d l = 1.5d s . The results are presented in non-dimensional form, where the length, time and velocity scalings (6.7a-d) are based on the small particle diameter d s and gravity g. The layer depth h is assumed to be 30d s . The simulations are performed in a three-dimensional cell that is periodic in the down and cross-slope directions, and has a fixed bed that is made rough with particles of diameter 1.2d s . The down and cross-slope dimensions are 20d s × 20d s . Figure 9 shows Tripathi & Khakhar's (2011) computed small particle concentration and downslope velocity for five different depth-averaged concentrations, ranging from pure small to pure large. For comparison, the bidisperse small particle concentration equation (2.19) is solved at steady state, assuming the functional forms (3.7) and (3.9) for the segregation velocity magnitude and diffusivity, i.e.

Comparison with the steady-state DEM solutions of
where A, B, C and E are non-dimensional constants andγ has been replaced by its equivalent strain-rate invariant, i.e. 2 D . Assuming that the downslope velocity and the small particle concentration are purely functions of the slope normal coordinate z, (2.19) can be integrated once with respect to z. Applying the no flux boundary condition (2.17) at the surface and/or base of the flow, the D d 2 dependence in the segregation and diffusive terms cancels out. As a result the final steady-state ODE for the concentration is independent of the shear rate, uncoupling it from the downslope momentum balance. The non-dimensional parameter C is primarily introduced to remove the pressure singularity at the free surface, and its measured value of C = 0.2712 makes very little difference to the shape of the concentration profile (Trewhela et al. 2021). If instead C is assumed to be zero, and the pressure is lithostatic (5.3), then the intrinsic grain density ρ * , gravity g and the slope angle ζ cancel out, leaving a non-dimensional steady-state ODE for the concentration which is dependent purely on the grain-size ratio R. This is separable, and can be integrated (Trewhela et al. 2021) to give the exact solution where K is a constant of integration and the coefficients λ 1 , λ 2 and λ 3 are 1)) .
(6.11a-c) For a given depth-averaged small particle concentration ϕ s the constant of integration K can be determined iteratively. Figure 9(a) shows the comparison between these exact There is also very good agreement at ϕ s = 50 % and 70 % using exactly the same non-dimensional constants A, B and E determined experimentally by Trewhela et al. (2021) and summarized in table 2. The theory therefore matches the stronger segregation at the surface than the base of the flow, as well as the gradient of the concentration profiles, without the need for any fitting parameters. This is strong evidence that Trewhela et al.'s (2021) theory captures the essence of the segregation process. It also contrasts with Gray & Chugunov's (2006) solution, where the ratio of segregation to diffusion is uniform with depth. The agreement between Trewhela et al.'s (2021) theory and Tripathi & Khakhar's (2011) DEM simulations is not as good at ϕ s = 30 %. The DEM results at ϕ s = 30 % look slightly odd, with a layer of almost pure small particles at the base of the cell and a much more diffuse profile higher up. It is therefore possible that, in this particular case, Tripathi & Khakhar's (2011) DEM simulations have not fully reached steady state. Figure 17 of Tripathi & Khakhar (2011) suggests that the friction in both their monodisperse and bidisperse systems was closely approximated by the classical μ(I) law (2.7), using the generalized inertial number (3.2) with a local average grain sized. To leading order, therefore, the macroscopic friction coefficients μ s and μ d , as well at the non-dimensional constant I 0 are the same for the large and small particles. Tripathi & Khakhar (2011) suggested that a good overall fit to the data was provided by μ s = tan(20.16 • ), μ d = tan(37.65 • ) and I 0 = 0.434. These values are, however, not good for the particular set of simulations shown in Tripathi & Khakhar's (2011) figure 9, and reproduced here in figure 9. To select a better fit, the Bagnold solution (5.4) has been non-dimensionalized using the scalings (6.7a-d) to givê and then fitted to the small particle velocity DEM data using a least squares fit. This determines the value of I ζ , which for the classical friction law (2.7) of Jop et al. (2006) is defined as The values of μ s , μ d and I 0 can therefore be modified, while still fitting the data, provided the same value of I ζ is obtained. There are an infinite number of combinations that will do this. This paper therefore assumes that the values of μ s and μ d are the same as Tripathi & Khakhar (2011), but that I 0 = 0.5106. To solve for the velocity profiles at other concentrations it is necessary solve the ODE (6.5) in non-dimensional form, i.e. (6.14) subject to a no-slip boundary condition at the base. In the ODE (6.14) the small particle concentration ϕ s = 1 − ϕ l is given by Trewhela et al.'s (2021) exact solution (6.10). The solutions are shown in figure 9(b). The 100 % small particle solution agrees extremely well with the Bagnold solution (6.12), which is not too surprising as the parameters μ s , μ d and I 0 have been chosen specifically to match this curve. The monodisperse large particle solution also has a Bagnold-like velocity profile, but the magnitude of the velocities are slightly underpredicted. In principle, the monodisperse small particle solution should be a factor R larger than the large particle solution. The fact that they are not, is an indication of either (i) the level of noise in Tripathi & Khakhar's (2011) system, or (ii) the basal roughness not scaling with the size of flowing particles (i.e. non-local effects).
Since the deviations of the intermediate solutions from the DEM data are of a similar order of magnitude, it is probably unwise to read too much into the precise comparisons. The predicted velocity profiles at intermediate concentrations monotonically decrease in magnitude as the small particle content decreases. This is broadly speaking what the DEM data show; however, the case ϕ s = 30 % the DEM solution is again anomalous, being below the case of pure large particles (Tripathi & Khakhar 2011). To be sure that this is real behaviour, rather than an anomaly, more precise DEM solutions are required that average over significantly more than the approximately 4000-12 000 particles used by Tripathi & Khakhar (2011).

Fully coupled rotating drum simulations
Particle segregation in non-circular rotating drums provides an ideal test case for the two-way coupled model, as the particles strongly segregate and diffuse in the near surface liquid-like avalanche, but not in the solid-like rotating body beneath. Computing the bulk flow field in a rotating drum is still a significant challenge. Indeed, recent segregation simulations in a circular rotating drum (Schlick et al. 2015) have prescribed the steady-state incompressible bulk velocity field based on fits to DEM simulations; a process that inherently relies on the steady-state nature and simple geometry of the circular drum. More complex models that do use a continuum approach to calculate the bulk flow in a circular drum (e.g. Liu, Gonzalez & Wassgren 2018, do so with rate-independent elasto-plastic constitutive laws which are prone to ill posedness (Schaeffer 1987). This paper goes considerably further, by using the partially regularized μ(I)-rheology and the recent segregation model of Trewhela et al. (2021) to simultaneously compute the fully two-way coupled bulk flow, segregation and diffusion in a square rotating drum.

Modelling the bulk flow, segregation and diffusion in a square rotating drum
It is useful to have two coordinate systems to simulate the flow in the drum. The first is a rectangular Cartesian coordinate system Oxz that is fixed and centred at the axis of rotation of the drum, which lies at the centre of the square. The z axis is aligned with the gravitational acceleration g, but points upwards in the opposite sense. A second coordinate system OXZ is inclined at an angle θ to Oxz and rotates with the drum. The axes are aligned with the drum walls, so the drum lies in the region −L ≤ X ≤ L, −L ≤ Z ≤ L, where 2L is the length of the walls. Initially, the OXZ axes coincide with Oxz and the concentrations of excess air, small particles and large particles are where L = 0.1 m and H = 0.04 m, implying a 70 % fill fraction with a 50 : 50 mix of large and small particles of diameters d l = 2 mm and d s = 1 mm, respectively. A fill fraction above 50 % is chosen so that an undisturbed core forms in the centre of the drum, consisting of material which never enters the avalanche (see e.g. Gray & Hutter 1997;Gray 2001). All the material is initially assumed to be in solid-body rotation where Ω is the rotation rate, the radial coordinate r = √ x 2 + z 2 and θ is the azimuthal unit vector. A constant rotation rate of Ω = −π/5 rad s −1 is specified, with the negative sign denoting clockwise rotation. This corresponds to one full revolution every 10 s, placing the flow at the upper end of rolling flow (Henein, Brimacombe & Watkinson 1983;     partially regularized friction law with parameters that match the steady-state DEM simulations of Tripathi & Khakhar (2011) shown in figure 9. The value of μ ∞ is chosen to ensure that the equations remain well-posed up to a maximum inertial number I max = 16.20, while I 1 is the minimum well-posed inertial number in the unregularized law. To handle the evolving free surface, the excess air segregates with a constant rate f al = f as from the large and small particles and does not diffuse with them. The air phase is assumed to segregate upwards in the direction of the unit vector k along the z-axis, which this time is aligned with gravity. Ding et al. 2002;Yang et al. 2008). This is also known as the continuous, or the continuously avalanching, regime (Rajchenbach 1990;Gray 2001) as a quasi-steady-state avalanche forms, with continuous erosion and deposition occurring with the solid rotating body of grains beneath. The frictional parameters are the same for the large and small particles, and the momentum coupling enters via the evolving local average particle sizē d in the generalized inertial number (3.2). The values of μ s , μ d and I 0 are the same as those used to fit the DEM simulations of Tripathi & Khakhar (2011) (see figure 9), and the theory is partially regularized ) by introducing a creep state at low inertial numbers and a linear friction regime at high inertial numbers. The values of all the frictional parameters are summarized in table 4, together with the particle sizes, and the air segregation and diffusion rates. The segregation of the particles is performed using the same non-dimensional constants A, B, C and E suggested by Trewhela et al. (2021) and summarized in table 2.
The velocity fieldû in the rotating frame is related to the velocity u in the fixed coordinate system by the relationû = u − Ωrθ . 7.2. Bulk velocity, pressure and inertial number in the square drum Initially the free surface is flat and the material is in solid-body rotation, see (7.1) and (7.2). The entire body of grains is therefore quasi-static in the moving frame, and the high viscosity cut-off (4.12) remains active until the free-surface inclination nears the static angle of friction ζ s = tan −1 (μ s ), (7.5) at which point the material near the free surface fails and avalanches downslope. After the initial failure and slump, a continuously avalanching regime rapidly establishes itself, as shown in figure 10 and the supplementary movies. There is a rapid liquid-like avalanche close to the free surface and a solid-like quasi-static region beneath ( figure 10a,b). The angle of the free surface remains close to the static value ζ s , but the position of the free surface subtly rises and falls as the finite volume of grains is incorporated into the constantly changing intersections with the shape of the drum during each quarter turn.
The flow therefore has a quasi-periodic pulsing behaviour, with peak surface velocities (at the centre and surface of the flow) changing in time, e.g. the peak free-surface velocity is faster at t = 80 s than at t = 78.75 s in figure 10(a,b). As shown in § 6.3, the variations in velocity with the flow composition are subtle, and do not provide a strong feedback on the bulk flow at this fill height. However, the experiments of Zuriguel et al. (2006) imply that for fill levels close to 50 % these subtle composition-dependent velocity changes can cause the formation of petal-like structures in the deposit, and so can be very important.
As the flow pulses, the surface avalanche becomes deeper directly beneath the region where the peak velocities are attained. The avalanche depth also changes along its length, reaching a peak near its centre. Typically the main flow is confined to a layer with a maximum thickness of 1.6 cm, which gives the rapid free-surface flow a shallow aspect ratio, consistent with the assumptions underpinning theoretical models for avalanches (e.g. Savage & Hutter 1989;Gray 2001;Gray & Edwards 2014). Close to the free surface the pressure is approximately lithostatic and is aligned with the free surface (figure 10c,d) as one might expect. However, lower down the pressure rises to much higher values and pulses as the overall volume of grains redistributes itself in the changing geometry of the walls that confine it. The base ten logarithm of the inertial number is shown in figure 10(e, f ) and also identifies the near surface region where the failure occurs. The flow is in a creep state for I < I 1 = 0.01886. This region lies significantly higher in the flow than the Newtonian viscous region, which is activated by the numerical regularization (4.12) at high viscosities. The dominant rheology in the simulation is therefore the granular rheology, which involves regions of both creep and dynamic motion.
For fill levels above approximately 55 % a solid core develops (Mounty 2007), within which particles are simply rotated around with the drum and undergo a small amount of creep when they are closer to the free surface. The remaining grains pass through both the solid-like and fluid-like regions. These particles are rotated around with the drum in the solid-body region, until they approach the near surface layers when they begin to creep downslope. As an individual particle is rotated further towards the free surface the creep becomes progressively stronger, until finally it avalanches downslope in the liquid-like surface avalanche. As more particles are entrained into this avalanche it becomes deeper and flows faster, so peak velocities are reached midway down the slope, after which particles are deposited from the base of the flow, and the avalanche thins and slows. An individual particle therefore accelerates downslope in the first half of the avalanche and decelerates after the midway point, before being deposited into the slowly creeping body of rotating grains beneath. Unlike a circular drum, where monodisperse particles form closed streamlines, the changing geometry of the confining walls adds considerable complexity to the problem. This is because the underlying particle trajectories become chaotic even for monodisperse flows (Hill et al. 1999;Ottino & Khakhar 2000). Particle-size segregation, however, introduces a strong organizational influence on the resulting patterns that form. 7.3. The particle-size distribution in the square drum The dynamics of the mixing and segregation process is shown in figure 11 and the supplementary online movie. As the drum rotates up to the static angle of friction there is no shear and hence no segregation or diffusion. However, as soon as the initial failure occurs, and the avalanche flows downslope, the particles begin to segregate with the large particles rising to the free surface and the small particles percolating downwards. The linear shear rate dependence in (6.8a,b) ensures that both the segregation and diffusion are confined to the thin avalanching layer close to the free surface, with the additional pressure dependence ensuring that segregation shuts off more rapidly than the diffusion with increasing depth (Golick & Daniels 2009;Fry et al. 2018;Trewhela et al. 2021). This effect is compounded by the fact that particles near the free surface travel the longest distance through the liquid-like avalanche, while those that are entrained at lower levels may move only a short distance before they are deposited back into the solid-like rotating body of grains beneath. As a result, after the first full rotation of the drum (at t = 10 s) the clearest segregation can be seen in the large particles that are able to rise to the surface and collect at the top of the flow before being deposited near the drum wall. The rest of the grains remain quite well mixed.
After approximately 3/4 of a drum revolution all the material that is able to pass through the surface avalanche has done so, and the material that is subsequently entrained is no longer homogenously mixed. The large particles that were deposited next to the drum wall are rotated around and re-entrained into the avalanche right at the back of the flow, where the avalanche is thinnest. There is therefore no need for them to segregate towards the surface again, as they are naturally re-entrained on trajectories that pass through the surface layers of the avalanche. The particles that lie closer to the drum core are naturally entrained onto paths that take them through lower regions of the avalanche. They therefore get another chance to segregate again each time they pass through the avalanche. This process can be seen slowly sharpening the segregation in the successive panels of figure 11. With increasing time, the large particle region adjacent to the drum wall thickens up and regions with high concentrations of small particles start to emerge, as the avalanche at the surface becomes progressively more inversely graded. Complete separation of the large and small grains does not occur, however, because of the diffusive remixing process in (2.19).
The combination of particle segregation and the rising and falling of the free-surface height as the drum rotates leads to the spontaneous formation of three lobes with high concentrations of small particles, that are oriented towards the corners of the drum. These lobes are interesting because they propagate around the drum faster than the drum rotates, with a period of approximately 7.5 s. The lobes are in qualitatively very good agreement with the experiments of Hill et al. (1999), Ottino & Khakhar (2000) and Mounty (2007), as shown in figure 12. The simulations also predict the formation of a central core within which the concentration is almost unchanged from its initial value. This core forms a shape that is almost square and lies at an angle of 45 • to the square drum walls, which is also qualitatively in agreement with the experiments of Hill et al. (1999), Ottino & Khakhar (2000) and Mounty (2007). However, the simulated central core is about half the diameter of that in experiment. The reason for this is that the surface avalanche is much deeper in the simulations than in the experiments. This is not necessarily a deficiency of the model. The experiments are performed in drums with a narrow gap between the sidewalls; as a result the avalanche is thinner and faster than in a drum with a wide cross-section (Jop et al. 2005) and hence the central core is larger.
FIGURE 11. Fully coupled simulation of a bidisperse granular mixture in a square rotating drum using the parameters in tables 2 and 4, where the drum walls are of length 0.2 m: (a) plots ϕ s at 10 s and (b-l) correspond to a further 6.25 s of rotation, or 5/8 of a full revolution, up to 78.75 s. A movie 6 is available in the online supplementary movies.
One of the important consequences of the D d 2 dependence in the segregation velocity magnitude and diffusivity in (6.8a,b) is that the time scale for segregation and diffusion to occur is proportional to h 2 /(γd 2 ). It therefore takes longer to segregate in a deeper avalanche, or with smaller average grain sizes. This makes direct comparison with the experiments difficult, as the depth and velocity of the surface avalanche is strongly influenced by sidewall friction. In principle, it is very easy to include the effect of sidewall friction in the simulations. However, the numerical method requires sufficient grid points to be located in the surface avalanche. For a regular grid this requires higher resolution throughout the drum, which dramatically increases the time necessary to produce grid converged results, and so this is not done here. Instead the grain sizes are made larger in order to get the pattern to form in the rotating drum simulations in a comparable time scale to that in the experiments.
The lobes do not appear to reach a quasi-periodic steady state, but have small protuberances that continue to evolve when tracking a particular lobe. This is also in accordance with experimental observation, as a non-homogeneous initial distribution of particles can lead to lobes of different sizes, which appear to persist indefinitely. Figure 11 also shows that over long periods of time the central core contracts towards the origin. This is due to the slow creep that occurs as the material in the solid-body region is rotated through the near surface creeping zone, where both segregation and diffusion can act over very long time scales. This creep can be minimized either by (i) introducing sidewall friction or (ii) by using constitutive equations with a static yield stress. However, both of these require additional physics to be included in the model.

Grid convergence
A grid convergence study was carried out for four different mesh refinements. Figure 13 shows the evolution over time of the integral I (t), defined as FIGURE 13. Integral I (t) of the small particle concentration difference between the current state and the state of one full revolution earlier, as a function of time and for four different mesh resolutions N to demonstrate grid convergence of the numerical solution.
difference between the current state and the state a full rotation period earlier. The maximal value of I (t) is unity. Numerical diffusion means that overall segregation is weaker when the flow is under resolved, and, at N = 50, regions of high concentration fail to coalesce. The increasing proximity of the curves with increasing grid resolution, and particularly the closeness of I (t) for N = 150 and N = 200, demonstrate numerical grid convergence. This relatively high number of grid points is required to properly resolve the shallow avalanche at the free surface of the flow where the segregation and diffusion predominantly occur. As noted earlier, higher grid resolutions will be necessary to resolve the thinner surface avalanches that develop in experiments with sidewall friction (Hill et al. 1999;Ottino & Khakhar 2000;Jop et al. 2005;Mounty 2007). The integral I (t) implies that the square drum does not approach a periodic quasi-steady solution, but settles down, after about four complete revolutions, to a state where dI (t)/dt is small and I (t) is non-zero. This represents a fully segregated mixture with time-dependent perturbations that propagate around the system.

Conclusions and discussion
This paper develops a general framework for simultaneously solving for the flow and segregation of polydisperse granular materials. At its heart lies the partially regularized incompressible μ(I)-rheology of  and the polydisperse segregation theory of , which is generalized here to allow for different diffusion rates and segregation directions between the constituents. The coupling between the models is crucial and can be very complex. Three primary coupling mechanisms are identified: (i) advection of the particle concentrations by the bulk velocity, (ii) feedback of the particle-size and/or frictional properties on the bulk flow field and (iii) influence of the shear-rate, pressure, gravity, particle size and particle-size ratio on the locally evolving segregation (Trewhela et al. 2021) and diffusion rates (Utter & Behringer 2004).
A general numerical method is developed to solve the resulting system of equations, which is implemented in OpenFOAM. In order to solve free-surface flow problems that commonly arise in both geophysical and industrial contexts, a new interface sharpening procedure is developed that uses the multi-component segregation theory to segregate excess air out of the granular material. The new method generates a sharp interface between the grains and the air and prevents the formation of mesh-dependent trapped air bubbles or air layers, which form with standard interface sharpening techniques (Rusche 2002;Weller 2008). In fluid flows, bubbles may be realistic, but in granular flows they are not, because the air can usually escape easily through the pore space. Bubble trapping in solid-like granular flows is a common problem, and the new segregation based approach to interface sharpening solves many of the issues when combining multiphase methods with granular flow theory and may be applicable to a wide range of problems.
The numerical method is used to investigate one-way coupled problems in § 5 and two-way coupled problems in § § 6 and 7. The advantage of investigating one-way coupled problems is that it allows the numerical method to be extensively tested against exact solutions for (i) concentration shock wave development (figure 3), (ii) steady-state Bagnold flow (figure 4a), (iii) steady-state concentration profiles (figure 4b) and (iv) the formation of steadily travelling flow fronts (figures 5 and 6). These simple one-way coupled simulations also qualitatively show how large and small particles are advected in a spatially and temporally evolving bulk flow field, allowing, for instance, the formation of a large rich flow front to develop (figure 3). When the large particles are more frictional than the fines (in § 6) the large rich flow front slows down, and a bulbous head develops (see figure 8) that is relevant for geophysical flows (Denissen et al. 2019).
To provide a quantitative test of the model, Trewhela et al.'s (2021) experimental scaling law for segregation is implemented together with a generalization of Utter & Behringer's (2004) diffusivity in § 6.3. Figure 9(a) shows very good agreement with Tripathi & Khakhar's (2011) DEM simulation data for the steady-state concentration profiles with depth, without the need for any fitting parameters. The frictional feedback arises through the use of the generalized inertial number (3.2), which is based on the average local grain sized defined in (3.1). For an inclined flow down a plane, this monotonically decreases the velocity at all heights as the proportion of large particles increases. This general trend is also seen in Tripathi & Khakhar's (2011) DEM data (figure 9b), although the fits are not precise. The fact that Tripathi & Khakhar's (2011) pure large and pure small simulations do not obey the Bagnold scaling precisely, suggests that more accurate DEM simulations are required to fully test the model, in particular, there may be an influence from the basal roughness, which does not change as the mean grain size changes, and their data at 30 % small particle concentration appear anomalous.
As a demonstration of the potential of the model, the fully coupled flow in a square rotating drum is computed in § 7. Such a configuration is a real challenge for current models, because the flow field is not steady and cannot easily be prescribed or approximated from DEM simulations. The numerical simulations (figures 10-12) show that the fully coupled model is able to compute the spatially evolving velocity, pressure and concentration fields as a function of time, and that a petal-like concentration pattern spontaneously forms in the rotating deposit, which is qualitatively very similar to that observed in the experiments of Hill et al. (1999), Ottino & Khakhar (2000) and Mounty (2007). Precise experimental comparison is not possible at this stage, however, because the experiments are strongly influenced by wall friction, making the free-surface avalanche thinner and faster than in the absence of sidewalls (Jop et al. 2005). Computations that include sidewall friction are possible, but will require finer meshes to resolve the segregation within the avalanche, and consequently will take much longer to run.
The examples investigated in this paper provide the briefest glimpse at what is possible within this powerful new theoretical and computational framework. There is still much work to be done to fully understand the feedbacks and how they can affect real world problems of practical interest. In some situations the feedback may be relatively subtle, i.e. the quantitative values for the velocity, pressure and concentrations are changed, but they don't have a big impact on the subsequent flow (see e.g. figures 10-12). However, in other situations these (sometimes small) quantitative changes can induce fundamental qualitative change in the solutions. A prime example of this is the formation of a bulbous large rich head (Gray & Ancey 2009;Denissen et al. 2019), which is calculated in two dimensions for the first time in § 6.2 and shown in detail in the supplementary movie to figure 8. In three dimensions this solution can become unstable to spanwise perturbations and instead generates a series of leveed flow fingers (Pouliquen et al. 1997;Woodhouse et al. 2012;Baker et al. 2016b) that are directly relevant to the self-channelization of snow avalanches, debris flows and volcanic pyroclastic flows (Pierson 1986;Iverson & Vallance 2001;Johnson et al. 2012;Rocha et al. 2019). Much less is known about the feedbacks between segregation and flow in industrial problems, but they most definitely occur (see e.g. Zuriguel et al. 2006). It is hoped that our new found understanding can also be exploited in future to improve, and control, the flowability of bulk solids as well as mitigate the worst effects of particle segregation.