2.1. Governing equations

Granular material is supplied from a point source onto a rough plane that is inclined at an angle $\zeta$ to the horizontal. A Cartesian coordinate system $Oxyz$ is defined with the origin $O$ located at the top centre of the plane. The $x$-axis points downslope, the $y$-axis points across the slope and the $z$-axis is the upward pointing normal. The unit normals are $\boldsymbol {e}_x$, $\boldsymbol {e}_y$ and $\boldsymbol {e}_z$, and the velocity $\boldsymbol {u}$ has components $(u,v,w)$ in each of these three directions, respectively. The equations are integrated through the flow depth $h$, measured normal to the plane, and the depth-averaged velocity $\bar {\boldsymbol {u}}$ has components

(2.1a,b)\begin{equation} \bar{u} = \frac{1}{h}\int_0^h u\,\mathrm{d}z, \quad \bar{v} = \frac{1}{h}\int_0^h v\,\mathrm{d}z, \end{equation}

in the downslope and cross-slope directions, respectively. The granular material is assumed to be incompressible with a constant-uniform bulk density $\rho$. It follows that the depth-averaged mass and momentum balance equations (Gray & Edwards Reference Gray and Edwards2014; Baker et al. Reference Baker, Barker and Gray2016a; Edwards et al. Reference Edwards, Viroulet, Kokelaar and Gray2017) are

(2.2)$$\begin{gather} \dfrac{\partial{h}}{\partial{t}} + \mbox{div}\left(h \bar{\boldsymbol{u}}\right) = 0, \end{gather}$$

(2.3)$$\begin{gather} \dfrac{\partial}{\partial{t}}\left(h\bar{\boldsymbol{u}}\right) + \mbox{div}\left(h \bar{\boldsymbol{u}} \otimes \bar{\boldsymbol{u}}\right) + \mbox{grad}\left(\frac{1}{2} h^2 g\cos{\zeta}\right) = hg\boldsymbol{S} \cos{\zeta} + \hbox{div}\left(\nu h^{3/2}\bar{\boldsymbol{D}}\right), \end{gather}$$

where $\otimes$ is the dyadic product and $g$ is the constant of gravitational acceleration. Note that the two-dimensional divergence and gradient operators, $\hbox {div}$ and $\hbox {grad}$, are used here to distinguish them from their three-dimensional counterparts, which are written in terms of the $\boldsymbol {\nabla }$ operator in § 5. In the source term on the right-hand side of (2.3), the non-dimensional net acceleration,

(2.4)\begin{equation} \boldsymbol{S} = \tan{\zeta}\,\boldsymbol{e_x} - \mu_b\boldsymbol{e}, \end{equation}

consists of the component of gravity pulling the avalanche downslope along the direction of the unit vector $\boldsymbol {e}_x$ and the effective basal friction $\mu _b$. The hysteretic effective basal friction law for $\mu _b$ is introduced in § 2.2. It is a generalization of Pouliquen's (Reference Pouliquen1999a) dynamic friction law, which is directly linked to the $\mu (I)$ rheology for granular flows (GDR-Midi 2004; Jop et al. Reference Jop, Forterre and Pouliquen2005, Reference Jop, Forterre and Pouliquen2006; Gray & Edwards Reference Gray and Edwards2014; Baker et al. Reference Baker, Barker and Gray2016a). If the avalanche is in motion, the direction of the friction, $\boldsymbol {e}$, opposes the motion, and when the grains are static, the friction opposes the downslope component of gravitational acceleration and the depth-averaged pressure gradient, which implies that

(2.5)\begin{equation} \boldsymbol{e} = \begin{cases} \dfrac{\bar{\boldsymbol{u}}}{|\bar{\boldsymbol{u}}|} & \text{if $|\bar{\boldsymbol{u}}|>0$},\\[10pt] \dfrac{\tan\zeta\,\boldsymbol{e}_x-\hbox{grad}\,h}{|\tan\zeta\,\boldsymbol{e}_x-\hbox{grad}\,h|} & \text{if $|\bar{\boldsymbol{u}}|=0$}. \end{cases} \end{equation}

The final term on the right-hand side of (2.3) arises from the inclusion of depth-averaged in-plane deviatoric stresses, which are neglected in most theories. Here, however, they play a crucial role in selecting a unique solution for the steady-state channel width and height. Gray & Edwards (Reference Gray and Edwards2014) and Baker et al. (Reference Baker, Barker and Gray2016a) derived the specific form for this second-order-gradient depth-averaged viscous term from the $\mu (I)$ rheology. It consists of a depth-averaged viscosity $\nu h^{1/2}/2$ that multiplies the depth-integrated strain-rate tensor

(2.6)\begin{equation} \bar{\boldsymbol{D}} = \tfrac{1}{2}\left(\bar{\boldsymbol{L}} + \bar{\boldsymbol{L}}^{T}\right), \end{equation}

where $\bar {\boldsymbol {L}}=\hbox {grad}\,\bar {\boldsymbol {u}}$ is the two-dimensional gradient of the depth-averaged velocity. The coefficient $\nu$ is determined from the $\mu (I)$ friction law (GDR-Midi 2004; Jop et al. Reference Jop, Forterre and Pouliquen2006) and an explicit formula will be given for it in § 2.3.

2.2. The hysteretic friction law

This paper uses the hysteretic friction law of Edwards et al. (Reference Edwards, Russell, Johnson and Gray2019), which is a generalization of Pouliquen & Forterre's (Reference Pouliquen and Forterre2002) original law. It consists of three cases that are known as the dynamic, intermediate and static friction regimes

(2.7)\begin{equation} \mu_b(h/d,Fr)=\left\{ \begin{array}{@{}ll@{}} \mu_{D}, & Fr \geqslant \beta_*, \\ \mu_{I}, & 0 < Fr \leqslant \beta_*,\\ \mu_{S}, & Fr=0. \end{array}\right. \end{equation}

These regimes are defined respectively by the functions

(2.8)$$\begin{gather} \mu_{D}=\mu_1 + \frac{\mu_2 - \mu_1}{1 + h \beta/(\mathcal{L}d(Fr+\varGamma))}, \end{gather}$$

(2.9)$$\begin{gather}\mu_{I}=\left(\dfrac{Fr}{\beta_*}\right)\left(\mu_1 + \frac{\mu_2 - \mu_1}{1 + h \beta/(\mathcal{L}d(\beta_*+\varGamma))} - \mu_3-\frac{\mu_2 - \mu_1}{1 + h /(\mathcal{L}d)}\right) + \mu_3 + \frac{\mu_2 - \mu_1}{1 + h /(\mathcal{L}d)}, \end{gather}$$

(2.10)$$\begin{gather}\mu_{S}=\textrm{min}\left( \mu_3 + \frac{\mu_2 - \mu_1}{1 + h /(\mathcal{L}d)}, |\tan \zeta\,\boldsymbol{e}_x-\hbox{grad}\,h| \right), \end{gather}$$

where $\mu _1=\tan \zeta _1$, $\mu _2=\tan \zeta _2$ and $\mu _3=\tan \zeta _3$ are the tangents of the friction angles $\zeta _1$, $\zeta _2$ and $\zeta _3$. These parameters are determined from experimental fits to the deposit depth $h_{stop}(\zeta )$ and initiation thickness $h_{start}(\zeta )$ as a function of $\zeta$. The fitting functions contain a frictional length scale $\mathscr {L}$, which Forterre & Pouliquen (Reference Forterre and Pouliquen2003) showed scales with the grain size $d$. This grain-size dependence has explicitly been included in (2.8)–(2.10) by defining the non-dimensional friction length scale $\mathcal {L}=\mathscr {L}/d$. The regime boundaries are determined by the Froude number $Fr$, which is the ratio of the flow speed to the gravity wave speed

(2.11)\begin{equation} Fr=\frac{|\bar{\boldsymbol{u}}|}{\sqrt{gh\cos\zeta}}. \end{equation}

If $Fr$ is above the minimum Froude number $\beta _*$ for steady-uniform flow, then the flow is in the dynamic regime, if the Froude number equals zero, it is in the static regime and if $Fr\in (0,\beta _*)$, then it is in the intermediate regime. The remaining non-dimensional constants $\beta$ and $\varGamma$ are determined from a pair of best fits to the steady-uniform-flow law for angular particles (Forterre & Pouliquen Reference Forterre and Pouliquen2003)

(2.12)\begin{equation} Fr = \beta \frac{h}{h_{stop}} - \varGamma. \end{equation}

Edwards et al. (Reference Edwards, Russell, Johnson and Gray2019) assumed that the minimum observable steady-uniform flow occurs at the same Froude number $Fr=\beta _*$ at all slope angles. It follows from (2.12) that the minimum steady-uniform-flow depth

(2.13)\begin{equation} h_* = \left(\frac{\beta_*+\varGamma}{\beta}\right) h_{stop}(\zeta), \end{equation}

which implies that $h_*$ is a multiple of $h_{stop}$. The frictional parameters must be chosen so that $h_{stop}< h_*< h_{start}$ for the hysteretic friction law to be well defined. A necessary, but not sufficient condition is that $\mu _1<\mu _3<\mu _2$. This paper adopts the parameters for sand used by Rocha et al. (Reference Rocha, Johnson and Gray2019), except for $\zeta _3$, which is increased from $31^\circ$ to $33^\circ$. This increases the maximum static friction and hence the thickness range $(h\in [h_{stop},h_{start}])$ with coexisting flowing and static states. Increasing $\mu _3$ gives the levees slightly greater stability than in the simulations of Rocha et al. (Reference Rocha, Johnson and Gray2019), which suppresses the periodically upslope propagating erosion-deposition waves seen in their movie 5. All the parameters (see table 1) are tightly constrained by existing experiments in the literature and are typical of those for sand flowing on a bed of glass beads. Rocha et al. (Reference Rocha, Johnson and Gray2019) also gave parameters for flows of glass beads, which also form levees even though $\varGamma =0$ (Félix & Thomas Reference Félix and Thomas2004). The values of $\mu _1$ and $\mu _3$ are, however, very close to one another, which leads to weak non-monotonic dependence of the friction $\mu$ on the Froude number $Fr$, as shown in figure 6(b) of Rocha et al. (Reference Rocha, Johnson and Gray2019). This leads to very weak levees. Indeed, Deboeuf et al. (Reference Deboeuf, Lajeunesse, Dauchot and Andreotti2006) have observed experimentally that levees of glass beads slowly creep outwards over very long periods of time, which is not the case for sand.

Table 1. Material properties for flows of sand with a mean diameter of $0.45$ mm. The values of $\zeta _1$, $\zeta _2$ and $\mathcal {L}$ were calculated from the $h_{stop}$ curve of Takagi et al. (Reference Takagi, McElwaine and Huppert2011) by Rocha et al. (Reference Rocha, Johnson and Gray2019). The values of $\beta$ and $\varGamma$ differ from those in Pouliquen & Forterre (Reference Pouliquen and Forterre2002) and Forterre & Pouliquen (Reference Forterre and Pouliquen2003) to account for the factor $\sqrt {\cos \zeta }$ in the Froude number (2.11). All the values are the same as those used by Rocha et al. (Reference Rocha, Johnson and Gray2019) except for $\zeta _3$, which is increased from $31^\circ$ to $33^\circ$ to give the sand slightly more stability in the hysteretic flow regime.

2.3. The coefficient $\nu$

Gray & Edwards (Reference Gray and Edwards2014) and Baker et al. (Reference Baker, Barker and Gray2016a) derived the depth-averaged viscous terms in (2.3) from the $\mu (I)$ rheology (GDR-Midi 2004; Jop et al. Reference Jop, Forterre and Pouliquen2006). In this equation, the coefficient

(2.14)\begin{equation} \nu = \frac{2}{9} \frac{\mathcal{L}d\sqrt{g}}{\beta} \frac{\sin{\zeta}}{\sqrt{\cos{\zeta}}} \left( \frac{\mu_2 - \tan{\zeta}}{\tan{\zeta} - \mu_1} \right), \end{equation}

is explicitly determined during the integration process in terms of the existing friction parameters, summarized in table 1, and the slope angle $\zeta$. Equation (2.14) is valid for slope inclinations in the range $\zeta \in [\zeta _1,\zeta _2]$, when steady-uniform flows are predicted to exist by the dynamic frictional law (2.8). In reality, the lowest observable steady-uniform flow occurs at $Fr=\beta _*$ and one might anticipate that the form of $\nu$ might change when the flow is in the intermediate (2.9) and static (2.10) regimes. However, in the computations presented in this paper, (2.14) is assumed to apply to all regimes.

2.4. Particle size and gravitational dependence

The governing equations can be made non-dimensional by introducing scalings based on a typical particle diameter $d$ and gravity $g$

(2.15a–e)\begin{equation} \left.\begin{array}{@{}c@{}} (x,y,z,h)=d\,(\tilde{x}, \tilde{y}, \tilde{z}, \tilde{h}),\quad t=\sqrt{d/g}\,\tilde{t},\\ \bar{\boldsymbol{u}}=\sqrt{gd}\,\tilde{\bar{\boldsymbol{u}}},\quad \bar{\boldsymbol{D}}=\sqrt{g/d}\,\tilde{\bar{\boldsymbol{D}}},\quad \nu=d\sqrt{g}\,\tilde{\nu}, \end{array}\right\} \end{equation}

where the tilde is used to indicate a non-dimensional variable. It follows that the non-dimensional form of the mass and momentum balances (2.2) and (2.3) are

(2.16)$$\begin{gather} \dfrac{\partial{\tilde{h}}}{\partial{\tilde{t}}} + \widetilde{\hbox{div}}\left(\tilde{h} \tilde{\bar{\boldsymbol{u}}}\right) = 0, \end{gather}$$

(2.17)$$\begin{gather} \dfrac{\partial}{\partial{\tilde{t}}}\left(\tilde{h}\tilde{\bar{\boldsymbol{u}}}\right) + \widetilde{\hbox{div}}\left(\tilde{h} \tilde{\bar{\boldsymbol{u}}} \otimes \tilde{\bar{\boldsymbol{u}}}\right) + \widetilde{\hbox{grad}}\left(\frac{1}{2} \tilde{h}^2\cos{\zeta}\right) = \tilde{h}\boldsymbol{S} \cos{\zeta} + \widetilde{\hbox{div}}\left(\tilde{\nu} \tilde{h}^{3/2}\tilde{\bar{\boldsymbol{D}}}\right), \end{gather}$$

respectively, where $\widetilde {\hbox {div}}$ and $\widetilde {\hbox {grad}}$ are the non-dimensionalized divergence and gradient operators and the non-dimensional coefficient in the viscous term

(2.18)\begin{equation} \tilde\nu = \frac{2}{9} \frac{\mathcal{L}}{\beta} \frac{\sin{\zeta}}{\sqrt{\cos{\zeta}}} \left( \frac{\mu_2 - \tan{\zeta}}{\tan{\zeta} - \mu_1} \right). \end{equation}

Moreover, the non-dimensional net acceleration $\boldsymbol {S}$, the unit vector $\boldsymbol {e}$ and the friction $\mu _b$, defined in (2.4), (2.5) and (2.7), respectively, are already written in dimensionless form, since $\bar {\boldsymbol {u}}/|\bar {\boldsymbol {u}}|=\tilde {\bar {\boldsymbol {u}}}/|\tilde {\bar {\boldsymbol {u}}}|$, $\hbox {grad}\,h=\widetilde {\hbox {grad}}\,\tilde {h}$ and $\mu _b(h/d,Fr)=\mu _b(\tilde {h},Fr)$. The non-dimensional system of equations is therefore independent of any additional non-dimensional groups involving $d$ or $g$ dependence. These non-dimensionalized equations are therefore scale and gravity independent. This remarkable result is due to Kokelaar et al. (Reference Kokelaar, Bahia, Joy, Viroulet and Gray2017), and implies that provided all the length scales (height, width and length) scale up with the particle size, the results will be the same. A self-channelized geophysical flow of large boulders, such as that shown in figure 2, may therefore be exactly equivalent to a small-scale experimental flow involving millimetre-sized grains (Rocha et al. Reference Rocha, Johnson and Gray2019). This scaling on the particle size relies on the fact that the frictional length scale $\mathscr {L}=\mathcal {L}d$ scales linearly on $d$. Forterre & Pouliquen (Reference Forterre and Pouliquen2003) showed that $\mathscr {L}=1.65d$ for glass beads and $\mathscr {L}=2.03d$ for sand. In this paper, the non-dimensional frictional length $\mathcal {L}=\mathscr {L}/d$ is assumed to be equal to $2$ (table 1).

Kokelaar et al. (Reference Kokelaar, Bahia, Joy, Viroulet and Gray2017) applied these ideas to large-scale dry granular flows on the Moon, which occur during crater wall collapses. The reduced gravity on the Moon implies that the velocities will be slower than an equivalent sized flow on Earth and the flows would therefore have taken correspondingly longer to run out. However, Kokelaar et al.'s (Reference Kokelaar, Bahia, Joy, Viroulet and Gray2017) figure 10 shows self-channelized fingered deposits in the Bessel crater in Mare Serenitatis (latitude 21.8$^\circ$N, longitude 17.9$^\circ$E) that are closely analogous to the small-scale segregation induced fingering experiments of Pouliquen et al. (Reference Pouliquen, Delour and Savage1997), Woodhouse et al. (Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012) and Baker et al. (Reference Baker, Johnson and Gray2016b). The presence of water in debris flows on Earth modifies the friction law and the scaling properties. However, the essential flow dynamics and morphological features remain closely similar as demonstrated by the USGS debris-flow-flume experiments (Johnson et al. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012).

3.1. Continuous release of grains on an inclined plane

A numerical simulation is first performed to capture the continuous release of sand from a point source onto a rough plane inclined at $\zeta =32^{\circ }$. The values of the critical thicknesses $h_{stop}$, $h_*$ and $h_{start}$ as well as $\nu$ are therefore constant, and their specific values are summarized in table 2. The computational domain consists of a rectangle covering the area $0\ {\rm m} \leqslant x \leqslant 2.5$ m and $-10\ {\rm cm} \leqslant y \leqslant 10$ cm, discretized to $2500\times 200$ finite volume cells (1 grid point per mm in both directions). The conservation equations (2.2) and (2.3) are solved numerically with a semi-discrete non-oscillatory central (NOC) scheme, which uses a generalized minmod limiter with $\theta =2$ (Kurganov & Tadmor Reference Kurganov and Tadmor2000). The time-stepping is performed with a second-order Runge–Kutta method, with the step size determined by a CFL (Courant–Friedrichs–Lewy) number of $0.225$ and limited to a maximum of $\Delta t = 10^{-4}$ s (LeVeque Reference LeVeque2002) to minimize the creep in the levees.

Table 2. Critical layer thicknesses $h_{stop}$, $h_*$ and $h_{start}$ as well as the coefficient $\nu$ (in the depth-averaged viscosity $\nu h^{1/2}/2$) for the material properties for sand in table 1 and the fixed slope angle of $\zeta =32.0^{\circ }$ used in the numerical computations.

The initial conditions at $t=0$ are $h=h_0=0.01$ mm (around 50 times smaller than the mean particle diameter) and $\bar {\boldsymbol {u}}=\boldsymbol {0}\ {\rm m}\ {\rm s}^{-1}$ everywhere. This thin static layer is required to mitigate numerical errors caused by the degeneracy of the governing equations at $h=0$. Following Rocha et al. (Reference Rocha, Johnson and Gray2019), a source term $S_{inflow}$ is included on the right-hand side of (2.2) to give

(3.1)\begin{equation} \dfrac{\partial{h}}{\partial{t}} + \hbox{div}\left(h \bar{\boldsymbol{u}}\right) = S_{inflow}, \end{equation}

which allows grains to be added in a small circular region of radius $r_0=2.5$ cm centred at $(x_0,y_0)=(15$ cm$,0)$ at a flow rate of $Q_m = 130\ {\rm g}\ {\rm s}^{-1}$, i.e.

(3.2)\begin{equation} S_{inflow}=\left\{\begin{array}{@{}ll@{}} \displaystyle\dfrac{3 Q_M}{{\rm \pi} \rho r_0^6} \left(r_0^2 - \left((x-x_0)^2 + y^2\right) \right)^2, & (x-x_0)^2 + y^2 \leqslant r_0^2,\\[10pt] 0, & (x-x_0)^2 + y^2 > r_0^2. \end{array}\right. \end{equation}

The inflowing particles move downslope and some of them eventually reach the downstream boundary $x=2.5$ m, where a free outflow condition is imposed by linear extrapolation of the values of $h$ and $h\bar {u}$ from the final two columns of interior cells. The conditions $h=h_0$ and $\bar {\boldsymbol {u}}=\boldsymbol {0}\ {\rm m}\ {\rm s}^{-1}$ are trivially satisfied at the top $x=0$ and sides $y= \pm 10$ cm of the domain where no flow reaches.

The results of this numerical simulation are plotted in figure 3 at times $t=7$, $14$, $21$ and $28$ s. The continuous release of grains is shown to quickly form a self-channelized flow with static levees. In the central flowing channel, which is of approximately constant thickness $9$ mm and uniform width $16$ cm, material is transported towards the flow front. When it reaches the front, the lateral confinement is lost and the flow spreads out, rapidly coming to rest and building up new sections of the levees just behind the flow front. The levees are approximately $2.5$ cm wide and have a thickness that starts at the flowing channel depth and decreases to zero. As a result, the front in figure 3 appears to propagate steadily downslope at constant velocity $u_F \simeq 0.07\ {\rm m}\ {\rm s}^{-1}$, suggesting that a two-dimensional travelling wave is generated.

Figure 3. Surface plots of the flow thickness $h$ for the numerical simulation of a release of particles at a flow rate of $Q_m = 130\ {\rm g}\ {\rm s}^{-1}$ on a plane inclined at $\zeta =32^{\circ }$, coloured by the downslope component of depth-averaged velocity $\bar {u}$ and shown at times (a) $t = 7$, (b) $14$, (c) $21$ and (d) $28$ s. The filled black circles indicate the non-zero source region where $S_{inflow}$ is given by (3.2) and the solid horizontal white lines indicate cross-slope flow thickness, mimicking the displacement of an experimental laser line. The online supplementary movie 1 available at https://doi.org/10.1017/jfm.2022.1089 shows the time dependent evolution of the flow.

3.2. Depth-averaged travelling frame simulation

For the travelling wave to be fully realised, it is useful to perform simulations in a frame moving at the front speed $u_F$. Using the change of coordinates

(3.3a,b)\begin{equation} \tau=t-t_0,\quad \xi = x - u_F t, \end{equation}

the mass and momentum conservation laws (2.2) and (2.3) can be written

(3.4)$$\begin{gather} \dfrac{\partial{h}}{\partial{\tau}} + \hbox{div}'\left(h \bar{\boldsymbol{u}}^{\prime}\right) = S'_{inflow}, \end{gather}$$

(3.5)$$\begin{gather} \dfrac{\partial}{\partial{\tau}}\left(h\boldsymbol{\bar{u}}^{\prime}\right) + \hbox{div}' \left(h \bar{\boldsymbol{u}}^{\prime} \otimes \bar{\boldsymbol{u}}^{\prime}\right) + \hbox{grad}'\left(\frac{1}{2} h^2 g\cos{\zeta}\right) = hg\boldsymbol{S} \cos{\zeta} + \hbox{div}'\left(\nu h^{3/2}\overline{\boldsymbol{D}^\prime}\right), \end{gather}$$

where $\hbox {div}'$ and $\hbox {grad}'$ are now the two-dimensional divergence and gradient operators in the moving frame $(\xi,y)$. The depth-averaged velocity in the moving frame is equal to $\bar {\boldsymbol {u}}^{\prime } = \bar {\boldsymbol {u}}-u_F\boldsymbol {e}_x$. The strain-rate tensor is unchanged in the moving frame, i.e. $\bar {\boldsymbol {D}}^\prime =\bar {\boldsymbol {D}}$. Note that the elimination of $\bar {\boldsymbol {u}}$ in the transformed momentum balance (3.5) relies on subtracting $u_F$ times the moving frame mass balance (3.4) from the transformed downslope momentum balance.

Since (3.4) and (3.5) have the same structure as the original mass and momentum balance equations (2.2) and (2.3), the numerical method described in § 3.1 can be used to compute the flow in the moving frame. This time, the computational domain is a rectangle covering the area $-0.8\ {\rm m} \leqslant \xi \leqslant 0.8$ m and $-10\ {\rm cm} \leqslant y \leqslant 10$ cm, discretized to $1600\times 200$ finite volume cells (with the same resolution of 1 grid point per mm in both directions as the continuous release simulation). The initial conditions at $\tau = 0$ are the $t = t_0 = 28$ s state of the continuous release simulation, which is shown in figure 3(d). At the downstream boundary, at $\xi = 0.8$ m, the precursor layer enters the domain with thickness $h=h_0$ and velocity $\bar {\boldsymbol {u}}^\prime =-u_F\boldsymbol {e}_x$, which is equivalent to zero velocity in the stationary frame. This condition is also trivially satisfied at the sides $y= \pm 10$ cm.

The upstream boundary condition is more complex. Physically, this is because some of the grains in the central channel are moving downslope towards the front, and therefore must enter the domain at $\xi =-0.8$ m, whilst other grains at the sides and base of the central channel, as well as in the static levees, are moving slower than the front and leave the domain. This has to be accounted for in a depth-averaged model by allowing an inflow along part of the upstream boundary and an outflow along the rest. To achieve this, a novel procedure has been developed. First, the net mass flux leaving the domain across the upstream boundary is calculated by evaluating the integral

(3.6)\begin{equation} Q_{out} ={-}\rho\int_{{-}W/2}^{{+}W/2} h\bar{u}^{\prime} \rvert_{\xi={-}0.8\ {\rm m}}\, \mathrm{d}y. \end{equation}

At steady state, the inflowing and outflowing grains should be in exact balance and $Q_{out}$ should equal zero. However, the initial state is not yet quite in steady state and therefore $Q_{out}$ is non-zero. To correct for the net loss/gain of particles across $\xi =-0.8$ m, a mass source term (similar to that applied in (3.1)) is applied across the first row of interior grid cells. The new source term in the moving frame simulations is therefore

(3.7)\begin{equation} S'_{inflow} = \dfrac{Q_{out}h\rvert_{\xi={-}0.8\ \mathrm{m}}}{\rho\Delta x \displaystyle\int\nolimits_{{-}W/2}^{{+}W/2} h\rvert_{\xi={-}0.8\ \mathrm{m}}\,\mathrm{d} y} \end{equation}

inside the first grid cell and zero otherwise, where $\Delta x = 10^{-3}$ m is the length of a grid cell in the downslope direction. The inflow is weighted so that more mass is supplied in thicker parts of the flow. It follows that when (3.7) is integrated over the domain $(-0.8,-0.8+\Delta x~\mathrm {m})\times (-W,W)$, it resupplies the material lost through the upstream boundary. An iterative procedure is required to refine the initial guess of $u_F=0.07065\ {\rm m}\ {\rm s}^{-1}$ to ensure convergence towards a steady state. Figure 4 shows the results of a numerical simulation with $u_F=0.07053\ {\rm m}\ {\rm s}^{-1}$. Figure 4(a,b) show that the initial thickness is almost indistinguishable from the final thickness 100 s later, implying that by $t=28$ s, the simulations in figure 3(d) were already close to steady state. The temporal evolution of $Q_{out}$ is shown in figure 4(c). By $\tau =100$ s, the simulation has converged on a steady state with $Q_{out}=0$, and this is now used to reconstruct the three-dimensional velocity field near the front of a self-channelized flow with static levees.

Figure 4. Colour maps of flow thickness $h$ in a frame moving with the front speed $u_F = 0.07053\ {\rm m}\ {\rm s}^{-1}$ plotted (a) at the initial time $\tau = 0$ s (equivalent to $t=t_0=28$ s of the continuous release simulation shown in figure 3d) and (b) after $\tau = 100$ s in the travelling frame simulation. (c) Mass flux $Q_{out}$ leaving the domain as a function of time $\tau$. This tends to zero for large times confirming that a travelling wave solution is achieved.

4.1. Reconstruction of the bulk flow velocity components

The depth-averaged theory in § 2 is able to quantitatively model the formation of self-channelized flows with static levees as shown by Rocha et al. (Reference Rocha, Johnson and Gray2019) and the simulations in § 3. Following Gray & Ancey (Reference Gray and Ancey2009) and Johnson et al. (Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012), the three-dimensional velocity field $\boldsymbol {u} = (u,v,w)$ is reconstructed from $\boldsymbol {\bar {u}}$ by assuming a velocity profile through the depth of the flow. For a general function $f(z/h)$, the horizontal components of the three-dimensional velocity field parallel to the slope are related to their depth-averaged counterparts by

(4.1)\begin{equation} (u,v) = f\left(\frac{z}{h}\right)(\bar{u},\bar{v}). \end{equation}

This paper assumes a Bagnold velocity profile (Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001; GDR-Midi 2004; Jop et al. Reference Jop, Forterre and Pouliquen2005; Gray & Edwards Reference Gray and Edwards2014; Baker et al. Reference Baker, Barker and Gray2016a) which implies

(4.2)\begin{equation} f\left(\frac{z}{h}\right)=\frac{5}{3}\left[1 - \left(1-\frac{z}{h}\right)^{{3}/{2}}\right]. \end{equation}

It is possible to reconstruct the three-dimensional velocity field using other profiles. For instance, Johnson et al. (Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012) reconstructed a qualitatively similar depth-averaged velocity field $\bar {\boldsymbol {u}}$ using (i) plug flow, (ii) simple shear and (iii) linear shear with basal slip. As shown in their figure 11, this has an important affect on the surface trajectories of particles and determines the degree to which they reach the front and are overturned, or are deposited at the surface of the static levees. For plug flow, there is no overturning, which is unrealistic, and simple shear has too much overturning compared with measurements of debris flows at the USGS flume. Johnson et al. (Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012) found that a combination of linear shear and basal slip was the most realistic. This is also consistent with the Bagnold profile used here (Baker et al. Reference Baker, Johnson and Gray2016b).

Saingier, Deboeuf & Lagrée (Reference Saingier, Deboeuf and Lagrée2016) showed that when there is velocity shear, non-unity values of the velocity shape factor $\chi =\overline {u^2}/(\bar {u})^2$ in the momentum transport terms lead to the formation of a precursor layer in depth-averaged avalanche models. Since a precursor layer is not observed in experiments, they inferred that the velocity profile must become plug-like at the front, i.e. $\chi \rightarrow 1$ as $h\rightarrow 0$. Many avalanche computations therefore assume $\chi =1$. Despite this, frontal overturning is often observed in experiments (Johnson et al. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012), and also occurs in non-depth-averaged numerical simulations using the $\mu (I)$ rheology (see figures 3 and 6 of Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021). This paper therefore takes the pragmatic view. The shape factors are assumed to be unity in the depth-averaged momentum balance (2.3), but Bagnold velocity profiles are assumed everywhere, which allows the flow to overturn at the front. Support for this approach comes from Saingier et al. (Reference Saingier, Deboeuf and Lagrée2016) who showed that at low Froude numbers ($Fr<1$), the Bagnold shape factor ($\chi =5/4$) barely changes the free surface shape. Since $Fr\simeq 0.3$ in the simulations in figure 3, this is therefore a good approximation. Open questions still remain about whether the velocity profile should also be Bagnold-like in the cross-stream direction, but it is a reasonable first approximation. In addition, this paper neglects all feedback of the evolving particle-size distribution on the local velocity profile (Rognon et al. Reference Rognon, Roux, Naaim and Chevoir2007) or the basal friction. This would require a fully coupled theory such as the depth-averaged one of Baker et al. (Reference Baker, Johnson and Gray2016b) or the non-depth-averaged one of Barker et al. (Reference Barker, Rauter, Maguire, Johnson and Gray2021).

Using the Bagnold profile (4.2), it follows that the horizontal components of the three-dimensional velocity field $(u^{\prime },v^{\prime })$ in a frame moving with the speed of the front (3.3a,b) are

(4.3)$$\begin{gather} u^{\prime} ={-}u_F + \frac{5}{3}(\bar{u}^{\prime} + u_F) \left[1 - \left(1-\frac{z}{h}\right)^{{3}/{2}}\right], \end{gather}$$

(4.4)$$\begin{gather}v^{\prime} = \frac{5}{3}\bar{v}^{\prime}\left[1 - \left(1-\frac{z}{h}\right)^{{3}/{2}}\right], \end{gather}$$

respectively. Figures 5 and 6 show the three-dimensional downslope and cross-slope velocity fields reconstructed from the two-dimensional travelling wave solutions in § 3.2 at $\tau =100$ s. In these figures, panels (a–d) show longitudinal cross-sections at equally spaced positions across the width of the flow (as indicated in panel h), whereas panels (e–h) show sections across the flow at regular downslope intervals (as indicated in panel d). The central plane (panel a) looks similar to a Pouliquen (Reference Pouliquen1999b) front, with a monotonically decreasing free surface, but the downslope free-surface velocity decreases towards the front indicating that there must be lateral motion. The highest downslope velocities are attained in the thickest part of the flow in the centre of the channel, and the velocity decreases to zero towards the bottom and sides of the flow. Within the levees (panel d and panels e–g), the downslope velocity $u$ is zero in the stationary frame, but in the moving frame, it appears to move upslope with velocity $u'=-u_F\simeq -0.07\ {\rm m}\ {\rm s}^{-1}$. The static/slow-moving region is broadly consistent with that inferred experimentally by Deboeuf et al. (Reference Deboeuf, Lajeunesse, Dauchot and Andreotti2006) in their figure 4, although a direct comparison with the steady state used here is not possible, due to the partial collapse of the levee walls and drainage of the channel as their flow was brought to rest. Note that Rocha et al. (Reference Rocha, Johnson and Gray2019) showed that the same depth-averaged model used here is able to capture this drainage and collapse (see their figures 14 and 15, and movies 5 and 6). Their model could therefore be used to reconstruct the full three-dimensional temporally evolving flow field, but that is beyond the scope of this paper.

Figure 5. Colour maps of the downslope component of the three-dimensional velocity field $\boldsymbol {u}$ reconstructed from the depth-averaged travelling wave computed in § 3.2 (at $\tau =100$ s) with front speed $u_F=0.07\ {\rm m}\ {\rm s}^{-1}$. The colour maps are shown at cross-slope positions (a) $y = 0$, (b) $\pm 2$, (c) $\pm 4$ and (d) $\pm 6$ cm, and downslope positions (e) $\xi = 0.1$, ( f) $0.3$, (g) $0.5$ and (h) $0.7$ m. The relative positions of the various cross-sections are indicated by the black tick marks and corresponding labels in panels (d,h). A single colour scale is used to show both $u$ and $u'=u-u_F$. The dashed white line in panel (e) is the no-net-flow line. It separates the region above, where the particles are moving towards the front ($u^{\prime }>0$), from the region beneath, where particles are moving away from the front ($u^{\prime }<0$). The depth-averaged flow parameters are summarized in tables 1 and 2, and a Bagnold velocity profile (4.2) is used in the reconstruction.

Figure 6. Colour maps of the cross-slope component $v=v'$ of the three-dimensional velocity field $\boldsymbol {u}$ reconstructed from the depth-averaged travelling wave computed in § 3.2 (at $\tau =100$ s) with front speed $u_F=0.07\ {\rm m}\ {\rm s}^{-1}$. The colour maps are shown at cross-slope positions (a) $y = 0$, (b) $-2$, (c) $-4$ and (d) $-6$ cm, and downslope positions (e) $\xi = 0.1$, ( f) $0.3$, (g) $0.5$ and (h) $0.7$ m. The relative positions of the various cross-sections are indicated by the black tick marks and corresponding labels in panels (d) and (h). Note that cross-slope velocity is antisymmetric about the centre line $y=0$. The depth-averaged flow parameters are summarized in tables 1 and 2, and a Bagnold velocity profile (4.2) is used in the reconstruction.

The central flowing channel is of almost constant depth at $\xi =0.1$ m (panel e) indicating that the flow is close to the steady-state solutions of Rocha et al. (Reference Rocha, Johnson and Gray2019) and becomes progressively more curved as the front is approached and the lateral motion advects the incoming grains into the slower moving parts of the flow and the static levees. The cross-slope velocity (figure 6) is antisymmetric about $y=0$, indicating that the oncoming flow is pushed out laterally on either side of the centre plane with the highest cross-slope speeds achieved along the free surface at the sides of the spreading front. The no-net-flow line, where $u'=0$, is shown in figure 5(e). Material that is above the no-net-flow line moves faster than the flow front ($u'>0$) and is advected downslope, while material that is below the no-net-flow line moves slower than the front ($u'<0$) and is transported back upstream. In the fixed laboratory frame, however, all the material is either moving downslope or is stationary.

An expression for the normal component $w'$ of the three-dimensional velocity in the moving frame is derived by integrating the incompressibility condition (upon which the depth-averaged theory is based, see e.g. Gray & Edwards Reference Gray and Edwards2014; Baker et al. Reference Baker, Barker and Gray2016a) with respect to $z$, subject to the condition that $w'=0$ at $z=0$, to give

(4.5)\begin{equation} w^{\prime}(\xi,y,z) ={-}\int_0^z \left( \dfrac{\partial{u^{\prime}}}{\partial{\xi}} + \dfrac{\partial{v^{\prime}}}{\partial{y}}\right) \mathrm{d}\hat{z}, \end{equation}

where $\hat {z}$ is a dummy variable. Substituting the horizontal velocity components (4.3) and (4.4) into (4.5), and using the steady-state mass balance equation in the moving frame (3.4) to simplify the result, implies that the normal velocity

(4.6)\begin{equation} w^{\prime} = f_1\left(\frac{z}{h}\right) u_F \dfrac{\partial{h}}{\partial{\xi}} - f_2\left(\frac{z}{h}\right) h\left[\dfrac{\partial{\bar{u}^{\prime}}}{\partial{\xi}} + \dfrac{\partial{\bar{v}^{\prime}}}{\partial{y}}\right], \end{equation}

where the functions

(4.7)$$\begin{gather} f_1\left(\frac{z}{h}\right)=\frac{1}{3}\left[2 - 5\left(1-\frac{z}{h}\right)^{{3}/{2}} + 3\left(1-\frac{z}{h}\right)^{{5}/{2}}\right], \end{gather}$$

(4.8)$$\begin{gather}f_2\left(\frac{z}{h}\right)=\frac{5}{3}\left[\frac{z}{h} - \left(1-\frac{z}{h}\right)^{{3}/{2}} + \left(1-\frac{z}{h}\right)^{{5}/{2}}\right]. \end{gather}$$

The reconstructed normal component of velocity for the two-dimensional travelling wave solution computed in § 3.2 is illustrated in figure 7, using the same format as in figures 5 and 6. This indicates that the flow moves rapidly downwards along the free surface of the front and that there is no normal velocity in the levees. Just like $u'$, $w'$ is symmetric about the centreline $y=0$.

Figure 7. Colour maps of the cross-slope component $w=w'$ of the three-dimensional velocity field $\boldsymbol {u}$ reconstructed from the depth-averaged travelling wave computed in § 3.2 (at $\tau =100$ s) with front speed $u_F=0.07\ {\rm m}\ {\rm s}^{-1}$. The colour maps are shown at cross-slope positions (a) $y = 0$, (b) $\pm 2$, (c) $\pm 4$ and (d) $\pm 6$ cm, and downslope positions (e) $\xi = 0.1$, ( f) $0.3$, (g) $0.5$ and (h) $0.7$ m. The relative positions of the various cross-sections are indicated by the black tick marks and corresponding labels in panels (d,h). The depth-averaged flow parameters are summarized in tables 1 and 2, and a Bagnold velocity profile (4.2) is used in the reconstruction.

4.2. Reconstruction of the free-surface velocity and particle trajectories

Figure 8(a) shows the velocity vectors in the stationary frame superimposed on top of a colour map of the downslope component of the free-surface velocity. Sufficiently far upstream, the static levees are close to their steady-state widths and there is a pronounced velocity profile across the central channel. Material is transported to the flow front along the central channel and as it nears the front, the levees decrease in width and the flow spreads out laterally and slows down. This is qualitatively similar to large-scale debris-flow experiments (see figure 9a of Johnson et al. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012), so this relatively simple small-scale simulation captures the essential features of much larger scale geophysical flows. Figure 8(b) shows the same solution, but in the moving frame, so it is now possible to see that relative to the moving front, the static levees are moving backwards in accordance with the large-scale experiments in figure 9(b) of Johnson et al. (Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012).

Figure 8. Plots of the surface values of (a) laboratory frame velocity $\boldsymbol {u}$ (white arrows), (b) moving frame velocity $\boldsymbol {u}^{\prime }$ (white arrows), and (c) moving frame streamlines (solid and dashed white lines) on top of colour maps of the depth-averaged downslope velocity component $\bar {u}$. The surface velocity is shown to be greatest along $y=0$ and the levees are stationary in the laboratory frame of reference (a) or equivalently propagating upslope with the front speed $u_F$ in the moving frame (b). Surface material flowing down the central channel is pushed around at the lateral extents of the snout where the magnitude of the cross-slope velocity is greatest. This is depicted by the streamlines in (c), where the outermost trajectories in the $y$-direction turn back on themselves (solid white lines) whilst the central ones (dashed white lines) are turned over by the flow front.

The motion of a particle at position $\boldsymbol { \xi }_p =(\xi _p,y_p,z_p)$ that is advected by the three-dimensional velocity field $(u^{\prime }, v^{\prime }, w^{\prime })$ in the moving frame of reference is determined by

(4.9a–c)\begin{equation} \dfrac{\mathrm{d} \xi_p}{\mathrm{d} \tau} = u^{\prime}(\boldsymbol{ \xi}_p), \quad \dfrac{\mathrm{d} y_p}{\mathrm{d} \tau} = v^{\prime}(\boldsymbol{ \xi}_p), \quad \dfrac{\mathrm{d} z_p}{\mathrm{d} \tau} = w^{\prime}(\boldsymbol{ \xi}_p). \end{equation}

These three-dimensional particle paths coincide with the streamlines of the bulk flow field. Figure 8(c) shows the reconstructed free-surface particle trajectories for the travelling wave solution in § 3.2. Particles in the centre of the channel are transported to the flow front and disappear from view as the flow overturns. Particles in the central channel that start closer to the static levees are advected downslope and are either recirculated back into the more slowly moving boundary layer adjacent to the levees, or deposited at the surface of the static levees and then move backwards away from the moving front. Again, this is exactly the behaviour that was observed in Johnson et al.'s (Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012) large-scale debris-flow experiments (see their figures 6 and 8). This underlines the fact that the travelling wave solution in § 3.2 and the reconstructed velocity field produces a bulk flow that is qualitatively very similar to self-channelized geophysical flows.

5.1. The bidisperse segregation equation

The flow is assumed to be composed of a bidisperse mixture of large and small particles that occupy volume fractions $\phi ^l$ and $\phi ^s$ per unit granular volume, respectively. The volume fractions (or concentrations) necessarily sum to unity

(5.1)\begin{equation} \phi^s + \phi^l = 1. \end{equation}

During shear, the large grains tend to rise above smaller grains. The segregation is governed by the segregation–advection–diffusion equation for the concentration of small particles (Gray & Thornton Reference Gray and Thornton2005; Gray & Chugunov Reference Gray and Chugunov2006; Gray Reference Gray2018)

(5.2)\begin{equation} \dfrac{\partial{\phi^s}}{\partial{t}} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left( \phi^s \boldsymbol{u} \right) + \boldsymbol{\nabla}\boldsymbol{\cdot}\left( f_{sl}\phi^s\phi^l \frac{\boldsymbol{g}}{|\boldsymbol{g}|}\right ) = \boldsymbol{\nabla}\boldsymbol{\cdot} \left( \mathcal{D}_{sl}\boldsymbol{\nabla}\phi^s \right), \end{equation}

where $\boldsymbol {\nabla }$ is the three-dimensional gradient operator, $f_{sl}$ is the segregation-velocity magnitude, $\mathcal {D}_{sl}$ is the diffusivity of the small and large particles, and $\boldsymbol {g}=g\sin \zeta \,\boldsymbol {e}_x-g\cos \zeta \,\boldsymbol {e}_z$ is the gravitational acceleration vector expressed in slope aligned coordinates. Equation (5.2) implies that the rate of change of the small particle concentration is balanced by an advective flux due to the bulk flow, a segregation flux that shuts off when a pure region of large or small particles is formed (i.e. when $\phi ^s=0$ or $1$) and a diffusive flux that smooths out the sharp concentration shocks that can otherwise form.

5.2. Expressions for the diffusivity and the segregation-velocity magnitude

There has been considerable recent progress in determining the functional form of $f_{sl}$ and $\mathcal {D}_{sl}$. For monodisperse flows of disks, Utter & Behringer (Reference Utter and Behringer2004) showed that the self diffusivity of the particles is of the form $\mathcal {D}= \mathcal {A} \dot {\gamma } d^2$, where $\mathcal {A}$ is a non-dimensional constant of proportionality, $d$ is the grain diameter and the shear rate

(5.3)\begin{equation} \dot{\gamma}=\sqrt{\tfrac{1}{2} \mathrm{tr}(\boldsymbol{D}^2)},\quad\hbox{where}\ \boldsymbol{D}=\tfrac{1}{2}\left(\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{\rm T}\right), \end{equation}

is the strain-rate tensor. Trewhela et al. (Reference Trewhela, Ancey and Gray2021) generalized this to bidisperse flows by replacing $d$ by the volume fraction weighted local mean grain diameter

(5.4)\begin{equation} \bar{d} = \phi^l d^l + \phi^s d^s, \end{equation}

where $d^l$ and $d^s$ are the diameters of the large and small particles, respectively. It follows that the bidisperse diffusivity

(5.5)\begin{equation} \mathcal{D}_{sl}= \mathcal{A} \dot{\gamma} \bar{d}^2. \end{equation}

Refractive index matched shear box experiments (Trewhela et al. Reference Trewhela, Ancey and Gray2021) have also been used to show that for grain-size ratios

(5.6)\begin{equation} R=d^l/d^s \end{equation}

close to unity, the segregation-velocity magnitude $f_{sl}$ scales as

(5.7)\begin{equation} f_{sl}= \frac{\mathcal{B}\rho_*g\dot{\gamma} \bar{d}^2(R-1)}{\mathcal{C}\rho_*g\bar{d}+p}, \end{equation}

where $\mathcal {B}$ and $\mathcal {C}$ are non-dimensional empirical constants, $\rho _*$ is the intrinsic density of the grains and $p$ is the pressure. In the theories of Gray & Edwards (Reference Gray and Edwards2014) and Baker et al. (Reference Baker, Barker and Gray2016a), which underly the depth-averaged model used to calculate the bulk flow field in §§ 2 and 3, the pressure is lithostatic, i.e.

(5.8)\begin{equation} p=\rho_*\varPhi g (h-z)\cos\zeta, \end{equation}

where $\varPhi$ is the constant solids volume fraction (GDR-Midi 2004). Substituting (5.8) into (5.7) implies that the segregation-velocity magnitude,

(5.9)\begin{equation} f_{sl}=\frac{ \mathcal{B}\dot{\gamma} \bar{d}^2(R-1)}{\mathcal{C}\bar{d}+\varPhi (h-z)\cos\zeta}, \end{equation}

is independent of the intrinsic particle density $\rho _*$ and the constant of gravitational acceleration $g$. Equation (5.9) also implies that $f_{sl}$ is linearly dependent on the shear rate $\dot {\gamma }$ and is approximately inversely proportional to the local particle depth $h-z$. All other conditions being equal, particles at the base of the flow will therefore segregate slower than particles near the free surface. The non-dimensional constant $\mathcal {C}$ was included by Trewhela et al. (Reference Trewhela, Ancey and Gray2021) to prevent the segregation-velocity magnitude from having a singularity at $z=h$. The values of the universal non-dimensional constants $\mathcal {A}$, $\mathcal {B}$ and $\mathcal {C}$ determined by Utter & Behringer (Reference Utter and Behringer2004) and Trewhela et al. (Reference Trewhela, Ancey and Gray2021) are summarized in table 3. For moderate particle-size ratios ($2< R<4$) and low small particle concentrations, Trewhela et al. (Reference Trewhela, Ancey and Gray2021) observed individual small particle intruders could percolate more rapidly through a shearing matrix of large grains than suggested by (5.9). This enhanced segregation rate was related to the transition to spontaneous percolation at large size ratios. Such effects are neglected in this paper, as are enhanced packing effects at large size ratios, which can reduce the segregation rate as mixtures of grains pack more densely together.

Table 3. The universal constants $\mathcal {A}$, $\mathcal {B}$ and $\mathcal {C}$ in the diffusivity (5.5) and the segregation-velocity magnitude (5.9) measured by Utter & Behringer (Reference Utter and Behringer2004) and Trewhela et al. (Reference Trewhela, Ancey and Gray2021), as well as the parameter values used in § 6. Note that for the simulations presented in this paper, $\mathcal {D}_{sl}$ is assumed to be zero.

5.3. Scale dependence of the segregation equation

It is also possible to show that the segregation equation (5.2) in conjunction with the diffusivity (5.5) and the segregation-velocity magnitude (5.9) are dependent on the typical particle diameter $d$. Introducing the general set of scalings

(5.10a–f)\begin{equation} \left.\begin{array}{@{}c@{}} (x,y,z,h,\bar{d},d^l,d^s)=d\,(\tilde{x}, \tilde{y}, \tilde{z}, \tilde{h},\tilde{\bar{d}}, \tilde{d}^l, \tilde d^s),\quad t=(d/U)\,\tilde{t},\quad \boldsymbol{g}=g\,\tilde{\boldsymbol{g}},\\ (\boldsymbol{u}, f_{sl})=U(\tilde{\boldsymbol{u}}, \tilde f_{sl}),\quad \dot{\gamma}=(U/d)\,\tilde{\dot{\gamma}},\quad \mathcal{D}_{sl}=(Ud)\,\tilde{\mathcal{D}}_{sl}, \end{array}\right\} \end{equation}

where $U$ is a typical velocity magnitude, it follows that the non-dimensional segregation equation

(5.11)\begin{equation} \dfrac{\partial{\phi^s}}{\partial{\tilde{t}}} + \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\left( \phi^s \tilde{\boldsymbol{u}} \right) + \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\left( \tilde{f}_{sl}\phi^s\phi^l \frac{\tilde{\boldsymbol{g}}}{|\tilde{\boldsymbol{g}}|}\right ) = \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \left( \tilde{\mathcal{D}}_{sl}\tilde{\boldsymbol{\nabla}}\phi^s \right), \end{equation}

the non-dimensional diffusivity

(5.12)\begin{equation} \tilde{\mathcal{D}}_{sl}= \mathcal{A} \tilde{\dot{\gamma}} \tilde{\bar{d}}^2 \end{equation}

and the non-dimensional segregation velocity

(5.13)\begin{equation} \tilde{f}_{sl}=\frac{ \mathcal{B}\tilde{\dot{\gamma}} \tilde{\bar{d}}^2(R-1)}{\mathcal{C}\tilde{\bar{d}}+\varPhi (\tilde{h}-\tilde{z})\cos\zeta} \end{equation}

are independent of any new non-dimensional groups. The non-dimensional segregation equation (5.11) is therefore scale invariant, and the typical particle diameter $d$ completely scales out of the problem. When $U=\sqrt {gd}$, the scalings (5.10a–f) are consistent with those for the bulk dry flow (2.15a–e) and the whole system scales with the typical particle diameter $d$ and gravity $g$. It is this fact that explains why the self-channelized fingered deposits in the Bessel crater in Mare Serenitatis (Kokelaar et al. Reference Kokelaar, Bahia, Joy, Viroulet and Gray2017) show very similar segregation to that observed in the small-scale experiments of Pouliquen et al. (Reference Pouliquen, Delour and Savage1997), Woodhouse et al. (Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012) and Baker et al. (Reference Baker, Johnson and Gray2016b), but with particles that have a much larger average grain diameter in a reduced gravity. Importantly, the typical particle size $d$ scales out of the equations even when $U$ does not scale as $\sqrt {gd}$. This implies that the solutions in §§ 6 and 7 are of direct relevance to general self-channelized flows that may have more complex scalings, e.g. the USGS debris-flow flume experiments (Johnson et al. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012) where the interstitial fluid may significantly influence the rheology (Iverson Reference Iverson1997; Denlinger & Iverson Reference Denlinger and Iverson2001; Iverson & Vallance Reference Iverson and Vallance2001; George & Iverson Reference George and Iverson2014).

5.4. Reduced segregation equation in the moving frame

Considerable insights into the effects of segregation have been obtained by constructing solutions to (5.2) in the absence of diffusion, which has the advantage that it allows the large and small particle paths to be tracked explicitly (Gray & Thornton Reference Gray and Thornton2005; Thornton, Gray & Hogg Reference Thornton, Gray and Hogg2006; Thornton & Gray Reference Thornton and Gray2008; Gray & Ancey Reference Gray and Ancey2009). In this paper, it is therefore assumed that the diffusivity

(5.14)\begin{equation} \mathcal{D}_{sl}=0. \end{equation}

In addition, the flow computed in §§ 3 and 4 is shallow, i.e. the typical length of the flow $L$ is much larger than a typical flow depth $H$. This can be exploited to simplify the segregation term. If typical downslope and cross-slope velocities are of order $U$, then incompressibility implies that slope normal velocities are of order $\varepsilon U$, where $\varepsilon =H/L$ is a small parameter. Introducing the scalings

(5.15a–c)\begin{equation} (x,y,z)=L(\hat{x}, \hat{y}, \varepsilon\hat{z}),\quad (u,v,w,f_{sl})=U(\hat{u},\hat{v}, \varepsilon\hat{w},\varepsilon\hat f_{sl}),\quad t=(L/U)\hat{t}, \end{equation}

where the segregation-velocity magnitude $f_{sl}$ is also assumed to scale as $\varepsilon U$, the non-dimensional segregation equation (5.2) becomes

(5.16) \begin{equation} \dfrac{\partial{\phi^s}}{\partial{\hat{t}}} + \dfrac{\partial}{\partial{\hat{x}}} \left(\phi^s (\hat{u}+\varepsilon\hat{f}_{sl}\phi^l\sin\zeta)\right)+ \dfrac{\partial}{\partial{\hat{y}}}\left(\phi^s \hat{v}\right) + \dfrac{\partial}{\partial{\hat{z}}}\left(\phi^s (\hat{w} -\hat{f}_{sl}\phi^l\cos\zeta)\right) = 0. \end{equation}

This shows that the downslope component of segregation (the $\sin \zeta$ term) is small compared with all the other terms and can therefore be neglected. Returning to dimensional variables and using the change of coordinates (3.3a,b), it follows that the reduced segregation equation in the moving frame is

(5.17) \begin{equation} \dfrac{\partial{\phi^s}}{\partial{\tau}} + \dfrac{\partial}{\partial\xi}\left(\phi^s u^{\prime}\right) + \dfrac{\partial}{\partial y}\left(\phi^s v^{\prime}\right) + \dfrac{\partial}{\partial z}\left(\phi^s \left(w^{\prime} - f_{sl}\phi^l \cos\zeta \right)\right) = 0, \end{equation}

where the physical domain $z=h(\xi,y)$ and the velocity field $\boldsymbol {u}^\prime (\xi,y,z)$ of the travelling wave (computed in § 3 and reconstructed in § 4) are now steady in the moving coordinates $(\xi,y,z)$.

6.1. Numerical method

The numerical solution of the advection-segregation equation (5.17) is obtained using the method of Kurganov & Tadmor (Reference Kurganov and Tadmor2000), extended straightforwardly to three dimensions in the domain $-0.4\ {\rm m} \leqslant \xi \leqslant 0.8\ {\rm m} \times -10\ {\rm cm} \leqslant y \leqslant 10\ {\rm cm} \times 0\leqslant z\leqslant 10$ mm. The flow thickness $h(\xi,y)$ and the velocity field $\boldsymbol {u}'(\xi,y,z)$ are assumed to be steady in time and given by the depth-averaged self-channelized flow simulations in § 3 and the reconstruction procedure in § 4. To accurately resolve steady solutions to (5.17), the numerical discretization must be well-balanced, in the sense that the constant solutions of the segregation equation $\phi ^s =0$ and $\phi ^s = 1$ are preserved exactly. Well-balanced discretizations often involve a balancing of discretized numerical fluxes with source terms to preserve exact solutions. Here there are no source terms, but instead, the numerically discretized form of the advection velocity field $\boldsymbol {u}'$ is required to be exactly divergence free and satisfy an appropriate discretization of $\boldsymbol {u}' \boldsymbol {\cdot } \boldsymbol {n}=0$ on $z=0$ and $z=h(\xi,y)$, so that no numerical flux enters or leaves the domain over these surfaces. This procedure is described in Appendix A.

6.2. Initial and boundary conditions

Initially, the flow is assumed to be entirely composed of small particles, i.e. $\phi ^s_0(\xi,y,z,0) = 1.0$ everywhere, and a mixture of large and small grains are fed into the domain across the upstream boundary at $\xi =-0.4$ m. There is no flux of large or small grains across $z=0$ and $z=h(\xi,y)$. By construction, $\boldsymbol {u}' \boldsymbol {\cdot } \boldsymbol {n}=0$ on $z=0, h(\xi,y)$, so this boundary condition reduces to the requirement that the small particle segregation flux normal to the boundary equals zero

(6.1)\begin{equation} \boldsymbol{F^s} \boldsymbol{\cdot} \boldsymbol{n}={-}f_{sl}\phi^s\phi^l\cos\zeta\,\boldsymbol{e}_z\boldsymbol{\cdot} \boldsymbol{n}=0,\quad \hbox{at} \ z=0,\, h(\xi,y). \end{equation}

The upstream boundary condition along $\xi =-0.4$ m is more complicated, because, in the moving frame, large and small grains flow in and out of the domain dependent on whether they lie above or below the no-net-flow line. The no-net-flow line is illustrated in figure 9. The inflowing grains are assumed to be inversely graded, with the large grains separated from the small grains beneath by a sharp interface $\eta (y)$. This is assumed to be a multiple of the inflow thickness

(6.2)\begin{equation} \eta(y)=\bar{\phi}^s_{in} h({-}0.4,y), \end{equation}

where $\bar {\phi }^s_{in}$ is a constant. Since the inflow thickness is almost constant across the flowing part of the channel, the interface height $\eta$ is almost constant, as illustrated in figure 9. The small particle concentration for the inflowing particles (above the no-net-flow line) is therefore

(6.3)\begin{equation} \phi^s({-}0.4,y,z,t)=\left\{\begin{array}{@{}ll@{}} 0, & z>\eta(y),\\ 1, & z\leqslant\eta(y). \end{array}\right. \end{equation}

For the out-flowing grains that lie below the no-net-flow line in figure 9, no boundary condition is required, and the concentration that develops within the domain is simply advected upstream out of the domain.

Figure 9. Colour map of the small particle concentration $\phi ^s$ plotted at $\xi =-0.3$ m and $\tau = 100$ s for an inflow concentration $\bar {\phi }^s_{in} = 0.8$ and $R = 1.49$ with $d^l = 0.61$ mm and $d^s = 0.41$ mm. The dashed white line is the no-net-flow line, where $u'=0$, or equivalently $u=u_F$. Above the no-net-flow line, the concentration is prescribed at the inflow at $\xi =-0.4$ m, with the large particles being sharply segregated from the small particles beneath. All the grains entering above the no-net-flow line are advected towards the flow front, segregate and eventually flow out of the domain below the no-net-flow line. The concentration below the no-net-flow line is determined by the downstream problem, and no boundary conditions are required. The white dotted lines indicate the levee-channel boundary.

Note that in figure 9, the large and small particles in the levees are not subject to further segregation because the shear rate $\dot {\gamma }=0$. Segregation does, however, continue in the mixed regions that lie in the central channel. The segregation rate is small, vanishingly so right at the margins, so it takes a very long time for the large particles to fully segregate. As a result of this and the finite domain size, some large particles flow out of the domain in an unstably stratified configuration. For example, all the mixed material in the central channel has further to segregate in figure 9. The steady-state results presented here are therefore dependent on the domain size, but only weakly so, since most of the large particles in the mixed regions would ultimately segregate to positions that lie below the no-net-flow line, so they do not affect the inflow. In reality, debris flows are only sustained in quasi-steady configurations for finite periods of time, so applying the upstream boundary condition at a finite distance from the front is more realistic than considering a semi-infinite avalanche.

Provided that the inflow is sufficiently far back, the outflowing grains on the centreline of the channel are small. It follows that $\bar {\phi }^s_{in}$ is actually the depth-averaged small particle concentration in the centre of the channel. The parameter $\bar {\phi }^s_{in}$ will therefore be used to describe the mixture composition. Note that $\bar {\phi }^s_{in}$ is not equal to the composition of the arrested mixture, because the large particles near the free surface are moving faster than the small particles beneath, and a thin layer of large particles may therefore produce a larger flux than a thicker layer of small grains (Gray & Ancey Reference Gray and Ancey2009; Wiederseiner et al. Reference Wiederseiner, Andreini, Epely-Chauvin, Moser, Monnereau, Gray and Ancey2011).

For the simulations in this section,

(6.4a,b)\begin{equation} \bar{\phi}^s_{in}=0.8,\quad \bar{\phi}^l_{in}=1-\bar{\phi}^s_{in}=0.2. \end{equation}

The volume fraction weighted mean particle diameter of the inflowing mixture,

(6.5)\begin{equation} \bar{d}_{in} = \bar{\phi}^l_{in} d^l + \bar{\phi}^s_{in} d^s, \end{equation}

is chosen to equal the particle diameter $d=0.45$ mm (table 1) used to compute the bulk flow field in § 3. The diameters of the large and small grains, $d^l$ and $d^s$, are still free to be chosen, allowing flows with different grain-size ratios $R$ to be investigated, whilst preserving $\bar {d}_{in}=0.45$ mm. For the simulations in this section, $d^l = 0.61$ mm and $d^s = 0.41$ mm, giving a grain-size ratio $R=1.49$.

6.3. Concentration profiles

The three-dimensional small-particle concentration $\phi ^s(\xi,y,z,\tau )$ is computed as a function of time $\tau$, subject to the initial and boundary conditions described in § 6.2. The temporal evolution of the depth-averaged small-particle concentration

(6.6)\begin{equation} \bar{\phi}^s(\xi,y,\tau) = \frac{1}{h} \int_0^h \phi^s(\xi,y,z,\tau)\,{\rm d}z, \end{equation}

towards a steady-state solution is shown in figure 10. Vertical downslope and cross-slope sections across the steady-state are shown in figure 11. Movie 2 shows the simultaneous temporal evolution of the small-particle concentration $\phi ^s(\xi,y,z,\tau )$ in each of the cross-sections, as well as the depth-averaged small-particle concentration $\bar {\phi }^s(\xi,y,\tau )$. The large and small particle trajectories on the centreline are superimposed on top of the small-particle concentration distribution in figure 12. Note that although the solutions are computed for $-0.4\ \mathrm {m}<\xi$, all these plots focus on the solution near the front in the range $0\leqslant \xi \le 0.8$ m.

Figure 10. Colour maps of depth-averaged small-particle concentration $\overline {\phi ^s}$ at (a) $\tau =5$, (b) $10$, (c) $15$, (d) $20$ and (e) $25$ s in the moving frame for $\bar {\phi }^s_{in} = 0.8$, $d^l = 0.61$ mm and $d^s = 0.41$ mm, which implies the grain-size ratio $R=1.49$. The white dotted lines in panel (e) indicate the levee–channel boundary. The online supplementary movie 2 shows the time dependent evolution of the flow towards steady state.

Figure 11. Colour maps of the small particle concentration $\phi ^s$ at time $\tau = 100$ s for $\bar {\phi }^s_{in} = 0.8$, $d^l = 0.61$ mm and $d^s = 0.41$ mm ($R = 1.49$) are plotted at cross-slope positions (a) $y = 0$, (b) $\pm 2$, (c) $\pm 4$ and (d) $\pm 6$ cm, and downslope positions (e) $\xi = 0.1$, ( f) $0.3$, (g) $0.5$ and (h) $0.7$ m. These locations are indicated by the black tick marks and corresponding labels in panels (d,h). The online supplementary movie 2 shows the time-dependent evolution towards the steady state in each of the sections.

Figure 12. Steady-state small-particle concentration $\phi ^s$ for $\bar {\phi }^s_{in} = 0.8$, $d^l = 0.61$ mm and $d^s = 0.41$ mm (${R = 1.49}$) showing an enlarged view ($\xi \in [0.2,0.8]$) of the breaking-size-segregation wave on the centreline $y=0$. The plot also shows the large-particle paths (solid white lines), the small-particle paths (dashed black lines) and the no-net-flow line (dashed white line).

At early times ($\tau <10$ s), the grains that flow into the domain at the upstream boundary (at $\xi =-0.4$ m) are already inversely graded and do not need to segregate vertically. The surface layer of large particles is therefore simply advected downstream by the bulk velocity field $\boldsymbol {u}'$, and consequently moves rapidly towards the front in the centre of the channel, spreading very little laterally, as shown in figure 10(a). It is only when this layer approaches the front, where the lateral velocity $v'$ becomes significant for $\xi > 0.5$ m (figure 6), that the large particles begin to spread laterally at the surface (figure 10b).

At approximately $\tau =10$ s, the leading edge of the free-surface layer of large particles reaches the front (figure 10b). Since the bulk flow wraps back on itself (figure 5), the large basal particles are transported backwards relative to the moving front. Physically, this corresponds to fast moving large grains at the free surface, flowing down over the front and entering a stationary or slow moving basal region, where the downstream velocity $u$ is much slower than the front speed $u_F$ (figure 5). Large basal particles in the centre of the channel are therefore overrun by finer grained material, and begin to segregate upwards again (see movie 2 panels (a–c) at $\tau =10\unicode{x2013}35$ s). In the moving frame, the segregating basal large particles initially appear to move backwards away from the front, because $u< u_F$ near the base, but they eventually reach a flow height where $u>u_F$ and start moving forwards towards the front again. This forms a recirculating loop. This process of large particle transport, segregation and recirculation generates a diffuse small-particle concentration region that is called a breaking size-segregation wave (Thornton & Gray Reference Thornton and Gray2008; Gray & Ancey Reference Gray and Ancey2009; Gray & Kokelaar Reference Gray and Kokelaar2010; Johnson et al. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012; Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021). In this three-dimensional flow, the breaking-size-segregation wave reaches a steady travelling state (figure 11) that moves downslope at the same speed as the front.

Thornton & Gray (Reference Thornton and Gray2008) derived an exact solution for the structure of the breaking-size-segregation wave in a two-dimensional flow field. They showed it consisted of two expansion fans and two concentration shocks that were arranged in a lens-like structure that travelled downslope at a speed that was equal to the depth average of the downslope velocity over the height of the wave. Since breaking waves do not occupy the full height of the flow, this implies that in two dimensions, they move slower than the depth-averaged speed of the flow. Above them, there is a layer of large grains that moves faster than the breaking wave, and these large particles are therefore transported forwards. In two dimensions, this allows a large-rich region to form ahead of the wave, and grow in time as increasingly more large grains accumulate there (Gray & Kokelaar Reference Gray and Kokelaar2010).

Figures 10(c–e) and 11 show that in a fully three-dimensional flow, a large-rich snout also develops, but it travels at the same speed as the breaking-size-segregation wave just behind it. For the breaking wave to move at the same speed as the flow front, there must be a physical mechanism that removes the large particles that are transported to the flow front. In two dimensions, Gray & Ancey (Reference Gray and Ancey2009) showed it was possible to do this by depositing the large frontal grains in a static layer at the base of the flow. In three dimensions, the frontal large particles can instead be removed by transporting them laterally out of plane. This is what happens in the self-channelized flows studied here.

On the centreline of the channel, $v^\prime =0$, and the small-particle concentration equation (5.17) can be written in the form

(6.7) \begin{equation} \dfrac{\partial{\phi^s}}{\partial{\tau}} + \dfrac{\partial}{\partial\xi}\left(\phi^s u^{\prime}\right) + \dfrac{\partial}{\partial z}\left(\phi^s \left(w^{\prime} - f_{sl}\phi^l \cos\zeta \right)\right) ={-}\phi^s \frac{\partial v^{\prime}}{\partial y}. \end{equation}

Since this equation is no longer dependent on cross-slope gradients of $\phi ^s$, the small-particle concentration $\phi ^s(\xi,0,z,\tau )$ can be solved for on the centreline provided $\boldsymbol {u}^\prime$ is prescribed. Johnson et al. (Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012) exploited this to compute the steady-state concentration on the centreline in a qualitatively similar flow field to that used here. This showed that the structure of the breaking-size-segregation wave was much more complicated than the two-dimensional solution of Thornton & Gray (Reference Thornton and Gray2008).

This paper goes much further by computing the three-dimensional time-dependent structure of the breaking-size-segregation wave in a self-channelized flow field. On either side of the centreline (figure 11b,c), the concentration distribution looks qualitatively similar to that on the centreline (figure 11a). However, the breaking wave extends progressively further upstream (figure 11c) as the slower moving sides of the central channel are approached. This is due to the reduced shear rate in the channel margins, which reduces the segregation rate (5.9). As a consequence, the basal large particles take longer to segregate upwards and then be transported towards the front again, which extends the length of the breaking wave. If the large particles are incorporated into the static levees, they do not segregate further at all (figure 11d–h).

In the centre of the channel, the breaking size-segregation wave has a high small-particle concentration eye, which sits just behind the flow front, on the no-net-flow line. A swirl of high large-particle concentration wraps underneath and around the eye, as shown in an enlarged view in figure 12, and progressively diminishes in strength. This characteristic breaking-wave swirl is present throughout the central channel (figure 11a–c). At the sides of the front, the bulk flow rotates the concentration distribution in the breaking-size-segregation wave through ninety degrees and transports it laterally, partially depositing it into the static levees (figure 11e–g). The levees are therefore not solely composed of large particles. Instead, there is a carapace of large particles that is wrapped around the top and bottom surfaces, and a central mixed core with higher concentrations of fine grains. This is consistent with the experimental observations of Kokelaar et al. (Reference Kokelaar, Graham, Gray and Vallance2014) who used a resin impregnation technique to make sections across the levees (see their figures 5–7).

Movie 2 shows how the small-particle concentration distribution is advected from the central channel into the levees, where it is frozen in to the static deposit. In particular, as the large particles reach the front and are advected sideways and deposited into the levee walls, a large-particle-rich deposited region propagates backwards relative to the moving front in the overhead views in figure 10(c–e). The large-rich carapace is fully developed by $\tau =25$ s. The steady-state development of the finer grained levee core takes considerably longer (approximately $\tau =50\unicode{x2013}60$ s). This is because the large particles, which form this structure, have been recirculated by the breaking-size-segregation wave and take longer to get there. Movie 2 shows this as a concentration swirl that progressively develops in the levee walls. It is important to realize that this swirl is not caused by motion of the large particles in the static levees, but just by their different arrival times in the static deposit and subsequent transport away from the front in the moving frame of reference. In the final, fully developed steady state, shown in figure 11, all the large and small particles that enter above the no-net-flow line at the upstream boundary at $\xi =-0.4$ m, exit the domain below the no-net-flow line (figure 9).

6.4. Particle trajectories

As in § 4.2 for the bulk flow, the large and small particle paths $\boldsymbol {\xi }^\nu _p =(\xi ^\nu _p,y^\nu _p,z^\nu _p)$ can be reconstructed by solving the particle-path equations

(6.8a–c)\begin{equation} \dfrac{\mathrm{d} \xi_p^\nu}{\mathrm{d} \tau} = u^\nu(\boldsymbol{ \xi}_p^\nu), \quad \dfrac{\mathrm{d} y_p^\nu}{\mathrm{d} \tau} = v^\nu(\boldsymbol{ \xi}_p^\nu), \quad \dfrac{\mathrm{d} z_p^\nu}{\mathrm{d} \tau} = w^{\nu}(\boldsymbol{ \xi}_p^\nu), \end{equation}

for each particle phase $\nu =l,s$. In the absence of diffusion, the bidisperse segregation equation (5.17) implies (Gray & Ancey Reference Gray and Ancey2009; Johnson et al. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012; Gray Reference Gray2018) that in the moving frame of reference, the large and small-particle velocity components are

(6.9a–c)\begin{gather} u^l=u',\quad v^l=v^\prime, \quad w^l = w^{\prime} + f_{sl}\phi^s \cos\zeta, \end{gather}

(6.10a–c)\begin{gather} u^s=u',\quad v^s=v^\prime, \quad w^s = w^{\prime} - f_{sl}\phi^l \cos\zeta, \end{gather}

respectively. Equations (6.8a–c) are integrated forwards in time, from the particle's initial position $\boldsymbol {\xi }^\nu _{p0} =(-0.4,y^\nu _{p0},z^\nu _{p0})$, above the no-net-flow line in figure 9, through the downstream flow field, until the particle exits the domain (below the no-net-flow line) across the $\xi =-0.4$ m plane.

The centreline $y=0$ is a somewhat pathological case, because $v^\prime =0$. The particle trajectories therefore stay on the centreline, even though mass is constantly being transported laterally out of plane by the velocity gradient term $\partial v^\prime /\partial y$ on the right-hand side of (6.7). The large and small particle trajectories on the centreline are superimposed on top of the small-particle concentration distribution in a zoomed-in view ($\xi \in [0.2,0.8]$) in figure 12. Large particles enter the flow adjacent to the free surface, and move parallel to it until they near the flow front. At the front, all the large particles are overrun, and then move backwards (in the moving frame of reference) in the slower moving flow below the no-net-flow line. As they do so, they segregate upwards into progressively faster moving parts of the flow, and when they eventually cross the no-net-flow line, they begin moving forwards towards the front again. This process of being overrun by the bulk flow, segregating upwards and recirculating forwards is repeated many times. However, unlike the two-dimensional depositing flow of Gray & Ancey (Reference Gray and Ancey2009), the large-particle paths do not form closed recirculating loops, but instead spiral in to an eye of high small-particle concentration that is centred on the no-net-flow line.

Large particles that initially move along the free surface are transported to the bottom of the flow and take longest to recirculate forwards. The free-surface large-particle trajectories therefore delineate the back of the breaking-size-segregation wave. Small particles enter the domain above the no-net-flow line, but beneath the sharp interface with the large particles (at $z=\bar {\phi }^s_{in}h$) imposed by the inflow conditions (6.2) and (6.3). All the small particle trajectories shown in figure 12 move towards the front, pass through the breaking-size segregation wave, exit beneath the no-net-flow line and then move backwards out of the domain. It is also possible (not shown) for small particles to be recirculated further upstream by the bulk flow field alone, without ever reaching the breaking-size-segregation wave. In both cases, small particles that start in the highest fastest moving part of the flow end up at or near the base, whereas those that start just above the no-net-flow line exit the domain just beneath it. However, it should be noted that the lateral velocity gradient $\partial v^\prime /\partial y$ term in (6.7) also continuously transports a mass of small particles laterally during the recirculation.

Figure 13(a) and movie 3 show selected three-dimensional particle paths for the bulk flow field, reconstructed using (4.9a–c). The particle paths are colour coded by their inflow position above the no-net-flow line in the $\xi =0$ plane to help identify them. From their initial location above the no-net-flow line, the particles are initially transported towards the flow front, but are then recirculated by the bulk flow field into regions of the flow that are moving slower than the front speed, and move backwards out of the domain below the no-net-flow line. Bulk particle paths near the centre of the channel move towards the front, where they are overrun by the bulk flow, whilst particles starting closer to the sides of the central channel are advected outwards and deposited into the levees before moving out of the frontal domain.

Figure 13. Three-dimensional particle trajectories starting and ending in the $\xi =0$ plane for (a) the bulk flow, (b) large particles and (c) small particles. The small particle concentration $\phi ^s$ is shown at the back of the flow in the $\xi =0$ plane and the area of the base of the flow at $z=0$ is filled in grey. The trajectories are coloured by their starting positions above the no-net-flow line to help identify them. These trajectory colours are used in figure 14 to show where inflowing particles starting at different locations are mapped to in the outflow. Movies 3–5 in the online supplementary material show flyby animations of the paths.

The start and end positions of the trajectories on the $\xi =0$ plane can be interpreted as a mapping. Figure 14(a) shows such a mapping using the same colour map that is used to identify the trajectories in figure 13(a). The colour map is defined on a regular grid that lies above the no-net-flow line. This is where the particles flow into the domain towards the front. The outgoing trajectories are shown on a deformed grid and colour map below the no-net-flow line. By comparing regions of the same colour above and below the no-net-flow line in figure 14(a), one can identify where inflowing material flows out of the domain. Material entering at the surface of the flow in the centre of the channel therefore flows out of the domain close to the base of the central channel, while material that starts closer to the sides of the central channel ends up in the levees. In particular, for starting positions that are sufficiently far off centre, particles starting on the surface remain at the surface and are deposited near the top of the levees.

Figure 14. Mapping between the inflow and outflow positions on the $\xi =0$ plane of the three-dimensional particles trajectories shown in figure 13 for (a) the bulk flow field, (b) large particles and (c) small particles. The colour scheme is the same as that used to identify the trajectories in figure 13. The initial positions of the particles are shown with a colour map and regular grid above the no-net-flow line (dashed white line). The corresponding outgoing positions are shown with the deformed grid and colour map. The steady-state concentration $\phi ^s(0,y,z)$ is shown in the background.

Figure 13(b) and movie 4 show selected three-dimensional large-particle paths starting from their inflow position above the concentration discontinuity near the free surface. Large particles in the centre of the flow are transported forwards to the front, where they are overridden by the bulk flow and appear to move backwards. However, they are able to segregate in the breaking-size-segregation wave and move up into progressively faster moving regions of the flow. This allows them to eventually move forwards and create looping trajectories as the grains are subsequently transported laterally into the levee walls. The mapping between the inflow and outflow positions is shown in figure 14(b). The recirculating loops generate a swirl-like structure in the core of the levee, as shown in the mapping in figure 14(b). Conversely, large particles that start closer to the sides of the central channel end up on the surface of the levee to form a large-rich carapace.

Figure 13(c) and movie 5 show selected three-dimensional small-particle paths that start between the no-net-flow line and the concentration discontinuity. They look qualitatively very similar to those of the bulk flow field, although particle segregation speeds their descent to the base of the flow. The small particles are all recirculated internally to the flow, so that none appear on the free surface. As shown in figure 14(c), small particles that start closer to the sides of the central channel tend to mix with the large grains and form the swirls that are deposited in the levee cores, consistent with the experimental observations of Kokelaar et al. (Reference Kokelaar, Graham, Gray and Vallance2014). These simulations also explicitly show that a pure layer of fine particles forms at the base of the flow, and that the lateral margins of the central channel and the internal levee wall are fines rich. As Kokelaar et al. (Reference Kokelaar, Graham, Gray and Vallance2014) pointed out, this fines lining could be very significant, because in fully coupled simulations (Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021, figure 7) a basal fine-grained layer significantly enhances the flow's speed.

7.1. Varying the particle-size ratio

To change the particle-size ratio $R=d^l/d^s$, the particle sizes need to be changed. As a result, the local volume-fraction-weighted average particle size $\bar {d}$, defined in (5.4), also changes, and this in turn affects the segregation-velocity magnitude (5.9). To minimize this secondary effect, the particle sizes are changed in such a way that the mean inflow particle diameter $\bar {d}_{in} = \bar {\phi }^l_{in} d^l + \bar {\phi }^s_{in} d^s = 0.45$ mm. This is the same value as that used in the simulations of the self-channelized flow in §§ 3 and 4, as well as for the small-particle concentration in § 6. The values of the particle diameters used in this section are summarized in table 4.

Table 4. The particle diameters corresponding to the assumed particle-size ratio in each panel of figure 15. The particle diameters are chosen so that when $\bar {\phi }^s_{in}=0.8$, the average particle diameter at the inflow $\bar {d}_{in} = \bar {\phi }^l_{in} d^l + \bar {\phi }^s_{in} d^s = 0.45$ mm.

Figure 15 shows the concentration distribution on the centreline $y=0$ and on the plane $\xi =0.1$ m for four values of the grain-size ratio $R$. The case $R=1.49$ is the same as in § 6 and serves as a benchmark for comparison. Figure 15(a–d) show that as the grain-size ratio is decreased, the breaking-size-segregation wave becomes progressively longer and extends further upstream. In addition, the swirl of large particles that wraps around the central eye of small particles, develops a more noticeable structure. Note that when $R=1.23$, the back of the breaking-size segregation wave extends out of the plotted domain in figure 15(d), but it is still within the computational domain $\xi \in [-0.4,0.8]$ m. Figure 15(e, f), show that when $R$ is decreased from $1.63$ to $1.49$, the width of the swirl of large grains that gets partially deposited in the levee core widens, and develops more structure. This is also true for $R=1.36$ and $R=1.23$, but it can not solely be deduced from figure 15(g,h), because the sections at $\xi =0.1$ m also cut across the back of the breaking-size segregation wave in the central channel. This widening of the partially deposited breaking-size-segregation waves implies that when $R$ is decreased, there are more fine grains in the levee core.

Figure 15. Colour maps of the small particle concentration $\phi ^s$ at time $\tau = 100$ s with $\bar {\phi }^s_{in} = 0.8$ and for particle-size ratios (a,e) $R=1.63$, (b, f) $R=1.49$, (c,g) $R=1.36$, (d,h) $R=1.23$. The average particle diameter at the inflow $\bar {d}_{in}=0.45$ mm. This constrains the particle diameters, which are summarized in table 4. Panels (a–d) show the $y=0$ centre plane, while panels (e–h) show the cross-section at $\xi = 0.1$ m. The complete temporal evolution for $R=1.63$, $1.36$ and $1.23$ are shown in movies 6–8. The case $R=1.49$ is the same as that in § 6.

Movies 6–8 show the full temporal evolution of the solution for the cases $R=1.63$, $1.36$ and $1.23$, using the same format as movie 2 for $R=1.49$. The surface layer of large particles still reaches the flow front at approximately $\tau =10$ s. However, as $R$ decreases, it takes longer for the large particles that are overrun to segregate upwards and then recirculate towards the front again. The time to reach steady state therefore increases, from approximately $50$ s for the case $R=1.63$ to $90$ s for the case $R=1.23$.

7.2. Varying the inflow concentration

The inflow concentration $\bar {\phi }^s_{in}$ is now varied, whereas the particle sizes and hence the particle-size ratio are held fixed at the same values as § 6. Figure 16 shows three new cases, for $\bar {\phi }^s_{in} = 0.9$, $0.7$ and $0.6$, respectively, as well as the existing case of $\bar {\phi }_s^{in} = 0.8$ as a benchmark for comparison. As $\bar {\phi }^s_{in}$ is decreased, the interface $z=\bar {\phi }^s_{in}h$ between the large and small particles at the inflow reduces in height, and the incoming layer of large particles gets thicker. Since the inflowing particles are inversely graded, they do not segregate and are simply transported by the bulk flow towards the front. The primary effect of reducing the inflow concentration is therefore to increase the thickness of the layer of large particles that develops at the free surface, as shown in figure 16(a,d). These large particles are then either advected and deposited onto the top of the levee to form a carapace, or overrun and recirculated before being deposited internally within the levees.

Figure 16. Colour maps of the small particle concentration $\phi ^s$ at time $\tau = 100$ s for depth-averaged small-particle inflow concentrations (a,e) $\bar {\phi }^s_{in} = 0.9$, (b, f) $\bar {\phi }^s_{in} = 0.8$, (c,g) $\bar {\phi }^s_{in} = 0.7$, (d,h) $\bar {\phi }^s_{in} = 0.6$. The particle diameters are $d^l = 0.61$ mm and $d^s = 0.41$ mm, which corresponds to a particle-size ratio $R=1.49$. Panels (a–d) show the $y=0$ centre plane, whereas panels (e–h) show the cross-section at $\xi = 0.1$ m. The complete temporal evolution for $\bar {\phi }^s_{in} = 0.9$, $0.7$ and $0.6$ are shown in movies 9–11. The case $\bar {\phi }_s^{in} = 0.8$ is the same as that in § 6.

For the case $\bar {\phi }^s_{in}=0.9$, the surface carapace of large particles on the levee walls gets very thin. As a result, the steady-state depth-integrated concentration (movie 9) has alternating bands of large/fine/large/fine material parallel to the downslope direction as one moves inwards from the levee margins towards the centre of the channel. These bands are due to the fact that a significant amount of large material that was in the breaking-size-segregation wave gets trapped close to the inner margin of the levee wall, and the carapace of large grains is thicker at both the levee margin and close to the interface between the levee and the central channel. This effect can also be seen for the original simulation in figures 10 and 11 for $\bar {\phi }^s_{in}=0.8$ with $R=1.49$ and movie 2, but seems less pronounced. It is also evident in movie 7, but is a weak effect in the other movies.

Movies 9–11 show the temporal evolution for the cases $\bar {\phi }^s_{in}=0.9$, $0.7$ and $0.6$. Once the large particles pass over the front at approximately $\tau =10$ s again, a breaking-size-segregation wave forms. As $\bar {\phi }^s_{in}$ decreases, there are more large particles that are overrun at the front, and the breaking wave is consequently richer in large particles, and the low concentration eye of fine grains diminishes in size and appears to move towards the front of the wave. Due to the increased thickness of the free-surface layer of large particles, the breaking size segregation wave gets shifted progressively upslope as $\bar {\phi }^s_{in}$ is reduced. The length of the breaking-size-segregation wave increases as the proportion of large grains is increased, but the change in length is not as strong as the dependence on the size ratio $R$. This reflects the fact that although the segregation-velocity magnitude increases with increasing local mean particle diameter $\bar {d}=\phi ^ld^l+\phi ^sd^s$, the normal component of the large particle segregation velocity (6.9a–c) decreases with increasing large particle concentration.

The increased depth of the incoming layer of large particles at the free surface also increases the thickness of the carapace of large particles that flows out of the domain at the sides of the main channel and in the cross-sections of the levees shown in figure 16(e–h). As the small-particle concentration $\bar {\phi }^s_{in}$ is decreased, the mixed central core in the levee wall diminishes in size and the levees become almost entirely composed of large particles (figure 16h). This imposes a limit on the ability of self-channelized flows to transport the incoming mixture in a steady manner. If there are too many large particles, they cannot be accommodated in the levee walls or in the narrow boundary layers on either side of the no-net-flow line (in the central moving channel) and instead, large particles will continue to accumulate at the flow front, as in the solutions of Gray & Ancey (Reference Gray and Ancey2009) and Gray & Kokelaar (Reference Gray and Kokelaar2010).