Segregation-induced finger formation in granular free-surface flows

Abstract

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.

1 Introduction

The process of particle size segregation, whereby mixtures of different sized particles separate into distinct grain-size classes during flow, can be very pronounced, 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; Gray & Ancey 2011; Marks, Rognon & Einav 2012; Gray & Ancey 2015) 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 (Johnson et al. 2012) 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).

Figure 1. Pyroclastic flow deposits from the eruption of Mount St Helens on July 22nd 1980 showing evidence of particle size segregation and finger formation during runout (Photo courtesy Dan Miller and USGS).

Experiments are carried out using a bidisperse mixture of spherical ballotini (white, $75{-}150~\unicode[STIX]{x03BC}\text{m}$ diameter) and angular carborundum (brown, $305{-}355~\unicode[STIX]{x03BC}\text{m}$ ) flowing down a plane inclined at $27^{\circ }$ , which is roughened by attaching a single layer of turquoise ballotini $(750{-}1000~\unicode[STIX]{x03BC}\text{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~\text{cm}~\text{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~\text{mm}~\text{s}^{-1}$ , compared to average bulk velocities of $30~\text{mm}~\text{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; Gray et al. 2015) 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 3 b and 4 b). This final phase can lead to unexpectedly long run-out distances, such as in figure 1, before eventually coming to rest and revealing the lubricating fine-grained levee lining (Kokelaar et al. 2014).

Figure 2. Experiments on a plane inclined at $27^{\circ }$ using 80 % ballotini (white, $75{-}150~\unicode[STIX]{x03BC}\text{m}$ ), 20 % carborundum (brown, $305{-}355~\unicode[STIX]{x03BC}\text{m}$ ) released from rest through a double gate system of inflow thickness 2 mm. The chute is roughened with turquoise ballotini $(750{-}1000~\unicode[STIX]{x03BC}\text{m})$ . Images show snapshots at approximate times $t=0.9~\text{s}$ , $t=2.6~\text{s}$ , $t=4.1~\text{s}$ , $t=6.0~\text{s}$ and $t=7.9~\text{s}$ . Supplementary movie 1 available online.

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~\unicode[STIX]{x03BC}\text{m})$ , 20 % brown carborundum $(305{-}355~\unicode[STIX]{x03BC}\text{m})$ is used and the inflow thickness is 2 mm.

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.

Figure 5. Experiments on a plane inclined at $27^{\circ }$ using monodisperse granular material consisting of 100 % ballotini $(75{-}150~\unicode[STIX]{x03BC}\text{m})$ released from rest through a double gate system of inflow thickness 2 mm. Images show snapshots at approximate times $t=0.4~\text{s}$ , $t=1.7~\text{s}$ , $t=3.0~\text{s}$ , $t=4.3~\text{s}$ and $t=5.7~\text{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.

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 at this slope inclination of 27 $^{\circ }$ 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 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 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 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. Gray & Edwards (2014) recently devised a strategy to achieve this using the $\unicode[STIX]{x1D707}(I)$ -rheology for dense granular flows (GDR MiDi 2004; da Cruz et al. 2005; Jop, Forterre & Pouliquen 2005, 2006). 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 $\unicode[STIX]{x1D707}(I)$ -rheology is able to predict this cutoff for a wide range of Froude numbers and slope angles without any fitting parameters (Gray & Edwards 2014).

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 Gray & Edwards’ (2014) depth-averaged $\unicode[STIX]{x1D707}(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.

2 A depth-averaged model for particle size-segregation

Consider a Cartesian coordinate system $Oxz$ with the $x$ -axis pointing downslope at an angle $\unicode[STIX]{x1D701}$ 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 $\unicode[STIX]{x1D719}\in [0,1]$ (so that the proportion of large particles is $1-\unicode[STIX]{x1D719}$ ), 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; Gray & Ancey 2011; Gajjar & Gray 2014),

(2.1) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\unicode[STIX]{x1D719}u)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}(\unicode[STIX]{x1D719}w)-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}(Q(\unicode[STIX]{x1D719}))=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left(D\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}\right), & & \displaystyle\end{eqnarray}$$

where the bulk velocity $\boldsymbol{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(\unicode[STIX]{x1D719})\geqslant 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(\unicode[STIX]{x1D719})=q\unicode[STIX]{x1D719}(1-\unicode[STIX]{x1D719})$ , (Gray & Thornton 2005) or skewed cubic $Q(\unicode[STIX]{x1D719})=q\unicode[STIX]{x1D719}(1-\unicode[STIX]{x1D719})(1-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719})$ , (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.

Figure 6. A schematic diagram of the coordinate axes $Oxz$ inclined at an angle $\unicode[STIX]{x1D701}$ to the horizontal, so that the $x$ -axis points downslope and the $z$ -axis is the upward pointing normal. The granular material lies between the base $z=b(x)$ and free surface $z=s(x,t)$ , giving a flow thickness $h(x,t)=s-b$ . At $z=l(x,t)$ there is an interface separating a layer of pure small particles ( $\unicode[STIX]{x1D719}=1$ ) of thickness $\unicode[STIX]{x1D702}(x,t)=l-b$ at the bottom of the flow from a layer of pure large particles ( $\unicode[STIX]{x1D719}=0$ ) lying on top.

The segregation equation (2.1) is subject to kinematic boundary conditions,

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle u_{b}\frac{\unicode[STIX]{x2202}b}{\unicode[STIX]{x2202}x}-w_{b}=0,\quad \text{at }z=b(x), & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}t}+u_{s}\frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}x}-w_{s}=0,\quad \text{at }z=s(x,t), & \displaystyle\end{eqnarray}$$

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) $$\begin{eqnarray}\displaystyle Q(\unicode[STIX]{x1D719})+D\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}=0,\quad \text{at }z=b(x)\text{ and }z=s(x,t). & & \displaystyle\end{eqnarray}$$

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

(2.5) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h\bar{\unicode[STIX]{x1D719}})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(h\overline{\unicode[STIX]{x1D719}u})-\left[\unicode[STIX]{x1D719}\left(\frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}t}+u\frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}x}-w\right)\right]_{b}^{s}=\left[Q(\unicode[STIX]{x1D719})+D\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}\right]_{b}^{s}, & & \displaystyle\end{eqnarray}$$

where

(2.6a,b ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D719}}=\frac{1}{h}\int _{b}^{s}\unicode[STIX]{x1D719}\,\text{d}z,\quad \overline{\unicode[STIX]{x1D719}u}=\frac{1}{h}\int _{b}^{s}\unicode[STIX]{x1D719}u\,\text{d}z, & & \displaystyle\end{eqnarray}$$

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 square-bracketed terms disappear and the depth-integrated segregation equation (2.5) reduces to

(2.7) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h\bar{\unicode[STIX]{x1D719}})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(h\overline{\unicode[STIX]{x1D719}u})=0. & & \displaystyle\end{eqnarray}$$

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

(2.8) $$\begin{eqnarray}\displaystyle \bar{u}(x,t)=\frac{1}{h}\int _{b}^{s}u(x,z,t)\,\text{d}z. & & \displaystyle\end{eqnarray}$$

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

(2.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}=\left\{\begin{array}{@{}ll@{}}0,\quad & l<z<s,\\ 1,\quad & b<z<l,\end{array}\right. & & \displaystyle\end{eqnarray}$$

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

(2.10) $$\begin{eqnarray}\displaystyle u(x,z,t)=\bar{u}(x,t)f(\hat{z}), & & \displaystyle\end{eqnarray}$$

where $\hat{z}=(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) $$\begin{eqnarray}\displaystyle \int _{0}^{1}f(\hat{z})\,\text{d}\hat{z}=1, & & \displaystyle\end{eqnarray}$$

to be consistent with the definition (2.8). Gray & Kokelaar (2010a , b ) used families of linear shear profiles given by

(2.12) $$\begin{eqnarray}\displaystyle f(\hat{z})=f_{L}(\hat{z})\equiv \unicode[STIX]{x1D6FC}+2(1-\unicode[STIX]{x1D6FC})\hat{z}, & & \displaystyle\end{eqnarray}$$

to derive their depth-averaged segregation equation, where the parameter $\unicode[STIX]{x1D6FC}\in [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,

(2.13) $$\begin{eqnarray}\displaystyle f(\hat{z})=f_{B}(\hat{z})\equiv {\textstyle \frac{5}{3}}(1-(1-\hat{z})^{3/2}), & & \displaystyle\end{eqnarray}$$

which can be derived as the steady uniform solution to the three-dimensional $\unicode[STIX]{x1D707}(I)$ -rheology for granular flows (e.g. GDR MiDi 2004; Gray & Edwards 2014). Substituting the inversely graded concentration (2.9) and velocity profile (2.10) into the flux integral in (2.6) gives

(2.14) $$\begin{eqnarray}\displaystyle \overline{\unicode[STIX]{x1D719}u}=\frac{1}{h}\int _{b}^{l}u\,\text{d}z=\bar{u}\int _{0}^{\bar{\unicode[STIX]{x1D719}}}f(\hat{z})\,\text{d}\hat{z}, & & \displaystyle\end{eqnarray}$$

which may then be inserted into the depth-integrated segregation equation (2.7) to give

(2.15) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h\bar{\unicode[STIX]{x1D719}})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(h\bar{\unicode[STIX]{x1D719}}\bar{u})-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(h\bar{u}G(\bar{\unicode[STIX]{x1D719}}))=0, & & \displaystyle\end{eqnarray}$$

where

(2.16) $$\begin{eqnarray}\displaystyle G(\bar{\unicode[STIX]{x1D719}})=\bar{\unicode[STIX]{x1D719}}-\int _{0}^{\bar{\unicode[STIX]{x1D719}}}f(\hat{z})\,\text{d}\hat{z}. & & \displaystyle\end{eqnarray}$$

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

(2.17) $$\begin{eqnarray}\displaystyle G(\bar{\unicode[STIX]{x1D719}})=G_{L}(\bar{\unicode[STIX]{x1D719}})\equiv (1-\unicode[STIX]{x1D6FC})\bar{\unicode[STIX]{x1D719}}(1-\bar{\unicode[STIX]{x1D719}}), & & \displaystyle\end{eqnarray}$$

as in Gray & Kokelaar (2010a ), and the Bagnold shear profile (2.13) giving

(2.18) $$\begin{eqnarray}\displaystyle G_{B}(\bar{\unicode[STIX]{x1D719}})\equiv {\textstyle \frac{2}{3}}(1-\bar{\unicode[STIX]{x1D719}})(1-(1-\bar{\unicode[STIX]{x1D719}})^{3/2}). & & \displaystyle\end{eqnarray}$$

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 $\unicode[STIX]{x1D6FC}=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-\bar{\unicode[STIX]{x1D719}})^{3/2}$ term in (2.18) results in complex values if round-off errors cause $\bar{\unicode[STIX]{x1D719}}$ to be slightly greater than unity. Though the linear profile (2.12) with $\unicode[STIX]{x1D6FC}=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.

Figure 7. (a) Plots of the linear (2.12) and Bagnold (2.13) shear profiles $f(\hat{z})$ . (b) The corresponding transport functions $G(\bar{\unicode[STIX]{x1D719}})$ given by (2.17) and (2.18) respectively. The value $\unicode[STIX]{x1D6FC}=1/7$ is chosen for the linear profiles so that the area under the curves in (b) is the same.

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, $\unicode[STIX]{x1D702}(x,t)=l-b$ , using the fact that $\unicode[STIX]{x1D702}=h\bar{\unicode[STIX]{x1D719}}$ , or the thickness of the large particle layer, $\unicode[STIX]{x1D705}(x,t)=h-\unicode[STIX]{x1D702}$ , as described in Gray & Kokelaar (2010a , b ). Here it shall be left in terms of the depth-averaged concentration of small particles $\bar{\unicode[STIX]{x1D719}}$ because this is more representative of what would actually be seen in overhead views of bidisperse experiments.

3 Segregation-mobility coupling

The large particle transport equation (2.15) may be solved for the depth-averaged concentration $\bar{\unicode[STIX]{x1D719}}$ for a prescribed flow thickness $h$ and bulk velocity $\bar{u}$ (e.g. Gray & Kokelaar 2010a , b ). In some cases $h$ and $\bar{u}$ can be inferred from experimental measurements (Johnson et al. 2012) 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 (Gray & Edwards 2014)

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(h\bar{u})=0, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h\bar{u})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\unicode[STIX]{x1D712}h\bar{u}^{2})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\frac{1}{2}gh^{2}\cos \unicode[STIX]{x1D701}\right)=ghS+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\unicode[STIX]{x1D708}h^{3/2}\frac{\unicode[STIX]{x2202}\bar{u}}{\unicode[STIX]{x2202}x}\right), & \displaystyle\end{eqnarray}$$

where $g$ is the constant of gravitational acceleration. The shape factor $\unicode[STIX]{x1D712}=\overline{u^{2}}/\bar{u}^{2}$ in (3.2) depends on the form of the velocity profile with depth. The Bagnold profile (2.13) gives a value $\unicode[STIX]{x1D712}=5/4$ but it shall be assumed that $\unicode[STIX]{x1D712}=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),

(3.3) $$\begin{eqnarray}\displaystyle S=\cos \unicode[STIX]{x1D701}\left(\tan \unicode[STIX]{x1D701}-\unicode[STIX]{x1D707}_{b}\text{sgn}(\bar{u})-\frac{\unicode[STIX]{x2202}b}{\unicode[STIX]{x2202}x}\right), & & \displaystyle\end{eqnarray}$$

where $\text{sgn}$ is the sign function and ensures friction always opposes the direction of motion. The effective basal friction coefficient $\unicode[STIX]{x1D707}_{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),

(3.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{b}(h,Fr,\bar{\unicode[STIX]{x1D719}})=\bar{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D707}_{b}^{S}(h,Fr)+(1-\bar{\unicode[STIX]{x1D719}})\unicode[STIX]{x1D707}_{b}^{L}(h,Fr), & & \displaystyle\end{eqnarray}$$

where

(3.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{b}^{S}(h,Fr)<\unicode[STIX]{x1D707}_{b}^{L}(h,Fr), & & \displaystyle\end{eqnarray}$$

are the basal friction coefficients for smooth small and frictional large particles, respectively, and are written as functions of thickness and Froude number,

(3.6) $$\begin{eqnarray}\displaystyle Fr=\frac{|\bar{u}|}{\sqrt{gh\cos \unicode[STIX]{x1D701}}}. & & \displaystyle\end{eqnarray}$$

It is assumed that the friction laws for the individual constituents are given by the dynamic friction law of Pouliquen & Forterre (2002),

(3.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{b}^{{\mathcal{N}}}(h,Fr)=\unicode[STIX]{x1D707}_{1}^{{\mathcal{N}}}+\frac{\unicode[STIX]{x1D707}_{2}^{{\mathcal{N}}}-\unicode[STIX]{x1D707}_{1}^{{\mathcal{N}}}}{(\unicode[STIX]{x1D6FD}^{{\mathcal{N}}}h)/({\mathcal{L}}^{{\mathcal{N}}}Fr)+1},\quad Fr>\unicode[STIX]{x1D6FD}^{{\mathcal{N}}}, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{N}}=S,L$ denotes small or large particles, respectively. The values $\unicode[STIX]{x1D707}_{1}^{{\mathcal{N}}}=\tan \unicode[STIX]{x1D701}_{1}^{{\mathcal{N}}}$ and $\unicode[STIX]{x1D707}_{2}^{{\mathcal{N}}}=\tan \unicode[STIX]{x1D701}_{2}^{{\mathcal{N}}}$ are constants, where angles $\unicode[STIX]{x1D701}_{1}^{{\mathcal{N}}}$ and $\unicode[STIX]{x1D701}_{2}^{{\mathcal{N}}}$ correspond to the minimum and maximum slope angles for which steady uniform flows are observed for a monodisperse material of constituent ${\mathcal{N}}$ . The length scales ${\mathcal{L}}^{{\mathcal{N}}}$ and dimensionless constants $\unicode[STIX]{x1D6FD}^{{\mathcal{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.

Table 1. Material parameters that will remain constant throughout this paper.

Strictly speaking, the individual basal friction laws (3.7) only hold providing $Fr>\unicode[STIX]{x1D6FD}^{{\mathcal{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.

The form of the final viscous term in the momentum equation (3.2) is motivated by the work done by Gray & Edwards (2014) for monodisperse flows, who used the $\unicode[STIX]{x1D707}(I)$ -rheology (GDR MiDi 2004; da Cruz et al. 2005; Jop et al. 2005, 2006) to incorporate more of the specific material properties into the depth-averaged governing equations. They showed that, to leading order, the $\unicode[STIX]{x1D707}(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 Gray & Edwards (2014) 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, $\unicode[STIX]{x1D708}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\unicode[STIX]{x2202}\bar{u}/\unicode[STIX]{x2202}x$ . Gray & Edwards (2014) were able to write the controlling coefficient $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}^{{\mathcal{N}}}$ explicitly in terms of the friction parameters of the monodisperse material as

(3.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D708}^{{\mathcal{N}}}=\frac{2{\mathcal{L}}^{{\mathcal{N}}}\sqrt{g}}{9\unicode[STIX]{x1D6FD}^{{\mathcal{N}}}}\frac{\sin \unicode[STIX]{x1D701}}{\sqrt{\cos \unicode[STIX]{x1D701}}}\left(\frac{\unicode[STIX]{x1D707}_{2}^{{\mathcal{N}}}-\tan \unicode[STIX]{x1D701}}{\tan \unicode[STIX]{x1D701}-\unicode[STIX]{x1D707}_{1}^{{\mathcal{N}}}}\right),\quad \unicode[STIX]{x1D701}_{1}^{{\mathcal{N}}}<\unicode[STIX]{x1D701}<\unicode[STIX]{x1D701}_{2}^{{\mathcal{N}}}. & & \displaystyle\end{eqnarray}$$

For the bidisperse flows being considered here it might be sensible to choose

(3.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}(\bar{\unicode[STIX]{x1D719}})=\bar{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D708}^{S}+(1-\bar{\unicode[STIX]{x1D719}})\unicode[STIX]{x1D708}^{L}, & & \displaystyle\end{eqnarray}$$

in an analogous manner to (3.4), where $\unicode[STIX]{x1D708}^{S}$ and $\unicode[STIX]{x1D708}^{L}$ are the coefficients for small and large particles and are given by (3.8). However, the coefficients $\unicode[STIX]{x1D708}^{S}$ and $\unicode[STIX]{x1D708}^{L}$ are only valid for slope angles $\unicode[STIX]{x1D701}_{1}^{{\mathcal{N}}}<\unicode[STIX]{x1D701}<\unicode[STIX]{x1D701}_{2}^{{\mathcal{N}}}$ , 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 $\unicode[STIX]{x1D707}(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) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D701}_{1}^{S}<\unicode[STIX]{x1D701}<\unicode[STIX]{x1D701}_{1}^{L}<\unicode[STIX]{x1D701}_{2}^{S}<\unicode[STIX]{x1D701}_{2}^{L}. & & \displaystyle\end{eqnarray}$$

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 $\unicode[STIX]{x1D708}>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.

4 Steady uniform flows

A simple solution to the system of equations (2.15), (3.1) and (3.2) is given by

(4.1a-c ) $$\begin{eqnarray}\displaystyle h=h_{0},\quad \bar{u}=\bar{u}_{0},\quad \bar{\unicode[STIX]{x1D719}}=\bar{\unicode[STIX]{x1D719}}_{0}, & & \displaystyle\end{eqnarray}$$

for constants $h_{0}>0$ , $\bar{u}_{0}>0$ , $\bar{\unicode[STIX]{x1D719}}_{0}\in [0,1]$ . 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,

(4.2) $$\begin{eqnarray}\displaystyle \tan \unicode[STIX]{x1D701}=\unicode[STIX]{x1D707}_{b}(h_{0},F,\bar{\unicode[STIX]{x1D719}}_{0}), & & \displaystyle\end{eqnarray}$$

where

(4.3) $$\begin{eqnarray}\displaystyle F=Fr_{0}=\frac{\bar{u}_{0}}{\sqrt{gh_{0}\cos \unicode[STIX]{x1D701}}}, & & \displaystyle\end{eqnarray}$$

is the steady uniform Froude number. Treating $h_{0}$ and $\bar{\unicode[STIX]{x1D719}}_{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

(4.4) $$\begin{eqnarray}\displaystyle AF^{2}+Bh_{0}F+Ch_{0}^{2}=0, & & \displaystyle\end{eqnarray}$$

where the coefficients are given by

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle A(\bar{\unicode[STIX]{x1D719}}_{0})=\bar{\unicode[STIX]{x1D719}}_{0}\unicode[STIX]{x1D707}_{2}^{S}+(1-\bar{\unicode[STIX]{x1D719}}_{0})\unicode[STIX]{x1D707}_{2}^{L}-\tan \unicode[STIX]{x1D701}, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle B(\bar{\unicode[STIX]{x1D719}}_{0})=\bar{\unicode[STIX]{x1D719}}_{0}(M^{S}\unicode[STIX]{x1D707}_{1}^{S}+M^{L}\unicode[STIX]{x1D707}_{2}^{S})+(1-\bar{\unicode[STIX]{x1D719}}_{0})(M^{S}\unicode[STIX]{x1D707}_{2}^{L}+M^{L}\unicode[STIX]{x1D707}_{1}^{L})-(M^{S}+M^{L})\tan \unicode[STIX]{x1D701}, & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle C(\bar{\unicode[STIX]{x1D719}}_{0})=(\bar{\unicode[STIX]{x1D719}}_{0}\unicode[STIX]{x1D707}_{1}^{S}+(1-\bar{\unicode[STIX]{x1D719}}_{0})\unicode[STIX]{x1D707}_{1}^{L}-\tan \unicode[STIX]{x1D701})M^{S}M^{L}, & \displaystyle\end{eqnarray}$$

with $M^{{\mathcal{N}}}=\unicode[STIX]{x1D6FD}^{{\mathcal{N}}}/{\mathcal{L}}^{{\mathcal{N}}}$ . For a slope angle in the range given by (3.10), it can be seen that $A(\bar{\unicode[STIX]{x1D719}}_{0})>0$ for all $\bar{\unicode[STIX]{x1D719}}_{0}\in [0,1]$ , whereas $C(\bar{\unicode[STIX]{x1D719}}_{0})>0$ for $\bar{\unicode[STIX]{x1D719}}_{0}<\bar{\unicode[STIX]{x1D719}}_{0}^{\ast }$ and $C(\bar{\unicode[STIX]{x1D719}}_{0})<0$ for $\bar{\unicode[STIX]{x1D719}}_{0}>\bar{\unicode[STIX]{x1D719}}_{0}^{\ast }$ , where

(4.8) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D719}}_{0}^{\ast }=\frac{\unicode[STIX]{x1D707}_{1}^{L}-\tan \unicode[STIX]{x1D701}}{\unicode[STIX]{x1D707}_{1}^{L}-\unicode[STIX]{x1D707}_{1}^{S}}. & & \displaystyle\end{eqnarray}$$

Consequently, the steady-state Froude number, found by taking the positive root of (4.4),

(4.9) $$\begin{eqnarray}\displaystyle F=h_{0}\left(\frac{-B+\sqrt{B^{2}-4AC}}{2A}\right), & & \displaystyle\end{eqnarray}$$

is only positive providing that $\bar{\unicode[STIX]{x1D719}}_{0}>\bar{\unicode[STIX]{x1D719}}_{0}^{\ast }$ , 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},\bar{\unicode[STIX]{x1D719}}_{0})$ , along with the regions where $\bar{\unicode[STIX]{x1D719}}_{0}<\bar{\unicode[STIX]{x1D719}}_{0}^{\ast }$ . In the pure small limit ( $\bar{\unicode[STIX]{x1D719}}_{0}=1$ ) the expression (4.9) reduces to that given in Gray & Edwards (2014),

(4.10) $$\begin{eqnarray}\displaystyle F=F(h_{0})=\frac{M^{S}h_{0}(\tan \unicode[STIX]{x1D701}-\unicode[STIX]{x1D707}_{1}^{S})}{\unicode[STIX]{x1D707}_{2}^{S}-\tan \unicode[STIX]{x1D701}}, & & \displaystyle\end{eqnarray}$$

which can also be derived from the more straightforward force balance, $\tan \unicode[STIX]{x1D701}=\unicode[STIX]{x1D707}_{b}^{S}(h_{0},F)$ . The corresponding steady uniform velocities $\bar{u}_{0}(h_{0},\bar{\unicode[STIX]{x1D719}}_{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.

Figure 8. Contour plots of the steady uniform Froude number $F(h_{0},\bar{\unicode[STIX]{x1D719}}_{0})$ , given by (4.9). The shaded regions represent where $\bar{\unicode[STIX]{x1D719}}_{0}<\bar{\unicode[STIX]{x1D719}}_{0}^{\ast }$ (given by (4.8)), meaning there are too many frictional large particles for steady uniform flow.

5 Linear stability analysis

5.1 Non-dimensionalisation

Assume the values $h_{0}$ and $\bar{\unicode[STIX]{x1D719}}_{0}$ are chosen such that a steady state $h=h_{0}$ , $\bar{\unicode[STIX]{x1D719}}=\bar{\unicode[STIX]{x1D719}}_{0}$ , $\bar{u}=\bar{u}_{0}(h_{0},\bar{\unicode[STIX]{x1D719}}_{0})>0$ exists with corresponding Froude number $F>0$ , as described in the previous section. It is then convenient to introduce the scalings

(5.1a-d ) $$\begin{eqnarray}\displaystyle h=h_{0}{\hat{h}},\quad \bar{u}=\bar{u}_{0}\hat{\bar{u}},\quad x=h_{0}\hat{x},\quad t=\frac{h_{0}}{\bar{u}_{0}}\hat{t}, & & \displaystyle\end{eqnarray}$$

where the hats denote dimensionless quantities. Note that the depth-averaged concentration $\bar{\unicode[STIX]{x1D719}}$ is already non-dimensional. Using these scalings, the governing equations may be written as

(5.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(h\bar{u})=0, & \displaystyle\end{eqnarray}$$

(5.3) $$\begin{eqnarray}\displaystyle & \displaystyle F^{2}h\left(\frac{\unicode[STIX]{x2202}\bar{u}}{\unicode[STIX]{x2202}t}+\bar{u}\frac{\unicode[STIX]{x2202}\bar{u}}{\unicode[STIX]{x2202}x}\right)+h\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}=h\left(\tan \unicode[STIX]{x1D701}-\unicode[STIX]{x1D707}_{b}(h,\bar{u},\bar{\unicode[STIX]{x1D719}})\right)+\frac{F^{2}}{R}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(h^{3/2}\frac{\unicode[STIX]{x2202}\bar{u}}{\unicode[STIX]{x2202}x}\right), & \displaystyle\end{eqnarray}$$

(5.4) $$\begin{eqnarray}\displaystyle & \displaystyle h\left(\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x2202}t}+\bar{u}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x2202}x}\right)-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(h\bar{u}G(\bar{\unicode[STIX]{x1D719}}))=0, & \displaystyle\end{eqnarray}$$

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 $\bar{u}>0$ . The basal friction coefficient is now written in terms of the non-dimensional variables as

(5.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{b}(h,\bar{u},\bar{\unicode[STIX]{x1D719}})=\bar{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D707}_{b}^{S}(h,\bar{u})+(1-\bar{\unicode[STIX]{x1D719}})\unicode[STIX]{x1D707}_{b}^{L}(h,\bar{u}), & & \displaystyle\end{eqnarray}$$

where the individual friction coefficients for each type of particle are given by

(5.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{b}^{{\mathcal{N}}}(h,\bar{u})=\unicode[STIX]{x1D707}_{1}^{{\mathcal{N}}}+\frac{\unicode[STIX]{x1D707}_{2}^{{\mathcal{N}}}-\unicode[STIX]{x1D707}_{1}^{{\mathcal{N}}}}{(\unicode[STIX]{x1D6FE}^{{\mathcal{N}}}h^{3/2})/\bar{u}+1}, & & \displaystyle\end{eqnarray}$$

for constants $\unicode[STIX]{x1D6FE}^{{\mathcal{N}}}=(\unicode[STIX]{x1D6FD}^{{\mathcal{N}}}h_{0})/({\mathcal{L}}^{{\mathcal{N}}}F)$ . The granular Reynolds number in (5.3) is defined as

(5.7) $$\begin{eqnarray}\displaystyle R=\frac{\bar{u}_{0}\sqrt{h_{0}}}{\unicode[STIX]{x1D708}}, & & \displaystyle\end{eqnarray}$$

and will typically take large values since the viscous terms are small in magnitude compared to the standard shallow-water contributions, as shown by Gray & Edwards (2014). The Reynolds number is a function of both steady uniform flow thickness $h_{0}$ and concentration $\bar{\unicode[STIX]{x1D719}}_{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$ , $\bar{u}=1$ , $\bar{\unicode[STIX]{x1D719}}=\bar{\unicode[STIX]{x1D719}}_{0}$ , and so the variables are perturbed about this base state,

(5.8) $$\begin{eqnarray}\displaystyle & \displaystyle h(x,t)=1+h_{1}(x,t),\quad |h_{1}|\ll 1, & \displaystyle\end{eqnarray}$$

(5.9) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{u}(x,t)=1+\bar{u}_{1}(x,t),\quad |\bar{u}_{1}|\ll 1, & \displaystyle\end{eqnarray}$$

(5.10) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D719}}(x,t)=\bar{\unicode[STIX]{x1D719}}_{0}+\bar{\unicode[STIX]{x1D719}}_{1}(x,t),\quad |\bar{\unicode[STIX]{x1D719}}_{1}|\ll 1. & \displaystyle\end{eqnarray}$$

The linearised governing equations then become

(5.11) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}h_{1}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}h_{1}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x}=0, & \displaystyle\end{eqnarray}$$

(5.12) $$\begin{eqnarray}\displaystyle & \displaystyle F^{2}\left(\frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x}\right)+\frac{\unicode[STIX]{x2202}h_{1}}{\unicode[STIX]{x2202}x}=-(h_{1}\unicode[STIX]{x1D707}_{h}+\bar{u}_{1}\unicode[STIX]{x1D707}_{\bar{u}}+\bar{\unicode[STIX]{x1D719}}_{1}\unicode[STIX]{x1D707}_{\bar{\unicode[STIX]{x1D719}}})+\frac{F^{2}}{R}\frac{\unicode[STIX]{x2202}^{2}\bar{u}_{1}}{\unicode[STIX]{x2202}x^{2}}, & \displaystyle\end{eqnarray}$$

(5.13) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D719}}_{1}}{\unicode[STIX]{x2202}t}-G_{0}\left(\frac{\unicode[STIX]{x2202}h_{1}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}\bar{u}_{1}}{\unicode[STIX]{x2202}x}\right)+(1-G_{0}^{\prime })\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D719}}_{1}}{\unicode[STIX]{x2202}x}=0, & \displaystyle\end{eqnarray}$$

where the simplified notation

(5.14a,b ) $$\begin{eqnarray}\displaystyle G_{0}\equiv G(\bar{\unicode[STIX]{x1D719}}_{0})={\textstyle \frac{6}{7}}\bar{\unicode[STIX]{x1D719}}_{0}(1-\bar{\unicode[STIX]{x1D719}}_{0}),\quad G_{0}^{\prime }\equiv G^{\prime }(\bar{\unicode[STIX]{x1D719}}_{0})={\textstyle \frac{6}{7}}(1-2\bar{\unicode[STIX]{x1D719}}_{0}), & & \displaystyle\end{eqnarray}$$

has been introduced to denote the steady-state values of the transport function and its derivative. The basal friction law contributes via terms

(5.15) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}_{h}=\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{b}}{\unicode[STIX]{x2202}h}\right|_{(1,1,\bar{\unicode[STIX]{x1D719}}_{0})}=-\frac{3}{2}\left(\bar{\unicode[STIX]{x1D719}}_{0}\unicode[STIX]{x1D6E4}^{S}+(1-\bar{\unicode[STIX]{x1D719}}_{0})\unicode[STIX]{x1D6E4}^{L}\right)<0, & \displaystyle\end{eqnarray}$$

(5.16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}_{\bar{u}}=\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{b}}{\unicode[STIX]{x2202}\bar{u}}\right|_{(1,1,\bar{\unicode[STIX]{x1D719}}_{0})}=\bar{\unicode[STIX]{x1D719}}_{0}\unicode[STIX]{x1D6E4}^{S}+(1-\bar{\unicode[STIX]{x1D719}}_{0})\unicode[STIX]{x1D6E4}^{L}>0, & \displaystyle\end{eqnarray}$$

(5.17) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}_{\bar{\unicode[STIX]{x1D719}}}=\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{b}}{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D719}}}\right|_{(1,1,\bar{\unicode[STIX]{x1D719}}_{0})}=\unicode[STIX]{x1D707}_{b}^{S}(1,1)-\unicode[STIX]{x1D707}_{b}^{L}(1,1)<0, & \displaystyle\end{eqnarray}$$

where the positive constants $\unicode[STIX]{x1D6E4}^{S}$ , $\unicode[STIX]{x1D6E4}^{L}$ , are given by

(5.18) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}^{{\mathcal{N}}}=\frac{\unicode[STIX]{x1D6FE}^{{\mathcal{N}}}(\unicode[STIX]{x1D707}_{2}^{{\mathcal{N}}}-\unicode[STIX]{x1D707}_{1}^{{\mathcal{N}}})}{(1+\unicode[STIX]{x1D6FE}^{{\mathcal{N}}})^{2}}, & & \displaystyle\end{eqnarray}$$

for ${\mathcal{N}}=L,S$ . Now seek normal mode solutions of the form

(5.19) $$\begin{eqnarray}\displaystyle (h_{1},\bar{u}_{1},\bar{\unicode[STIX]{x1D719}}_{1})=(H,U,\unicode[STIX]{x1D6F7})\text{e}^{\unicode[STIX]{x1D70E}t}\text{e}^{\text{i}kx}, & & \displaystyle\end{eqnarray}$$

Figure 9. Plots of the growth rates $\unicode[STIX]{x1D70E}_{M}(k)$ for both the inviscid and viscous equations. (a) When $F<F_{c}$ both growth rates are positive for all wavenumbers $k$ , meaning that perturbations grow in time and the base state is unstable. In the large wavenumber limits the inviscid curves tend to a positive constant (5.30) (dash-dot lines), whereas the viscous values decay to zero according to the asymptotics (5.34), giving an internal maximum $\unicode[STIX]{x1D70E}_{Max}$ at finite wavenumber $k_{M}$ . (b) For $F>F_{c}$ the inviscid values tend to the constant (5.29) as $k\longrightarrow \infty$ . In the viscous case there is an internal maximum as well as a cutoff wavenumber $k_{c}$ , above which all perturbations are stable. The different cases are established by fixing the depth-averaged concentration $\bar{\unicode[STIX]{x1D719}}_{0}=0.8$ , giving a critical Froude number $F_{c}=1.94$ , and varying the flow thickness $h_{0}$ . The values $h_{0}=2~\text{mm}$ and $h_{0}=3~\text{mm}$ give corresponding Froude numbers $F=1.57$ and $F=2.36$ for (a) and (b) respectively. The coefficient in the effective viscosity is set to be $\unicode[STIX]{x1D708}=0.001~\text{m}^{3/2}~\text{s}^{-1}$ , but qualitatively similar behaviour is found for all positive values.

Figure 10. Plots of the growth rates $\unicode[STIX]{x1D70E}_{M}(k)$ at the critical Froude number $F=F_{c}$ . The inviscid growth rates grow unboundedly $\unicode[STIX]{x1D70E}_{M}\propto k^{1/2}$ as $k\longrightarrow \infty$ according to the scaling (5.31). The viscous values are also unstable for all wavenumbers but $\unicode[STIX]{x1D70E}_{M}$ remains bounded and decays to zero like $1/k^{4}$ as $k\longrightarrow \infty$ (5.35), giving an internal maximum $\unicode[STIX]{x1D70E}_{Max}=\max (\unicode[STIX]{x1D70E}_{M})$ at $k=k_{M}$ . The parameters used are $\unicode[STIX]{x1D708}=0.001~\text{m}^{3/2}~\text{s}^{-1}$ and $\bar{\unicode[STIX]{x1D719}}_{0}=0.8$ , giving a critical Froude number $F_{c}=1.94$ at flow thickness $h_{0}=2.47~\text{mm}$ .

where, for temporal stability analysis, the wavenumber $k$ is real and $\unicode[STIX]{x1D70E}(k)=\unicode[STIX]{x1D70E}_{R}(k)+\text{i}\unicode[STIX]{x1D70E}_{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

(5.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63C}\boldsymbol{W}=\unicode[STIX]{x1D70E}\boldsymbol{W}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{W}=(H,U,\unicode[STIX]{x1D6F7})^{\text{T}}$ and the matrix $\unicode[STIX]{x1D63C}$ is given by

(5.21) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63C}=\left(\begin{array}{@{}ccc@{}}-\text{i}k & -\text{i}k & 0\\ \displaystyle -\frac{\unicode[STIX]{x1D707}_{h}}{F^{2}}-\frac{\text{i}k}{F^{2}} & \displaystyle -\frac{\unicode[STIX]{x1D707}_{\bar{u}}}{F^{2}}-\text{i}k-\frac{k^{2}}{R} & \displaystyle -\frac{\unicode[STIX]{x1D707}_{\bar{\unicode[STIX]{x1D719}}}}{F^{2}}\\ \text{i}kG_{0} & \text{i}kG_{0} & \text{i}k(G_{0}^{\prime }-1)\end{array}\right). & & \displaystyle\end{eqnarray}$$

Equation (5.20) has non-zero solutions for $\boldsymbol{W}$ if and only if $|\unicode[STIX]{x1D63C}-\unicode[STIX]{x1D70E}\unicode[STIX]{x1D644}|=0$ , (where $\unicode[STIX]{x1D644}$ is the $3\times 3$ identity matrix) leading to a cubic characteristic polynomial

(5.22) $$\begin{eqnarray}\displaystyle f(\unicode[STIX]{x1D70E})\equiv f_{0}+f_{1}\unicode[STIX]{x1D70E}+f_{2}\unicode[STIX]{x1D70E}^{2}+\unicode[STIX]{x1D70E}^{3}=0, & & \displaystyle\end{eqnarray}$$

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 $\unicode[STIX]{x1D70E}^{(1)}$ , $\unicode[STIX]{x1D70E}^{(2)}$ , $\unicode[STIX]{x1D70E}^{(3)}$ with corresponding real parts $\unicode[STIX]{x1D70E}_{R}^{(1)}$ , $\unicode[STIX]{x1D70E}_{R}^{(2)}$ , $\unicode[STIX]{x1D70E}_{R}^{(3)}$ . The growth rate $\unicode[STIX]{x1D70E}_{M}$ is then given by the maximum of these values,

(5.23) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}_{M}=\max (\unicode[STIX]{x1D70E}_{R}^{(1)},\unicode[STIX]{x1D70E}_{R}^{(2)},\unicode[STIX]{x1D70E}_{R}^{(3)}). & & \displaystyle\end{eqnarray}$$

If $\unicode[STIX]{x1D70E}_{M}(k)<0$ for all values of $k$ then all perturbations decay exponentially in time and the base state is linearly stable. On the other hand, if there exists a wavenumber such that $\unicode[STIX]{x1D70E}_{M}(k)>0$ then this perturbation grows in time and the steady flow is linearly unstable. Example plots of $\unicode[STIX]{x1D70E}_{M}(k)$ are shown in figures 9 and 10 for both the new viscous equations and the inviscid equivalent ( $\unicode[STIX]{x1D708}=0$ ). This inviscid case corresponds to a one-dimensional version of Woodhouse et al.’s (2012) model, although they used the original exponential form of the basal friction law (Pouliquen 1999a ) as opposed to (3.7). In the viscous regime, the coefficient $\unicode[STIX]{x1D708}$ is set to be $0.001~\text{m}^{3/2}~\text{s}^{-1}$ for all stability results shown on figures 9 and 10, but qualitatively similar behaviour is found for different values of $\unicode[STIX]{x1D708}>0$ .

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 9 b).

An important result in Woodhouse et al. (2012) was the discovery of unbounded growth rates, $\unicode[STIX]{x1D70E}_{M}\longrightarrow \infty$ as $k\longrightarrow \infty$ (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

(5.24) $$\begin{eqnarray}\displaystyle \frac{\text{d}x}{\text{d}t}=\unicode[STIX]{x1D706}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D706}$ is the characteristic wavespeed. For the non-dimensional mass and momentum balance equations (5.2), (5.3) there are two wavespeeds given by

(5.25a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}^{(1)}=\bar{u}+\frac{1}{F}\sqrt{h},\quad \unicode[STIX]{x1D706}^{(2)}=\bar{u}-\frac{1}{F}\sqrt{h}, & & \displaystyle\end{eqnarray}$$

(see, for example, Courant & Hilbert 1962) whereas for the large particle transport equation (5.4) the wavespeed is

(5.26) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}^{(3)}=\bar{u}(1-G^{\prime }(\bar{\unicode[STIX]{x1D719}})). & & \displaystyle\end{eqnarray}$$

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,\bar{u},\bar{\unicode[STIX]{x1D719}})=(1,1,\bar{\unicode[STIX]{x1D719}}_{0})$ , and noting that $G_{0}^{\prime }<1$ , it can be seen that $\unicode[STIX]{x1D706}^{(1)}>0$ and $\unicode[STIX]{x1D706}^{(3)}>0$ . The other wavespeed $\unicode[STIX]{x1D706}^{(2)}$ is positive for $F>1$ and negative for $F<1$ . Most importantly, at the critical Froude number

(5.27) $$\begin{eqnarray}\displaystyle F=F_{c}\equiv \frac{1}{|G_{0}^{\prime }|}, & & \displaystyle\end{eqnarray}$$

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}^{\prime }$ , 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\gg 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.

5.3 Inviscid high wavenumber asymptotics

For completeness, we begin by considering the inviscid case, $\unicode[STIX]{x1D708}=0$ . The real part of the three roots is found to behave like

(5.28) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{R}^{(1)}\sim \frac{(FG_{0}^{\prime }-1)(F\unicode[STIX]{x1D707}_{h}-\unicode[STIX]{x1D707}_{\bar{u}})+FG_{0}(1-F)\unicode[STIX]{x1D707}_{\bar{\unicode[STIX]{x1D719}}}}{2F^{2}(FG_{0}^{\prime }-1)}, & \displaystyle\end{eqnarray}$$

(5.29) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{R}^{(2)}\sim \frac{-(FG_{0}^{\prime }+1)(F\unicode[STIX]{x1D707}_{h}+\unicode[STIX]{x1D707}_{\bar{u}})+FG_{0}(1+F)\unicode[STIX]{x1D707}_{\bar{\unicode[STIX]{x1D719}}}}{2F^{2}(FG_{0}^{\prime }+1)}, & \displaystyle\end{eqnarray}$$

(5.30) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{R}^{(3)}\sim \frac{G_{0}(1-G_{0}^{\prime })\unicode[STIX]{x1D707}_{\bar{\unicode[STIX]{x1D719}}}}{F^{2}(G_{0}^{\prime })^{2}-1}, & \displaystyle\end{eqnarray}$$

in the high wavenumber limit (see § A.1), meaning that the growth rate $\unicode[STIX]{x1D70E}_{M}$ tends to a constant value as $k\longrightarrow \infty$ . From the definition (5.17) and the assumption (3.5) that large particles experience greater frictional forces, it follows that $\unicode[STIX]{x1D707}_{\bar{\unicode[STIX]{x1D719}}}<0$ . Using this, and the fact that $G_{0}^{\prime }<1$ by (5.14), it can be seen that that $\unicode[STIX]{x1D70E}_{R}^{(3)}>0$ for $F<F_{c}$ , where the critical Froude number is defined by (5.27). The growth rate $\unicode[STIX]{x1D70E}_{M}$ is therefore always positive for $F<F_{c}$ , meaning perturbations grow in time and the base state is unstable as $k\longrightarrow \infty$ (figure 9 a). For $F>F_{c}$ the third root $\unicode[STIX]{x1D70E}_{R}^{(3)}$ is now stable in the asymptotic limit and one must investigate the sign of the other roots $\unicode[STIX]{x1D70E}_{R}^{(1)}$ and $\unicode[STIX]{x1D70E}_{R}^{(2)}$ . For the parameters given in table 1 at least one of (5.28) or (5.29) is positive for all regions of $(h_{0},\bar{\unicode[STIX]{x1D719}}_{0})$ space, meaning that $\unicode[STIX]{x1D70E}_{M}>0$ and the base state is again unstable (figure 9 b).

At the critical Froude number $F=F_{c}$ the third root $\unicode[STIX]{x1D70E}_{R}^{(3)}$ and either $\unicode[STIX]{x1D70E}_{R}^{(1)}$ or $\unicode[STIX]{x1D70E}_{R}^{(2)}$ become infinite, depending on the sign of $G_{0}^{\prime }$ . In these singular cases, the leading-order growth rate is instead given by the distinguished limit

(5.31) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}_{R}^{(\pm )}\sim \pm {\textstyle \frac{1}{2}}|G_{0}G_{0}^{\prime }(1-G_{0}^{\prime })\unicode[STIX]{x1D707}_{\bar{\unicode[STIX]{x1D719}}}|^{1/2}k^{1/2}, & & \displaystyle\end{eqnarray}$$

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

(5.32) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{R}^{(1)}\sim -\frac{1}{R}k^{2}, & \displaystyle\end{eqnarray}$$

(5.33) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{R}^{(2)}\sim -\frac{R}{F^{2}}, & \displaystyle\end{eqnarray}$$

(5.34) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{R}^{(3)}\sim \frac{R^{2}G_{0}(1-G_{0}^{\prime })(F^{2}(G_{0}^{\prime })^{2}-1)\unicode[STIX]{x1D707}_{\bar{\unicode[STIX]{x1D719}}}}{F^{4}(G_{0}^{\prime })^{2}}\frac{1}{k^{2}}, & \displaystyle\end{eqnarray}$$

for $k\gg 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}^{\prime }<1$ , $\unicode[STIX]{x1D707}_{\bar{\unicode[STIX]{x1D719}}}<0$ , and $F^{2}(G_{0}^{\prime })^{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\longrightarrow \infty$ , and instead decays to zero like $1/k^{2}$ (figure 9 a). This is significant because it means that the growth rate $\unicode[STIX]{x1D70E}_{M}(k)$ will take its maximum value $\unicode[STIX]{x1D70E}_{Max}=\max (\unicode[STIX]{x1D70E}_{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 9 b and 11 b). 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.

Figure 11. (a) Plots of the maximum growth rate $\unicode[STIX]{x1D70E}_{Max}$ against steady-state Froude number $F$ . The inviscid values are infinite at $F=F_{c}$ , whereas $\unicode[STIX]{x1D70E}_{Max}$ remains bounded for all Froude numbers when viscosity is included. (b) 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 $\unicode[STIX]{x1D708}$ lead to smaller $\unicode[STIX]{x1D70E}_{Max}$ and $k_{c}$ . In both plots the depth-averaged concentration is fixed at $\bar{\unicode[STIX]{x1D719}}_{0}=0.8$ and $F$ is varied by changing the steady-state thickness $h_{0}$ .

At the critical Froude number the expression (5.34) reduces to zero (since $F^{2}(G_{0}^{\prime })^{2}=1$ ) and a new asymptotic analysis leads to

(5.35) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}_{R}^{(3)}\sim 2R^{3}G_{0}^{2}(G_{0}^{\prime })^{2}(1-G_{0}^{\prime })^{2}\unicode[STIX]{x1D707}_{\bar{\unicode[STIX]{x1D719}}}^{2}\frac{1}{k^{4}}, & & \displaystyle\end{eqnarray}$$

in the high wavenumber limit. The expression (5.35) is strictly positive, meaning the root is unstable for $k\gg 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\longrightarrow \infty$ (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 $\unicode[STIX]{x1D70E}_{Max}$ as a function of steady uniform Froude number $F$ . For the inviscid equations $\unicode[STIX]{x1D70E}_{Max}\longrightarrow \infty$ as $F\longrightarrow 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 $\unicode[STIX]{x1D708}=0.001~\text{m}^{3/2}~\text{s}^{-1}$ , and qualitatively similar behaviour is found for all positive choices of $\unicode[STIX]{x1D708}$ . The specific value does not affect whether the equations are well posed, but it does have an influence on the stability of steady uniform flows. Increasing $\unicode[STIX]{x1D708}$ , and consequently reducing the Reynolds number $R$ , leads to lower maximum growth rates $\unicode[STIX]{x1D70E}_{Max}$ (figure 11 a) 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 $\unicode[STIX]{x1D708}$ (figure 11 b).

6 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). Woodhouse et al. (2012) were 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.

6.1 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. Gray & Edwards 2014; 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 $\unicode[STIX]{x1D707}(I)$ -rheology of Gray & Edwards (2014) to two dimensions. Using their work as a base and introducing a cross-slope coordinate $y$ and depth-averaged velocity $\bar{v}$ , the two-dimensional segregation-mobility feedback equations (reverting to dimensional variables) are

(6.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(h\bar{\boldsymbol{u}})=0, & \displaystyle\end{eqnarray}$$

(6.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h\bar{\boldsymbol{u}})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(h\bar{\boldsymbol{u}}\otimes \bar{\boldsymbol{u}})+\unicode[STIX]{x1D735}\left(\frac{1}{2}gh^{2}\cos \unicode[STIX]{x1D701}\right)=gh\boldsymbol{S}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D708}h^{3/2}\bar{\unicode[STIX]{x1D63F}}), & \displaystyle\end{eqnarray}$$

(6.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h\bar{\unicode[STIX]{x1D719}})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(h\bar{\boldsymbol{u}}\bar{\unicode[STIX]{x1D719}})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(h\bar{\boldsymbol{u}}G(\bar{\unicode[STIX]{x1D719}}))=0, & \displaystyle\end{eqnarray}$$

where $\bar{\boldsymbol{u}}=(\bar{u},\bar{v})$ is the depth-averaged velocity vector, $\unicode[STIX]{x1D735}=(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x,\unicode[STIX]{x2202}/\unicode[STIX]{x2202}y)$ denotes the gradient operator in $(x,y)$ space, ‘ $\boldsymbol{\cdot }$ ’ is the dot product and $\otimes$ the dyadic product. Assuming no basal topography, the source terms $\boldsymbol{S}=(S_{x},S_{y})$ are given by

(6.4a,b ) $$\begin{eqnarray}\displaystyle S_{x}=\cos \unicode[STIX]{x1D701}\left(\tan \unicode[STIX]{x1D701}-\unicode[STIX]{x1D707}_{b}\frac{\bar{u}}{|\bar{\boldsymbol{u}}|}\right),\quad S_{y}=-\cos \unicode[STIX]{x1D701}\left(\unicode[STIX]{x1D707}_{b}\frac{\bar{v}}{|\bar{\boldsymbol{u}}|}\right), & & \displaystyle\end{eqnarray}$$

where $|\bar{\boldsymbol{u}}|=(\bar{u}^{2}+\bar{v}^{2})^{1/2}$ is the magnitude of the velocity. The basal friction law $\unicode[STIX]{x1D707}_{b}(h,Fr,\bar{\unicode[STIX]{x1D719}})$ is defined by (3.4) and (3.7) using the more general Froude number definition

(6.5) $$\begin{eqnarray}\displaystyle Fr=\frac{|\bar{\boldsymbol{u}}|}{\sqrt{gh\cos \unicode[STIX]{x1D701}}}. & & \displaystyle\end{eqnarray}$$

The coefficient $\unicode[STIX]{x1D708}>0$ remains a free parameter for slope angles in the range (3.10), and the depth-integrated strain-rate tensor is

(6.6) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D63F}}={\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}\bar{\boldsymbol{u}}+(\unicode[STIX]{x1D735}\bar{\boldsymbol{u}})^{\text{T}}). & & \displaystyle\end{eqnarray}$$

The generalised large particle transport equation (6.3) can be derived in a similar manner to the one-dimensional case (2.15), providing assumptions are made about the vertical variation of the transverse velocity. Whilst the downslope velocity magnitudes are typically much larger than their lateral counterparts, the thinness of the flow means that material responds to shear comparably in both directions, and thus the shear 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 $\unicode[STIX]{x1D6FC}=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.

6.2 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 (Johnson et al. 2012), granular roll waves (Gray & Edwards 2014; 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\longrightarrow 0$ . To get around this potential problem a minimum flow thickness smaller than one particle diameter, $h_{min}=10^{-4}~\text{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

(6.7a-d ) $$\begin{eqnarray}\displaystyle h(x,y,0)=0,\quad \bar{u}(x,y,0)=0,\quad \bar{v}(x,y,0)=0,\quad \bar{\unicode[STIX]{x1D719}}(x,y,0)=0, & & \displaystyle\end{eqnarray}$$

are prescribed. Simulations are carried out on a computational domain $L_{x}\times L_{y}$ discretized over $N_{x}\times N_{y}$ grid cells, giving a spatial resolution of $N_{x}/L_{x}\times 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

(6.8a-d ) $$\begin{eqnarray}\displaystyle h(0,y,t)=(1-\text{e}^{-10t})h_{0},\quad \bar{u}(0,y,t)=\bar{u}_{0},\quad \bar{v}(0,y,t)=0,\quad \bar{\unicode[STIX]{x1D719}}(0,y,t)=\bar{\unicode[STIX]{x1D719}}_{0}+\unicode[STIX]{x1D6FF}\bar{\unicode[STIX]{x1D719}}_{p}(y), & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

are enforced. These correspond to the steady uniform flow of § 4, with the thickness prefactor $(1-\text{e}^{-10t})$ used to ensure a smooth transition from an empty slope. Note that the simulations represent a continuous inflow (akin to figures 3 a, 4 a) as opposed to a finite release of material (figures 3 b, 4 b). To introduce transverse variation into the system the inflow concentration is perturbed with a periodic function

(6.9) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D719}}_{p}(y)=\sin (2\unicode[STIX]{x03C0}y/L_{y}), & & \displaystyle\end{eqnarray}$$

of magnitude $\unicode[STIX]{x1D6FF}=8\times 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.

Figure 12. Numerical solutions of the system of PDEs (6.1)–(6.3) showing the depth averaged concentration $\bar{\unicode[STIX]{x1D719}}$ at times $t=1.9~\text{s}$ (non-dimensional $\hat{t}=198$ ),  $t=3.3~\text{s}$ ( $\hat{t}=343$ ),  $t=5.1~\text{s}$ ( $\hat{t}=530$ ),  $t=6.2~\text{s}$ ( $\hat{t}=645$ ) and $t=7.5~\text{s}$ ( $\hat{t}=780$ ). A large-rich region quickly develops at the flow front $(t=1.9~\text{s})$ , starts to become unstable $(t=5.1~\text{s})$ and develops into fingers $(t=6.2~\text{s})$ which elongate and coarsen over time. It is clear that lateral levees bounding the fingers consist predominantly of large particles. The colour scheme has been chosen to mimic experiments, with turquoise representing the region where the chute is empty, $h<h_{min}$ . Parameters are $h_{0}=2~\text{mm}$ , $\bar{\unicode[STIX]{x1D719}}_{0}=0.8$ , $\bar{u}_{0}=0.208~\text{ms}^{-1}$ , $F=1.57$ , $\unicode[STIX]{x1D708}=0.01~\text{m}^{3/2}~\text{s}^{-1}$ , $N_{x}/L_{x}=N_{y}/L_{y}=750~\text{m}^{-1}$ . Note that the axes limits $L_{x}=1.5~\text{m}$ , $L_{y}=0.5~\text{m}$ correspond to non-dimensional values of $\hat{L}_{x}=750$ and $\hat{L}_{y}=250$ respectively. Supplementary movie 3 available online.

Figure 13. Numerical solutions of the system of PDEs (6.1)–(6.3) showing the speed $|\bar{\boldsymbol{u}}|=(\bar{u}^{2}+\bar{v}^{2})^{1/2}$ at times $t=1.9~\text{s}$ , $t=3.3~\text{s}$ , $t=5.1~\text{s}$ , $t=6.2~\text{s}$ and $t=7.5~\text{s}$ . Once the flow breaks up into fingers the lateral levees are close to stationary, with areas of much faster flow down the central, channelised regions representing the more mobile interior. Parameters are the same as in figure 12, meaning that a dimensional velocity of $|\bar{\boldsymbol{u}}|=0.208~\text{m}~\text{s}^{-1}$ would correspond to a non-dimensional value of $|\hat{\bar{\boldsymbol{u}}}|=1$ . Supplementary movie 3 available online.

Figure 14. Numerical solutions of the system of PDEs (6.1)–(6.3) showing the flow thickness $h$ . Simulations are shown at times $t=1.9$ s, $t=3.3$ s, $t=5.1~\text{s}$ , $t=6.2~\text{s}$ and $t=7.5~\text{s}$ . The propagating front breaks up into a series of fingers going to zero thickness at the boundaries. A region of thicker flow follows behind the main front. Parameters used are the same as figures 12 and 13, meaning a dimensional thickness of $h=2~\text{mm}$ represents ${\hat{h}}=1$ in non-dimensional terms. Supplementary movie 3 available online.

Figures 12–14 and supplementary movie 3 show the resulting contours of the flow height $h$ , depth-averaged concentration $\bar{\unicode[STIX]{x1D719}}$ and speed $|\bar{\boldsymbol{u}}|$ . 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 ( $\bar{\unicode[STIX]{x1D719}}=1$ ). Assuming no transverse variation, the monodisperse equations admit a steady travelling front solution (Gray & Edwards 2014), 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 point in the simulations. Whilst this is not a major problem for the fingers that form from a continuous source (figures 3 a, 4 a), as the inflow is cutoff and the flow wanes the static regions become more significant and lead to distinct fingers separated by grain-free regions (figures 3 b, 4 b). In order to bring material to rest and prevent 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 (Johnson et al. 2012; 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 3 a). 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.

Figure 15. Numerical simulations of depth-averaged concentration $\bar{\unicode[STIX]{x1D719}}$ at time $t=7.5~\text{s}$ for different grid resolutions and domain widths $L_{y}$ . Black dotted lines denote the maximum and minimum front position, as defined by (6.10) and (6.11) respectively. The final results are not identical, but the width and downslope extent of the fingers is similar for all runs. Other parameters used are the same as in figures 12–14.

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 $\bar{\unicode[STIX]{x1D719}}$ is plotted at time $t=7.5~\text{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 $\unicode[STIX]{x1D6EC}$ 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~\text{s}$ . However, running the simulations for longer can lead to merging events, making precise definitions more difficult and meaning that $\unicode[STIX]{x1D6EC}$ will also evolve over time.

Figure 16. The effect of changing the grid resolution on (a) the mean finger wavelength $\unicode[STIX]{x1D6EC}$ at time $t=7.5~\text{s}$ and (b) the front elongation $x_{l}$ (given by (6.10)–(6.12)) at times $t=2.5~\text{s}$ , $t=5.0~\text{s}$ and $t=7.5~\text{s}$ . Both tend to roughly constant values for large enough numbers of grid points.

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

(6.10) $$\begin{eqnarray}\displaystyle & \displaystyle x_{f}^{+}(t)=\max \{x:h(x,y,t)<h_{min}\}, & \displaystyle\end{eqnarray}$$

(6.11) $$\begin{eqnarray}\displaystyle & \displaystyle x_{f}^{-}(t)=\min \{x:|\bar{\unicode[STIX]{x1D719}}(x,y,t)-\bar{\unicode[STIX]{x1D719}}_{0}|>\bar{\unicode[STIX]{x1D719}}_{0}/2\}, & \displaystyle\end{eqnarray}$$

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) $$\begin{eqnarray}\displaystyle x_{l}(t)=x_{f}^{+}(t)-x_{f}^{-}(t). & & \displaystyle\end{eqnarray}$$

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.

Figure 17. The effect of changing the viscosity coefficient $\unicode[STIX]{x1D708}$ on (a) the mean finger wavelength $\unicode[STIX]{x1D6EC}$ at a given time $t=7.5~\text{s}$ and (b) the front elongation $x_{l}$ . The grid resolution is fixed at $N_{x}/L_{x}=N_{y}/L_{y}=500~\text{m}^{-1}$ .

Figure 18. Plots of minimum and maximum front position, $x_{f}^{-}$ and $x_{f}^{+}$ , and finger length $x_{l}$ as functions of time for different values of $\unicode[STIX]{x1D708}$ . All travel at approximately constant velocities until the onset of finger formation, which occurs at earlier times for smaller viscosities.

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 $\unicode[STIX]{x1D708}=0.01~\text{m}^{3/2}~\text{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 $\unicode[STIX]{x1D708}$ ranging from 0.005 to 0.05 to investigate what effect this has on the physical characteristics of the fingers. The results are shown on figure 17, where it can be seen that increasing $\unicode[STIX]{x1D708}$ 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 to a sudden increase in the rate of elongation of the frontal region. This transition point happens at earlier times for lower values of $\unicode[STIX]{x1D708}$ , 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.

7 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 Gray & Edwards (2014) to incorporate higher-order terms into their depth-averaged $\unicode[STIX]{x1D707}(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 $\unicode[STIX]{x1D707}(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. 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 $\unicode[STIX]{x1D708}$ . 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 $\unicode[STIX]{x1D708}$ is currently treated as a free parameter, since fingers only form at slope angles where the depth-averaged $\unicode[STIX]{x1D707}(I)$ -rheology of Gray & Edwards (2014) needs to be regularised. Choosing $\unicode[STIX]{x1D708}$ 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.