2.1. Mechanics

We focus on the limit in which the channel is sufficiently thin, or the piston sufficiently slow, that the inertia of the fluid–solid mixture is unimportant in the macroscopic momentum balance. In this limit, which we quantify more formally in § 2.3, conservation of linear momentum can be written as

(2.5) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma } = \rho g \boldsymbol{e}_z, \end{align}

where $\boldsymbol{\sigma }$ is the total stress in the mixture and $\rho =\rho _{\!f} \phi _{\!f} + \rho _s \phi _s$ is the local volume-averaged density of the mixture, and where we have adopted the convention that tensile stresses are positive. Going forward, we assume that the ratio of typical vertical to horizontal length scales is sufficiently small that we can adopt a lubrication-type approximation, in which $w \ll u$ and $\partial _x \ll \partial _z$ . We discuss the validity of these assumptions in §§ 5 and 6.

2.1.1. Upper layer (fluid)

The upper layer, $h\lt z\lt b$ , comprises a single incompressible Newtonian fluid of viscosity $\eta _{\!f}$ , for which the lubrication approximation suggests that the fluid pressure $p_{\!f}$ is hydrostatic, $p_{\!f}(x,z,t) = P_{\!f}(x,t) + \rho _{\!f} g (b - z)$ , such that

(2.6) \begin{align} \frac {\partial \sigma _{xz}}{\partial z} = \frac {\partial P_{\!f}}{\partial x}, \quad \textrm {with} \quad \sigma _{xz} = \eta _{\!f} \frac {\partial u_{\!f}}{\partial z}, \end{align}

where $P_{\!f}$ is the unknown fluid pressure at $z = b$ . We impose a no-slip boundary condition on the flow at the top wall, $u_{\!f} = 0$ on $z = b$ , as is standard. The appropriate boundary condition at the interface $z = h$ is less obvious, but for simplicity, we follow previous models of bedload transport and impose no tangential slip between the overlying fluid and the granular solid by requiring that $u_{\!f} = u_s$ at $z=h^+$ (e.g. Aussillous et al. Reference Aussillous, Chauchat, Pailha, Médale and Guazzelli2013). Note that this choice does not constrain the fluid velocity on the granular side of the interface, where in general we expect a non-zero Darcy flow (i.e. fluid flow relative to the solid grains), and neither does it prevent fluid flow across the interface. Note also that in reality, some degree of tangential slip is expected at the boundary between an external fluid flow and a porous medium (e.g. Le Bars & Worster Reference Le Bars and Worster2006). Many associated slip boundary conditions have been proposed, most famously the so-called Beavers–Joseph condition (Beavers & Joseph Reference Beavers and Joseph1967; Le Bars & Worster Reference Le Bars and Worster2006). We revisit this point briefly below, at the end of § 2.1.2.

Imposing these two no-slip conditions, we arrive at an expression for the velocity of the fluid in the upper layer:

(2.7) \begin{align} u_{\!f}(x, z, t) = \frac {(b - z)(h - z)}{2 \eta _{\!f}} \frac {\partial P_{\!f}}{\partial x} + \left (\frac {b - z}{b - h}\right )u_s(z = h). \end{align}

The shear stress exerted on the granular layer by this flow is then

(2.8) \begin{align} \sigma _{xz}(z = h) = \eta _{\!f} \left [ \frac {h - b}{2 \eta _{\!f}} \frac {\partial P_{\!f}}{\partial x} - \frac {u_{s}(z = h)}{b - h} \right ]\!. \end{align}

Equation (2.7) suggests that the flow in the upper layer is a Poiseuille–Couette flow, combining a parabolic profile from the background pressure gradient with a linear shear profile from the motion of the granular layer relative to the top wall.

2.1.2. Lower layer (fluid–solid mixture)

For the lower layer, we adopt Terzaghi’s principle (Terzaghi Reference Terzaghi1936) by decomposing the total stress into the fluid pressure and the effective solid stress $\boldsymbol{\sigma }'$ :

(2.9) \begin{align} \boldsymbol{\sigma } = \boldsymbol{\sigma }' - p_{\!f} \boldsymbol{I}. \end{align}

Following Jop et al. (Reference Jop, Forterre and Pouliquen2006), we then decompose $\boldsymbol{\sigma }'$ into its isotropic and deviatoric components:

(2.10) \begin{align} \boldsymbol{\sigma }' = -p_s \boldsymbol{I} + \boldsymbol{\tau }, \end{align}

where $p_s = -(1/3)\, \textrm {tr}(\boldsymbol{\sigma }')$ and $\textrm {tr}(\boldsymbol{\tau }) = 0$ . Further, we impose that the fluid flows through the porous, granular material according to Darcy’s law,

(2.11) \begin{align} (1 - \phi _s^*) (\boldsymbol{v}_{\!f} - \boldsymbol{v}_s) = -\frac {k}{\eta _{\!f}} \left ( \boldsymbol{\nabla }\! p_{\!f} - \rho _{\!f} g \boldsymbol{e}_z \right )\!, \end{align}

where $k$ is the permeability of the granular skeleton, which we expect to scale with the square of the grain size: $k\sim {}d^2$ .

Combining (2.5), (2.9) and (2.10) under the lubrication approximation, we arrive at

(2.12) \begin{align} \frac {\partial }{\partial z} (p_s + p_{\!f}) = -\rho g \quad \textrm {and} \quad \frac {\partial \sigma _{xz}'}{\partial z} = \frac {\partial }{\partial x} (p_s + p_{\!f}), \end{align}

along with $\sigma _{xx}' = \sigma _{zz}' = -p_s$ . The vertical component of (2.11) suggests that the fluid pressure is hydrostatic; further requiring that it be continuous across $z=h$ leads to $p_{\!f} = P_{\!f} + \rho _{\!f} g (b - z)$ . It then follows that $p_s = P_s(x, t) + \phi _s^*\, \Delta \rho\, g (h - z)$ , where $\Delta \rho = \rho _s - \rho _{\!f}$ , and $P_s$ is the as-yet-unspecified solid pressure at the interface $z = h$ . As such,

(2.13) \begin{align} \sigma _{xz}' = \sigma _{xz} &= \sigma _{xz}(z = h) - \left ( h-z \right ) \left ( \frac {\partial P}{\partial x} + \phi _s^*\, \Delta \rho\, g \frac {\partial h}{\partial x} \right )\!, \end{align}

where $P = P_{\!f} + P_s$ . Further, the horizontal component of (2.11) becomes

(2.14) \begin{align} (1 - \phi _s^*) (u_{\!f} - u_s) = - \frac {k}{\eta _{\!f}} \frac {\partial P_{\!f}}{\partial x}. \end{align}

Note that continuity of normal stress requires that $P_s$ must vanish when the solid does not bridge the gap (i.e. when $h \lt b$ ), such that $P_s \ne 0$ only in the bridged region (i.e. where $h \equiv b$ ). Correspondingly, conservation of mass requires that the fluid and solid in the bridged region move at the same speed, thus $\partial _x P_{\!f} = 0$ in the bridged region.

Imposing no tangential slip between the granular solid and the fluid in the upper layer, as we did above, implies a jump in tangential fluid velocity between the free flow in the upper layer and the Darcy flow in the lower layer (i.e. across the interface). In reality, we anticipate a narrow transition region at the interface (on the scale of one grain/pore, or $\sim {}d\sim \sqrt {k}$ ) over which the fluid velocity smoothly adjusts (see e.g. Beavers & Joseph Reference Beavers and Joseph1967; Le Bars & Worster Reference Le Bars and Worster2006). We expect the difference between applying a no-slip condition and a more elaborate slip condition to scale with $\sqrt {k}$ , which is typically extremely small. As noted above, this boundary condition does not constrain flow across the interface; rather, the normal component of the fluid velocity is dictated by conservation of mass, as is standard in lubrication problems.

Note that the incomplete model above already provides some insight into the flow in the far field, where the grains are undisturbed at uniform depth $h_0$ , and the fluid flows both above them and through them. For stationary grains ( $u_s=0$ ), the flow in the upper layer is a simple Poiseuille flow (see (2.7)), and the uniform Darcy flow through the lower layer follows from (2.14). The total flow rate is $\dot {\delta }b$ , as determined by the piston, with a fraction $(b-h_0)^3/[(b-h_0)^3+12 h_0 k]$ flowing in the upper layer, and a fraction $12 h_0 k/[(b-h_0)^3+12 h_0 k]$ flowing through the grains in the lower layer.

2.1.3. Granular friction and wall friction

As noted above, we adopt a $\mu (I)$ -type rheology for the internal mechanics of the granular material (Cassar et al. Reference Cassar, Nicolas and Pouliquen2005; Jop et al. Reference Jop, Forterre and Pouliquen2005; Pouliquen et al. Reference Pouliquen, Cassar, Jop, Forterre and Nicolas2006). In these model constructions, the internal mechanical interactions between grains are taken to obey a Coulomb-like friction model:

(2.15) \begin{align} |\sigma _{xz}'| = \mu p_s \quad \textrm {if } \left | \sigma _{xz}' \right | \geqslant \mu _s p_s, \end{align}

(2.16) \begin{align} \frac {\partial u_s}{\partial z} = 0 \quad\textrm {otherwise}, \end{align}

where $\mu$ is a dimensionless, rate-dependent coefficient of kinetic friction, as discussed further in § 2.2, and $\mu _s$ is a dimensionless constant coefficient of static friction. If the magnitude of the shear stress is sufficiently small, then (2.16) indicates that grains will resist internal rearrangement and the solid phase will act locally as a rigid body. If the shear stress is sufficiently large, then the solid phase will yield such that $|\partial _z u_s| \gt 0$ and the solid velocity will satisfy (2.15), with the sign of $\sigma _{xz}'$ given by the sign of the shear rate $\partial _z u_s$ . As a result, the vertical profile $u_s(z)$ may feature multiple coupled yielded and unyielded regions separated by free boundaries. We consider the various possible flow scenarios in § 2.1.4. We specify the granular rheology in more detail in § 2.2.

At the bottom surface of the granular layer, we assume that the mechanical interactions between the grains and the wall obey a Coulomb-like friction model:

(2.17) \begin{align} \sigma _{xz}'(z = 0) \leqslant \mu _w p_s(z = 0), \end{align}

such that

(2.18) \begin{align} u_s(z = 0) = 0 \quad \textrm {if } \quad \sigma _{xz}'(z = 0) \lt \mu _w p_s(z = 0), \end{align}

(2.19) \begin{align} u_s(z = 0) \gt 0 \quad \textrm {if } \quad \sigma _{xz}'(z = 0) = \mu _w p_s(z = 0), \end{align}

where $\mu _w$ is the dimensionless constant coefficient of friction between the grains and the wall. Note that we assume that $u_{s} \geqslant 0$ , such that $\sigma _{xz}'(z=0)\geqslant 0$ , and thus that friction resists relative motion. To generate a single, continuous problem in $x$ , we unify these conditions in a regularised form in § 2.3.

The appropriate boundary condition at the top surface of the granular layer depends on whether or not the grains are in contact with the top wall. If the grains have bridged the gap (i.e. if $h \equiv b$ ), then we similarly apply Coulomb friction at the top wall:

(2.20) \begin{align} \sigma _{xz}'(z=b) \geqslant -\mu _w p_s(z=b) \quad \mathrm{if}\ h\equiv b, \end{align}

such that

(2.21) \begin{align} u_s(z=h) = 0 \quad \textrm {if } \quad \sigma _{xz}'(z=b) \gt -\mu _w p_s(z=b) \quad \mathrm{and} \quad h\equiv b, \end{align}

(2.22) \begin{align} u_s(z=h) \gt 0 \quad \textrm {if } \quad \sigma _{xz}'(z=b) = -\mu _w p_s(z=b) \quad \mathrm{and} \quad h\equiv b, \end{align}

where the differences in sign relative to (2.17)–(2.19) reflect the different direction of the normal vector, again ensuring that friction resists relative motion for $u_s\geqslant 0$ . If the grains have not bridged the gap between the plates (i.e. if $h\lt b$ ), then we instead impose that the shear stress is continuous with the upper fluid layer via (2.8). To generate a single, continuous problem in $x$ , we also unify these conditions in a regularised form in § 2.3.

2.1.4. Velocity profile

Before giving specific details of the granular rheology, we first present in figures 2 and 3 a selection of example velocity profiles that can occur in the transition zone ( $h\lt b$ ) and in the bridged region ( $h\equiv {}b)$ , respectively. In both cases, the qualitative nature of the profile depends centrally on whether the coefficient of wall friction, $\mu _w$ , is larger or smaller than the static coefficient of internal friction, $\mu _s$ .

Figure 2.

Examples of the different possible stress and velocity configurations in the unbridged state ( $h \lt b$ and $P_s\equiv {}0$ ), neglecting the fluid phase ( $\partial _x{p_{\!f}}=0$ and $\sigma '_{xz}(z=h) = 0$ ). (a) Possible shear-stress profiles $\sigma _{xz}'(z)$ (red and blue lines), which must be linear from (2.13), together with the effective yield stress of the granular layer $\pm \mu _s p_s$ , where $p_s = \phi _s^* g\, \Delta \rho\, (h - z)$ (solid yellow); the shear stress lies below the yield stress in the shaded region. Different colours correspond to whether the wall friction is less than the internal static friction, $\mu _w = \mu _{w1} \lt \mu _s$ (blue), or greater, $\mu _w = \mu _{w2} \gt \mu _s$ (red), with corresponding velocity profiles shown in (b,c) and (d,e), respectively. In each case, there are two possibilities: (b,d) the shear stress on the bottom wall does not reach the slip stress, $|\sigma _{xz}'(z=0)|\lt \mu _w p_s(0)$ , and the grains do not slip there, or (c,e) it does and they do. Where the shear stress is below the internal yield stress, $|\sigma _{xz}'|\lt \mu _sp_s$ , the granular layer behaves as a rigid plug, as in (b,c).

Figure 3.

Examples of the different possible stress and velocity configurations in the bridged state ( $h \equiv b$ , $P_s\gt 0$ and $\sigma '_{xz}(z=h) \lt 0$ ), neglecting the fluid phase ( $\partial _x{p_{\!f}}=0$ ). (a) Possible shear-stress profiles $\sigma _{xz}'(z)$ (red and blue lines), together with the effective yield stress of the granular layer $\pm \mu _s p_s$ , where $p_s = P_s + \phi _s^* g\, \Delta \rho\, (b-z)$ (solid yellow); the shear stress lies below the yield stress in the shaded region. If $\mu _w = \mu _{w1} \lt \mu _s$ (blue), then the granular layer does not yield anywhere and must move as a rigid plug, as in (b). If $\mu _w = \mu _{w1} \gt \mu _s$ (red), then different flow profiles are possible depending on whether the shear stress attains the slip stress on either the upper or lower wall, as sketched in (c,d,e), each corresponding to the equivalent line style in (a).

We first consider the unbridged portions of the granular layer (i.e. the transition region and the far field) in figure 2, for which $h \lt b$ and $P_s \equiv 0$ . For simplicity, we neglect the role of the fluid by taking $\sigma _{xz}'(z = h) = 0$ and $\partial _xp_{\!f}=0$ . Recalling that the shear stress at any horizontal location increases linearly with depth (2.13), we find that four types of qualitative profile arise. When $\mu _w \lt \mu _s$ (blue curves), the shear stress can never grow large enough for the granular layer to yield internally: $|\sigma _{xz}'| \lt \mu _s p_s$ and $|\partial _z u_s| \equiv 0$ for all $z$ . As a result, the grains do not deform, and either remain stationary (figure 2 b) or move as a rigid body (figure 2 c), depending on whether $\sigma _{xz}'(z = 0)$ is less than or equal to $\mu _w p_s(z=0)$ , respectively. By contrast, when $\mu _w \gt \mu _s$ (red curves), the granular layer yields over its entire thickness, again with either no slip (figure 2 d) or slip (figure 2 e) on the base, depending on the size of the basal stress. Note that analogous profiles have been observed in the case of viscoplastic fluids with wall slip (Lawal et al. Reference Lawal, Kalyon and Yilmazer1993).

We next consider the bridged portions of the granular layer, in figure 3, for which $h\equiv 0$ , $P_s\gt 0$ , and $\sigma _{xz}'(z=h= b)$ is dictated by friction with the top wall. Note that we again neglect the role of the fluid. As in the unbridged portions (figure 2), there are different cases depending on whether $\mu _w \lt \mu _s$ or $\mu _w \gt \mu _s$ . In the former case (blue, figure 3 b), the grains remain everywhere unyielded and slide over both boundaries as a rigid plug (because of the moving piston, the bridged region cannot remain stationary). If $\mu _w \gt \mu _s$ (red, figures 3 c–e), then there are more possibilities. The grains may be able to slip, or not, on both the upper and lower walls, depending on the local shear and normal stresses, and in all cases there may be some unyielded interior region where the shear stress falls below the effective yield stress (figure 3 a). Across the entire parameter space, we expect the solid velocity to be smaller near the bottom wall than near the top wall because $p_s$ increases from top to bottom due to gravity.

To accommodate the scenarios shown in figures 2 and 3, we seek in general a piecewise solid velocity of the form

(2.23) \begin{align} u_s=\begin{cases} u_{s1}(x,z,t) & \text{if} \ 0\leqslant z \leqslant h_{y1},\\[3pt] u_{s2}(x,t) & \text{if} \ h_{y1} \lt z \leqslant h_{y2},\\[3pt] u_{s3}(x,z,t) & \text{if} \ h_{y2} \lt z \leqslant h,\\ \end{cases} \end{align}

such that the granular layer is unyielded ( $\partial _z u_s = 0$ ) in an interior region $h_{y1} \lt z \lt h_{y2}$ and yielded ( $|\partial _z u_{s}| \gt 0$ ) outside that region, with continuity of velocity requiring that $u_{s1} = u_{s2}$ at $z = h_{y1}$ and $u_{s2} = u_{s3}$ at $z = h_{y2}$ . The locations $h_{y1}$ and $h_{y2}$ are determined by solving $\sigma _{xz}' = \pm \mu _s p_s$ to give

(2.24) \begin{align} h_{y1} = h - \frac {\sigma _{xz}'(z = h) - \mu _s P_s}{\partial _x P + \phi _s^*\, \Delta \rho\, g (\partial _x h + \mu _s) } \quad \textrm {and} \quad h_{y2} = h - \frac {\sigma _{xz}'(z = h) + \mu _s P_s}{\partial _x P + \phi _s^*\, \Delta \rho\, g (\partial _x h - \mu _s) }. \end{align}

In the yielded regions ( $0 \lt z \lt h_{y1}$ and $h_{y2} \lt z \lt h$ ), we combine (2.13) and (2.15) with the appropriate granular rheology from § 2.2 to obtain ordinary differential equations (ODEs) for $u_{s1}$ and $u_{s3}$ .

2.2. Granular rheology

In the region where the granular layer has yielded, $|\sigma _{xz}'| \gt \mu _s p_s$ , we adopt a $\mu (I)$ -type rheology to describe its evolution (Cassar et al. Reference Cassar, Nicolas and Pouliquen2005; Jop et al. Reference Jop, Forterre and Pouliquen2005; Pouliquen et al. Reference Pouliquen, Cassar, Jop, Forterre and Nicolas2006). These rheologies relate the effective viscosity of the solid phase, $\mu$ , to either the inertial number $I$ , or the viscous number $I_v$ , both of which are defined below. The specific choice depends on the magnitude of the Stokes number $\textit{St} \sim I^2/I_v \sim \rho _s d^2\, |\partial _z u_{s}|/ \eta _{\!f}$ , which compares inertial forces to viscous forces on the scale of the grains. Viscous drag from the fluid is expected to dominate inertia at the grain scale for $\textit{St}\ll {}1$ , whereas grain-scale inertia is expected to be dominant for $\textit{St}\gg {}1$ . As such, different rheologies are appropriate in these two regimes. Some effort has been made to combine the two regimes into a single rheology (Boyer et al. Reference Boyer, Guazzelli and Pouliquen2011; Baumgarten & Kamrin Reference Baumgarten and Kamrin2019; Rettinger et al. Reference Rettinger, Eibl, Rüde and Vowinckel2022; Tapia et al. Reference Tapia, Ichihara, Pouliquen and Guazzelli2022), but we treat them separately for simplicity.

In § 3, we consider the ‘dry’ case where $\textit{St}\gg {}1$ , in which limit the fluid phase is negligible. In this case, we follow Iordanoff & Khonsari (Reference Iordanoff and Khonsari2004) and Da Cruz et al. (Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005) by taking

(2.25) \begin{align} \mu (I) = \mu _s + \left (\frac {\mu _2 - \mu _s}{I + I_0}\right ) I, \end{align}

where $I_0$ and $\mu _2$ are constants. The inertial number $I$ is

(2.26) \begin{align} I = \sqrt {\frac {d^2 \rho _s}{p_s}} \left | \frac {\partial u_s}{\partial z} \right |\!, \end{align}

where, under the lubrication approximation, $|\partial u_s/\partial z|$ is the leading-order contribution to the magnitude of the rate-of-strain tensor. Equation (2.25) is an empirical fit to experiments and numerical simulations, and is defined such that $\mu \to \mu _s$ for $I\ll I_0$ and $\mu \to \mu _2$ for $I\gg {}I_0$ .

In § 4, we consider the ‘wet’ case where $\textit{St}\ll {}1$ , in which limit fluid viscosity dominates granular inertia. In this case, we follow Boyer et al. (Reference Boyer, Guazzelli and Pouliquen2011) by taking

(2.27) \begin{align} \mu (I_v) = \mu _s + \left (\frac {\mu _2 - \mu _s}{I_v + I_0} \right )I_v + I_v + \frac {5}{2}\phi _s^*\sqrt {I_v}, \end{align}

where $I_0$ and $\mu _2$ are again constants, the appropriate values of which are not necessarily the same as in (2.25). The viscous number $I_v$ is

(2.28) \begin{align} I_v = \frac {\eta _{\!f}}{p_s} \left | \frac {\partial u_s}{\partial z} \right | \!, \end{align}

having again taken $|\partial _z u_s|$ to be the leading-order contribution to the magnitude of the rate-of-strain tensor. The last two terms on the right-hand side of (2.27) mean that $\mu (I_v)$ will increase without bound as $I_v$ increases.

By substituting either (2.25) and (2.26) or (2.27) and (2.28) into the frictional-rheology model in (2.15), we arrive at a nonlinear ODE for $u_s$ , with boundary conditions described in more detail below.

Note that both the $\mu (I)$ and $\mu (I_\nu )$ models presented here have been shown to be mathematically ill-posed in their general, time-dependent, tensorial forms (Barker et al. Reference Barker, Schaeffer, Bohorquez and Gray2015, Reference Barker, Gray, Schaeffer and Shearer2023); specifically, they can be unstable to arbitrarily small-wavelength perturbations for some ranges of $I$ and $I_\nu$ . This problem can be resolved by carefully accounting for compressibility in each case (Barker et al. Reference Barker, Schaeffer, Shearer and Gray2017, Reference Barker, Gray, Schaeffer and Shearer2023). In the present lubrication framework, the momentum balance (2.12) is quasi-steady, which avoids any direct issues of ill-posedness. In the future, it could be instructive to consider a compressible version of the present model that remains formally well-posed even away from the lubrication limit.

2.3. Non-dimensionalisation and regularisation

To non-dimensionalise our model, we use $b$ , $b/\dot {\delta }$ , $\dot {\delta }$ and $\phi _s^*\, \Delta \rho\, g b$ as characteristic spatial, temporal, velocity and pressure scales. The key non-dimensional parameters are then

(2.29) \begin{align} \theta = \frac {\eta _{\!f} \dot {\delta }}{b^2 \phi _s^*\, \Delta \rho\, g}, \quad K = \frac {k}{b^2} \quad \textrm {and} \quad \textit{Fr} = \frac {\dot {\delta } d}{b \sqrt {b g}}, \end{align}

as well as $\mu _w$ , $\mu _s$ , and the other constants in the $\mu (I)$ and $\mu (I_v)$ constitutive laws. The parameter $K$ is the non-dimensional permeability. Noting that $\eta _{\!f} \dot {\delta } / b$ is the characteristic shear stress from the upper fluid layer, the parameter $\theta$ can be viewed as a modified Shields parameter that measures the strength of fluid shear relative to the weight of the granular bed. The modified Froude number $\textit{Fr}$ measures the strength of inertial forces relative to gravitational forces, scaled by $d/b$ . The ratio $\textit{Fr}^2/\theta \equiv \textit{St}$ recovers the Stokes number, which, as noted above, compares the relative importance of inertial and viscous forces on the scale of the grains. In the ‘wet’ granular limit ( $\textit{St} \ll 1$ ), we find that $\textit{Fr}$ drops out of the problem; in the ‘dry’ granular limit ( $\textit{St} \gg 1$ ), $\textit{Fr}$ becomes the main control parameter.

In § 2.1, we neglected the inertia of the fluid–solid mixture in the macroscopic momentum balance (2.5). The present scalings suggest that the relative importance of macroscopic inertia is governed by the square of the macroscopic Froude number, $\mathcal{F}^2 \equiv \dot \delta ^2/(bg) = (\textit{Fr} \,b/d)^2 = \textit{St} \, \theta (b/d)^2$ . Since the grains are much smaller than the gap, $\mathcal{F}$ is always larger than $\textit{Fr}$ , meaning that neglecting macroscopic inertia is reasonable only when $\textit{Fr}^2 \ll 1$ . Nonetheless, grain-scale inertia will control the granular rheology when $\textit{St} = \textit{Fr}^2/\theta \gg 1$ (i.e. in the ‘dry’ limit). Together, these criteria suggest that macroscopic inertia can be neglected, even when the local dynamics are controlled by grain-scale inertia, when $\theta \ll \textit{Fr}^2 \ll 1$ . In the ‘wet’ limit (i.e. $\textit{St}\ll {}1$ ), $\textit{Fr}$ drops out and can be freely chosen to ensure that inertial effects are negligible. Note that we have also neglected the inertia of the overlying fluid layer, which is reasonable when the fluid Reynolds number $\textit{Re}_{\!f} = \rho _{\!f} \dot \delta b/\eta _{\!f}$ is sufficiently small.

From this point on, we retain our original notation, but all quantities are dimensionless unless explicitly noted otherwise. The evolution equation for the height of the granular layer, (2.2), becomes

(2.30) \begin{align} \frac {\partial h}{\partial t} + \frac {\partial }{\partial x} (\bar {u}_s h) = 0,\quad \textrm {with} \quad \bar {u}_s h = \int _0^{h_{y1}} u_{s1} \, \text{d}z + (h_{y2} - h_{y1}) u_{s2} + \int _{h_{y2}}^h u_{s3} \, \text{d}z, \end{align}

and

(2.31) \begin{align} h_{y1} = h - \frac {\sigma _{xz}'(z = h) - \mu _s P_s}{\partial _x (P + h) + \mu _s } \quad \textrm {and} \quad h_{y2} = h - \frac {\sigma _{xz}'(z = h) + \mu _s P_s}{\partial _x (P + h) - \mu _s }, \end{align}

recalling that $P = P_{\!f} + P_s$ , where $P_{\!f}$ and $P_s$ are the dimensionless fluid pressure at $z=1$ and the dimensionless solid pressure at $z=h$ , respectively. To solve for the velocity profile $u_s(z)$ , we equate (2.13) for the shear stress with (2.15) for the rheology to obtain

(2.32) \begin{align} \sigma _{xz}'(z=h) - (h - z) \left ( \frac {\partial P}{\partial x} + \frac {\partial h}{\partial x} \right ) &= \mu p_s, \quad 0 \leqslant z \lt h_{y1},\\[-12pt]\nonumber \end{align}

(2.33) \begin{align} \frac {\partial u_s}{\partial z} &= 0, \quad h_{y1} \leqslant z \lt h_{y2}, \\[-12pt]\nonumber \end{align}

(2.34) \begin{align} \sigma _{xz}'(z=h) - (h - z) \left ( \frac {\partial P}{\partial x} + \frac {\partial h}{\partial x} \right ) &= -\mu p_s, \quad h_{y2} \leqslant z \lt h, \\[10pt] \nonumber \end{align}

where $p_s = P_s + (h-z)$ is the dimensionless solid pressure, and noting that the velocity profile $u_s(z)$ must be continuous in $z$ . The shear stress at the top of the granular layer, $\sigma _{xz}'(z=h)$ , is given below. As discussed in § 2.2, we use a frictional rheology model with $\mu = \mu (I)$ for $\textit{St} \gg 1$ (2.25) and $\mu = \mu (I_v)$ for $\textit{St} \ll 1$ (2.27), where the inertial and viscous numbers are

(2.35) \begin{align} I = \frac {\textit{Fr}}{\sqrt {p_s}} \left |\frac {\partial u_s}{\partial z} \right | \quad \textrm {and} \quad I_v = \frac {\theta }{p_s} \left |\frac {\partial u_{s}}{\partial z}\right |\!, \end{align}

respectively, and the associated expressions for $\mu$ are given in (2.25) and (2.27).

The shear stress $\sigma _{xz}^\prime (z=0)$ is determined by wall friction via (2.18) and (2.19). Similarly to internal granular friction, a challenge with wall friction is that the horizontal profile $u_s(x,z=0)$ may feature multiple coupled sliding and non-sliding regions. To avoid the need to track these regions explicitly, we unify (2.18) and (2.19) in a dimensionless, regularised form as

(2.36) \begin{align} \sigma _{xz}'(z=0) = \mu _w(P_s+h)\tanh \left [ \frac {u_s(z=0)}{\varepsilon } \right ]\!, \end{align}

where $\varepsilon \ll 1$ is a dimensionless regularisation parameter. Equation (2.36) provides a single continuous expression for $\sigma _{xz}'(z = 0)$ over the entire $x$ -domain, ensuring that $|\sigma _{xz}'(z = 0)| \lt \mu _w(P_s+h)$ when $u_s(z = 0) \lesssim {} \varepsilon \ll 1$ , while converging rapidly to (2.19) for larger values of $u_s(z=0)$ . This type of regularisation is common in the numerical simulation of viscoplastic fluids, because it avoids the need to introduce internal free boundaries between yielded and unyielded regions (e.g. Balmforth et al. Reference Balmforth, Frigaard and Ovarlez2014). Adopting it here for wall friction allows us to solve a single, continuous problem in $x$ , which is convenient for the quasi-one-dimensional flow studied here and even more-so for complex, quasi-two-dimensional frictional-flow patterns (i.e. Sandnes et al. Reference Sandnes, Flekkøy, Knudsen, Måløy and See2011). Some care should be taken with regularisations of this form to avoid the potential misidentification of yielded regions (e.g. Putz et al. Reference Putz, Frigaard and Martinez2009); we illustrate the convergence of our solutions with decreasing $\varepsilon$ in Appendix A.

The shear stress $\sigma _{xz}^\prime (z=h)$ is determined locally by wall friction if the grains are in contact with the top wall (i.e. (2.21) and (2.22)) for $h\equiv {}1$ ), or by fluid shear otherwise (i.e. (2.8) for $h\lt 1$ ). These scenarios can be captured concisely by introducing an indicator function $\chi (h)$ that transitions from $\chi (h\lt 1)=0$ to $\chi (h\equiv 1)=1$ , such that

(2.37) \begin{align} \sigma _{xz}'(z = h) = -\chi \mu _w P_s \tanh \left [ \frac {u_s(z = h)}{\varepsilon } \right ] + (1 - \chi ) \left [ \frac {h - 1}{2} \frac {\partial P_{\!f}}{\partial x} - \frac {\theta u_{s}(z = h)}{1 - {h}} \right ]\!, \end{align}

where the first term on the right-hand side is regularised wall friction and the second is viscous shear.

In (2.37), the term proportional to $\theta$ is inversely proportional to $1-h$ , which illustrates an awkward feature of the model. As discussed in § 2.1.1, the upper fluid layer satisfies no-slip conditions on both the stationary top wall and the moving granular layer. However, the thickness of this fluid layer vanishes as the granular layer approaches and then touches the top wall (i.e. as $h\to {}1$ , the bridging transition). As a result, the shear stress required to drive relative motion between the granular layer and the wall – which requires shearing the vanishing fluid layer – can become arbitrarily large as $h\to 1$ , driving the system towards a no-slip condition between the granular layer and the wall. For a system of two fluid layers, this smooth approach toward no-slip between the lower layer and the wall is physically appropriate. In the present scenario, however, the lower layer should instead transition a frictional condition when it touches the upper wall; this frictional condition features an upper limit on the shear stress, allowing for sliding (i.e. slip) when that limit is reached. Hence the lower layer could experience a jump in velocity and shear stress upon contacting the top wall unless this transition is regularised. In reality, the top of the granular layer is a permeable surface with grain-scale roughness. As a result, we expect fluid flow in the upper layer to transition smoothly into flow through the granular porous medium, and the shear stress at the top surface of the granular layer to transition smoothly from fluid shear to wall friction, as the irregular granular surface gradually makes contact with the top wall.

To reflect this smooth transition, we regularise the implementation of (2.37) in two ways. First, we place a non-zero minimum on the thickness of the fluid layer (i.e. on $1-h$ ) as it appears in the final term of (2.37), which reflects the fact that this construction breaks once the gap becomes comparable to the grain size. Specifically, we use the dimensionless permeability $K$ as a floor on $1-h$ in the last term in (2.37) in order to link the regularisation to the grain size ( $\sqrt {K}$ would scale more directly with the grain diameter, but would also regularise more aggressively). Second, we replace the sharp indicator function $\chi (h)$ with a smoothed step function of the form $\chi = h^n$ with $n\gg 1$ , thus regularising the bridging transition. The latter regularisation has the added benefit of preserving a continuous problem in $x$ , avoiding the need to track the free boundary between bridged and unbridged regions. We show in Appendix A, and specifically figure 17(c), that our results are insensitive to the value of $n$ provided that it is large enough.

Equations (2.32)–(2.34) lead to a nonlinear first-order ODE for $u_s(z)$ that also depends on $P_s$ and $\partial _x P_{\!f}$ . Expressions for the latter two quantities follow from evaluating (2.13) at $z=0$ :

(2.38) \begin{equation} \sigma _{xz}'(z = h) - \sigma _{xz}'(z = 0) = h \left ( \frac {\partial P}{\partial x} + \frac {\partial h}{\partial x} \right )\!. \end{equation}

Combined with (2.36) and (2.37), this relationship determines $\partial _x P_{\!f}$ in the unbridged region (where $P_s = 0$ ), and $P_s$ in the bridged region (where $\partial _x P_{\!f} = 0$ ), recalling that $P = P_s + P_{\!f}$ .

Finally, the ODE for $u_s$ requires a boundary condition. The appropriate constraint follows from global conservation of mass, (2.4), which becomes

(2.39) \begin{align} 1 = h \bar {u}_s - \frac {\bar {k}}{\theta } \frac {\partial P_{\!f}}{\partial x} + \frac {1 - h}{2} u_s(z = h), \end{align}

where $\bar {k} = K h + (1 - h)^3/12$ is the dimensionless gap-averaged permeability.

2.4. Numerical details – boundary and initial conditions

To solve the model presented in § 2.3, we simplify computation by shifting to a frame $X=x-t$ that moves with the piston. With this change of frame, (2.30) becomes

(2.40) \begin{align} \frac {\partial h}{\partial t} + \frac {\partial }{\partial X} \left [ (\bar {u}_s - 1) h \right ] = 0. \end{align}

We solve the problem from initial condition $h(X,0)=h_0$ and subject to boundary conditions $u_s = u_{\!f} = 0$ at $X=0$ (grains and fluid travelling with the impermeable piston) and $\partial h/\partial X = 0$ at $X=L_X$ (grains undisturbed in the far field), where the domain size $L_X$ is chosen to be sufficiently far ahead of the transition zone that its value has no impact on the results. We discretise the numerical domain $0 \leqslant X \leqslant L_X$ into $1+L_X/\Delta {X}$ nodes of equal spacing $\Delta X = 5 \times 10^{-3}$ . We then solve (2.40) for $h$ by approximating the spatial derivative using a second-order, essentially non-oscillatory scheme and integrating in time using the forward Euler method with fixed time step $\Delta t = 10^{-5}$ .

Evolving (2.40) in time requires that we compute $\bar {u}_s$ at every discrete point $X_j$ via (2.31)–(2.39) at each time step. To do so, we discretise the vertical domain at each point $X_j$ , $0 \leqslant z \leqslant h(X_j)$ , into $N_z=250$ points of equal spacing $\Delta {z}_j=h(X_j)/(N_z-1)$ . We then approximate $\partial _z u_s$ and $\partial _X P_s$ at each point $X_j$ as $(\partial _z u_s)_j\approx (u_s(X_j,z_{i+1})-u_s(X_j,z_i))/\Delta z_j$ and $(\partial _X P_s)_j \approx (P_s(X_{j+1})-P_s(X_{j}))/\Delta X$ , respectively, noting that $P_s=0$ in the unbridged region. Thus (2.31)–(2.39) become a coupled, nonlinear system of equations in $N_z + 6$ unknowns: $P_s(X_{j})$ , $(\partial _X P_{\!f})_j$ , $\sigma _{xz}'(X_j,0)$ , $\sigma _{xz}'(X_j,h(X_j))$ , $h_{y1}(X_j)$ , $h_{y2}(X_j)$ and $u_s(X_j,z_i)$ for $i = 1$ to $N_z$ . We solve this system of equations via Newton’s method using the solution at the previous time step as the initial guess. We then compute $\bar {u}_s(X_j)$ by averaging $u_s$ over appropriate regions in $z$ .

For the regularisation parameters in (2.36) and (2.37), we fix $\varepsilon = 10^{-3}$ and $n = 200$ ; our results are insensitive to these specific values (see Appendix A).

3.1. Bridged region

Since $h\equiv 1$ in the bridged region by definition, (2.30) becomes unnecessary and (2.39) simplifies to $\bar {u}_s = 1$ . However, as discussed in § 2.1.2, the solid velocity is non-trivial due to intra-granular friction and wall friction. As illustrated in the first row of figure 5, when $\textit{Fr}$ is sufficiently small, the grains will yield over the entire gap around the start of the bridged region ( $h_{y1} = h_{y2} = 1$ ), and the grains do not slide over the bottom plate. As the length of the bridged region grows, the grains stop sliding over the top plate, resulting in the velocity profile illustrated in figure 3(e). When $\mu _w \lt \mu _2$ (figure 5 a), as $\textit{Fr}$ increases, we find that the grains begin to slide over both the top and bottom plates. Further, for large $\textit{Fr}$ , the grains appear to resist internal rearrangement such that $|\partial _z u_s| \ll 1$ across the entire gap. In comparison, when $\mu _w \gt \mu _2$ (figure 5 b), we find that the grains do not slide over the bottom or top plate for all $\textit{Fr}$ , and the velocity profile does not appear to change significantly as $\textit{Fr}$ is varied.

To better understand the influence of the Froude number and the wall friction parameter, we compute the velocity profile of the grains at $z = 1$ and $0$ over a range of values of $\textit{Fr}$ with $\mu _s \lt \mu _w = 0.55 \lt \mu _2$ , as shown in figures 6(a) and 6(b). The velocity of the grains in contact with the top plate (figure 6 a) is fastest at the start of the bridged region, and monotonically decreases to $0$ as the length of the bridged region grows, after which the grains are stationary. Further, the length of the region where grains slide over the top plate increases with $\textit{Fr}$ . For the grains in contact with the bottom plate (figure 6 b), $\textit{Fr}$ must be sufficiently large or the grains will be stationary over the entire bridged region. When $\textit{Fr}$ becomes sufficiently large, the grains will begin to slide, which corresponds to the point where $\mu (I) = \mu _w$ , as shown in figure 6(c). However, unlike grains in contact with the top plate, grains will be stationary at the start of the bridged zone, resulting in a region where grains slide sandwiched between regions where the grains are stationary. As $\textit{Fr}$ increases further, this region where $u_s(z = 0) \gt 0$ will ‘spread’ and the grains at the start of the bridged zone will begin to move. In comparison, when $\mu _w = 0.85$ , grains do not slide over either plate for all $\textit{Fr}$ . Figure 6(c) indicates that while $\mu (I)$ (denoted by dashed lines for $\mu _w = 0.85$ ) increases with $\textit{Fr}$ as it does for $\mu _w = 0.55$ , it will always be smaller than $\mu _w$ , in this case, because it is bounded above by $\mu _2 = 0.7 \lt \mu _w$ .

Figure 6.

Velocity of the grains on the (a) top and (b) bottom plates with $\mu _w = 0.55$ and (light to dark blue) $\textit{Fr} = 0.01$ , 0.04, 0.06, 0.1 and 0.15. (c) The friction factor $\mu (I)$ evaluated at $z = 0$ for $\mu _w = 0.55$ (solid) and $\mu_w=0.85$ (dashed). The dashed black lines denote $\mu _s = 0.4$ , and the two values of $\mu _w$ being considered.

Figure 6. Long description

Panel A: A line graph shows the velocity of grains on the top plate. The x-axis is labeled with x-tilde ranging from -5 to 2, and the y-axis is labeled with u sub s(z equals 1) ranging from 0 to 1.5. The graph includes multiple lines in varying shades of blue representing different values. An arrow labeled Fr is also present. Panel B: Another line graph displays the velocity of grains on the bottom plate. The x-axis is labeled with x-tilde ranging from -5 to 2, and the y-axis is labeled with u sub s(z equals 0) ranging from 0 to 1.5. This graph also includes multiple lines in varying shades of blue. Panel C: A line graph illustrates the friction factor evaluated at z equals 0. The x-axis is labeled with x-tilde ranging from -5 to 2, and the y-axis is labeled with mu(I(z equals 0)) ranging from 0.4 to 0.8. The graph features solid and dashed lines in blue and black, with the dashed black lines denoting specific values.

As the length of the bridged region grows and $P_s$ becomes sufficiently large, for all $\textit{Fr}$ and $\mu _w \gt \mu _s$ , the grains will not be able to slide over both the top and bottom plates. Further, as indicated by figure 7(a), $h_{y1}$ and $h_{y2}$ evolve from 0 to 1, respectively, towards 1. As such, the velocity profile of the grains will tend to $u_s \approx 1$ , except for the increasingly thin regions around $z = 0$ and 1. Such a profile is not observed in experiments, as discussed further in § 5, and we anticipate that our continuum approximation will break down when $h_{y1}$ and $1 - h_{y2}$ become sufficiently small. When the velocity of the grains in contact with the walls is of order $\varepsilon$ , figure 7(b) shows that $\tanh (u_s/\varepsilon )$ at $z = 0$ and 1 will tend to $\mu _s / \mu _w$ , which still satisfies $u_s \ll 1$ . Thus from (2.38), it follows that $\partial _x P_s \sim -2\mu _s P_s$ , indicating that pressure builds up exponentially through the bridged region, as illustrated in figure 7(c). By comparison, when $\mu _w \lt \mu _s$ , we have the trivial velocity profile $u_s = 1$ , $\tanh (u_s / \varepsilon ) \sim 1$ and $\partial _x P_s \sim -2\mu _w P_s$ resulting in the same exponential behaviour, but not controlled by the wall friction factor. This exponential relationship between the length of the bridged region and the normal stress at the piston is consistent with previous DEM simulations (Marks et al. Reference Marks, Sandnes, Dumazer, Eriksen and Måløy2015).

Figure 7.

(a) The location of the yield layers $h_{y1}$ (dashed curve) and $h_{y2}$ (solid curve) determined by (2.31). (b) The hyperbolic tangent of the solid velocity at $z = 1$ (solid curve) and $z = 0$ (dashed line), which are used to compute the regularised shear stress according to (2.36) and (2.37), recalling that $\varepsilon = 10^{-3}$ . (c) Natural log of $|\partial _x P_s|$ computed from (2.38). Solutions are computed with $\textit{Fr} = 0.1$ and $\mu _w = 0.55$ (blue) and $\mu_w=0.85$ (red).

3.2. Transition zone

We now turn our attention to the behaviour of the transition zone, which, as illustrated in figure 1(c), is the wedge-shaped region with $h \lt 1$ and $P_s = 0$ . In this region, given that the impact of the fluid on the grains is negligible, it follows that $\sigma _{xz}(z = h) \sim 0$ , and the (non-dimensional) shear stress, (2.13), reduces to

(3.1) \begin{align} \sigma _{xz}' = (z - h) \frac {\partial h}{\partial x}. \end{align}

As illustrated in figure 2, when $\sigma _{xz}(z = h) = 0$ , the entire solid layer will be plugged ( $h_{y1} = 0$ and $h_{y2} = h$ ) if $\mu _w \lt \mu _s$ . If $\mu _w \gt \mu _s$ , then the solid layer will yield $(h_{y1} = h_{y2} = h)$ if $|\partial _x h| \gt \mu _s$ , or be plugged otherwise. As such, we can reduce (2.23) to

(3.2) \begin{align} u_s &= \begin{cases} u_{s}(z = 0) & \textrm {if } \ 0 \leqslant z \lt h_y, \\[3pt] u_{sy}(x,z,t) & \textrm {if } \ h_y \leqslant z \lt h, \end{cases} \end{align}

where

(3.3) \begin{align} h_y = \begin{cases} 0 & \textrm {if } \ |\partial _x h| \gt \mu _s, \\[3pt] h & \textrm {otherwise}. \end{cases} \end{align}

As such, we arrive at the so-called Bagnold profile (Andreotti et al. Reference Andreotti, Forterre and Pouliquen2013) in the yielded region

(3.4) \begin{align} u_{sy} = u_s(z = 0) + \frac {2 I_0}{3 \textit{Fr}} \frac {\partial _x h + \mu _s}{\partial _x h + \mu _2} \left [ (h - z)^{3/2} - (h - h_y)^{3/2} \right ] \!, \end{align}

where

(3.5) \begin{align} 0 = \frac {\partial h}{\partial x} + \mu _w \tanh \left ( \frac {u_s(z = 0) }{\varepsilon } \right ) \!. \end{align}

Further, by substituting the inertial number $I = \textit{Fr} \ \partial _z u_s / \sqrt {h-z}$ into (2.25), we arrive at

(3.6) \begin{align} I = -\frac {I_0 (\partial _x h + \mu _s)}{\partial _x h + \mu _2}, \end{align}

which combined with (2.25) reduces to $\mu (I) = -\partial _x h$ .

Recall that representative solutions were shown in figure 5(a) with $\mu _w \lt \mu _2$ , and in figure 5(b) with $\mu _2 \lt \mu _w$ . We find that as $\textit{Fr} \to 0$ (first row), the height of the granular pile tends to a line of slope $-\mu _s$ . Further, the grains on the bottom plate are stationary such that the velocity profile takes the form illustrated in figure 2(b). When $\mu _w \lt \mu _2$ (figure 5 a), as $\textit{Fr}$ increases, the length of the transition zone decreases and the grains begin to slide along the bottom plate. For sufficiently large $\textit{Fr}$ (third row), the grains appear to resist internal rearrangement such that $|\partial _z u_s| \ll 1$ . When $\mu _2 \lt \mu _w$ (figure 5 b), again the slope of the transition zone decreases with increasing $\textit{Fr}$ . However, unlike in figure 5(a), the grains remain stationary on the bottom plate for all $\textit{Fr}$ . Analogous behaviour was observed in the bridged region discussed previously. For the case where $\mu _w \lt \mu _s$ (not shown), the grains do not yield; as such, solutions are independent of $\textit{Fr}$ and the granular pile takes the shape of a wedge with $\partial _x h = -\mu _w$ .

We compute the maximum magnitude of the slope of $h$ as a function of $\textit{Fr}$ over a range of values of $\mu _w$ (figure 8 a), as well as the corresponding length $\lambda$ of the transition zone (figure 8 b), where we define the latter as the distance between the last bridged point and the first undisturbed point. We determine these values by integrating our full time-dependent model for a sufficiently long time such that the shape and velocity of the transition zone become constant. We find that for all $\mu _w \gt \mu _s$ , $\max |\partial _x h|$ – and in turn $\max |\mu (I)|$ – will tend to $\mu _s$ as $\textit{Fr} \to 0$ . This behaviour can be explained by taking $\textit{Fr} \to 0$ with $u_{s1} \sim {O}(1)$ in (3.4), which implies that $\partial _x h \sim -\mu _s$ . When $\mu _s \lt \mu _w \lt \mu _2$ , $\max |\partial _x h|$ (and $\max |\mu (I)|$ ) will monotonically increase with $\textit{Fr}$ from $\mu _s$ and plateau at $\mu _w$ . Once $\max |\partial _x h|$ has plateaued, the shear stress is sufficiently large on the base such that the grains will begin sliding, $u_s(z = 0) \gt 0$ . Given that the inertial number $I \sim \textit{Fr} \ \partial _z u_s$ does not increase above $\mu _w$ once grains begin sliding along the base, it follows that $\partial _z u_s \to 0$ as $\textit{Fr}$ becomes large, resulting in the plug-like flow observed in the third row of figure 5. Further, taking the limit that $\textit{Fr}$ is large, (3.4) and (3.5) imply that $u_s \sim u_s(z = 0)$ and $\partial _x h \sim -\mu _w$ . When $\mu _w \gt \mu _2$ , we find that $\max |\partial _x h|$ will monotonically increase from $\mu _s$ to $\mu _2$ with increasing $\textit{Fr}$ , while $\lambda$ will decrease from $\lambda \to 1/(2\mu _s)$ to $1/(2\mu _2)$ . As $\max |\partial _x h| \lt \mu _w$ for all $\textit{Fr}$ , the shear stress at the base is never sufficiently large for the grains to slide over the bottom plate and remain stationary. As such, solutions are independent of $\mu _w$ for all $\mu _w \gt \mu _2$ . In the case where $\mu _w \lt \mu _s$ , as discussed above, we find that $\max |\partial _x h| = \mu _w$ and $\lambda = 1/(2\mu _w$ ) for all $\textit{Fr}$ . Note that it is no longer reasonable to neglect the macroscopic inertia of the fluid–solid mixture once $\textit{Fr} \gtrsim O(1)$ , as discussed in § 2.3, so inertial effects should play a role for the solutions with the largest $\textit{Fr}$ here.

Figure 8.

(a) The maximum magnitude of the slope of $h$ , recalling that $I = -\partial _x h$ , computed with $h_0 = 0.5$ , $\mu _s = 0.4$ and $\mu _w = 0.3, 0.4, \ldots, 0.7$ . (b) The corresponding length $\lambda$ of the transition region.

Although $h$ will tend to a straight line in the limits such that $\textit{Fr}$ becomes sufficiently small or large, the shape of the granular pile is not trivial for intermediate values of $\textit{Fr}$ . To better understand how the shape of the transition zone varies, we solve for $h$ over a range of values of $\textit{Fr}$ (figure 9 a), and extract the corresponding slope (figure 9 b). For small $\textit{Fr}$ , the magnitude of the slope of $h$ tends to $\mu _s$ , as expected. As $\textit{Fr}$ increases, $h$ becomes slightly more nonlinear and $|\partial _x h|$ non-uniformly increases, where the slope is steeper at the start of the transition zone and decreases to 0 by the end. The magnitude of the slope increases with $\textit{Fr}$ until it reaches $\mu _w$ , which corresponds to the point at which the grains will slide along the bottom plate. We note that the shape of the transition zone is in good qualitative agreement with previous results. Specifically, Marks et al. (Reference Marks, Sandnes, Dumazer, Eriksen and Måløy2015), who performed DEM simulations of bulldozing in a confined geometry, showed that the shape of the transition zone is approximately linear with a gradient that steepens as the equivalent of the Froude number increases. Similarly, Sauret et al. (Reference Sauret, Balmforth, Caulfield and McElwaine2014), who investigated a continuum model of dry bulldozing in an unconfined geometry, showed that profiles are primarily linear in shape but can become slightly ‘curvy’, particularly around the leading edge, with increasing $\textit{Fr}$ .

Figure 9.

(a) The height of the granular pile for $h_0 = 0.5$ , $\mu _w = 0.6$ and $\textit{Fr}$ = { $0.001$ , $0.017$ , $0.036$ , $0.06$ , $0.1$ , $0.17$ , $1$ }. (b) Corresponding magnitude of the slope of $h$ . Dotted black lines denote values of $\mu _s$ , $\mu _w$ and $\mu _2$ .

We repeat this process for a fixed $\textit{Fr}$ across a range of values of $h_0$ (figure 10). For larger values of $h_0$ (darker blue curves), $|\partial _x h|$ monotonically decreases over the transition zone. For smaller values of $h_0$ (blue curves), $|\partial _x h|$ is non-monotonic such that it initially increases and reaches a peak, then decreases to 0. For sufficiently small $h_0$ , $\max |\partial _x h|$ rises above $\mu _w$ , suggesting that there will be a region where the grains slide along the base ( $|\sigma _{xz}'| \gt \mu _w p_s$ at $z = 0$ ) sandwiched between regions where the grains are stationary, $u_s(z = 0) = 0$ .

Figure 10.

(a) The height of the granular pile for $\textit{Fr} = 0.06$ , $\mu _s = 0.4$ , $\mu _w = 0.6$ and $h_0$ = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7}. (b) Corresponding magnitude of the slope of $h$ . Dotted black lines denote the values of $\mu _s$ , $\mu _w$ and $\mu _2$ .