Multiple solutions for granular flow over a smooth two-dimensional bump

Geophysical granular flows, such as avalanches, debris flows, lahars and pyroclastic flows, are always strongly influenced by the basal topography that they flow over. In particular, localised bumps or obstacles can generate rapid changes in the flow thickness and velocity, or shock waves, which dissipate significant amounts of energy. Understanding how a granular material is affected by the underlying topography is therefore crucial for hazard mitigation purposes, for example to improve the design of deflecting or catching dams for snow avalanches. Moreover, the interactions with solid boundaries can also have important applications in industrial processes. In this paper, small-scale experiments are performed to investigate the flow of a granular avalanche over a two-dimensional smooth symmetrical bump. The experiments show that, depending on the initial conditions, two different steady-state regimes can be observed: either the formation of a detached jet downstream of the bump, or a shock upstream of it. The transition between the two cases can be controlled by adding varying amounts of erodible particles in front of the obstacle. A depth-averaged terrain-following avalanche theory that is formulated in curvilinear coordinates is used to model the system. The results show good agreement with the experiments for both regimes. For the case of a shock, time-dependent numerical simulations of the full system show the evolution to the equilibrium state, as well as the deposition of particles upstream of the bump when the inflow ceases. The terrain-following theory is compared to a standard depth-averaged avalanche model in an aligned Cartesian coordinate system. For this very sensitive problem, it is shown that the steady-shock regime is captured significantly better by the terrain-following avalanche model, and that the standard theory is unable to predict the take-off point of the jet. To retain the practical simplicity of using Cartesian coordinates, but have the improved predictive power of the terrain-following model, a coordinate mapping is used to transform the terrain-following equations from curvilinear to Cartesian coordinates. The terrain-following model, in Cartesian coordinates, makes identical predictions to the original curvilinear formulation, but is much simpler to implement.


78
S. Viroulet and others

Introduction
Granular free-surface flows are encountered in many geophysical and industrial processes. They occur at a wide range of scales from table-top flows in a rotating drum (Gray 2001) to geophysical mass flows, such as avalanches (Savage & Hutter 1989), debris flows (Iverson & Denlinger 2001;Johnson et al. 2012) and pyroclastic flows (Branney & Kokelaar 1992;Mangeney et al. 2007). Despite differences in scale of several orders of magnitude, all of these flows may be considered shallow, after the initial release, with typical flow thicknesses being much smaller than in-plane length scales, and they can therefore be modelled with a relatively simple depth-averaged avalanche theory. Such shallow-water-type equations were first used to model snow avalanches (Grigorian, Eglit & Iakimov 1967), before Savage & Hutter (1989) presented one of the earliest formal derivations. They used a Mohr-Coulomb internal rheology and a constant Coulomb basal friction law, which introduced additional source terms into the momentum equation as well as an earth pressure coefficient multiplying the depth-averaged pressure, that switched between 'active' and 'passive' values dependent on whether the flow was dilatational or compressional. These equations were then applied to the motion of a finite mass of granular material on a slope with variable topography by using a slope fitted curvilinear coordinate system (Savage & Hutter 1991;Greve & Hutter 1993). Since their seminal work, this theory has been used in many different situations and extended to two-dimensional flows (downslope and transverse) over complex topography Bouchut & Westdickenberg 2004;Doyle, Hogg & Mader 2011). Gray, Tai & Noelle (2003) simplified the Savage & Hutter (1989) model to a conservative shallow-water-like structure with source terms, by assuming that the in-plane deviatoric stresses were negligible and hence that the effective earth pressure coefficient was equal to unity. This was significant because the resulting hyperbolic system of equations admitted discontinuous solutions with abrupt changes in material thicknesses and velocities, features that drastically change the overall properties of the flow and represent a key challenge in natural hazard mitigation, as well as optimisation of industrial processes. Such 'shock waves' have been extensively studied for a single layer of Newtonian fluid, whether it be direct applications to geophysical events such as tidal bores (Stoker 1949;Chanson 2009) or underwater dune formation (Fourriere, Claudin & Andreotti 2010;Andreotti et al. 2012), or more fundamental studies into the stability of hydraulic jumps in the presence of a single obstacle (Lawrence 1987;Baines & Whitehead 2003;Defina & Susin 2003) or in a wavy channel (Wierschem & Aksel 2004). High-speed granular free-surface flows exhibit shock waves (Gray et al. 2003) that are very similar to those observed in shallow water (Rouse 1938;Ippen 1949), gas dynamics (Ames Research & Staff 1953) and flows in collapsible tubes (Shapiro 1977). These shocks may be generated in a number of different ways, for example using an obstacle (Hakonardottir et al. 2003;Hakonardottir & Hogg 2005;Cui & Gray 2013) or a contraction (Aker & Bokhove 2008) in the flow. Studies into granular flows past obstacles in particular (e.g. Tai et al. 1999;Gray et al. 2003;Cui, Gray & Johannesson 2007;Johannesson et al. 2009) provide useful insight for the design of deflecting or catching dams for geophysical hazards, where understanding the shock structure is crucial. A recent review article on the applications of the theory to geophysical flows and small-scale experiments is given in Delannay et al. (2017).
Incorporating the effects of topography as well as the presence of erodible material, which may cause the flow to bulk up, remain key challenges in developing accurate Multiple solutions for granular flow over a bump 79 mathematical models of geophysical mass flows. Recent studies have shown that the topography, in particular centrifugal forces related to curvature effects, have a major impact on the seismic signals that are generated by an avalanche (Favreau et al. 2010;Levy et al. 2015). These signals can be used to determine the flow history and properties such as its velocity and rheology for example (Brodsky, Evgenii Gordeev & Kanamori 2003). Since understanding and quantifying the rheology of a geophysical flow remains one of the main difficulties in predicting their behaviour, it is of crucial importance to accurately quantify the topography effects in those inverted seismic signals which represent an interesting method for direct measurement of geophysical events (Farin, Mangeney & Roche 2014;Levy et al. 2015). Another key aspect for geophysical flows is the existence of erodible material in the path of the flow, which may often have been left by a previous event. Laboratory experiments have shown that the presence of an erodible granular substrate can drastically change the behaviour of the flow by enhancing the runout distance Farin et al. 2014) or giving rise to erosion-deposition waves (Edwards & Gray 2015) that could have an important effect on the amount of material transported.
In some flow configurations the presence of obstacles can lead to grain-free regions, as well as shocks, when material flows around an obstruction (Gray et al. 2003;Cui & Gray 2013), and it is still possible to capture this behaviour within the depth-averaged framework of Gray et al. (2003). In other cases, a build-up of material on the upstream side of a barrier may eventually lead to overtopping (Faug, Beguin & Chanut 2009;Chanut, Faug & Naaim 2010) and, if the oncoming velocities are sufficiently high, particles losing contact with the base and becoming airborne to form a granular jet (Hakonardottir et al. 2003;Faug 2015). Beyond the take-off point, the depth-averaged governing equations are no longer directly applicable, but Hakonardottir et al. (2003) showed that individual particle trajectories may be accurately predicted using a ballistic approach.
This paper investigates the dynamics of dense granular flows down an inclined chute, with a smooth bump acting as an obstacle to the flow. Interestingly, two contrasting steady-state regimes are found to coexist, for the same chute angle and mass flux, with only the initial distribution of grains on the chute determining which is selected. If the chute is initially empty, the flow leaving the hopper is able to accelerate sufficiently so that it takes off and forms a granular jet when moving over the bump. However, placing erodible particles in front of the obstacle reduces the energy of the oncoming flow and can generate a shock upstream of the bump. A depth-averaged avalanche theory (Savage & Hutter 1991;Gray et al. 1999) in a curvilinear coordinate system that follows the topography of the bump is able to accurately predict both the take-off point of the jet (which may then subsequently be modelled as a series of ballistic trajectories) as well as the shock position and thickness profile at steady state. Time-dependent numerical solutions of this terrain-following avalanche model are used to show the transient dynamics before the system reaches equilibrium, as well as the deposition of a static deposit on the bump once the inflow ceases. A standard depth-averaged avalanche theory in an aligned Cartesian coordinate system (Savage & Hutter 1989;Gray et al. 2003) is also presented and compared to the curvilinear terrain-following model for the shock at steady state. For this highly sensitive problem, this demonstrates that it is important to take into account the local changes in the slope inclination and the basal curvature in the depth-averaged momentum balance (Savage & Hutter 1991;Gray et al. 1999), rather than using a fixed inclination and topography height gradients to describe the effect of flowing over the bump. For many geophysical applications, 80 S. Viroulet and others however, it is not convenient to set-up terrain-following curvilinear coordinates. In § 7.2 it is therefore shown how to use a coordinate transformation to convert the curvilinear terrain-following model of Savage & Hutter (1991) and Gray et al. (1999) into Cartesian coordinates, i.e. it is explicitly shown how to produce a Cartesian terrain-following avalanche model.

Small-scale experiments
The experimental set-up consists of a 1.8 m long smooth Perspex chute, inclined at an angle θ to the horizontal, with rigid side walls 5 cm apart. The base of the chute incorporates a smooth bump extending across the width of the channel, described by a hyperbolic secant profile with a maximum height of 4.75 cm and with its centre located 43 cm downstream of the inflow. The bump has a characteristic length scale of 4 cm, which implies that ∼90 % of its amplitude change occurs over a downstream distance of 12 cm and its wavelength is therefore approximately 24 cm. The experiments are performed with spherical glass beads of diameter 600-800 µm, which are coloured blue and white to aid visualisation. Grains are loaded into a hopper at the top of the channel and released from rest using a double gate system consisting of a fixed gate to control the initial flow thickness h 0 and another movable barrier to control the release time. For all of the experiments presented here the gate height remains constant at h 0 = 1.5 cm; qualitatively similar behaviour has been observed when different values are used.
Two different types of initial condition are implemented in experiments. In the first, the chute is cleared of all downstream particles before the gate is opened, so that the granular material flows down a smooth, empty channel. In the second, static particles (of the same type) are placed slightly upstream of the bump, and the oncoming flow from the hopper then travels over a partially erodible bed. These different initial conditions evolve to two dramatically different stable steady-state regimes, which shall be referred to as the 'jet' and the 'shock', respectively. Figures 1 and 2 show the time evolution of these two types of flows, and figure 3 shows the final steady-state behaviour. Supplementary movies of these flows are available online at https://doi.org/10.1017/jfm.2017.41.
An initially empty chute usually leads to the formation of a jet. As soon as the gate is opened, the grains flow out of the hopper and accelerate downstream. For slope angles θ 35 • , they reach a sufficiently high velocity to detach from the base and become airborne as they flow over the bump. Once they have passed this take-off point, the grains approximately follow ballistic trajectories, before landing at a point downstream of the bump (see figure 1 and movie 1).
A jet may still form when only a small mass of particles is placed in front of the obstacle. In this case, the oncoming material has enough momentum to entrain the erodible bed into the bulk flow, which then takes off as before. However, adding more static particles can lead to the formation of a steady shock upstream of the bump (see figure 2 and movie 2). When sufficient mass is added, the erodible layer provides enough resistance to drastically decrease the bulk velocity of the accelerating flow from the inlet, without itself being pushed over the top of the bump. To conserve mass, the flow thickness must consequently increase, and this sharp transition in flow height and velocity is referred to as a granular shock. Some material may also be scattered into the air during initial impact with the stationary material, as seen in figure 2. The shock propagates upstream until it reaches a steady-state position. It is also possible to generate shocks using alternative dissipative mechanisms. One such Multiple solutions for granular flow over a bump 81 FIGURE 1. Snapshots of an experiment showing the time evolution of the jet to the steady state. As the oncoming material flows over the top of the bump it is able to detach from the base and follow ballistic trajectories, before landing and coming into contact with the chute once again. The experiment is performed at a constant slope angle θ = 39 • with pictures taken at approximate times t = 0.3; 0.6; 0.9 and 4.0 s. Note that the images have been slightly rotated to maximise space. The bump height of 4.75 cm acts as a scale. The time-dependent evolution is shown in supplementary movie 1, which is available online.

82
S. Viroulet and others FIGURE 2. Snapshots of an experiment showing the time-dependent evolution of the shock towards steady state. As the oncoming material from the inflow collides with the layer of static particles upstream of the bump there is a sharp decrease in bulk velocity and associated increase in flow thickness. This shock propagates upstream until it reaches an equilibrium position. The experiment is performed at a constant slope angle θ = 39 • with pictures taken at times t = 0; 0.4; 0.7; 1.0; 1.5 and 4.0 s. Note that the images have been slightly rotated to maximise space. The bump height of 4.75 cm acts as a scale. A supplementary movie is available online.
Multiple solutions for granular flow over a bump 83 FIGURE 3. Experimental results showing the two different steady states possible for the same slope angle and different initial conditions. Both experiments have been performed at θ = 40 • with the same inflow conditions. Note that the images have been slightly rotated to maximise space. The bump height of 4.75 cm acts as a scale. approach is demonstrated in supplementary movie 3, where an initially empty chute leads to the formation of a jet, and a rigid obstacle is temporarily placed into the path of the flow. This again leads to a shock that eventually settles down to an equilibrium position. Movie 3 also shows that the shock position remains stable to perturbations in the flow. After reaching steady state, the flow is again obstructed downstream of the shock. This momentarily causes the shock to migrate upstream, but as soon as the obstacle is removed the shock relaxes back to its steady-state position. Similarly, sweeping away small amounts of the slower moving material in the shock causes it to temporarily move downstream before returning to its original position. However, sweeping away a larger proportion of the slowly moving grains can lead to complete remobilisation and transition back to the jet regime.
When an initial deposit of static particles is used to generate a steady shock, a critical mass of stationary material is required to sufficiently reduce the momentum of the flowing grains. This critical mass depends on the inclination angle of the chute. Several experiments have been performed with varying slope angles and masses of erodible particles to determine the necessary conditions for the formation of a steady shock. The results are summarised in figure 4, which shows a phase diagram indicating whether a shock or a jet is formed, depending on the slope angle and mass of static grains. As expected, more particles are needed to generate a shock Also shown are the extreme slope angle regions, where it is not possible to keep any particles at rest (high slope angles) or a spontaneous shock forms and propagates upstream to the gate (low angles).
when the slope angle is higher. For slope angles of 34 • or lower, the flow never becomes fast enough to pass over the bump, and a shock is spontaneously generated even when there are no static particles. However, the shock does not reach a steady state and keeps propagating upwards until it reaches the gate and the flow stops completely. Contrastingly, for steep slopes θ 41 • , the friction on the smooth base is not sufficient to keep any particles placed in front of the obstacle at rest. They roll over the bump and a jet always forms. Assuming there are enough particles in front of the bump to create a shock in the first place, the position of the shock does not depend on the initial mass. Note that the phase diagram of figure 4 is specific to this precise experimental set-up. In general, the critical mass and maximum and minimum slope angles for steady-shock formation will depend on the geometry and position of the bump, the inflow gate height and the frictional properties of the chute and specific granular material used. However, it is expected that qualitatively similar behaviour will be found for all configurations.

Depth-averaged terrain-following avalanche theory
The flow is modelled using a depth-averaged avalanche theory that is based on the work of Savage & Hutter (1991) and Gray et al. (1999), and which uses a curvilinear coordinate system to follow the basal topography. This is referred to throughout this paper as the terrain-following avalanche model. In order to construct the curvilinear coordinates, an inclined Cartesian coordinate system OXZ is first defined, with the unit vector i aligned with the downslope direction of the chute, which lies at a constant angle θ to the horizontal. The unit vector k is perpendicular to the chute and points Multiple solutions for granular flow over a bump 85 upwards. The shape of the rigid topography on which the avalanche flows is then described by where the coefficients are all in metres. This is illustrated in figure 5(a) and forms a reference surface for the terrain-following curvilinear coordinate system Oxz, where the position vector of a point is given by the distance x along this reference surface and the height z in the direction normal to the base (see Gray et al. 1999, and figure 6 for details). For a small increment in the down chute coordinate dX, there is a small increment in the topography height db and by Pythagoras' theorem the associated increment in the curvilinear coordinate dx satisfies dx 2 = dX 2 + db 2 . In differential form this implies which can be integrated to give a mapping between the downstream distance in curvilinear and Cartesian coordinates, i.e.
The curvilinear coordinate is therefore just the arc length of the basal topography. As a result, the curvilinear distance is slightly longer than the Cartesian distance as shown in figure 5(b). A key advantage of the terrain-following model (Savage & Hutter 1991;Gray et al. 1999) is that, as the avalanche flows over the topography, it experiences changes in the local slope inclination ζ with the downstream coordinate x. The angle is computed from the definition of the topography (3.1) by where θ is the angle of the Cartesian coordinates to the horizontal. If the chute is inclined at an angle of θ = 39 • then the local slope inclination ζ ranges between approximately 8 • and 70 • as shown in figure 5(c). The local variation of the inclination angle ζ has a very important influence on the net balance between the downslope acceleration and the resistance due to basal friction in the terrain-following model. The final element of the curvilinear coordinate system is that it also takes account of changes in curvature κ, which is defined as where the factor (3.6) Note that the final result in (3.5) follows from differentiating arctan. Once the local inclination ζ and the curvature κ have been computed from the prescribed topography (3.1), the mapping ( . This forms a reference surface for the slope-fitted curvilinear coordinate system Oxz, with the x-axis following the reference surface and being locally inclined at an angle ζ (x) to the horizontal. The z-axis is in the direction of the local normal, which is at an angle ζ to the gravity acceleration vector.
that the curvilinear coordinates introduce a singularity corresponding to where adjacent z-axes intersect , which occurs at a height z = 1/κ. For the flow over a smooth bump considered here, κ is never large enough that the thickness of the avalanche is greater than the local radius of curvature 1/κ of the topography, so there is no overlap of coordinates. Following Gray et al. (1999) the depth-averaged mass and momentum balance equations for the terrain-following avalanche model in curvilinear coordinates are (3.8) respectively, where g represents the gravitational acceleration, h is the flow thickness measured normal to the base,ū is the depth-averaged velocity in the direction of the reference surface and µ is the effective friction coefficient, which may, in general, depend on the flow variables. Note that the system (3.7)-(3.8) implicitly assumes that 88 S. Viroulet and others the density is constant and uniform, i.e. the flow is incompressible. In addition, the earth pressure coefficient (Savage & Hutter 1989) is assumed to be unity, based on the scaling analysis of Gray & Edwards (2014). The shape factor arises in the depth integration of the momentum transport equation to account for the fact that the product of the depth-averaged velocity is not, in general, equal to the average of the velocity product. Particle image velocimetry (PIV) measurements of the accelerating region through the side walls suggest that the vertical velocity may be approximated by a linear profile of the form u =ū(α + 2(1 − α)z/h), with a shear parameter α ≈ 0.87 representing a large degree of basal slip. This gives a corresponding shape factor χ = 1.006. Although the presence of side wall friction may influence these measurements (Baker, Barker & Gray 2016) they are consistent with surface and basal velocity measurements performed by Faug et al. (2015) in the centre of their channel, which suggest χ = 1.009 if a Bagnold profile with basal slip is used to calculate χ. Similarly small deviations of χ away from unity are also obtained in the discrete element simulations of Brodu, Richard & Delannay (2013). For this reason the velocity profile is approximated here by plug flow, which corresponds to a value χ = 1. The curvature term κχhū 2 on the right-hand side of (3.8) is a key feature of the depth-averaged terrain-following model, and describes the enhancement or reduction of the basal pressure due to centrifugal forces as the avalanche flows over topography, in the same way as one is pushed down harder or thrown up into the air on a roller coaster.
With the assumption that χ = 1, the characteristics of the governing equations (3.7) and (3.8) are given in (x, t) space by where the characteristic wave speeds are λ ± =ū ± gh cos ζ .
( 3.11) This motivates the definition of the granular Froude number as the ratio of the flow speed to the speed of surface gravity waves, and it can be seen that for supercritical flow (Fr > 1) both characteristics move in the same direction, whereas in the subcritical regime (Fr < 1) they travel in opposite directions.
3.1. Steady solution As observed during the experiments, the flow can evolve into a steady state, and in this steady state the depth-averaged mass and momentum balance equations, (3.7) and (3.8), reduce to The mass balance equation (3.13) can be integrated directly, subject to the boundary condition that the avalanche velocity and thickness at the inflow areū 0 and h 0 , respectively, to give a relation between the flow thickness and the velocity By expanding out the derivatives, noting that the slope angle ζ is x-dependent in this curvilinear system, and dividing equation (3.14) by h, the depth-averaged momentum balance becomes where (3.13) has been used to simplify (3.16) and it has been assumed thatū > 0. Substituting κ = −∂ζ /∂x for the curvature and using (3.15) to replace the depthaveraged velocity allows (3.16) to be written as a single ordinary differential equation (ODE) governing the evolution of the flow thickness, For a given basal friction coefficient µ and inflow conditions (h 0 ,ū 0 ), equation (3.17) can be integrated numerically subject to h = h 0 at x = 0. This gives the thickness profile h(x) as material leaves the hopper, which is identical for both the jet and shock regimes. The depth-averaged velocity profile can then be recovered from the mass balance equation (3.15). In the experiments material accelerates and thins as it comes out of the inflow gate, and this behaviour is recovered in the numerical solutions if dh/dx 0 at x = 0. The curvature terms in (3.17) are negligible at the inflow, and since tan ζ > µ for accelerating flows the numerator is always positive at x = 0. Examining the denominator, it can be seen that (3.17) only admits accelerating solutions if Fr 0 1, where is the inflow Froude number and ζ 0 is the inclination angle at x = 0. For the rest of this paper it will be assumed that the Froude number at the inflow equals unity, i.e. Fr 0 = 1, which corresponds to an infinite gradient at the inflow. This is consistent with the notion that particles are fully mobilised along the inside wall of the hopper and are free to move downwards, i.e. there is no dead zone adjacent to the hopper wall. As these particles exit the hopper they are therefore moving down normal to the chute, before they are accelerated downstream. Note that it is possible to choose other values Fr 0 > 1 with less steep thickness profiles (see e.g. Cui & Gray 2013), but further justification for setting Fr 0 = 1 will be given in § 3.2. With this choice there is a singularity in (3.17) at the origin. A power series expansion can be used to integrate directly from the inflow position. To achieve this, the curvature terms are neglected and the ODE is written in the alternative form where the function

S. Viroulet and others
satisfies f (h 0 ) = 0. Expanding about this point by assuming allows the leading-order behaviour of the ODE (3.19) near the inflow to be written .
( 3.22) Here, the relationshipū 0 = √ gh 0 cos ζ 0 has been used to simplify the derivative f (h 0 ), which is a positive constant. Equation (3.22) can be solved, subject to h = 0 at x = 0, to give . (3.23) Choosing the negative root to ensure a thinning, accelerating flow, the thickness is calculated near the hopper by the algebraic expression and for positions further downstream it is computed numerically using the full ODE (3.17). In the complete absence of curvature terms, and a constant basal friction coefficient µ, it is also possible to construct an exact solution by integrating (3.17) directly and solving a cubic equation forū (see Cui & Gray 2013, p. 320).
3.2. Basal and wall friction For a given basal topography, the only parameter left in (3.7), (3.8) is the effective friction µ. An efficient way to determine this coefficient on a smooth slope is to measure the surface velocity of the accelerating flow and match it to the one-dimensional theory (see Cui & Gray 2013). A similar approach is adopted here by using a high-speed camera (Teledyne DALSA Falcon2 with macro lens) to capture images of the top of the flow at a frame rate of 2000 frames per second (fps). The resolution is sufficient to capture details of circa 100 µm, and individual grains are seen to travel approximately one grain diameter between frames for the highest slope angles. These images are processed with the PIV software PIVlab for MATLAB (Thielicke & Stamhuis 2014) to calculate the surface velocity profiles u s (y) (where y is the transverse coordinate across the chute) at different downstream locations X. The camera is always aligned perpendicular to the local slope, so the measured velocity corresponds to the curvilinear value, and the velocities are averaged over the time window that the images cover.
The results of this PIV are shown in figure 7 for slope angles θ = 38 • and θ = 40 • . All of the profiles show clear curvature across the chute, with wall friction effects leading to noticeably smaller velocities towards the lateral boundaries (Brodu et al. 2013;Faug et al. 2015;Baker et al. 2016). This confirms that any measurements made from the side may not be representative of the overall flow, and also suggests that the effective friction coefficient should account for friction at the wall as well as the base. Following Taberlet et al. (2003), an effective friction coefficient that models both the basal and side wall contributions can be expressed as FIGURE 7. Time-averaged surface velocity profiles u s (y) at different downslope locations X, calculated using high-speed camera images and PIV from the top of the flow.
where µ b and µ w are constant coefficients for the basal and wall friction, respectively, and W is the chute width. This is a friction law designed for smooth frictional beds, as opposed to a rough bed friction law, such as those of Pouliquen (1999) and Pouliquen & Forterre (2002), where a layer of particles is glued to the bed. The thickness-dependent wall friction term will be negligible for a thin accelerating flow, but might become significant after the shock where the flow thickness is comparable to the channel width. Note that expression (3.25) only has empirical support for h/W < 1.5 (Taberlet et al. 2003), but the flows studied in this paper remain in this regime. Substituting (3.25) into the ODE (3.17) allows the thickness profile, and subsequently the depth-averaged velocity, to be calculated for material leaving the gate. Figure 8 shows the resulting velocity profiles, which have been converted to Cartesian downslope distances X, for ease of comparison, by inverting the integral mapping (3.3). Material accelerates downslope as it leaves the hopper, is slowed down by the topography and then accelerates once more after passing the bump. All of the experimental measurements are taken upstream of the take-off point for the formation of the jet, which implies that the equations for dense granular flows still apply. For calibration, the PIV results are averaged in the cross-slope direction (since the model is one-dimensional). Assuming plug flow, the surface velocities correspond to the depth-averagedū, and choosing µ b = tan(23 • ) and µ w = tan(7.5 • ) gives a good fit for both slope angles, as shown in figure 8 (solid lines). Note that this agreement provides justification for the choice Fr 0 = 1 at the inflow, since trying to model the experimental results using larger values Fr 0 > 1 give unphysical values for the friction parameters and/or less satisfactory fits to the data points. Also shown in figure 8 are the theoretical profiles in the absence of wall friction (dashed lines). In this thin accelerating regime this additional friction term is almost negligible, but it is expected to play a more significant role in the thicker flows near the shock and hence is included in the effective friction coefficient µ. A not-to-scale schematic of the basal topography is included for both slope angles, with shaded areas representing where µ b > tan ζ .

Formation of a jet downstream of the obstacle
When there are no (or too few) grains initially upstream of the bump, a jet usually appears downstream of the obstacle (see figures 1, 3(a) and 4, as well as movie 1). When the flow is released from the hopper it accelerates downslope before being directed upwards (relative to the chute base) and detaching from the topography at a take-off point, where all the particles lose contact with the base and proceed to follow a parabolic trajectory. The take-off point corresponds to where the normal traction at the base of the avalanche vanishes. Following Gray et al. (1999), this normal traction in curvilinear coordinates is expressed by where σ b is the Cauchy stress at the base, ρ is the flow density (assumed constant) and h, ζ , κ andū are all x-dependent. Setting (4.1) to zero, the normal traction vanishes whenū 2 = − g cos ζ κ . (4.2) The depth-averaged velocity along the chute can be calculated by integrating the ODE (3.17) to determine h(x) and then using the mass balance equation (3.15) to determineū, as shown in § 3.2 and figure 8. The computed velocity can then be used to determine the take-off point x = x J , where condition (4.2) is satisfied. The corresponding flow thickness h J = h(x J ) and depth-averaged velocityū J =ū(x J ) can then be extracted and used as initial conditions for the jet, whose motion is computed by using Newton's second law (assuming the only body forces are due to gravity) to determine the trajectory of airborne particles. Note that it is much easier to compute the jet dynamics in the aligned Cartesian coordinate system than in the curvilinear system. In Cartesian variables the initial take-off point at the free surface of the flow (X JS , Z JS ) is found from the curvilinear coordinates by the projections

4)
Multiple solutions for granular flow over a bump 93 where X(x J ) represents the Cartesian abscissa of the take-off point x J , ζ J = ζ (x J ) is the inclination of the x-axis at the take-off point and β J = θ − ζ J is the angle between the z and Z axes at this same position. Throughout this paper the free surface of the non-airborne avalanche is reconstructed with a similar projection method to (4.3)-(4.4) at a general point (x, h). Assuming the dominant velocity is purely downslope at x = x J (i.e. the flow is tangential to the topography at the take-off point), the initial take-off velocity has components (U J , W J ) in the downslope and normal directions, respectively, where (4.6) The trajectories (X(t), Z(t)) are then governed by the second-order ODEs which may be integrated subject to the initial conditions (4.3)-(4.6) to give X(t) = 1 2 gt 2 sin θ + U J t + X JS , (4.9) (4.10) Solving (4.9) for time t and substituting into (4.10) implies that Z(X) = − 1 2 gT(X) 2 cos θ + W J T(X) + Z JS , (4.11) where (4.12) and the positive root is chosen in (4.12) to ensure that X > X JS for positive times. Equations (4.11), (4.12) determine the trajectory of the top surface of the jet once the initial positions (4.3) and (4.4) are determined. Trajectories for particles starting at heights between zero and h J may also be constructed in a similar manner, but lie close to the surface jet trajectory, due to the shallowness of the flow and the assumption of plug flow prior to take-off. Figure 9 shows a comparison between the theoretical and experimental surface jet trajectories, at two different slope angles θ = 38 • and θ = 40 • . In both cases, the (dashed) paths given by (4.11) and (4.12) show fairly good agreement with experiments for the centre of the jet. There is, however, significant vertical spreading of the airborne material, with the effective flow thickness becoming much greater compared to the upstream region. One possible explanation for this is that there is vertical shear through the depth of the avalanche prior to take-off, since this would imply that the upper particles have greater take-off velocities and follow higher trajectories. However, as already noted, the amount of shear is low, and a more probable reason for the vertical spreading is the effect of wall friction. The theory is quasi-one-dimensional, with a single mean value of velocity used at each cross-slope position in the flow. However, as shown in figure 7, particles at the walls 94 S. Viroulet and others The solid lines are calculated using the measured surface velocities at the centre of the channel for the top of the jet, and the dash-dotted lines are using the wall velocity for the lower trajectory (see figure 7). The red dots represent the surface take-off points (x J , h J ). Note that the images have been slightly rotated to maximise space. The bump height of 4.75 cm acts as a scale.
are slower than the mean and those in the centre are faster. The images in figure 9 are taken from the side of the chute through clear walls, so that the upper trajectories in figure 9 actually correspond to the fastest particles in the centre. Similarly, the lower paths are the slower moving wall particles. By keeping the take-off position x J predicted by (4.2) and instead imposing the surface velocity measured at the middle of the flow for X = 0.40 m (near the take-off point), much better agreement is found for the surface trajectories at both slope angles (see solid lines in figure 9). In the same way, using the surface velocity measured near the walls gives a significantly improved fit for the lower paths (see dash-dotted lines in figure 9). Hence, much of the apparent spreading of the jet can be accounted for by including wall friction effects. The underestimation of the top surface of the jet may also be explained by the fact that the surface velocities are slightly underestimated by the theory near the take-off point (especially for θ = 40 • ), as seen in figure 8.

Formation of a steady shock upstream of the obstacle
As observed during the experiments, introducing a static deposit of particles in front of the obstacle can cause a steady shock to form before the bump instead of a jet behind it ( figure 3 and movie 2). This shock separates a thin fast moving (supercritical) upstream region from a slower (subcritical) regime with an associated increase in thickness on the downstream side. Experimentally, there is a rapid, but smooth, transition across the shock, with a width that is controlled by the viscosity of the system (Whitham 1974). The mean shock position (referred to as simply shock position from this point onwards) is defined to be the point of maximal free-surface gradient dh/dx in this zone. Conversely, in the hyperbolic inviscid depth-averaged theory (3.7)-(3.8) the shock represents a mathematical discontinuity and therefore has an infinitesimally small width and unique downslope position. Despite this important difference, the model is able to predict the key features of the position of the shock and corresponding thickness profile.
The jump conditions govern the conservation of mass and momentum when the thickness and velocity are discontinuous. They are derived from the integral form of the conservation laws using a limiting argument (see e.g. Chadwick 1976). For the case of a one-dimensional stationary shock at x = x s , the general form (e.g. Gray et al. 2003; reduces to where the local chute inclination angle ζ s = ζ (x s ) and the jump brackets f = f + − f − denote the difference between the enclosed function on the forward (+) and rearward (−) sides of the shock. Note that the source terms do not enter into the momentum jump (5.2), since their integral vanishes as the control volume is shrunk down onto the shock (Chadwick 1976). Written out explicitly, equations (5.1) and (5.2) are giving two expressions relating the unknown quantities h ± ,ū ± . The shock position x s is also unknown at this stage. The mass jump condition (5.1) implies that hū is the same on either side of the shock and, using (3.15), it follows that the volume flux is equal to h 0ū0 everywhere on the chute, i.e. h 0ū0 = h −ū− = h +ū+ = hū. (5.5) This can be used to eliminate the velocitiesū ± in the momentum jump condition (5.4) and it also follows that the ODE (3.17) governs the thickness on both the upstream and downstream sides of the shock. To ensure that sufficient information is propagated into the shock in order to determine its position, the causality condition (see e.g. Ockendon et al. 2004) states that three families of characteristics must travel into the discontinuity. Note this is also equivalent to the Lax entropy condition (Lax 1957), which ensures that the vanishing viscosity solution of the equations is selected. Assuming thatū is positive, it follows from the definition of the characteristics (3.10)-(3.11) that this occurs if and only if the upstream side is supercritical (Fr > 1) and the downstream side is subcritical (Fr < 1), giving a range of permissible shock positions. The ODE (3.17) can then be integrated forward in space from x = 0 up to x = x s , giving h − andū − , before applying the jump conditions (5.3) and (5.4) to calculate h + andū + . The downstream region of the flow is found by integrating (3.17) with initial condition h = h + at x = x s .

S. Viroulet and others
O Z X FIGURE 10. Diagram showing a sketch of the solution structure of the shock. There is a supercritical flow (Fr > 1) out of the hopper, which transitions to a subcritical flow (Fr < 1) across a shock located at x s in curvilinear and X s in Cartesian coordinates, respectively. At x 1 and X 1 the flow transitions smoothly back from a subcritical to a supercritical flow of height h 1 . The inflow height is denoted h 0 .
Even after enforcing the causality condition, there is a non-unique way of choosing the shock position x s . This is in direct contrast to the experimental set-up, which selects a specific shock position for a given slope angle. There is a very simple mathematical resolution of this issue. Immediately after the shock the flow is subcritical, so the solution can either continue to decelerate and thicken, which is not what is observed, or it can accelerate. If it accelerates, then the Froude number will increase, and when it reaches unity the denominator of the ODE (3.17) will be zero and hence the thickness gradient becomes infinite (dh/dx → ∞), which is again unphysical. The only exception to this is when the numerator of (3.17) is also zero, which implies that the gradient is undefined, or at least must still be determined by L'Hôpital's rule. In this latter case it is possible to produce a smooth transition from subcritical to supercritical flow as the avalanche accelerates down the chute, as sketched in figure 10 and shown in figure 11. There is therefore a critical point, x = x 1 (say), where the Froude number is equal to one and hence the flow speed is equal to the gravity wave speed.
The critical point is analogous to the sonic point in gas dynamics (e.g. Laney 1998) and plays a vital role in selecting the correct shock position, in a parallel way to flows in collapsible tubes (Shapiro 1977). Denoting the flow thickness and velocity at x 1 by h 1 andū 1 , respectively, and recalling that Fr 0 = 1 at the inflow as well as at the critical point, it follows that where (5.6) represents the conservation of mass (3.15), and the slope angles at the inflow and the critical point are ζ 0 = ζ (0) and ζ 1 = ζ (x 1 ), respectively. Using (5.6) and (5.7) the thickness at the critical point is which can be substituted into the numerator of (3.17) to give h 3 1 g cos ζ 1 (tan ζ 1 − µ(h 1 )) − h 1 µ(h 1 )κ 1 (h 0ū0 ) 2 − 1 2 gh 4 1 κ 1 sin ζ 1 = 0. FIGURE 11. Comparison of the experimental and theoretical (white line) free-surface profile for the steady shock at a slope angle θ = 39 • . Note that the image has been slightly rotated to maximise space. The bump height of 4.75 cm acts as a scale.
unity and both the numerator and the denominator of the ODE (3.17) are equal to zero. The gradient of the solution at the critical point therefore has to be determined either by a Taylor series expansion or by L'Hôpital's rule (see appendix A for details). However, once the position and gradient of the critical point have been determined, the ODE (3.17) can be integrated both upstream and downstream of x = x 1 , to construct a solution that transitions smoothly from subcritical flow (for x < x 1 ) to supercritical flow (for x > x 1 ) as sketched in figure 10. The final part of the problem is to connect the smoothly varying solution through the critical point to the supercritical solution that emerges from the hopper described in § 3. By construction the mass jump condition (5.1) is satisfied everywhere, so it only remains to find the shock position x = x s where the momentum jump condition (5.2) is satisfied, which can again be solved for using a standard numerical root finding method. There are in fact two solutions; however, only the one that lies furthest upstream is stable to small perturbations and this is the one that is observed in experiments. Figure 11 shows a comparison of the thickness profile for the experiment and theory, at a slope angle θ = 39 • . The shock position is in relatively good agreement with the experiment. There are, however, some discrepancies in the flow thickness, especially in the shock region, where, as already noted, the flow is non-shallow and has a rapid (but smooth) transition in thickness. This is in contrast with the discontinuity predicted by the shallow-water model. It is possible to construct smooth shocks by including a depth-averaged version of the µ(I)-rheology (Gray & Edwards 2014) into the theory as in Edwards & Gray (2015). However, this is based on the rough-bed friction law (Pouliquen 1999), which is inconsistent with the Coulomb-type law (3.25) and thus not included here. Another important effect that is missing is the dilatation of the flow. Faug et al. (2015) recently observed that, for avalanches on steep slopes on a smooth base, the supercritical flow becomes more dilute as it leaves the gate and accelerates downstream, whereas after the shock the solids volume fraction is increased in the slower moving regime. This would potentially lead to an under prediction of the flow 98 S. Viroulet and others thickness in the accelerating region and an over prediction of the thickness after the shock.
In addition, a much greater degree of shear is observed in the thick region, immediately behind the shock, compared to the accelerating regions, with particles at the bottom becoming close to static, especially at lower angles of inclination. Long exposure photographs indicate that all the grains are in fact in motion, but there is certainly much greater shear through the flow, than in the region upstream of the shock. This shear is not accounted for by the theory, since plug flow is assumed everywhere. In the thick region the velocity profile at the side wall is better approximated by an exponential function (see e.g. Wiederseiner et al. 2011) which has a corresponding shape factor For the value of Λ = 3.24 found by Wiederseiner et al. (2011) this gives a shape factor χ = 1.75, which is no longer close to unity. It is therefore surprising that the theory is able to accurately describe such regions. There are two places in (3.8) where χ appears. The first is in the momentum transport term and the second is in the centrifugal correction to the basal pressure. It is useful to define two non-dimensional parameters P 1 = χhū 2 1 2 gh 2 cos ζ = 2χFr 2 , and P 2 = κχ hū 2 gh cos ζ = χ κhFr 2 . (5.12a,b) The first is the ratio of the momentum transport terms to the depth-averaged pressure gradient in the momentum balance (3.8), whilst the second is the ratio of the centrifugal pressure correction to the basal pressure in the source terms. The product κh is always less than unity, otherwise the avalanche is thicker than the intersection point of the curvilinear coordinates. Moreover the Froude number typically lies between 0.19 < Fr < 0.25 in the thick slowly moving region of the flow, as shown in figure 13. Hence the Froude number squared lies in the range 0.03 to 0.06, and the ratios P 1 and P 2 are very small unless χ is very large. It follows that even if there was a fairly large deviation of χ away from unity, its effect in the thick, slowly moving region would be small. The simple depth-averaged model therefore provides a good approximation in the slowly moving region despite the shape factor being unrepresentative of the shear profile. Conversely, in the thin fast moving regions, where momentum transport terms and centrifugal effects can dominate, there is considerable slip at the base of the avalanche, so that χ = 1 (as assumed) is a good approximation. Remarkably, then, despite the simplicity of the model there is an excellent agreement between the experimental and theoretical shock position for a wide range of slope angles, using the same friction coefficients obtained in § 3.2, as shown in figure 12.

Time and spatially dependent numerical simulations
The full governing equations (3.7)-(3.8) are now solved numerically to show the transient evolution to the equilibrium shock state described in the previous section, as well as the formation of a static deposit on the bump when the inflow ceases. The hyperbolic system is solved using the shock-capturing non-oscillatory central scheme of Kurganov & Tadmor (2000), whose semi-discrete formulation is combined with a second-order Runge-Kutta time integrator.
In order to solve the system, the depth-averaged equations (3.7), (3.8) together with the friction law (3.25) are written in terms of conserved variables in vector form as The critical inflow Fr = 1 at x = 0 means that two inflow conditions must be specified (e.g. Weiyan 1992), which are the same as those for the steady ODE, i.e. where the function h init (x) is prescribed. Figure 14 shows a simulation of an avalanche as it flows down a plane inclined at θ = 39 • and hits a static deposit that has been left by a previously computed avalanche. In fact, the final deposit that is computed at the end of this simulation, which is shown in figure 15, is indistinguishable from the initially assumed deposit, so the solution is quasi-periodic in nature.
As the avalanche first impinges on the deposit, at approximately 0.22 s, it generates an erosion shock that accelerates the thin static deposit on the bump into a slowly moving thick layer of material. This shock, which is similar to those observed by Edwards & Gray (2015), propagates rapidly downslope and reaches the end of the static deposit at approximately 0.52 s, and dissipates. At the same time as the static material is being mobilised, the fast moving upstream avalanche produces a shock at the rear of the thick slowly moving layer, which also moves downslope, but at a much slower rate than the erosive shock at the front. At approximately 0.68 s this shock reverses direction and slowly starts to propagate upstream, reaching the steady-state equilibrium position at approximately t = 5 s, with no further propagation of the shock occurring before the final image in figure 14 at t = 8 s.
It is very important to note that if the static layer of grains were replaced by a rigid smooth topography with the same shape, the flow would automatically develop into the jet solution, which would be even stronger due to the slightly steeper slope. In the absence of the dissipation mechanism, provided by the erodible grains, the 102 S. Viroulet and others FIGURE 15. Temporal evolution of the free-surface height once the inflow ceases at slope angle θ = 39 • . Note that the images have been slightly rotated to maximise space, the aspect ratio is 1:1 and the bump height of 4.75 cm acts as a scale. The final image corresponds to t steady = 15 s. A supplementary movie is available online. oncoming avalanche tip would simply propagate smoothly across the new topography in the same way as it does upstream and downstream of the erodible layer. The erodible material is therefore vital for the brief formation of the downslope-facing erosive shock shown at t = 0.3 and t = 0.4 s in figure 14. This dissipates a lot of energy by mobilising the static grains, and, crucially, thickens and slows the oncoming flow allowing the upslope-facing shock to be triggered. This very fine balancing of different physical effects is what makes this experiment such a sensitive and interesting one. Figure 15 shows another sequence of simulation images, which start from the steady-shock solution shown in figure 14. At t = 10 s the inflow is instantaneously closed off, with the fast moving tail of the avalanche reaching the shock at t = 10.24 s. At this point the steep shock collapses, with some material being pushed upslope until 10.50 s and the avalanche then slowly thinning as the residual material, that is able to creep over the bump, slowly flows downslope. By t = 14 s a static deposit is formed on the upstream side of the bump. This material is sufficient to trigger any subsequent oncoming avalanche into the steady-shock regime. As previously mentioned it is essentially this deposit that provided the initial condition (6.5a,b) for the beginning of the simulations shown in figure 14. A supplementary movie showing the full time-dependent development of the solution from the impingement of the avalanche onto the static deposit to the formation of the static deposit at the end is available online.
The numerical computations are compared to experiments in figure 16, where the free-surface profiles are overlaid on the experimental photos. The slope angle is again θ = 39 • , but this time the initial condition is set by the actual free-surface profile of the static grains in the experiment. As the avalanche impinges on the static deposit some of the particles are thrown into the air during the mobilisation process, rather than generating the simulated erosive shock on the downstream side of the slowly moving material. However, this anomaly rapidly propagates through the deposit and dissipates. The rearward shock between the slowly moving, recently mobilised, grains and the rapid upstream avalanche has already started to move upslope by t = 0.8 s in both the experiments and the computation. The experimental shock wave is diffuse and probably somewhat upstream of the computed shock, but it is in reasonable agreement. As the experimental shock attains its equilibrium position, described in § 5, the numerical shock catches up with it, so that at steady state, the mean position of the diffuse shock and the jump are remarkably close. This is a particularly strong vindication of the theory given how sensitively dependent the shock position is on the basal friction and the chute inclination. It should be noted, however, that when the grains flow over the bump, the avalanche is on the verge of detaching from the base at t = 0.8 s, and the flow is quite dilute on the downstream side. As the backward-propagating shock establishes itself, it becomes more efficient at slowing the grains, and the dilute region on the lee side of the bump is reduced in size, although there is some dilation even at steady state, which occurs at times greater than 4 s.

Comparison between standard and terrain-following approaches
It is also possible to derive a simpler depth-averaged theory in aligned Cartesian coordinates OXZ for modelling the shock evolution and steady-state position (see e.g. Savage & Hutter 1989;Gray et al. 2003;Gray & Edwards 2014). This is referred to as the standard depth-averaged avalanche model throughout this paper. This theory 104 S. Viroulet and others FIGURE 16. Temporal evolution of the theoretical free-surface profile (red solid lines) to the exact steady-state solution (green dashed line) for a shock at slope angle θ = 39 • . Note that the images have been slightly rotated to maximise space, the aspect ratio is 1:1 and the bump height of 4.75 cm acts as a scale. assumes that the slope is inclined at a fixed mean angle and the effects of terrain only enter the equations via a topography height gradient in the source term, i.e. there is no local variation in the slope angle and there are no centrifugal effects. For this rather sensitive experiment, the terrain-following model (see e.g. Savage & Hutter 1991;Gray et al. 1999) described in § § 3-6 is much more successful at predicting the steady-state shock position than the standard model. In addition, the standard approach is also unable to predict the take-off point of the jet. These results therefore provide compelling evidence of the advantages of the terrain-following approach.
7.1. The standard depth-averaged avalanche model Using a Cartesian coordinate system OXZ aligned at a constant angle θ to the horizontal (as shown in figure 6), the flow thickness H(X, T) and depth-averaged velocityŪ(X, T) are now measured perpendicular and parallel, respectively, to the X-axis. The depth-integrated mass and momentum balance equations (see e.g. Savage & Hutter 1989;Gray et al. 2003;Gray & Edwards 2014) are then where the shape factor χ has been assumed to equal unity. The source term S is composed of the downslope component of gravity, Coulomb basal friction, wall friction and topography gradients, and is equal to Note that in the standard approach, the effect of topography arises purely from the downslope component of the basal pressure acting on the bump in the basal shear stress (see e.g. Gray et al. 2003;Gray & Edwards 2014). This is, therefore, a much less sophisticated treatment of the topography than the terrain-following approach, where both the local inclination angle varies as a function of x and the basal pressure accounts for centrifugal forces. It has, nevertheless, proved to be very useful for simulating complicated flows past obstacles, where shock waves are generated (see e.g. Gray et al. 2003;Aker & Bokhove 2008;Cui & Gray 2013). Following the same approach as in § 3.1, it follows that integrating the steady-state depth-integrated mass balance equation (7.1) in X implies that where the thickness and depth-averaged velocity at the inflow are H 0 andŪ 0 , respectively. Moreover, the mass jump condition, HŪ = 0, implies that the volume flux per unit width is equal to H 0Ū0 everywhere. The steady-state depth-integrated momentum balance (7.2) then implies that the ODE governing the thickness upstream and downstream of the shock is

S. Viroulet and others
where it has implicitly been assumed thatŪ > 0. Note that the denominator in (7.5) takes the same form as in the terrain-following system (3.17), and the numerator is similar except with a topography gradient replacing the curvature term. In the standard model, the Froude number is defined as (7.6) and this is again assumed to be equal to unity at the inflow. Hence, the inflow velocitȳ U 0 = √ gH 0 cos θ. At the critical point, which lies at X = X 1 , the Froude number is also equal to unity and the denominator in the ODE (7.5) is equal to zero. It follows that the thickness at the critical point, H 1 , is the same as the thickness at the inflow, i.e.
In order to ensure the thickness gradients in the ODE (7.5) remain bounded at X = X 1 the numerator must also be zero, which implies that Since the topography b = b(X) is given by (3.1) and H 1 = H 0 , equation (7.8) can be solved numerically to find the position of the critical point X = X 1 .
In the same way as in the terrain-following solution, the ODE (7.5) is integrated backwards, and forwards, away from the critical point to construct a solution for H that varies smoothly between subcritical (X < X 1 ) and supercritical flow (X > X 1 ). A smoothly varying supercritical solution can also be constructed from the inflow and the two portions are pieced together by applying the momentum jump condition HŪ 2 + 1 2 gH 2 cos θ = 0, (7.9) to find the position of the shock X = X s . Figure 17 shows the computed shock position as a function of the slope angle θ. The inflow conditions and wall friction µ w = tan(7.5 • ) are kept at the same constant values as in § 3, but different basal friction coefficients µ b are plotted (black dashed lines). For each angle, it is possible to choose µ b so that the shock position is predicted by the Cartesian theory, for example the value µ b = tan(27.5 • ) gives good agreement at θ = 39 • (dash-dotted line in figure 17). However, the same parameters are unable to reproduce the results at other slope angles, whereas the terrain-following approach matches the shock position across the whole range for the same friction coefficients (solid line). This is highlighted in figure 18, which shows a comparison between the thickness profiles obtained with both theories and the experiments. At 39 • there is little difference between the two theories and both fit the experimental data equally well. However, at 37 • the terrain-following system gives a significantly improved fit, with the standard model predicting a shock position far upstream of the actual location. Note that the basal friction value µ b = tan(27.5 • ) required for fitting the standard theory at θ = 39 • is unphysically high and not representative of that calculated using PIV (µ b = tan(23 • )). When using the experimentally measured value (µ b = tan(23 • )), the prediction of the shock position is far away from the actual values at low slope angles. At higher angles (θ = 36 • ) for µ b = tan(23 • ), the standard theory is unable to predict any shock solutions at all. The measured friction values are not high enough to be able to satisfy the momentum shock condition in the standard system. This is another Curvilinear terrain-following Cartesian terrain-following FIGURE 17. Shock position X s as a function of the slope angle θ. The symbols are experimental data, and the solid red line shows the terrain-following theory, computed in curvilinear coordinates with the measured µ b = tan(23 • ) and µ w = tan(7.5 • ) from PIV. The different black dashed lines denote the standard Cartesian theory, calculated with fixed wall friction µ w = tan(7.5 • ) and varying basal friction µ b , with the value µ b = tan(27.5 • ) chosen as the best fit to the shock position at θ = 39 • (dash-dotted line). The green dashed line shows the terrain-following avalanche model computed independently in Cartesian coordinates, which, as expected, exactly reproduce the curvilinear results.
advantage of the terrain-following model since shocks are observed experimentally for a wider range of angles than the standard theory predicts.
7.2. Coordinate transformation of the terrain-following model For practical applications the topography is often defined by a geographical information system (GIS) model in Cartesian coordinates with a vertical axis aligned with gravity. In order to benefit from the improved predictive power of the terrain-following model it would therefore be useful to express it in a Cartesian coordinate system. As a specific example of this, it is now shown how the terrain-following model can be transformed from the curvilinear Oxz coordinates to the aligned Cartesian coordinates OXZ used by the standard model. This can be achieved by making the simple transformation of variables (7.10a,b) which implies that the derivatives where the normalisation factor ∆ b , defined in (3.6), is repeated here for completeness (7.12) Since ∆ b is a function of X only, it can be incorporated inside the time derivative to write the terrain-following mass and momentum balances (3.7)-(3.8) in conservative form and in Cartesian coordinates as 2 gh 2 cos ζ = ∆ b hS terrain-following , (7.14) where the terrain-following source term is S terrain-following = g sin ζ −ū |ū| µ(g cos ζ + χ κū 2 ) . (7.15) It is important to recognise that the transformation (7.10a,b) does not change the model, i.e. it is still the terrain-following theory of Savage & Hutter (1991) and Gray et al. (1999). It is just expressed in a Cartesian coordinate system. For GIS applications, the same transformation will also work, except that now the topography height (3.1) must be given relative to Cartesian coordinates aligned with gravity, rather than the inclined coordinates used by the standard model. This idea is not new. Bouchut et al. (2003) and Bouchut & Westdickenberg (2004) have investigated various general transformations of the Savage & Hutter (1991) model, including into an arbitrary system of final coordinates. A simplified treatment is provided in § 3.1 and appendix A of Mangeney et al. (2007), but the new features of the transformed model are not really exploited fully, since the paper is focused on levee formation on an inclined plane, where the standard and terrain-following models are identical.
7.3. The steady-shock solution for the terrain-following model in Cartesian coordinates The steady-shock problem of § 5 is now considered again to explicitly demonstrate that the transformation of variables (7.10a,b) does not change the resultant predictions of the terrain-following model. In general, the unsteady jump conditions for the depthaveraged mass and momentum balances, (7.13) and (7.14), are hū − ∆ b hV n = 0, (7.16) hū 2 − ∆ b hūV n + 1 2 gh 2 cos ζ = 0, (7.17) where V n = dX/dT is the shock speed expressed in Cartesian variables. Since the curvilinear shock velocity v n = ∆ b V n , it follows that the Cartesian curvilinear jump conditions (7.16) and (7.17) are precisely the same as the general unsteady curvilinear jump conditions. In particular, when the shock speed is equal to zero, these also reduce to the steady curvilinear jump conditions (5.1) and (5.2) used in the shock solution in § 5. Assuming steady fully developed flow, the mass balance (7.13) can be integrated subject to the condition that h = h 0 andū =ū 0 at the inflow to show that hū = h 0ū0 , (7.18) which, by virtue of the steady mass jump condition, is valid everywhere. Similarly, assuming that the shape factor χ is equal to unity, the steady-state depth-averaged momentum balance (7.14) can be expanded using (7.18) and (3.5) to show that the ODE governing the flow thickness is

S. Viroulet and others
To confirm this, the steady-state shock solution has been computed independently, using the terrain-following model in Cartesian coordinates. In figure 17 the shock position as a function of inclination angle, computed using Cartesian coordinates, is shown by the green dashed line. It lies precisely on top of the red solid line, which shows the solution of the terrain-following model computed in curvilinear coordinates. The results of the terrain-following model are, as one would expect, therefore completely independent of the coordinate system that the equations are solved in.
Using the terrain-following model in Cartesian coordinates does, however, have the advantage that if the topography is defined in Cartesian coordinates, it does not have to be converted into curvilinear variables before the model can be used. Moreover, the results no longer have to be mapped back from curvilinear coordinates in order to plot them. It should be noted, however, that both the thickness and the velocity are the curvilinear ones, which are implicitly assumed to lie normal and tangential to the local topography. As a result, in order to reconstruct the free-surface position it is still important to project the thickness using the equations X surface = X − h sin(θ − ζ (X)), (7.21) Z surface = b(X) + h cos(θ − ζ (X)).
( 7.22) These are in fact the same formulae as (4.3) and (4.4), which were used in § 4 to project the free surface at the take-off point back into Cartesian coordinates. It is therefore possible to have the improved accuracy of the terrain-following model, (7.13)-(7.15), while still using a simple Cartesian coordinate system.

Conclusion
In this paper small-scale experiments have been used to show that the flow of a dry granular material over variable topography may exhibit two very different types of behaviour depending on the initial conditions. On an initially empty chute, the avalanche accelerates as it leaves the hopper and may reach fast enough velocities to take-off as it flows over a smooth bump, forming a detached airborne 'jet' (figure 1 and supplementary movie 1). In contrast, placing a sufficient mass of static, erodible particles slightly upstream of the obstacle reduces the energy of the impacting flow, generating a sharp decrease in velocity and associated increase in flow thickness, or 'shock' that propagates upstream until it reaches a steady-state position (figure 2 and movie 2). Once a shock is created, its equilibrium position does not depend on the amount of material placed in front of the obstacle and is stable to perturbations in the base flow (see supplementary movie 3). A depth-averaged terrain-following avalanche model (Savage & Hutter 1991;Gray et al. 1999) that is formulated in slope fitted curvilinear coordinates is used to model these two types of behaviour. Based on experimental measurements of the upstream accelerating region, which is present in both the jet and shock regimes, a friction coefficient has been introduced in the equations that takes both basal and wall friction into account. The theory is then used to calculate the take-off position, corresponding to the point in the flow where the normal traction at the base of the avalanche vanishes. The trajectory of the jet is then calculated by assuming that the particles follow a series of ballistic trajectories. The mean particle path is found to be in good agreement with experiments, and the vertical spreading of the jet is largely accounted for by including wall friction effects (figure 9), although details of the shear profile through the depth of the avalanche may also be required for more accurate modelling.
For the steady-state shock regime, material accelerates and thins as it leaves the hopper, and the shock represents a transition from supercritical (Fr > 1) to subcritical (Fr < 1) flow. Downslope of the bump, the flow accelerates to become supercritical once more, and consequently passes through a critical point, where Fr = 1. At this point the denominator of the ODE (3.17), which governs the free-surface height, is equal to zero. The critical point therefore plays a vital role in determining the position of the steady-state shock, since smooth solutions that pass through it are only possible if the numerator of the ODE (3.17) is also zero. The mass and momentum jump conditions, (5.1) and (5.2), can then be used to find a unique stable position for the shock that connects the smooth solution out of the hopper with the smooth solution through the critical point, as sketched in figure 10 and shown in figure 11. Despite the simplicity of the theory, figure 12 shows that it is able to match the mean position of the experimentally measured shock location for a wide range of inclination angles θ using the previously determined friction coefficients. The largest discrepancy between the solution and the experiment occurs close to the shock, which is smoothed out in experiments rather than being a discontinuity. This is due to high shear in this region, so that the nonlinear granular viscous effects (GDR-MiDi 2004;Jop, Forterre & Pouliquen 2006;Gray & Edwards 2014;Baker et al. 2016) become important. Smaller discrepancies may also be due to the flow dilating as it accelerates down the slope and then rapidly compressing as it flows over the bump .
Full time-dependent numerical simulations using a high resolution shock-capturing method are shown in figures 14 and 15 as well as in a supplementary online movie. These snapshots and animations highlight why this is such a sensitive and interesting problem. In particular, the formation of an erosive shock (second and third snapshots in figure 14) as the avalanche propagates over the erodible layer is of great importance. It mobilises the static grains into a thick slowly moving layer, while dissipating significant amounts of energy from the oncoming flow. As a result, the thin high-speed flow out of the hopper collides with the rear of the slowly moving thick layer and forms a second shock, that initially moves downstream, before reversing and propagating upslope to its steady-state location. The existence of the erosive shock is therefore critical in allowing the avalanche to self-trigger into the steady-shock state.
The depth-averaged terrain-following theory formulated in curvilinear coordinates is compared to a standard avalanche model (Savage & Hutter 1989;Gray et al. 2003;Gray & Edwards 2014). This uses an aligned Cartesian coordinate system and the topography only enters the equations through the source term, through gradients in topography height. This is a much less sophisticated treatment of the topography than the terrain-following theory, which includes local changes in the slope inclination and centrifugal effects. The standard theory is unable to predict the take-off point of the jet, which falls out naturally in the terrain-following model. In addition, the steady-state shock position is much more accurately captured using topography-following coordinates. Whilst it is possible to choose friction parameters so that the standard theory matches the experiments for a given inclination angle, accurate fitting cannot be obtained across the whole slope angle range (figure 17 ) and the parameter values required to fit at a specified slope angle do not correspond to those calculated using PIV. This is in contrast to the terrain-following system, where the same values measured from the accelerating regime can be used to predict the shock position for all slope angles, as well as the jet trajectories. The local variation in the inclination of the terrain-following coordinates together with the centrifugal effects in the source terms therefore play an important role in achieving a more accurate flow representation.
Both of the regimes investigated in this paper could have important implications for the understanding and simulation of natural flows. The formation of an airborne jet, characterised by vanishing normal basal forces as material flows over complex topography, could be an important design consideration for flow defences and could, in principal, be inferred from seismic measurements (e.g. Favreau et al. 2010;Levy et al. 2015). The fact that the shock is generated as the flow moves over an erodible layer of particles is also significant because, in a geophysical context, such a deposit could be left behind from a previous event. Consequently, the characteristics of an initial flow and any subsequent flows may be drastically different, even if the source conditions are similar. This is consistent with the work of Moretti et al. (2012), who showed that the presence of erodible material on the slope strongly affected the resulting dynamics.
When modelling geophysical flows it is advantageous to be able to use a Cartesian coordinate system, so that GIS data can be used to define the topography height above a fixed datum. In § 7.2 a coordinate transformation is therefore introduced to convert the terrain-following model of Gray et al. (1999) from curvilinear to Cartesian coordinates. This is particularly useful in this paper, since the topography height b(X), the local slope angle ζ (X) and the curvature κ(X), are all defined in Cartesian coordinates. If the Cartesian formulation of the terrain-following model is used, then all of these variables no longer need to be converted into functions of the curvilinear coordinate x, nor do the subsequent results need to be transformed back to Cartesians using the transformation (3.3). It must be stressed, however, that the coordinate system used by the terrain-following model of Gray et al. (1999) does not change the results, i.e. it does not matter whether curvilinear or Cartesian coordinates are used. All the results for the jet, in § 4, as well as the steady and unsteady shock, in § § 5 and 6, will be precisely the same when computed in Cartesian coordinates. To demonstrate this, the steady shock has been computed independently using the terrain-following model in Cartesian coordinates and the results are shown by the dashed green curve in figure 17. As expected, the shock position as a function of the chute inclination angle, is precisely the same as in curvilinear formulation (solid red line). It is important to note that when using the terrain-following model the thickness h and the velocityū are still implicitly assumed to lie normal and tangential to the local topography, and so the simple transformations (7.21) and (7.22) need to be applied when reconstructing the free-surface profile.
In order to solve for the steady-state free-surface profile using the terrain-following model it is useful to define the gradient of the solution at the critical point. This can Downloaded from https://www.cambridge.org/core. Multiple solutions for granular flow over a bump 113 be done by considering the ODE (3.17), which is restated here, as dh dx = h 3 g cos ζ (tan ζ − µ) − hµκ(h 0ū0 ) 2 − 1 2 gh 4 κ sin ζ h 3 g cos ζ − (h 0ū0 ) 2 .
(A 1) The critical point is the point x = x 1 , where the thickness h = h 1 and the Froude number is equal to unity. It has the special property that both the numerator N = h 3 g cos ζ (tan ζ − µ) − hµκ(h 0ū0 ) 2 − 1 2 gh 4 κ sin ζ , (A 2) and the denominator D = h 3 g cos ζ − (h 0ū0 ) 2 , (A 3) of equation (A 1) are zero, and hence that the gradient is undefined. It can, however, be calculated using L'Hôpital's rule. Since the topography (3.1) is defined in Cartesian coordinates, it is easiest to take the limit in X-coordinates, i.e.

S. Viroulet and others
with h = h 1 and substituting them into (A 4) yields a quadratic equation for dh/dX at X = X 1 , which can be solved for the gradient at the critical point. This is the gradient required when solving the terrain-following model in Cartesian coordinates, as in § 7.2. To convert the gradient to curvilinear coordinates used in § 5 one needs to divide by the normalisation factor ∆ b , i.e. dh/dx = (1/∆ b ) dh/dX.