2.1. Local equations, constitutive relations and boundary conditions

Figure 1.

Channelized granular flows over erodible beds: (a) variable definitions and velocity distribution for a typical granular column; (b) evolving profiles for a generic space- and time-varying flow; (c) the seesaw channel configuration used to produce transient uniform flows.

In an open channel of constant width $W$ , we consider a dense flow of dry grains over an erodible granular substrate (figure 1 a). At a given location $x$ and time $t$ , grains move downslope at elevation-dependent longitudinal velocity $u(z)$ , over depth $h$ between two evolving surfaces: the free surface $\widetilde{z}$ and a basal interface ${\displaystyle\mathop{z}\limits_{{\sim}}}$ below which granular shearing ceases. Granular motions are assumed to be nearly parallel to the angle of repose ${\it\alpha}$ , and coordinates $x,z$ are taken parallel and perpendicular to this direction. The flow is assumed to be uniform across the channel width. Neglecting solid fraction variations, the equation for granular mass balance is

(2.1) $$\begin{eqnarray}\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}=0,\end{eqnarray}$$

where $u$ and $w$ are the granular velocity components in the $x$ and $z$ directions, and the transverse velocity component $v$ in the $y$ direction is neglected. For smooth walls, the approximation that the velocity profile is uniform across the channel width is supported by three-dimensional ${\it\mu}(I)$ calculations (Jop et al. Reference Jop, Forterre and Pouliquen2007) and discrete particle simulations (Chou et al. Reference Chou, Lee, Chung and Hsiau2012). In the shallow flow approximation, momentum balance in the longitudinal direction yields the width-averaged equation of motion

(2.2) $$\begin{eqnarray}{\it\rho}\left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+w\frac{\partial u}{\partial z}\right)={\it\rho}g_{_{\Vert }}-\frac{\partial {\it\sigma}}{\partial x}+\frac{\partial {\it\tau}}{\partial z}-\frac{2{\it\tau}_{W}}{W},\end{eqnarray}$$

where ${\it\rho}=c_{S}{\it\rho}_{S}$ is the mean density of the bulk phase, ${\it\sigma}={\it\rho}g_{_{\bot }}(\,\widetilde{z}-z)$ is the normal stress or pressure, assumed lithostatic, $g_{_{\Vert }}=g\sin {\it\alpha}$ and $g_{_{\bot }}=g\cos {\it\alpha}$ are the parallel and perpendicular components of the gravitational acceleration, ${\it\tau}_{W}={\it\mu}_{W}{\it\sigma}$ is the shear stress along the sidewalls (with ${\it\mu}_{W}$ the wall friction coefficient) and ${\it\tau}$ is the internal shear stress. For the latter, we assume the viscoplastic dense granular flow rheology (da Cruz et al. Reference da Cruz, Emam, Prochnow, Roux and Chevoir2005)

(2.3) $$\begin{eqnarray}{\it\tau}={\it\tau}_{0}+{\it\tau}_{1}={\it\mu}_{0}{\it\sigma}+{\it\chi}D({\it\rho}{\it\sigma})^{1/2}\dot{{\it\gamma}},\end{eqnarray}$$

expressed as the sum of two components: a plastic yield stress ${\it\tau}_{0}$ proportional to ${\it\sigma}$ (with ${\it\mu}_{0}=\tan {\it\alpha}$ their constant ratio) and a viscous stress ${\it\tau}_{1}$ that varies linearly with the shear rate $\dot{{\it\gamma}}=\partial u/\partial z$ , with effective viscosity ${\it\chi}D({\it\rho}{\it\sigma})^{1/2}$ dependent on grain diameter $D$ , normal stress ${\it\sigma}$ and dimensionless rheological coefficient ${\it\chi}$ . It is a linearized version of the ${\it\mu}(I)$ rheology proposed by Jop et al. (Reference Jop, Forterre and Pouliquen2005), relating the internal friction coefficient ${\it\mu}={\it\tau}/{\it\sigma}$ to the inertial number $I=\dot{{\it\gamma}}D/({\it\sigma}/{\it\rho}_{S})^{1/2}$ . This rheology is known to fail in the quasi-static and collisional limits, corresponding respectively to low and high inertial numbers, but is well supported for slow dense flows at intermediate inertial numbers (Jop et al. Reference Jop, Forterre and Pouliquen2007; Henann & Kamrin Reference Henann and Kamrin2013). It should be noted that the constitutive relations for wall friction and internal shear stress are written assuming $u>0$ and $\dot{{\it\gamma}}>0$ . The wall friction coefficient ${\it\mu}_{W}$ , yield ratio ${\it\mu}_{0}$ and rheological coefficient ${\it\chi}$ are presumed to take constant values dependent on material properties. The equations must further be supplemented by kinematic and dynamic boundary conditions along the free surface,

(2.4a−c ) $$\begin{eqnarray}\frac{\partial \widetilde{z}}{\partial t}+\widetilde{u}\frac{\partial \widetilde{z}}{\partial x}=\widetilde{w},\quad \widetilde{{\it\sigma}}=0,\quad \widetilde{{\it\tau}}=0,\end{eqnarray}$$

and along the basal yield surface,

(2.5a−c ) $$\begin{eqnarray}{\displaystyle\mathop{ u}\limits_{{\sim}}}=0,\quad {\displaystyle\mathop{w}\limits_{{\sim}}}=0,\quad {\displaystyle\mathop{{\it\tau}}\limits_{{\sim}}}={\it\mu}_{0}{\displaystyle\mathop{{\it\sigma}}\limits_{{\sim}}}.\end{eqnarray}$$

Following Liggett (Reference Liggett1994), tildes above and below variables denote values sampled along the free surface and the basal interface respectively.

2.2. Exact solutions for equilibrium heap flows

The above equations admit exact solutions for equilibrium heap flows (Jop et al. Reference Jop, Forterre and Pouliquen2005; Berzi & Jenkins Reference Berzi and Jenkins2008). These are steady uniform flows of constant surface inclination ${\it\beta}>{\it\alpha}$ , obtained by alimenting a long open channel at a constant upstream feed rate and waiting for the flows to become stationary. In this case, the equation of motion (2.2) can be integrated to obtain the equilibrium shear stress profile ${\it\tau}({\it\eta})$ , where ${\it\eta}=\widetilde{z}-z$ is the depth below the free surface. The depth ${\displaystyle\mathop{{\it\eta}}\limits_{{\sim}}}=h$ where the basal yield condition ${\displaystyle\mathop{{\it\tau}}\limits_{{\sim}}}={\it\mu}_{0}{\displaystyle\mathop{{\it\sigma}}\limits_{{\sim}}}$ is attained determines the equilibrium flow depth

(2.6) $$\begin{eqnarray}h=h_{eq}=(W/{\it\mu}_{W})S,\end{eqnarray}$$

where $S=-\partial \widetilde{z}/\partial x=\tan ({\it\beta}-{\it\alpha})$ is the slope excess beyond the angle of repose. Once the stress profile is known, the rheology (2.3) yields an ODE for $u({\it\eta})$ that can be solved for the velocity profile,

(2.7) $$\begin{eqnarray}u(\widehat{{\it\eta}})=\left({\textstyle \frac{7}{3}}-{\textstyle \frac{35}{6}}\widehat{{\it\eta}}^{3/2}+{\textstyle \frac{7}{2}}\widehat{{\it\eta}}^{5/2}\right)\overline{u},\end{eqnarray}$$

where $\widehat{{\it\eta}}={\it\eta}/h=(\,\widetilde{z}-z)/h$ is a normalized depth coordinate and the depth-averaged velocity $\overline{u}$ is given by

(2.8) $$\begin{eqnarray}\overline{u}=\overline{u}_{eq}=\frac{4}{35}\frac{g_{_{\bot }}^{1/2}}{{\it\chi}D}\left(\frac{W}{{\it\mu}_{W}}\right)^{3/2}S^{5/2}.\end{eqnarray}$$

For small inertial number $I$ , these results agree with the more general velocity profile derived by Jop et al. (Reference Jop, Forterre and Pouliquen2005) from their nonlinear ${\it\mu}(I)$ rheology. Applied to flows without wall friction over a rough rigid floor, the same equations would yield instead the Bagnold velocity profile (see e.g. Gray & Edwards Reference Gray and Edwards2014). For equilibrium heap flows, exact solutions can thus be found for the flow layer depth $h$ , depth-averaged velocity $\overline{u}$ and velocity profile $u(z)$ , dependent on the excess inclination $S$ and other parameters. The profile shape (2.7) was earlier proposed by Takahashi (Reference Takahashi1991) for equilibrium debris flows. The profile is sigmoidal in shape and features zero shear rates at the base and top, with a maximum shear rate at $\widehat{{\it\eta}}=1/3$ . The shape matches well velocity profiles measured in debris flow experiments (Tsubaki, Hashimoto & Suetsugi Reference Tsubaki, Hashimoto and Suetsugi1982). For dry granular flows, experiments show that the velocity profiles deviate from this shape, but support the scalings implied by (2.6) and (2.8).

2.3. Depth-integrated equations

To solve the equations away from equilibrium, numerical solutions can be sought using solvers developed for free-surface flows. This has been pursued by Kamrin (Reference Kamrin2010) using an elasto-plastic formulation solved by finite elements, by Lagrée et al. (Reference Lagrée, Staron and Popinet2011) using a volume-of-fluid Navier–Stokes solver and by Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2014) using an arbitrary Lagrangian–Eulerian scheme. To make the problem more tractable and the behaviour more transparent, we wish to work instead with depth-integrated equations. Our objective is to derive reduced equations governing the evolution of depth $h(x,t)$ , basal elevation ${\displaystyle\mathop{z}\limits_{{\sim}}}(x,t)$ and mean velocity $\overline{u}(x,t)$ , for arbitrary transient and non-uniform flows (figure 1 b). To deduce the full velocity field $u(x,z,t)$ from these variables, we will assume that the equilibrium shape of the velocity profile (2.7) remains applicable away from equilibrium.

For non-equilibrium flows, the mass balance equation (2.1) can be integrated over depth to yield the continuity equation

(2.9) $$\begin{eqnarray}\frac{\partial \widetilde{z}}{\partial t}+\frac{\partial }{\partial x}(h\overline{u})=0.\end{eqnarray}$$

Likewise, integration over depth of the local equation of motion (2.2) yields, after invoking the boundary conditions and the Leibniz integral rule,

(2.10) $$\begin{eqnarray}\frac{\partial }{\partial t}(h\overline{u})+\frac{\partial }{\partial x}\left(h\overline{u^{2}}\right)=g_{_{\Vert }}h-g_{_{\bot }}h\frac{\partial \widetilde{z}}{\partial x}-\frac{1}{{\it\rho}}{\displaystyle\mathop{{\it\tau}}\limits_{{\sim}}}-\frac{2}{{\it\rho}W}h\overline{{\it\tau}}_{W}.\end{eqnarray}$$

Here and elsewhere, overbars denote depth averages, defined for dummy variable ${\it\xi}$ by $\overline{{\it\xi}}=1/h\int _{{\displaystyle\mathop{z}\limits_{{\sim}}}}^{\widetilde{z}}{\it\xi}(z)\text{d}z$ . Assuming that the velocity profile remains similar to the equilibrium profile even for non-equilibrium flows, the shape function (2.7) can be used to work out the momentum flux $h\overline{u^{2}}=\int _{{\displaystyle\mathop{z}\limits_{{\sim}}}}^{\widetilde{z}}u^{2}\text{d}z=(77/48)h\overline{u}^{2}$ . Substituting this result and the constitutive relations for ${\it\tau}$ and ${\it\tau}_{W}$ , we obtain the depth-integrated momentum balance equation in the explicit form

(2.11) $$\begin{eqnarray}\frac{\partial }{\partial t}(h\overline{u})+\frac{\partial }{\partial x}\left(\frac{77}{48}h\overline{u}^{2}\right)=-g_{_{\bot }}h\frac{\partial \widetilde{z}}{\partial x}-\frac{{\it\mu}_{W}}{W}g_{_{\bot }}h^{2},\end{eqnarray}$$

where the terms $g_{_{\Vert }}h$ and ${\displaystyle\mathop{{\it\tau}_{0}}\limits_{{\sim}}}/{\it\rho}$ have cancelled each other out thanks to our choice of reference inclination ${\it\alpha}$ such that $\tan {\it\alpha}={\it\mu}_{0}$ . Equations similar to (2.9) and (2.11) have earlier been proposed to describe flows over rigid boundaries (Savage & Hutter Reference Savage and Hutter1991; Gray & Edwards Reference Gray and Edwards2014) and erodible beds (Tai & Kuo Reference Tai and Kuo2008; Lê & Pitman Reference Lê and Pitman2010).

Over rigid boundaries, depth-integrated equations obtained in this way provide a complete macroscopic description of the dynamics of open-channel flows. Over erodible beds, however, the description suffers from two blind spots. First, evolution of the two variables $\widetilde{z}={\displaystyle\mathop{ z}\limits_{{\sim}}}+h$ and $q=h\overline{u}$ does not suffice to determine the three degrees of freedom ${\displaystyle\mathop{z}\limits_{{\sim}}}$ , $h$ and $\overline{u}$ . Recent works have sought to address this problem by complementing (2.9) and (2.10) with an additional closure relation for the erosion rate $e=-\partial {\displaystyle\mathop{z}\limits_{{\sim}}}\,/\partial t$ (Tai & Kuo Reference Tai and Kuo2008; Lê & Pitman Reference Lê and Pitman2010; Iverson Reference Iverson2012). In order to deduce rates of erosion or entrainment $e$ from first principles, we wish to avoid this tactic. The second problem with (2.10) is that, after depth integration, only the basal internal shear stress ${\displaystyle\mathop{{\it\tau}}\limits_{{\sim}}}$ remains. Because the shear rate vanishes at the base, moreover, so does the rate-dependent component of the shear stress there, ${\displaystyle\mathop{{\it\tau}_{1}}\limits_{{\sim}}}=0$ . Hence nowhere in (2.9) and (2.10) does ${\it\tau}_{1}$ play any role. As a result, the equations cannot recover the equilibrium mean velocity (2.8), controlled by the rate-dependent shear stress via the rheological parameter ${\it\chi}$ .

To address these problems, we propose to complement (2.9) and (2.10) with an additional evolution equation expressing the balance of kinetic energy. We derive this equation by integrating (2.2) once more, this time weighted by the local velocity $u(z)$ , to obtain

(2.12) $$\begin{eqnarray}\frac{\partial }{\partial t}\left(\frac{1}{2}h\overline{u^{2}}\right)+\frac{\partial }{\partial x}\left(\frac{1}{2}h\overline{u^{3}}\right)=g_{_{\Vert }}h\overline{u}-g_{_{\bot }}h\overline{u}\frac{\partial \widetilde{z}}{\partial x}-\frac{1}{{\it\rho}}h\overline{{\it\tau}\dot{{\it\gamma}}}-\frac{2}{{\it\rho}W}h\overline{{\it\tau}_{W}u},\end{eqnarray}$$

where boundary conditions and the Leibniz integral rule have again been invoked. In this equation, the left-hand terms give the rate of change of the depth-integrated kinetic energy, and the divergence of the kinetic energy flux. The right-hand terms are the gravitational production, the work of the normal stress gradient and the rates of energy dissipation by internal and wall shear stresses. It should be noted that (2.2) could be multiplied by any weight function ${\it\omega}(z)$ prior to integration over depth. The alternative weight function ${\it\omega}(z)=z-(\,\widetilde{z}-(1/2)h)$ , for example, could be used to derive an equation for the balance of angular momentum about mid-depth (Steffler & Jin Reference Steffler and Jin1993). As in Richard & Gavrilyuk (Reference Richard and Gavrilyuk2013), however, the kinetic energy equation produced by the choice ${\it\omega}(z)=u(z)$ turns out to be particularly tractable and straightforward to interpret. Upon substituting the constitutive relations and evaluating integrals using the velocity profile (2.7), all terms of (2.12) can be explicitly worked out to obtain

(2.13) $$\begin{eqnarray}\frac{\partial }{\partial t}\left(\frac{77}{96}h\overline{u}^{2}\right)+\frac{\partial }{\partial x}\left({\it\kappa}h\overline{u}^{3}\right)=-g_{_{\bot }}h\overline{u}\frac{\partial \widetilde{z}}{\partial x}-\frac{35}{9}{\it\chi}Dg_{\bot }^{1/2}\frac{\overline{u}^{2}}{h^{1/2}}-\frac{5}{9}\frac{{\it\mu}_{W}}{W}g_{_{\bot }}h^{2}\overline{u},\end{eqnarray}$$

where ${\it\kappa}=(342\,853/233\,376)\approx 1.469$ . This one additional balance equation turns out to solve the different problems. By tracking the kinetic energy $k=(77/96)h\overline{u}^{2}$ in addition to the momentum $q=h\overline{u}$ and surface elevation $\widetilde{z}$ , we can retrieve all three degrees of freedom $\overline{u}=(96/77)k/q$ , $h=q/\overline{u}$ and ${\displaystyle\mathop{z}\limits_{{\sim}}}=\widetilde{z}-h$ . The rate-dependent shear stress ${\it\tau}_{1}$ , moreover, now intervenes via the depth-averaged dissipation term $\overline{{\it\tau}_{1}\dot{{\it\gamma}}}=(35/9){\it\rho}{\it\chi}Dg_{\bot }^{1/2}\overline{u}^{2}h^{-3/2}$ . As a result, application of (2.11) and (2.13) to equilibrium heap flows yields the expected depth (2.6) and mean velocity (2.8). The extended equations thus admit the correct equilibrium state.