Self-channelisation and levee formation in monodisperse granular flows

Dense granular flows can spontaneously self-channelise by forming a pair of parallel-sided static levees on either side of a central flowing channel. This process prevents lateral spreading and maintains the flow thickness, and hence mobility, enabling the grains to run out considerably further than a spreading flow on shallow slopes. Since levees commonly form in hazardous geophysical mass flows, such as snow avalanches, debris flows, lahars and pyroclastic flows, this has important implications for risk management in mountainous and volcanic regions. In this paper an avalanche model that incorporates frictional hysteresis, as well as depth-averaged viscous terms derived from the $\unicode[STIX]{x1D707}(I)$ -rheology, is used to quantitatively model self-channelisation and levee formation. The viscous terms are crucial for determining a smoothly varying steady-state velocity profile across the flowing channel, which has the important property that it does not exert any shear stresses at the levee–channel interfaces. For a fixed mass flux, the resulting boundary value problem for the velocity profile also uniquely determines the width and height of the channel, and the predictions are in very good agreement with existing experimental data for both spherical and angular particles. It is also shown that in the absence of viscous (second-order gradient) terms, the problem degenerates, to produce plug flow in the channel with two frictionless contact discontinuities at the levee–channel margins. Such solutions are not observed in experiments. Moreover, the steady-state inviscid problem lacks a thickness or width selection mechanism and consequently there is no unique solution. The viscous theory is therefore a significant step forward. Fully time-dependent numerical simulations to the viscous model are able to quantitatively capture the process in which the flow self-channelises and show how the levees are initially emplaced behind the flow head. Both experiments and numerical simulations show that the height and width of the channel are not necessarily fixed by these initial values, but respond to changes in the supplied mass flux, allowing narrowing and widening of the channel long after the initial front has passed by. In addition, below a critical mass flux the steady-state solutions become unstable and time-dependent numerical simulations are able to capture the transition to periodic erosion–deposition waves observed in experiments.

Dense granular flows can spontaneously self-channelise by forming a pair of parallel-sided static levees on either side of a central flowing channel. This process prevents lateral spreading and maintains the flow thickness, and hence mobility, enabling the grains to run out considerably further than a spreading flow on shallow slopes. Since levees commonly form in hazardous geophysical mass flows, such as snow avalanches, debris flows, lahars and pyroclastic flows, this has important implications for risk management in mountainous and volcanic regions. In this paper an avalanche model that incorporates frictional hysteresis, as well as depth-averaged viscous terms derived from the µ(I)-rheology, is used to quantitatively model self-channelisation and levee formation. The viscous terms are crucial for determining a smoothly varying steady-state velocity profile across the flowing channel, which has the important property that it does not exert any shear stresses at the levee-channel interfaces. For a fixed mass flux, the resulting boundary value problem for the velocity profile also uniquely determines the width and height of the channel, and the predictions are in very good agreement with existing experimental data for both spherical and angular particles. It is also shown that in the absence of viscous (second-order gradient) terms, the problem degenerates, to produce plug flow in the channel with two frictionless contact discontinuities at the levee-channel margins. Such solutions are not observed in experiments. Moreover, the steady-state inviscid problem lacks a thickness or width selection mechanism and consequently there is no unique solution. The viscous theory is therefore a significant step forward. Fully time-dependent numerical simulations to the viscous model are able to quantitatively capture the process in which the flow self-channelises and show how the levees are initially emplaced behind the flow head. Both experiments and numerical simulations show that the height and width of the channel are not necessarily fixed by these initial values, but respond to changes in the supplied mass flux, allowing narrowing and widening of the channel long after the initial front has passed by. In addition, below a critical mass flux the steady-state solutions become unstable and time-dependent numerical simulations are able to capture the transition to periodic erosion-deposition waves observed in experiments.

Introduction
Self-channelisation and levee formation can occur in a wide range of geophysical mass flows that take place in volcanic and mountainous regions throughout the world. These include highly mobile and destructive pyroclastic flows (Rowley, Kuntz & MacLeod 1981;Wilson & Head 1981;Branney & Kokelaar 1992;Calder, Sparks & Gardeweg 2000;Jessop et al. 2012), water-saturated lahars (Vallance 2000;Vallance & Iverson 2015) and debris flows (Sharp & Nobles 1953;Costa & Williams 1984;Pierson 1986;Iverson 1997;Major 1997) as well as wet snow avalanches (Jomelli & Bertran 2001;Ancey 2012;Bartelt et al. 2012;Schweizer, Bartelt & van Herwijnen 2014). Although these flows vary greatly in composition, a unifying feature is that, on shallow slopes, they spontaneously form static parallel-sided levees that bound a central flowing channel. The static levees prevent lateral spreading, allowing deeper flows to be sustained for longer than a spreading flow, and thereby increase the flow's mobility and its eventual run-out distance.
Although the mobility of grains is important for both industrial and geophysical applications, modelling the self-channelisation process still presents a significant theoretical challenge. This is firstly because it is a particularly subtle test of the granular rheology, since the flow spontaneously selects its own width, rather than being laterally unconfined or having the width imposed by side walls. Secondly, it also raises fundamental questions about how static and flowing regions can coexist, which is a longstanding issue in modelling granular materials. Although this paper is motivated by complex geophysical granular flows, it is focussed on determining the simplest possible formulation that will capture the levee formation process in quasi-monodisperse dry granular flows (Félix & Thomas 2004;Deboeuf et al. 2006;Takagi, McElwaine & Huppert 2011).
In order to narrow down the physical mechanisms required for levee formation, Félix & Thomas (2004) performed small-scale experiments with 300-400 µm spherical dry glass beads that were steadily released from a point source onto an inclined plane roughened with a glued layer of 425-600 µm glass beads, to ensure that there was no slip at the base. For illustration, a similar experiment using 160-200 µm red sand on a bed of turquoise ballotini (750-1000 µm) is shown in figures 1 and 2 as well as supplementary movie 1 (online) available at https://doi.org/10.1017/jfm.2019.518. Both of these sets of experiments are dry and have a very narrow range of particle sizes, but, as Félix & Thomas (2004) showed, static levees still form. This suggests that neither interstitial fluid nor particle-size segregation are essential to the self-channelisation process, although they may strongly enhance its effects (Pouliquen, Delour & Savage 1997;Pouliquen & Vallance 1999;Félix & Thomas 2004;Goujon, Dalloz-Dubrujeaud & Thomas 2007;Iverson et al. 2010;Woodhouse et al. 2012;Kokelaar et al. 2014;Baker, Johnson & Gray 2016b).
A self-channelised leveed flow has three distinct phases of motion (Félix & Thomas 2004) as illustrated in figure 2(a-c) respectively and supplementary movie 1. As grains flow down the incline they form a fully mobilised head at the front of the flow, which propagates at a constant speed and lies downslope of the levees. The downslope velocity is greatest in the surface layers in the centre of the channel (Jop, Forterre & Pouliquen 2005;Kokelaar et al. 2014;Baker, Barker & Gray 2016a). This region therefore transports grains towards the flow head, but when the grains reach the front they lose their confinement and spread out. Some of the grains are overrun by the flow itself, while others are advected to the sides and come to rest, extending the static levees , which advance downslope at the back of the head (figures 1 and 2a).

Levees
Flow head Channel FIGURE 1. An oblique view of a self-channelised leveed flow of 160-200 µm red sand on a plane roughened with a layer of 750-1000 µm turquoise spherical ballotini and inclined at ζ = 34 • ± 0.1 • . The granular material is released from a funnel at the top of the chute, and a self-channelised flow rapidly develops, which moves down the plane at constant speed. The grains spontaneously select a flowing channel width W = 3.60 ± 0.05 cm and the total width of the levees and the channel is W total = 5.00 ± 0.05 cm. A long shutter time is used to blur the moving grains and thus highlight the parallel-sided static levees where the grains are in sharp focus. The steady-state mass flux Q M = 8.6 ± 0.1 g s −1 is measured as the grains flow off the end of the chute. A movie showing the time-dependent evolution is available in the online supplementary material (movie 1).
The flow within the head is slightly deeper than that in the channel (Félix & Thomas 2004), so as the head passes by there is a small decrease in flow thickness and the levees stabilise, selecting a width for the channel (figure 2b). For sufficiently wide flows the central channel is of almost constant thickness (Félix & Thomas 2004), but in narrower channels the flow may have a pronounced cross-slope gradient in the free surface (Félix & Thomas 2004;Takagi et al. 2011), with the deepest flow in the centre of the channel. This gradient may be an indication that normal stress differences are an important component of the underlying granular rheology (McElwaine, Takagi & Huppert 2012). Once the parallel-sided levee-channel morphology has been established, it persists until the inflow ceases and the channel then partially drains to leave pronounced levee walls on either side of a thinner deposit within the channel (figure 2c). The thickness of these levees is similar to the thickness of the flow that generated them, especially for wide channels.
The width of the flowing channel is not necessarily set at the point where the levees are emplaced immediately behind the head. Figure 3(a-c) shows an overhead view of the flow front advancing from left to right down the slope and depositing levees behind the flow head. A subsequent internal surge in the flow (figure 3d-f ) then pushes material over the top of one of the levee walls and re-mobilises it. This allows the central channel to become slightly wider and a new levee is formed, further out, which maintains the self-channelisation of the flow. This process of levee-bank overtopping allows the channel to adjust to its final steady-state width after the initial passage of the flow head. For the angular red sand particles used in this paper, the final steady-state levee width is achieved within a few tens of seconds. Smaller fluctuations around the steady state constantly erode and deposit small amounts of material on the inside of the levee walls (see movie 2) without remobilising the entire levee. The levees are therefore in active equilibrium with the flow, constantly undergoing minor readjustments that ensure that the width of the channel reflects its current inflow mass flux.
Particle shape strongly influences the formation of self-channelised leveed flows. In flows of spherical glass ballotini, the levees are considerably less stable (resistant to erosion) than they are for sand, and the levees creep outwards over a time scale that is much longer than that of their initial emplacement. Deboeuf et al. (2006) found that for flows of 300-400 µm spherical glass beads on a rough plane made of 200 µm sandpaper, the channel widened and the flow thickness decreased very slowly until an asymptotic state was approached after very long times (>1 h). Similar experiments for even longer times (Takagi et al. 2011) indicated that the flow widened and eventually became intermittent, with waves further widening the levees, implying that a steadystate width was never attained.
In contrast, Takagi et al. (2011) showed that flows of 300-600 µm angular sand particles, on a bed made of the same sand, rapidly established a steady-state leveed channel for sufficiently high mass fluxes. The flow had a fixed channel width and thickness, as well as a well-defined downslope surface velocity profile across the channel. As the mass flux was increased, Takagi et al. (2011) found that the flow thickness stayed approximately constant, but the width of the channel increased. If the mass flux was decreased below a critical threshold, however, an unsteady regime developed in which regular pulses of material flowed down the static channel as a series of erosion-deposition waves (Daerr 2001;Börzsönyi, Halsey & Ecke 2005;Clément et al. 2007;Börzsönyi, Halsey & Ecke 2008;Takagi et al. 2011;Edwards & Gray 2015;Edwards et al. 2017).
A less well studied aspect of self-channelised leveed flow is the changing shape of the static levees as the inflow mass flux is altered (figure 4). For an inflow mass flux of Q M = 17.0 ± 0.1 g s −1 a stable central channel forms with a width W = 4.40 ± 0.05 cm (figure 4a). When the flux is reduced to Q M = 7.7 ± 0.1 g s −1 the channel narrows to a new steady-state width W = 2.70 ± 0.05 cm (figure 4b). The old levee is left in situ and, as the flowing region retreats into the centre of the channel, very small leveed structures are left on the inside of the old levee walls. For these small-scale experiments these features are barely one grain-diameter thick, but they can be still be seen in figure 4(b) due to the oblique lighting. An immediate consequence of this  (a-f ) A series of overhead images at t = 3, 5, 7, 9, 10 and 11 s showing (a-c) the initial formation of the static levees and (d-f ) their subsequent reworking by levee-bank overtopping. The mass flux Q M = 8.6 ± 0.1 g s −1 , ζ = 34 • ± 0.1 • and the final width of the flow (moving channel plus levees) is W total = 5.00 ± 0.05 cm. Note that at 9 s a wave begins to propagate down the channel and continuously overtops the levee bank on one side, forming a new levee that is slightly further out. This allows the flowing channel width to readjust to its steady-state value W = 3.60 ± 0.05 cm. The downslope flow direction is from left to right. A movie showing the time-dependent evolution is available in the online supplementary material (movie 1).
FIGURE 4. Sequence of overhead images of a self-channelised flow of 160-200 µm red sand on a rough plane inclined at ζ = 34 • ± 0.1 • . The aperture of the funnel supplying the grains to the top of the chute is changed in order to reduce the mass flux from (a) its initial value of Q M = 17.0 ± 0.1 g s −1 to (b) 7.7 ± 0.1 g s −1 and then (c) back up to 27.0 ± 0.1 g s −1 . As the mass flux is reduced, the flowing channel width narrows from (a) W = 4.40 ± 0.05 cm to (b) W = 2.70 ± 0.05 cm by mass accretion to the inside of the levee walls. The old levees are left in situ and record the fact that a higher flux once propagated down the channel. When (c) the mass flux is subsequently increased again the levee walls are pushed out by levee-bank overtopping and the old levees are fully remobilised to form a new channel of width W = 5.70 ± 0.05 cm. The downslope flow direction is from left to right. A movie showing the narrowing and widening of the channel is available in the online supplementary material (movie 3). observation is that the shape of the material outside of the central flowing channel is not unique and at least partially records the history of the flow. If, on the other hand, the inflow mass flux is increased to Q M = 27.0 ± 0.1 g s −1 a new stable width W = 5.70 ± 0.05 cm rapidly develops by levee-bank overtopping (figure 4c) and the history, preserved in the old stacked levees, is erased. The significance of this process is that it may allow information about the time history of natural geophysical flows to be inferred from the deposit, for example from sequences of subtly stacked levees that occur in pyroclastic deposits (Rowley et al. 1981;Wilson & Head 1981;Branney & Kokelaar 1992;Calder et al. 2000).
The fact that static levees are deposited on an inclined plane alongside flowing grains led Félix & Thomas (2004) to suggest that levee formation is related to frictional hysteresis (Daerr & Douady 1999;Pouliquen 1999;Pouliquen & Forterre 2002;Edwards & Gray 2015;Edwards et al. 2017). A simple example of this hysteresis is when a steady uniform flow is brought to rest on a slope of angle ζ , it leaves a deposit of thickness h = h stop (ζ ), but a layer of this thickness will not start to flow again until the inclination angle is increased to ζ = ζ start (h). Since the inverse function h start (ζ ) is greater than h stop (ζ ), there are a range of thicknesses over which static and flowing layers can coexist. The underlying cause of this behaviour is a non-monotonic relationship between the flow velocity and the basal friction coefficient, as expressed by the friction law of Pouliquen & Forterre (2002). This phenomenological law is applicable to spherical grains and combines a flow rule for steady uniform flows (Pouliquen 1999) with measurements of ζ start (h) to describe the friction of both static and flowing layers. Mangeney et al. (2007) used Pouliquen & Forterre's (2002) non-monotonic empirical friction law to model the self-channelisation process within the framework of classical depth-averaged avalanche equations (e.g. Grigorian, Eglit & Iakimov 1967;Savage & Hutter 1989;Gray, Wieland & Hutter 1999;Heinrich, Piatanesi & Hebert 2001;Gray, Tai & Noelle 2003;Mangeney-Castelnau et al. 2003;Pitman et al. 2003). Mangeney et al.'s (2007) numerical simulations exhibited both a central flowing channel and parallel-sided static/very-slowly creeping margins that were similar to those seen in experiments. The computations explicitly demonstrated how the material flowing down the central channel spread out laterally at the flow front, slowed down and then deposited to form static lateral margins behind the head. These simulations therefore implied that a heterogenous rheology, due to interstitial fluid and/or size segregation, was not essential for modelling levee formation. Despite the ground-breaking nature of this work, the computed flow thickness in the fully developed flowing channel was only very slightly greater than h stop , the minimum thickness for steady uniform flow in the friction law of Pouliquen & Forterre (2002). The channel was therefore significantly thinner and wider than those observed experimentally (Félix & Thomas 2004) as Mangeney et al. (2007) acknowledged themselves. Moreover, the downslope velocity in Mangeney et al.'s (2007) simulations was plug-like across the channel, which contradicted the experimental observation that the downslope velocity decreases smoothly to zero at the levee-channel interfaces across a shear band (Félix & Thomas 2004;Deboeuf et al. 2006;Takagi et al. 2011).
In this paper a depth-integrated model for self-channelised flows is developed, taking into account contributions of the in-plane deviatoric stress, which lead to depth-averaged viscous-like terms (Gray & Edwards 2014;Baker et al. 2016a). This model is referred to as the viscous depth-averaged model in this paper to distinguish it from classical inviscid shallow-water-like depth-averaged avalanche models, which do not have viscous second-order gradient terms in their depth-averaged momentum balances. A steady-state equation of motion is derived, which describes the balance of stresses that cause shear bands (e.g. Pouliquen & Gutfraind 1996;Schall & van Hecke 2010) to form, and shows that these stresses provides the vital missing physics that sets the steady-state thickness and width of the whole channel. This stress balance results from the interaction of two physical mechanisms, namely frictional hysteresis (Daerr & Douady 1999;Pouliquen & Forterre 2002;Mangeney et al. 2007;Edwards et al. 2017Edwards et al. , 2019 and lateral viscous stresses (Baker et al. 2016a). Steady-state theoretical predictions of the width, thickness and velocity profile across the channel (as functions of mass flux and slope angle) are in good quantitative agreement with laboratory experiments performed with spherical glass beads (Félix & Thomas 2004), and also sand (Takagi et al. 2011). Importantly, this paper clearly demonstrates that the inviscid theory is a singular limit of the viscous case, with degenerate non-unique steady longitudinally uniform states that have a range of channel thicknesses and widths for the same mass flux. In addition, two-dimensional time-dependent numerical solutions of this viscous avalanche model are able to explicitly compute how channels with the correct steady-state thickness and width are established dynamically, as well as allowing more complex unsteady flows to be investigated.

Governing equations
The depth-averaged theory of Edwards et al. (2017) for erosion-deposition waves is used here to model the formation of self-channelised leveed flows on non-erodible slopes.
2.1. Depth-averaged model Granular material is supplied with constant mass flux from a funnel onto a rough plane inclined at an angle ζ to the horizontal. A rectangular Cartesian coordinate system Oxyz is defined with the x-axis pointing downslope, the y-axis in the cross-slope direction and z the upward perpendicular to the plane (figure 5). The granular material is assumed to be incompressible with constant bulk density ρ and velocity components u = (u, v, w) in the downslope, cross-slope and normal directions, respectively. The depth-averaged velocity fieldū = (ū,v) in the downslope and cross-slope directions is then defined as where h is the avalanche thickness. Depth averaging the incompressibility condition and applying kinematic boundary conditions (Savage & Hutter 1989)   . A schematic diagram of the coordinate system Oxyz inclined at an angle ζ to the horizontal, where the x-axis points downslope, the y-axis is the cross-slope direction and the z-axis is the upward normal to the chute. For a steady uniform self-channelised flow, the cross-slope thickness profile is h(y), whilst the velocity field is given by u(y, z). The moving central channel has width W and is bounded on either side by parallel static levees. where g is the constant of gravitational acceleration and ν is a coefficient in the depthaveraged kinematic viscosity νh 1/2 /2, discussed in § 2.3. Implicit in the momentum transport terms of these equations is the assumption that the shape factor (u 2 /ū 2 ) is unity (Pouliquen 1999 which arise from the gravitational force that pulls the grains downslope and the effective basal friction. The direction of the frictional force is determined by the vector where i is the unit vector in the downslope direction and ∇ is the two-dimensional gradient operator. This ensures that the friction opposes the motion whenū is nonzero, and opposes the resultant force due to gravity and the depth-averaged pressure gradient when the material is stationary. Equations (2.2)-(2.4) were originally derived by assuming that the flow was shallow and depth averaging the mass and momentum balance equations, assuming the µ(I)rheology (GDR-MiDi 2004;Jop, Forterre & Pouliquen 2006) and a no-slip condition at the base (Gray & Edwards 2014;Baker et al. 2016a). To leading and first order in the aspect ratio, this process yielded the classical inviscid depth-averaged avalanche equations (given by setting ν = 0), with an effective basal friction law corresponding to that measured empirically by Pouliquen & Forterre (2002) for steady uniform flows (which is referred to as the dynamic frictional regime in § 2.2). The depth-averaged viscous terms (i.e. those terms multiplied by ν) emerge at second order, by retaining the in-plane normal and shear stresses and approximating them using the leading-order lithostatic pressure distribution and Bagnold velocity profile through the flow depth. This explicitly determines the coefficient ν in terms of known frictional parameters, see § 2.3, rather than it being a fitting parameter. The viscous terms are usually very small, but the fact that they are the highest gradient terms makes their inclusion a singular perturbation of the equations. As a result there are some problems where they play a crucial role. These include (i) obtaining the correct cutoff frequency of roll waves (Forterre 2006;Gray & Edwards 2014), (ii) generating cross-stream velocity profiles in channels (Baker et al. 2016a) and (iii) producing well-posed models of segregation-induced fingering (Baker et al. 2016b). It shall be shown in this paper that they are also play a vital role in the selection of the height, width and velocity profile across a monodisperse leveed channel.

The effective basal friction law
The effective coefficient of basal friction µ b encodes information about the µ(I)rheology (GDR-MiDi 2004;Jop et al. 2006) as well as the hysteretic behaviour of the granular material (Daerr & Douady 1999;Pouliquen & Forterre 2002;Mangeney et al. 2007;Edwards et al. 2017Edwards et al. , 2019. This hysteretic behaviour is described by a non-monotonic friction law (Pouliquen & Forterre 2002;Forterre & Pouliquen 2003;Edwards et al. 2017Edwards et al. , 2019   In Pouliquen & Forterre (2002) the dynamic friction µ D is based on the empirical flow rule for glass beads Fr = βh/h stop (Pouliquen 1999), where β is a constant obtained experimentally. By combining this with measurements of h stop , Pouliquen (1999) derived a dynamic friction law that was an increasing function of Fr/h and was valid for Fr β. Angular sand grains obey a more general flow rule (Forterre & Pouliquen 2003) where Γ = 0.84, and β = 0.71 is considerably higher than the value β = 0.15 for glass beads (Forterre & Pouliquen 2003). Note that these values have been corrected from the previously published values to account for the factor √ cos ζ in the Froude number (2.8). The flow rule (2.9) implies that a steady uniform flow at h = h stop has a Froude number Fr = β − Γ that is negative for the parameters for sand, a contradiction since Fr 0 from (2.8). This led Edwards et al. (2017) to suggest that (2.9) is valid only for Fr β * , where β * is both strictly positive and greater than β − Γ .
The flow rule (2.9) can be used to eliminate h stop in the reciprocal form of the empirical h stop curve (Pouliquen & Forterre 2002) to derive the dynamic friction law where Fr β * and the parameters µ 1 = tan ζ 1 , µ 2 = tan ζ 2 and L are fitted to measurements of the h stop (ζ ) curve (Pouliquen 1999;Pouliquen & Forterre 2002;Forterre & Pouliquen 2003). The angles ζ 1 and ζ 2 are the minimum and maximum angles, respectively, at which steady uniform flow is observed.
The friction force for static material (Fr = 0) is exactly that which is required to keep the material stationary, up to a maximum static friction coefficient. This coefficient therefore takes the form where the first argument to the min function is the friction required to balance the other forces acting on the static layer and the second argument is the maximum static friction. The maximum static friction is obtained by fitting measurements of h start (ζ ), and introduces a further friction-law parameter µ 3 = tan ζ 3 , where ζ 3 is the minimum angle at which an infinitely deep static layer will start to flow spontaneously (Daerr & Douady 1999;Pouliquen & Forterre 2002). The friction in the intermediate regime 0 < Fr < β * is a power-law interpolation between the maximum static friction and the minimum dynamic friction at Fr = β * and takes a much more complicated form for angular particles than for spherical ballotini (Pouliquen & Forterre 2002). The general form that accounts for spherical and irregular particles is given by Edwards et al. (2017) as where κ is a constant interpolation power. In the absence of experimental measurements for slow flows 0 < Fr < β * , Pouliquen & Forterre (2002) suggested that κ = 10 −3 . However, this value produces an extremely rapid decrease in the friction as the Froude number is increased from zero, so sharp that it cannot be resolved numerically using standard double-precision floating-point numbers . Edwards et al. (2017) instead suggest that κ needs to be at least 10 −1 to create a robust metastable state, in which a static layer, just slightly thinner than h start (i.e. within the hysteretic region), remains stationary unless it is perturbed. In this paper it is assumed for simplicity, as in Edwards et al. (2017Edwards et al. ( , 2019, that κ = 1, which is consistent with the experimental results of Russell et al. (2019) on retrogressive failures. The parameter β * sets the Froude number of the transition between the dynamic and intermediate regimes. It also defines the thickness h * of the minimum steady uniform flow, which lies between h stop (ζ ) and h start (ζ ). In Pouliquen & Forterre (2002) β * = β (and thus h * = h stop ), while Edwards et al. (2017) assumed that h * was half-way between h stop and h start . This paper follows Edwards et al. (2019) where Λ is a constant for all slope angles. Since steady uniform flows satisfy the empirical flow rule (2.9) it follows that the transition Froude number is constant (2.14) In general, Λ must be greater than unity and chosen so that h * < h start , to ensure that there is hysteresis in flows of thickness h ∈ [h * , h start ]. The new approach  has the advantage that it defines h * over the complete range of steady uniform flow angles ζ ∈ [ζ 1 , ζ 2 ], rather than just in the range [ζ 3 , ζ 2 ], as in Edwards et al. (2017). It also allows the intermediate friction to be accessed even for materials, such as sand, for which Γ > β.
It is important to note that the extended law proposed by Edwards et al. To enable quantitative comparison between theory and experiment, the material parameters used throughout the paper (table 1) correspond as closely as possible to the sand used by Takagi et al. (2011) and the glass beads used by Félix & Thomas (2004). For the experiments with sand, the angles in the friction law, ζ 1 and ζ 2 , are obtained by fitting h stop (ζ ) to the experimental data of Takagi et al. (2011) using the functional form proposed by Pouliquen & Forterre (2002) (figure 6a), which leads to friction coefficients in the form of (2.10)-(2.12). The characteristic length scale L and the parameters β and Γ are taken from measurements of Forterre & Pouliquen (2003) for a similar flow of sand on a rough bed of the same material. The angle ζ 3 is determined from measurements of h start (ζ ), but these measurements are not reported by Takagi (Takagi et al. 2011, blue markers) and glass beads (Félix & Thomas 2004, red markers), with best-fit h stop curves using the friction-law parameters calibrated from these experiments (table 1). (b) The two friction laws for sand (blue curve) and glass beads (red curve) as a function of the Froude number Fr at fixed thickness and inclination angle. Parameters h = 8 mm and ζ = 32 • for sand and h = 3.2 mm and ζ = 25 • for glass beads are chosen to match those in the respective experiments.  (2004), but with a smaller L to quantitatively fit the h stop measurements of Félix & Thomas (2004). Note that the values for β and Γ differ from those in Pouliquen & Forterre (2002) and Forterre & Pouliquen (2003) to account for the factor √ cos ζ in the definition of the Froude number (2.8).
measurements of similar flows of sand over a rough bed. For the experiments with glass beads performed by Félix & Thomas (2004), the parameters are chosen to be the same as the ones used by Mangeney et al. (2007) to simulate these flows, although here a value of L smaller than that of Mangeney et al. (2007) is used to quantitatively match the experimental h stop curve (figure 6a). The resultant friction laws are plotted in figure 6(b) as a function of the Froude number. The frictional hysteresis of the angular sand particles (µ(Fr = 0) − µ(Fr = β )) is much greater than that of glass beads, which makes the sand levees significantly stronger.

Depth-averaged kinematic viscosity
The second-order depth-averaged viscous-like terms in the momentum balances (2.3)-(2.4) contain a parameter ν in the depth-averaged kinematic viscosity νh 1/2 /2, for which Gray & Edwards (2014) This is a function of the slope angle ζ and parameters that are already known from the effective basal friction law, so no new parameters are introduced into the theory. It is well-defined provided ζ ∈ [ζ 1 , ζ 2 ], i.e. in the range of angles where steady uniform flows develop. Outside this range, some form of regularisation is required to ensure that it does not become negative. In this paper, the two-dimensional viscous terms derived by Baker et al. (2016a) are used for flows in all three frictional regimes, even though their original derivation (Gray & Edwards 2014) implicitly assumed the flow was in the dynamic frictional regime.

Fully developed self-channelised flow
This section presents exact steady-state solutions for the height, width and downslope velocity profile across the central flowing channel. The depth-averaged granular viscosity provides the vital mechanism for producing a smoothly varying velocity profile across the channel, allowing appropriate boundary conditions to be imposed at the levee-channel interfaces. These boundary conditions together with an integral mass flux constraint are then able to determine a unique equilibrium channel thickness and width, completely independently of the flow head dynamics. In the absence of viscosity the equations are not closed, so there is no unique solution and inviscid theories therefore lack a crucial thickness and width selection mechanism.
Outside of this region, in the levees, there are multiple static states (as demonstrated experimentally in figure 4) because the static friction encoded in (2.11) can take any value between zero and its maximum value. While there is not a unique solution for the width of the levees, the minimum levee width necessary to support a central flowing channel of a given thickness, can be determined by assuming all the static grains are at the maximum static friction.
3.1. Steady-state depth-averaged equations in the flowing channel Consider a fully developed steady-state parallel-sided leveed flow that has a central flowing channel of width W and is supplied by a mass flux Q M . It follows that the flow is independent of time t and the downstream coordinate x. In the central channel the depth-averaged mass balance (2.2) then reduces to which is subject to a condition of no flow across the levee-channel interfaces, v = 0 at y = ±W/2. This can be integrated to show that the flow thickness is constant where the constant H is, as yet, unknown. The downslope component of the depthaveraged momentum balance (2.3) then reduces to whereū is assumed to be strictly positive. This is a second order autonomous ordinary differential equation ( and no lateral shear stress at both channel-levee interfaces. In the fully developed steady state, the mass flux of grains entering the chute is equal to the mass flux flowing down the central channel at any downstream location. Hence,ū(y) is subject to the integral constraint The thickness H, width W and the depth-averaged downslope velocityū(y) across the flowing channel, are direct predictions of the model (3.5)-(3.8). To visualise the solution for the entire flow including the levees, two additional assumptions can be made, which are described in appendices A and B. The first assumption is of a constant velocity profile through the depth of the flow, which allows the downslope velocity u(y, z) to be reconstructed from the depth averaged downslope velocityū(y). The second assumption is that the static levees are everywhere on the brink of yield (Hulme 1974;Balmforth, Burbidge & Craster 2001), which allows the thickness profile of the static levee to be computed. These assumptions are not intrinsic to the model, but extend its predictions to allow comparison with experimental measurements of the surface velocity and the combined width of the flow and levees, respectively.

Inviscid solutions
In the inviscid case, the coefficient ν = 0 and the steady-state equation of motion (3.5) reduces from a second-order ODE to an algebraic balance between friction and gravity only, µ b = tan ζ . This reduction in order is indicative of a singular perturbation problem in which, as is the case here, the inviscid ν = 0 model differs qualitatively from the model with any non-zero viscosity ν > 0. Since the thickness H is constant across the channel, solving µ b (H,ū) = tan ζ with µ b in the dynamic regime (3.4) gives a constant velocityū across the channel, where (3.10) and the thickness H is still a free parameter. The reduction of order in the inviscid system means that both tangential velocity boundary conditions (3.6) and (3.7) must be relaxed. Instead, since the material in the central channel moves at constant speed u steady and the grains in the levee are static, there will be contact discontinuities at the levee-channel boundaries y = ±W/2 (i.e. jumps in the downslope velocity parallel to the margins, see e.g. Chadwick 1976) as shown in figure 7. This is not what is observed in experiment (Félix & Thomas 2004;Deboeuf et al. 2006;Takagi et al. 2011), where there are smoothly varying velocity profiles across the channel, rather than frictionless surfaces sliding past one another.
Assuming that there is a single channel, the integral constraint (3.8) together with the velocity field (3.9) provide one equation expressing the flowing channel width W as a function of the thickness H, .
( 3.11) With no further boundary conditions imposed in the inviscid model, there is no equation to determine the thickness H. The only additional constraint on the inviscid solution is that the depth-averaged pressure across the contact discontinuities must be equal (Chadwick 1976), i.e. the depth-averaged pressure in the levee is just sufficient to balance the depth-averaged pressure in the central channel at y = ±W/2. This is equivalent to the condition that the thickness must be continuous across the contact discontinuities. Since static and moving grains must coexist for the same thickness, it follows that H must lie in the metastable range h * H h start . These inequalities are not sufficient, however, to determine H uniquely, and hence (3.9) and (3.11) do not determineū(H) and W(H) uniquely either. For a given mass flux Q M , the inviscid avalanche model has an infinite set of solutions that are parameterised by the H ∈ [h * , h start ]. Figure 7 shows three cases for Q M = 12 g s −1 and a slope angle ζ = 25 • using the parameters for glass beads from table 1. The widest and slowest moving central channel (W = 9.42 cm,ū steady = 2.93 cm s −1 ) occurs when H = h * = 2.9 mm, while the narrowest and fastest moving flow (W = 5.73 cm,ū steady = 3.94 cm s −1 ) is obtained when h = h start = 3.54 mm. It follows, that even a relatively modest range of thicknesses (H ∈ [2.9, 3.54] mm) can lead to a substantial range of channel widths. In this case the widest flowing channel is 1.64 times wider than the narrowest. Interestingly, the total width, including the minimal static levees on either side of the central flowing channel, is almost the same for both cases (i.e. W total (h * ) = 12.61 cm, W total (h start ) = 12.14 cm) as shown in figure 7(a,c). The relative insensitivity of the total channel width arises because the static levees need to be much wider to support the central flowing channel as its thickness approaches h start . This is due to the fact that dh/dy → 0 in (B 4) as h → h start , i.e. the levee becomes less stable (in the sense that it finds it harder to support thickness gradients) as the thickness approaches h start . Note that the total width of the channel is a non-monotonic function of H, since W total = 11.27 cm when H = 0.5(h * + h start ) as shown in figure 7(b).
A prime example of a physical system that has an infinite number of steady states are the static levees. As shown in figure 4, these can have arbitrary free-surface shapes that are dependent on the history of the flow that emplaced them. Having an infinite number of steady states for the inviscid flowing channel is not, however, physically realistic, because there is strong experimental evidence (Félix & Thomas 2004;Deboeuf et al. 2006;Takagi et al. 2011) that the thickness H and width W are functions of the imposed mass flux Q M . Since the model is not able to uniquely determine the central channel thickness H, and therefore (3.9) and (3.11) cannot be used to calculate W(H) andū(H) uniquely either, the inviscid theory is missing an important physical mechanism for selecting the channel thickness, width and velocity profile.

Viscous solutions for the central flowing channel
A solution to the steady-state equation of motion (3.5) will now be constructed that includes the effects of lateral viscous stresses. The system to be solved comprises of the second-order ODE (3.5) with two unknown parameters H and W. Four boundary conditions are therefore required to close the system. These are provided by the boundary conditions at the levee-channel interface (3.6) and (3.7) (which, due to the symmetry of the solutions as y → −y, provide three independent conditions) and by the integral mass flux condition (3.8). The boundary value problem, associated boundary conditions and the flux constraint therefore provide a closed system of equations, which determines the thickness, width, and velocity profile of a self-channelised flow.
Much of the solution can be obtained algebraically if the problem is written in terms of the flow thickness H, and the mass flux Q M is obtained as a function of H using (3.8). This relationship Q M (H) can then be inverted to numerically solve for the thickness at a specified mass flux. The symmetry implies that it is sufficient to solve the problem in the half-domain y ∈ [−W/2, 0] and it is therefore convenient to define a new coordinate system that is centred at the levee-channel interfacê (3.12) Immediately adjacent to the levees there is a slowly moving region that lies in the intermediate frictional regime. Recalling that κ = 1, the transformed ODE (3.5) with the intermediate friction law (2.12) takes a particularly simple form, where the coefficients are thickness dependent, Solving (3.13) subject to the boundary conditions (3.6) and (3.7) it follows that in the slowly moving layer adjacent to the levee. If the Froude number Fr < β * everywhere in the channel, the solution (3.16) is valid everywhere and the channel width is equal to one period of the cosine function, i.e. W = 2π/ √ a. These solutions are, however, not the ones that are observed physically.
Instead, the velocity increases until Fr = β * and then the flow transitions across to the dynamic frictional regime. From the definition of the friction law (2.7) and the Froude number (2.8), it follows that this occurs whenū is equal to a transition velocitȳ u transition = β * gH cos ζ .
( 3.17) From (3.16) this transition occurs when and, from (3.16)-(3.18), the velocity gradient at this point is dū dŷ whereū steady is the steady uniform flow velocity defined in (3.9), Since H is constant,ū steady , c and d are constant, and hence (3.20) is an autonomous second-order ODE that can be solved by making the substitution p = dū/dŷ and using the chain rule to give p dp dū =ū −ū steady cū + d .
( 3.22) Integrating (3.22) subject to the boundary conditions thatū =ū centre and p = dū/dŷ = 0 on the symmetry line, and taking the positive root, implies that  Valid solutions are only found for certain values ofū centre , and since the centreline velocity is a function of the thickness, there should be a range of possible values for H. In fact, it is observed that in order to obtain a self-channelised solution the thickness is limited to a finite interval, i.e. H ∈ [H min , H max ]. The first boundary can be found by assuming that the central steady flow is at the transition between dynamic and intermediate friction, i.e. that the centreline velocityū centre =ū transition . In this case, the right-hand side of (3.24) is zero, so dū/dŷ|ŷ transition = 0. Using (3.19) it follows thatū transition (H) = 2b(H)/a(H), which can be solved to determine an upper bound for H H max (ζ ). The second limit, H min , is found by noticing thatū centre is a decreasing function of H, and that there is a value of the thickness, for which the self-channelised orbit in the phase space becomes homoclinic. For this value of H trajectories starting at the origin of the phase plane go directly to the saddle point given by the equilibrium point between dynamic friction and gravity (dashed green line in figure 8a). Hence, the lower boundary H min (ζ ) is determined by settinḡ u centre =ū steady in (3.24). For the friction-law parameters and slope angle calibrated to the experiments of Takagi et al. (2011) (table 1), H min ≈ 7.9 mm and H max ≈ 9.2 mm, which are significantly deeper than both h stop ≈ 5 mm and h * ≈ 6.7 mm. On the other hand, for the glass beads parameters (Félix & Thomas 2004;Mangeney et al. 2007) (table 1) the range is much narrower, H min ≈ 3.152 mm and H max ≈ 3.166 mm, but, nonetheless, is still deeper than h stop ≈ 2.16 mm and h * ≈ 2.9 mm.
By choosing H ∈ [H min , H max ] a solution to (3.24) for the centreline velocityū centre is guaranteed to exist and can be found by numerical root finding techniques. Once this is given, the separable ODE (3.23) can then be solved by quadrature, i.e.
By construction the velocity gradient in (3.23) is zero whenū =ū centre , so it follows that the integrand in (3.25) is singular asū →ū centre . However, asū centre →ū steady there is a large central region of the flow where the downslope velocity profile is flat and the integral (3.25) is difficult to evaluate. It is therefore useful to linearise the righthand side of (3.23) aboutū =ū steady and then solve to obtain the approximate solution u =ū steady − (ū steady −ū centre ) cosh ŷ − W/2 cū steady + d .
(3.26) By integrating (3.25) toū =ū centre (1 − ε), where ε 1, and then using the approximation (3.26) it is possible to accurately determine the solution close toū steady and hence the half-width of the channelŷ centre = W/2. The solutions for the velocity profile and channel width W are shown in figure 8(b) for a range of flow thicknesses H. This relationship is inverted numerically to find the solution corresponding to a given mass flux Q M . For a fixed slope angle, flows with larger mass flux are wider and faster moving (but, somewhat counter-intuitively, are thinner) than those with a smaller mass flux.
3.4. Shear-band structure adjacent to the levees The physical mechanism that prevents erosion of the levees is the granular viscosity, since it allows steady-state solutions to be constructed that do not exert any shear stress on the inside of the levee walls, even though material is flowing down the central channel. It is therefore of interest to understand the boundary layer structure that allows this.
For As the mass flux is increased, the centreline velocity approaches the steady uniform flow solutionū steady and the velocity profile flattens, forming a central region of material of uniform thickness and approximately uniform depth-averaged velocity.
As H → H min the orbit becomes homoclinic, directly connecting the origin and the saddle point, and the period of the orbit, channel width W and mass flux Q M all diverge. The centreline velocityū centre tends to the steady uniform flow velocityū steady , given by (3.9). Equation (3.9) therefore defines an upper bound for the velocity in a self-channelised flow, as observed experimentally (Takagi et al. 2011). Close to this limit, an increase in mass flux primarily increases the channel width, with little effect on the flow thickness or speed. At constant mass flux, an increase in slope angle decreases the flow thickness (since H min (ζ ) is decreasing in ζ ), which, counter-intuitively, decreases the flow velocity.
Away from the saddle point the orbits are almost identical, so the boundary layers at the sides of the flow collapse on top of one another when plotted in the leveecentred coordinateŷ (figure 9c). These boundary layers have a shear-band structure (Pouliquen & Gutfraind 1996;Fenistein & van Hecke 2003;Schall & van Hecke 2010) that seamlessly connects the static material in the levee with the central region. The boundary layer has two length scales associated with it: an inner boundary layer where the intermediate friction is active, and an outer layer that is in the dynamic friction regime. The exact intermediate solution (3.16) implies that the inner boundary layer is of width while the approximate solution (3.26) suggests that the outer boundary layer has a typical width W outer = cū steady + d. (3.28) For the parameters given in table 1 for sand, the inner boundary layer width W inner ≈ y transition (H min ) = 7.16 mm, whilst the outer dynamic boundary layer is approximately 8.52 mm wide. These are shown by the pink and light green shaded regions in figure 9(c) and provide a good order of magnitude estimate for the boundary layer width. The governing ODE (3.5) describes a balance between its three terms, the forces due to the slope-tangential component of gravity, basal friction and lateral viscous stress. Hysteresis in the basal friction results in two distinct force balances, depending on the local flow velocity. In figure 9(d) forces are plotted as a function of the levee-centred cross slope directionŷ. The magenta solid line represents the resultant force of gravity and friction, whilst the blue curve shows how viscous forces vary across the channel for Q M = 150 g s −1 . All the forces are non-dimensionalised. Right next to the levees, where the flow is slow, the basal friction is greater than the downslope component of gravity, leading to an upslope resistive force. In this region, viscous stresses act in the downslope direction, balancing the friction and sustaining motion. As the velocity increases towards the centre, friction decreases and eventually is balanced by gravity at the intermediate equilibrium (the green circle in figure 9c,d).
Since this is an inflection point d 2ū /dy 2 = 0 and viscous forces also vanish. In the faster flow away from the levee wall, d 2ū /dy 2 < 0 and the viscous stress instead acts in the upslope direction, slowing the flow. In the absence of viscous stresses the flow in this region would accelerate to the steady uniform flow velocityū steady (3.9). In the central region, where the profile flattens (ŷ > 5 cm), viscous contributions are extremely small and the flow is basically governed by the balance between gravity and friction explaining why the velocity approachesū steady . Therefore, lateral viscous stresses provide the mechanism that connects the flow in these two regions and allow the interface between static and flowing granular layers to be modelled.

Comparison with experiments
The steady-state viscous theory is now compared with the experimental results of Takagi et al. (2011), for sand, and Félix & Thomas (2004) for glass beads. The width W and thickness H of the central flowing channel are direct results from the depth-averaged theory and provide the strongest comparisons with the experiments. The theory also solves for the depth-averaged velocity profile across the channel u(y), but this is not directly comparable with the surface velocities u s (y) measured in Takagi et al.'s (2011) experiments, because there is shear through the depth of the avalanche. The additional assumptions necessary to make a comparison are summarised in appendix A. Similarly Félix & Thomas (2004) report the total width of the flow W total , rather than the width of the flowing channel W, so a comparison requires additional assumptions about the minimal width of the levees as detailed in appendix B.

Takagi et al.'s (2011) experiments with sand
The predicted flow thickness H is plotted in figure 10(a) as a function of the mass flux Q M for the slope angle 32 • used by Takagi et al. (2011). The solution consists of two regimes. At low mass fluxes, only the intermediate friction is active and the solution is given by (3.16), while above a critical threshold (Q min = 5.045 g s −1 ) the solutions also have a region of dynamic friction. These will be referred to as the intermediate and dynamic regimes. In the intermediate regime, the thickness H ∈ [H max , h start ] is a rapidly decreasing function of the mass flux Q M and is not physically realised. In the dynamic regime, the thickness H ∈ [H min , H max ] ≈ [7.9, 9.2] mm is almost constant and asymptotes to H min as Q M → ∞. This is significantly above h stop = 5 mm and h * = 6.7 mm. These observations are consistent with the experiments of Takagi et al. (2011) with 300-600 µm angular sand particles on a 32 • slope who found that maximum flow thickness stayed approximately constant at a value of H = 8.3 ± 0.4 mm as the mass flux was increased as shown in figure 10(a). Notably, if κ = 10 −3 is chosen, as in Pouliquen & Forterre (2002), instead of κ = 1, both H min and H max are only slightly (∼0.1 %) greater than h * = 6.7 mm, in contradiction with the experimental measurements.  where the velocitiesū steady andū centre are given by (3.9) and by solving (3.24), respectively. The approximate solution is shown by the magenta line in figure 10(b) and lies very close to the experimental data, as well as the physical solution. The approximation (3.29) therefore provides a useful formula for the channel width. The maximum depth-averaged velocity increases initially with the mass flux and then saturates above Q M = 50 g s −1 . This trend is broadly in line with what Takagi et al. (2011) observed for the maximal surface velocity (figure 11a). The near-constant ratio of the surface velocity to the depth-averaged velocity provides information about shear through the depth of the avalanche. The best fit for the data implies u s /ū = 2.35, which lies significantly above the ratio u s /ū = 1.67 for Bagnold flow. This is consistent with the discrete element method (DEM) simulations of Silbert, Landry & Grest (2003) and non-local theory of Kamrin & Henann (2015), which predict a transition from Bagnold flow to exponential-like velocity profiles with depth, as the flow gets close to h * , even in the complete absence of side walls. The ratio of the surface velocity to depth-averaged velocity in figure 11 strongly motivates the use of a concave exponential velocity profile with u s /ū = 2.35 to reconstruct the non-depthaveraged velocity in this paper, which is discussed in more detail in appendix A.  Figure 11(b) shows a comparison between the predicted surface velocity across the channel, assuming u s = 2.35ū, and that derived from particle image velocimetry (PIV) measurements of Takagi et al. (2011) for a flow of sand on a 32 • slope 2 m down from the source. The steady-state theory predicts a slightly narrower channel with a slightly faster maximum surface velocity than that observed in the experiments. This may be an indication that the depth-averaged viscosity νh 1/2 is slightly larger than that assumed in the theory, which is derived for thicker flows that have well-developed Bagnold profiles. Takagi et al. (2011) observed that there was a minimum mass flux below which the steadily flowing leveed channel was replaced by an unsteady sequence of avalanches triggered at equal time intervals as the grains piled up near the source and failed periodically. An illustrative movie showing the equivalent experiment with the red sand used in this paper is available online (movie 4). This will be investigated further in § 4.4.

Félix & Thomas's (2004) experiments with glass beads
The solution for glass beads has the same structure as that for sand, consisting of an unphysical intermediate regime for low fluxes and a physically realised dynamic regime above Q min = 3.6 g s −1 . The small range of H ∈ [H min , H max ] ≈ [3.152, 3.166] mm implies that the predicted channel thickness is almost independent of the mass flux Q M and is also significantly larger than both h stop ≈ 2.16 mm and h * ≈ 2.9 mm. The agreement between H and the experimental thickness data of Félix & Thomas (2004) is very good, as shown in figure 12(a). As with sand, the channel thickness H is in agreement with experiments when κ = 1, but if κ = 10 −3 is chosen the thickness is under-predicted, with both H min and H max extremely close to h * ≈ 2.9 mm. Félix & Thomas (2004) also reported the total width of the flow as a function of the mass flux. This makes direct comparison with experiments more difficult, not least because frictional hysteresis implies that the levee widths are not unique. However, by assuming that all points in the levee are on the brink of yield (Hulme 1974;Balmforth et al. 2001) a unique solution for the minimal levee width can be obtained, as described in detail in appendix B. Adding two minimal levees widths to the flowing channel width W then gives a good approximation for the minimum total width W total of the channel. Figure 12(b) shows a comparison of both W total and W to the total width measured by Félix & Thomas (2004) as a function of the mass flux Q M . For relatively low mass fluxes most of the experimental data lie close to the minimum total width W total predicted by the solutions in which the dynamic friction is active. However, for larger values of Q M the width appears to be over predicted by the steady-state theory. Félix & Thomas (2004) do not report the time or position at which their steady-state data were collected, but their chute was relatively short (2 m) and they only collected data for 80 s. The most likely explanation for the apparent discrepancy at high mass fluxes is therefore not that the steady-state predictions are wrong, but that the experimental data were not collected far enough downstream or after long enough times. For instance, Deboeuf et al. (2006) performed experiments at Q M = 25 g s −1 on a 3 m chute and found that although the thickness was close to steady state within 100 s the width of the channel adjusted on a much longer time scale and was still slowly widening at 3500 s (see their figure 2). Even longer experiments by Takagi et al. (2011) found that the levee margins for glass beads eventually became unstable after 70-90 min.  Félix & Thomas (2004). The theoretical levee profiles assume a perfect balance between gravity, pressure gradient and maximum static friction as described in appendix B.
F. M. Rocha, C. G. Johnson and J. M. N. T. Gray Deboeuf et al.'s (2006) and Takagi et al.'s (2011) observations suggests that for weakly hysteretic materials, such as glass beads, it can take a long time and/or distance for the width of the flow to reach steady state and that steady state may itself be unstable. This is borne out by a comparison between the experiments of Félix & Thomas (2004) and the full steady-state thickness profiles (including the levees and the central channel) in figure 12(c). All of the experimental flows have thicknesses that are close to their steady-state value and the widths also agree for low mass fluxes Q M = 3.7 and 8.0 g s −1 . However, as the mass flux is increased to Q M = 17.5 and 38 g s −1 the width of the channel is seen to be far from steady state. As stated above, this apparent lack of agreement at high mass fluxes is most likely due to the experiments having not been run for long enough and/or the measurement position not being far enough downstream.
3.6. Reconstruction of the smoothly varying velocity field Figure 13 shows a reconstruction of the steady-state downslope velocity for glass beads in a cross-section through the channel for Q M = 3.7, 8.0, 17.5, 38 g s −1 . These solutions assume an exponential velocity profile through the avalanche depth, which is appropriate for thin flows close to h * (see appendix A). It is striking that the velocity varies smoothly across the whole channel, with a maximum at the surface and centre of the flow, as well as boundary layers adjacent to the levee-channel interfaces providing a seamless connection to the static levees. This is in stark contrast to the inviscid model, which has a uniform velocity profile across the channel with contact discontinuities (velocity jumps) at the levee-channel margins as shown in figure 7. These plots also highlight that for a given mass flux Q M and inclination angle ζ , the viscous model has a unique solution for the height, width and velocity, rather than an infinite number of steady states, parameterised by H ∈ [h * , h start ], in the inviscid case. Moreover the height of the flows is approximately the same and the minimal levees are therefore closely similar. It should be noted that in the bottom left and right corners of the central channel the flow is moving very slowly. This is consistent with Deboeuf et al.'s (2006) and Kokelaar et al.'s (2014) colour change observations, which demonstrated that the material in these regions was either static or very slowly moving. Such linings of the channel do not occur in the inviscid model as shown in figure 7.

Time-dependent numerical simulations
To investigate the evolution towards the steady state, fully time and spatially dependent numerical solutions to the system of conservation laws (2.2)-(2.4) are now computed.

Numerical method
Numerical solutions are calculated using a high-resolution semi-discrete nonoscillatory central (NOC) scheme for convection-diffusion equations (Kurganov & Tadmor 2000). This method has proved its ability to solve similar systems of conservation laws for erosion-deposition waves (Edwards & Gray 2015;Edwards et al. 2017), segregation-induced finger formation (Baker et al. 2016b) and bi-disperse roll waves (Viroulet et al. 2018). The equations are solved in conservative form, Self-channelisation and levee formation in monodisperse granular flows  (4.2a,b) and diffusive fluxes where an additional volume source term S inflow (x, y) has been added to the right-hand side of the depth-averaged mass balance equation (2.2). This provides a simple way of modelling the experimental supply of grains onto the chute from a funnel (see § 1).
Grains are added to the domain in a small circular region, described by whereS is a normalisation constant, R = 2.5 cm is the radius of the circular inflow and x 0 = 15 cm is the downstream position of its centre. Defining polar coordinates with an origin at the centre of the circle by x = x 0 + r cos θ and y = r sin θ, the total volume flux of material entering the chute is Since the total mass flux is Q M = ρQ, (4.8) the normalisation factor in (4.6) isS = 3Q M πρR 6 . (4.9) The numerical solutions start with initial conditions of an empty chute, h = 0, m = 0 at t = 0. An outflow boundary condition is implemented at the bottom of the chute by extrapolating the solution to a row of ghost cells beyond the computational domain (LeVeque 2002). The flow does not reach the top (x = 0) or sides (y = ±L y /2) of the domain and so a boundary condition h = 0 is trivially satisfied here.
The two-dimensional numerical domain (L x , L y ) = (2.5, 0.3) m is discretised into a rectangular grid with (n x , n y ) = (1000, 300) grid points, respectively. Numerical fluxes are computed by a generalised minmod limiter with θ = 2 (Kurganov & Tadmor 2000) and the spatially discretised equations are integrated in time using a secondorder Runge-Kutta method, with time steps determined by a CFL (Courant-Friedrichs-Lewy) number of 0.225, and a maximum step size t = 10 −4 s (LeVeque 2002), which is required to minimise creep of the levees. Numerical errors caused by the degeneracy in the conservative form of the equations when h = 0 are mitigated by introducing a minimum thickness threshold (set to 10 −5 m or approximately one tenth of a grain diameter), below which the thickness and momentum are set to zero.

Formation and partial drainage of a self-channelised flow
The downslope velocityū and thickness h of a typical numerical solution are shown in figures 14(a) and 15(a) respectively (supplementary movies 5 and 6). For the parameters of sand at a slope angle ζ = 32 • , the mass flux Q M = 85 g s −1 is large enough for a self-channelised avalanche to form. After an initial adjustment from the source condition in the uppermost ∼1 m, the simulated flow spontaneously selects its own steady-state channel thickness, width and cross-slope velocity profile. A flow head is rapidly established and propagates down the plane at constant speed (figure 14a), laying down a steady-state levee-channel morphology behind it once some initial transients have dissipated. Figure 14(d) shows the downslope velocity at x = 1.75 m. At t = 32 s the whole width of the flow is mobile at x = 1.75 m, but, as the front passes by, the sides of the flow rapidly solidify and by t 60 s the velocity profile across the channel is rapidly tending to the steady-state solution computed directly from (3.5). The velocity therefore is almost constant in the centre of the channel and smoothly connects to the levee-channel interfaces on either side across two boundary layers. The thickness in this flowing section is constant, as shown in figure 15(b) and by the thickness profile lines across the channel at x = 0.5, 1, 1.5 and 2 m in figures 14 and 15.
By t = 60 s (figures 14b and 15b) the velocity profile and constant channel thickness H = 7.9 mm are in good agreement with the steady-state predictions in § 3.3 (figures 14d and 15d). Once established the width of the flowing channel remains constant in the numerical simulations. The outer edges of the levees have some residual creep, which is less than 1 % of the flow velocity in the channel, for the time step of 10 −4 s used here. Importantly this creep velocity scales with the time step and so the levees converge to a static state under time step refinement.
The time scale for establishment of the steady state is comparable to that in the experiments with 160-200 µm red sand (see supplementary movie 1), but is much quicker than Takagi et al.'s (2011) experiments which took up to 20 min. This discrepancy may be due to a difference in bed roughness: whereas Takagi et al. and a deposition wave moves rapidly downslope, thinning and then depositing material that was previously flowing in the central channel (figures 14c and 15c). This process is illustrated by cross-slope thickness profiles at x = 1.75 for t = 100, 105.3, 110 s in figure 15(d). At t = 100 s the channel thickness is still equal to the theoretically predicted steady-state value H = 7.9 mm and the surface gradient of the levee walls is below the maximum static friction, so they are stable and support the central channel.
At t = 105.3 s, in the middle of the deposition wave, the flow wanes and only a region near the centre of the channel remains flowing. This deposition wave leaves behind a static deposit with super-elevated levee walls that have a maximum thickness close to the original flow thickness H and a deposit in the central channel that lies just above h stop (figure 15d, t = 110 s). As noted by Félix & Thomas (2004) and Mangeney et al. (2007), the flow thickness H and channel width W are preserved in the deposit, allowing field observations to be used to estimate the inflow mass flux required to sustain the channel during flow. This provides an important constraint on the time evolution of the geophysical processes that created the levees, even when the flow event was not observed directly.

Narrowing and widening of the central flowing channel
The response of the steady fully developed solution to changes in the inflow mass flux is now investigated. The initial state is chosen to be the same as in figures 14(b) and 15(b) for an inflow mass flux Q M = 85 g s −1 and slope angle ζ = 32 • . The mass flux is then reduced to Q M = 70 g s −1 at t = 100 s as shown in movies 7 and 8. The channel does not narrow uniformly as shown in figure 16(a,b). An expanding wave travels rapidly downslope, which separates a downstream region, where the channel is still at the steady state for Q M = 85 g s −1 , and an upstream region where the flow attains the new steady state for Q M = 70 g s −1 . The front section of this wave travels slightly faster than the rear, so the flowing regions pinch off and separate from one another. As a result there is a brief period where the grains come to rest in the channel, before they are remobilised by the second half of the wave as shown in figure 16(a,b). This behaviour is also observed in the experiments with red sand as shown in movie 3. The wave eventually propagates out of the domain to leave a parallel-sided leveed channel with a width W = 10.3 cm that is narrower than the original width W = 12.1 cm for Q M = 85 g s −1 . The thickness of the channel before and after the change in flux are H before = 7.9145 mm and H after = 7.904 mm, respectively, which are virtually the same as shown in figure 16(e). Once the channel has established itself the computed downslope velocities at x = 1.75 m are almost identical to the steady-state theory in § 3.3 as shown in figure 16( f ). During the process of channel narrowing, grains are deposited on the inside of the pre-existing levee walls, but there is no change to the existing levees, so the total width of the channel does not change, as shown in figure 16(b).
Once the new fully developed steady state for Q M = 70 g s −1 has been established the inflow flux is increased to Q M = 100 g s −1 at t = 150 s as shown in movies 7 and 8. A wave again propagates down the channel, which adjusts its width. This time the wave remobilises the levees on both sides, broadening the central channel and then re-solidifying to form levees that are slightly further out as shown in figure 16(c,d). This process is inherently unsteady and so the outer margins of the new levees are no longer perfectly parallel. The initial wavefront, associated with the change in mass flux, grows in size as it erodes the static material in the pre-existing wide levees and deposits a narrower levee behind it. As a result it moves faster than  the steady uniform flow behind it and therefore detaches from it, producing a brief period where the grains come to rest before the new steady-state flow establishes itself. The computed steady-state thickness and velocity profiles across the channel are shown in figures 16(e) and 16( f ), respectively. These show that, as the mass flux is increased, the width increases to W = 13.7 cm, the thickness barely changes and the velocity profile lies very close to the steady-state value in § 3. All traces of the previous history of the flow that was recorded in the earlier levees are therefore destroyed, which is consistent with the experiments in figure 4.

Unsteady periodic avalanching regime
When the mass flux is reduced to Q M = 16 g s −1 grains pile up near the source, fail periodically and then come to rest again by upslope propagating stopping waves. The initial phases of this solution are very complicated with a sequence of pile collapses that produce two avalanches that propagate down either side of a central ridge, as shown in the online supplementary movies 9 and 10. However, after approximately 95 s the collapses become more structured, producing a single erosion-deposition wave that propagates downslope over a thin deposit that has been laid down by previous flows, as shown in figures 17 and 18. Erosion-deposition waves continuously erode material at their leading edge and deposit material at the rear as they propagate downslope (Edwards & Gray 2015;Edwards et al. 2017). The waves therefore do not represent a single body of material that is released from the collapse, but a constantly changing body of grains. From t = 95-300 s the static deposit has an approximately elliptical shape that becomes progressively elongated with increasing time, with the source centred at the upstream focus point. The erosion-deposition waves are able to propagate easily on top of this deposit and grow in size in the first half of the ellipse expanding laterally to span its width. However, as the width of the elliptical deposit narrows, the avalanche overtops the sides and the material on the static bed travels a short distance before it stops. As a result, as the avalanche continues to propagate downslope, it is progressively deposited at the sides and decreases in size, before it eventually overtops the deposit front and then rapidly comes to rest. The front and sides of the deposit are therefore progressively advanced in a sequence of steps as each avalanche passes by.
The movies also show that as the erosion-deposition wave propagates downslope it detaches from the flowing material at the source and a retrogressive (upslope propagating) stopping wave forms that travels up to the source and brings the grains to rest. As material then builds up at the source it eventually overcomes the static friction and collapses to form the next erosion-deposition wave. In this way a quasi-periodic sequence of avalanches is formed that transport the grains downslope in a series of steps. This process is visualised in figure 19, which shows the thickness at a fixed position x = 1.25 m as a function of time. The front first reaches this position at approximately 90 s. During the passage of the first five or six waves the deposit thickness builds up to just over 6 mm, which lies between h stop = 5 mm and h * = 6.7 mm. The presence of this erodible layer allows the avalanches to propagate easily downslope with a total thickness of approximately 8 mm. The erodible layer reduces the apparent friction experienced by the grains, in the sense that an avalanche of the same volume released on the same slope, but without the erodible layer, would rapidly come to rest. Each of the individual waves have the same shape as the one-dimensional erosion-deposition waves modelled by Edwards & Gray (2015) and develop into a periodic wave train with a period of approximately 2) m that was divided into (n x , n y ) = (400, 300) grid points, this implies that the downslope velocity dropped from the steady uniform value to zero over (n y /L y )W boundary ≈ 3.3 grid cells. Mangeney et al. (2007) used friction parameters corresponding to those in the experiments of Félix & Thomas's (2004), but with a thickness length scale L = 0.8 mm that was larger than the experimental best fit of L = 0.65 mm, leading to a numerical h stop about 25 % larger than that in experiments. Despite this, their simulated flow thickness of H ≈ 1.03h stop = 2.754 mm was still considerably thinner than the 3.12 mm and 3.37 mm measured experimentally by Félix & Thomas (2004) for Q M = 8 g s −1 and Q M = 17.5 g s −1 , respectively. Mangeney et al. (2007) acknowledged this discrepancy, and concluded that '[in] particular, the absence of lateral dissipation between the flowing mass and the quasi-static shoulders likely leads to underestimation of the flowing thickness and the thickness of the deposit in the central channel'.
As shown in this paper, lateral dissipation (viscosity) is required to produce a unique steady channel thickness, but a hysteretic friction law with κ = 1 (a gradual transition 632 F. M. Rocha, C. G. Johnson and J. M. N. T. Gray from static to dynamic friction coefficient) is also required, in order for this unique thickness to be significantly greater than both h stop and h * , as observed in experiments.
As has been seen in § 3, the depth-averaged viscous terms in (2.3)-(2.4) are a singular perturbation of the model that produces qualitatively different solutions. Rather than the unique solution found when the viscous terms are present, in the inviscid case (ν = 0) there are an infinite family of steady-state solutions with constant channel thicknesses in the metastable range [h stop , h start ], plug-like downslope velocity profiles across the slope and contact discontinuities at the channel edge ( § 3.2). Mangeney et al.'s (2007) numerical simulations are consistent with a flow that is approaching one of these steady-state inviscid solutions. Indeed, Mangeney et al. (2007) derived equations for the velocity and width of the flow (their equations (24) and (25)), which are directly equivalent to the steady-state equations (3.9) and (3.11) for glass beads (Γ = 0). However, they did not draw out the fact that the steady-state inviscid system of equations is not uniquely determined, nor why the observed thickness H ≈ 1.03h stop is so close to the thinnest possible steady-state solution.
In the absence of any numerical viscosity, the initial conditions and time-dependent evolution of the inviscid equations could, in principle, select the observed channel width from the many possible steady-state solutions. In reality, this time dependence does not appear to be the mechanism setting the channel width, since Deboeuf et al. (2006) found that the width was independent of past changes in flux and Takagi et al. (2011) showed that an initial layer of erodible particles covering the chute did not influence the equilibrium flowing channel width.
Based on solutions of Mangeney et al.'s (2007) inviscid equations and quantitative reproduction of their observations made using the numerical scheme of this paper, there appear to be two reasons why the numerical solutions of Mangeney et al. (2007) select a flow thickness very close to h stop . These relate to the Mangeney et al.'s (2007) use of an extremely rapid transition from static to flowing friction (κ = 10 −3 ) and the interaction of this with small but unavoidable errors in numerical schemes.
Firstly, as pointed out by Edwards et al. (2019), the initial sharp decrease in the friction as Fr increases from zero, when κ = 10 −3 , is extremely difficult to represent in standard floating-point numerical computations. For instance, at the smallest non-zero value the Froude number (∼10 −16 , based on IEEE-754 double precision round-off error of O(1) numerical quantities in the finite-volume scheme) the friction coefficient drops to a value that is well below the maximum static friction and is insufficient to bring to rest any flowing material thicker than 1.02h stop = 2.727 mm. This suggests that when Mangeney et al. (2007) stopped their numerical simulations (when H ∼ 1.03h stop ) the flow was still slowly spreading, but that it was very close to the range of thicknesses where inviscid steady-state solutions could exist.
Secondly, the singular nature of the viscous terms in the governing equations means that even a small amount of numerical viscosity results in the solution of the viscous equations, with the viscosity coefficient in this case determined artificially by the grid resolution and details of the numerical scheme. As noted previously, when κ = 10 −3 the unique steady viscous solution for a given flux lies in a very narrow range of thicknesses [H min , H max ] that is slightly (∼0.1 %) greater than the thinnest steady uniform flow (i.e. h stop , for the friction law used by Mangeney et al. 2007). It is therefore possible that Mangeney et al.'s (2007) inviscid solutions are approaching a steady state of the viscous equations, but this relies entirely on the mesh-size dependent and numerical scheme dependent viscosity.
In contrast, the combination in this paper of physically derived viscous terms and the friction law of Edwards et al. (2019) (with κ = 1) provides (i) a unique steady channel thickness and width in good agreement with experiments, and (ii) a timedependent model that unambiguously approaches this unique steady state at late times.

Summary of results
This paper shows that two physical processes are needed to quantitatively predict self-channelisation and levee formation in monodisperse granular flows using depth-averaged avalanche models. These are (i) a non-monotonic effective basal friction law (Pouliquen & Forterre 2002;Félix & Thomas 2004;Deboeuf et al. 2006;Mangeney et al. 2007;Edwards & Gray 2015;Edwards et al. 2017Edwards et al. , 2019 and (ii) depth-averaged lateral viscous stresses (Gray & Edwards 2014;Baker et al. 2016a). The non-monotonic friction models the hysteretic behaviour that allows static and flowing regions to coexist, while the viscous terms are vital to close the system of equations and produce smooth transitions between the static and flowing regions that do not exert shear stresses on the levee walls.
The two-dimensional time-dependent simulations in figures 14 and 15 show that the resultant equations (2.2)-(2.15) model the formation of a central flowing channel bounded by parallel-sided static levees. When material flowing down the channel reaches the front, it loses its lateral confinement, spreads out and slows down. The grains that are still in the centre of the flow are then overrun, while the particles that migrate to the sides come to rest and build a new section of the levee wall, advancing it downslope (Félix & Thomas 2004;Johnson et al. 2012). In this way the flow front is able to propagate downslope as a steadily travelling wave, continuously laying down a pair of levees at the back of the flow head. The final width of the channel is not always set immediately behind the head; an equilibrium width of the central channel is established when the flow thickness and depth-averaged velocity profile are such that there is no net erosion or deposition at the levee-channel boundaries. Once the channel has formed, the levees prevent the flow spreading laterally and thereby maintain the flow's thickness, and hence mobility, leading to enhanced run-out. The flow can widen by levee-bank overtopping in response to an increase in mass flux, and retreat into the centre of the channel if the flux is reduced, as shown experimentally (figure 4, movie 3) and in the numerical simulations (figure 16, movies 7 and 8).
To understand the force balance leading to the leveed channel, a steady-state problem is formulated in § 3.1, which shows that across the width W of a parallel-sided flowing channel, the thickness h is equal to a constant H and the downslope velocity profileū(y) satisfies the second-order ODE (3.5). For a given mass flux Q M , this can be solved as a boundary value problem for W, H andū(y), subject to zero velocity (3.6) and zero shear stress (3.7) boundary conditions at the levee-channel interfaces, as well as the integral constraint (3.8) on the mass flux. The viscous terms allow the static shear-stress-free levees to be connected to the central flow, where steady uniform flow dominates, across a viscous boundary layer or shear band (see figure 9c). There is a unique solution to the problem for all values of the mass flux Q M although the character of the solution changes at Q = Q min . For Q M < Q min the solutions are purely in the intermediate frictional regime (2.12), whilst for Q M > Q min , the solutions also have regions of dynamic friction (2.10). It is these latter solutions that are physically realised.
The solutions to the boundary value problem (see figures 8-9) are in excellent quantitative agreement with the experimental measurements of the height and width of the central flowing channel (Takagi et al. 2011) as shown in figure 10 for sand. These solutions show that, as the mass flux Q M is increased, the thickness of the central flowing channel remains nearly constant and its width increases approximately linearly. The constant thickness H in the channel is significantly above both h stop and h * . Takagi et al. (2011) also made measurements of the surface velocity profile across the channel. These are harder to compare, because the velocity shear through the depth of the avalanche necessarily implies that the surface velocity is not equal to the depth-averaged velocity calculated by the theory. Figure 11 shows that the depth-averaged velocityū and the surface velocity measurements u s have the same functional behaviour with the mass flux Q M , and that u s /ū ≈ 2.35. This large ratio of surface velocity to depth-averaged velocity is consistent with the discrete element method simulations of Silbert et al. (2003) and non-local theory of Kamrin & Henann (2015), which predict a transition from Bagnold flow to exponential-like velocity profiles as the flow gets close to h * .
The model is also in very good agreement with the thickness measurements of Félix & Thomas (2004) for glass beads as shown in figure 12 for a wide range of mass fluxes. However, the predictions of the steady-state minimal total width of the channel (see § 3.5.2 and appendix B) do not agree with those measured by Félix & Thomas (2004) at high mass fluxes. The most likely explanation for this discrepancy is that either Félix & Thomas's (2004) experimental chute was not long enough, or the flow was not observed for a sufficiently long time, to measure a true steady state. This is certainly consistent with the observations of Deboeuf et al. (2006) and Takagi et al. (2011) who found that for weakly hysteretic materials, such as glass beads, the width of the flow adjusted on very long time scales (∼70-90 min).
Below a critical mass flux the steady-state solutions go unstable in experiments (Félix & Thomas 2004;Takagi et al. 2011) as shown for red sand in movie 4. This also occurs in the numerical solutions (figures 17, 18 and 19), which produce a series of finite pulses of flowing material separated by intervals during which the grains are completely static. These are another example of erosion-deposition waves, which have previously been observed on erodible beds (Daerr & Douady 1999;Börzsönyi et al. 2005;Clément et al. 2007;Edwards & Gray 2015;Edwards et al. 2017) and have been simulated numerically by Edwards et al. (2017) using the same system of equations used in this paper. The governing equations (2.2)-(2.15) are therefore extremely versatile and have the potential to explain many flows with qualitatively very different behaviour.
The viscous terms play a critical role in setting the velocity profile and hence the thickness and width of the self-channelised flow. In the absence of viscosity, the system of equations (2.2)-(2.5) reduce to the classical avalanche equations, which are hyperbolic in character (see e.g. Grigorian et al. 1967;Savage & Hutter 1989;Gray et al. 1999Gray et al. , 2003Mangeney et al. 2007). As shown in § 3.2, in this inviscid model the constant channel thickness H implies that the velocity is equal to the steady-uniform value (3.9) across the whole channel, with contact discontinuities at the levee-channel interfaces where the downslope velocity jumps from the steady-uniform valueū steady to zero. The thickness of the flow H is not determined, and may lie anywhere in the range [h * , h start ], where h * > h stop is the minimum observable thickness of a steady uniform flow . As a result, even for the same mass flux Q M , there are a corresponding range of possible widths (see figure 7), which implies that inviscid theories lack the physical mechanism to set the steady-state thickness and width of the flow.
In practical numerical computations there is usually some mesh-size-dependent numerical diffusion. It follows that numerical computations with Mangeney et al.'s (2007) inviscid model will appear to behave somewhat like the theory presented in this paper (see 4.5). The important distinction is that it is the mesh-size-dependent numerical diffusion that selects the solution, rather than the physically based depth-averaged viscosity (Gray & Edwards 2014;Baker et al. 2016a) derived from the µ(I)-rheology (GDR-MiDi 2004;Jop et al. 2006). The viscous depth-averaged theory presented in this paper is therefore a significant advance over previous models.
Despite the theory's clear success in simulating many features of the flow, one subtle aspect that is not captured correctly is the shape of the free surface. Experimental measurements of self-channelised flows (Félix & Thomas 2004;Deboeuf et al. 2006;Takagi et al. 2011) show that, instead of having a constant depth across the channel, the maximum thickness occurs in the centre, where the material moves faster. The curvature in the free surface appears to coincide with regions of high shear, i.e. for narrow channels, with a parabolic-like velocity profile, experiments suggest that the free surface is curved across the whole channel. For wider channels, the free surface is almost flat in the central steady uniform region, with curvature observed only at the boundary layers at the sides (Félix & Thomas 2004). This observation may be important in the development of new constitutive equations for granular flows that go beyond the µ(I)-rheology (see e.g. Kamrin  Another possible explanation is that the low velocity regions, on the inside walls at the base levees, jam and become quasi-static rather than inertial (Deboeuf et al. 2006;Kokelaar et al. 2014). This would change the height of the effective topography that the avalanche was flowing over and potentially allow cross-slope variations in thickness. Understanding how to model this quantitatively would probably require a compressible granular theory (see e.g. Barker et al. 2017;Heyman et al. 2017;Schaeffer et al. 2019) that could span the quasi-static and inertial regimes (Chialvo, Sun & Sundaresan 2012). There is, however, a long way to go before either of these effects can be included in depth-averaged avalanche models, and the existing theory appears to be able to make accurate predictions without them.

Implications for the interpretation of levee-channel deposits
The results of this paper provide a theoretical framework for the interpretation of leveed deposits from geophysical mass flows such as debris flows and pyroclastic density currents. Specifically, the predictions of the thickness and width of a self-channelised leveed flow in § 3.3 may be inverted to infer information about the mass flux or rheology of the original flow from measurements of the deposit.
A key observation is that the thickness of a levee is similar to the thickness of the flow that created it, and is nearly independent of the mass flux of that flow (figures 10a and 12a). A difference in thickness between levees in successive flows therefore suggests a change in flow composition and rheology, not simply mass flux. The width of the static channel between the levees is roughly proportional to the mass flux of the originating flow (figures 10b and 12b). Very wide levees, with a sequence of nested parallel levees inside the main channel, are strongly indicative of a decrease in mass flux prior to deposition (figure 4b). However, an increase in flux with time may not be reflected in the deposit, due to erosion of the existing levees (figure 4c).
Interpretation of deposits with levee-channel morphology is complicated by the existence of other levee-forming flows. A disturbance to an erodible layer of grains on an inclined plane may lead to a localised erosion-deposition wave, in which the static layer is continuously eroded at the flow front and deposited a short time later at the back of the avalanche ). This flow leaves behind a pair of static parallel-sided levees that are raised above the surface of the original erodible layer, located either side of a depression or trough. The leveed deposit resulting from this erosion-deposition wave therefore looks similar to that produced by a continuous flux of grains (figure 2c), which might make the two cases difficult to distinguish in geophysical flow deposits. A key difference between these flows is the origin of the material that forms the levees. The levees described in this paper are composed of grains provided by a continuous source upslope, whereas those generated by the erosion-deposition waves of Edwards et al. (2017) are composed of grains that are continuously eroded at the front of the avalanche.
Field observations (Wilson & Head 1981;Bartelt et al. 2012) suggest that the margins of the flow may be differentially fluidised during emplacement of levees, which implies that the rheology is inhomogeneous (Iverson & Vallance 2001;Iverson 2003). In debris flows, the evolving pore fluid pressure (Iverson 1997;Iverson & Denlinger 2001) and particle-size distribution due to segregation (Gray 2018) significantly alter the bulk flow by generating drier bouldery margins that resist the motion and form pronounced levees . Although these effects present additional modelling challenges, the heterogeneous particle-size distribution of a deposit also potentially provides considerably more information about the flow that created it than the deposit geometry alone. The incorporation of such heterogeneity into the model for levee formation presented here is, therefore, a promising avenue for future research.
provides a good fit to the observations in this paper, where λ is an adjustable parameter. Evaluating this at the free surface implies that the surface velocity u s is related to the depth-averaged velocityū by the relation The surface velocity data for sand (Takagi et al. 2011) shown in figure 11 suggest that u s = 2.35ū, which implies λ = 2.05. This is within 15 % of the value u s = 2.7ū (corresponding to λ = 2.45) measured experimentally in a thin unconfined flow of glass beads .