Formation of levees, troughs and elevated channels by avalanches on erodible slopes

Snow avalanches are typically initiated on marginally stable slopes with a surface layer of fresh snow that may easily be incorporated into them. The erosion of snow at the front is fundamental to the dynamics and growth of snow avalanches and they may rapidly bulk up, making them much more destructive than the initial release. Snow may also deposit at the rear, base and sides of the flow and the net balance of erosion and deposition determines whether an avalanche grows or decays. In this paper, small-scale analogue experiments are performed on a rough inclined plane with a static erodible layer of carborundum grains. The static layer is prepared by slowly closing down a flow from a hopper at the top of the slope. This leaves behind a uniform-depth layer of thickness $h_{stop}$ at a given slope inclination. Due to the hysteresis of the rough bed friction law, this layer can then be inclined to higher angles provided that the thickness does not exceed $h_{start}$ , which is the maximum depth that can be held static on a rough bed. An avalanche is then initiated on top of the static layer by releasing a fixed volume of carborundum grains. Dependent on the slope inclination and the depth of the static layer three different behaviours are observed. For initial deposit depths above $h_{stop}$ , the avalanche rapidly grows in size by progressively entraining more and more grains at the front and sides, and depositing relatively few particles at the base and tail. This leaves behind a trough eroded to a depth below the initial deposit surface and whose maximal areal extent has a triangular shape. Conversely, a release on a shallower slope, with a deposit of thickness $h_{stop}$ , leads to net deposition. This time the avalanche leaves behind a levee-flanked channel, the floor of which lies above the level of the initial deposit and narrows downstream. It is also possible to generate avalanches that have a perfect balance between net erosion and deposition. These avalanches propagate perfectly steadily downslope, leaving a constant-width trail with levees flanking a shallow trough cut slightly lower than the initial deposit surface. The cross-section of the trail therefore represents an exact redistribution of the mass reworked from the initial static layer. Granular flow problems involving erosion and deposition are notoriously difficult, because there is no accepted method of modelling the phase transition between static and moving particles. Remarkably, it is shown in this paper that by combining Pouliquen & Forterre’s (J. Fluid Mech., vol. 453, 2002, pp. 133–151) extended friction law with the depth-averaged $\unicode[STIX]{x1D707}(I)$ -rheology of Gray & Edwards (J. Fluid Mech., vol. 755, 2014, pp. 503–544) it is possible to develop a two-dimensional shallow-water-like avalanche model that qualitatively captures all of the experimentally observed behaviour. Furthermore, the computed wavespeed, wave peak height and stationary layer thickness, as well as the distance travelled by decaying avalanches, are all in good quantitative agreement with the experiments. This model is therefore likely to have important practical implications for modelling the initiation, growth and decay of snow avalanches for hazard assessment and risk mitigation.

Snow avalanches are typically initiated on marginally stable slopes with a surface layer of fresh snow that may easily be incorporated into them. The erosion of snow at the front is fundamental to the dynamics and growth of snow avalanches and they may rapidly bulk up, making them much more destructive than the initial release. Snow may also deposit at the rear, base and sides of the flow and the net balance of erosion and deposition determines whether an avalanche grows or decays. In this paper, small-scale analogue experiments are performed on a rough inclined plane with a static erodible layer of carborundum grains. The static layer is prepared by slowly closing down a flow from a hopper at the top of the slope. This leaves behind a uniform-depth layer of thickness h stop at a given slope inclination. Due to the hysteresis of the rough bed friction law, this layer can then be inclined to higher angles provided that the thickness does not exceed h start , which is the maximum depth that can be held static on a rough bed. An avalanche is then initiated on top of the static layer by releasing a fixed volume of carborundum grains. Dependent on the slope inclination and the depth of the static layer three different behaviours are observed. For initial deposit depths above h stop , the avalanche rapidly grows in size by progressively entraining more and more grains at the front and sides, and depositing relatively few particles at the base and tail. This leaves behind a trough eroded to a depth below the initial deposit surface and whose maximal areal extent has a triangular shape. Conversely, a release on a shallower slope, with a deposit of thickness h stop , leads to net deposition. This time the avalanche leaves behind a levee-flanked channel, the floor of which lies above the level of the initial deposit and narrows downstream. It is also possible to generate avalanches that have a perfect balance between net erosion and deposition. These avalanches propagate perfectly steadily downslope, leaving a constant-width trail with levees flanking a shallow trough cut slightly lower than the initial deposit surface. The cross-section of the trail therefore represents an exact redistribution of the mass reworked from the initial static layer. Granular flow problems involving erosion and deposition are notoriously difficult, because there is no accepted method of modelling the phase transition between static and moving particles. Remarkably, it is shown in this paper that by combining Pouliquen & Forterre's (J. Fluid Mech., vol. 453, 2002, pp. 133-151) extended friction law with the depth-averaged µ(I)-rheology

Introduction
Large-scale field measurements of snow avalanches (Sovilla, Sommavilla & Tomaselli 2001;Sovilla & Bartelt 2002;Sovilla, Burlando & Bartelt 2006) indicate that the dominant mechanism for entrainment and growth of an avalanche is via frontal ploughing into a layer of fresh snow. Sovilla et al. (2006) showed that over eighteen events the avalanche mass typically increased by a factor of four, which had a significant effect on the flow dynamics as well as the run-out distance. The primary limiting factors on how much snow can be entrained are the structure of the snowpack and the amount of available material. Virtually all of the entrainment occurred at, or very close to, the flow front, with basal erosion by the body of the flow being much less significant. The process of frontal entrainment allows the avalanche to bulk up, and hence become much more destructive than the initial release. As well as eroding material at the front, snow may also be deposited at the rear, base and sides, and there is a subtle balance which determines whether there is overall growth or decay in the total avalanche mass. Deposition can also produce key qualitative changes to the flow dynamics, such as the formation of static levees at the flanks of the avalanche, which leave a trail in the deposit as shown in figure 1. This image also shows that between the levees there is a trough that has been incised down below the initial height of the snow cover, and which contains material that had previously been mobilized by the avalanche. The flow may also rapidly stop if the slope inclination changes and the avalanche decays away through mass deposition (Bartelt, Buser & Platzer 2007).
In order to investigate the effects of erosion and deposition experimentally, we perform small-scale analogue experiments on a rough inclined slope using carborundum particles of 280-350 µm. For angles below ζ 2 = 47.5 • , the maximum angle for steady uniform flow, a static layer of erodible particles forms as the grains come to rest at a thickness h stop (Pouliquen 1999b). Due to the hysteresis in the rough bed friction law Pouliquen & Forterre 2002), once the grains have stopped the slope can be inclined to a higher angle before the particles begin to flow again, and thicker layers of particles are stable provided their depth lies below h start > h stop . In the experiments, an avalanche is initiated on top of the erodible layer by releasing a finite mass of carborundum particles from a small aspect ratio hollow cylinder. Figure 2 shows a typical flow on a slope inclined at ζ = 35.2 • to the horizontal and covered with a uniform-depth layer of erodible grains of thickness h stop . A long exposure time is selected so that the moving grains appear slightly blurred and the static grains can be seen in sharp focus. The avalanche has a crescentic  shape, with a steep erosive front and a tail that decays more gradually in height. As the avalanche flows downslope, it continuously erodes particles at the front and deposits them at the rear, sides and base, to form a trail behind the flow, with an incised trough, that lies below the initial static deposit height and which is flanked by static raised levees. These are subtle features that can be seen in figure 2 due to the use of oblique lighting. In this case, the avalanche has organized itself into a three-dimensional, steadily travelling erosion-deposition wave (Edwards & Gray 2015) that can propagate indefinitely along the slope, so long as the erodible layer ahead of the wave does not change depth and the inclination remains the same. The formation of an incised trough and levees is closely similar to what was observed in the real snow avalanche in figure 1. This suggests that there is a strong link between the analogue experiments and real field-scale observations. Permanent-form solitary waves on erodible beds were first observed by Daerr (2001). Börzsönyi, Halsey & Ecke (2005 made further observations and postulated a Burgers-type model, but unfortunately solitary n-wave-type solutions decay in height with increasing time (Whitham 1974;Edwards & Gray 2015), so this approach cannot describe erosion-deposition waves. Clément et al. (2007) developed a partial fluidization theory with two equations to describe the evolution of the flow thickness and an order parameter that governs the amount of fluidization. This model does have a family of solitary wave solutions, but the structure of the equations is radically different to those one might expect as a generalization of the shallow-water-like avalanche models that are conventionally used to model geophysical mass flows Formation of granular levees, troughs and elevated channels 281 FIGURE 2. Oblique view of an erosion-deposition wave travelling downslope on a plane inclined at an angle ζ = 35.2 • to the horizontal. The right-hand side of the wavefront appears brighter due to oblique illumination from the downslope end of the plane. A long time exposure has been used so that moving grains are blurred and the static regions are sharply in focus. The width of this avalanche is approximately 8.5 cm across the wave crest. Behind the elevated front of the flow lateral levees are deposited on either flank and between them a trough that is slightly beneath the level of the original deposit is left behind on the erodible bed.
Early models for snow avalanches (Briukhanov et al. 1967;Grigorian et al. 1967;Eglit & Demidov 2005) recognized the importance of frontal entrainment and solved for the motion of the front by treating it as a shock wave separating a finite-depth solid-like erodible layer from a rapidly moving liquid-like avalanche. Mass and momentum jump conditions were then used to solve for the entrainment of mass and the speed of front propagation. This approach requires the moving interface between the static snow cover and the avalanche to be tracked as part of any computation. To avoid this, Eglit (1983) proposed a model with more gradual erosion at the base of the avalanche, with the rate being proportional to the flow speed, which is highest at the front. Most current avalanche models try to incorporate erosion and deposition by introducing mass and momentum source terms into the depth-averaged mass and momentum balances (e.g. Bouchaud et al. 1994;Douady, Andreotti & Daerr 1999;Gray 2001;Doyle et al. 2007;Tai & Kuo 2008;Gray & Ancey 2009;Iverson 2012;Tai & Kuo 2012) and sometimes to augment them with an energy balance equation (e.g. Bouchut et al. 2008;Capart, Hung & Stark 2015). These approaches are notoriously difficult since additional non-trivial empirical relations are required to close the models and there may be slow creep deep beneath the avalanche (Komatsu et al. 2001), which makes the interface between the flowing and static regions hard to define.
Recently Edwards & Gray (2015) have developed a one-dimensional system of depth-averaged equations to model erosion-deposition waves, which combines the extended basal friction law of Pouliquen & Forterre (2002) with the depth-averaged µ(I)-rheology of . The avalanche is assumed to be either moving or static throughout its entire depth, with the basal friction law essentially determining which regions are in motion. Although this approximation is crude, in that it does not explicitly resolve the basal entrainment or deposition, Edwards & Gray (2015) showed that it was able to predict accurately the amplitude, wavelength and coarsening dynamics of erosion-deposition waves that spontaneously form in a long channel from a continuous inflow. They observed that typical waves had a ratio of peak height to static layer depth of 2.6 and a typical mobile wavelength of 59 cm, which are both very similar to the waves observed by Takagi, McElwaine & Huppert (2011) in their low-inflow-rate experiments.
A key element of the  theory is the introduction of a depth-averaged viscous term, which represents a singular perturbation to the equations. Strong support for this approach is provided by the fact that the depth-averaged µ(I)-rheology is also able to predict the cutoff frequency of roll waves (Forterre & Pouliquen 2003;Forterre 2006) without the need for any additional fitting parameters, since the granular viscosity is fully determined ). This contrasts with inviscid avalanche models that incorrectly predict the growth of granular roll waves at all frequencies above the critical Froude number, Fr c > 2/3. Baker, Barker & Gray (2016a) have generalized the depth-averaged µ(I)-rheology to two dimensions by using the principle of frame invariance. Their extended model is able to capture the depth-averaged downslope velocity profile that develops across a fully developed chute flow with either wall slip or zero velocity at the rough side walls. Again, this is something that the inviscid theory is unable to capture. Moreover, for sufficiently wide channels the reconstructed two-dimensional downslope velocity profile across the channel compares favourably with steady-state simulations using the full (non-depth-averaged) µ(I)-rheology (Jop, Forterre & Pouliquen 2006), as well as with experimental measurements of the surface velocity.
In this paper we focus on using the Baker et al. (2016a) extension to the depth-averaged µ(I)-rheology to model three-dimensional erosion-deposition waves. Early observations of such waves started with  who triggered two types of avalanche behaviour on a rough bed with an erodible layer of grains. The static layer was prepared by pouring glass beads down a slope covered with a velvet cloth, which left a deposit of thickness h stop on a slope of inclination ζ stop . They observed that this static layer could be inclined to a higher angle ζ start before it would spontaneously start flowing again, i.e. there is hysteresis. For a small increase in the slope angle, between ζ stop and ζ start , a small perturbation to the static layer at a single point created an avalanche that propagated downslope entraining mass and widening as it flowed. This left a triangular pattern in the deposit, which led  to name these phenomena 'triangular' avalanches. They also showed that on less stable slopes the avalanching region propagated upwards as well, so that eventually the whole erodible layer was set in motion. Pouliquen (1999a) measured h stop as a function of the inclination angle ζ , and performed experiments to determine the thickness h and the depth-averaged velocitȳ u during steady uniform flow. He found that all the experimental data collapsed with the scaling Fr = βh/h stop , where the constant of proportionality β = 0.136. Note that Pouliquen (1999a) defined the Froude number as Fr = |ū|/ √ gh, where |ū| was the depth-averaged speed and g was the constant of gravitational acceleration. Using this relation to substitute for h stop , Pouliquen (1999a) determined an empirical friction law for rough beds provided Fr > β.  showed that although the rough bed friction law looks like a basal friction law it actually encodes the leading-order behaviour of the µ(I)-rheology through the depth of the granular avalanche. Pouliquen & Forterre (2002) combined Pouliquen's (1999a) dynamic rough bed friction law with the  concept that once the grains were stopped (i.e. Fr = 0) flow could not start again until the higher angle ζ start was exceeded. For the intermediate regime 0 < Fr < β there was no further experimental information, so Pouliquen & Forterre (2002) suggested a power law extrapolation between the two laws. Pouliquen & Forterre (2002) applied their extended friction law to the release of a finite mass of granular material on a rough surface and found that an inviscid shallowwater avalanche model was able to predict accurately the position and time at which it stopped, as well as the overall spreading. They also performed the same experiment on a bed of erodible grains. In this configuration they found that the experimental avalanche spread out, and formed a drop shape similar to that observed by Daerr (2001), with a sharp front and gradual decrease in height in the tail that propagated steadily downslope, eroding and depositing material as it flowed. On the other hand, their inviscid simulations predicted that the flow did not approach a steadily travelling wave. Instead, the computed solution spread out considerably more, so that it was thinner and faster than the experiment with a large almost constant thickness region. This may be an indication that the computation was partially locking onto the steady uniform flow solution to the intermediate flow regime found by Edwards & Gray (2015). In the light of the successful prediction of the formation of two-dimensional steadily travelling erosion-deposition waves (Edwards & Gray 2015), using the depthaveraged µ(I)-rheology , this paper investigates whether their model is able to more accurately predict what happens in three dimensions.

Experimental set-up and methods
To investigate the formation of three-dimensional erosion-deposition waves, smallscale experiments are performed to highlight some key qualitative features not previously described. The experimental set-up is shown in figure 3 and consists of a plane inclined at an angle ζ to the horizontal, which is roughened by attaching a monolayer of spherical glass beads of diameter 750-1000 µm using double-sided sticky tape. A coordinate system Oxyz is defined with the x-axis pointing downslope, the y-axis across the slope, with y = 0 at the midpoint of the width of the plane, and the z-axis pointing normal to the rough plane at z = 0. The avalanche is assumed to have depth-averaged velocity componentsū andv in the downslope x and cross-slope y directions, respectively. The inclined plane is 1.5 m long by 0.5 m wide and has a hopper and gate mechanism at the top, which allows the bed to be coated with a static layer of 280-350 µm diameter carborundum particles prior to the start of the experiments. For a slope at an angle ζ 0 to the horizontal, the initial static layer is of thickness h 0 = h stop (ζ 0 ), i.e. the same thickness as that measured by  and Pouliquen (1999a). Once the grains have stopped the hysteresis of the friction law allows the chute to be inclined to a new angle ζ < ζ start that is steeper than ζ 0 = ζ stop without the grains becoming mobile again.
The release point is centred at (x, y) = (0, 0), which is 80 cm from the end of the chute. A finite volume of 280-350 µm carborundum particles is then placed on top of the static layer by filling a hollow cylinder of radius R = 1.4 cm and height h c = 1.6 cm. The flow is released by raising the cylinder, which causes the particles inside to spread from the downslope edge, like a small inclined column collapse (Mangeney et al. 2010), before forming an avalanche that travels down the plane whilst eroding the stationary layer at the wavefront and depositing particles behind. Figure 2 shows an oblique view of one of these avalanches as it travels downslope. In the trail left behind the flowing avalanche there are small static levees and a shallow trough, made visible by illuminating obliquely from the downstream end of the chute.
To capture some of the subtle features of the flow, a high-speed camera (Teledyne DALSA Genie HM1400) is used to photograph it from overhead at a rate of 100 f.p.s. with oblique illumination from the downstream end of the chute. A space-time plot is constructed by extracting the middle row of pixels (along y = 0) from each image and combining them into a single plot with elapsed timet along the abscissa and distance x along the ordinate. Note that elapsed timet is offset from the real time t by a different unknown amount for each set of images. In addition, a Micro-Epsilon scanCONTROL 2700-100 laser profile sensor is used to acquire thickness data aligned along the y-axis at three different downstream positions x L = 12 cm, x L = 30 cm and x L = 50 cm. The laser line measures the distance of the bulk flow particles away from the sensor at a frequency of 100 Hz and to an accuracy of ±0.2 mm (approximately a grain diameter) by laser triangulation. The thickness profile h of the avalanche in the z-direction can then be calculated along the laser line by measuring the distance between the sensor and the bed before it is coated with the static layer of carborundum particles. Space time contour plots of the thickness are then built up for each position x L , although the experiment has to be repeated for each location. Experiments were carried out for various slope inclination angles and static layer depths, and three qualitatively different behaviours were observed, which are described in greater detail in the following sections.
Formation of granular levees, troughs and elevated channels 285 x y R O z FIGURE 3. A wooden plane with a layer of 750-1000 µm spherical glass beads stuck to the surface is inclined at an angle ζ to the horizontal. A coordinate system Oxyz is centred at the release point with the x-axis pointing downslope, the y-axis across the slope and the z-axis is aligned with the upward pointing normal. The chute is prepared with a constant-depth layer of 280-350 µm diameter carborundum particles of thickness h 0 = h stop (ζ 0 ) and a finite volume of the same grains is released on top of this layer from a hollow cylinder of radius R = 1.4 cm and height h c = 1.6 cm centred at the origin (x, y) = (0, 0). Note that the hysteresis of the basal friction allows the inclination angle ζ to be greater than the preparation angle ζ 0 of the constant-depth layer provided ζ < ζ start .

Growing avalanches
In order to generate an avalanche that has a net increase in mass as it propagates downslope the plane is first inclined at an angle of ζ = 34.4 • and the bed is coated with a static layer of carborundum that has a thickness h 0 = h stop (34.4 • ) ≈ 2.2 mm. Without disturbing the grains in the static layer, the slope angle is then increased to a higher angle of ζ = 35.2 • . The deposit layer then is thicker than the usual h stop at 35.2 • by virtue of the hysteresis in the friction law in the range ζ ∈ [ζ stop , ζ start ]. The cylinder is held on top of the layer, filled with grains and the avalanche is then released.
Overhead images of the avalanche on the inclined plane and a space-time plot constructed from them are shown in figure 4. There is also a movie of the flow in the online supplementary material available at https://doi.org/10.1017/jfm.2017.309 (movie 1). After the initial release the avalanche rapidly develops a sharp front that erodes material and a tail that decreases in height much more gently, and from which grains are continually being deposited. Overall there is net erosion and the avalanche grows in size, leaving behind it a trough formed from the previously mobilized grains, which is thinner than the initially prepared layer. The vertical grey lines on the space-time plot in figure 4( f ) indicate stationary grains and the flow front appears as a strong line across the plot due to the oblique lighting from the downstream end of the chute. This indicates that the avalanche travels at a near constant wavespeed of u w ≈ 0.26 m s −1 although there is evidence of it beginning to slightly accelerate towards the end of the plane. Data from the laser profilometer are shown in figure 5 as (a) contours of thickness h in the (t, y)-plane and (b) as thickness along the centre of the plane y = 0 at each of the distances x L . After the initial collapse the avalanche propagates downslope and both the width W and the wave peak height h w continually increase to the end of the plane. The measured values are W ≈ 11 cm, h w ≈ 5.6 mm at x L = 30 cm and W ≈ 14 cm, h w ≈ 5.9 mm at x L = 50 cm (figure 5). Since the overall length is approximately the same there is a net increase in avalanche mass and a shallow trough forms behind the avalanche. The trough widens with downstream distance as shown in the overhead photos in figure 4(e), which makes it analogous to the 'triangular avalanches' first described by . The deposited layer depth in the trough, denoted by h deposit , continually decreases to a minimum value of h deposit ≈ 1.7 mm at x L = 50 cm. This is close to the static layer thickness h stop (35.2 • ) ≈ 1.7 mm.

Decaying avalanches
Next the plane is inclined at an angle of ζ = 34.1 • to the horizontal and coated with a static layer of thickness h 0 = h stop (34.1 • ) ≈ 2.5 mm. A sequence of overhead images of the avalanche and a space-time plot constructed from them are shown in figure 6. A movie is also available in the online supplementary material (movie 2).  The avalanche decreases in width as it travels downslope, depositing an elevated channel behind it with levees along its flanks, before it finally comes to rest.
It can be seen that after the initial release of grains the avalanche travels downslope for a short distance before it stops abruptly. The nearly straight diagonal line in the space-time plot between 0.15 x + R 0.45 m indicates that typical wavespeeds are u w ≈ 0.13 m s −1 and that there is a very rapid deceleration in the final stages before the avalanche abruptly comes to rest at x = 55 cm downslope. As the avalanche flows down the inclined plane it erodes mass at the front, but deposits far more mass behind it, to form an elevated channel that lies above the height of the undisturbed deposit and is flanked by levees on either side. To make this precise, figure 7, shows (a) thickness contours in the (t, y)-plane and (b) the thickness h along the centre of the plane y = 0 at each of the x L locations. The grains spread to a maximum width of W ≈ 8.5 cm at x L = 12 cm, before reducing to a width of W ≈ 5.5 cm when the avalanche reaches x L = 50 cm (figure 7a). The narrowing of the deposit with increasing downstream distance can also be clearly seen in figure 6(e). The wave peak height h w ≈ 5.1 mm throughout the flow. However, the deposit thickness in the elevated channel increases significantly above the static layer thickness h stop ≈ 2.5 mm, to h deposit ≈ 2.7 mm at x L = 30 cm and h deposit ≈ 3.3 mm at x L = 50 cm (figure 7b). The deposit is able to remain static by virtue the hysteresis in the friction law Pouliquen & Forterre 2002) since h stop < h deposit < h start . The fact that the wavespeed u w is almost constant during most of the flow may be related to the fact that the wave peak height is the same at each of the downstream locations

Solitary steadily travelling avalanches
Finally, the plane is inclined at an angle of ζ = 35.2 • and the bed is coated with a static layer of thickness h 0 = h stop (35.2) ≈ 1.6 mm. A series of overhead images of the avalanche propagating down the plane and a space-time plot constructed from them are shown in figure 8. A movie is also available in the online supplementary material (movie 3). After the initial release the avalanche rapidly develops into a steadily travelling wave, so that byx = x + R ≈ 15 cm it remains at a constant width and has a near constant wavespeed of u w ≈ 0.20 m s −1 for the remaining length of the experimental plane. This suggests that the frontal erosion and deposition at the base, sides and tail of the avalanche are in exact balance, i.e. it has developed into a steadily travelling solitary wave similar to the two-dimensional erosion-deposition waves observed by Edwards & Gray (2015). In this case, however, the wave has a well defined three-dimensional crescentic shape and it leaves behind static levees and a trough incised below the height of the undisturbed deposit. This is precisely shown by the laser profilometer thickness data in figure 9. After the initial release from the cylinder, the grains spread to a maximum width of W ≈ 8.5 cm, which is maintained as the avalanche travels downslope (figure 9a). The wave peak height increases from the initial stages up to h w ≈ 4.3 mm at x L = 30 cm and afterwards remains nearly constant, as in the previous experiment. In this case however, the deposit layer in the avalanche trough also remains approximately constant at h deposit ≈ 1.5 mm, which is, crucially, very slightly thinner than the static layer thickness, h stop ≈ 1.7 mm (figure 9b). The deficit in grains in the trough is made up by the formation of very small static levees, which are slightly higher than the original erodible layer depth h 0 as shown in figure 9(a). Since this flow is at the same angle as the growing avalanche, shown in figure 4, it is clear that it is possible to switch between growing, decaying and steady avalanches by changing h 0 , at least for a certain range of slope angles.
3. Measurements of the rough bed friction law Pouliquen (1999a) performed laboratory experiments for flows of spherical glass beads on a rough chute and found an empirical friction law that was valid for steady uniform flows at various slope angles. This law was derived from a steady uniform flow rule that was later found by Forterre & Pouliquen (2003) to take a slightly different form for flows of sand on a rough chute. An extension of the Pouliquen (1999a) friction law was made by Pouliquen & Forterre (2002) to include flows outside of the dynamic regime, when the Froude number was lower than the minimum required for a steady uniform flow to be possible, down to Fr = 0 for static material. This comprised a law that took a different form in each of these dynamic, static and intermediate regimes. They considered two critical slope inclination angles as functions of the flow thickness h, namely ζ stop (h) and ζ start (h), where ζ stop (h) is the slope angle at which a steady uniform flow leaves a deposit of thickness h and ζ start (h) is the angle at which a layer of thickness h is mobilized. The thickness of a deposit left by a steady uniform flow at an inclination angle ζ is denoted by h stop (ζ ), which is the inverse function of ζ stop (h), and the thickness of a static layer that is mobilised when the inclination is increased to an angle ζ is denoted by h start (ζ ), which is the inverse function of ζ start (h).
The procedures for determining the critical slope angles and the empirical steady uniform flow law have been repeated here for the experimental set-up described above, namely for a flow of 280-350 µm carborundum particles on a rough bed of spherical glass beads. The critical slope angle curves, ζ stop (h) and ζ start (h) are determined by fits to experimental data, which comprise measurements of the respective thicknesses, h stop (ζ ) and h start (ζ ), at various angles ζ with the laser profilometer. Following Pouliquen & Forterre (2002), the experimental fits take the functional form tan ζ stop,start = tan ζ 1,3 + tan ζ 2 − tan ζ 1 1 + h/L . (3.1) Henceforth, the critical thickness curves take the inverse functional form There is no steady-uniform flow for inclination angles ζ < ζ 1 , which is the asymptote of the curve ζ stop (h) for large h, and the flow is accelerated for ζ > ζ 2 = ζ stop (0). The third critical angle ζ 3 is the asymptote of the curve ζ start (h) for large h. The parameter L (having the dimensions of a length) is the characteristic depth of flow over which a transition between the angles ζ 1 and ζ 2 occurs in the friction law and as such it is  dependent on the properties of the grains and on the bed roughness. The fits of these parameters to the experimental data with the functions (3.1) are shown in figure 10. The values of these and all of the other previously obtained or yet to be introduced parameters, whose values remain fixed throughout this paper, are given in table 1.
3.1. Dynamic regime For steady uniform flows the empirical flow relation between the ratio of the flow thickness h to h stop (ζ ) and the Froude number was found by Forterre & Pouliquen (2003) to be Measurements of the flow thickness are made with the laser profilometer, whilst the depth-averaged downslope velocityū, which is equal to |ū|, is determined by the speed of a front travelling down the plane from a space-time plot of overhead images (measured in the same way as the wavespeeds in § 2). Here the constants β = 0.63 and Γ = 0.40 are best fits to the measured data, as shown in figure 11. The values of these constants were previously found to be β = 0.65 and Γ = 0.77 for flows of sand on a rough bed of the same material (Forterre & Pouliquen 2003) and β = 0.136/ √ cos ζ and Γ = 0 for flows of spherical glass beads on a rough bed of the same material (Pouliquen 1999a).
With the critical flow thickness functions (3.1) and empirical flow rule (3.4) known, a friction law for the dynamic regime can now be determined. In a steady uniform flow there is a balance between gravity and the shear stress at the bed, which implies that the force balance may be written as (Pouliquen 1999a) µ = tan ζ . (3.5) By defining the tangent of the critical stopping angle as the friction coefficient for the static layer is found through the steady uniform flow relation (3.5) and the empirical law (3.4) to be Formation of granular levees, troughs and elevated channels

295
However the empirical law (3.4) and therefore the friction law (3.7) are only valid for flows in the dynamic regime, where steady uniform flows are possible. According to Pouliquen & Forterre (2002), a flow is in the dynamic regime if h h stop (ζ ), which is equivalent to Fr β in the case of glass beads, since the Froude number offset Γ = 0. For different granular materials, such as sand and the carborundum particles used here, which have higher values of β, the transition between the dynamic and intermediate regimes is more complicated. Edwards & Gray (2015) found that allowing a dynamic regime for h h stop (ζ ) in the same way as for glass beads, with β = 0.65 and Γ = 0, leads to the formation of erosion-deposition waves that deposit grains in the wave troughs with a layer of thickness h deposit that is significantly lower than h stop (ζ ). This implies that for a steady uniform flow to leave a layer of thickness h stop the transition between dynamic and intermediate friction regimes must occur at higher Froude number. Furthermore, for materials such as sand, where β < Γ , the minimum Froude number for the transition to the intermediate regime, i.e. Fr = β − Γ by (3.4), is negative, implying that all flows can only be in the dynamic regime. To overcome these problems, a friction law transition point h * (ζ ) is introduced at a thickness between h stop (ζ ) and h start (ζ ). This transition thickness is defined as where 0 a 1 is a parameter whose value is chosen (a = 0.5 here) such that numerical simulations of the h stop experimental measurement procedure correctly produce deposit layers with thickness h deposit ≈ h stop . Changing the value of a does not qualitatively change the possible behaviours, but it changes h * and therefore the flow thickness at which the transition between regimes occurs. Defining β * to be the Froude number at h = h * , and substituting the flow rule (3.4) into (3.8) gives the transition Froude number between the dynamic and intermediate regimes as Nevertheless the offset in the friction law will be used in conjunction with the depth-averaged µ(I)-rheology  in this paper. This is justified by the fact that for shallow flows the additional depth-averaged viscous term, where some modification may be necessary, is already much smaller than the leading-order terms in the downslope momentum balance.

Static regime
Defining the tangent of the critical starting angle as µ start (h) = tan(ζ start (h)), (3.10) Pouliquen & Forterre (2002) showed that a knowledge of the functions µ stop (h) and µ start (h) is sufficient to define an empirical friction law for the whole range of possible In general, the static friction may lie below µ start and simply balance the forces in order to prevent motion. For both the inviscid and depth-averaged µ(I)-rheology this is just a balance between effective basal friction, gravity pulling it down the slope and pressure gradients. It follows that µ = |tan ζ i − ∇h|, where i is the unit vector in the x-direction.

Intermediate regime
In the intermediate friction regime when 0 < Fr < β * , the friction coefficient µ is given by a power law extrapolation between the friction laws in the static and dynamic friction regimes as (3.12) The power of the extrapolation chosen by Pouliquen & Forterre (2002) was κ = 10 −3 , in order to give a rapid decrease of the friction coefficient close to Fr = 0. They found their results were insensitive to changes in κ up to a value of 10 −2 . Both of these values are extreme. As a result the decrease in friction from the static value is concentrated in a narrow boundary layer that lies close to Fr = 0, and for much of the range the friction coefficient is essentially constant. In fact, these values of κ are so extreme, that, to the accuracy of machine precision, the boundary layer is essentially invisible in numerical simulations. This may explain why Pouliquen & Forterre (2002) found that when they simulated a finite mass release on an erodible layer it had difficulty coming to rest, and ran out considerably further than in experiments. For the simulations presented here a value of κ = 1 is used instead. This gives stationary material in the hysteretic region of parameter space much greater stability to small disturbances, as well as slowly moving material, above h stop , a greater propensity to come to rest. In general, the value of κ should be at least of order 10 −1 to prevent static layers that are thinner than h start from moving without being disturbed and no greater than order 1, to allow them to flow if forced with the steady uniform flow velocity.
Formation of granular levees, troughs and elevated channels

Governing equations
The erosion-deposition waves observed experimentally in §2 are now modelled using Baker et al.'s (2016a) two-dimensional generalization of Gray & Edwards' (2014) depth-averaged µ(I)-rheology combined with the modification of Pouliquen & Forterre's (2002) friction law discussed in § 3. The theory is therefore a two-dimensional depth-averaged generalization of the erosion-deposition wave model of Edwards & Gray (2015).

Depth-averaged equations with viscous dissipation
The theory does not explicitly resolve the interface between static and flowing layers in the normal direction to the chute or the erosion and deposition rates between them. Instead, the avalanche thickness h and depth-averaged velocityū = (ū,v) are defined over the entire layer. With this representation the depth-averaged mass and momentum balance equations are ∂h ∂t + ∂ ∂x (hū) + ∂ ∂y (hv) = 0, (4.1) where χ = u 2 /ū 2 is the shape factor, S x,y are source terms and V x,y are viscous-like terms. The µ(I)-rheology implies that for steady uniform flow a Bagnold velocity profile develops (see e.g. GDR-MiDi 2004;) and the resulting shape factor χ = 5/4. However, non-unity values of the shape factor change the characteristic structure of the inviscid equations, and causes problems when handling grain-free regions (Hogg & Pritchard 2004), so virtually all avalanche models (e.g. Grigorian et al. 1967;Savage & Hutter 1989;Gray et al. 1999;Pouliquen 1999b;Pouliquen & Forterre 2002;Gray et al. 2003) assume, as we do here, that χ = 1. This is a reasonable assumption since the dominant momentum balance is between source terms and pressure gradients, whilst the acceleration terms are small. Furthermore, 298 A. N. Edwards, S. Viroulet, B. P. Kokelaar and J. M. N. T. Gray Viroulet et al. (2017) showed that variation in the shape factor away from unity has little effect for flows in which the Froude number is less than one. The source terms consist of the downslope component of gravity and the effective basal friction µ between the avalanche and the rough plane, which opposes the direction of motion.  showed that to leading order the inviscid avalanche equations emerge naturally from depth averaging the µ(I)-rheology (GDR-MiDi 2004;Jop et al. 2006) with the basal friction given by the dynamic basal friction law of Pouliquen & Forterre (2002). The depth-averaged µ(I)-rheology Baker et al. 2016a), used here, differs from the standard inviscid equations by the inclusion of in-plane deviatoric stresses, which give rise to depth-averaged viscous-like terms in the momentum balances (4.2)-(4.3). In two dimensions these are respectively, where the coefficient ν in the effective viscosity νh 1/2 /2 is explicitly determined in the depth-integration process and is given by (4.8) with the parameters L, β, ζ 1 and ζ 2 arising from Pouliquen & Forterre's (2002) basal friction law described in § 3. The value of ν for the various slope inclinations is summarized in table 2. Strong evidence for these viscous terms is provided by the fact that, unlike the inviscid avalanche model, they predict a downslope velocity profile across a channel with rough side walls (Baker et al. 2016a). The viscous terms are also vital in regularizing depth-averaged models of segregation-induced fingering (Pouliquen, Delour & Savage 1997;Pouliquen & Vallance 1999;Woodhouse et al. 2012;Baker, Johnson & Gray 2016b). In addition, the functional dependence on the slope inclination ζ and thickness h to the three-halves power is crucial in matching the cutoff frequency of roll waves without any additional fitting parameters . It should be noted, however, that the theory is only well posed for slope inclinations in the range of steady uniform flow, i.e. for ζ ∈ [ζ 1 , ζ 2 ], reflecting the fact that the full µ(I)-rheology is ill posed for high and low inertial numbers (Barker et al. 2015). Outside of this range of angles some form of regularization is therefore required, i.e. the negative viscosity for ζ / ∈ [ζ 1 , ζ 2 ] must be prevented. Edwards & Gray (2015) found that by using Pouliquen & Forterre's (2002) extended basal friction law, which has additional static and intermediate regimes, it was possible to construct steadily travelling erosion-deposition wave solutions with static and mobile regions. In principle, the form of the viscous term should change when the friction law is in the intermediate or static regimes, i.e. for Fr < β * . However,  it is not obvious how to achieve this, since there is no longer a steady uniform flow solution to determine the velocity and pressure profiles in the depth-integration process. Neither is it obvious how the offset Γ for angular particles (e.g. sand and carborundum) can be included in the viscosity. However, since the viscous terms represent a singular perturbation to the equations (i.e. they are the highest-order derivatives, but most of the time they are very small), the simple form in (4.6)-(4.7) with the viscous coefficient (4.8) are assumed to apply to all regimes in this paper.

Numerical method
The standard depth-averaged avalanche equations (e.g. Gray et al. 2003) represent a system of hyperbolic equations that require high-resolution shock capturing numerical methods (e.g. Nessyahu & Tadmor 1990) to solve them. Although our problems are still convection dominated the inclusion of the depth-averaged µ(I)-rheology Baker et al. 2016a) changes the system into a set of convection-diffusion equations. This paper therefore uses the closely related semi-discrete high-resolution non-oscillatory central schemes of Kurganov & Tadmor (2000) together with a second-order Runge-Kutta method, with a step size of 10 −4 s, for their time evolution. The scheme has been tested against the one-dimensional travelling wave solutions of Edwards & Gray (2015) and the value of the time step has been chosen to minimize the small variations in the static layer thickness due to overshoots, while maximizing the efficiency of the computation. The travelling wave results have also been validated by using an alternative third-order Runge-Kutta adaptive step method (Medovikov 1998) for the time evolution, although this method is much more computationally expensive.
In order to solve the system, the depth-averaged equations (4.1)-(4.7) together with the friction law (3.15) are written in conservative vector form as where w = (h, m, n) T is the vector of conservative variables h, m = hū and n = hv, and the subscripts denote partial derivatives. The resulting convection flux functions f and g, and the source term vector S = (0, S x , S y ) T are (4.12) For the problems of interest here the flows are subcritical (Fr < 1) everywhere, which requires two inflow conditions and one outflow condition to be specified for the hyperbolic system (e.g. Weiyan 1992). The introduction of diffusive terms means that the system is now parabolic and an extra boundary condition is required, which is applied here at the outflow end of the domain. Periodic boundary conditions are imposed at the cross-slope boundaries, although there is no flow across them in any of the following simulations.

Numerical simulations on a two-dimensional plane
The two-dimensional system of conservative equations (4.9)-(4.12), with the friction law (3.15), are solved numerically for the experimental configurations described in § 2. A rectangular computational domain is defined in the range −5 x 95 cm, −10 cm y 10 cm and is discretized with 1000 grid points in the x-direction and 100 grid points in the y-direction. The initial conditions are of a cylindrical mass of radius R = 1.4 cm and height h c = 1.6 cm centred at (x, y) = (0, 0) on top of a layer of thickness h = h 0 , with m = n = 0 everywhere. This will be referred to as the standard release. As the column of grains collapses it sets the previously static particles immediately downstream of the cylinder into motion to form an avalanche. At the downstream boundary x = 0.95 m free outflow is imposed via a linear extrapolation of the values of h and m from the final two columns of interior cells. At the upstream boundary x = −0.05 m, a stationary inflow is imposed with h = h 0 and m = 0.

Growing avalanches
Simulations are performed to model the growing avalanche shown in figure 4. The slope angle is set to ζ = 35.2 • and a deep deposit layer of thickness h 0 = 2.2 mm is prescribed, which is equivalent to h stop (34.4 • ) = 2.2 mm on a shallower slope. The results are plotted as thickness contours in the (x, y)-plane at various times in figure 12. The region that is in motion is shaded grey and the dashed contour at 2.1 mm approximately delineates the outer boundary of the static trough deposit that forms behind the avalanche, which is thinner than the undisturbed erodible layer. A movie is available in the online supplementary material (movie 4). As the avalanche propagates downslope it erodes the thick layer in front of it and progressively grows in size by depositing less material behind it and by broadening laterally. By x L = 30 cm it is approximately 14 cm wide, which is roughly 25 % wider than the experimental avalanche, but it has qualitatively the correct behaviour. As a result, the maximal extent of the trough has a triangular shape, which is consistent with the experimental results of . As the avalanche bulks up, the wave accelerates slightly, although over the length of the chute its wavespeed remains close to its mean value of u w = 0.31 m s −1 , which is slightly greater than the experimental value of 0.26 m s −1 . In order to compare the numerical simulations with the experimental results shown in figure 5, the flow thickness contours in the (t, y)-plane and the temporal evolution of the flow thickness on the centreline are plotted in figure 13. The temporal evolution of the downslope velocity on the centreline and the final static layer depth across the plane at x L = 30 cm are also shown. The plots in figure 13(b,c) show that the avalanche develops into a classical erosion-deposition wave soon after the initial collapse, with the highest velocities occurring during the peak of the wave, which has a sharp front and a more gradual decay of both height and velocity in the tail. Outside of the wave the velocity is zero. The avalanche peak height at x L = 50 cm is h w = 5.8 mm, which is only slightly less than the experimental value of approximately 5.9 mm. As the wave propagates downslope it erodes more grains from the front than it deposits in the tail, and the incised trough becomes progressively deeper further downstream, reaching h deposit = 1.6 mm at x L = 50 cm. This is slightly lower than the experimental value of 1.7 ± 0.2 mm. The cross-section in figure 13(d) shows that after the avalanche has passed by there is a strongly incised trough, which has a base that is 0.6 mm lower than the initial deposit thickness. At the sides are some very small levees, but they are not big enough to balance the mass that has been lost, so that this channel is typical of a growing avalanche. Overall the model produces a slightly wider flowing region, with a shorter thinner avalanche than is observed experimentally. However, the important features of the experimental flow are captured in the simulations and the results are qualitatively comparable, as they produce an avalanche that is continually growing in width and in wave peak height.
There are several possible reasons for the apparently increased width of the numerical simulations as compared to the experiments. One possibility is that the friction law (3.15) is still not perfect, and that by modifying it one might be able to get closer to the experiments. Another possibility is that the relatively crude representation of the deposition mechanism is at fault, and that to get better agreement the interface between flowing and non-flowing material really does have to be resolved in the vertical. Unfortunately there is no accepted way of doing this yet. A final possibility is that the initial collapse from the cylinder may not be captured well in the depth-averaged framework. Certainly a few of the grains are stranded at the experimental release point and do not, as assumed, propagate down the plane. This would reduce the mass that should be used in the simulations and hence also produce a width reduction. Better agreement between theory and experiment could also potentially be achieved by commencing the simulations once the avalanche is shallow. However, this has not been done here, because it would require us to prescribe the unknown thickness and velocity fields at this later time.

Decaying avalanches
With the slope at ζ = 34.1 • and the static layer thickness at h 0 = 2.5 mm, equal to the theoretical h stop for 34.1 • , the standard release again generates a three-dimensional erosion-deposition wave that propagates downslope. This time, however, the contour plots in figure 14 show that there is net deposition and by t = 2.5 s the avalanche has completely stopped. The shape of the deposit and the final stopping position of the head, at x = 53 cm, are very close to the experiments shown in figure 6, where  the maximum distance reached by the avalanche is 55 cm. The space-time plots in figure 15(a,b) indicate that the avalanche is again wider than in the experiments in figure 7, but that the model otherwise accurately captures the deposition of mass. At x L = 30 cm the deposit thickness in the centre of the channel h deposit = 2.8 mm, which is very close to the experimental value of 2.7 mm. The cross-section in figure 15(d) shows that the avalanche is leaving behind an elevated channel, which is everywhere higher than the undisturbed deposit and which is flanked by small levees. This cross-slope profile is typical of a depositing flow, and as one progresses downstream the elevated channel narrows. During the flowing phase the avalanche has an almost constant wavespeed u w = 0.15 m s −1 , which is moderately larger than the experimental value of approximately 0.13 m s −1 . However, figure 15(b,c) shows that as the avalanche propagates downslope both the peak thickness and the magnitude of the velocities decrease. At x L = 30 cm the peak height h w = 4.8 mm, which is slightly less than the experiments where the thickness is approximately 5.0 mm. This suggests that the model is again allowing for a wider spread of material into a shorter and thinner avalanche than is observed experimentally, but that all the main attributes are qualitatively captured and the quantitative comparison is still good.

Steady avalanches
Perhaps most interestingly, with the slope at ζ = 35.2 • and the static layer thickness h 0 = 1.7 mm (equal to h stop (35.2 • ) and the experimental value of h 0 ≈ 1.7 mm) after the initial collapse the avalanche appears to approach a steadily travelling wave as shown in figure 16. In order for this to occur there has to be a perfect balance between the erosion of material at the front of the avalanche and the deposition that occurs behind it. The wavespeed is calculated by measuring the distance travelled by the wave peak in the time after the initial spreading and before it leaves the computational domain, and is found here to be u w = 0.21 m s −1 , which is only slightly greater than the value of 0.20 m s −1 found in experiments. The space-time plots in figure 17(a-c) indicate that after the initial collapse the avalanche quickly attains a constant width W ≈ 11 cm, which is 30 % wider than the experiments in figure 8. The avalanche peak height at x L = 30 cm is h w = 4.0 mm, which is slightly less than the experimental value of approximately 4.3 mm. As the avalanche moves farther downstream its peak height and maximum velocity continue to evolve, albeit very slowly. Behind the wave, the deposit in the centre of the trough approaches a constant thickness h deposit = 1.6 mm, which is slightly lower than the initial layer thickness, and close to the experimental value of 1.5 mm. The cross-slope profile in figure 17(d) shows that levees form on either flank of the shallow trough. The form of this channel is therefore typical of a steady travelling erosion-deposition wave, in which the initial erodible layer ahead of the wave is reworked into a pair of levees and a trough behind it, whose cross-sectional area exactly equates to that in the undisturbed deposit.

Numerical simulations of the avalanche in a travelling frame
The survival of the steadily propagating avalanche in § 5.3 suggests that there is a three-dimensional travelling wave solution to the governing equations. This will now be investigated further by solving the problem for much longer times in a frame moving with the wave. 6.1. Governing equations in a travelling frame of reference In a frame moving with constant speed u w the transformation of coordinates Formation of granular levees, troughs and elevated channels implies that the depth-averaged mass and momentum conservation laws (4.1)-(4.3) transform to where S ξ = S x and S η = S y are given by (4.4) and (4.5) respectively, and the diffusive terms are In order to solve the travelling wave equations (6.2)-(6.6) numerically, they are written in terms of the conserved variables w = (h, m, n) T = (h, hū, hv) T as The resulting convection flux functions arê with g given in (4.10). The source term vector is S = (0, S ξ , S η ) T = (0, S x , S η ) T and the diffusive flux functions are given bŷ (6.10) 6.2. Numerical simulation The conservative form of the governing equations (6.7)-(6.10) with the source terms (4.4)-(4.5) are solved numerically on a rectangular computational domain (0 ξ 50 cm, −10 η 10 cm), which is discretized over 500 grid points in the ξ -direction and 100 grid points in the η-direction. The initial conditions are taken from the results of the numerical simulation of § 5.3 in the region 45 x 95 cm downstream at t = 3.7 s (figure 16e). Free outflow is imposed at the upstream boundary ξ = 0 m via a linear extrapolation of the values of h and m from the first two columns of interior cells. At the downstream boundary ξ = 0.5 m, a stationary inflow is imposed with h = h 0 and m = 0, such that the initial avalanche effectively has an infinitely long static layer of thickness h 0 ahead of it. The slope angle ζ = 35.2 • and h 0 = 1.7 mm. Based on the measured wavespeed in the previous simulations in § 5.3 the speed of the travelling frame is set to u w = 0.21 m s −1 .
The initial condition and the computed state after 20 s are shown in figure 18(a,b). A movie is also available online (movie 7). These show that the initial and final states are virtually unchanged. The front position has moved slightly to the left, which simply implies that the imposed speed of the travelling frame was slightly greater than the actual wavespeed of the avalanche. However, there is nothing to suggest that after a longer period of computation, and the exact wavespeed of the avalanche being known, that it would not continue to propagate unchanged in this domain forever. It is rather remarkable that the presence of an erodible layer of grains can change the apparent friction of the slope and can thus produce indefinitely long run out. This is potentially very important for quantifying the risks posed by natural hazards when there are slopes covered with easily erodible material, such as in the case of snow avalanches.
The steady travelling wave solution is the three-dimensional equivalent of the two-dimensional erosion-deposition waves observed by Edwards & Gray (2015). The final computed state shown in figure 18(b) compares very favourably with the experimentally observed deposit contour shown in red. The avalanche is wider and shorter than the experimental avalanche, but the quantitative comparison is reasonably good and the qualitative features of a levee and shallow trough are very close to what are observed in the simulations. One interesting feature in the experiments that one does not see in the numerical computations is that the tail of the moving part of the avalanche has a pronounced 'V'-shape. This may be an indication that the avalanche is eroding down to different depths within the channel. To help visualize the computed wave figure 19 shows a three-dimensional surface plot of an erosion-deposition wave with an exaggerated vertical height and colour scale used to reveal the shape of the wave (in red), the trough (in blue), the levees (in yellow) and the undisturbed deposit (in green).

Conclusions
This paper shows that when a finite mass of carborundum particles is released on a rough plane covered with an erodible layer of the same grains, there are three qualitatively different behaviours that can be observed dependent on the slope inclination and the deposit thickness. For sufficiently shallow slopes the avalanche propagates for a finite distance before coming to rest, as one may intuitively expect (Pouliquen & Forterre 2002). At a steeper slope angle, the avalanche has been observed to either propagate steadily with constant wavespeed and shape, or, when released on a thicker erodible layer, to accumulate additional particles (bulking) and therefore grow and accelerate. The steady travelling avalanche is the three-dimensional equivalent of the two-dimensional erosion-deposition waves observed by Edwards & Gray (2015). Provided that both the static layer in front of the flow and the inclination angle do not change, this wave can travel arbitrarily long distances, which is an intriguing observation. Behind it, the wave leaves a shallow deposited trough that is flanked by levees on either side as illustrated in figure 19. Since erosion and deposition are in exact balance, the total volume of material in the levees and the trough is the same as the equivalent volume in the undisturbed layer ahead of the wave. For the growing avalanches the avalanche essentially leaves behind an eroded trough of previously mobilized material that widens with increasing downstream distance . Conversely, for decaying avalanches, there is net mass loss, and a channel is left behind that is elevated above the initially static deposit height and narrows with downstream distance, terminating in a static deposited head. Given the notorious difficulty in modelling erosion-deposition problems, it is remarkable that all three of these behaviours can be modelled by combining the extended rough bed friction law of Pouliquen & Forterre (2002) with a twodimensional generalization (Baker et al. 2016a) of the depth-averaged µ(I)-rheology of . The theory is crude, in that the interface between the static and moving grains is not resolved in the vertical direction. Instead, the flow is either static or moving through its entire depth, which is a good approximation for the shallow erodible beds considered here. An important element of the model is that the static, dynamic and intermediate regimes in the rough bed friction law (3.15) have been modified slightly to give greater stability to static layers that form in the hysteretic region between h stop and h start . In particular, the power law transition starts at a higher Froude number β than is normal and the subsequent transition between the static and dynamic states is much more gradual (κ = 1) than in the original formulation (Pouliquen & Forterre 2002) where κ was equal to 10 −3 . The resulting model is able to capture qualitatively all of the key behaviour observed in the experiments, and the quantitative agreement is also good, even though the avalanches tends to be shorter and wider than observed. The growing, decaying and steady avalanches will stop if they flow out onto a run-out plane. The same equations can be used to compute the arrest, however, care must be taken to regularize the coefficient ν in the viscosity for angles ζ < ζ 1 .
Both the bulking up of snow slab avalanches through frontal ploughing (Sovilla et al. 2001;Sovilla & Bartelt 2002;Sovilla et al. 2006) as well as their sudden arrest in the run-out zone are strongly affected by the net erosion and deposition in the flow. In particular, the existence of an easily erodible layer confers much greater mobility to the flow than one might otherwise expect (Mangeney et al. 2010;Iverson et al. 2011) and run-out distances can be significantly enhanced. It is therefore anticipated that the relatively simple model for erosion and deposition, developed in this paper, may find important application in predicting the risk posed by a wide range of geophysical mass flows on erodible slopes. A key challenge for the future is how to determine the h stop and h start curves for real geophysical materials where the fixed bed may not be clearly defined.