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. Reference Hutter, Siegel, Savage and Nohguchi1993; Pouliquen & Forterre Reference Pouliquen and Forterre2002). 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. Reference Pouliquen, Cassar, Jop, Forterre and Micolas2006) 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 Reference Gray and Edwards2014; Baker et al. Reference Baker, Barker and Gray2016a).

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. Reference Edwards, Russell, Johnson and Gray2019),

(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. (Reference Russell, Johnson, Edwards, Viroulet, Rocha and Gray2019) and Edwards et al. (Reference Edwards, Russell, Johnson and Gray2019), 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 (Reference Pouliquen and Forterre2002) and Edwards et al. (Reference Edwards, Russell, Johnson and Gray2019). 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. Reference Jop, Forterre and Pouliquen2006) by Gray & Edwards (Reference Gray and Edwards2014), 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. Reference Baker, Barker and Gray2016a; Rocha et al. Reference Rocha, Johnson and Gray2019). The viscous terms fundamentally change the system's structure from hyperbolic (Savage & Hutter Reference Savage and Hutter1989; Gray et al. Reference Gray, Tai and Noelle2003; Hákonardóttir & Hogg Reference Hákonardóttir and Hogg2005; Hogg et al. Reference Hogg, Gray and Cui2005; Cui et al. Reference Cui, Gray and Jóhannesson2007; Johnson & Gray Reference Johnson and Gray2011; Cui & Gray Reference Cui and Gray2013; Viroulet et al. Reference Viroulet, Baker, Edwards, Johnson, Gjaltema, Clavel and Gray2017; Edwards et al. Reference Edwards, Viroulet, Johnson and Gray2021) to hyperbolic-parabolic. The inclusion of the viscous terms is vital to properly model the formation of levees (Rocha et al. Reference Rocha, Johnson and Gray2019). 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. Reference Gray, Tai and Noelle2003; Hákonardóttir & Hogg Reference Hákonardóttir and Hogg2005; Gray & Cui Reference Gray and Cui2007; Pudasaini et al. Reference Pudasaini, Hutter, Hsiau, Tai, Wang and Katzenbach2007; Cui & Gray Reference Cui and Gray2013). 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 Reference Gray and Cui2007), 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. Reference Edwards, Viroulet, Kokelaar and Gray2017; Rocha et al. Reference Rocha, Johnson and Gray2019), 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. Reference Rocha, Johnson and Gray2019). 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 Reference Edwards and Gray2015; Edwards et al. Reference Edwards, Viroulet, Kokelaar and Gray2017; Rocha et al. Reference Rocha, Johnson and Gray2019; Viroulet et al. Reference Viroulet, Edwards, Johnson, Kokelaar and Gray2019; Edwards et al. Reference Edwards, Viroulet, Johnson and Gray2021). 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. (Reference Edwards, Russell, Johnson and Gray2019). 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 Reference Kurganov and Tadmor2000). 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 Reference Edwards and Gray2015; Baker et al. Reference Baker, Johnson and Gray2016b; Edwards et al. Reference Edwards, Viroulet, Kokelaar and Gray2017; Rocha et al. Reference Rocha, Johnson and Gray2019). 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 Reference LeVeque2002). 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 Reference Forterre and Pouliquen2003; Gray & Edwards Reference Gray and Edwards2014; Edwards & Gray Reference Edwards and Gray2015; Viroulet et al. Reference Viroulet, Baker, Rocha, Johnson, Kokelaar and Gray2018). 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. Reference Edwards, Russell, Johnson and Gray2019). 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. Reference Edwards, Viroulet, Kokelaar and Gray2017). 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. Reference Viroulet, Edwards, Johnson, Kokelaar and Gray2019; Edwards et al. Reference Edwards, Viroulet, Johnson and Gray2021). 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. Reference Edwards, Russell, Johnson and Gray2019) 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. Reference Viroulet, Edwards, Johnson, Kokelaar and Gray2019; Edwards et al. Reference Edwards, Viroulet, Johnson and Gray2021) 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 Reference Hadeler and Kuttler1999; Pauli & Gioia Reference Pauli and Gioia2007; Colombo, Guerra & Monti Reference Colombo, Guerra and Monti2012)

(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 Reference Sneddon1957; Stavroulakis & Tersian Reference Stavroulakis and Tersian2004). 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 Reference Stavroulakis and Tersian2004; Pauli & Gioia Reference Pauli and Gioia2007) 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 Reference Chadwick1976; Cui & Gray Reference Cui and Gray2013)

(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. (Reference Hogg, Gray and Cui2005). 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.