**Abstract** Granular materials composed of a mixture of grain sizes are notoriously prone to segregation during shaking or transport. In this paper, a binary mixture theory is used to formulate a model for kinetic sieving of large and small particles in thin, rapidly flowing avalanches, which occur in many industrial and geophysical free-surface flows. The model is based on a simple percolation idea, in which the small particles preferentially fall into underlying void space and lever large particles upwards. Exact steady-state solutions have been constructed for general steady uniform velocity fields, as well as time-dependent solutions for plug-flow, that exploit the decoupling of material columns in the avalanche. All the solutions indicate the development of concentration shocks, which are frequently observed in experiments. A shock-capturing numerical algorithm is formulated to solve general problems and is used to investigate segregation in flows with weak shear.

## 1. Introduction

The blending or separation of grains and powders is of considerable practical importance to the pharmaceutical, bulk chemical, mining and food industries, which process huge quantities of granular materials. US production alone is estimated at a trillion kilograms per year (Shinbrot & Muzzio 2000). In some processes, such as those in the mineral processing industry, size segregation is a desired and useful effect (Wills 1979). In many others, where the aim is to mix two cohesionless granular materials into a consistent blend (Metcalfe *et al*. 1995; Shinbrot *et al*. 1999; Gray 2001), size segregation can be a source of considerable frustration. Incorrect blending or separation can degrade the quality of a product, and in some situations it can have serious safety implications.

There are a number of mechanisms for the segregation of dissimilar grains in granular flows (Bridgwater 1976), including inter-particle percolation, convection (Ehrichs *et al*. 1995), inertia, collisional condensation (Jenkins 1998), differential air drag, clustering (Mullin 2000) and ordered settling. This paper focuses on kinetic sieving (Savage & Lun 1988), which is the dominant mechanism in dense granular free-surface flows. The basic idea is that as grains avalanche downslope, the local void ratio fluctuates, and the small particles fall into gaps that open up beneath them because they are simply more likely to fit into the available space than the large ones. The small particles, therefore, migrate towards the bottom of the flow and lever the large particles upwards by force imbalances. In dry frictional cohesionless flows, this process is so efficient that zones with 100% large and 100% small particles develop, which are separated by a sharp concentration jump (Savage & Lun 1988; Vallance & Savage 2000). However, in more energetic collisional flows this mechanism is less efficient, because of diffusive mixing, and the sharp transitions are smoothed out (Jenkins & Yoon 2001).

Kinetic sieving within shallow granular avalanches is the key mechanism for the formation of layers, which have high concentrations of large particles overlying high concentrations of fines. In geology this is known as ‘inverse’ or ‘reverse’ grading (Middleton & Hampton 1976), and distinguishes deposits from granular flows, such as rockfalls, lahars (Vallance 2000), debris flows and pyroclastic flows (Iverson & Vallance 2001), from ‘normally’ graded deposits that are formed as large particles drop out of a water suspension before the fines. If an inversely graded layer is sheared, large particles tend to migrate towards the front of the flow and smaller ones towards the rear. This local size distribution can have a controlling influence on the dynamics of the bulk flow. For instance, in some flows, the larger particles at the front experience greater friction with the underlying topography than the rest of the material, which can lead to lateral instability and the formation of lobes and fingers (Pouliquen *et al*. 1997). These are frequently observed in pyroclastic lobate deposits (Iverson & Vallance 2001) and in the formation of levees in debris flows (Iverson 1997).

Many free-surface flows in industry and geophysics, interact closely with a solid body of grains beneath. In industrial rotating tumblers, for instance, a series of avalanches at the free-surface can leave inversely graded deposits that form amazingly ordered and beautiful patterns in some regimes (Gray & Hutter 1997, 1998) and chaotic patterns in others (Hill *et al*. 1999). Figure 1 shows patterns produced by rotating a thin circular cylinder, which is partially filled with a mixture of large and small grains. As the drum is rotated, avalanches are released intermittently. The passage of each avalanche leaves an inversely graded stripe parallel to the free-surface, which is then rotated and buried by subsequent avalanches. For large fill heights, a Catherine wheel pattern is formed, while lower fill heights generate a palm-leaf pattern, which has arms that radiate outwards in the opposite sense. Similar deposition mechanisms are responsible for the formation of stratification and segregation patterns in heaps and silos (Williams 1968; Makse *et al*. 1997; Gray & Hutter 1997; Baxter *et al*. 1998), which are as equally familiar to those handling bulk solids in industry as they are to geologists (Fineberg 1997).

Although kinetic-sieving is of fundamental importance in many industrial, geophysical and pattern-forming processes, there has been little theoretical work to date. Savage & Lun (1988) used statistical mechanics and information entropy ideas to derive the first theory, which was able to predict the steady-state particle size distribution in a steady uniform flow. The model has many good features, but one of its weaknesses is that it predicts segregation even in the absence of gravity. The percolation of small particles in the kinetic sieving mechanism is, however, a fundamentally gravity driven process, so kinetic sieving models should not predict segregation if gravity is ‘turned off’. Note that segregation may occur in the absence of gravity, but it is driven by spatial gradients in the fluctuation energy of dissimilar grains (Jenkins & Mancini 1987; Xu *et al*. 2003). The problem with Savage & Lun's (1988) model stemmed from particles being instantaneously transferred between layers in their capture mechanism. Here individual constituent momentum balances are used to model the percolation of the large and small particles, which provides a natural way of introducing gravity into the theory. The resulting equations are closely linked to Savage & Lun's (1988) model, but are able to compute both the temporal and spatial evolution of the particle size distribution in any incompressible three-dimensional granular free-surface flow.

## 2. Mixture framework and conservation laws

A simple two constituent theory is formulated in which the interstitial pore space is incorporated into the volume fractions of the grains. The granular material is therefore assumed to be a bi-disperse mixture of ‘large’ and ‘small’ particles, which occupy volume fractions *ϕ*^{l} and *ϕ*^{s}, per unit mixture volume, respectively. With these definitions, the large and small particle granular volume fractions must sum to unity:(2.1)Mixture theory (e.g. Truesdell 1984; Morland 1992) assumes that every point in the material is simultaneously occupied by both phases, so that overlapping partial densities, *ρ*^{ν}, partial velocities, **u**^{ν}, and partial pressures, *p*^{ν}, can be defined for each constituent *ν*. The constituent letters l and s will be used for the ‘large’ and ‘small’ particles throughout this paper. Each phase must satisfy individual balance laws for the conservation of mass(2.2)and momentum(2.3)where ⊗ is the dyadic product and *ρ*^{ν}* g* is the gravitational force. The interaction force,

*β*^{ν}, is the force exerted on phase

*ν*by the other constituent. By definition, the sum over the two constituents is equal to zero,

*β*^{l}+

*β*^{s}=0. The bulk density,

*ρ*, and the bulk pressure,

*p*, are defined as the sum of the partial densities and partial pressures:(2.4)A key feature of mixture theory is how partial variables are related to their physical, or intrinsic, counterparts. In standard mixture theory, the partial and intrinsic velocity fields are identical, while the other fields are related by linear volume fraction scalings(2.5)where the superscript * denotes an intrinsic variable. Note, that in this simple two constituent formulation, the interstitial pore space is assumed to be incorporated into each phase, so the intrinsic densities,

*ρ*

^{ν*}, are in fact the mean solids fraction times the bulk solid density. This also implicitly assumes that the mean solids fraction is constant.

## 3. The particle size segregation model

Let O*xyz* be a coordinate system with the *x*-axis pointing down a chute inclined at an angle *ζ* to the horizontal, the *y*-axis across the chute, and the *z*-axis being the upward pointing normal as shown in figure 2. The large and small particles are assumed to have the same constant density,(3.1)which is necessarily equal to the bulk density, *ρ*. If the normal acceleration terms are negligible, the sum of the momentum balance components (2.3) over large and small constituents implies(3.2)where *g* is the constant of gravitational acceleration. Since *ρ* is constant and the free-surface is traction free, equation (3.2) can be integrated through the avalanche depth, *h*, to show that the bulk pressure is hydrostatic(3.3)The key idea behind the kinetic sieving model is that while the small particles percolate through the matrix, they support less of the overburden pressure than they should, and the large particles must therefore carry proportionately more of the load. A new pressure scaling is therefore introduced(3.4)where the factors *f*^{l} and *f*^{s} determine the proportion of the hydrostatic load which is carried by the large and small particles. This departs from standard mixture theory, which assumes as in equation (2.5) that the pressures scale linearly with the volume fraction. Specific forms for these factors will be proposed later, but it is worth noting that equation (2.4) implies that they must sum to unity, *f*^{l}+*f*^{s}=1. Experimental observations of the kinetic sieving process suggest that there is an analogy with the percolation of fluids through porous solids, and Darcy's law therefore motivates an interaction drag of the form (e.g. Morland 1992)(3.5)where *c* is the coefficient of inter-particle drag and is the bulk velocity. The second term on the right-hand side is simply a resistance term, while the first term combines with ∇(*f*^{ν}*p*) in the momentum balance to leave *f*^{ν}∇*p*. This ensures, as in Darcy's law, that the percolation process is driven by intrinsic rather than partial pressure gradients.

The large and small particle percolation velocities are assumed to be of the same order of magnitude as the normal bulk velocity, but much smaller than typical bulk downstream velocities. To reflect this, the constituent velocities in the down and cross-slope directions are assumed to be equal to the bulk down and cross-slope velocity components,(3.6)The normal constituent velocities, *w*^{ν}, are obtained by substituting equations (3.3)–(3.5) into the normal component of the momentum balance equation (2.3), with the assumption that the normal acceleration terms are negligible, to give(3.7)The significance of the pressure scalings *f*^{ν} is now clear. If *f*^{ν}>*ϕ*^{ν}, then the particles will rise, if *f*^{ν}<*ϕ*^{ν}, the particles will fall, and if *f*^{ν}=*ϕ*^{ν} there will be no motion relative to the bulk normal flow. The functions *f*^{ν} must satisfy the constraint that when only one type of particles are present, they must support the entire load, i.e.(3.8)The simplest non-trivial functions that satisfy the constraints (3.8) and the condition that *f*^{l}+*f*^{s}=1 are(3.9)where the first terms on the right-hand side are the usual volume fraction scalings and the second terms represent a perturbation of non-dimensional magnitude *B*. Note, that since *ϕ*^{l}+*ϕ*^{s}=1, the perturbation is zero whenever a pure phase is present. When the functions (3.9) are substituted into equation (3.7) they imply that the large and small particle percolation velocities relative to the bulk are(3.10)respectively, where(3.11)is the mean segregation velocity. Therefore large particles move up through the matrix at a velocity proportional to the volume fraction of small particles, while the small particles drain down at a velocity proportional to the volume fraction of large particles. In both cases, the segregation stops when a 100% concentration of that constituent is reached.

An equation to compute the volume fraction *ϕ*^{s} can be formulated by substituting equations (3.6) and (3.10) into the mass balance equation (2.2), for the small particles to give(3.12)The bulk flow * u*=(

*u*,

*v*,

*w*) will be prescribed in this paper, but, in general, the velocity can be computed by solving an existing avalanche model, such as the shallow water type theories of Grigorian

*et al*. (1967), Kulikovskii & Eglit (1973), Eglit (1983) and Gray

*et al*. (2003), or the Mohr–Coulomb models of Savage & Hutter (1989, 1991) and Gray

*et al*. (1999). Both sets of models assume that the granular avalanche is incompressible(3.13)with constant uniform density and a hydrostatic pressure distribution through their depth. This is consistent with the bulk density and pressure fields assumed in this model, and the incompressibility condition is recovered by summing the mass balance equations (2.2) over the large and small constituents and using equation (3.1).

### (a) Comparison with Savage and Lun's theory

Savage & Lun (1988) used a statistical mechanics approach to derive a segregation equation to compute steady-state particle size distributions in a uniform steady two-dimensional flow field. Direct comparison with equation (3.12) is difficult, as Savage and Lun formulated their theory in terms of the layer number density ratio, *η*, and the particle diameter ratio, *σ*, instead of volume fractions. However, an analysis of the definitions of the partial densities in their equations (3.25) and (3.27) reveals that(3.14)In addition, comparing (3.12) with their equations (5.3) and (5.4), it is clear that and are simply the large and small particle percolation velocities relative to the bulk. Savage & Lun (1988) define these velocities in their equations (6.4) and (6.3), which using (3.14) can be expressed as(3.15)where their mean segregation velocity,(3.16)is a function of the diameter of the large particles, *D*_{a}, the shear rate, d*u*/d*z*, and the non-dimensional volume averaged velocities, and . These are in turn dependent on a further six variables. A comparison of equations (3.10) and (3.15) shows that the two theories have the same structure, differing only in their definitions of the segregation velocities, *q* and *q*^{SL}, in equations (3.11) and (3.16). A significant advantage of the new theory, is that the segregation velocity is dependent on the normal component of gravity, *g* cos *ζ*, which automatically defines the direction for segregation and ensures that the kinetic sieving process does not act in the absence of gravity. Ultimately, it may be possible to incorporate some of the other more complex dependencies of *q*^{SL} into *q*, but, in the present paper, the case of constant *q* is investigated as it is the simplest mathematical structure that leads to segregation.

### (b) Non-dimensionalization

Avalanche models all exploit the shallowness of the flow to derive a system of depth-averaged mass and momentum equations for the thickness and the mean downslope velocity. Anticipating that the bulk flow will be computed using such models, the variables are non-dimensionalized by the standard avalanche scalings:(3.17)where *U* is a typical downslope velocity magnitude, and the typical avalanche length, *L*, is much larger than the typical thickness, *H*. Dropping the tildes and the superscript s, the segregation equation (3.12) becomes(3.18)where the non-dimensional segregation number(3.19)is the ratio of the mean segregation velocity to typical magnitudes of the normal bulk velocity, *w*. The non-dimensional form of the incompressibility condition (3.13) can be used to simplify the conservation form of the segregation equation (3.18) to(3.20)which, when * u* is given, is a classical first order quasi-linear equation for the volume fraction of small particles. When

*S*

_{r}≡0, there is no segregation and equation (3.20) reduces to the

*tracer equation*. This has been used by Gray & Tai (1998) and Gray

*et al*. (2000) to model the formation of stratification patterns in a pre-segregated bi-disperse mixture. The non-dimensional segregation number,

*S*

_{r}, determines the strength of the segregation. Strong segregation is usually observed (Savage & Lun 1988) when there is a significant gradient in the downslope velocity through the avalanche depth, but segregation can also occur, over longer distances, in avalanches with weak shear. Avalanche models usually assume a

*plug-flow*regime with uniform down- and cross-slope velocity profiles, but vertical structure can easily be incorporated by the inclusion of

*shape factors*(e.g. Savage & Hutter 1989) in the depth averages of

*u*

^{2},

*uv*and

*v*

^{2}in the momentum transport terms. Both strong and weak shear can be generated in laboratory experiments, and it is therefore of interest to see what effect they have on the resulting particle size distribution.

## 4. Steady-state segregation in steady uniform flows

Steady-state segregation generated by a homogeneous inflow of particles in a simple shear field was first investigated by Savage & Lun (1988) who derived exact solutions for the flow and compared them to a series of detailed laboratory experiments. A similar set of problems may be formulated for the new theory, if the velocity(4.1)in an avalanche of constant unit height. Assuming that there are no cross-slope gradients, the conservation form of segregation equation (3.18) reduces to(4.2)which must be solved subject to the condition that a homogeneous mixture of concentration *ϕ*_{0} enters the chute at *x*=0,(4.3)and there is no normal flux of particles out through the free-surface or the base(4.4)

### (a) Characteristics

The conservative segregation equation (4.2) can be rewritten as a simple first order quasi-linear equation by expanding out the derivatives to give(4.5)Solutions to equation (4.5) may be constructed by the method of characteristics, which implies that *ϕ* is equal to the constant *ϕ*_{λ} on the characteristic curve given by(4.6)As *ϕ*_{λ} is constant and *u* is a function of *z*, this is a separable equation that can easily be integrated once the velocity field is known. Solutions for general velocity fields can, however, be constructed by defining a depth-integrated velocity coordinate(4.7)which increases monotonically with increasing *z*, provided *u*≥0 and any points with zero velocity are isolated. Since the avalanche velocity magnitude is set by the scalings (3.17), we may assume without loss of generality that *ψ*=1 at the free-surface *z*=1. In these coordinates (4.6) reduce to(4.8)which is independent of the prescribed velocity field, and has the solution(4.9)where (*x*_{λ}, *ψ*_{λ}) is the initial position of the characteristic. In transformed variables, all characteristics are therefore straight lines, the gradients of which are set by the non-dimensional segregation parameter and the small particle concentration. For regions with 100% concentration of large particles, the characteristics propagate downwards with gradient −*S*_{r}; for regions with 100% concentration of small particles, the characteristics propagate upwards with gradient *S*_{r} and for particles entering the domain with concentration *ϕ*_{0}, the characteristics have gradient *S*_{r}(2*ϕ*_{0}−1), which may point either upwards or downwards, depending on the value of *ϕ*_{0}. The characteristics are illustrated for an inflow concentration of 60% in figure 3. The position of the characteristics in physical space can be calculated by inverting transformation (4.7) once *u*(*z*) is prescribed.

### (b) Segregation jump condition

Savage & Lun (1988) observed the development of sharp concentration jumps in their laboratory experiments. Once such shocks appear, the segregation field equations (4.2) and (4.5), are no longer valid and instead a jump condition (e.g. Chadwick 1976) must be applied across the discontinuous interface. The jump condition can be derived from an integral version of the conservative form of segregation equation (4.2)(4.10)Assuming that there is a jump in *ϕ* at *z*=*s*(*x*), equation (4.10) becomes(4.11)where the plus and minus superscripts denote evaluation of the limits on either side of the discontinuity. Interchanging the order of differentiation and integration by applying Leibniz's rule (Abramowitz & Stegun 1970),(4.12)where the jump bracket is the difference of the enclosed quantity on the forward and rearward sides of the shock. Shrinking the domain [*L*_{1},*L*_{2}] onto the shock by taking the limits *L*_{1}→*s*^{−} and *L*_{2}→*s*^{+} yields the jump condition(4.13)where *s*′=d*s*/d*x*. This can be rearranged to give(4.14)which is the ordinary differential equation for the height of the shock. Using depth-integrated velocity coordinates, defined in equation (4.7), this can be transformed to(4.15)which is again independent of the assumed velocity profile.

### (c) Shock solutions in mapped coordinates

Savage & Lun's (1988) experiments showed that a layer of 100% fine particles developed near the base and that there was a sharp concentration jump between this region and the homogeneous inflowing mixture. There is a simple physical explanation for this. As the mixture flows into the chute, small particles drain down through the matrix, and in turn lever large particles towards the surface. The net effect is that the local volume fraction of small particles remains at the inflow concentration *ϕ*_{0} throughout most of the flow, just as the solution by the method of characteristics suggests. However, at the lower boundary there are no more large particles to be levered up, and the no-normal flux condition, equation (4.4), implies that a region of 100% fine particles develops at the base. This layer becomes progressively thicker downstream, because no large particles can be supplied from the pure phase.

The position of the lower shock, which separates the fine particles from the homogeneous mixture, can be computed from (4.15) by substituting *ϕ*^{+}=*ϕ*_{0} and *ϕ*^{−}=1 and integrating subject to the boundary condition that *ψ*=0 at *x*=0. This implies(4.16)where the subscript 1 is used to denote the lower shock. In depth-integrated velocity coordinates, this grows linearly with downstream distance, as shown in figure 3, but, in general, when mapped back to physical coordinates, it will describe a curve. Specific results for a range of velocity profiles will be analyzed at the end of this section, once the complete solution has been constructed.

Savage & Lun (1988) also observed the development of a similar concentration shock near the free-surface, between a region of 100% coarse particles and the homogeneous mixture. This time the pure phase of large particles is generated because there are no more small particles to fall down through the matrix. Using exactly the same arguments as above, the position of the top shock, *ψ*_{2}, can be computed by substituting *ϕ*^{+}=0 and *ϕ*^{−}=*ϕ*_{0} into equation (4.15) and integrating subject to the boundary condition *ψ*_{2}=1 at *x*=0, which implies(4.17)This falls linearly from the free-surface and meets the bottom shock at *x*=1/*S*_{r} at a height *ψ*=*ϕ*_{0} in depth-integrated velocity variables. When the two shocks merge a third shock is formed, between the pure phases of small and large particles, creating a *triple-point*. The shock position is determined by substituting *ϕ*^{+}=0 and *ϕ*^{−}=1 into equation (4.15) and integrating to give(4.18)The solution consists of three domains of constant concentration, which are separated by straight shocks in depth-integrated velocity coordinates (*x*, *ψ*). This is illustrated in figure 3. The full structure is now clear. At *x*=0, the homogeneous mixture enters the domain and the initial concentration, *ϕ*_{0}, is swept into the triangular region adjacent to the *ψ*-axis by the characteristics. At the base and free-surface, there are no more large or small particles to propagate through the domain, and pure phases of small and large particles develop. Within these pure regions no further segregation takes place, and the particles move downstream along constant height trajectories by virtue of (3.10). The small particles have characteristics that propagate upwards and the large ones that propagate downwards, and these eventually intersect, either with one another, or, with the characteristics from the homogeneous domain, to generate the three shocks.

### (d ) Physical solutions

The beauty of the depth-integrated velocity coordinates is that the solution given by equations (4.16)–(4.18) is valid for all velocity fields that yield a well-defined mapping from depth-integrated to physical coordinates. This is guaranteed provided *u*(*z*)≥0, and any points where the velocity is zero are isolated. In this paper, the solution for several different physical velocity fields will be constructed by solving integral (4.7) to obtain *ψ*=*ψ*(*z*). The first set of velocity fields all have linear profiles(4.19)where the parameter *α* is used to generate plug-flow for *α*=1, simple shear for *α*=0 and linear shear with basal slip, for intermediate values. The integral (4.7) implies that the depth-integrated velocity coordinate(4.20)which has the property that at the free-surface *ψ*(1)=1, as required. This can easily be inverted to give the position of the shocks in physical space(4.21)The exact solutions are illustrated in figure 4, for a series of linear velocity profiles generated by parameter values of *α*=0, 1/2 and 1, an initial concentration *ϕ*_{0}=50 and 30% and for a segregation number *S*_{r}=1. Plug-flow is the simplest case, as the physical and transformed coordinates are identical. The shocks and characteristics are therefore all straight lines, as described above. For 50% inflow concentration, the third shock lies at *z*=0.5 to create an inversely graded layer of large particles overlying small particles. At 30% inflow concentration the solutions are similar, but the third shock lies at *z*=0.3. In the case of linear shear with slip (*α*=1/2), the shocks and characteristics are mapped to physical space using the full quadratic mapping defined in equation (4.21). The solution, therefore, has the same basic structure as plug-flow, but the upper and lower shocks are now curved and the third concentration shock is higher up, to reflect the fact that there is a greater mass flux close to the free-surface than at the base.

### (e) Comparison of simple shear with Savage and Lun

The case of simple shear is particularly interesting as it can be compared to the results of Savage & Lun (1988). The inverse mapping (4.21) reduces to , so that the shocks are simply(4.22)which are illustrated in the lower middle panels of figure 4. The upper and lower shocks have square root profiles, with the lower one having an infinite gradient at the origin. The third shock is again straight, but it is significantly higher than the cases of plug-flow and linear shear with slip.

Savage & Lun's (1988) solutions look superficially very similar to those constructed here, but were based on an incorrect assumption. Instead of using a shock condition, they assumed that the fall line (or shock) marking the 0% fines region, was given by the characteristic curve emanating from the top left corner of the domain, i.e.(4.23)They then used the fact that, at any position along the chute, the depth-integrated flux of small particles must be equal to the inflow flux,(4.24)which, upon evaluation of the integrals, implies that the bottom shock(4.25)With the upper and lower shocks given by equations (4.25) and (4.23), Savage & Lun's (1988) theory predicts that complete segregation occurs at(4.26)Using definitions (3.14), it is easy to show that equations (4.23), (4.25) and (4.26) are equivalent to equations (6.11)–(6.14) of Savage & Lun (1988). The upper and lower shocks are steeper than those in the correct solution, and therefore the distance for complete segregation is shorter. However, the conservation of the depth-integrated flux of small particles in equation (4.24) ensures that the final height of the third shock is correct. In the dilute limit, there is only a small difference between the two solutions, as figure 5 shows, so there is still very good agreement between the new theory and the experiments of Savage & Lun (1988) and Vallance & Savage (2000).

### (f ) A velocity field that scales with the thickness to the 3/2

Steady uniform flows are observed over a range of inclination angles on bumpy chutes and exhibit an interesting scaling of the mean velocity with the avalanche thickness to the 3/2 (Vallance 1994). Pouliquen (1999) has developed a new drag law for rough slopes based on this scaling and recently Silbert *et al*. (2001) have performed molecular dynamics simulations, which also show the same scaling law. Based on Bagnold's (1954) constitutive law for the grain inertia regime in which the shear stress is proportional to the square of the strain-rate, , Silbert *et al*. (2001) derived an expression for the downslope velocity profile,(4.27)which agrees well with their simulations. Scaling the flow depth on the avalanche thickness, *h*, and the velocity magnitude *U*=2(*ρg* sin *ζ*)^{1/2}*h*^{3/2}/(5*A*), so that the depth-integrated flux is unity at the free-surface, the non-dimensional downslope velocity becomes(4.28)Integral (4.7) implies that the depth-integrated velocity coordinate(4.29)This transformation cannot be inverted to give explicit relations for the shocks, but it is a simple matter to produce contour plots of the results, which are shown in the bottom two panels of figure 4. They look quite similar to the case of simple shear, as the bottom shock has an infinite gradient at the origin, but the third shock is not as high. Note that Louge (2003) has extended Silbert *et al*.'s (2001) theory to determine the constant of proportionality *A*^{2}.

## 5. Time-dependent segregation in steady uniform plug-flows

Segregation by kinetic sieving is usually associated with flows that are strongly sheared through their depth. However, many geophysical scale (Dent *et al*. 1998) and laboratory scale (Savage & Hutter 1989; Keller *et al*. 1998; Eckart *et al*. 2003) avalanches on smooth slopes have relatively blunt downstream velocity profiles with slip at the base. Segregation can still occur in these flows, provided the particle size difference is large enough and the agitation strong enough for the matrix to dilate sufficiently for percolation to take place. Considerable insight into segregation in this weak shear limit is provided by the case of plug-flow.

### (a) Segregation in independent columns

For uniform plug-flow in a domain of unit height the segregation equation (3.18) reduces to(5.1)where the transport velocity, *u*_{0}, can, without loss of generality, be assumed to be unity by virtue of the scalings in equation (3.17). As the velocity is independent of depth, considerable simplification can be achieved by transforming to a frame moving downstream with speed *u*_{0}. Using the change of coordinates(5.2)the segregation equation (5.1) reduces to(5.3)in the moving frame, where the primes are now dropped. Since this equation is independent of *ξ*, it implies that particle size segregation in a fixed moving column of granular material is completely independent of the segregation taking place in adjacent columns. This uncoupling is very useful for the construction of exact time-dependent solutions.

A comparison of equations (5.3) and (4.2) shows that the segregation problem in a moving column has exactly the same structure as the steady-state problem with *u*≡1, except that the spatial coordinate is now replaced by time. The solutions constructed in § 4 can therefore be applied directly here. It is a simple matter to show that an initially homogeneous column of material with concentration *ϕ*_{c} at *t*=*t*_{c} generates three straight shocks(5.4)which separate the homogeneous mixture and the pure phases in an exactly analogous manner to the problems in §4. These solutions are effectively illustrated in the top two panels of figure 4. The only difference is that the *x*-axis must now be replaced by the *t*−*t*_{c} axis.

### (b) General time-dependent solutions for plug-flow

A full time-dependent solution can be constructed for plug-flow by using the simple column solution (5.4) in a series of adjacent columns moving downstream at speed *u*_{0}. The coordinate *ξ* will be used to uniquely identify each column, by labelling them using their initial position(5.5)It follows from equation (5.2) that at a general time *t*, column *ξ* has position(5.6)Columns with positive *ξ* are therefore initially within the avalanche, while columns with negative *ξ* first enter the avalanche, at *x*=0, at time −*ξ*/*u*_{0}. The *transition point*, *ξ*=0, between the columns initially within the avalanche and those that subsequently enter it, is transported downstream with constant speed *u*_{0}, and has position(5.7)Since the columns are independent of one another, different values of the constants *t*_{c} and *ϕ*_{c} can be chosen in each column, and they may, therefore, be considered to be functions of *ξ*. A single function, *Φ*, is, therefore, used to parameterize both the initial and boundary conditions. Assuming that the homogeneous inflow concentration varies as *Φ*(*t*), and that the initial mixture has a uniform concentration equal to *Φ*(0), the column parameters *t*_{c} and *ϕ*_{c} are(5.8)Substituting these into solutions (5.4) implies that the three shocks are simply time-dependent in the region that was initially in the chute(5.9)but have both space and time dependence in the region that flows into the chute(5.10)Shock solutions (5.9) and (5.10) determine the general solution for the small particle concentration for any time-dependent function of the inflow concentration, *Φ*.

Figure 6 shows contour plots of the small particle concentration at a series of time intervals for *S*_{r}=1, *u*_{0}=1 and a constant initial and inflow concentration(5.11)Initially, the chute is filled entirely with a homogenous uniform mixture with a concentration of 50%. As soon as the system is released, the particles segregate in the interior, creating two spatially uniform shocks that separate the pure phases from the mixture. These two shocks eventually meet at *t*=1 to create a stationary third shock at *z*=1/2, which separates the large particles from the small ones beneath. As this process is taking place, new material is fed into the chute and segregates. Substituting equation (5.11) into equation (5.10), it is clear that these shocks are straight lines that are independent of time,(5.12)The transition point, therefore, marks the divide between the steady-state solution and the region of transient adjustment to the initial conditions, and it propagates downstream with speed *u*_{0}. As the initial conditions are uniform and equal to the inflow concentration, there is no further change to the solution after the transition point passes the triple-point at *x*=1 and *t*=1. The solution, therefore, reaches the same steady-state illustrated in the top left panel of figure 4.

In physical experiments the material entering the avalanche is supplied from a hopper, and once flow starts, inhomogeneities often develop, causing the inflow concentration to vary as a function of time. Exact solutions for the plug-flow regime can be constructed for this case using the method above. Let us suppose that the variation in the homogeneous inflow concentration is parameterized by(5.13)The solution is illustrated in figure 7 at a sequence of time intervals. The structure in front of the transition point is identical to that in the first problem, with two uniform shocks generated at the boundaries that propagate inwards and meet at *z*=1/2 at time *t*=1, to leave a stationary shock. However, behind the transition, in the domain controlled by the boundary conditions, the sinusoidal variation propagates into the domain and distorts the shocks, which are fully time and space dependent. At *t*=1, the two shocks meet at *x*=1 to form a triple point, which oscillates up and down with time. For *t*>1, the third shock becomes time-dependent in the region controlled by the boundary conditions and remains straight in the region determined by the initial conditions, as shown in the right middle panel of figure 7.

## 6. Computational method and results

High-resolution shock-capturing numerical methods can be used to solve equation (3.18) with prescribed velocity fields that include shearing effects. These methods have a long history starting from the classic papers of Godunov (1959), Van Leer (1979), Harten (1983) and Yee (1987), and there are now a wide range of textbooks on these powerful schemes (e.g. LeVeque 2002; Godlevski & Raviart 1996; Toro 1997). A simple high-order modified TVD Lax-Friedrichs (TVDLF) scheme is used here, which has the advantage that it does not require a knowledge of the approximate or exact Riemann solvers for the system. The method was first suggested by Yee (1989) and was extended and improved by Tóth & Odstrcil (1996).

### (a) The modified total variational diminishing Lax Friedrichs method

Godunov type operator spitting (e.g. LeVeque 2002) is used to reduce the multi-dimensional segregation equation, equation (3.18), into a series of one-dimensional problems. A regularly spaced mesh is constructed and the discretized concentration is defined at each grid point *j* and at time-step *n*. The cell interfaces lie mid-way between grid points and are denoted by the subscript *j*+1/2. The TVDLF method starts with a Hancock half-step to achieve second order temporal accuracy:(6.1)where Δ*t* is the size of the full-time-step, Δ*x* is the grid-size and *F* is the flux in the appropriate direction. The limited difference is computed using a *superbee* slope limiter:(6.2)where *a*=sgn(Δ*ϕ*_{j+1/2}) and Δ*ϕ*_{j+1/2}=*ϕ*_{j+1}−*ϕ*_{j}. To compute the full-step, upwinded left and right states, denoted *ϕ*^{L} and *ϕ*^{R}, are constructed at the half-step(6.3)and are used to calculate the mean flux at the cell interfaces(6.4)With these definitions the full-step is(6.5)where the dissipative limiter, *Ψ*^{LR}, is equal to in the original TVDLF method (Yee 1989). Here, we follow Tóth & Odstrcil (1996) who multiplied Yee's limiter by the global Courant number to obtain a less diffusive scheme. Since the maximum global wave-speed is equal to two, this implies(6.6)The modified TVDLF method is used to solve equation (3.18) in [0, 2]×[0, 1] using 300 grid points in each direction. The inflow concentration is prescribed as a function of time at *x*=0 and the velocity * u* is given throughout the domain. Along the upper and lower boundaries at

*z*=0, 1, the no normal flux condition, equation (4.4), is imposed and at

*x*=2 the material flows freely out of the domain. Initially, the concentration is assumed to be spatially uniform, and the equations are integrated forward in time using a time-step satisfying the Courant–Friedrichs–Levy (CFL) condition. The method is accurate, robust and has been extensively tested against the exact solutions derived in §§4 and 5.

### (b) Numerical results

The effect of shear is analyzed by computing numerical solutions to the problems outlined in §5 in a shearing and translating flow with *α*=0.8. Figure 8 shows the results for the case when the initial and inflow concentrations are 50%, which are comparable to those in figure 6. The solution in the uniform region determined by the initial conditions is identical to the plug-flow solution in §5. The most important effect of shear is that the transition from initial to boundary condition controlled solutions occurs at different times at different levels in the flow. This can clearly be seen in the position of the kinks in the shocks at *t*=0.7 in the lower left panel. The top transition, therefore, intersects with the lower horizontal shock at *z*=1/2 at *t*=1, but, because of the shear, it does not link with the bottom transition, which intersects with the top shock a short time later. This creates a third shock that has three clearly defined regions, a straight portion at the steady-state height defined in equation (4.18), a transitions zone, and another straight section at *z*=1/2, as can be seen at *t*=1.5. The mismatch in height is swept downstream and out of the domain, so that locally the solution is close to steady-state at *t*=2. The shear causes the mismatch interface to steepen very slowly and it eventually breaks in finite time far downstream.

Figure 9 shows the results when the inflow concentration is controlled by equation (5.13), which varies sinusoidally about a background concentration of 50%. The results are similar to those in figure 7, but the concentration contours in the inflow region tip over with downstream distance in response to the shear. When the two shocks meet a triple-point is formed, which moves from side to side, as well as up and down. This introduces oscillations into the boundary controlled section of the third shock, which are advected downstream and slowly steepen and break.

## 7. Discussion and conclusions

This paper uses binary mixture theory to derive a simple kinetic sieving model for the segregation of large and small particles in shallow granular avalanches. The model reduces to a single first order quasi-linear conservation equation (3.18), for the volume fraction of small particles. In order to solve it, the bulk velocity field in a shallow three-dimensional incompressible granular free-surface flow must either be prescribed or computed using existing avalanche models (e.g. Grigorian *et al*. 1967; Kulikovskii & Eglit 1973; Eglit 1983; Savage & Hutter 1989, 1991; Gray *et al*. 1999, 2003). A significant advantage of this theory is that the segregation velocity is explicitly dependent on gravity. This sets an orientation for the direction of segregation and ensures that the kinetic sieving process does not drive segregation in the absence of gravity.

Exact steady-state concentration solutions have been derived for general steady uniform velocity fields, by using a concentration jump condition, equation (4.13) and a coordinate mapping, equation (4.7). All solutions consist of three shocks that separate the inflowing mixture from pure phases of the large and small particles. Sufficiently far downstream, complete segregation occurs and an inversely graded layer is obtained with the large particles separated from the small ones beneath by a concentration jump. Specific solutions have been constructed for general linear velocity fields with basal slip as well as a more complex velocity profile that obeys the thickness to the power 3/2 scaling (Vallance 1994; Pouliquen 1999; Silbert *et al*. 2001). The solutions are in close agreement with Savage & Lun's (1988) original laboratory experiments, as well as those of Vallance & Savage (2000). Fully time-dependent solutions have also been constructed for the plug-flow regime by exploiting the decoupling of material columns as they are advected downstream. These solutions yield considerable insight into segregation in geophysical flows, where there is strong slip at the base. The shock-capturing modified TVDLF method can be used to compute solutions to segregation equation (3.18) with any incompressible three-dimensional velocity field, and for general initial and boundary conditions. It will, therefore, be useful for calculating numerical solutions to more general problems in future.

## Acknowledgments

This research was supported by the Royal Society through grant 22919 and by the EPSRC through a Doctoral Training Account and an Advanced Research Fellowship (GR/S50052/01 and GR/S50069/01).

## Footnotes

As this paper exceeds the maximum length normally permitted, the authors have agreed to contribute to production costs.

- Received March 8, 2004.
- Accepted October 4, 2004.

- © 2005 The Royal Society