Segregation-induced finger formation in granular free-surface flows

Geophysical granular flows, such as landslides, pyroclastic flows and snow avalanches, consist of particles with varying surface roughnesses or shapes that have a tendency to segregate during flow due to size differences. Such segregation leads to the formation of regions with different frictional properties, which in turn can feed back on the bulk flow. This paper introduces a well-posed depth-averaged model for these segregation-mobility feedback effects. The full segregation equation for dense granular flows is integrated through the avalanche thickness by assuming inversely graded layers with large particles above fines, and a Bagnold shear profile. The resulting large particle transport equation is then coupled to depth-averaged equations for conservation of mass and momentum, with the feedback arising through a basal friction law that is composition dependent, implying greater friction where there are more large particles. The new system of equations includes viscous terms in the momentum balance, which are derived from the $\unicode[STIX]{x1D707}(I)$ -rheology for dense granular flows and represent a singular perturbation to previous models. Linear stability calculations of the steady uniform base state demonstrate the significance of these higher-order terms, which ensure that, unlike the inviscid equations, the growth rates remain bounded everywhere. The new system is therefore mathematically well posed. Two-dimensional simulations of bidisperse material propagating down an inclined plane show the development of an unstable large-rich flow front, which subsequently breaks into a series of finger-like structures, each bounded by coarse-grained lateral levees. The key properties of the fingers are independent of the grid resolution and are controlled by the physical viscosity. This process of segregation-induced finger formation is observed in laboratory experiments, and numerical computations are in qualitative agreement.

Segregation-induced finger formation in granular free-surface flows 169 with experiments showing rapid vertical segregation into regions of nearly pure small and large particles (Savage & Lun 1988;Vallance & Savage 2000;Golick & Daniels 2009). When this is combined with periodic deposition it can lead to the formation of striking alternating stratified layers (Gray & Hutter 1997;Makse et al. 1997;Gray & Ancey 2009) in heaps as well as petal-like patterns in rotating drums (Hill et al. 1999;Gray & Chugunov 2006;Zuriguel et al. 2006). For dense granular flows, the dominant physical mechanisms driving segregation are thought to be kinetic sieving and squeeze expulsion (Middleton 1970;Savage & Lun 1988;van der Vaart et al. 2015). As a polydisperse material is sheared, smaller particles are more likely to be able to percolate down through cavities that open up, which in turn exerts an upward force on the larger particles. Several models have been proposed to capture this effect (e.g. Bridgwater, Foo & Stephens 1985;Savage & Lun 1988;Dolgunin & Ukolov 1995;Gray & Thornton 2005;Gray & Chugunov 2006;Marks, Rognon & Einav 2012; which all have a similar structure and describe the evolving particle size distribution for a given bulk flow. A recent review can be found in Gray, Gajjar & Kokelaar (2015).
Field studies (e.g. Pierson 1986;Iverson 2003;Lube et al. 2007) have provided strong evidence for the occurrence of particle size segregation in geophysical flows. In particular, debris flow deposits show self-organisation into leveed channels, with large particles being vertically segregated to the free surface, sheared to the flow front and then shouldered aside into coarse-grained static regions (Félix & Thomas 2004;Johnson et al. 2012). The finer material forms a lining on the inside wall of these lateral levees (Kokelaar et al. 2014), which reduces the friction in the channel and enhances the mobility of the mixed interior. Experiments at the United States Geological Survey (USGS) debris flow flume in Oregon, USA  as well as smaller-scale laboratory investigations (Deboeuf et al. 2006;Goujon, Dalloz-Dubrujeaud & Thomas 2007;Kokelaar et al. 2014) have been able to reproduce these feedback effects, with runout distances for a bidisperse material being greater than for either type of particle in pure phase. A related phenomenon is the formation of segregation-induced fingering instabilities in granular free-surface flows (Pouliquen, Delour & Savage 1997;Pouliquen & Vallance 1999;Aranson, Malloggi & Clement 2006;Malloggi et al. 2006;Woodhouse et al. 2012). These studies can be motivated by field observations of geophysical flows advancing as a series of lobate structures, for example the pyroclastic currents following the Mount St Helens eruption in July 1980 (figure 1).
Experiments are carried out using a bidisperse mixture of spherical ballotini (white, 75-150 µm diameter) and angular carborundum (brown, 305-355 µm) flowing down a plane inclined at 27 • , which is roughened by attaching a single layer of turquoise ballotini (750-1000 µm) to the base with double-sided tape (figure 2 and supplementary movie 1 available at https://doi.org/10.1017/jfm.2016.673). Initially well-mixed material is released from rest using a double gate system with an inflow height of 2 mm. As it flows down the slope, the large particles are segregated to the surface and preferentially sheared to the front. This front becomes unstable due to greater frictional forces and splits into a number of different channels, or fingers, with the internal structure of each finger resembling that of a single leveed channel.
The time scales associated with this instability are relatively short, with the early traces of fingers beginning to appear after approximately 1 s. However, fingers only develop after the segregation of large particles to the free surface and subsequent accumulation at the front, and hence the fingering time scale must necessarily be slower than that of particle size segregation. There have been several attempts to calculate this segregation rate, for example in large-scale experimental debris flows, where Johnson et al. (2012) found large particles rising at approximately 3.5 cm s −1 , or 1 % of the typical bulk downslope velocity. In laboratory experiments of dry glass beads, Wiederseiner et al. (2011) measured percolation rates of 1.5 mm s −1 , compared to average bulk velocities of 30 mm s −1 . This ratio is consistent with the discrete element model (DEM) simulations used by Staron & Phillips (2014) to calculate segregation time scales. Such segregation rates suggest that these thin flows (less than 2 mm, or approximately 10 particle diameters) rapidly segregate before the onset of the fingering instability.
There is an important distinction to be made between two different finger formation regimes that occur for different inflow conditions. For the experiments shown in figure 2 and movie 1, large quantities of granular material are loaded into the hopper, meaning grains are supplied at a constant flux for the entire observed duration. The resulting fingers are bounded by coarse-rich levees, and also have regions of pure carborundum at the rear of each channel wall, which are eroded by oncoming material from the inflow. This erosion process is particularly apparent in movie 1, where it can be seen that the 'large particle islands' creep downslope in a series of discrete surges. These islands move more slowly than the flow front, leading to finger elongation, although levees of adjacent fingers typically remain in contact. Figure 3(a) shows a close up of the experimental frontal zone, and a schematic of this behaviour is given by figure 4(a). The continuous inflow regime is representative of early fingering instability experiments (Pouliquen et al. 1997;Pouliquen & Vallance 1999), whereas later work (Woodhouse et al. 2012; used only a finite amount of material in the hopper. In this case the initial onset of finger formation is identical, but as the inflow stops and the supply wanes, erosion of the large particle islands ceases. These stationary regions then act as a barrier between adjacent channels, preventing contact and allowing distinct separated fingers to form as the remainder of the material lengthens the pre-existing fingers (figures 3b and 4b). This final phase can lead to unexpectedly long run-out distances, such as in figure 1, before eventually Segregation-induced finger formation in granular free-surface flows 171 5 cm FIGURE 2. Experiments on a plane inclined at 27 • using 80 % ballotini (white, 75-150 µm), 20 % carborundum (brown, 305-355 µm) released from rest through a double gate system of inflow thickness 2 mm. The chute is roughened with turquoise ballotini (750-1000 µm). Images show snapshots at approximate times t = 0.9 s, t = 2.6 s, t = 4.1 s, t = 6.0 s and t = 7.9 s. Supplementary movie 1 available online.
coming to rest and revealing the lubricating fine-grained levee lining (Kokelaar et al. 2014).
Continuously supplied experiments are also conducted using a monodisperse flow of small ballotini (figure 5 and supplementary movie 2). There are some small irregularities as the front advances, most likely due to imperfections of the inflow layer and on the channel bed, as well as the formation of roll waves, but the same finger structures do not form and propagation is approximately uniform across the slope. This is consistent with the work of Pouliquen (1999b), who showed that a monodisperse granular front flows with a constant velocity and well-defined shape on a rough inclined plane. The process of finger formation is therefore driven by particle size segregation. Note that experiments using pure large particles do not flow FIGURE 3. Close ups of the experimental flow fronts for (a) a continuous supply of particles from the inflow gate and (b) a finite release of granular material, where the supply has already been cutoff. In both cases a bidisperse mixture of 80 % white ballotini (75-150 µm), 20 % brown carborundum (305-355 µm) is used and the inflow thickness is 2 mm. at this slope inclination of 27 • because the angular carborundum in pure phase is too resistive. This highlights another key component of the instability mechanism, which requires the larger particles to have a higher effective friction coefficient than the smaller ones. In natural flows, the interstitial pore pressure is dissipated more rapidly through large particles, meaning that large-particle-rich regions experience greater frictional forces, even if the particles themselves are not more angular like in the experiments shown here (Iverson 1997;Johnson et al. 2012). The equivalent experiments have also been carried out using a bidisperse mixture of different sized spheres and a frontal instability does still form, although the resulting fingers have FIGURE 4. Schematic illustrating the difference between the initial onset of finger formation and fully developed fingers. (a) A continuous supply of material from the inflow gate causes the large particles at the back of the levees to be slowly eroded and move downstream. The front of the fingers propagates faster, meaning they lengthen over time, and adjacent fingers remain in contact with each other. (b) When the inflow is cutoff the regions at the rear of the levees come to rest and all remaining material flows down the pre-established channels. This leads to elongated distinct fingers with grain-free zones in between, which will eventually arrest as the flow wanes. In both diagrams shaded regions correspond to coarse-rich areas and dotted lines denote extent of the fingers at an earlier time.
weaker, less stable levee walls. In this case the geometrical properties of the two spherical species are the same, but the large particles are slightly more resistive due to their interaction with the bed roughness (Goujon, Thomas & Dalloz-Dubrujeaud 2003). On the other hand, the fingering instability does not form in experiments using rough small particles and smooth large grains, where it is found that the larger particles shear off the top of the fines, which are deposited on the chute without the formation of fingers. The above observations suggest that any theoretical model should account for both the bulk flow and the effect of particle size segregation, in particular the relative frictional differences. Pouliquen & Vallance (1999) proposed a model for these segregation-mobility feedback effects in bidisperse granular flows based on their experimental work. Depth-averaged mass and momentum balance equations were coupled to the depth-averaged concentration (representing the distribution of large and small particles) through a basal friction law that was weighted according to the evolving composition. However, this work did not explicitly model the size-segregation process, instead prescribing an initial concentration distribution and allowing it to be advected with the bulk flow. The work of Gray & Kokelaar (2010a,b) in depth integrating previous three-dimensional segregation equations (e.g. Gray & Thornton 2005) allowed the development of fully coupled avalanche-segregation models. This was exploited by Woodhouse et al. (2012), where the coupling was 174 J. L. Baker, C. G. Johnson and J. M. N. T. Graỹ 5 cm FIGURE 5. Experiments on a plane inclined at 27 • using monodisperse granular material consisting of 100 % ballotini (75-150 µm) released from rest through a double gate system of inflow thickness 2 mm. Images show snapshots at approximate times t = 0.4 s, t = 1.7 s, t = 3.0 s, t = 4.3 s and t = 5.7 s. Note the time scales are shorter than the equivalent bidisperse experiments (figure 2) as pure small particles travel faster. Supplementary movie 2 available online.
achieved through a concentration-dependent version of Pouliquen's (1999a) friction law. This model was able to capture the qualitative features of spontaneous leveed finger formation, but the authors showed that, at a critical concentration, the equations were mathematically ill posed in the sense of Joseph & Saut (1990), i.e. a linear stability analysis produced unbounded growth rates in the high wavenumber limit. The critical Froude number at which this occurred corresponded to where one of the characteristics of the shallow water equations coincided with that from the large particle transport equation (Gray & Kokelaar 2010a,b) and the system loses strict hyperbolicity. Consequently, at a specific concentration any numerical grid-scale noise grows unboundedly as the grid size tends to zero and the ill posedness manifests Segregation-induced finger formation in granular free-surface flows 175 itself in the form of grid-dependent simulations, with the number of fingers being governed by the numerical viscosity.
The Woodhouse et al. (2012) model suggests that additional physics is required to regularise the depth-averaged governing equations.  recently devised a strategy to achieve this using the µ(I)-rheology for dense granular flows (GDR MiDi 2004;da Cruz et al. 2005;Jop, Forterre & Pouliquen 2005. To leading order, they showed that this three-dimensional constitutive law only contributed via an effective basal friction, equivalent to the dynamic friction law for rough beds (Pouliquen 1999a;Pouliquen & Forterre 2002), and the depth-averaged equations reduce to a standard hyperbolic avalanche model (e.g. Gray, Tai & Noelle 2003). Using the steady uniform Bagnold velocity and lithostatic pressure profiles (GDR MiDi 2004) they were able to include the gradient of the depth-averaged in-plane deviatoric stress into the downstream momentum balance. These higher-order viscous terms represent a singular perturbation to the system and in many situations they can be neglected. However, strong evidence for their inclusion is provided by roll waves, where the standard shallow water avalanche equations are unable to predict the cutoff frequency observed in experiments (Forterre & Pouliquen 2003). With viscous terms, the depth-averaged µ(I)-rheology is able to predict this cutoff for a wide range of Froude numbers and slope angles without any fitting parameters .
In addition, Edwards & Gray (2015) showed that the extra terms play a crucial role in the formation of steadily propagating erosion-deposition waves on erodible beds. Baker, Barker & Gray (2016) recently proposed a two-dimensional extension of the equations to account for lateral variation, and applied the model to steady uniform channel flows. The generalised viscous terms give rise to downslope velocities with cross-slope profiles, another physical feature not captured by classical shallow-water models. These very promising results for monodisperse flows suggest that  depth-averaged µ(I)-rheology could provide the dissipative mechanism to regularise the depth-averaged segregation-mobility feedback equations. This paper therefore describes how to generalise their work into a bidisperse set-up, and shows that the resulting model is mathematically well posed. A two-dimensional (downslope and lateral) extension of the system of equations, based on the work of Baker et al. (2016), admits numerical solutions showing the formation of fingering instabilities on an inclined plane, with the key finger characteristics being independent of the grid resolution and controlled by the newly introduced physical viscosity.

A depth-averaged model for particle size-segregation
Consider a Cartesian coordinate system Oxz with the x-axis pointing downslope at an angle ζ to the horizontal and the z-axis being the upward pointing normal (figure 6). A bidisperse mass of granular material is assumed to lie between a free surface at z = s(x, t) and rigid base at z = b(x), so that the flow thickness is h(x, t) = s − b. Denoting the volume fraction of small particles as φ ∈ [0, 1] (so that the proportion of large particles is 1 − φ), the evolving concentration distribution can be modelled by a general segregation-diffusive-remixing equation (e.g. Bridgwater 1976; Savage & Lun 1988;Dolgunin & Ukolov 1995;Gray & Chugunov 2006;Gajjar & Gray 2014), where the bulk velocity u has components (u, w) in the downslope and normal directions respectively. The first three terms on the left-hand side represent the advection of the concentration with the bulk flow, whereas the fourth term accounts for vertical segregation. The flux function Q(φ) 0 satisfies Q(0) = Q(1) = 0 to ensure the segregation mechanism shuts off in the monodisperse limits. Different functional forms for Q have been proposed, including a simple quadratic Q(φ) = qφ(1 − φ), (Gray & Thornton 2005) or skewed cubic Q(φ) = qφ(1 − φ)(1 − γ φ), (Gajjar & Gray 2014;van der Vaart et al. 2015) the latter being motivated by experimental observations of asymmetric segregation, which has also been found from discrete particle method simulations (Tunuguntla, Bokhove & Thornton 2014). The exact dependence will not be important in this paper. The right-hand side of (2.1) represents diffusive remixing, where the diffusivity D may, in general, depend on the flow variables. The segregation equation (2.1) is subject to kinematic boundary conditions, where subscripts b and s denote evaluation of the velocity field at the base and free surface respectively. In addition, there is no flux of either large or small particles across the boundaries, i.e.
(2.4) Segregation-induced finger formation in granular free-surface flows Following Gray & Kokelaar (2010a,b), the segregation-diffusive-remixing equation (2.1) may be integrated through the avalanche thickness using Leibniz' rule (Abramowitz & Stegun 1970) to interchange the order of differentiation and integration, giving are the depth-averaged small particle concentration and small particle flux respectively. The kinematic and no-flux boundary conditions (2.2)-(2.4), ensure that the squarebracketed terms disappear and the depth-integrated segregation equation (2.5) reduces to The model is closed by deriving expressions relating the depth-averaged concentration flux to the depth-averaged downslope velocity, the latter being defined analogously to (2.6) asū Since bidisperse flows have been observed to rapidly segregate into inversely graded layers (Gray & Hutter 1997;Gray & Ancey 2009), Gray & Kokelaar (2010a,b) suggested using a concentration profile representing a layer of pure small particles lying on top of a layer of pure large particles, where z = l(x, t) denotes the height of the separating interface. In addition, the bulk velocity is assumed to take the form whereẑ = (z − b)/h is the rescaled vertical coordinate and f is the vertical shear profile, which should be an increasing function to ensure surface velocities are greater than those at the base and should also satisfy (2.11) to be consistent with the definition (2.8). Gray & Kokelaar (2010a,b) used families of linear shear profiles given by to derive their depth-averaged segregation equation, where the parameter α ∈ [0, 1] controls the relative amount of shear and basal slip. These were also employed by Johnson et al. (2012) to reconstruct the full velocity field at the USGS flume. Whilst simple linear profiles capture the basic features of the flow, a more physically accurate choice is the Bagnold velocity profile, which can be derived as the steady uniform solution to the three-dimensional µ(I)-rheology for granular flows (e.g. GDR MiDi 2004;. Substituting the inversely graded concentration (2.9) and velocity profile (2.10) into the flux integral in (2.6) gives which may then be inserted into the depth-integrated segregation equation (2.7) to give (2.16) The first two terms in (2.15) represent advection of the depth-averaged concentration with the bulk flow, and the third term captures the preferential shearing of the large particles to the flow front (the minus sign implies that fines are transported to the rear). For this reason it is referred to as the 'large particle transport equation' and is a more general version of that derived by Gray & Kokelaar (2010a,b) and Woodhouse et al. (2012). The form of the 'transport function' G depends on the choice of shear profile, with the linear shear profile (2.12) leading to the quadratic as in Gray & Kokelaar (2010a), and the Bagnold shear profile (2.13) giving The functions (2.17) and (2.18) have similar forms, with both satisfying G(0) = G(1) = 0, meaning the concentration is simply advected at the same speed as the bulk flow in both of the monodisperse limits. The Bagnold transport function (2.18) is skewed slightly towards smaller concentrations of small particles. However, the difference is relatively small (<7 % of the maximum amplitude) and (2.18) may be approximated using a quadratic of the form (2.17) (figure 7). A value α = 1/7 is chosen to ensure that the total area under the two curves, and hence the mean transport rate across all different concentrations, is the same, and such a fitted quadratic for G shall be assumed throughout this paper. This makes subsequent computations more straightforward, since the (1 −φ) 3/2 term in (2.18) results in complex values if round-off errors causeφ to be slightly greater than unity. Though the linear profile Segregation-induced finger formation in granular free-surface flows (2.12) with α = 1/7 is qualitatively different to the Bagnold shear (2.13) due to the non-zero basal slip velocity, the remainder of this work does not distinguish between the velocity at different vertical positions, meaning this simplification is appropriate when dealing with depth-averaged quantities.
Note the similar structure of the original segregation equation (2.1) and the large particle transport equation (2.15), with the vertical segregation in the former being replaced by lateral segregation in the latter. Also note that it is possible to reformulate (2.15) in terms of the small particle layer thickness, η(x, t) = l − b, using the fact that η = hφ, or the thickness of the large particle layer, κ(x, t) = h − η, as described in Gray & Kokelaar (2010a,b). Here it shall be left in terms of the depth-averaged concentration of small particlesφ because this is more representative of what would actually be seen in overhead views of bidisperse experiments.

Segregation-mobility coupling
The large particle transport equation (2.15) may be solved for the depth-averaged concentrationφ for a prescribed flow thickness h and bulk velocityū (e.g. Gray & Kokelaar 2010a,b). In some cases h andū can be inferred from experimental measurements ) but typically they are unknown and need to be solved for as part of the problem. Furthermore, it is expected that the evolving concentration distribution will feed back on the bulk motion and this coupling should be built into the model. The equations representing conservation of mass and momentum for the bulk flow are where g is the constant of gravitational acceleration. The shape factor χ = u 2 /ū 2 in (3.2) depends on the form of the velocity profile with depth. The Bagnold profile (2.13) gives a value χ = 5/4 but it shall be assumed that χ = 1 for simplicity, since non-unity values change the characteristic structure of the inviscid equations and cause problems near zero-thickness regions (Hogg & Pritchard 2004). This is common across the granular flow literature (Grigorian, Eglit & Iakimov 1967;Savage & Hutter 1989;Gray, Wieland & Hutter 1999;Pouliquen & Forterre 2002), even though it is formally inconsistent with the sheared velocity profile. The source terms S are due to a combination of gravity, effective basal friction and changes in basal topography (e.g. Gray et al. 2003), where sgn is the sign function and ensures friction always opposes the direction of motion. The effective basal friction coefficient µ b provides a mechanism to incorporate segregation-mobility feedback effects into the governing equations. As noted in § 1, the different species of particle have different frictional properties, and for fingers to develop it is required that the larger particles experience greater resistance to motion. This is accounted for by taking a concentration-weighted sum (e.g. Pouliquen & Vallance 1999;Woodhouse et al. 2012), are the basal friction coefficients for smooth small and frictional large particles, respectively, and are written as functions of thickness and Froude number, It is assumed that the friction laws for the individual constituents are given by the dynamic friction law of Pouliquen & Forterre (2002), where N = S, L denotes small or large particles, respectively. The values µ N 1 = tan ζ N 1 and µ N 2 = tan ζ N 2 are constants, where angles ζ N 1 and ζ N 2 correspond to the minimum and maximum slope angles for which steady uniform flows are observed for a monodisperse material of constituent N . The length scales L N and dimensionless constants β N are found empirically, and may depend on both the granular material and bed composition. These constants are estimated for the laboratory set-up of figures 2-5, and are given in table 1, along with the other parameters that are kept constant in this paper.
Strictly speaking, the individual basal friction laws (3.7) only hold providing Fr > β N . For slower flows the extended law of Pouliquen & Forterre (2002) should be implemented, which accounts for arresting and static regions (see e.g. Johnson & Gray 2011;Edwards & Gray 2015). For simplicity it shall be assumed that (3.7) is valid everywhere for both types of particle. The implications of this assumption will be discussed in § § 6 and 7.
Segregation-induced finger formation in granular free-surface flows Material parameters that will remain constant throughout this paper.
The form of the final viscous term in the momentum equation (3.2) is motivated by the work done by  for monodisperse flows, who used the µ(I)-rheology (GDR MiDi 2004;da Cruz et al. 2005;Jop et al. 2005Jop et al. , 2006 to incorporate more of the specific material properties into the depth-averaged governing equations. They showed that, to leading order, the µ(I)-rheology only contributes via the basal friction coefficient, which is equivalent to (3.7). The resulting shallow-water-like equations are similar to those that have been successfully used in many granular flow configurations (Grigorian et al. 1967;Pouliquen 1999b;Gray et al. 2003). Higher-order viscous terms were introduced using the steady-state Bagnold velocity profile and lithostatic pressure distribution to derive an expression for the depth-averaged in-plane deviatoric stress, which  then wrote in the same form as in (3.2) using the relationship between the depth-averaged Bagnold velocity and flow thickness. In this formulation, νh 1/2 /2 may be interpreted as the kinematic viscosity, which acts, in the depth-integrated momentum balance equation, on the gradient term h∂ū/∂x.  were able to write the controlling coefficient ν = ν N explicitly in terms of the friction parameters of the monodisperse material as For the bidisperse flows being considered here it might be sensible to choose in an analogous manner to (3.4), where ν S and ν L are the coefficients for small and large particles and are given by (3.8). However, the coefficients ν S and ν L are only valid for slope angles ζ N 1 < ζ < ζ N 2 , where steady uniform flows are possible. Outside of this range the coefficient of viscosity is negative, and therefore the monodisperse depth-averaged theory is ill posed and must be regularised. This reflects the underlying ill posedness of the µ(I)-rheology (Barker et al. 2015). In order to get levee and finger formation the slope angle must be such that large particles in pure phase are brought to rest, whilst small particles and mixtures may still flow, i.e. (3.10) In this range the coefficient of viscosity for large particles is undefined, and it is not currently clear how to extend (3.8) to all slope angles. Instead of using (3.8) and (3.9), a constant bulk value ν > 0 is imposed in this paper, which may now be considered as a free parameter. The effect of changing this constant will be investigated and discussed. The large particle transport equation (2.15), together with the mass and momentum balances (3.1), (3.2), define a fully coupled system for the flow thickness and depth-averaged velocity and concentration. Segregation-mobility feedback effects are achieved through the effective basal friction in the momentum equation (3.2), with higher concentrations of large particles resulting in greater friction. From the monodisperse expressions it is known that the viscous terms are typically small in magnitude compared to the standard shallow-water contributions. The importance of these terms should not be underestimated, however, as they represent a singular perturbation to the inviscid equations (Woodhouse et al. 2012) which are ill posed at a critical Froude number. It will be shown here that the inclusion of viscosity is sufficient to regularise the equations.

Steady uniform flows
A simple solution to the system of equations (2.15), (3.1) and (3.2) is given by This represents a steady, fully developed flowing layer. Upon substitution into the governing equations, conservation of mass (3.1) and the large particle transport equation (2.15) are automatically satisfied.
Assuming there are no topography gradients, the momentum equation (3.2) reduces to a force balance between gravity and basal friction, is the steady uniform Froude number. Treating h 0 andφ 0 as known control parameters, equation (4.2) can be solved for F as a function of thickness and concentration. Substituting the friction law (3.4) and (3.7) into the force balance (4.2) leads to the quadratic equation where the coefficients are given by with M N = β N /L N . For a slope angle in the range given by (3.10), it can be seen (4.8) Consequently, the steady-state Froude number, found by taking the positive root of (4.4),  (4.9). The shaded regions represent whereφ 0 <φ * 0 (given by (4.8)), meaning there are too many frictional large particles for steady uniform flow.
is only positive providing thatφ 0 >φ * 0 , meaning steady uniform flow is not possible if there are too many frictional large particles. Figure 8 shows a contour plot of the two-parameter family of steady states F(h 0 ,φ 0 ), along with the regions whereφ 0 <φ * 0 . In the pure small limit (φ 0 = 1) the expression (4.9) reduces to that given in , which can also be derived from the more straightforward force balance, tan ζ = µ S b (h 0 , F). The corresponding steady uniform velocitiesū 0 (h 0 ,φ 0 ) may be recovered from the Froude number (4.9) using the relation (4.3). As a final point, the inclusion of higher-order terms into the momentum balance (3.2) does not change the steady-state values derived above, allowing direct comparisons to be made with the inviscid equations in subsequent sections.

Non-dimensionalisation
Assume the values h 0 andφ 0 are chosen such that a steady state h = h 0 ,φ =φ 0 ,ū = u 0 (h 0 ,φ 0 ) > 0 exists with corresponding Froude number F > 0, as described in the previous section. It is then convenient to introduce the scalings where the hats denote dimensionless quantities. Note that the depth-averaged concentrationφ is already non-dimensional. Using these scalings, the governing equations may be written as where the hats have been dropped for brevity and conservation of mass (5.2) has been used to simplify the momentum and large particle transport equations, (5.3) and (5.4) respectively. In (5.3) it has also implicitly been assumed thatū > 0. The basal friction coefficient is now written in terms of the non-dimensional variables as where the individual friction coefficients for each type of particle are given by and will typically take large values since the viscous terms are small in magnitude compared to the standard shallow-water contributions, as shown by . The Reynolds number is a function of both steady uniform flow thickness h 0 and concentrationφ 0 .

5.2.
Linearised equations and the characteristic polynomial By construction, there is a one-parameter family of steady-state solutions to (5.2)-(5.4) given by h = 1,ū = 1,φ =φ 0 , and so the variables are perturbed about this base state, The linearised governing equations then become

12)
Segregation-induced finger formation in granular free-surface flows where the simplified notation has been introduced to denote the steady-state values of the transport function and its derivative. The basal friction law contributes via terms where the positive constants Γ S , Γ L , are given by (5.18) for N = L, S. Now seek normal mode solutions of the form where, for temporal stability analysis, the wavenumber k is real and σ (k) = σ R (k) + iσ I (k) is complex. Substituting this ansatz into the linearised equations (5.11)-(5.13) allows the system to be written as an eigenvalue problem where W = (H, U, Φ) T and the matrix A is given by (5.21) Equation (5.20) has non-zero solutions for W if and only if |A − σ I| = 0, (where I is the 3 × 3 identity matrix) leading to a cubic characteristic polynomial where the coefficients f 0 , f 1 and f 2 are functions of the wavenumber k and base state, and are given in appendix A. Equation (5.22) can be solved to give three roots σ (1) , The growth rate σ M is then given by the maximum of these values, All of the cases considered in figures 9 and 10 show regions of positive growth rate for small k, meaning these base states are unstable to small wavenumber perturbations. However, there are major differences at larger values of k, depending on the Froude number and whether the viscous or inviscid equations are being considered. In the inviscid regime, the growth rates remain positive for all wavenumbers and are increasing functions of k, whereas the viscous growth rates take their maximum at a finite value of k and, in some cases, have a cutoff wavenumber above which perturbations are stable ( figure 9b).
An important result in Woodhouse et al. (2012) was the discovery of unbounded growth rates, σ M −→ ∞ as k −→ ∞ (figure 10), at a critical Froude number, meaning the inviscid system of equations was ill posed at this specific point in parameter space (Joseph & Saut 1990). Although the model was well posed almost everywhere, numerical simulations of the fingering instability showed that the width of the fingers was grid dependent. This was because the computations always encompassed a line of points where the depth-averaged concentration was equal to the critical concentration for ill posedness, and refining the computational domain introduced increasingly unstable small wavelength perturbations in these regions. The critical regime can be related to the characteristics of the inviscid equations, which are given in (x, t) space by the lines where λ is the characteristic wavespeed. For the non-dimensional mass and momentum balance equations (5.2), (5.3) there are two wavespeeds given by Note that it is only possible to separate the two sets of characteristics in this way because the segregation-mobility coupling arises through the source terms of the momentum equation, which do not contribute to the characteristic structure. Evaluating at the steady state (h,ū,φ) = (1, 1,φ 0 ), and noting that G 0 < 1, it can be seen that λ (1) > 0 and λ (3) > 0. The other wavespeed λ (2) is positive for F > 1 and negative for F < 1. Most importantly, at the critical Froude number 27) the large particle transport equation's characteristic wavespeed (5.26) coincides with one of those from the shallow water equations (5.25), depending on the sign of G 0 , and the system loses strict hyperbolicity. This is the point at which unbounded growth rates are found in the linear stability analysis for the inviscid model. It is therefore vital to check that the growth rate for this new viscous model remains bounded for all wavenumbers by conducting an asymptotic analysis of the characteristic polynomial (5.22) for k 1. The reader is referred to appendix A for the full details, but a summary of the key results is given in the following sections.

Inviscid high wavenumber asymptotics
For completeness, we begin by considering the inviscid case, ν = 0. The real part of the three roots is found to behave like (5.30) in the high wavenumber limit (see § A.1), meaning that the growth rate σ M tends to a constant value as k −→ ∞. From the definition (5.17) and the assumption (3.5) that large particles experience greater frictional forces, it follows that µφ < 0. Using this, and the fact that G 0 < 1 by (5.14), it can be seen that that σ (3) R > 0 for F < F c , where the critical Froude number is defined by (5.27). The growth rate σ M is therefore always positive for F < F c , meaning perturbations grow in time and the base state is unstable as k −→ ∞ (figure 9a). For F > F c the third root σ (3) R is now stable in the asymptotic limit and one must investigate the sign of the other roots σ (1) R and σ (2) R . For the parameters given in table 1 at least one of (5.28) or (5.29) is positive for all regions of (h 0 ,φ 0 ) space, meaning that σ M > 0 and the base state is again unstable (figure 9b).
At the critical Froude number F = F c the third root σ (3) R and either σ (1) R or σ (2) R become infinite, depending on the sign of G 0 . In these singular cases, the leading-order growth rate is instead given by the distinguished limit Segregation-induced finger formation in granular free-surface flows 189 for large k. The positive root in (5.31) will give unbounded growth rates proportional to k 1/2 (figure 10), which is the same as the one-dimensional result found in Woodhouse et al. (2012) and means that the inviscid model is ill posed in the sense of Joseph & Saut (1990).
5.4. Viscous high wavenumber asymptotics Returning to the viscous equations (finite R > 0), the leading-order behaviour of the three roots are for k 1 (see § A.2). In this high wavenumber limit, the first two of these are stable for all Froude numbers, whereas the third root is stable for F > F c and unstable for F < F c . To see this, note that G 0 < 1, µφ < 0, and F 2 (G 0 ) 2 may be written as (F/F c ) 2 . Whilst the sign is consistent with the inviscid case (5.30), the third root (5.34) does not tend to a constant value as k −→ ∞, and instead decays to zero like 1/k 2 (figure 9a). This is significant because it means that the growth rate σ M (k) will take its maximum value σ Max = max (σ M ) at a finite wavenumber k = k M , corresponding to the most unstable mode. Furthermore, for finite values of R > 0 the base state is always stable in the high wavenumber limit when F > F c . This leads to a cutoff wavenumber k = k c , and all perturbations with k > k c are stable (figures 9b and 11b).
Neither of these features are present in the inviscid regime, where the maximum growth rates are given by the asymptotic limits (5.28), (5.29) or (5.30). It is possible to derive analytical expressions for the cutoff wavenumber k c , details of which are given in § A.3. At the critical Froude number the expression (5.34) reduces to zero (since F 2 (G 0 ) 2 = 1) and a new asymptotic analysis leads to in the high wavenumber limit. The expression (5.35) is strictly positive, meaning the root is unstable for k 1. Intuitively, this result agrees with the inviscid analysis for the critical regime, since these are in some sense the most unstable cases. However, the growth rate crucially remains bounded for all values of k when viscous terms are included, and in fact decays like 1/k 4 as k −→ ∞ ( figure 10). This ensures that the model is well-posed, even in the previously problematic critical regime. These important results are highlighted in figure 11(a), which shows the maximum growth rates σ Max as a function of steady uniform Froude number F. For the inviscid equations σ Max −→ ∞ as F −→ F c , whereas the viscous curves remain bounded for all Froude numbers. As already noted, the linear stability calculations shown on figures 9 and 10 are computed using a fixed coefficient in the effective viscosity ν = 0.001 m 3/2 s −1 , and qualitatively similar behaviour is found for all positive choices of ν. The specific value does not affect whether the equations are well posed, but it does have an influence The cutoff wavenumber k c as a function of F, which only exists for the viscous regime, providing that F > F c . Larger values of the coefficient ν lead to smaller σ Max and k c . In both plots the depth-averaged concentration is fixed atφ 0 = 0.8 and F is varied by changing the steady-state thickness h 0 . on the stability of steady uniform flows. Increasing ν, and consequently reducing the Reynolds number R, leads to lower maximum growth rates σ Max (figure 11a) and therefore has a stabilising effect on the base state. The most unstable mode k M and cutoff wavenumber k c (when it exists) decrease with increasing ν (figure 11b).

Two-dimensional numerical simulations for a propagating front
Having eliminated the possibility of unbounded growth rates as a source of ill posedness, the capability of the new governing equations to model physical systems may be tested. One such system is segregation-induced fingering instabilities that develop as a front of bidisperse material propagates down an inclined plane (see figures 2, 3, movie 1 and Pouliquen et al. 1997;Pouliquen & Vallance 1999;Woodhouse et al. 2012, for example). In an initially homogeneous mixture large particles are rapidly segregated to the surface, where the higher velocity shears them to the front. Here they are overrun, resegregated upwards and recirculated by the bulk to form a coarse-rich flow head. This head experiences enhanced frictional forces and may break down into a series of 'finger-like' structures. Each finger is bounded by static or slowly moving levees of coarse material that channelise the finer, more mobile interior. This phenomenon is important because it is an example of a segregation-mobility feedback effect, with the instability being suppressed in experiments using monodisperse material (see figure 5) Segregation-induced finger formation in granular free-surface flows 191 able to produce numerical simulations showing spontaneous leveed finger formation and elongation, but their results were grid dependent due to the ill posedness discussed in the previous sections.

Generalised equations
Many granular flow configurations have little variation in the lateral direction, allowing them to be modelled by one-dimensional (depth-averaged) equations governing the evolution in time and downslope direction (e.g. Edwards & Gray 2015). However, the instability that occurs as a bidisperse mixture propagates down a plane is predominantly in the transverse direction, meaning that the lateral dependence needs to be explicitly accounted for in the model. The governing equations must therefore be extended to two spatial dimensions. Baker et al. (2016) examined a similar problem for monodisperse material and generalised the one-dimensional depth-averaged µ(I)-rheology of  to two dimensions. Using their work as a base and introducing a cross-slope coordinate y and depth-averaged velocityv, the two-dimensional segregation-mobility feedback equations (reverting to dimensional variables) are is the depth-averaged velocity vector, ∇ = (∂/∂x, ∂/∂y) denotes the gradient operator in (x, y) space, '·' is the dot product and ⊗ the dyadic product. Assuming no basal topography, the source terms S = (S x , S y ) are given by S x = cos ζ tan ζ − µ bū |ū| , S y = − cos ζ µ bv |ū| , (6.4a,b) where |ū| = (ū 2 +v 2 ) 1/2 is the magnitude of the velocity. The basal friction law µ b (h, Fr,φ) is defined by (3.4) and (3.7) using the more general Froude number definition Fr = |ū| √ gh cos ζ . (6.5) The coefficient ν > 0 remains a free parameter for slope angles in the range (3.10), and the depth-integrated strain-rate tensor is D = 1 2 (∇ū + (∇ū) T ). (6.6) The generalised large particle transport equation ( profile is taken to be the same (see Woodhouse et al. 2012). The transport function G then remains unchanged and is given by (2.17) with α = 1/7. Note that, for simplicity, the linear stability analysis conducted in § 5 is for the one-dimensional equations and does not strictly apply to their generalised form (6.1)-(6.3). One cannot guarantee that these will also be well posed without further calculations. Woodhouse et al. (2012) carried out a full two-dimensional stability analysis for their inviscid equations and found that the most unstable mode was usually in the downslope direction. In particular, the unbounded growth rates were only apparent in the large k x (downslope wavenumber) limit. The equivalent analysis has been carried out for the viscous equations (see appendix B), and unstable modes are only found for non-zero downslope perturbations. In all cases, the growth rates remain bounded for all wavenumbers, meaning the generalised equations are well posed.

Numerical solutions
The system of coupled partial differential equations (PDEs) (6.1)-(6.3) is now solved numerically using the shock-capturing central scheme of Kurganov & Tadmor (2000). This is second order in space and has a semi-discrete formulation, allowing it to be combined with a time stepper of choice. In this case a Runge-Kutta-Chebyshev adaptive step method (Medovikov 1998) is employed. The scheme is well suited to this particular problem because it is easily generalised to multiple spatial dimensions and is capable of handling convection-diffusion equations. It has previously been used to solve similar systems governing granular flows, for example a two-dimensional breaking size-segregation wave , granular roll waves Razis et al. 2014) and erosion-deposition waves (Edwards & Gray 2015). To calculate the numerical fluxes, the scheme requires the specification of a flux limiter. Here, the weighted essentially non-oscillatory (WENO) limiter detailed in Noelle (2000) is chosen. In order to utilise the numerical method of Kurganov & Tadmor (2000), the governing equations must be written in conservative form. This introduces numerical singularities in both the convective and diffusive fluxes as h −→ 0. To get around this potential problem a minimum flow thickness smaller than one particle diameter, h min = 10 −4 m, is introduced and the fluxes are set to zero whenever h < h min .
For numerical simulations of a front of granular material advancing down an inclined plane, initial conditions of an empty slope h(x, y, 0) = 0,ū(x, y, 0) = 0,v(x, y, 0) = 0,φ(x, y, 0) = 0, (6.7a−d) are prescribed. Simulations are carried out on a computational domain L x × L y discretized over N x × N y grid cells, giving a spatial resolution of N x /L x × N y /L y points per metre. Periodic boundary conditions are specified in the y direction, and at the upstream boundary (x = 0) the inflow conditions h(0, y, t) = (1 − e −10t )h 0 ,ū(0, y, t) =ū 0 ,v(0, y, t) = 0,φ(0, y, t) =φ 0 + δφ p (y), (6.8a−d) are enforced. These correspond to the steady uniform flow of § 4, with the thickness prefactor (1 − e −10t ) used to ensure a smooth transition from an empty slope. Note that the simulations represent a continuous inflow (akin to figures 3a, 4a) as opposed Segregation-induced finger formation in granular free-surface flows 193 to a finite release of material (figures 3b, 4b). To introduce transverse variation into the system the inflow concentration is perturbed with a periodic function φ p (y) = sin (2πy/L y ), (6.9) of magnitude δ = 8 × 10 −5 . Simulations are halted before material reaches the downstream boundary and so, since the equations degenerate when h = 0, no boundary conditions are required at this end of the domain. Figures 12-14 and supplementary movie 3 show the resulting contours of the flow height h, depth-averaged concentrationφ and speed |ū|. It can be seen that the governing equations are able to qualitatively reproduce many of the key phenomena present in the experimental set-up. The concentration plots ( figure 12) show that a large-particle-rich front quickly develops and subsequently breaks up into a series of fingers. These are bounded by coarse-grained lateral levees, and elongate as time progresses. As the fingers emerge from the uniform front the large particles are shouldered aside into slow moving coarse-particle-rich levees, and the large particle regions at the finger tips are dramatically reduced in size. In practice a breaking size-segregation wave will form just behind the front (Thornton & Gray 2008;Gray & Ancey 2009;Johnson et al. 2012;Gajjar et al. 2016) which in the depth-averaged segregation theory (Gray & Kokelaar 2010a,b) is replaced by a depth-averaged concentration shock. The existence of a breaking size-segregation wave in the physical experiments, as well as diffusion, makes the transition at the large-rich front much more diffuse than predicted in the simulations.
Note that the coarse-rich front that forms before the onset of finger formation is not a steady travelling wave solution to the system of equations, because large particles continue to accrue at the flow front (Gray & Kokelaar 2010a). A steadily travelling base state must include a mechanism to counteract this accumulation, which Gray & Ancey (2009) achieved by depositing the coarser grains arriving at the front onto the underlying substrate. In the absence of basal deposition the large particles may be removed from the frontal region through lateral advection towards the levees, and Johnson et al. (2012) showed that this gave rise to a steady travelling segregation profile on the flow centreline, for a specified flow thickness and velocity field. We now believe that a single leveed finger can propagate steadily downslope, and anticipate that such a travelling wave solution to equations (6.1)-(6.3) would be an appropriate base state, although this has not been confirmed. This behaviour is in direct contrast to the pure small limit (φ = 1). Assuming no transverse variation, the monodisperse equations admit a steady travelling front solution , which can be found analytically (Gray & Ancey 2009) or numerically (Pouliquen 1999b) and is unaffected by the new viscous terms. Two-dimensional computations suggest that this base state does not break down into fingers, which is consistent with the monodisperse experiments (figure 5), although roll waves may form and travel to the flow front.
From figure 13 it can be seen that once the fingers form the material travelling down the centre of the channels is moving approximately twice as fast as the steady uniform flow behind. This can be explained by examining the coarse-rich lateral levees, where the flow approaches zero velocity at the margins between adjacent fingers. To account for these reduced flux sections, the central regions must increase their flux, through a combination of flow thickening and increased speed, in order to conserve material from the constant inflow. There is an important discrepancy in the slow moving zones between experiments and numerics. In the laboratory, some of the material in the levees is completely stationary, but this is not achieved at any merging and coarsening of fingers in the final stages of flow, the extended friction law of Pouliquen & Forterre (2002) needs to be generalised to bidisperse material and implemented. The flow thickness plots ( figure 14) show that the region immediately behind the flow front is elevated above the steady uniform inflow depth. This is consistent with observations of bulbous flow fronts in both experiments and geophysical events Kokelaar et al. 2014). There are also regions at the rear of the fingers that are significantly higher than the rest of the material. Comparing with figures 12 and 13, these correspond to the large particle islands that are also seen in small-scale experiments. For the continuous inflow regime simulated here, they move slowly downstream, which matches the experimental erosion processes (figure 3a).
Whilst there are no quantitative comparisons made at this stage, typical time scales for the onset of finger formation are roughly consistent. The channel widths do not correspond to the imposed inflow perturbation, meaning they are set by the system itself, and are of the correct order of magnitude (around 5 cm for the simulations shown). Further discussion about the finger characteristics follows in the next section. 6.3. Finger characteristics: numerical versus physical viscosity As well as showing that the new equations remove the possibility of unbounded growth rates, it also needs to be checked that the resulting simulations are grid convergent. Both the experiments and the non-linear equations exhibit high sensitivity to the initial conditions, which is highlighted numerically by the apparent random nature of the resulting fingers. Even though these arise through the integration of deterministic PDEs, numerical round-off error, subsequently magnified by the underlying instability of the equations, is sufficient to break the symmetry of the inflow (6.8) and (6.9). Changing the grid resolution will change these conditions, and one can therefore not expect to obtain identical results when running simulations with different numbers of grid points. Attention is instead given to the wavelength of fingers (cross-slope distance between two adjacent channels) and their elongation (downslope distance between front and rear of leveed region). Woodhouse et al. (2012) showed that, for the inviscid equations, the wavelength got progressively smaller and the fingers more elongated with decreasing mesh size, suggesting that numerical viscosity was a controlling factor for channel characteristics.
Computations are carried out for different grid resolutions and domain widths L y , and example results are shown on figure 15. The depth-averaged concentrationφ is plotted at time t = 7.5 s. As expected, the results are not perfectly identical for varying numbers of grid points, with the position and shape of the fingers changing slightly between runs. This is actually a desirable property, since no two experiments are identical and so it is good that the numerics exhibit a similar degree of sensitivity. However, like the experiments, the numerical finger width and their elongation stays approximately the same, which is in direct contrast to Woodhouse et al. (2012). At a given time it is straightforward to calculate the mean finger wavelength Λ by dividing the domain width L y by the number of fingers, and this is shown on figure 16(a) for the different computations. There is some variation when too few cells are used, but the wavelength converges at sufficiently high resolutions and is independent of the domain width. Note that counting the number of fingers from plots such as figure 15 is fairly intuitive since the individual channels are well defined at this time, t = 7.5 s. However, running the simulations for longer can lead to merging events, making precise definitions more difficult and meaning that Λ will also evolve over time.
To estimate the elongation of the fingers, a maximum and minimum front position, x + f and x − f respectively, are defined for each time as x + f (t) = max {x : h(x, y, t) < h min }, (6.10) x − f (t) = min {x : |φ(x, y, t) −φ 0 | >φ 0 /2}, (6.11) where (6.11) is chosen to capture the sharp transition between the mixed inflow and pure large region at the back of the fingers. The lines where x = x + f (t) and x = x − f (t) are shown for reference on figure 15. These values can be used to calculate the length, or elongation, of the fingers as (6.12) Figure 16(b) shows the variation of x l with changing numbers of grid points at different times t. Again, there is slight variation at low resolutions but it eventually converges to roughly the same value at all times. This confirms that the numerical viscosity is no longer controlling the finger properties, which is major progress in modelling spontaneous finger formation. The wavelength and downslope extent of the fingers is now being set by the equations themselves, in particular the newly introduced physical viscosity. All the computations thus far were carried out with a coefficient in the effective viscosity ν = 0.01 m 3/2 s −1 . This may be unphysically high, but was chosen to ensure that it would outweigh any numerical diffusion in checking grid convergence. Since it is being treated as a free parameter in this paper, simulations were also conducted with values of ν ranging from 0.005 to 0.05 to investigate what effect this has on the Segregation-induced finger formation in granular free-surface flows  physical characteristics of the fingers. The results are shown on figure 17, where it can be seen that increasing ν typically leads to fewer fingers (larger wavelength) that are less elongated (smaller x l ). One interpretation is that higher material viscosities suppress finger formation. This is investigated further on figure 18, where the front positions x − f , x + f and finger lengths x l are plotted as functions of time t for different viscosities. During the initial uniform propagation, both the minimum and maximum front positions move at constant velocities, with x + f travelling faster than x − f leading to a steady growth in x l . At the onset of finger formation, there is an increase in the speed of x + f as the more mobile interior breaks through the resistive front. Interestingly, the maximum extent of the flow then travels faster than the prescribed inflow, which could have important implications from a hazards mitigation perspective. At the same time there is a corresponding deceleration of x − f as the slow moving large particle islands form, and combined with the acceleration of the finger tips this leads 200 J. L. Baker, C. G. Johnson and J. M. N. T. Gray to a sudden increase in the rate of elongation of the frontal region. This transition point happens at earlier times for lower values of ν, which is consistent with the linear stability results in § 5 where higher viscosities lead to slower maximum growth rates. It is also consistent with the idea that the frontal break down is closely tied to the region of steady uniform flow immediately behind.

Conclusions
Polydispersity in granular flows can significantly alter the overall flow characteristics, with particle size segregation leading to the formation of regions with different sized grains. If these grains also have different frictional properties this then feeds back on the bulk flow and can produce rich behaviour that is not seen in monodisperse flows, such as self-organisation into coarse-grained levees and segregation-induced fingering instabilities. This paper has presented a fully coupled depth-averaged model for such segregation-mobility feedback effects. The evolving particle distribution is governed by a large particle transport equation, which is derived from a full segregation equation by assuming perfect vertical inverse grading and a shear profile through the avalanche. It is shown that using a physically motivated Bagnold profile leads to qualitatively similar transport functions to previous models (Gray & Kokelaar 2010a,b;Woodhouse et al. 2012), providing the correct shear parameter is chosen. This Bagnold assumption is also consistent with that used by  to incorporate higher-order terms into their depth-averaged µ(I)-rheology for monodisperse flows. These second-order terms are generalised for the bidisperse regime considered here, giving a viscous model with coupling achieved via a concentration-weighted effective basal friction.
Linear stability calculations of the steady uniform base state highlight the significance of the higher-order terms, which were not present in Woodhouse et al. (2012). The growth rates now remain bounded everywhere, whereas for the inviscid equations there is a critical Froude number that gives rise to unbounded growth rates in the high wavenumber asymptotic limit. This means that the introduction of viscosity regularises the problem and ensures that the governing equations are mathematically well posed. Perhaps this is not such a surprising result, but the advantage of the particular viscosity formulation used in this paper is that the structure of the higher-order terms is physically motivated, with the parameters completely determined by the µ(I)-rheology and no additional fitting parameters required, at least in the monodisperse regime.
The system is then generalised to two spatial dimensions in an analogous way to Baker et al. (2016) (for monodisperse flows) and numerical simulations of a propagating flow front are presented. The equations are able to predict the formation of a coarse-rich flow front that develops instabilities and breaks up into a series of finger-like structures, elongating over time. The lateral boundaries of these fingers consist of large particle levees, where material is travelling significantly slower than the finer-grained interior, which itself moves faster than the steady uniform supply. There is also a noticeable speed-up at the tip of the flow as the fingers first emerge. At the rear of the levees are clusters of coarse grains that migrate slowly downslope, consistent with the erosion processes observed during continuous inflow experiments. Unlike the physical system, there are no positions in the simulation domain where the velocities reach precisely zero. This is not problematic for the onset of finger formation considered in this paper, but is expected to become more significant when examining the final run-out of flows after the supply of particles has been stopped.
Segregation-induced finger formation in granular free-surface flows

201
In this case the static regions are important in preventing lateral levee spreading and forming grain-free regions between distinct fingers. To bring material to rest and lock in the structure of the channel walls the extended friction law of Pouliquen & Forterre (2002) appears to be necessary.
The ill posedness of the inviscid model (Woodhouse et al. 2012) manifested itself as grid-dependent simulations, with increasing numbers of computational cells leading to larger numbers of fingers that spread over greater downslope distances. This suggests that the finger characteristics were being determined by numerical viscosity. Solutions of the new equations do not suffer from the same problem, with both the finger wavelength and their elongation converging for sufficiently high resolutions. The sensitivity of the system means that there is still some variety in the exact shape and position of the channels, but this is not necessarily undesirable since the laboratory experiments are also highly sensitive to the initial conditions. The channels in consecutive experimental runs will always form in slightly different places, although the characteristic width and time scales will be similar.
This paper is only a preliminary investigation into finger formation, focusing on establishing the well posedness and grid convergence of the system as opposed to detailed analysis of finger growth rates and wavelengths. Indeed, such analysis is complicated by a number of factors, including the sensitivity described above, the temporal evolution of the base state (which is uniform in y but evolving in x and t), nonlinear coarsening, and the interactions between instabilities of the front and instabilities of the steady uniform flow behind. The linear stability analysis of § 5 is presented in one spatial dimension to simplify and aid visualisation, and, whilst the two-dimensional version (appendix B) of a uniform base state is important for checking well posedness in the general case, it is difficult to relate this to the numerical computations of a non-uniform base state.
Despite these difficulties, the qualitative properties of the fingers in the simulations are now controlled by the physical viscosity via the effective coefficient ν. Increasing this leads to larger wavelength fingers that do not spread as far in the downslope direction. Higher viscosity flows also take longer to break down into leveed channels, which is consistent with linear stability calculations of the steady uniform base state. The value of ν is currently treated as a free parameter, since fingers only form at slope angles where the depth-averaged µ(I)-rheology of  needs to be regularised. Choosing ν so that the numerical simulations match experimental data may provide insight into precisely how to achieve this regularisation. It may also be possible to calibrate the coefficient of viscosity from other bidisperse experiments, for example using the cutoff frequency for the roll-wave instability. This appendix provides derivations of the asymptotic results presented in § § 5.3 and 5.4, as well as analytical expressions for the cutoff wavenumber for instability in the viscous regime. Firstly, the characteristic polynomial (5.22) is given again for completeness, It is useful to expand the coefficients in powers of the wavenumber k, A.1. Inviscid asymptotics Firstly consider the inviscid case, so that f 0,4 = f 1,3 = f 2,2 = 0. An asymptotic expansion is sought of the form σ ∼ σ 0 k p + σ 1 k q + σ 2 k r + σ 3 k s + · · · (A 8) for k 1, where the exponents p > q > r > s · · · are to be determined. The terms σ 0 , σ 1 , σ 2 , σ 3 , . . . are all strictly O(1) and are calculated in increasing order until a non-zero real part is found. This then determines the leading-order growth rate of the root. Substituting the ansatz (A 8) into the characteristic polynomial (A 1) gives the dominant balance p = 1 and the leading-order behaviour at O(k 3 ) is with associated solutions These are all purely imaginary, so the next order in the expansion is considered by setting q = 0 which gives, at O(k 2 ), the linear equation, Segregation-induced finger formation in granular free-surface flows

203
Rearranging for σ 1 and substituting in the expressions (A 10) gives the real values The growth rate of all three roots therefore tends to a constant as k −→ ∞, with the value being determined by (A 12)-(A 14). At the critical Froude number F = F c = 1/|G 0 |, two of the above roots become infinite because the coefficient of σ 1 in (A 11) is zero (meaning the denominator in (A 12)-(A 14) degenerates). In this case the alternative dominant balance q = 1/2 is chosen, which gives at O(k 2 ), Substituting in σ 0 = −i(1 − G 0 ), which is the same for both critical roots, gives and the corresponding real parts In this critical regime these two roots scale with k 1/2 in the large wavenumber limit, and the positive choice in (A 17) grows unboundedly. By definition (Joseph & Saut 1990), this means that the inviscid equations are ill posed at the critical Froude number.
A.2. Viscous asymptotics Repeating the same process for the viscous equations, two possible dominant balances are found for the leading-order behaviour. Choosing p = 2 gives, at O(k 6 ), which has non-zero real solution The growth rate of the first root is therefore negative and decays like k 2 for large k.
The other roots are found by letting p = 1, giving leading-order behaviour at O(k 4 ) determined by the quadratic which has purely imaginary solutions The next-order dominant balance is achieved by setting q = 0, which gives, at O(k 3 ), the linear equation Substituting in for σ 0 = σ (2) 0 from (A 21) gives the real solution Hence, the second root is also stable in the large wavenumber limit, with growth rates tending to the negative constant (A 23). However, using σ 0 = σ (3) 0 gives σ 1 = 0. In this case the alternative scaling q = −1 is chosen and the next-order behaviour is at O(k 2 ) and given by which leads to the imaginary solution The next order is at O(k) and achieved by setting r = −2. This gives the linear equation with real solution The third growth rate therefore decays to zero like O(k −2 ) for k 1. It may be stable or unstable, depending on the sign of (A 27), but remains bounded for all wavenumbers. Note that (A 27) degenerates to zero at F = F c . In this case one must instead choose r = −3, giving at O(1), (if 1,1 + 2f 2,0 σ 0 + f 2,2 σ 1 )σ 1 + (if 1,3 + 2f 2,2 σ 0 )σ 2 = 0, (A 28) with associated imaginary solution The next highest order is achieved by setting s = −4, which gives at O(k −1 ), ( f 1,2 σ 2 + 2if 2,1 σ 0 σ 2 + if 2,1 σ 2 1 + 3σ 0 σ 2 1 + 3σ 2 0 σ 2 ) + (if 1,3 + 2f 2,2 σ 0 )σ 3 = 0, (A 30) with real solution This is positive, meaning the third root will be unstable at the critical Froude number in the large wavenumber limit. However, the choice s = −4 means that it will decay to zero according to k −4 , crucially remaining bounded for all values of k. The viscous equations therefore remain well posed.
Segregation-induced finger formation in granular free-surface flows 205 A.3. Viscous cutoff wavenumber The cutoff wavenumber for instability, k = k c , occurs when the growth rate is purely imaginary, say σ = iσ c . Substituting into the characteristic polynomial (A 1) and equating real and imaginary parts gives Equation (A 32) is a quadratic in σ c which can be solved to give the relation where the functions C ± (k c ) are defined as . (A 35) Substituting (A 34) into the imaginary parts equation (A 33) and neglecting the trivial root k c = 0 gives the cubic equation with associated real solutions Now, equation (A 35) can be rearranged, squared and factorised to give Two solutions of (A 38) for k c are given by which are independent of the value of C. However, since µū > 0 and the cutoff wavenumber must be real these can immediately be discarded. The other roots are given by When C = C 3 the denominator in (A 40) is zero, meaning there are no additional roots in this case. Substituting in the other values of C given by (A 37) leads to the solutions Written out explicitly for reference, these coefficients are f 1,4,0 = 1 2F 2 R (1 + 3F 2 (2G 0 − 3)), f 1,5,0 = (2 − G 0 ) 2R 2 , f 2,0,0 = tan ζ F 4 µū, f 2,0,2 = 1 2F 2 R (tan ζ + 2(R + µū)), f 2,0,4 = 1 2R 2 , f 2,1,0 = 1 F 2 ((3 − G 0 )(µū + tan ζ ) − µ h + G 0 µφ), f 2,1,2 = 3(3 − G 0 ) 2R , f 2,2,0 = 1 2F 2 R (µū + 2 tan ζ + 2R(1 + 3F 2 (G 0 − 2))), f 2,2,2 = 1 R 2 , B.1. High downslope wavenumber asymptotics First consider the high downslope wavenumber limit, k x 1, for a fixed cross-slope wavenumber k y 0. This closely resembles the one-dimensional analysis conducted in § A.2, and is a viscous analogue of the results presented in Woodhouse et al. (2012). An asymptotic expansion is sought of the form σ ∼ σ 0 k p x + σ 1 k q x + σ 2 k r x + σ 3 k s x + · · · (B 6) for k x 1 and exponents p > q > r > s > · · · to be determined, alongside the strictly O(1) values σ 0 , σ 1 , σ 2 , σ 3 , . . . Omitting the details for simplicity, the leading-order 208 J. L. Baker, C. G. Johnson and J. M. N. T. Gray growth rates of the four roots are found to be for k x 1. Note the similarities to the one-dimensional case, with the asymptotic behaviour being independent of the transverse wavenumber k y and remaining bounded for all k x . The extra spatial dimension simply introduces an extra root, here referred to as σ (2) R , which is stable in the asymptotic limit. At the critical Froude number σ (2) R in (B 7) degenerates to zero, and the leading-order growth rate in this case is instead given by for k x 1. This is always positive, meaning that the root is unstable, but remains bounded even in this critical regime, decaying like 1/k 4 x as k x −→ ∞. Consequently, when considering the asymptotic behaviour (B 7), (B 8) of all of the roots, it can be seen that the maximum growth rate σ is bounded above for all parameter values.
The real part of σ (1) and σ (±) is negative for all values of k y , meaning these roots are stable. Consequently, for k x = 0, the maximum growth rate is given by σ M ≡ 0 and hence the base state is neutrally stable. One requires non-zero downslope perturbations for disturbances to grow in time, so the subsequent calculations assume k x > 0 and seek an asymptotic expansion of the form σ ∼ σ 0 k p y + σ 1 k q y + · · · . (B 10) Following the same procedure as in appendix A, the leading-order growth rates of the four roots is given by for k y 1. The first three of these are stable and analogous to the downslope results (B 7), whereas the fourth root is unstable for non-zero k x in this asymptotic limit. For a fixed downslope wavenumber, the growth rate remains bounded as k y −→ ∞ (decaying to zero) and hence so too does the maximum growth rate σ M .
B.3. High wavenumber asymptotics In order to completely establish the well posedness of the system, it also needs to be checked that the growth rates remain bounded when both the downslope and crossslope wavenumbers tend to infinity, i.e. k x 1 and k y 1. The analysis conducted in § § B.1 and B.2 ensures it is sufficient to consider the curves k x = k, k y = mk a , (B 12a,b) for k 1, where the constant m > 0 and exponent a > 0. In this case an asymptotic expansion is sought of the form σ ∼ σ 0 k p + σ 1 k q + σ 2 k r + · · · , (B 13) where σ 0 , σ 1 , σ 2 , . . . and the parameters p > q > r > · · · may now depend on a. For 0 < a < 1 the leading-order growth rates of the first three roots is given by for k 1, which are all stable. The growth rate of the final root is more complicated to calculate, but note that the first two terms in the asymptotic expansion (B 13) are given by the imaginary expressions for k 1, and hence the exponent of k in the leading-order real growth rate must therefore be less than −1. This growth rate will therefore decay to zero in the large wavenumber asymptotic limit, remaining bounded even if the root is unstable. Similar behaviour is found for other values of the parameter a, with the first three roots for a = 1 having leading-order growth rates for k 1, and the expansion of the fourth root begins with the imaginary terms For the final case a > 1, it can be shown that the roots have growth rates for k 1, alongside the asymptotic behaviour σ (4) ∼ −i(1 − G 0 )k + iRG 0 (1 − G 0 )µφ F 2 G 0 m 2 1 k 2a−1 + · · · (B 19) for the fourth root, and hence the real growth must necessarily remain bounded for all wavenumbers. This concludes the two-dimensional linear stability analysis, where all Downloaded from https://www.cambridge.org/core. IP address: 54.70.40.11, on possible variants of the asymptotic limit |k| = (k 2 x + k 2 y ) 1/2 −→ ∞ have been covered. Whilst the full algebraic expressions have sometimes been omitted for simplicity, it has been shown that, in all cases, the growth rates remain bounded and therefore that the two-dimensional governing equations are well posed.
All of the above stability analysis has been validated numerically, using timedependent solutions of the fully nonlinear equations as described in § 6.2.