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

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

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

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

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

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 (Reference Gray and Edwards2014). 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,

The linearised governing equations then become where the simplified notation

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

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

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

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

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

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

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,

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 (Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012) model, although they used the original exponential form of the basal friction law (Pouliquen Reference Pouliquen1999a ) 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. (Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012) 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 Reference Joseph and Saut1990). Although the model was well posed almost everywhere, numerical simulations of the fingering instability showed that the width of the fingers was grid dependent. This was because the computations always encompassed a line of points where the depth-averaged concentration was equal to the critical concentration for ill posedness, and refining the computational domain introduced increasingly unstable small wavelength perturbations in these regions. The critical regime can be related to the characteristics of the inviscid equations, which are given in $(x,t)$ space by the lines

where $\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

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

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

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

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

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. (Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012) and means that the inviscid model is ill posed in the sense of Joseph & Saut (Reference Joseph and Saut1990).

5.4 Viscous high wavenumber asymptotics

Returning to the viscous equations (finite $R>0$ ), the leading-order behaviour of the three roots are

for $k\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

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.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 Reference Gray and Edwards2014; Edwards & Gray Reference Edwards and Gray2015). 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. (Reference Baker, Barker and Gray2016) examined a similar problem for monodisperse material and generalised the one-dimensional depth-averaged $\unicode[STIX]{x1D707}(I)$ -rheology of Gray & Edwards (Reference Gray and Edwards2014) 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

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

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

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

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. Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012). 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. (Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012) 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 (Reference Kurganov and Tadmor2000). 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 Reference Medovikov1998) 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. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012), granular roll waves (Gray & Edwards Reference Gray and Edwards2014; Razis et al. Reference Razis, Edwards, Gray and van der Weele2014) and erosion–deposition waves (Edwards & Gray Reference Edwards and Gray2015). 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 (Reference Noelle2000) is chosen. In order to utilise the numerical method of Kurganov & Tadmor (Reference Kurganov and Tadmor2000), 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

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

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

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 Reference Thornton and Gray2008; Gray & Ancey Reference Gray and Ancey2009; Johnson et al. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012; Gajjar et al. Reference Gajjar, van der Vaart, Thornton, Johnson, Ancey and Gray2016) which in the depth-averaged segregation theory (Gray & Kokelaar Reference Gray and Kokelaar2010a ,Reference Gray and Kokelaar 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 Reference Gray and Kokelaar2010a ). A steadily travelling base state must include a mechanism to counteract this accumulation, which Gray & Ancey (Reference Gray and Ancey2009) 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. (Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012) 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 Reference Gray and Edwards2014), which can be found analytically (Gray & Ancey Reference Gray and Ancey2009) or numerically (Pouliquen Reference Pouliquen1999b ) 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 (Reference Pouliquen and Forterre2002) 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. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012; Kokelaar et al. Reference Kokelaar, Graham, Gray and Vallance2014). 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. (Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012) 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. (Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012). 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

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

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.