Subcritical and supercritical granular flow around an obstacle on a rough inclined plane

Abstract

A blunt obstacle in the path of a rapid granular avalanche generates a bow shock (a jump in the avalanche thickness and velocity), a region of static grains upstream of the obstacle, and a grain-free region downstream. Here, it is shown that this interaction is qualitatively altered if the incline on which the avalanche is flowing is changed from smooth to rough. On a rough incline, the friction between the grains and the incline depends on the flow thickness and speed, which allows both rapid (supercritical) and slow (subcritical) steady uniform avalanches. For supercritical experimental flows, the material is diverted around a blunt obstacle by the formation of a bow shock and a static dead zone upstream of the obstacle. Downslope, a grain-free vacuum region forms, but, in contrast to flows on smooth beds, static levees form at the boundary between the vacuum region and the flow. In slower, subcritical, flows the flow is diverted smoothly around the dead zone and the obstacle without forming a bow shock. After the avalanche stops, signatures of the dead zone, levees and (in subcritical flows) a deeper region upslope of the obstacle are frozen into the deposit. To capture this behaviour, numerical simulations are performed with a depth-averaged avalanche model that includes frictional hysteresis and depth-averaged viscous terms, which are needed to accurately model the flowing and deposited regions. These results may be directly relevant to geophysical mass flows and snow avalanches, which flow over rough terrain and may impact barriers or other infrastructure.

1. Introduction

When a granular avalanche impacts a blunt obstacle, material is diverted around the obstruction resulting in changes to the speed and direction of the flow. For rapid flows on smooth beds, this includes the formation of a bow shock upstream of a static region of grains that is deposited at the front face of the obstacle (Gray, Tai & Noelle 2003). Downslope of the obstacle, the diverted material expands into the space behind the obstacle forming a grain-free region. The speed and thickness of an avalanche play a key role in determining its dynamics and deposits. Shocks and bores can occur whenever the speed of the avalanche exceeds the local propagation speed of the surface gravity waves for that material (Savage 1979; Brennen, Sieck & Paslaski 1983; Savage 1984; Gray et al. 2003; Gray & Cui 2007; Pudasaini et al. 2007; Cui & Gray 2013; Faug et al. 2015; Mejean, Faug & Einav 2017). Such flows where the speed of the avalanche exceeds the wave speed are called supercritical flows, and correspondingly where the avalanche speed is less than the wave speed the flow is subcritical. The Froude number is defined as the ratio of the depth-averaged flow speed to the speed of propagation of gravity waves

(1.1)\begin{equation} {\textit{Fr}} = \frac{|{\boldsymbol{\bar{u}}}|} {\sqrt{g h \cos \zeta}}, \end{equation}

where ${\boldsymbol{\bar{u}}}$ is the depth-averaged flow velocity, $h$ the flow thickness on a slope inclined to an angle $\zeta$ and $g$ is the acceleration due to gravity. The Froude number, therefore, determines whether the flow is supercritical (${ {Fr}}>1$), critical (${ {Fr}}=1$) or subcritical (${ {Fr}}<1$). This is analogous to the Mach number in gas dynamics, which is defined as the ratio of the speed of a compressible flow of air to the relevant wave propagation speed, i.e. the speed of sound (Johnson 2020).

Gray et al. (2003) conducted experiments of supercritical granular flows impacting a pyramidal obstacle fixed to a smooth bed. They found that the shape of the obstacle plays a key role and analysed two cases where the flow impacted directly onto a vertex, or onto the face of a triangular pyramid. For a pointed obstacle, oblique shocks form, whereas for a blunt obstacle, a bow shock forms upstream of a region of static deposited grains (termed a dead zone) that also contributes to the flow diversion. These ideas have been expanded upon with studies looking into the properties of the shock and how it changes depending on the angle of the impacted face (Hákonardóttir & Hogg 2005; Gray & Cui 2007) and the obstacle shape (Vreman et al. 2007; Cui & Gray 2013). The diversion of material results in an increase to the flow speed and thickness around the obstacle and on the lee side a grain-free region is formed as the material expands. This protected grain-free region is referred to as a granular vacuum, by analogy to the similar feature in gas dynamics (Gray et al. 2003; Hogg, Gray & Cui 2005; Gray & Cui 2007; Cui & Gray 2013). Through these studies we have gained a good understanding of supercritical granular flows on smooth beds.

Initial studies for obstacle interaction on a rough bed were conducted by Faug, Naaim & Naaim-Bouvet (2004) who looked into the effect of slope angle and fence height on the retention of material from a flow impacting a fence on a roughened slope. Rough beds allow steady uniform flows to develop for a wide range of flow depths and slope inclination angles. This is not the case for smooth beds without sidewalls, where the gravitational and frictional source terms in the depth-averaged downslope momentum balance imply a net accelerative or decelerative background motion. In particular, rough beds allow subcritical, as well as supercritical, granular flow around a blunt obstacle to be studied. For supercritical flow on a rough bed, a bow shock and potentially a static dead zone form upstream of the obstacle, as shown in figure 1(a). This is similar to what one would expect on a smooth bed, although the shape and size of the dead zone differ. On the lee side of the obstacle, instead of the fastest flow being located at the edge of the vacuum region (as is the case on a smooth bed), levees of static material form around the vacuum region. This is a new feature for flows impacting obstacles and is the opposite behaviour to that on the smooth bed.

Figure 1. Panel (a) shows an upstream supercritical flow (${ {Fr}} \approx 2.35$, $\zeta = 35^\circ$) of $d=300 - 400$ ${\rm \mu}$m soft masonry sand impacting a rectangular obstacle of dimension $3 \times 8 \times 3$ cm on a rough bed made from 750–1000 ${\rm \mu} $m turquoise glass ballotini. The shutter speed is $1/8$ s and the photographs are contrast enhanced. This procedure highlights the difference between static and moving regions of grains through motion blur. For supercritical upstream flow, a bow shock can be observed upstream of a static dead zone, while on the lee side of the obstacle a grain-free vacuum region opens up that is bounded by static levees. On the downstream face of the obstacle the magnetic fixings are visible, but they do not influence the flow as they are entirely contained in the vacuum region. Panel (b) shows an upstream subcritical flow (${ {Fr}} \approx 0.77$, $\zeta = 32^\circ$). In this case a bow shock does not form, but there is an upstream region that decelerates the flow and smoothly deflects the oncoming grains around the static dead zone and the obstacle. On the lee side the levees are wider and thicker than in the supercritical case. The sand, basal surface and obstacle dimensions are the same in this and subsequent experiments. Movies showing the time-dependent evolution are available in the online supplementary material.

For a subcritical inflow on a rough bed (figure 1b) there is no bow shock, but there is an upstream region that decelerates the flow and allows it to smoothly deflect around the static dead zone and the obstacle. The static deposits of material, both in the levees bounding the vacuum region and in the upstream dead zone, are larger than their supercritical counterparts. In both supercritical and subcritical flows, the levees provide an important buffer zone between the flow and vacuum region. The vacuum boundary is therefore very stable to perturbations in the main flow, and is able to extend far downstream, as shown in figure 2(a). Unlike a smooth bed, where all the grains (except those in the dead zone) flow off the inclined plane when the inflow ceases (Gray et al. 2003), on a rough bed a static deposit of grains forms in the region that was occupied by flowing grains, as shown in figure 2(a,b). This is related to the rough-bed friction law (Pouliquen 1999; Pouliquen & Forterre 2002; Edwards et al. 2017, 2019), which is much more complicated than the simple Coulomb friction used in many depth-averaged numerical simulations on smooth beds (Savage & Hutter 1989; Gray et al. 2003; Cui & Gray 2013). As well as a thin coating of grains, the subcritical flow leaves behind a flat-topped pile adjacent to the obstacle and a central ridge of material that extends far upstream (figure 2a).

Figure 2. Panels (a,b) show photographs of the final deposits of an upstream subcritical flow (${ {Fr}} \approx 0.77$, $\zeta = 32^\circ$) once the flow has subsided. The flat-topped pile above the obstacle is shown in panel (a) and the raised levees in the obstacle wake in panel (b). The orientation of the photographs is such that in panel (a) the upstream face of the obstacle is shown, whereas in panel (b) the view is changed to show the rear downstream face.

There is no universal constitutive law for granular flows that captures the complete behaviour observed in experiments and geophysical flows. Due to this there are a variety of approaches to modelling avalanches. Flow-obstacle interaction has been studied using discrete particle method simulations and compared with experiments undertaken on smooth beds (see, e.g. Teufelsbauer et al. 2009, 2011). This approach has also been used to examine the impact pressure on the obstacle itself for both rigid and flexible bodies (Teufelsbauer et al. 2011; Zhan et al. 2019). Discrete methods have the disadvantage that scaling up numerical experiments becomes very computationally expensive, and, as such, to consider the global effect of an avalanche impacting an obstacle, continuum methods are required. Recent advances in our understanding of the constitutive behaviour of dry granular flows with the so called $\mu (I)$-rheology (where $\mu$ is the friction and $I$ is the inertial number) (GDR-MiDi 2004; Jop, Forterre & Pouliquen 2006; Barker et al. 2015; Barker & Gray 2017; Schaeffer et al. 2019; Rauter, Barker & Fellin 2020) make continuum simulations of avalanches possible (Lagrée, Staron & Popinet 2011; Baker, Barker & Gray 2016a; Barker et al. 2021). However, these laws do not contain the non-local effects (Kamrin & Koval 2012; Kamrin & Henann 2015) that are necessary to model the more complex flow behaviour on rough beds, such as the formation of static levees during the flow (figure 1) and the coating of the inclined plane with a layer of grains at the end of the experiment (figure 2).

Avalanches typically have a shallow aspect ratio (i.e. the flow depth is much less than its length and width), which allows them to be successfully modelled using a depth-averaged approach (Grigorian, Eglit & Iakimov 1967; Eglit 1983; Savage & Hutter 1989; Gray, Wieland & Hutter 1999; Gray & Edwards 2014). Typically the system of equations resembles those of the shallow water or St Venant equations of fluid mechanics, with additional depth-averaged source terms to represent the effect of gravity driving the grains down the slope, basal topography gradients and frictional resistance to motion. These equations remove one spatial dimension from the problem and have the advantage that they scale up easily from laboratory experiments to full-scale geophysical flows (Kokelaar et al. 2017). The constitutive behaviour is usually encapsulated in the frictional source term. Savage & Hutter (1989) assumed a constant Coulomb friction $\mu$ on a slope inclined at an angle $\zeta$ to the horizontal. As a result, the flow is accelerative if $\mu <\tan \zeta$, decelerative if $\mu >\tan \zeta$ and non-accelerative when $\mu =\tan \zeta$. This is a reasonably good approximation for short smooth beds (Cui & Gray 2013; Viroulet et al. 2017), but over longer distances the avalanche does not tend to a terminal velocity, which is unrealistic. As a result, velocity-dependant basal-drag laws have been introduced in snow-avalanche models, such as the RAMMS model, which is widely used for hazard zoning in Switzerland (Salm 1993; Bouchet et al. 2004; Naaim et al. 2004; Christen, Kowalski & Bartelt 2010; Bartelt et al. 2015).

More complicated frictional models are required for dry granular flows on rough beds, where the static material is metastable for thicknesses in the range $h\in [h_{stop}(\zeta ), h_{start}(\zeta )]$ (Daerr & Douady 1999). This hysteretic behaviour implies that at the same inclination angle $\zeta$ and thickness $h$, there can be coexisting regions of static and flowing material. In particular, when a steady uniform flow is brought to a rest (by ceasing the inflow) it leaves a deposit of thickness $h = h_{stop}(\zeta )$ on a slope of angle $\zeta$. The slope angle can then be increased until $\zeta _{start}(h)$ before the grains are remobilised, where $\zeta _{start}(h)$ is the inverse function of $h_{start}(\zeta )$. Pouliquen & Forterre (2002) showed that this behaviour could be captured by making the friction $\mu$ a non-monotonic function of the Froude number and flow thickness. This law, which was originally developed for spherical glass ballotini, has been extended to model the hysteresis of angular grains, such as carborundum and sand (Edwards et al. 2017, 2019).

A rich variety of flow features, such as roll waves (Forterre & Pouliquen 2003; Gray & Edwards 2014; Razis et al. 2014; Viroulet et al. 2018), erosion-deposition waves (Edwards & Gray 2015), retrogressive failures (Daerr & Douady 1999; Russell et al. 2019), as well as levees, troughs and elevated channels (Edwards et al. 2017; Rocha, Johnson & Gray 2019; Viroulet et al. 2019; Edwards et al. 2021) can be quantitatively modelled with these rough-bed friction laws. Many of these problems require the additional inclusion of depth-averaged in-plane viscous terms in the avalanche equations (Gray & Edwards 2014; Baker et al. 2016a; Kanellopoulos 2021). These second-order gradient terms represent a singular perturbation to the equations and are necessary to predict (i) the cut-off frequency of roll waves (Forterre 2006; Gray & Edwards 2014), (ii) the downslope velocity profile across channels (Baker et al. 2016a), (iii) the height and width of self-channelised flows (Rocha et al. 2019) and (iv) segregation induced flow fingers (Woodhouse et al. 2012; Baker, Johnson & Gray 2016b). In this paper it is shown that the qualitative features of granular flow around obstacles on a rough bed, such as shock and levee formation, as well as the static pile and vacuum shape, can be quantitatively predicted using the same model.

The interaction of large-scale dense avalanches, such as debris flows, pyroclastic flows and snow avalanches, with obstacles in their path can lead to major changes in the speed and direction of the avalanche. Due to the hazardous nature of such avalanches, many attempts have been made to divert or control the flow to protect people and infrastructure using defences such as dams and barriers, or in some cases using natural features such as maintained forests, on common avalanche paths (see, e.g.  Cui, Gray & Jóhannesson 2007; Olschewski et al. 2012; Luong, Baker & Einav 2020). Historical examples of avalanche defences date back as far as the 17th century when a church in the sub-parish of Frauenkirch near Davos, after being destroyed by an avalanche in 1602, was rebuilt with a spaltkeil (i.e.  a wedge-shaped structure used to divert the avalanche around the building) (Gray et al. 2003). Further examples include a 4-m-high and 80-m-long avalanche-deflection wall which was built in Leukerbad around the same period (Pudasaini & Hutter 2007). Through time, mitigation strategies have become more varied. This can be observed on the Walmendinger Horn in the Kleinwalsertal/Mittelberg where arrays of avalanche protection barriers have been made using dry stone walls, gabions and earthworks that date from various points in the early 20th century (Drexel & Macher 2018). A more modern example of this tradition are the protective measures above Siglufjör$\eth$ur in northern Iceland which include multiple collection dams directly above the town and a series of barriers on the steeper slopes above (Arnalds et al. 2001), or at the Taconnaz site in Chamonix, France where series of breaking dams are used to slow the avalanche (Naaim, Faug & Naaim-Bouvet 2003). The features observed in such real-world flows are strikingly similar to those in dense dry granular avalanches indicating a strong correlation between the two (Sovilla et al. 2008; Köhler, McElwaine & Sovilla 2018), where large-scale experiments on snow avalanches impacting obstacles have been conducted with features showing many similarities to small-scale granular avalanches. This connection has influenced the development and study of avalanche protection design (see, e.g. Naaim et al. 2003; Cui et al. 2007; Barbolini et al. 2009). A recent example is that of the Cabane Arpitettaz, shown in figure 3, where boulders have been placed upstream of the cabin to divert oncoming avalanches around it. The shape of the protective structure is similar to the maximal static dead zone that would naturally form against a rigid obstacle, as will be shown in § 7. Overall this indicates that there is a long history of developing effective mitigation measures for avalanche hazards. This study extends our understanding of the avalanche-obstacle interaction to situations where the terrain has significantly more complicated frictional dependence than in previous studies.

Figure 3. Photograph of a snow avalanche defence built against the Cabane Arpitettaz near Zinal, Switzerland. The geometry of boulders built up against the structure resembles a maximal pile (§ 7), which is a stable shape that effectively diverts incoming snow away from the building.

2. Governing equations

Consider a slope inclined at an angle $\zeta$ to the horizontal and let $Oxyz$ be a Cartesian coordinate system with the $x$-axis orientated in the downslope direction, the $y$-axis in the cross-slope direction and the $z$-axis pointing in the upward normal direction to the slope, as shown in figure 4. The obstacle lies within this domain and is included in calculations through the applied boundary conditions, discussed further in § 4. In these coordinates the depth-averaged mass and momentum conservation laws for the avalanche thickness $h(x,y,t)$ and the depth-averaged velocity ${\boldsymbol{\bar{u}}}(x,y,t)=(\bar {u}(x,y,t),\bar {v}(x,y,t))$ are

(2.1)$$\begin{gather} \frac{\partial h}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} ( h {\boldsymbol{\bar{u}}})= 0, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial}{\partial t} ( h {\boldsymbol{\bar{u}}} ) + \boldsymbol{\nabla} \boldsymbol{\cdot} ( h {\boldsymbol{\bar{u}}} \otimes {\boldsymbol{\bar{u}}} ) + \boldsymbol{\nabla}\left( \frac{1}{2} g h^2 \cos \zeta \right) = g h {\boldsymbol{{S}}} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\nu h^{{3}/{2}} {\boldsymbol{\bar{D}}}), \end{gather}$$

where $\boldsymbol {\cdot }$, $\otimes$ and $\boldsymbol {\nabla } = (\partial / \partial x, \partial / \partial y)$ are the two-dimensional dot product, dyadic product and gradient operators, respectively. The shape factor has been implicitly assumed to be equal to unity, in keeping with many depth-averaged avalanche models (see, e.g. Hutter et al. 1993; Pouliquen & Forterre 2002). The mass equation (2.1) is derived under the approximation that the avalanche is incompressible. While this is not true of very rapid granular flows, the typical inertial numbers in the experiments presented here are small enough that the change in volume fraction $\phi (I)$ (Pouliquen et al. 2006) between static and flowing regions is relatively small, $\sim$5 %.

Figure 4. Experimental sketch of a granular flow (yellow) impacting a blunt obstacle (grey) on an inclined plane with a (turquoise) rough bed. The chute is inclined at an angle $\zeta$ to the horizontal and the flow is released from an upstream hopper. The flow is steady and uniform prior to impacting the obstacle, which is attached to the chute with neodymium magnets through the depth of the chute. The dashed line shows the approximate position of the bow shock for a supercritical inflow, while the solid lines on the chute show the approximate particle paths of individual grains as they flow around the obstacle.

The source term ${\boldsymbol{{S}}}$ includes contributions from gravity and the effective friction

(2.3)\begin{equation} {\boldsymbol{{S}}} = \sin \zeta {\boldsymbol{{i}}} - \mu \cos \zeta {\boldsymbol{{e}}}, \end{equation}

where $\mu$ is the friction coefficient and ${\boldsymbol{{i}}}$ is a unit vector in the downslope direction. The frictional force is oriented along the direction ${\boldsymbol{{e}}}$, which is defined as

(2.4)\begin{equation} {\boldsymbol{{e}}} = \begin{cases} \displaystyle \frac{{\boldsymbol{\bar{u}}}}{|{\boldsymbol{\bar{u}}}|}, & |{\boldsymbol{\bar{u}}}| > 0, \\ \displaystyle\frac{\tan\zeta {\boldsymbol{{i}} } - \boldsymbol{\nabla} h}{|\tan\zeta {\boldsymbol{{i}} } - \boldsymbol{\nabla} h| }, & |{\boldsymbol{\bar{u}}}| = 0. \end{cases} \end{equation}

It follows that when $|\bar {\boldsymbol{{u}}}| = ({\bar u}^2 + {\bar v}^2)^{1/2}$ is non-zero, the frictional force opposes the direction of motion, and when the material is stationary, it opposes the resultant force due to gravity and the depth-averaged pressure gradient. The final term on the right-hand side of (2.2) is the depth-averaged viscous term, where $\bar {\boldsymbol{{D}}} = ( \boldsymbol {\nabla }\bar {\boldsymbol{{u}}} + (\boldsymbol {\nabla }\bar {\boldsymbol{{u}}})^\textrm {T})/2$ is the depth-averaged strain-rate tensor and the coefficient $\nu$ determines the depth-averaged kinematic viscosity $\nu h^{1/2}/2$ (Gray & Edwards 2014; Baker et al. 2016a).

2.1. The effective friction law

The non-monotonic friction law describes the hysteretic behaviour of the grains and consists of dynamic, intermediate and static regimes (Edwards et al. 2019),

(2.5)\begin{equation} \mu(h,{{Fr}}) =\begin{cases} \mu_{D}, & {{Fr}} \geq \beta_*,\\ \mu_{I}, & 0 < {{Fr}} < \beta_*, \\ \mu_{S}, & {{Fr}} = 0, \end{cases} \end{equation}

which are defined by the functions

(2.6)$$\begin{align} \mu_D &= \mu_1 + \frac{\mu_2 - \mu_1}{1 + h\beta/(\mathscr{L}({{Fr}} + \varGamma))},\end{align}$$

(2.7)$$\begin{align}\mu_I &= \left(\frac{{{Fr}}}{\beta_*}\right)^\kappa\left(\mu_1 + \frac{\mu_2 - \mu_1}{1 + h\beta/(\mathscr{L}(\beta_{{\ast}} + \varGamma))} - \mu_3 - \frac{\mu_2 - \mu_1}{1 + h/\mathscr{L}} \right)\nonumber\\ &\quad + \mu_3 + \frac{\mu_2 - \mu_1}{1 + h/\mathscr{L}},\end{align}$$

(2.8)$$\begin{align}\mu_S &= \min\left(|\tan\zeta \boldsymbol{{i}} - \boldsymbol{\nabla} h|,\ \mu_3 + \frac{\mu_2 - \mu_1}{1 + h/\mathscr{L}}\right). \end{align}$$

Following Russell et al. (2019) and Edwards et al. (2019), the frictional transition in the intermediate regime is assumed to be linear, i.e. $\kappa =1$ in (2.7). Measurements are made to determine the material parameters for $d=300 - 400$ ${\rm \mu} $m masonry sand using the methodology described by Pouliquen & Forterre (2002) and Edwards et al. (2019). The parameters $\beta$, $\beta _*$ and $\varGamma$ are determined by fits to the empirical flow law

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

as shown in figure 5(a), while the friction angles $\mu _1 = \tan \zeta _1$, $\mu _2 = \tan \zeta _2$, $\mu _3 = \tan \zeta _3$, and the frictional-length scale $\mathscr {L}$ are found by experimental fits to the curves

(2.10a,b)\begin{align} h_{stop}(\zeta) =\mathscr{L}\left(\frac{\tan\zeta_2-\tan\zeta_1}{\tan\zeta-\tan\zeta_1}-1\right),\quad h_{start}(\zeta)=\mathscr{L}\left(\frac{\tan\zeta_2-\tan\zeta_1}{\tan\zeta-\tan\zeta_3}-1\right), \end{align}

as shown in figure 5(b). The minimum dynamic friction is assumed to occur at a constant Froude number ${ {Fr}}=\beta _*$, which by (2.9) implies that the corresponding transition thickness $h_*=(\beta _*+\varGamma )h_{stop}/\beta$ is a multiple of $h_{stop}(\zeta )$, as shown in figure 5(b). The experimental parameters for the friction law are summarised in table 1. While the quantitative results in this paper are specific to this friction law for masonry sand, qualitatively similar results are observed experimentally and predicted theoretically when glass spheres of diameter 125–160 ${\rm \mu} \text {m}$ are used, with quite different friction parameters $\zeta _1=21.27^\circ$, $\zeta _2=33.89^\circ$, $\zeta _3=25.3^\circ$, $\beta =0.143$, $\beta _*=0.19$, $\varGamma =0$ and $\mathscr {L}=0.2351$ mm.

Figure 5. Panel (a) shows the experimental measurements of the Froude number ${ {Fr}}$ as a function of $h/h_{stop}$. A linear best fit to this data determines the gradient $\beta$ and the intercept $\varGamma$. The lowest Froude number for which a steady uniform flow exists determines $\beta _*$ (horizontal yellow line). Panel (b) shows the experimental data for the deposit layer thickness $h_{stop}(\zeta )$ and the corresponding activation thickness $h_{start}(\zeta )$. The best-fit curves (black and blue lines) to the functional forms (2.10) determine the parameters $\zeta _1$, $\zeta _2$, $\zeta _3$ and $\mathscr {L}$. The transition between the intermediate and dynamic friction regimes occurs at $h_*$, which is shown in yellow.

Table 1. Experimentally determined friction law parameters for masonry sand of grain size $d$.

2.2. The parameter $\nu$ in the depth-averaged kinematic viscosity

The parameter $\nu$ in the depth-averaged kinematic viscosity was derived from the $\mu (I)$-rheology (GDR-MiDi 2004; Jop et al. 2006) by Gray & Edwards (2014), assuming the flow was in the dynamic frictional regime. It is assumed to apply to all the frictional regimes in the simulations presented here. The value of $\nu$ is a function of the slope angle and parameters that are already known from the rough-bed friction law,

(2.11)\begin{equation} \nu = \frac{2}{9} \frac{\mathscr{L} \sqrt{g}}{\beta} \frac{\sin \zeta}{\sqrt{\cos \zeta}} \left(\frac{\tan \zeta_2 - \tan \zeta}{\tan \zeta - \tan \zeta_1} \right), \end{equation}

and it is well defined provided $\zeta \in [\zeta _1,\zeta _2]$. The viscous terms are a singular perturbation to the theory, only becoming significant in rapid shear zones and shocks, which become smooth transitions rather than sharp jumps (see, e.g. Baker et al. 2016a; Rocha et al. 2019). The viscous terms fundamentally change the system's structure from hyperbolic (Savage & Hutter 1989; Gray et al. 2003; Hákonardóttir & Hogg 2005; Hogg et al. 2005; Cui et al. 2007; Johnson & Gray 2011; Cui & Gray 2013; Viroulet et al. 2017; Edwards et al. 2021) to hyperbolic-parabolic. The inclusion of the viscous terms is vital to properly model the formation of levees (Rocha et al. 2019). However, the system is still convection dominated, so understanding the hyperbolic problem is still useful. In one spatial dimension, the method of characteristics implies that the characteristics $\lambda =\bar u\pm (gh\cos \zeta )^{1/2}$. Assuming the downslope velocity is positive, and substituting for $\bar u$ from the definition of the Froude number (1.1) implies that $\lambda =(gh\cos \zeta )^{1/2}({ {Fr}}\pm 1)$. It follows that both characteristics will point downslope for supercritical flows with ${ {Fr}}>1$, while for subcritical flows (with ${ {Fr}}<1$), one characteristic propagates upslope and the other downslope. Normal shocks require information from two upstream characteristics and one downstream one, to determine the jump in thickness and velocity, as well as the speed to the shock. As a result, steady shocks are only observed for supercritical flows.

3. Small-scale experimental flows around a blunt obstacle

The experiments are carried out on an inclined plane which is $1.65$ m long and $0.58$ m wide. It is roughened with a monolayer of 750–1000 ${\rm \mu} $m turquoise spherical glass beads that are stuck to the surface using double-sided tape. The flows of 300–400 ${\rm \mu} $m soft masonry sand are fed onto the slope from a hopper with an adjustable gate, which allows the flux to be controlled by a combination of the gate height and slope angle. Flow thicknesses are measured using a laser profilometer (Micro-Epsilon scanCONTROL 2700 - 100) and depth-averaged velocities of steady uniform flows are determined by measuring the speed of a steadily propagating flow front. A $3\times 8\times 3$ cm metal obstacle is fixed to the base using strong neodymium magnets that are placed on either side of the chute. The obstacle is located $0.7$ m downslope from the hopper, and the results are insensitive to this position. Sequences of experimental photographs of the avalanches are captured using an overhead camera positioned perpendicular to the chute. Long exposure times and contrast enhancement are used to highlight the flowing and static regions.

3.1. Supercritical experiment on a rough inclined plane

The sequence of experimental photographs in figure 6 illustrates the behaviour of a typical supercritical flow for a gate height of $7$ mm and a slope angle of $35^\circ$. The first image (at $t=0$ s) shows the sand propagating downslope towards the obstacle from the hopper (figure 6a). A supercritical steady uniform flow with ${ {Fr}} \approx 2.35$ is rapidly established on the chute prior to impact with the obstacle. This contrasts with experiments on smooth beds, where the avalanche is usually accelerating, or decelerating, and the upstream Froude numbers are typically in the range ${ {Fr}}=3$–$5$ (Gray et al. 2003; Hákonardóttir & Hogg 2005; Gray & Cui 2007; Pudasaini et al. 2007; Cui & Gray 2013). As the avalanche hits the obstacle (figure 6b–d) a static dead zone forms on the front face and propagates upslope together with a detached bow shock which deflects the oncoming material around the dead zone and the obstacle. The bow shock and the static dead zone are close to their steady-state positions by $t=2.27$ s. The height and velocity of the flow change rapidly at the bow shock, which is approximately parabolic in shape and lies some distance upstream of the triangular static deposit. The change in height can be seen as a shadow in the photographs in figure 6, and the change in direction is indicated by the streak lines created by the long exposure. A close-up view of the upstream supercritical flow structure, with schematic annotations, is shown in figure 7(a). Figure 8(a) shows thickness contours taken from a laser scan of the upstream region for supercritical flow with the same inflow parameters. The surface of the flow during steady state was constructed from successive cross-slope laser thickness profiles taken at half-centimetre intervals in the $x$-direction. The parabolic shape of the bow shock is shown as a jump in thickness.

Figure 6. Sequence of overhead photographs for a gate height of $7$ mm on a slope angle of $35^\circ$, which generates a supercritical steady uniform flow with ${ {Fr}} \approx 2.35$. The downslope direction is from top to the bottom of each image. Panels (a–h) correspond to $t = 0$, 0.91, 2.27, 7.74, 8.65, 10.93, 14.57 and $20.49$ s, respectively. The exposure time is $0.3$ s and the photographs are contrast enhanced in order to highlight the moving and static areas of grains via motion blur. Panel (d) shows the flow in steady state, prior to the inflow being stopped and the formation of a static deposit (e–h). Movies showing the evolution, both with and without contrast enhancement, are available in the online supplementary material.

Figure 7. Annotated photographs showing the (a) supercritical and (b) subcritical steady-state flow structure upstream of the obstacle. Streamlines are shown in blue and arrows indicate the flow direction and relative velocity magnitude. The static dead zone is shaded bright yellow. In panel (a) the bow shock is shown as a parabolic red line. Upstream of it the flow is steady and uniform. As grains cross the shock, the flow rapidly changes in thickness and velocity, and is deflected outwards to flow around the static dead zone and the obstacle. In the subcritical flow shown in panel (b) a shock does not form and instead the changes are more gradual. The flow has four domains that are broadly arranged in a St Andrew's cross-like structure (red lines) that is centred at the apex of the dead zone. The dead zone lies in the lower quadrant of the cross. In the upper quadrant the oncoming steady uniform flow increases slightly in thickness and is decelerated (indicated by the decreasing length of the arrows). As the grains cross the red lines into the remaining two quadrants, the flow is rotated and accelerated outwards around the dead zone and the obstacle. Note that in both cases the dead zone is not perfectly triangular, but slightly scalloped.

Figure 8. Photographs with overlaid experimental thickness contours showing the (a) supercritical (${ {Fr}} = 2.35$, $\zeta = 35^\circ$) and (b) subcritical (${ {Fr}} = 0.77$, $\zeta = 32^\circ$) steady-state flow structure upstream of the obstacle. Panels (c) and (d) show the final deposits from the supercritical and subcritical inflows, respectively. Experimental thickness contours are calculated from a laser data scan, where successive cross-slope profiles at half-centimetre intervals in the $x$-direction are used to construct the surface of the flow. In panel (a) the bow shock is shown alongside the flat-topped pile structure, which lengthens and rounds in the final deposit shown in (c). In panel (b) the contours show a raised region upstream of the dead-zone apex in the shape of an hourglass indicating the ‘St Andrew's cross’ structure, described in figure 7. Panel (d) shows the final subcritical deposit, with a central ridge that extends far upstream.

The bow shock may also be viewed as a detached oblique shock (Gray & Cui 2007), since the static region acts like an obstacle itself. Here, however, the flow can interact with the dead zone, which may change its shape by eroding or depositing material. The static material supported by the obstacle is due to the multi-valued static friction balancing the pressure and the gravitational body forces. While the approximately triangular outline of the dead zone (when viewed from above) is consistent across repeated experiments, its thickness is not. Deposited granular material can collapse off the pile due to the material surface exceeding the angle of repose. In addition, small fluctuations in the flow can lead to the avalanche rolling over the static region and depositing further material onto the existing pile surface. In this case the pile remains fairly flat topped from its initial formation and does not reach the maximal size that can be supported on the inclined plane. Although not uniquely defined, the flat-topped nature of the pile structure can be observed in the experimental thickness contours in figure 8(a,c). Once the steady state has been established (figure 6d) the static region no longer changes shape or size. As the close-up in figure 7(a) shows, the steady-state dead zone is not perfectly triangular in shape, but slightly scalloped.

Figure 6(b,c) show that on the lee side of the obstacle the oncoming flow wraps around and forms a grain-free granular vacuum that reaches steady state by $t=7.74$ s (figure 6d). The granular vacuum is bounded by static levees (Edwards et al. 2017; Rocha et al. 2019), which are a new feature of the flow. This contrasts strongly with smooth beds, where the highest velocities are attained on the vacuum boundary, as the oncoming flow expands laterally into the space behind the obstacle. Here, as material is deflected by the bow shock and the dead zone, it forms two jets of diverted material on either side of the obstacle that are raised in height. As these pass the rear face, the flow expands to fill the space just like on smooth beds. However, when the thickness of the flow falls below $h_{stop}$ it has a strong propensity to deposit, unless it is driven downslope by gradients of thickness or velocity. At the levee-flow boundary, both the velocity and the viscous shear stress are zero at steady state (Rocha et al. 2019). The levees start from the rear face of the obstacle and continue for the entire length of the grain-free region. They curve slightly inwards towards the centre of the plane before intersecting again further down the slope forming a cusp. Downstream of this point the combined levees continue to exist until far downstream, when erosion eventually eradicates them and the avalanche tends back to a steady uniform flow. Note, that sufficiently far away in the $y$-direction, the obstacle does not influence the avalanche at all, and it remains at steady uniform flow down the entire chute.

After $t=7.74$ s the material supply ceases, and the flow thins and decelerates. This runout period is shown in figure 6(e–h). As the material runs out, the dead zone and levees relax and reach the final static deposit on the plane. For steeper slope angles, the collapsed material from the dead zone and the levees forms small erosion-deposition waves which propagate off the inclined plane in discrete wave pulses (Edwards & Gray 2015; Edwards et al. 2017; Rocha et al. 2019; Viroulet et al. 2019; Edwards et al. 2021). By $t=20.49$ s the flow comes to a rest and forms a final deposit. This consists of a thin layer of material close to $h_{stop}$ over most of the plane (figure 6h), but with remnants of the dead zone forming a static pile upstream of the obstacle and raised levees remaining in situ on either side of the grain-free granular vacuum on the lee side. Laser thickness contours of a final deposit are shown in figure 8(c). Comparison between figures 8(a) and 8(c) shows the shortening and rounding of the dead-zone structure as it relaxes to the final deposit.

3.2. Subcritical experiment on a rough inclined plane

As opposed to smooth beds, where subcritical flows only occur briefly as the avalanche decelerates in the runout zone, rough beds allow the study of sustained subcritical flows. Figure 9 shows a sequence of experimental images for a 7 mm gate height on a $32^\circ$ slope. A steady uniform flow with ${ {Fr}} \approx 0.77$ rapidly establishes itself and propagates downslope as shown in figure 9(a). As the material hits the front face of the obstacle, grains are deposited and an interface between the flowing and static material propagates upslope. A flat-topped triangular static dead zone is formed (figure 9b–d), which is slightly longer than before, but no bow shock is formed. Instead, a region of upstream influence develops, which can be more clearly seen in figures 1(b) and 8(b) and in the supplementary movies (available at https://doi.org/10.1017/jfm.2021.1074). At steady state the flow develops four regions that form a St Andrew's cross structure that is centred at the apex of the dead zone, as shown schematically in figure 7(b) and as a raised region in the experimental thickness contours. The static dead zone is located in the lower quadrant, while in the upper quadrant the oncoming material from the hopper is slowed down and increases slightly in height. In the two domains on either side, the flow deflects to flow parallel to the dead zone, and then accelerates around the obstacle. The regions of acceleration and deceleration are localised just upstream of the obstacle and do not extend back to the hopper, as shown in the online supplementary movie. The St Andrew's cross structure is a novel feature of the flow that has not been reported before.

Figure 9. Sequence of overhead photographs for a gate height of $7$ mm on a slope angle of $32^\circ$, which generates a subcritical steady uniform flow with ${ {Fr}} \approx 0.77$. The downslope direction is from top to the bottom of each image. Panels (a–h) span the times $t = 0$, 0.91, 1.83, 4.56, 8.2, 10.94, 14.58 and $23.69$ s, respectively. Panel (e) is in steady state, prior to the flow waning ( f,g) and forming a static deposit (h). Movies showing the evolution, both with and without contrast enhancement, are available in the online supplementary material.

On the lee side of the obstacle, the flow expansion and the formation of a grain-free vacuum region bounded by static levees are very much like the supercritical case. Here, however, the levees are considerably wider and thicker. This is partly because $h_{stop}$ and $h_{start}$ are larger on the lower inclination slope, and partly because the flow is slower and thicker than the supercritical case, so the levees need to be more substantial to prevent lateral spreading. The overall shape of the vacuum region is similar to the supercritical case, although the cusp is slightly further downstream. Figure 9( f–h) show the runoff of material and the formation of a final static deposit after the inflow ceases. At this slope angle, discrete erosion-deposition waves do not form and instead a deposition front simply propagates downstream. This is analogous to the travelling deposition wave described by Edwards et al. (2019). As the region of upstream influence relaxes, it largely drains away, leaving only a ridge that extends beyond the static pile left by the dead zone (see figures 1b, 2 and 8d). The levees also relax slightly and lose some material as the deposition wave propagates through the system. Over the remaining parts of the chute that were occupied by flowing material, the $h_{stop}$ layer is slightly thicker than the supercritical case, because of the reduced slope angle.

4. Numerical method

To investigate the system further, the depth-averaged equations (2.1)–(2.2) are solved with a high-resolution semi-discrete non-oscillatory central scheme for convection–diffusion equations (Kurganov & Tadmor 2000). In particular, this scheme can handle the non-standard viscous dissipation terms in the theory, and has been extensively tested against exact solutions and experiments (Edwards & Gray 2015; Baker et al. 2016b; Edwards et al. 2017; Rocha et al. 2019). A well-balanced discretisation is used for the basal friction source terms (2.4) and (2.8) so that, in regions of static grains that are below yield, the friction exactly balances the other forces and thereby keeps the material stationary.

The numerical domain $(L_x,L_y) = (2,0.6)$ m is discretised into a rectangular grid with $(nx,ny) = (1200,360)$ grid points, respectively. The initial conditions specify a small stationary precursor layer, of thickness 0.1 mm, over the domain, which ensures stability of the simulations by counteracting numerical errors caused by the degeneracy in the equations when the thickness approaches zero. Two constraints determine the numerical time step. For stable integration of the hyperbolic equations, a maximum time step given by a CFL (Courant–Friedrichs–Lewy) condition is used, with a CFL number of $0.25$ (LeVeque 2002). Accurate integration of the source terms in static regions requires a time step no greater than $1\times 10^{-4}$ s. The smaller of these two values is chosen to be the numerical time step. This maximal step size is required to minimise the creep observed for both the pile above the obstacle and the levees. For the first part of the simulation, a steady uniform flow of a specified Froude number is prescribed at the inflow boundary at $x = 0$. Periodic boundary conditions are applied on the sides of the domain at $y = -0.3$ and $y = 0.3$ and outflow conditions are specified at $x = 2.0$ m. A no-penetration boundary condition $\bar {\boldsymbol{{u}}} \boldsymbol {\cdot } \boldsymbol{{n}}=0$ (where $\boldsymbol{{n}}$ is the normal vector to the obstacle boundary) is applied along the boundaries of the $3 \times 8$ cm obstacle, which is specified within the domain. Ghost cells at the obstacle boundary are used to facilitate this by applying matching conditions for the thickness at the boundary and specifying the required fluxes. This simulates an infinitely tall impermeable obstacle within the numerical domain. To mimic the drainage regime in the experiments, the height and downslope flux parameters at the inflow are reduced exponentially over $0.5$ s after $t=45$ s. The exponential reduction is a simple representation of the waning flux as the last material leaves the hopper. The sustained flux up to $t=45$ s allows the flow to establish a steady state prior to transitioning to the drainage regime. The resulting stopping wave propagates downstream over the domain, after which all material is static. The final deposits can then be analysed and compared with the experimental counterparts.

5. Numerical simulations of the supercritical flow and deposit

The numerical method allows the thickness $h$, as well as the velocity components $\bar u$ and $\bar v$, to be computed for a Froude number ${ {Fr}} = 2.35$ flow impacting the obstacle on a $35^\circ$ slope. This corresponds to the same conditions as those in the supercritical experiment in § 3.1. The time-dependent behaviour of each of the fields is shown in figures 10–12. Streamlines calculated from the instantaneous velocity components, $\bar {u}$ and $\bar {v}$, are overlaid on the downslope velocity plots and give an indication of the direction of the flow of material. In addition, there is a movie available online, which shows the simultaneous time-dependent evolution of all the fields, together with the Froude number ${ {Fr}}$.

Figure 10. Numerical simulations of the avalanche thickness $h$ on a $35^\circ$ slope and with ${ {Fr}} = 2.35$. The relative times $t = 0$, 0.9, 2.3 and $7.7$ s in panels (a–d) match the experiment in figure 6. The drainage sequence is triggered at $45$ s. In panels (e–h) the times match features from the experiments at $t = 44.7$, $44.9$, $50.3$ and $58.3$ s, respectively. The $3\times 8$ cm obstacle is highlighted in grey and is located in the same position as in the experiment. Grain-free regions, where the thickness of the flow is below a grain diameter, are shown in purple. A movie showing the time-dependent evolution is available in the supplementary material online.

The sequence of plots in figures 10–12 exhibit similar features to the experiments. The first panel in each plot (figures 10–12a) shows a steadily propagating front approaching the obstacle. Figures 10–12(b–d) show the transient formation of the bow shock, dead zone, vacuum region and static levees. Material impacts the front face of the obstacle and comes to rest to form a dead zone, while a curved bow shock forms upstream. This structure then propagates upslope until it reaches steady state. The thicker slower moving deflected material after the shock accelerates and forms a jet as it flows around either side of the dead zone and the obstacle. The dead zone is wedge shaped and has a flat top, which is qualitatively similar to the experiments in both the method of formation and the final shape. The agreement in the time scales for the formation of the features is good, and the simulation reaches a steady state in a similar length of time. Likewise, the observation that the shock is detached from the obstacle and the dead zone is similar to the experimental case.

Figure 11. Numerical simulations of the depth-averaged downslope velocity $\bar u$ on a $35^\circ$ slope and with ${ {Fr}} = 2.35$. The times are the same as in figure 10. Streamlines, calculated from the instantaneous velocity fields, are overlaid in black. Purple is used to highlight grain-free regions and the obstacle is shown in grey. A movie showing the time-dependent evolution is available in the supplementary material online.

Figure 12. Numerical simulations of the depth-averaged cross-slope velocity $\bar v$ on a $35^\circ$ slope and with ${ {Fr}} = 2.35$. The times are the same as in figure 10. The speckled colouring is observed within the deposit regions, due to a small residual creep in the numerical method and the small threshold for the contour scale. Purple is used to highlight grain-free regions and the obstacle is shown in grey. A movie showing the time-dependent evolution is available in the supplementary material online.

As the grains spread out into the space on the lee side of the obstacle, the depth-averaged pressure causes the flow to expand laterally and form a grain-free granular vacuum. The friction law (2.5) implies that material which has a thickness below $h_{start}$ can come to rest, and static levees, therefore, form adjacent to the grain-free vacuum region. Similar to the experiments, the levees curve inwards before intersecting far downstream. Once the levees have combined, the intersection point retreats upstream, shrinking the vacuum region until a steady state is reached. The levees persist downstream of the grain-free region and are only slowly eroded with increasing downstream distance. By $t=7.7$ s the flow is close to steady state (figures 10–12d). Small fluctuations do, however, imply that there can be some residual lateral motion in the dead zone. These small fluctuations trigger the formation of subtle waves further downstream, which can be clearly seen in the cross-slope velocity field in the online supplementary movie. These are likely to develop into roll waves on a long enough chute (Forterre & Pouliquen 2003; Gray & Edwards 2014; Edwards & Gray 2015; Viroulet et al. 2018). Unlike flows on smooth beds, the online movie indicates that the Froude number is close to its steady uniform value over the entire domain, except immediately downstream of the shock, in the dead zone and in the levees.

The final sequence of figures 10–12( f–g) show the drainage as the inflow stops at $t=45$ s. The material is brought to rest via a stopping wave that propagates downstream and leaves a deposit that has a thickness close to $h_{stop}$ (Edwards et al. 2019). As the depth-averaged pressure diminishes, the dead zone and the levees relax and reduce in size. The material that is activated in this process creates finite pulses of flowing material next to the levees. These are similar to a finite release of material onto an erodible bed, which due to the high slope angle breaks into discrete erosion-deposition waves (see Edwards et al. 2017). This is captured in figures 10–12(g), and is particularly nicely visualised in contours of $\bar v$ at $t=50 - 60$ s in the online supplementary movie. Note that, as the material comes to rest, the shape of the pile formed by the collapsing dead zone changes, allowing material to be deposited further upstream. This in turn results in a longer pile shape. Figures 10–12(h) show the final deposit where the whole domain is considered to be static. Here the levees bounding the vacuum region are raised and the static pile upstream of the obstacle is reduced in size and becomes more rounded. This corresponds to the experimental case where the pile shape relaxes and elongates.

Figure 13 shows a three-dimensional reconstruction of the steady state and final deposit, as well as cross-slope thickness profiles, from both the simulation and the corresponding experiment, at fixed intervals downslope to show the shapes of the features in more detail. The experimental cross-slope laser profiles shown are downsampled in order to remove noise and highlight the shape of the features. An animation of the three-dimensional reconstruction is available in the online supplementary material. In particular, this shows that no waves are able to propagate upslope past the bow shock. At steady state (figure 13c,e,g,i), the bow shock has a diffuse profile and the static dead zone has a distinct wedge shape with a flat top. On the lee side, the levees are most pronounced close to the obstacle and diminish in height with increasing downstream distance. The profiles of the final deposit in figure 13(d, f,h,j) show that there is a significant decrease in the pile and levee heights, and there is evidence of mass shedding during the propagation of the discrete erosion-deposition waves (Viroulet et al. 2019; Edwards et al. 2021). The final shape of the collapsed pile is much more rounded than that in the experiment, which maintains the flat top. Figure 14 shows a direct comparison between the thickness contours and the shape of the experimental features on the $xy$-plane, for both the steady state and the final deposit. The main difference between the steady-state simulation and experimental cases is the standoff distance of the shock, which is under-predicted by the simulation. In terms of the qualitative features, the prediction is good for both the dead-zone shape and size at steady state as well as its collapsed form in the final deposit. The computed vacuum region and the static levees are in excellent quantitative agreement with the experiments.

Figure 13. Three-dimensional surface reconstruction for (a) a steady-state supercritical flow with ${ {Fr}} = 2.35$ on a $35^\circ$ slope at $t = 45$ s and (b) the deposit at $t = 60$ s. Lighting effects are added to highlight the shape of the features. The cross-slope profiles above the obstacle (c)–( f) are brightness graded in green and correspond to the lines in panels (a) and (b), respectively. Panels (e) and ( f) show experimental laser thickness profiles at the same upslope intervals (and corresponding colours), under the same inflow and slope angle conditions. The results are downsampled to remove noise. The red graded profiles in panels (g) and (h) are the corresponding cross-sections below the obstacle, and likewise panels (i) and (j) show the corresponding laser thickness profiles. The figures show the rounding of the flat-topped pile upstream of the obstacle and the raised levees in the deposit. A movie showing the time-dependent evolution is available in the supplementary material online.

Figure 14. Comparison of the computed thickness contours with the flow features in the experimental photographs for supercritical flow with ${ {Fr}}=2.35$ on a $35^\circ$ slope. Panel (a) shows the steady state at $t = 7.74$ s and (b) shows the final deposit. The purple line indicates the edge of the vacuum region in the simulation, where the interior is less than a grain diameter.

6. Numerical simulations of the subcritical flow and deposit

Figures 15–17 show a subcritical flow with ${ {Fr}}=0.77$ impacting the obstacle on a $\zeta = 32^\circ$ slope, which are the same conditions as in the experiments in § 3.2. In figures 15–17(a) the pre-impact front is steadily propagating and laterally uniform as in the experimental case. The steady-state dead zone and the vacuum region then form in figures 15–17(b–e). As in the supercritical case, the material comes to rest on the front face of the obstacle and a deposition front propagates upslope until an upstream wedge-shaped flat-topped static dead zone is formed. Its size and time scale for formation are in good agreement with the experiments in figure 9, and it extends slightly further upstream than the supercritical case. For subcritical flows, there is no bow shock to slow and deflect the flow. Instead, the flow decelerates upstream of the pile. This is consistent with the St Andrew's cross structure observed in experiments and shown in the annotated close-up photographs in figure 7(b). The region of upstream influence is most easily seen in the cross-slope velocity plots in figure 17(b–e) and the associated movie (available online). A deflection wave propagates outwards and upslope, and attains a maximum steady-state upstream position $x\approx 0.4$ m. This lies significantly higher upslope than the apparent changes in thickness and downslope velocity. Nevertheless, this wave can also be seen propagating upslope in the movie of the three-dimensional reconstructed thickness (see figure 18). This wave diminishes in height with increasing upstream distance, so there is no influence on the hopper outflow, and, hence, the results are independent of the obstacle location, provided it is far enough away from the hopper.

Figure 15. Numerical simulations of the flow thickness $h$ on a $32^\circ$ slope and with ${ {Fr}} = 0.77$. Panels (a–e) match their experimental counterparts in figure 9 at the relative times $t = 0$, 0.9, 1.8, 4.6 and $8.2$ s. The drainage sequence is triggered $45$ s after the simulation is initiated. In panels ( f–h) the times are chosen to match features from the experiments at $t = 46.9$, $47.7$ and $50.5$ s, respectively. Purple is used to highlight grain-free regions and the obstacle is shown in grey. Contour lines are overlaid using a darker colourmap to highlight the shapes of the features. A movie showing the time-dependent evolution is available in the supplementary material online.

Figure 16. Numerical simulations of the depth-averaged downslope velocity $\bar u$ on a $32^\circ$ slope and with ${ {Fr}} = 0.77$. The times are the same as in figure 15. Streamlines, calculated from the instantaneous velocity fields, are overlaid in black. Purple is used to highlight grain-free regions and the obstacle is shown in grey. A movie showing the time-dependent evolution is available in the supplementary material online.

Figure 17. Numerical simulations of the depth-averaged cross-slope velocity $\bar v$ on a $32^\circ$ slope and with ${ {Fr}} = 0.77$. The times are the same as in figure 15. The speckled colouring is observed within the deposit regions, due to a small residual creep in the numerical method and the small threshold for the colour scale. Purple is used to highlight grain-free regions and the obstacle is shown in grey. A movie showing the time-dependent evolution is available in the supplementary material online.

Figure 18. Three-dimensional surface reconstruction for (a) a steady-state subcritical flow with ${ {Fr}} = 0.77$ on a $32^\circ$ slope at $t = 45$ s and (b) the deposit at $t = 60$ s. Lighting effects are added to highlight the shape of the features. The cross-slope profiles above the obstacle (c) and (d) are brightness graded in green and correspond to the lines in panels (a) and (b), respectively. Panels (e) and ( f) show experimental laser thickness profiles at the same upslope intervals (and corresponding colours), under the same inflow and slope angle conditions. The results are downsampled to remove noise. The red graded profiles in panels (g) and (h) are the corresponding cross-sections below the obstacle, and likewise panels (i) and (j) show the corresponding laser thickness profiles. The figures show the rounding of the flat-topped pile upstream of the obstacle and the raised levees in the deposit. A movie showing the time-dependent evolution is available in the supplementary material online.

Sufficiently far away from the obstacle in the $y$-direction the flow remains steady and uniform. As a result, the material that was deflected to the sides, upstream of the obstacle, flows past it with a slightly increased thickness and velocity. Once past the obstacle, the material expands into the free space and forms a granular vacuum that is bounded by a static levee, in the same way as in the supercritical case, but the levee is much more substantial in height and width. The levees curve inwards and intersect downstream resulting in a shape, which is in good agreement with the experiment. Figures 15–17(e) show the steady-state solution close to the obstacle, although the front continues to propagate downslope as shown in the movie online. This time, no roll waves are triggered, although in principle they could develop. The subsequent figures 15–17 show the drainage of the material leading to the final deposit. A deposition wave (Edwards et al. 2019) propagates down the chute leaving a final deposit that is close to $h_{stop}(32^\circ )$ over much of the chute. Unlike in the experiments, the flat-topped static dead zone is not perfectly preserved in the final deposit, but is progressively eroded from the sides, shrinking it in size, as shown in figure 18(d, f). As the thicker region of upstream influence also relaxes, it leaves a ridge upstream of the collapsed dead zone, just like the experiment in figures 2 and 8. The grains that are released from the region of upstream influence and the dead zone form two strips of flowing material on either side of the obstacle that pinch off from the main part of the draining flow at $t=49$ s, and form secondary mass shedding waves (Viroulet et al. 2019; Edwards et al. 2021) that deposit adjacent to the levees, which themselves are well preserved in the deposit. The deposition sequence is particularly well captured in the three-dimensional thickness reconstruction movie in the online supplementary material. Figure 19 shows a direct overlay of the thickness contours onto the experimental photographs for the subcritical flow in steady state and for its final deposit. The length of the computed steady-state dead zone is greater than that in the experiment, however, the length of the final static pile is in good agreement. The shape and size of the levees are also generally in good agreement with one another, although the cusp where the two levees meet is slightly further downstream in the experiments.

Figure 19. Comparison of the computed thickness contours with the flow features in the experimental photographs for subcritical flow with ${ {Fr}}=0.77$ on a $32^\circ$ slope. Panel (a) shows the steady state at $t = 8.2$ s and (b) shows the final deposit. The purple line indicates the edge of the vacuum region in the simulation, where the interior is less than a grain diameter.

7. The shape of the static dead zone

The dead-zone length is defined as the distance from the tip of the static dead zone to the front face of the obstacle measured along the downstream coordinate $x$. Figure 20(a) shows the numerically simulated dead-zone length for a wide range of slope angles, obstacle widths and Froude numbers. At low slope angles, the maximum pile length is obtained at low Froude numbers, and the length reduces substantially as the flow depth and, hence, the Froude number is increased. As the slope angle is raised, this effect is progressively reduced until by $\zeta =37^\circ$ the pile length is relatively insensitive to the Froude number. In order to explain this behaviour, it is instructive to consider simplified solutions for the maximal static pile that can be supported by the obstacle.

Figure 20. The numerical (a) dead-zone length, (c) stagnation-point thickness and (e) corner thickness as a function of the inclination angle, obstacle width and Froude number. Panels (b), (d) and ( f) show the comparison to the approximate theoretical solution. The value of $\mu _S=\mu _3$. The values from laboratory experiments are shown with the marker indicated by $w_{\text {e}}$.

7.1. Eikonal theory for the shape of a static pile on an empty plane

Let $OXYZ$ be a rectangular Cartesian coordinate system with the origin $O$ lying at the base of the middle of the upstream face of the obstacle. The $Z$-axis is aligned with gravity, the horizontal $X$-axis lies in the same plane as the $x$-axis and the $Y$-axis points across the slope in the same direction as the $y$-axis. The height $Z=H(X,Y)$ of the maximal pile satisfies the Eikonal equation (Hadeler & Kuttler 1999; Pauli & Gioia 2007; Colombo, Guerra & Monti 2012)

(7.1)\begin{equation} |\boldsymbol{\nabla} H| = \sqrt{\left(\frac{\partial H}{\partial X}\right)^2 + \left(\frac{\partial H}{\partial Y}\right)^2 } = \mu_S, \end{equation}

where the static friction $\mu _S$ is assumed constant for simplicity. The Eikonal equation can conveniently be written in the form

(7.2a–c)\begin{equation} F=\frac{1}{2}(p^2+q^2-\mu^2_S), \quad\text{where}\ p=\frac{\partial H}{\partial X}\quad \text{and}\quad q=\frac{\partial H}{\partial Y}, \end{equation}

and solved using Charpit's method (see, e.g. Sneddon 1957; Stavroulakis & Tersian 2004). This solution is obtained along characteristic curves $(X(\tau ),Y(\tau ))$ that are parameterised by $\tau$, and satisfy the ordinary differential equations (ODEs)

(7.3a,b)\begin{equation} \frac{\textrm{d}X}{\textrm{d}\tau}=\frac{\partial F}{\partial p}=p,\quad \frac{\textrm{d}Y}{\textrm{d}\tau}=\frac{\partial F}{\partial q}=q. \end{equation}

The values of $p$ and $q$ satisfy the ODEs

(7.4a,b)\begin{equation} \frac{\textrm{d}p}{\textrm{d}\tau}={-}\frac{\partial F}{\partial X}-p\frac{\partial F}{\partial H}=0,\quad \frac{\textrm{d}q}{\textrm{d}\tau}={-}\frac{\partial F}{\partial Y}-q\frac{\partial F}{\partial H}=0, \end{equation}

and the free-surface height satisfies

(7.5)\begin{equation} \frac{\textrm{d}H}{\textrm{d}\tau}=p\frac{\partial F}{\partial p}+q\frac{\partial F}{\partial q}=p^2+q^2=\mu_S^2. \end{equation}

The system of five (7.3)–(7.5) are known as Charpit's ODEs, and are subject to the initial conditions

(7.6a–e)\begin{equation} X=X_0(s),\quad Y=Y_0(s),\quad p=p_0(s),\quad q=q_0(s),\quad H=H_0(s), \end{equation}

defined along a curve parameterised by $s$ and where $p_0^2(s)+q_0^2(s)=\mu _S^2$ by (7.2a–c). In classical Eikonal problems (Stavroulakis & Tersian 2004; Pauli & Gioia 2007) an auxiliary condition is used to determine the boundary conditions for $p_0(s)$ and $q_0(s)$ on a boundary prescribed by $x_0(s)$ and $y_0(s)$. Here the pile selects its own boundary, which needs to be solved for as part of the problem.

For an empty chute and a constant value of the static friction, the maximal static pile solution is generated by two expansion fans centred at $X=0$ and $Y=\pm w/2$, where $w$ is the width of the obstacle. The obstacle height is implicitly assumed to be sufficient to support the static pile. Since equations (7.4) are equal to zero, it follows that $p$ and $q$ are equal to their initial values

(7.7a,b)\begin{equation} p=p_0,\quad q=q_0={\mp}\sqrt{\mu_S^2-p_0^2}, \end{equation}

where $p_0$ now parameterises the characteristic fan. This in turn allows (7.3) and (7.5) to be integrated, subject to the initial conditions (7.6) to give the solution in parametric form

(7.8a–c)\begin{equation} X=p_0\tau,\quad Y={\mp}\sqrt{\mu_S^2-p_0^2}\,\tau\pm\frac{w}{2},\quad H=\mu_S^2\tau. \end{equation}

The outermost characteristics of the expansion fans are determined by the intersection of the free-surface height (7.8c) and the inclined plane

(7.9)\begin{equation} Z={-}\mu X. \end{equation}

Substituting for $X$, from (7.8a), and cancelling $\tau$, implies the minimum $p_0=-\mu _S^2/\mu$. Similarly, the innermost characteristic is given by the characteristic that intersects with the obstacle along

(7.10)\begin{equation} Z=X/\mu, \end{equation}

which implies that the maximum value of $p_0=\mu _S^2\mu$. The boundaries of the pile are therefore characteristics. The static pile is shown in figure 21(a,d) and is parameterised by characteristics $p_0\in [-\mu _S^2/\mu,\mu _S^2\mu ]$, which emanate from $X=0$, $Y=\pm w/2$. The characteristics from each fan intersect along the pile ridge line, which lies along $Y=0$. The characteristics are therefore parameterised by $\tau \in [0,w/(2\sqrt {\mu _S^2-p_0^2})]$. Eliminating $\tau$ between (7.8a,b) implies that the characteristics are straight lines,

(7.11)\begin{equation} Y={\mp}\sqrt{\mu_S^2-p_0^2}\,\frac{X}{p_0}\pm\frac{w}{2}, \end{equation}

and the projections of these lines onto the free surface are also straight lines whose gradients are maximal, by construction. Eliminating $p_0$ and $\tau$ in (7.8), the free-surface height $H$ can be written as a function of $X$ and $Y$, i.e.

(7.12)\begin{equation} H=\mu_S \sqrt{X^2+(|Y|-w/2)^2}. \end{equation}

Figure 21(a,d) shows a reconstruction of the maximal static pile supported by the obstacle on an empty chute assuming that $\mu _S=\mu _3$ and mapped back to the $Oxyz$ coordinate system. During the experiments, the dead zone does not achieve its maximal normal height, but typically has a flat top that is roughly parallel to the plane. This truncation does not affect the characteristics that generate the solution up to this point, provided the truncated surface provides enough support from the obstacle to keep the pile static. The simplified model is therefore still relevant for calculating the shape of the dead zone.

Figure 21. Three-dimensional reconstruction (a) and overhead view (d) in $Oxyz$ coordinates of the maximal static pile solution (light yellow) for a grey obstacle of width $w=8$ cm on a $32^\circ$ slope (mauve). The black solid lines are of maximal slope $|\boldsymbol {\nabla } H|=\mu _S=\mu _3=0.62032$, the blue lines show the apex of the pile and the red lines show the intersection with the inclined plane. Panels (b,c,e, f) show the modified solution when the static pile is submerged beneath a layer of thickness $h_{stag}$ and the effective width of the obstacle $w_{effective}$ is increased, so that the static pile has thickness $h_{stag}$ at the corners of the obstacle at $x=x_{corner}$ and $y=\pm w/2$. Panels (b,e) show the shape of the collapsed static pile, when there is a static stabilizing layer of thickness $h_{stag}=h_{stop}=1.68$ mm and $w_{effective}=8.4$ cm, whereas panels (c, f) show the dead-zone shape when it is submerged by a layer of thickness $h_{stag} = 5.6$ mm and $w_{effective} = 9.4$ cm. This submergence depth is equivalent to an upstream flow of thickness $h_0=4.3$ mm and ${ {Fr}}_0 = 0.77$.

7.2. Partial submergence of the pile and increased width

During flow, the static pile is partially submerged by flowing material, which provides additional support at the sides and modifies its shape. To model this effect it is necessary to make a good approximation for the run-up height onto the pile. Assuming that the obstacle is impacted by a steady uniform flow of thickness $h_0$ and depth-averaged velocity $\bar u_0$, the furthest upstream point where the velocity is zero, is a stagnation point, and it lies on the central ridge of the static pile. The stagnation-point height $h_{stag}$ determines how much of the static pile is submerged by the flow. For supercritical flow, the steady uniform flow on the central streamline first passes through a normal shock, which is approximated by the steady-state inviscid mass and momentum jump conditions (see, e.g. Chadwick 1976; Cui & Gray 2013)

(7.13a,b)\begin{equation} [\![ h\bar u]\!]=0,\quad [\![ h\bar u^2+\tfrac{1}{2}gh^2\cos\zeta]\!]=0, \end{equation}

where the jump brackets $[\![ \chi ]\!]=\chi _1-\chi _0$ are the difference of the enclosed quantity on the forward and rearward sides of the shock. These imply that the downstream thickness and velocity are

(7.14a,b)\begin{equation} h_1=\frac{h_0}{2}(\sqrt{1+8{{Fr}}_0^2}-1),\quad \bar u_1=\frac{h_0\bar u_0}{h_1}, \end{equation}

where ${ {Fr}}_0$ is the upstream Froude number. Note that when ${ {Fr}}_0$ equals unity, the thickness $h_1=h_0$. On the central streamline, downstream of the shock the lateral velocity $\bar v$ is zero. Using the mass balance (2.1) it follows that the inviscid steady-state downstream momentum balance (2.2) reduces to

(7.15)\begin{equation} \frac{\partial}{\partial x}\left(\frac{\bar u^2}{2}+hg\cos\zeta\right)=g\cos\zeta(\tan\zeta-\mu). \end{equation}

Assuming that the source term on the right-hand side of (7.15) can be neglected, (7.15) can be integrated to relate the stagnation-point height to the downstream conditions after the shock, i.e.

(7.16)\begin{equation} h_{stag}=h_1+\frac{\bar u_1^2}{2g\cos\zeta}. \end{equation}

For subcritical flows, (7.16) also gives an approximation for the stagnation-point thickness, but since there is no shock $h_1=h_0$ and $\bar u_1=\bar u_0$. Figure 20(c) shows the simulated stagnation-point height as a function of the inclination angle and Froude number. These heights are very closely modelled by the approximate relation (7.16) for $h_{stag}$ as shown in figure 20(d).

The height that the flow runs up on either side of the static pile is also close to $h_{stag}$ all the way along the dead-zone/flow boundary. In particular, the run-up height as the flow passes the upstream corners of the obstacle $h_{corner}$ is well approximated by $h_{stag}$, as shown in figure 20(e, f). The static pile height at the obstacle corner is therefore not zero (as assumed in the static pile solution on an empty plane in § 7.1), but attains a finite thickness $h_{stag}$ when material is flowing around the obstacle. In order to approximate the length of the dead zone, it is necessary to modify the exact solution (7.12) for the shape of the dead zone, by replacing the width of the obstacle $w$ with its effective width $w_{effective}$, i.e.

(7.17)\begin{equation} H=\mu_S \sqrt{X^2+(|Y|-w_{effective}/2)^2}. \end{equation}

An expression for $w_{effective}$ can be obtained by assuming that $H=H_{stag}=h_{stag}/\cos \zeta$ on the edges of the obstacle, which by (7.10) lie at $X_{corner}=\mu H_{stag}$ and $Y=\pm w/2$. Substituting these expressions into (7.17) implies that

(7.18)\begin{equation} w_{effective}=w+2 H_{stag}\sqrt{\frac{1}{\mu_S^2}-\mu^2}. \end{equation}

A simple prediction of the resulting dead-zone shape can be obtained by intersecting its height (7.17) with an inclined plane that is parallel to the base and offset by a vertical distance $H_{stag}$, i.e. the plane

(7.19)\begin{equation} Z={-}\mu X+H_{stag}. \end{equation}

The resulting shape is given by

(7.20)\begin{equation} Y={\pm}\frac{w_{effective}}{2}\mp\sqrt{\left(\frac{-\mu X+H_{stag}}{\mu_S}\right)^2-X^2}. \end{equation}

For an empty chute, when $H_{stag}=0$, this is simply a triangle, as shown in figure 21(d) in $Oxyz$ coordinates. As $H_{stag}$ increases in depth, this shape becomes progressively shorter and more scalloped (figure 21e, f).

The furthest upstream point of the dead zone is obtained by intersecting the modified free-surface height (7.17) with the submerging plane (7.19) along the $Y=0$ axis, and solving the resulting quadratic for

(7.21)\begin{equation} X_{min}=\frac{2\mu H_{stag}-\sqrt{4\mu_S^2H_{stag}^2+(\mu^2-\mu_S^2)\mu_S^2w_{effective}^2}}{2(\mu^2-\mu_S^2)}. \end{equation}

The submerging plane (7.19) intersects the obstacle at

(7.22)\begin{equation} X_{max}=\frac{\mu}{1+\mu^2}H_{stag}. \end{equation}

It follows that the length of the pile in slope oriented coordinates is

(7.23)\begin{equation} L_{dead zone}=\frac{X_{max}-X_{min}}{\cos\zeta}. \end{equation}

Figure 20(b) shows that (7.23) does indeed provide a good collapse of the numerically computed dead-zone lengths for a wide range of Froude numbers, slope inclinations and obstacle widths. This simplified exact solution therefore provides considerable insight into what sets the shape and size of the dead zone.

8. Vacuum length

Immediately downstream of the obstacle, the grains flowing downslope either side of the obstacle expand inwards, eventually coming into contact with one another at the obstacle centreline and closing the vacuum region. This expansion, and consequently the length of the vacuum region, can be approximated using an asymptotic argument for the flow far downstream of the obstacle, similar to that in Hogg et al. (2005). This exploits the fact that, far downstream, variation in the downslope direction occurs over much larger length scales than in the cross-slope direction.

In this section the $Oxyz$ coordinate system defined in § 2 is used, but the origin is moved to the lower right corner of the downstream face of the obstacle, so that the obstacle lies in the quadrant $x < 0$, $y < 0$. The shape of the expansion on this side of the obstacle depends on the basal friction, and the depth and the speed of the flow at a given slope angle. For simplicity, the flow is assumed to be inviscid ($\nu =0$), and, hence, the steady conservation equations reduce to

(8.1)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} ( h {{\boldsymbol{\bar{u}}}}) = 0, \end{gather}$$

(8.2)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} ( h {{\boldsymbol{\bar{u}}}} \otimes {{\boldsymbol{\bar{u}}}} ) + \boldsymbol{\nabla}\left( \frac{1}{2} g h^2 \cos \zeta \right) = gh\left(\sin\zeta{\boldsymbol{{i}}}-\mu\cos \zeta \frac{{\boldsymbol{\bar{u}}}}{|{\boldsymbol{\bar{u}}}|}\right). \end{gather}$$

The downslope momentum balance equation (8.2) is further simplified by taking the basal friction coefficient to be a constant. This is chosen to be equal to the tangent of the slope angle ($\mu =\tan \zeta$) to reflect the fact that the flow adopts a near-uniform downstream speed $\bar u_0$ and thickness $h_0$ far downslope.

An asymptotic solution can be developed based on the assumption that the downslope component of velocity is greater than the lateral component, i.e. $| \bar v / \bar u | \equiv \epsilon \ll 1$. The governing equations (8.1)–(8.2) admit a set of distinguished scalings representing flow far downstream,

(8.3a–e)\begin{equation} \bar u=\bar u_0+O(\epsilon), \quad \bar v=O(\epsilon), \quad h=O(1), \quad x=O(1/\epsilon^{2}), \quad y=O(1/\epsilon), \end{equation}

under which the inertial and acceleration terms of the governing equations become negligible. At leading order, the downslope components of friction and gravity are in balance, while the cross-slope momentum equation reduces to a balance of cross-slope pressure gradients and basal friction, $h \partial h / \partial y= -\mu h \bar v/\bar u_0$. Using this cross-slope equation to write $\bar v$ in terms of $\partial h / \partial y$, the mass equation then becomes

(8.4)\begin{equation} \frac{\partial h}{\partial x} = \frac{1}{\mu} \frac{\partial}{\partial y}\left(h\frac{\partial h}{\partial y} \right). \end{equation}

As $y \to \infty$, the boundary condition $h \to h_0$ is applied, describing the flow approaching the thickness of the steady uniform flow far from the obstacle. The vacuum region is represented by the condition $h = 0$ for $y \leq y_b(x)$, where the unknown function $y_b(x)$ represents the boundary of the vacuum region, at which a kinematic condition $\bar v = \bar u_0\, \textrm {d}y_b/\textrm {d} x$ applies.

Transforming to similarity variables,

(8.5a,b)\begin{equation} h(x,y) = h_0 \mathcal{H}(\eta),\quad \eta = \sqrt{\frac{\mu}{h_0}} \frac{y}{\sqrt{x}}, \end{equation}

reduces the parabolic partial differential equation (8.4) to a second-order ODE,

(8.6)\begin{equation} \eta \mathcal{H}' + (\mathcal{H}^2)'' = 0, \end{equation}

where the prime indicates a derivative with respect to $\eta$. In these variables, the boundary condition of matching to the far-field flow becomes $\mathcal {H} \rightarrow 1$ as $\eta \rightarrow \infty$. The vacuum boundary at $y=y_b(x)$ becomes $\eta =\eta _0$, where $\eta _0$ is a constant to be determined and, at this point, $\mathcal {H}(\eta _0) = 0$ and the kinematic condition implies $\mathcal {H}'(\eta _0) = -\eta _0/2$. These three boundary conditions allow the second-order ODE (8.6) to be integrated numerically, and determine $\eta _0 \approx -1.2385$.

The vacuum boundary far downstream of the obstacle is then at

(8.7)\begin{equation} y = y_b(x) \equiv \eta_0 \sqrt{\frac{h_0}{\mu}} \sqrt{x}, \end{equation}

which starts at the downstream corner of the obstacle $y_b=0$ at $x=0$. Finding the intersection of this boundary with the obstacle centreline, which is at $y = -w / 2$ in the coordinate system used in this section, gives the length of the vacuum region,

(8.8)\begin{equation} x = x_v = \frac{\mu}{4\eta_0^2} \frac{w^2}{h_0}. \end{equation}

This solution is shown against the numerical results in figure 22 for a variety of obstacle widths and inflow conditions. For low Froude number values, the vacuum length is under-predicted by the theory with the assumption of constant Coulomb friction. This suggests that the levees at the edge of the vacuum region, which are of the greatest volume for low Froude number values, play a role in extending the vacuum region length. Empirically a better collapse can be shown by keeping the $w^2 / h$ dependence, but also including a ${ {Fr}}^{-1/4}$ dependence in place of the Coulomb friction coefficient $\mu$. This is shown in the inset plot in figure 22. In this plot it is shown that the vacuum length is well predicted for all cases. However, the physical explanation for the ${ {Fr}}^{-1/4}$ scaling is unclear.

Figure 22. The numerical vacuum region length against a candidate function for vacuum length from the theory. The inset plot shows the numerical vacuum length against an empirical scaling. The colours indicate the Froude number of the respective numerical results and marker symbols indicate various obstacle widths. The $w_{e}$ marker indicates the value from laboratory experiments.

9. Discussion and conclusions

This paper investigates the granular flow past a blunt obstacle on a rough inclined plane. Rough beds differ from smooth beds in that the friction allows steady uniform flows to develop over a wide range of inclination angles and Froude numbers. This enables the study of subcritical as well as supercritical flows past an obstacle, and it gives rise to a number of effects not seen on smooth beds.

For supercritical flows, a bow shock and static dead zone form upstream of the obstacle, and a grain-free vacuum region develops on the lee side, in a similar manner to supercritical flows on smooth beds (figures 1 and 6). However, as the vacuum forms, a static levee develops along its boundary during the flow. This is in complete contrast to smooth beds, where the largest velocities in the expansion are attained directly on the vacuum boundary. When the flow wanes, the flow gradually slows down and stops, leaving a layer of particles on the chute whose thickness is close to $h_{stop}$ (Pouliquen & Forterre 2002). In the final stages of the stopping process, the dead zone partially collapses and the material triggers a sequence of erosion-deposition waves that propagate downslope (Edwards & Gray 2015; Edwards et al. 2017; Viroulet et al. 2019; Edwards et al. 2021). All of these features are captured by a depth-averaged avalanche model (2.1)–(2.2) that uses a non-monotonic hysteretic friction law (2.5)–(2.8) (Edwards et al. 2019) and includes depth-averaged in-plane viscous terms (Gray & Edwards 2014; Baker et al. 2016a). The simulation results are shown in figures 10–14 as well as in the supplementary online movies.

Sustained subcritical flows do not usually develop on smooth beds. The rough bed experiments, shown in figures 1, 2 and 9, therefore allow a range of new phenomena to be observed for the first time. For subcritical flows, a bow shock does not form upstream of the obstacle, and instead gravitational waves propagate upstream, which convey the presence of the obstacle. This leads to a large upstream region of influence, where the oncoming flow is gradually deflected to either side of the obstacle. The experiments suggest that a St Andrew's cross type structure forms above the obstacle (see the schematic diagram figure 7b and experimental thickness contours figure 8b). In the lower quadrant (adjacent to the obstacle) a static dead zone forms, in the upper quadrant the flow decelerates in the downslope direction and spreads laterally, and in the two quadrants on either side the flow accelerates and flows parallel to the sides of the erodible dead zone. On the lee side, the flow expands again to form a grain-free vacuum that is bounded by static levees, which are more substantial than those formed during supercritical flow. All of these features can be visualised in the contrast enhanced movies that are available online. As the flow wanes, a $h_{stop}$ layer is deposited on most of the chute, but a raised ridge also forms upstream of the dead zone as shown in figure 2. This ridge develops as a consequence of the support that is provided by the obstacle, which extends far upstream when the chute inclination is close to the angle of repose. Numerical simulations of subcritical flow (figures 15–18) are able to capture all these features, and show evidence of the St Andrew's cross structure that exists upstream of the dead-zone apex.

In § 7 the Eikonal equation is used to develop a simplified model for the height of the maximal static pile that can be supported by the obstacle. In addition, by partially submerging this pile with a uniform thickness layer of flowing grains, and accounting for the support that this material provides (see figure 21), approximate solutions are obtained for the dead zone's scalloped shape (7.20) and its length (7.23). Figure 20(b) shows that this simplified model provides a good collapse of the numerically simulated dead-zone lengths for a wide range of obstacle widths, slope inclination angles and Froude numbers that span both the supercritical and subcritical regimes. Similarly, in § 8 an approximate similarity theory is developed, which provides a reasonable approximation for the length of the grain-free vacuum region in the lee of the obstacle.

It is anticipated that the results of this work may find application to a range of geophysical mass flows in the natural environment. These include snow avalanches, debris flows and dense volcanic pyroclastic avalanches, which pose a risk to people and infrastructure in mountainous and volcanic regions. In particular, in recent years global warming has led to the increased frequency of wet snow avalanches (Pielmeier et al. 2013; Naaim et al. 2016), which move at much slower speeds than their drier counterparts. Despite their slow speeds, wet snow avalanches impose very large loads on obstacles in their path and can be highly destructive. These flows readily form self-channelised flows with static levees in their runout zones (Gray & Kokelaar 2010; Bartelt et al. 2012), which suggests that a simple model might characterise their frictional behaviour with a rough-bed friction law of the form (2.5)–(2.8).