2.1. The segregation-advection-diffusion equation and boundary conditions

Consider a bidisperse mixture of large and small particles flowing down a chute inclined at an angle $\zeta$ to the horizontal. Define a Cartesian coordinate system $Oxyz$ with the $x$ -axis pointing downslope, the $y$ -axis pointing in the cross-slope direction and the $z$ -axis being the upward pointing normal. Unit vectors $\boldsymbol i$ , $\boldsymbol j$ and $\boldsymbol k$ are defined in each of these directions. All the particles have the same intrinsic density $\rho _*$ and the solids volume fraction $\varPhi$ is assumed to be constant. The volume fraction of large and small particles per unit granular volume are defined as $\phi ^l$ and $\phi ^s$ , respectively. It follows that their sum is equal to unity

(2.1) \begin{align} \phi ^s+\phi ^l=1. \end{align}

The particle concentrations $\phi ^l$ and $\phi ^s$ satisfy the segregation-advection-diffusion equations (Gray & Thornton Reference Gray and Thornton2005; Gray & Chugunov Reference Gray and Chugunov2006; Gray & Ancey Reference Gray and Ancey2011; Gray Reference Gray2018)

(2.2) \begin{align} \frac {\partial \phi ^l}{\partial t}+{ \boldsymbol{\nabla }}\,{\boldsymbol{\cdot }}\,\big(\phi ^l\boldsymbol{u} \big)+{ \boldsymbol{\nabla }}\,{ \boldsymbol{\cdot }}\, \big(f_{ls}\phi ^l\phi ^s\boldsymbol{e} \big) &= {\boldsymbol {\nabla }}\,{\boldsymbol{\cdot }}\,\big(D{ \boldsymbol{\nabla }}\phi ^l \big), \\[-12pt] \nonumber \end{align}

(2.3) \begin{align} \frac {\partial \phi ^s}{\partial t}+{\boldsymbol{\nabla }}\,{ \boldsymbol{\cdot }}\,\big(\phi ^s\boldsymbol{u} \big)-{ \boldsymbol{\nabla }}\,{ \boldsymbol{\cdot }}\,\big(f_{ls}\phi ^s\phi ^l\boldsymbol{e} \big) &= {\boldsymbol {\nabla }}\,{ \boldsymbol{\cdot }}\,\big(D{ \boldsymbol{\nabla }}\phi ^s \big), \\[10pt] \nonumber \end{align}

where $\boldsymbol {\nabla}$ is the gradient operator, $\boldsymbol {\cdot}$ is the dot product, ${\boldsymbol u}=(u,v,w)=u{\boldsymbol i}+v\kern-1pt{\boldsymbol j}+w{\boldsymbol k}$ is the bulk velocity field, $f_{ls}$ is the segregation velocity magnitude, ${\boldsymbol e}=(e_x,e_y,e_z)$ is the unit vector in the direction in which large particles segregate and $D$ is the isotropic particle diffusion coefficient. The segregation direction will be discussed at greater length in § 4, but it is usually aligned in the opposite direction to the gravitational acceleration vector $\boldsymbol g$ . When (2.2)–(2.3) are summed over all constituents, the summation condition (2.1) implies that the bulk velocity field $\boldsymbol u$ is incompressible

(2.4) \begin{align} {\boldsymbol{\nabla }}\,{ \boldsymbol{\cdot }}\,{\boldsymbol u}=0. \end{align}

The system (2.2)–(2.3) can be reduced to a single equation for the small-particle concentration $\phi ^s$ (say), by using the summation condition (2.1) to eliminate $\phi ^l$ , to give

(2.5) \begin{align} \frac {\partial \phi ^s}{\partial t}+{ \boldsymbol{\nabla }}\,{ \boldsymbol{\cdot }}\,(\phi ^s{\boldsymbol u})-{\boldsymbol {\nabla }}\,{\boldsymbol {\cdot }}\,\big(f_{ls}\phi ^s(1-\phi ^s){\boldsymbol e} \big) = {\boldsymbol {\nabla }}\,{\boldsymbol {\cdot }}\,(D{\boldsymbol{\nabla }}\phi ^s). \end{align}

The small-particle concentration therefore evolves in time due to advection by the bulk flow, downward segregation parallel to the direction $\boldsymbol e$ and isotropic diffusion. This is subject to a no-flux condition on the boundaries of the domain. Assuming that the bulk flow satisfies kinematic conditions on the free surface $z=s(x,y,t)$ and base $z=b(x,y,t)$ of the flow (see e.g. Savage & Hutter Reference Savage and Hutter1989; Gray et al. Reference Gray, Tai and Noelle2003), i.e. that

(2.6) \begin{align} \frac {\partial s}{\partial t}+u \frac {\partial s}{\partial x}+v \frac {\partial s}{\partial y}-w&=0,\quad \hbox{on}\quad z=s(x,y,t) , \\[-12pt] \nonumber \end{align}

(2.7) \begin{align} \frac {\partial b}{\partial t}+u \frac {\partial b}{\partial x}+v \frac {\partial b}{\partial y}-w&=0,\quad \hbox{on}\quad z=b(x,y,t) , \\[10pt] \nonumber \end{align}

the no-flux condition for the small particles reduces to

(2.8) \begin{align} f_{ls}\phi ^s(1-\phi ^s){\boldsymbol e}\,{ \boldsymbol{\cdot }}\,{\boldsymbol n}+D{\boldsymbol {\nabla }}\phi ^s\,{\boldsymbol {\cdot }}\,{\boldsymbol n}=0,\quad \text{on}\quad z=s(x,y,t),\,b(x,y,t) , \end{align}

where ${\boldsymbol n}=(-\partial z/\partial x,-\partial z/\partial y,1)$ is the normal to the free surface or base of the flow (see e.g. Gray & Chugunov Reference Gray and Chugunov2006; Gray & Ancey Reference Gray and Ancey2011).

2.2. Downslope uniformity

The straight chute experiments of Pearse et al. (Reference Pearse, Johnson and Gray2026) suggest that the stripes are only weakly dependent on the downslope coordinate $x$ . This motivates the simplification of the problem by assuming downslope uniformity of the particle concentration $\phi ^s(y,z,t)$ . Note that this assumption will be made for the theory developed here and the single-cell computations made in §§ 3 and 5, but is not made when developing the general depth-averaged theory in § 4 or when modelling the spatial development of the stripes along the chute in § 7. Throughout this paper the velocity is assumed to be independent of the downslope coordinate and independent of time, i.e. ${\boldsymbol u}={\boldsymbol u}(y,z)$ as shown in figure 2. As a result the incompressibility condition (2.4) reduces to

(2.9) \begin{align} \frac {\partial v}{\partial y}+\frac {\partial w}{\partial z}=0, \end{align}

and the segregation-advection-diffusion equation (2.5) becomes

(2.10) \begin{align} \frac {\partial \phi ^s}{\partial t}+\frac {\partial }{\partial y}(\phi ^s v)+\frac {\partial }{\partial z}(\phi ^s w)-\frac {\partial }{\partial z}\left (q\phi ^s(1-\phi ^s) \right ) = \frac {\partial }{\partial y}\left (D\frac {\partial \phi ^s}{\partial y} \right ) +\frac {\partial }{\partial z}\left (D\frac {\partial \phi ^s}{\partial z}\right ), \end{align}

where the normal component of the small-particle segregation velocity is

(2.11) \begin{align} q=f_{ls}\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{k}. \end{align}

In general, the segregation velocity $q$ is linearly dependent upon the shear rate and the average particle diameter, and has further dependencies on the particle-size ratio, intrinsic particle density, pressure and the particle concentrations, while the diffusion coefficient $D$ is linearly dependent on the shear rate and has a quadratic dependence on the average particle diameter (Savage & Lun Reference Savage and Lun1988; Utter & Behringer Reference Utter and Behringer2004; Gray & Thornton Reference Gray and Thornton2005; Trewhela et al. Reference Trewhela, Ancey and Gray2021; Maguire et al. Reference Maguire, Barker, Rauter, Johnson and Gray2024). Since gradients of the secondary flow are much smaller than the downslope velocity gradient through the depth of the flow, the shear rate $\dot \gamma =2\lVert {\unicode{x1D63F}}\rVert$ is approximately equal to $|\partial u/\partial z|$ . The downslope velocity profile therefore plays a dominant role in setting both the segregation rate and the diffusion, but does not appear explicitly in the equations.

2.3. Assumed form of the secondary vortices

The root cause of the secondary vortices is unclear and still a subject of active research. As discussed in the introduction, likely explanations are second-normal-stress differences (McElwaine et al. Reference McElwaine, Takagi and Huppert2012; Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2021; Kim & Kamrin Reference Kim and Kamrin2023; Levay et al. Reference Levay, Claudin, Somfai and Börzsönyi2025; Gadal et al. Reference Gadal, Johnson and Gray2026) and/or the formation of a dilute Leidenfrost layer at the base of the avalanche which leads to instability and basal slip (Campbell Reference Campbell1989; Forterre & Pouliquen Reference Forterre and Pouliquen2001; Eshuis et al. Reference Eshuis, van der Weele, van der Meer and Lohse2005; Börzsönyi et al. Reference Börzsönyi, Ecke and McElwaine2009; Wang et al. Reference Wang, Jing, Kwok, Sobral, Weinhart and Thornton2025). Following Pearse et al. (Reference Pearse, Johnson and Gray2026), the secondary vortices are parameterised in a simple way in order to solve for the associated segregation problem. The vortices are assumed to be rectangular in shape with incompressible bulk velocity field with transverse components

(2.12) \begin{align} v &=\varOmega \sin {\left (\frac {\pi y}{W}\right )}(2z-H),\quad 0 \leqslant z \leqslant H, \quad mW\leqslant y\leqslant nW, \\[-12pt] \nonumber \end{align}

(2.13) \begin{align} w &=-\frac {\varOmega \pi }{W}\cos {\left (\frac {\pi y}{W}\right )}z(z-H),\quad 0 \leqslant z \leqslant H,\quad mW\leqslant y\leqslant nW, \\[10pt] \nonumber \end{align}

where $\varOmega$ determines the rotation rate, $H$ is the constant vortex height and $W$ the constant vortex width, as shown in figure 2. A typical flow has $n{-}m$ counter-rotating vortices, where $m$ and $n\gt m$ are integers. Each stripe measured from centre to centre of adjacent large–rich bands consists of two counter-rotating vortices. Note that this velocity field is constructed to be incompressible and satisfy the no normal flux condition ${\boldsymbol u}\,{ \boldsymbol{\cdot }}\,{\boldsymbol n}=0$ on the cell boundaries. The cross-slope velocity (2.12) is non-zero at the base, i.e. basal slip is assumed. This is consistent with the experiments in § 6, which show that there is a significant amount of basal slip in the downslope direction. The slip could be facilitated by the formation of a low-density Leidenfrost layer at the base of the flow.

2.4. Non-dimensionalisation

In order to investigate the competing effects of segregation, secondary flow and diffusion, this paper assumes that the segregation velocity $q$ and the diffusion coefficient $D$ are constant. Introducing the non-dimensional variables, denoted by a tilde, via the scalings

(2.14) \begin{align} (y,z)=H(\tilde {y},\tilde {z}),\quad (v,w)=\varOmega H(\tilde {v},\tilde {w}),\quad t=(H/q)\tilde {t}, \end{align}

where $\varOmega H$ is a typical rotational velocity magnitude and $H/q$ is a typical time scale for segregation to occur. These scalings imply the incompressibility condition (2.9) becomes

(2.15) \begin{align} \frac {\partial \tilde {v}}{\partial \tilde {y}}+\frac {\partial \tilde {w}}{\partial \tilde {z}}=0, \end{align}

and the non-dimensional segregation-advection-diffusion equation (2.10) becomes

(2.16) \begin{align} \frac {\partial \phi ^s}{\partial \tilde {t}}+\frac {1}{\varLambda }\frac {\partial }{\partial \tilde {y}}(\phi ^s \tilde {v})+\frac {1}{\varLambda }\frac {\partial }{\partial \tilde {z}}(\phi ^s \tilde {w})-\frac {\partial }{\partial \tilde {z}}\big ( \phi ^s(1-\phi ^s) \big ) = \frac {1}{\textit{Pe}}\left ( \frac {\partial ^2\phi ^s}{\partial \tilde {y}^2} +\frac {\partial ^2 \phi ^s}{\partial \tilde {z}^2} \right ), \end{align}

where the non-dimensional parameters

(2.17a,b) \begin{align} \varLambda =\frac {q}{\varOmega H},\quad \hbox{and}\quad \textit{Pe}=\frac {qH}{D}, \end{align}

are the ratio of a typical segregation velocity to a typical vortex velocity and the ratio of a typical segregation velocity to a typical diffusive velocity, respectively. The scalings also imply that the no-flux boundary condition in (2.8) along the surface and base of the flow, when ${\boldsymbol n}=\pm {\boldsymbol k}$ , becomes

(2.18) \begin{align} \phi ^s(1-\phi ^s)+\frac {1}{\textit{Pe}}\frac {\partial \phi ^s}{\partial \tilde {z}}=0, \quad \tilde {z}=0,1, \end{align}

while on the sides of the flow (when ${\boldsymbol n}=\pm {\boldsymbol j}$ ) it reduces to the zero gradient condition

(2.19) \begin{align} \frac {\partial \phi ^s}{\partial \tilde {y}}=0, \quad \tilde {y}=\tilde {W}m,\tilde {W}n, \end{align}

provided that the segregation direction $\boldsymbol e$ does not have a component in the cross-slope direction, i.e. ${\boldsymbol e}\,{ \boldsymbol{\cdot }}{\boldsymbol j}=0$ . In addition an initial condition is required to solve for the small-particle concentration

(2.20) \begin{align} \phi ^s(\tilde {y},\tilde {z},0)=\phi ^s_0(\tilde {y},\tilde {z}), \quad (\tilde {y},\tilde {z}) \in [\tilde {W}m,\tilde {W}n]\times [0,1]. \end{align}

The non-dimensional scalings (2.14) imply that the velocity field (2.12)–(2.13) becomes

(2.21) \begin{align} \tilde {v} &=\sin {\left (\frac {\pi \tilde {y}}{\tilde {W}}\right )}(2\tilde {z}-1),\quad 0 \leqslant \tilde {z} \leqslant 1,\quad \tilde {W}m \leqslant \tilde {y} \leqslant \tilde {W}n, \\[-12pt] \nonumber \end{align}

(2.22) \begin{align} \tilde {w} &=-\frac {\pi }{\tilde {W}}\cos {\left (\frac {\pi \tilde {y}}{\tilde {W}}\right )}\tilde {z}(\tilde {z}-1),\quad 0 \leqslant \tilde {z} \leqslant 1,\quad \tilde {W}m \leqslant \tilde {y} \leqslant \tilde {W}n, \\[10pt] \nonumber \end{align}

where $\tilde {W}$ is the non-dimensional width of the cell. In this paper it is assumed that $\tilde {W}=1.85$ , so that the velocity field is the same as that used by Pearse et al. (Reference Pearse, Johnson and Gray2026) in the absence of diffusion, and the vortex aspect ratio matches their straight chute experiments. Figure 7 of Pearse et al. (Reference Pearse, Johnson and Gray2026) shows a velocity vector plot, contour plots of the velocity components and the streamfunction for the velocity field above.

3.1. The numerical method

Numerical solutions of the segregation-advection-diffusion equation (2.16) are obtained using the high-resolution central scheme of Kurganov & Tadmor (Reference Kurganov and Tadmor2000) with the Runge–Kutta–Chebyshev adaptive time-stepping method of Medovikov (Reference Medovikov1998). This time stepper is an efficient way of accounting for the stiffness of the system due to the presence of diffusion. The velocity field (2.21)–(2.22) is prescribed in a single vortical cell $(\tilde {y},\tilde {z})\in [0,1.85]\times [0,1]$ . All simulations begin from an inversely graded initial state

(3.1) \begin{align} \phi ^s_0= \begin{cases} 0, \quad \tilde {z}\gt 0.7,\\[4pt] 1, \quad \tilde {z}\leqslant 0.7, \end{cases} \end{align}

to reflect the 30 : 70 mix of large (green) and small (white) particles used in the experiments, and the no-flux boundary conditions (2.18) and (2.19) are applied at the edges of the domain.

3.2. Time-dependent solutions with diffusion

The evolution of a bidisperse mixture from an inversely graded initial distribution is shown in figure 3 for a single clockwise-rotating vortex cell. The flow has parameter values $(\textit{Pe},\varLambda )=(50,1)$ , which implies that the diffusion is relatively weak compared with the segregation, but that the vortical motion is comparable in magnitude. After $\tilde {t}=0$ , the sharp horizontal interface separating the large particles above, from the small particles below, is diffused and rotated by the velocity field (2.21)–(2.22). As in the case with no diffusion (see figure 8 of Pearse et al. Reference Pearse, Johnson and Gray2026), the small particles concentrate in the upwelling part of the flow (on the left of the cell) and the large particles accumulate in the downwelling region (on the right of the cell), as shown in figure 3(b). The diffuse interface steepens and breaks at approximately $\tilde {t}=1.4$ non-dimensional units, as shown in movie 2 in the online supplementary material, to form a diffuse breaking-size-segregation wave (figure 3 c). This grows in height and width until it spans the full depth of the flow, and then oscillates backwards and forwards until it reaches steady state by approximately $\tilde {t}=20$ non-dimensional time units. Although many features of the evolution towards steady state are similar to those in the absence of diffusion, the diffusion results in the interface between large and small particles being much less sharp, which may reflect real flows more accurately.

Figure 3.

The evolution of the small-particle concentration $\phi ^s$ at times (a) $\tilde {t}=0$ , (b) $\tilde {t}=1$ , (c) $\tilde {t}=2$ , (d) $\tilde {t}=3$ , (e) $\tilde {t}=4$ and steady state (f) $\tilde {t}=20$ for a 30 : 70 mix of large and small particles with parameter values $(\textit{Pe},\varLambda )=(50,1)$ . The solution is shown in a single clockwise-rotating secondary vortex cell with velocity field (2.21)–(2.22) applied. Movie 2 shows the complete time-dependent solution. The contour lines are shown in the range $\phi ^s\in [0.1,0.9]$ in 0.08 increments.

Figure 3. Long description

The time evolution of a 30:70 mixture of large-green material and small-white material in a single clockwise-rotating vortical cell. Initially, large-green material sits above small-white material with a horizontal interface separating the two phases. By one non-dimensional time unit, the interface is diffused and has rotated to form an s-shape with small material concentrated in the upwelling region (left side of cell) and large material concentrated in the downwelling region (right side of cell). By two non-dimensional time units, the interface has steepened to vertical and breaks to form a breaking-size-segregation wave separarying the pure phases of small and large material. By three non-dimensional time units, the breaking wave has grown to be the full depth of the flow after which it oscillates backwards and forwards until it reaches steady state.

At steady state the diffuse breaking-size-segregation wave separates an almost pure upwelling region of small particles (on the left) from an almost pure downwelling region of large particles (on the right). The development of the breaking wave is crucial for the formation of longitudinal stripes that extend all the way through the flow depth, because it allows large particles that are advected under small particles to segregate back up towards the surface, and conversely allows small particles that are advected over large grains to percolate down towards the base. Pearse et al. (Reference Pearse, Johnson and Gray2026) showed that in the absence of diffusion, at steady state, the large- and small-particle paths form two interpenetrating co-rotating sub-vortices in the projected cross-slope plane, that allow stable stripes to be maintained. The full three-dimensional particle paths (taking into account the downslope motion) develop very complicated trajectories that depend strongly on the mixing ratio of the two species.

Figure 4.

Perspective view of the time-dependent formation of two diffuse longitudinal stripes of large and small particles for $(\textit{Pe},\varLambda )=(50,1)$ and a 30 : 70 mix. The solution spans four counter-rotating vortex cells. The small-particle concentration is shown at non-dimensional times (a) $\tilde {t}=0$ , (b) $\tilde {t}=1$ , (c) $\tilde {t}=2$ , (d) $\tilde {t}=3$ , (e) $\tilde {t}=4$ and (f) $\tilde {t}=20$ when the solution is close to steady state. The colour bar is the same as that in figure 3. Movie 3 shows the full time-dependent evolution. The contour lines are shown in the range $\phi ^s\in [0.1,0.9]$ in 0.08 increments.

Figure 4. Long description

The three-dimensional time evolution of a 30:70 mixture of large-green material and small-white material across four counter-rotating vortical cells. The initially inversely graded flow produces a green surface which persists at one non-dimensional time unit even though the counter-rotating vortices have rotated the internal interface between large and small particles to form a diffuse wavy interface. Between 2 to 4 non-dimensional time units, small particles appear at the surface to form compositional stripes. The small white bands are centred on the upwelling part of the flow, whereas the green large particle bands are concentrated in the downwelling regions. The wave then oscillates backwards and forwards until a diffuse steady-state is approached at twenty non-dimensional time units.

To aid understanding, a three-dimensional perspective view of the time-dependent computations in figure 3 is plotted in figure 4 for four counter-rotating vortex cells by using the symmetry of the solution. This shows how the evolving particle-size distribution in the $(\tilde {y},\tilde {z})$ -plane manifests itself at the free surface, which is the only thing that is visible in experiments. The initially inversely graded flow at $\tilde {t}=0$ produces a green surface which persists at $\tilde {t}=1$ even though the counter-rotating vortices have rotated the internal interface between large and small particles to form a diffuse wavy interface (figure 4 b). It is only between $\tilde {t}=2$ and 4 non-dimensional time units that small particles appear at the surface to form compositional stripes. The small white stripe is centred on the upwelling part of the flow, whereas the green large-particle bands are concentrated in the downwelling regions. The wave then oscillates backwards and forwards until a diffuse steady state is approached at $\tilde {t}=20$ , as can be seen in movie 3 in the online supplementary material. With these parameters the approximate 30 : 70 mixing ratio is reflected in the width of the steady-state large- and small-particle bands.

Figure 5.

Steady-state solutions for the small-particle concentration $\phi ^s$ for non-dimensional parameters $\textit{Pe}=\infty ,500,50,20,5,1$ and $\varLambda =\infty ,100,10,3/2,1,0.5$ for a 30 : 70 mix of large and small particles. Increasing $ \textit{Pe}$ corresponds to reducing the amount of diffusion compared with segregation, whereas increasing $\varLambda$ corresponds to reducing the strength of the secondary flow relative to segregation. Each panel consists of a single clockwise-rotating vortex cell $[0,1.85]\times[0,1]$ with velocity field (2.21)–(2.22).

Figure 5. Long description

A total of 36 steady-state solutions in a single clockwise-rotating vortical cell spanning the full (Pe,Lambda) parameter space combining 6 values of Pe (1-infinity) and 6 values of Lambda (0.5-infinity). Lambda is on the x-axis and Pe is on the y-axis. The top row (Pe=infinity) corresponds to the case of no diffusion with sharp bands of large and small particles separated by a breaking-size-segregation wave whose width decreases as Lambda tends to infinity. When there’s no circulation (far right column) inversely-graded distributions are produced. The second row has relatively weak diffusion. Strong vortical motion (left side) produces similar solutions to the row above. As the vortical motion diminishes (move across to the right) wavy interfaces are produced. As diffusion increases (move further down y axis), the solutions in the second row become gradually more diffuse. When diffusion is too strong (bottom row), near-uniform distributions are produced.

3.3. The steady-state parameter space

Steady-state solutions have been calculated for a wide range of values of the non-dimensional parameters $ \textit{Pe}$ and $\varLambda$ , and are shown in figure 5 for a 30 : 70 mix of large and small particles. The top row corresponds to the case investigated by Pearse et al. (Reference Pearse, Johnson and Gray2026), when there is no diffusion and $\textit{Pe}=\infty$ . Provided that there is some vortical motion (i.e. $\varLambda \ne \infty$ ) the flow has a pure band of large particles in the downwelling region and a pure band of small particles that is concentrated in the upwelling section, which are separated by a breaking-size-segregation wave that spans the full depth of the flow (Thornton & Gray Reference Thornton and Gray2008; Gray & Ancey Reference Gray and Ancey2009; Edwards et al. Reference Edwards, Rocha, Kokelaar, Johnson and Gray2023; Pearse et al. Reference Pearse, Johnson and Gray2026). As the strength of the vortical motion decreases, which is equivalent to increasing $\varLambda$ , the width of the breaking-size-segregation wave diminishes, so that as $\varLambda \rightarrow \infty$ it tends to a vertical shock separating large- and small-particle bands. The case when $\varLambda =\infty$ corresponds to no vortical motion, and is special. In this case, there is no rotational motion to steepen the interface between large and small particles, so it remains horizontal and separates the large particles above from the small particles below.

The right-hand column of figure 5 shows the effect of increasing the amount of diffusion in the absence of vortical motion, which is equivalent to decreasing $ \textit{Pe}$ and holding $\varLambda =\infty$ . These solutions are all inversely graded (i.e. large particles are concentrated near the surface and small particles are concentrated near the base), but the sharp shock that develops in the absence of diffusion is progressively diffused. These particle-size distributions are equivalent to the exact solution derived by Gray & Chugunov (Reference Gray and Chugunov2006) in their § 3.

The solutions in the second row of figure 5 have relatively weak diffusion. When the vortical motion is strong, i.e. $\varLambda$ is low, the solutions look very similar to those in the absence of diffusion in the top row ( $\textit{Pe}=\infty$ ). However, as the vortical motion diminishes, the lateral diffusion can compete against the overturning effect of secondary flow, and a steady state develops with a relatively sharp ‘S’-shaped interface between the large and small particles (see the case $\textit{Pe}=500$ , $\varLambda =100$ ). The ability of the lateral diffusion to balance the overturning effect of the secondary flow becomes progressively stronger as the strength of the diffusion is increased (i.e. $ \textit{Pe}$ is decreased). Indeed for the case $\textit{Pe}=20$ and $\varLambda =10$ , the diffusion is strong enough to prevent small particles reaching the surface of the flow, so in an experiment one would not see the formation of stripes even though there was an underlying vortical flow. Note that, as the diffusion is increased, the ‘S’-shaped interface between large and small particles becomes progressively more diffuse. Indeed, if the diffusion is too strong (e.g. for $\textit{Pe}=1$ in the bottom row) then it produces an almost uniform particle-size distribution throughout the cell.

5.1. Basic assumptions

The theory is now reduced to the case in which there is no basal topography, $b=0$ , and the free surface lies at a constant height $h=s=H$ . Furthermore, the stripes and the velocity field are assumed to be uniform in the downslope direction, so that the $x$ -derivatives are zero, and the diffusion coefficient $D$ is assumed to be constant. In this case (4.4) reduces to

(5.1) \begin{align} \frac {\partial }{\partial t}(H\overline {\phi ^s}) + \frac {\partial }{\partial y}(H\overline {\phi ^s v}) = \frac {\partial }{\partial y}\left (D \frac {\partial }{\partial y}(H\overline {\phi ^s})\right )\!. \end{align}

Following Gray & Kokelaar (Reference Gray and Kokelaar2010a ,Reference Gray and Kokelaar b ), and based on the simulations in § 3 that show there is a region of parameter space where the interface between large and small particles is relatively sharp, the small-particle concentration may reasonably be approximated by

(5.2) \begin{align} \phi ^s= \begin{cases} 0, \hspace {0.3cm} \eta \lt z\leqslant H,\\ 1, \hspace {0.3cm} 0\leqslant z\leqslant \eta ,\\ \end{cases} \end{align}

where $\eta$ is the inversely graded interface height between large and small particles, as sketched in figure 6. It follows that the depth-averaged concentration, defined in (4.3a), satisfies

(5.3) \begin{align} H\overline {\phi ^s}=\eta . \end{align}

An approximation for $\overline {\phi ^s v}$ can be found by substituting the small-particle concentration (5.2) and the cross-slope velocity (2.12) into the integral (4.3c ) to give

(5.4) \begin{align} H\overline {\phi ^s v}=\varOmega \sin \left (\frac {\pi y}{W}\right )(\eta ^2-H\eta ). \end{align}

Substituting the integrals (5.3) and (5.4) into (5.1) implies that the depth-averaged segregation equation for the interface height $\eta$ is

(5.5) \begin{align} \frac {\partial \eta }{\partial t} + \frac {\partial }{\partial y}\left (\varOmega \sin \left (\frac {\pi y}{W}\right )(\eta ^2-H\eta )\right ) = \frac {\partial }{\partial y}\left (D \frac {\partial \eta }{\partial y}\right ). \end{align}

The quadratic structure of the small-particle flux on the left-hand side, implies that this equation will seek to segregate the large and small particles in the $y$ direction, but as opposed to earlier theories, there is a diffusive flux on the right-hand side that can prevent the interface from steepening up to the vertical (Gray & Kokelaar Reference Gray and Kokelaar2010a ,Reference Gray and Kokelaar b ; Woodhouse et al. Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012; Baker et al. Reference Baker, Johnson and Gray2016).

5.2. Non-dimensionalisation

The depth-averaged model is non-dimensionalised using the scalings in (2.14) and

(5.6) \begin{align} \eta =H\tilde {\eta }, \end{align}

to give the non-dimensional depth-averaged segregation-diffusion equation

(5.7) \begin{align} \frac {\partial \tilde {\eta }}{\partial \tilde {t}}+\frac {1}{\varLambda }\frac {\partial }{\partial \tilde {y}}\left ( \sin \left ( \frac {\pi \tilde {y}}{\tilde {W}}\right )(\tilde {\eta }^2-\tilde {\eta }) \right )=\frac {1}{\textit{Pe}}\left (\frac {\partial ^2 \tilde {\eta }}{\partial \tilde {y}^2}\right )\!, \end{align}

where $ \textit{Pe}$ and $\varLambda$ are defined in (2.17). This is the depth-averaged equivalent to the two-dimensional segregation-diffusion equation (2.16). It is subject to a zero gradient boundary condition at the cell interfaces

(5.8) \begin{align} \frac {\partial \tilde {\eta }}{\partial \tilde {y}}=0\quad \hbox{at}\qquad \tilde {y}=0,\tilde {W}, \end{align}

and an initial condition for the interface position

(5.9) \begin{align} \tilde {\eta }(\tilde {y},0)=\tilde {\eta }_0(\tilde {y}), \quad \tilde {y}\in [0,\tilde {W}], \end{align}

where $\tilde {\eta }_0$ is the initial interface height.

Figure 7.

Comparison of the one-dimensional interface height $\tilde {\eta }$ (magenta line) between large and small particles, and the full two-dimensional solution for the concentration $\phi ^s$ in a single clockwise-rotating vortex cell, with the velocity field (2.21)–(2.22), and for a 30 : 70 mix. Panels (a)–(e) show the time-dependent evolution for $(\textit{Pe},\varLambda )=(500,100)$ , while panels (f)–(j) show the solution for $(\textit{Pe},\varLambda )=(500,3/2)$ . The steady-state exact solutions are also shown for (e) $\varXi =2.94$ and (j) $\varXi =196.29$ (thicker semi-transparent light blue line). Movies 4 and 5 show the complete time-dependent evolution. The contour lines are shown in the range $\phi ^s\in [0.1,0.9]$ in 0.08 increments.

Figure 7. Long description

(a) Inversely-graded distribution with large-green material above small-white material with a horizontal interface. (b-e) The interface becomes increasingly wavy until a steady state is reached with large-green particles concentrated in the down-welling region (right) and small white grains concentrated in the upwelling region (left). The depth-averaged and exact solution matches the interface excellently. (f) The initially inversely-graded particle-size distribution becomes wavy in (g). (h) The interface steepens and breaks, but the breaking wave isn't yet the full height of the cell so there is still large-green material near the free-surface. (i) The breaking-wave is now the full height of the cell and separates sharp bands of small-white material (left) from large-green material (right). (j) Sharp bands of small and large material through the flow depth seperated by a breaking-size-segregation wave. The depth-averaged solution replaces the breaking wave with a shock.

5.3. Depth-averaged numerical solutions for the interface

Depth-averaged numerical solutions for the interface are obtained using the same shock-capturing numerical method as § 3.1. Computations are performed for a single depth-averaged vortex cell $\tilde {y}\in [0,1.85]$ with no-flux boundary conditions (5.8) applied at each end of the domain. For consistency, a 30 : 70 mix of large and small particles is used, that is assumed to have a uniform inversely graded distribution across the cell, i.e. the initial interface height

(5.10) \begin{align} \tilde {\eta }_0=0.7,\quad \tilde {y}\in [0,1.85]. \end{align}

To investigate the depth-averaged model’s ability to capture the sharp interface, the case when $(\textit{Pe},\varLambda )=(500,100)$ is studied. Figure 7(a–e) shows the time-dependent evolution of the one-dimensional interface $\tilde {\eta }$ and the corresponding two-dimensional particle-size distribution calculated using the full two-dimensional model developed in § 3. Initially the interface between the large and small particles starts off horizontal, but it is gradually rotated by the secondary flow, so that large particles concentrate in the downwelling region on the right, and the small particles concentrate in the upwelling region on the left. The diffusion is weak compared with the segregation, so the interface in the two-dimensional simulation stays sharp, and is well approximated by the one-dimensional depth-averaged theory throughout the entire time-dependent development, which can be seen in movie 4. The lateral diffusion is, however, strong enough to prevent the interface steepening up to the vertical and forming a breaking-size-segregation wave. The depth-averaged theory, and the full two-dimensional theory, therefore slowly approach a steady state with an ‘S’-shaped interface by $\tilde {t}=500$ non-dimensional time units. The exact steady-state solution is calculated in § 5.4 and is in excellent agreement with both computations. At steady state there is a thin layer of large particles on the surface of the upwelling region, but they are mainly concentrated in the downwelling region. In contrast, at steady state there is still quite a thick layer of small particles at the base of the downwelling region.

Figure 7(f–j) shows a similar simulation for $(Pe,\varLambda )=(500,3/2)$ . The relative strength of the diffusion to segregation is the same, but the secondary recirculation is significantly stronger. In this case the initial motion of the interface is also excellently captured by the depth-averaged theory, but at approximately $\tilde {t}=2$ the interface breaks to form a diffuse breaking-size-segregation wave (Thornton & Gray Reference Thornton and Gray2008; Gray & Ancey Reference Gray and Ancey2009; Gray & Kokelaar Reference Gray and Kokelaar2010a ,Reference Gray and Kokelaar b ; Pearse et al. Reference Pearse, Johnson and Gray2026), which is not captured in the one-dimensional model. The two-dimensional simulations show that the breaking wave grows until it reaches the full height of the cell, and then it slowly oscillates backwards and forwards before settling down to a steady state. Despite the formation of the breaking wave, the one-dimensional model accurately captures the position of the sharp inversely graded interface, away from the breaking wave, during the subsequent evolution. In this case, the large-particle layer at the surface of the flow is so thin in the upwelling region that it can barely be seen, and the large particles concentrate throughout the entire depth of the downwelling region. Although the depth-averaged theory cannot capture the structure of the breaking-size-segregation wave, it does accurately capture its approximate position. The full time-dependent evolution of both models are shown in movies 4 and 5.

5.4. Steady-state solution for the ‘S’-shaped interface

At steady state, an exact one-dimensional depth-averaged solution can be derived for the interface $\tilde {\eta }$ . In this case the depth-averaged segregation-diffusion equation becomes

(5.11) \begin{align} \frac {1}{\varLambda }\frac {{d}}{\mathrm{d} \tilde {y}} \left ( \sin {\left (\frac {\pi \tilde {y}}{\tilde {W}}\right )}(\tilde {\eta }^2-\tilde {\eta }) \right )=\frac {1}{\textit{Pe}}\left ( \frac {{d}^2\tilde {\eta }}{\mathrm{d} \tilde {y}^2} \right ), \end{align}

and since $ \textit{Pe}$ and $\varLambda$ are constant, this can be integrated once with respect to $\tilde {y}$ . Applying the boundary condition (5.8) at $\tilde {y}=0$ (say), yields the ordinary differential equation (ODE)

(5.12) \begin{align} \frac {\mathrm{d}\tilde {\eta }}{\mathrm{d}\tilde {y}}=\frac {\textit{Pe}}{\varLambda }\sin {\left (\frac {\pi \tilde {y}}{\tilde {W}}\right )}(\tilde {\eta }^2-\tilde {\eta }). \end{align}

This is separable, and can be integrated to give the exact solution

(5.13) \begin{align} \tilde {\eta }=\displaystyle \frac {1}{1+A \exp \left (-\varXi \cos \left (\frac {\pi \tilde {y}}{\tilde {W}}\right )\right )} , \end{align}

where $A$ is a constant of integration and the non-dimensional parameter

(5.14) \begin{align} \varXi =\frac {\textit{Pe} \tilde {W}}{\pi \varLambda }=\frac {HW\varOmega }{\pi D} \end{align}

sets the non-dimensional length scale of the $S$ -shaped interface. It can be thought of as a ratio of a typical circulation velocity to the velocity induced by lateral diffusion, and is independent of the segregation rate $q$ . Moreover, since $ \textit{Pe}$ and $\varLambda$ appear as a ratio, it follows that the steady-state theory has the same interface shape for the same value of $\textit{Pe}/\varLambda$ . For instance, the shape of the interface for ( $\textit{Pe},\varLambda )=(500,100)$ is exactly the same as for ( $\textit{Pe},\varLambda )=(50,10)$ , as illustrated in figure 5, although in the latter case the particle-size distribution is more diffuse in the full model, due to the lower value of $ \textit{Pe}$ . Note that low values of $\varXi$ generate smooth wavy inversely graded interfaces that do not approach the free surface, whereas high values of $\varXi$ correspond to solutions in which stripes will be visible at the free surface, as shown in figure 7(e, j), respectively.

The constant of integration $A$ in (5.13) is determined iteratively, by enforcing that the total volume of small grains in the secondary vortex is conserved, i.e. that

(5.15) \begin{align} \int _0^{\tilde {W}} \tilde {\eta }\,\mathrm{d}\tilde {y} = \tilde {W}\tilde {\eta }_0. \end{align}

As already shown in § 5.3, the exact solution (5.13) accurately predicts the final steady-state interface position provided that the secondary recirculation and the diffusion are not too strong compared with the segregation rate (e.g. for $(\textit{Pe},\varLambda )=(500,100)$ in figure 7 e). Importantly, this shows that, even if secondary vortices develop in the flow, it does not necessarily mean that compositional stripes will be seen at the surface of the experiment, because the overturning effect can be balanced by lateral diffusion. Many apparently mundane flows may therefore have more complex underlying motion that is hard to spot. Note also that, if the recirculation becomes too strong, then the sharp interface between large and small particles becomes unstable and forms a breaking-size-segregation wave, which is not captured by the depth-averaged theory (e.g. for $(\textit{Pe},\varLambda )=(500,3/2)$ in figure 7 f–j). Nevertheless, the exact steady-state solution (5.13) is still able to approximate the position of the breaking wave, and accurately captures the interface outside of its immediate vicinity. The depth-averaged theory may therefore still provide useful predictions even outside the region of parameter space for which it was originally derived.

5.5. Solving for the time-dependent evolution in the absence of diffusion

It is also possible to derive an exact solution for the evolving interface in the absence of diffusion. In this case the depth-averaged segregation equation (5.7) reduces to

(5.16) \begin{align} \frac {\partial \tilde {\eta }}{\partial \tilde {t}}+\frac {1}{\varLambda }\frac {\partial }{\partial \tilde {y}}\left (\sin \left ({\frac {\pi \tilde {y}}{\tilde {W}}}\right )(\tilde {\eta }^2-\tilde {\eta })\right )=0. \end{align}

In quasi-linear form this is equivalent to

(5.17) \begin{align} \frac {\partial \tilde {\eta }}{\partial \tilde {t}}+\frac {1}{\varLambda }\sin \left ({\frac {\pi \tilde {y}}{\tilde {W}}}\right )(2\tilde {\eta }-1)\frac {\partial \tilde {\eta }}{\partial \tilde {y}}=-\frac {1}{\varLambda }\frac {\pi }{\tilde {W}}\cos \left ({\frac {\pi \tilde {y}}{\tilde {W}}}\right )(\tilde {\eta }^2-\tilde {\eta }), \end{align}

which can be solved by the method of characteristics. Introducing the variable $\lambda$ to parameterise the characteristic curve $(\tilde {t}(\lambda ),\tilde {y}(\lambda ))$ , then the rate of change of the interface along the characteristic curve is

(5.18) \begin{align} \frac {\mathrm{d}\tilde {\eta }}{\mathrm{d}\lambda }=\frac {\partial \tilde {\eta }}{\partial \tilde {t}}\frac {\mathrm{d}\tilde {t}}{\mathrm{d}\lambda } +\frac {\partial \tilde {\eta }}{\partial \tilde {y}}\frac {\mathrm{d}\tilde {y}}{\mathrm{d}\lambda }. \end{align}

Comparing (5.17) with (5.18) implies the characteristic equations are

(5.19) \begin{align} \frac {\mathrm{d}\tilde {t}}{\mathrm{d}\lambda } &=1, \\[-12pt] \nonumber \end{align}

(5.20) \begin{align} \frac {\mathrm{d}\tilde {y}}{\mathrm{d}\lambda }&=\frac {1}{\varLambda }\sin \left ( \frac {\pi \tilde {y}}{\tilde {W}} \right )(2\tilde {\eta }-1), \\[-12pt] \nonumber \end{align}

(5.21) \begin{align} \frac {\mathrm{d}\tilde {\eta }}{\mathrm{d}\lambda }&=-\frac {1}{\varLambda }\frac {\pi }{\tilde {W}}\cos \left ( {\frac {\pi \tilde {y}}{\tilde {W}}} \right )(\tilde {\eta }^2-\tilde {\eta }). \\[10pt] \nonumber \end{align}

Since $\lambda$ does not appear on the right-hand side of this system, the chain rule and (5.19) can be used to eliminate it from (5.20) and (5.21), to generate the reduced system of two ODEs

(5.22) \begin{align} \frac {\mathrm{d}\tilde {y}}{\mathrm{d}\tilde {t}}&=\frac {1}{\varLambda }\sin \left ({\frac {\pi \tilde {y}}{\tilde {W}}}\right )(2\tilde {\eta }-1), \\[-12pt] \nonumber \end{align}

(5.23) \begin{align} \frac {\mathrm{d}\tilde {\eta }}{\mathrm{d}\tilde {t}}&=-\frac {1}{\varLambda }\frac {\pi }{\tilde {W}}\cos \left ({\frac {\pi \tilde {y}}{\tilde {W}}}\right )(\tilde {\eta }^2-\tilde {\eta }), \\[10pt] \nonumber \end{align}

which are subject to the initial conditions

(5.24) \begin{align} \tilde {y}=\tilde {y}_0, \qquad \hbox{and}\qquad \tilde {\eta }=\tilde {\eta }_0, \qquad \hbox{at}\qquad \tilde {t}=0. \end{align}

The chain rule can also be used to eliminate $\tilde {t}$ , and the resulting separable ODE can be integrated to show that

(5.25) \begin{align} \sin \left (\frac {\pi \tilde {y}}{\tilde {W}}\right )(\tilde {\eta }^2-\tilde {\eta })=\psi _0, \end{align}

where $\psi _0$ is a constant of integration. Since the streamfunction for the non-dimensional velocity field (2.21)–(2.22), evaluated at the interface $\tilde {\eta }$ , is

(5.26) \begin{align} \psi =\sin {\left (\frac {\pi \tilde {y}}{\tilde {W}}\right )}(\tilde {\eta }^2-\tilde {\eta }), \end{align}

equation (5.25) is equivalent to

(5.27) \begin{align} \psi =\psi _0. \end{align}

The characteristic curves therefore correspond to the closed streamlines of the bulk flow field. To solve for the position of the interface along the streamline at any given time, it is actually easier to solve the whole system (5.22)–(5.23) numerically. The initial interface position (5.24) is therefore parameterised by a series of points, whose subsequent motion is calculated by solving (5.22)–(5.23). Figure 8 and movie 6 show a comparison of the time-dependent evolution of the interface, solved by the method of characteristics, and the full two-dimensional solution for $\phi ^s$ for parameters $(Pe,\varLambda )=(\infty ,100)$ . The depth-averaged solution accurately picks off the steepening of the sharp interface up to approximately $\tilde {t}=140$ non-dimensional time units, when the interface in the non-depth-averaged simulation breaks to form a very thin breaking-size-segregation wave. At subsequent times the interface position can still be calculated by the method of characteristics, but the solution is multi-valued. A shock, or discontinuity in the interface height, then develops, which connects two smooth single-valued branches of the solution.

Figure 8.

Comparison of the one-dimensional interface height $\tilde {\eta }$ (black line) calculated using the method of characteristics, and the full two-dimensional solution for $\phi ^s$ in a single clockwise-rotating vortex cell, with the velocity field (2.21)–(2.22), and for a 30 : 70 mix and $(\textit{Pe},\varLambda )=(\infty ,100)$ . The solutions are shown at (a) $\tilde {t}=0$ , (b) $\tilde {t}=100$ , (c) when the shock (red line) first forms at $\tilde {t}=140$ , (d) $\tilde {t}=150$ , (e) $\tilde {t}=170$ , (f) $\tilde {t}=250$ , (g) steady state $\tilde {t}=510$ and (h) when the lobes cross the fitted shock as they continue to be swept round after steady state, $\tilde {t}=650$ . After the shock has formed the lobes of the interface are shown as dashed lines. Movie 6 shows the complete time-dependent solution.

Figure 8. Long description

(a) Inversely-graded distribution with large-green material above small-white material with a horizontal interface. (b)-(c) The interface becomes increasingly waved until it breaks, after which the depth-averaged solution (black line) requires a shock to be fitted (red line) between the two lobes of the depth-averaged solution. As the shock grows until it is the full flow depth, the lobes continue to be swept round by the bulk. The depth-averaged solution accurately captures the two-dimensional interface.

The growth and dynamics of the shock are controlled by a jump condition (Whitham Reference Whitham1974; Chadwick Reference Chadwick1976; Gray & Kokelaar Reference Gray and Kokelaar2010a ,Reference Gray and Kokelaar b ). The jump condition is derived from the conservative form of the depth-averaged segregation equation (5.16), and in this case is

(5.28) \begin{align} -v_n\left[\!\left[\tilde {\eta }\right ]\!\right ]+ {\left [\!\left [\frac {1}{\varLambda }\sin \left (\frac {\pi \tilde {y}}{\tilde {W}}\right )(\tilde {\eta }^2-\tilde {\eta }) \right ]\!\right ]} =0, \end{align}

where $v_n$ is the shock propagation speed and the jump bracket $[\![ f ]\!]=f_+-f_-$ is the difference of the enclosed quantity on the forward and rearward sides of the shock. This implies the shock speed is

(5.29) \begin{align} v_n=\frac {1}{\varLambda }\sin {\frac {\pi \tilde {y}}{\tilde {W}}}(\tilde {\eta }_++\tilde {\eta }_- -1)=\frac {d\tilde {y}}{d\tilde {t}}, \end{align}

which requires a knowledge of $\tilde{\eta} _+$ and $\tilde{\eta} _-$ across the shock, making it difficult to use directly. Instead, a shock fitting procedure is used to determine the shock position. This ensures that the lobes formed by the multi-valued characteristic solution, and which are cutoff by the shock, are of equal area (Whitham Reference Whitham1974). The shock forms at the point where the interface becomes vertical, and grows in height as the secondary vortex continues to steepen the interface between large and small particles. By $\tilde {t}=510$ non-dimensional time units, the shock reaches the full height of the cell and a steady state develops. Throughout the entire time-dependent evolution, the depth-averaged theory closely matches the position of the interface calculated with the full two-dimensional segregation theory.

7.1. A two-dimensional depth-averaged segregation model

For simplicity, it is assumed that the flow is of constant depth $H$ , is sharply inversely graded, with the particle-size distribution given by (5.2), and hence $H\overline {\phi ^s}=\eta$ , just as in (5.3) of the one-dimensional model developed in § 5.1. If the diffusion coefficient is also constant, then (4.4) reduces to

(7.1) \begin{align} \frac {\partial \eta }{\partial t}+ \frac {\partial }{\partial x}(H\overline {\phi ^s u}) + \frac {\partial }{\partial y}(H\overline {\phi ^s v}) = \frac {\partial }{\partial x}\left (D \frac {\partial \eta }{\partial x}\right )+\frac {\partial }{\partial y}\left (D \frac {\partial \eta }{\partial y}\right )\!. \end{align}

This is subject to an initial condition

(7.2) \begin{align} \eta (x,y,t)=\eta _0(x,y),\quad \hbox{for}\quad (x,y)\in [0,\mathscr{L}\,]\times [-y_{{wall}},y_{{wall}}] , \end{align}

where $\mathscr{L}$ is the length of the chute and the sidewalls are located at $y=\pm y_{{wall}}$ . At the inflow ( $x=0$ ), a constant interface height $\eta _{\mathrm{o}}$ is prescribed across the chute

(7.3) \begin{align} \eta (0,y,t)=\eta _{\mathrm{o}},\quad \hbox{in}\quad y\in [-y_{{wall}},y_{{wall}}], \end{align}

and along the sidewalls the no-flux condition implies that

(7.4) \begin{align} \frac {\partial \eta }{\partial y}=0,\quad \hbox{at}\quad y=\pm y_{{wall}}. \end{align}

The grains are allowed to freely flow out of the chute at $x=\mathscr{L}$ , which is equivalent to the condition

(7.5) \begin{align} \frac {\partial \eta }{\partial x}=0,\quad \hbox{at}\quad x=\mathscr{L}. \end{align}

In order to close the model, it remains to determine expressions for the depth-averaged flux of small particles in the downslope and cross-slope directions.

7.2. The downslope velocity profile and flux function

Following Gray & Kokelaar (Reference Gray and Kokelaar2010a ,Reference Gray and Kokelaar b ), a simple linear velocity profile is prescribed for the downslope velocity through the depth of the flow, i.e.

(7.6) \begin{align} u=\overline {u}\left (\alpha +2(1-\alpha )\frac {z}{H}\right ), \end{align}

where the depth-averaged downslope velocity $\overline {u}=\overline {u}(y)$ is still able to vary across the channel. This should be viewed as a first-order Taylor series approximation to the true velocity profile. The parameter $\alpha \in [0,1]$ sets the amount of shear through the avalanche depth, with $\alpha =1$ corresponding to plug flow and $\alpha =0$ representing simple shear. In the centre of the channel, (7.6) implies that the downslope surface and depth-averaged velocities, $u^s_{\textit{centre}}$ and $\bar {u}_{\textit{centre}}$ , are related by the formula

(7.7) \begin{align} u^s_{\textit{centre}}=\overline {u}_{\textit{centre}}(2-\alpha ), \end{align}

which can be used to infer the value of

(7.8) \begin{align} \alpha =2-\frac {u^s_{\textit{centre}}}{\bar {u}_{\textit{centre}}}\simeq 0.727 \end{align}

from the experimental measurements in § 6. This value of $\alpha$ implies that there is shear through the flow depth as well as a significant amount of basal slip, which is consistent with the formation of a low solids volume fraction basal Leidenfrost layer (Wang et al. Reference Wang, Jing, Kwok, Sobral, Weinhart and Thornton2025).

Substituting the downstream velocity profile (7.6) into the integral for the downstream flux of small particles (4.3b ) implies that

(7.9) \begin{align} H\overline {\phi ^s u}=\overline {u}\left (\eta -(1-\alpha )\left (\eta -\frac {\eta ^2}{H}\right )\right ). \end{align}

Despite its simplicity, the downslope velocity profile (7.6) and the associated downslope flux of small particles (7.9) capture many essential features of the flow, and have been used to model the accumulation of large particles at avalanche fronts, as well as segregation-induced flow fingering in fully coupled depth-averaged models (Gray & Kokelaar Reference Gray and Kokelaar2010a ,Reference Gray and Kokelaar b ; Woodhouse et al. Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012; Baker et al. Reference Baker, Johnson and Gray2016).

The experimental measurements of the surface velocity, shown in figure 10, suggest that the sides of the flow move slower than the centre. To parameterise this, the downslope surface velocity is assumed to take the form

(7.10) \begin{align} u^s(y)= \begin{cases} u^s_{\textit{centre}}+(u^s_{{wall}}-u^s_{\textit{centre}})\left (\frac {y+y_{{BL}}}{-y_{{wall}}+y_{{BL}}} \right )^2, & -y_{{wall}}\leqslant y\lt -y_{{BL}},\\[4pt]u^s_{\textit{centre}}, & -y_{{BL}}\leqslant y \leqslant y_{{BL}},\\[4pt]u^s_{\textit{centre}}+(u^s_{{wall}}-u^s_{\textit{centre}})\left (\frac {y-y_{{BL}}}{y_{{wall}}-y_{{BL}}} \right )^2, & y_{{BL}}\lt y \leqslant y_{{wall}}, \end{cases} \end{align}

where the surface wall-slip velocity is equal to $u^s_{{wall}}$ at $y=\pm y_{{wall}}$ and the central plug of material moves downslope with constant surface velocity $u^s_{\textit{centre}}$ . The flow is symmetric about $y=0$ and has two boundary layers that smoothly connect to the central velocity at $y=\pm y_{\textit{BL}}$ . This approximation is illustrated in figure 10(d) for the case $y_{{BL}}=5W$ . It neglects the subtle variations in the downslope surface velocity across the slope, which correlate with the flow bands/ridges, but is a reasonable first approximation. The reconstructed two-dimensional downslope velocity field $u$ is shown using dimensional variables in figure 11(b).

Figure 11.

(a) Comparison between the two-dimensional steady-state solution for the small-particle concentration $\phi ^s(y,z)$ (contour plot) and the one-dimensional depth-averaged model interface height $\eta (y)$ (magenta line) at $t=10$ s for the experimental parameters in table 1. This is equivalent to the case $(\textit{Pe},\varLambda )=(49.04,2.0155)$ . The two solutions are computed with a prescribed velocity field ${\boldsymbol u}(y,z)$ , which is independent of the downslope coordinate $x$ . Panels (b, c, d) show contour plots of the downslope, cross-slope and slope-normal velocity components $(u,v,w)$ , respectively. The contour lines in (a) are shown in the range $\phi ^s\in [0.1,0.9]$ in 0.08 increments.

Figure 11. Long description

(a) Comparison between the two-dimensional steady-state solution (contour plot) and the one-dimensional depth-averaged model interface height (magenta line) at ten seconds for the experimental parameters in table 1, equivalent to the case (Pe,Λ)= (49.04,2.0155). The figure shows 10 cells of width W and a wider boundary-layer cell adjacent to each sidewall. The depth-averaged solution accurately predicts the interface between small and large replacing the diffuse breaking-size-segregation waves with a shock in the small particle concentration. The solutions are computed with a prescribed velocity field which is independent of the downslope coordinate x. Panels (b,c,d) show contour plots of the downslope, cross-slope and slope normal velocity components respectively. The downslope velocity is largest at the free-surface and decreases as you approach the base and sidewalls. The cross-slope and normal velocity combine to form the alternating direction vortical velocity cells.

7.3. The cross-slope flux function

Figure 10 shows that the sidewall boundary layers have an upwelling band of fines-rich material immediately adjacent to the wall, and that the ends of the boundary layers are approximately at the centre of the first large–rich band. Moreover, in the centre of the flow, where the downslope surface velocity is plug-like, there are five complete stripes measured from centre to centre of each large–rich band. This implies that in the centre of the flow there must be ten counter-rotating secondary vortices, which are flanked on either side by a wider secondary vortex in each wall boundary layer. The cross-slope velocity field can therefore be parameterised by

(7.11) \begin{align} v= \begin{cases} -\varOmega _{{BL}}\sin {\left (\frac {\pi }{W_{{BL}}}(y+5W)\right )}(2z-H), & -y_{{wall}}\leqslant y\lt -5W,\\[12pt]\varOmega \sin {\left (\frac {\pi y}{W}\right )}(2z-H), & -5W\leqslant y \leqslant 5W,\\[12pt]-\varOmega _{{BL}}\sin {\left (\frac {\pi }{W_{{BL}}}(y-5W)\right )}(2z-H), & 5W\lt y \leqslant y_{{wall}}, \end{cases} \end{align}

where $W_{{BL}}$ is the width of the sidewall boundary layer and $\varOmega _{{BL}}$ is its rotation rate. Note that the prescribed cross-slope velocity in the central channel is exactly the same as in (2.12). The cross-slope velocity in the two boundary layers has a similar functional form, but is implicitly assumed to rotate in the opposite direction to the adjacent cell (hence the minus sign). The experiments in § 6 suggest that the edge of the downslope velocity boundary layers is the same as the edge of the (wider) sidewall secondary vortices, i.e. that

(7.12) \begin{align} y_{{BL}}=5W, \end{align}

and hence the width of the boundary layer is

(7.13) \begin{align} W_{{BL}}=y_{{wall}}-5W. \end{align}

In order to determine $\varOmega _{{BL}}$ it is necessary to consider the slope-normal velocity $w$ . Since the velocity field is incompressible and satisfies (2.9), the normal velocity can be reconstructed by differentiating (7.11) with respect to $y$ and integrating with respect to $z$ , subject to the boundary condition that $w=0$ at $z=0$ , to give

(7.14) \begin{align} w= \begin{cases} \frac {\varOmega _{{BL}}\pi }{W_{{BL}}}\cos {\left (\frac {\pi }{W_{{BL}}}(y+5W)\right )}z(z-H), & -y_{{wall}}\leqslant y\lt -5W,\\[12pt]-\frac {\varOmega \pi }{W}\cos {\left (\frac {\pi y}{W}\right )}z(z-H), & -5W\leqslant y \leqslant 5W,\\[12pt]\frac {\varOmega _{{BL}}\pi }{W_{{BL}}}\cos {\left (\frac {\pi }{W_{{BL}}}(y-5W)\right )}z(z-H), & 5W\lt y \leqslant y_{{wall}}. \end{cases} \end{align}

The normal velocity $w$ must vary smoothly across the interfaces between the boundary layers and the central region. Evaluating (7.14) at $y=\pm 5W$ it follows that the counter-rotating vortices in the boundary layers are consistent with those in the central flow, provided

(7.15) \begin{align} \varOmega _{{BL}}=\left (\frac {W_{{BL}}}{W}\right )\varOmega . \end{align}

Since the boundary layers are wider than the secondary vortices in the central part of the flow, they must rotate faster to be compatible (i.e. to prevent a contact discontinuity in $w$ forming at $y=\pm 5W$ ). The components of the reconstructed velocity field are shown in figure 11.

Finally, multiplying the cross-slope velocity (7.11) by the small-particle concentration (5.2) and integrating through the depth, implies that the cross-slope flux function is

(7.16) \begin{align} H\overline {\phi ^s v}= \begin{cases} -\varOmega _{{BL}}\sin {\left (\frac {\pi }{W_{{BL}}}(y+5W)\right )}(\eta ^2-H\eta ), & -y_{{wall}}\leqslant y\lt -5W,\\[12pt]\varOmega \sin {\left (\frac {\pi y}{W}\right )}(\eta ^2-H\eta ), & -5W\leqslant y \leqslant 5W,\\[12pt]-\varOmega _{{BL}}\sin {\left (\frac {\pi }{W_{{BL}}}(y-5W)\right )}(\eta ^2-H\eta ), & 5W\lt y \leqslant y_{{wall}}. \end{cases} \end{align}

Note that, in the centre of the flow, (7.16) is identical to the cross-slope flux function (5.4) derived in § 5. Before using the two-dimensional depth-averaged theory (7.1) to model the stripe formation down the chute, it is instructive to compute the depth-averaged solution in the chute geometry shown in figure 11 assuming spatial uniformity in the downslope direction $x$ , and to compare it against the two-dimensional non-depth-averaged theory of § 2. Both theories still need a number of experimental parameters to be determined.

7.4. Experimental parameter estimation

In order to solve the full segregation-advection-diffusion equation (2.10), the segregation rate $q$ and the diffusivity $D$ are approximated using the experimental scaling laws of Utter & Behringer (Reference Utter and Behringer2004) and Trewhela et al. (Reference Trewhela, Ancey and Gray2021). In its generalised form, Utter & Behringer’s (Reference Utter and Behringer2004) law implies that

(7.17) \begin{align} D= \mathcal{A}\dot \gamma \bar d^2, \end{align}

where $\mathcal{A}$ is a universal constant specified in table 1, $\dot \gamma$ is the shear rate and $\bar {d}$ is the average local particle size. The downslope flow is much larger than the secondary flow, so to leading order the shear rate

(7.18) \begin{align} \dot \gamma =\left |\frac {\partial u}{\partial z}\right |=2(1-\alpha )\frac {\bar {u}}{H}, \end{align}

which is constant through the flow depth, provided the linear downslope velocity profile (7.6) is assumed. The constant $\alpha$ is determined in (7.8). The average local particle size

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

where $d^l$ and $d^s$ are the diameters of the large and small particles, respectively. This is approximated using the values of the particle diameters and concentrations given in table 1. It follows that the diffusivity in the central channel is approximately equal to

(7.20) \begin{align} D=4.1607\times 10^{-7}\,\mathrm{m}^2\,\mathrm{s}^{-1}. \end{align}

In the simulations this constant value is assumed throughout the entire flow, even in the sidewall boundary layers, where the shear rate and slip may be modified.

Trewhela et al.’s (Reference Trewhela, Ancey and Gray2021) experimental scaling law for the segregation rate

(7.21) \begin{align} q=\frac {\mathcal{B}\dot \gamma \bar d^2\mathcal{F}}{\mathcal{C}\bar d+\varPhi (H-z)}, \end{align}

captures the observed shear-rate and pressure dependence through the flow depth, as well as the particle-size ratio $R=d^l/d^s$ dependence in the factor

(7.22) \begin{align} \mathcal{F}=(R-1)+\mathcal{E}\phi ^l(R-1)^2. \end{align}

In this law $\mathcal{B}$ , $\mathcal{C}$ and $\mathcal{E}$ are universal constants that are given in table 1. In order to find a representative constant value for $q$ , the scaling law is evaluated at $z=H/2$ and $\bar {d}$ is approximated by a global average (as above). Substituting the values given in table 1 implies that the segregation velocity

(7.23) \begin{align} q=0.0076\,\mathrm{m}\,\mathrm{s}^{-1}. \end{align}

The large and small particles are therefore able to segregate from one another at a maximum velocity of 7.6 mm s $^{-1}$ . This is quite fast, but is a consequence of the thinness of the flow.

Measuring the cross-slope surface velocity is difficult, because it is so much smaller than the downslope velocity. However, it is possible to use the depth-averaged segregation equation (7.1) to infer an approximate value. At steady state, the dominant balance is between the downslope flux gradient and the cross-slope flux gradient

(7.24) \begin{align} \frac {\partial }{\partial x} (H\overline {\phi ^s u})\sim \frac {\partial }{\partial y}(H\overline {\phi ^s v}). \end{align}

The experiments in § 6 suggest that order-unity changes in the downslope direction occur on a length scale $L=2.3$ m, but order-unity changes across the slope are of the same order as the flow depth $H$ . Moreover, assuming that the downslope flux $H\overline {\phi ^s u}$ in (7.9) can be approximated by $\bar {u}\eta$ and the cross-stream flux $H\overline {\phi ^s v}$ in (7.16) by $\varOmega (H\eta -\eta ^2)$ , then (7.24) implies that

(7.25a,b) \begin{align} \frac {\bar {u}\eta }{L}\sim \frac {\varOmega (H\eta -\eta ^2)}{H},\quad \Rightarrow \quad \varOmega \sim \frac {\bar {u}H}{L(H-\eta )}. \end{align}

Substituting the measured value of $\bar {u}$ from (6.1) and assuming that the interface lies at approximately $\eta =0.7H$ then it follows from (7.25b ) that the rotation rate in the centre of the channel is approximately

(7.26) \begin{align} \varOmega =1.389\,\mathrm{s}^{-1}. \end{align}

This implies that maximum surface cross-slope velocities are of the order of $\varOmega H\sim 3.8$ mm s $^{-1}$ . It follows from (7.15) that the faster rotation rate of the secondary vortices in the sidewall boundary layers is

(7.27) \begin{align} \varOmega _{\textit{BL}}= 3.8892\,\mathrm{s}^{-1}, \end{align}

implying a maximum surface velocity of 10.5 mm s $^{-1}$ . This is consistent with the bi-disperse roll wave experiments of Viroulet et al. (Reference Viroulet, Baker, Rocha, Johnson and Gray2018), where fines-rich boundary layers also formed. In particular, Viroulet et al.’s (Reference Viroulet, Baker, Rocha, Johnson and Gray2018) movie 1 shows that, as each large–rich roll wave front overtakes the material in front of it, large particles are pushed towards the sidewall, but are then advected back towards the centre of the channel prior to the arrival of the next wave. Finally, using their definitions (2.17a,b ) and (5.14) the non-dimensional parameters are

(7.28) \begin{align} \textit{Pe}=49.04,\quad \varLambda = 2.0155,\quad \varXi = 14.3455. \end{align}

All the dimensional and non-dimensional numbers needed to simulate the experimental flow, in both the depth-averaged and non-depth-averaged models, are summarised in table 1.

7.5. Comparison between the one- and two-dimensional representations of the experiment

In this section the particle-size distribution is assumed to be spatially uniform in the downslope direction. It follows that the depth-averaged model reduces to a one-dimensional problem across the chute, while the non-depth-averaged theory reduces to a two-dimensional problem in a $(\tilde {y},\tilde {z})$ section. The flow is assumed to be initially sharply segregated in an inversely graded configuration (large particles above small ones) with an interface located at

(7.29) \begin{align} \eta _0(y)=0.7 H. \end{align}

This is equivalent to an initial particle-size distribution of the form

(7.30) \begin{align} \phi ^s_0(y)= \begin{cases} 0,\quad z\gt 0.7 H,\\ 1,\quad z\leqslant 0.7 H, \end{cases} \end{align}

to again reflect the 30 : 70 mix of large (green) and small (white) particles used in the experiment. The solutions are computed in dimensional variables to ensure that they can easily be compared with the experimental observations. The results are shown in figure 12 and movie 8. Focusing first on the two-dimensional theory, at $t=0.5$ s the interface between large and small particles is still very sharp, but it has been rotated by the secondary flow to create a wavy interface. The large particles concentrate in the downwelling sections and small particles concentrate in the upwelling sections. The amplitude of the wavy interface grows with increasing time, with parts of the interface becoming progressively steeper until they break at $t\simeq 1$ s to form a series of small diffuse breaking-size-segregation waves (Thornton & Gray Reference Thornton and Gray2008; Gray & Ancey Reference Gray and Ancey2009; Edwards et al. Reference Edwards, Rocha, Kokelaar, Johnson and Gray2023). These grow in height and width until they occupy the whole depth of the flow, and white stripes of small particles become visible at the surface at $t\simeq 1.5$ s. The breaking-size-segregation waves in the boundary layers are larger and wider than those in the centre of the flow, which is a consequence of the increased rotation rate $\varOmega _{{BL}}$ of the secondary vortices. The central diffuse breaking waves rapidly attain a steady state, as shown in movie 8. However, the diffusion is not strong enough to immediately damp out the oscillation of the larger breaking waves in the sidewall boundary layers, and they oscillate for some time. A total of six stripes develop across the chute, reflecting the underlying secondary flow, and in agreement with the experimental results in figures 9 and 10. The stripes adjacent to the sidewall boundary layers are wider than those in the central channel, which matches the observations in figure 9(b). Throughout the entire time-dependent evolution, the one-dimensional depth-averaged theory agrees remarkably well with the interface position predicted by the full theory. The depth-averaged model cannot capture the structure of the diffuse breaking-size-segregation waves, but instead represents them by a sharp transition that accurately picks off their approximate location. The depth-averaged theory therefore provides a very useful reduced representation of the much more complex two-dimensional small-particle concentration field, and has the advantage that it can be used to efficiently model the development of the stripes down the entire chute.

Figure 12.

(a–f) Comparison between the two-dimensional time-dependent solutions for the small-particle concentration $\phi ^s(y,z,t)$ (contour plot) and the one-dimensional depth-averaged model interface $\eta (y,t)$ (magenta line) at a sequence of times and for the experimental parameters in table 1. This is equivalent to the case $(\textit{Pe},\varLambda )=(49.04,2.0155)$ . The flow assumes ten central secondary vortices of width $W=0.5$ cm and a boundary layer on either side with a vortex cell of width $W_{{BL}}=1.4$ cm, as illustrated in figure 11. The boundaries of the secondary vortices are indicated by the tick marks along the $y$ axis. Movie 8 shows an animation of the time-dependent evolution of the solution. The contour lines are shown in the range $\phi ^s\in [0.1,0.9]$ in 0.08 increments.

Figure 12. Long description

(a-f) Comparison of the two-dimensional time-dependent solutions for the small-particle concentration (contour plot) and the one-dimensional depth-averaged model interface (magenta line) at a sequence of times and for the experimental parameters in table 1.The flow assumes ten central secondary vortices of width W and a wider boundary layer cell adjacent to each side wall. Starting from an inversely-graded profile in (a) the velocity field rotates the interface which becomes increasingly wavy in (b)-(c) before it steepens to vertical and breaks in (d). The breaking waves oscillate back and forth (e) before reaching steady state in (f) where the boundary layer cells have slightly wider breaking waves than the interior cells. The breaking waves seperate pure regions of small and large material forming 6 large-green stripes at the free-surface with boundary layer stripes slightly wider than interior.

7.6. Depth-averaged solutions for the stripe development along the chute

The two-dimensional depth-averaged segregation equation (7.1), with the flux functions (7.9) and (7.16), is solved in the domain $(x,y)\in [0,2.84\,\mathrm{m}]\times [-y_{{wall}},y_{{wall}}]$ . The domain is initially filled with material of constant height $H$ that is moving downslope and simultaneously undergoing secondary recirculation. The velocity field $\boldsymbol u$ is assumed to be steady and fully developed in the downstream direction, with velocity components $(u,v,w)$ given by (7.6), (7.11) and (7.14), respectively, and illustrated in figure 11. This is implicitly encoded in the flux functions (7.9) and (7.16).

Material is continuously supplied at the inflow at $x=0$ , and is assumed to be uniformly distributed across the chute. Since the inflow is implicitly assumed to be inversely graded, with all the large particles on top of the small ones, the small-particle layer needs to be thicker in order to transport the correct flux down the chute. Assuming the downslope velocity profile (7.6), the downslope small-particle flux

(7.31) \begin{align} F^s=\int _0^H \phi ^s u\,\mathrm{d}z=\overline {u}\left (\alpha \eta +(1-\alpha )\frac {\eta ^2}{H}\right )\!, \end{align}

and the total flux

(7.32) \begin{align} F=\int _0^H u\,\mathrm{d}z=H\overline {u}. \end{align}

The correct proportion of small grains in the flowing mixture is achieved when the ratio of the small-particle flux $F^s$ to the total flux $F$ is equal to the concentration of the static mixture of grains $\phi ^s_{\textit{static}}$ , i.e. $F^s/F=\phi ^s_{\textit{static}}$ . Substituting (7.31) and (7.32) into this expression yields the quadratic equation

(7.33) \begin{align} (1-\alpha )\left (\frac {\eta }{H}\right )^2+\alpha \left (\frac {\eta }{H}\right )-\phi ^s_{\textit{static}}=0, \end{align}

which can be solved to determine the corrected inflow height

(7.34) \begin{align} \eta _{\mathrm{o}}(y)=\left (\frac {-\alpha +\sqrt {\alpha ^2+4(1-\alpha )\phi ^s_{\textit{static}}}}{2(1-\alpha )}\right )H. \end{align}

For the 30 : 70 mix the static concentration $\phi ^s_{\textit{static}}=0.7$ . It follows that for the velocity profile (7.6) with $\alpha$ given by (7.8), the corrected inflow height is

(7.35) \begin{align} \eta _{\mathrm{o}}(y)=0.751\,H, \end{align}

which is only a 7 % correction, because the downslope velocity has a lot of basal slip. The material in the chute at $t=0$ is also assumed to have the same initial interface height, i.e.

(7.36) \begin{align} \eta _0(x,y)=0.751\,H. \end{align}

Note that the depth-averaged theory implicitly assumes that the bidisperse mixture is sharply inversely graded as it enters the chute. This may not be the case, despite some segregation probably occurring as the grains are sheared towards the hopper exit. However, balancing the downslope flux of small particles with the segregation term in (2.5) suggests that the dominant balance is

(7.37) \begin{align} \frac {\partial }{\partial x}(\phi ^s u) \sim \frac {\partial }{\partial z} \big(q \phi ^s(1-\phi ^s)\big). \end{align}

If vertical gradients occur on the length scale of the flow depth $H=0.0027$ m, a typical depth-averaged downslope velocity $\bar {u}=$ 0.9584 m s $^{-1}$ and a typical segregation velocity $q=0.0076$ m s $^{-1}$ then the approximate downslope distance $L_{\textit{seg}}$ for complete vertical segregation through the flow depth is

(7.38) \begin{align} L_{\textit{seg}}\sim \frac {H\bar {u}}{q}=34\,\hbox{cm}. \end{align}

The segregation into inversely graded layers will therefore have occurred in the first 12 % of the chute, and may well have completed earlier if there is segregation prior to the grains exiting the hopper. Moreover, even if segregation is continuing to occur as the grains are flowing downslope, they are already being rotated by the secondary flow, so large and small particles will begin to concentrate in the downwelling and upwelling regions almost immediately.

Figure 13.

Contour plot showing the depth-averaged interface height $\eta (x,y)$ at $t=8$ s for the experimental parameters in table 1 and the velocity field in figure 11. The simulation domain is split into three regions (a) 0–1.04 m, (b) 0.91–1.95 m and (c) 1.8–2.84 m consistent with experiments shown in figure 9. The colour bar is skewed to make the colours representative of surface observations in a real flow. A time-dependent movie of the simulation is shown in movie 9 and figure 14 shows a series of cross-sections across the chute.

Figure 13. Long description

Contour plot showing the depth-averaged interface height at eight seconds for the experimental parameters in table 1 and the velocity field in figure 11. The simulation domain is split into three regions (a) 0-1.04 m, (b) 0.91-1.95 m and (c) 1.8-2.84 m consistent with experiments shown in figure 9. The colourbar is skewed to make the colours representative of surface observations in a real flow. In (a) material entering the chute is a uniform shade of green, with white stripes appearing near the sidewalls at x=0.5 m which thicken to their maximum width x=1.2 m . Some surface stripes begin to form between x=0.8-1 m. (b) The surface stripes gradually widen and sharpen to their steady-state width by x=2.3 m (c) and continue until material leaves the end of the chute at x=2.84 m.

Figure 14.

(a–g) A series of $(y,z)$ sections reconstructed from the steady-state depth-averaged solution $\eta (x,y)$ at $t=8$ s shown in figure 13 and computed with the experimental parameters in table 1. Each section shows the interface position $\eta (y)$ , which lies along the interface between large particles (shaded green) and small particles (shaded white). The boundaries of the secondary vortices are indicated by the tick marks along the $y$ axis. Movie 10 shows the time-dependent development of the solution in each cross-section.

Figure 14. Long description

Cross-sections of the interface solution through the flow depth at different downstream distances reconstructed from the steady-state depth-averaged solution in figure 13. (a) At 0 m, the boundary condition enforces a inversely-graded inflow. (b)-(c) Between 0.2-0.5 m, the interface becomes wavy with the free-surface remaining green. (d) At 0.7 m, the interface becomes more wavy, small-white material has broken onto the free surface by the sidewalls but the surface remains green elsewhere. (e)-(f) Between 1-1.5 m, the interfaces become steeper, with small-white material breaking onto the free-surface forming stripes that widen and sharpen with downstream distance. (g) By the end of the chute, the stripes have reached steady-state with 6 distinct large-green stripes and 7 distinct small-white stripes both of which are slightly wider in the boundary layer.

Movie 9 shows a top down animated contour plot of the time-dependent solution along the chute, and the steady state is shown in figure 13. In these plots the colour bar is chosen to reflect what a viewer sees from above. Note that when there is more than a monolayer of large green particles at the experimental flow surface, the surface appears green, regardless of how deep the large–small-particle interface $\eta$ lies beneath it. To approximate this effect, the colour bar is therefore skewed towards green. In particular, the gradient in the colour bar, from green to white, is chosen so that it occurs when the interface height $\eta$ is within one large-particle diameter of the surface. It is also important to stress that in the experiment there is initially no material in the chute, so the time-dependent simulation and movie 7 are not directly comparable. However, the steady-state solution shown in figure 13 is comparable to figure 9.

Material entering the chute is a uniform shade of green, as in the experiments, and a white stripe begins to emerge adjacent to each sidewall at approximately $x=0.5$ m. At this point the material in the centre of the channel still appears uniformly green. However, the reconstructed cross-section of the small particle-size distribution across the channel at $x=0.5$ m, shown in figure 14(c), indicates that the secondary recirculation has already caused a wavy interface to form across the channel. This remains unseen both in the experiment and in the simulated top down view in figure 13. The white stripes adjacent to the sidewalls thicken with downstream distance and are approximately fully formed by 1.2 m downslope, after which they do not noticeably change width down the whole of the rest of the channel. All of this behaviour is in good agreement with the experiment, although the white stripes in the experiment appear to start to form slightly higher up the chute and reach their steady-state width somewhat earlier.

The first vestiges of the stripes that form in the central channel begin to emerge in the simulations at $x=0.8-1.0$ m in figure 13(a). This is approximately the same distance that higher and lower concentrations of large and small particles are seen in the experiments. The cross-section at $x=1$ m shows that the interface $\eta$ is very close to the free surface in the fines-rich upwelling regions, but it has not fully penetrated through to the surface (figure 14 e). Equally, in the large–rich downwelling sections, large particles have still not reached the base of the flow. By $x=1.5$ m these thin layers become progressively thinner (see figure 14 f) and the stripes become visible at the surface (figure 13 b) in accordance with the experiment. A total of six large–rich bands then form in the channel, as a direct consequence of the twelve counter-rotating vortices across the width of the flow (figure 11) and in accordance with the conceptual picture in figure 2. With increasing downstream distance the bands sharpen up and achieve their steady-state width approximately 2.3 m downslope, as suggested by the scaling argument in (7.25).

Movie 10 shows an animation of the time-dependent development of the stripes in each of the cross-sections shown in figure 14. Each of these cross-sections looks a little like those illustrated in figure 12, but here they are cross-sections at a fixed time and a series of different locations, rather than being at a fixed location and at variable times. In addition, movie 10 looks considerably simpler than movie 8, because the depth-averaged theory cannot represent the structure of the breaking-size-segregation waves. Nevertheless this simplified reconstruction of the particle-size distribution from the depth-averaged model captures many essential features of the development of the stripes.

In movie 10, the section at $x=0.2$ m starts off with a horizontal interface (reflecting the inversely graded initial conditions) and starts to develop a series of waves. However, by $t=0.2$ s the waves have already reached their maximum steady-state amplitude and do not develop further. In this wavy configuration the interface $\eta$ never reaches the free surface, and so the surface of the flow will appear green in the experiments. At $x=0.5$ m the waves become more pronounced and stop developing at approximately $t=0.5$ s. At the sides of the flow the interface $\eta$ is sufficiently close to the surface that the upwelling fine material becomes visible in the top down reconstruction in figure 13. In the centre of the channel, however, there is still no sign of the internal waviness in the surface view.

The sections at $x=0.7$ and $x=1$ m reach steady state in the centre of the channel at approximately $t=0.7$ and $t=1$ s, respectively, reflecting the fact that the depth-averaged downslope velocity $\bar {u}_{\textit{centre}}\simeq 0.9584$ m s $^{-1}$ in the centre of the channel. Note that there is some continued evolution of the interface in the sidewall boundary layers after this time, because the downslope velocity is lower there. At $x=1$ m the stripes start to become visible in the top down view, and they become progressively stronger with increasing downstream distance (see figure 13 and movie 9). At $x=1.5$ and $x=2.84$ m the initial evolution of the interface looks like that in the time-dependent movie of the spatially uniform problem in figure 12 and movie 8. This is because the initial condition (7.36) implies that the initial evolution is precisely uniform in the downslope direction $x$ , and it is not until $t=1.5$ and $t=2.84$ s that the influence of the upstream boundary condition reaches these locations, respectively. Interestingly, after this time both profiles adjust slightly from the solutions that would be generated by a purely spatially uniform problem, as the lateral diffusion pushes some of the large grains that were initially in the sidewall vortices into the neighbouring internal vortex cells. Overall the depth-averaged model (7.1) does a remarkably good job of capturing both the qualitative and quantitative development of stripes observed in experiments in figure 9. In particular, it predicts the formation of alternating longitudinal bands of pure large and pure small particles that extend all the way through the flow depth. This is consistent with what Pearse et al. (2026) observed in the curved chute deposit, shown in their figure 3.