2.1. Mixture framework

Continuum mixture theory (Truesdell Reference Truesdell1984; Morland Reference Morland1992) provides a framework to describe the motion of multi-phase materials and postulates that each spatial point is simultaneously occupied by all the phases. Debris flows may be composed of grains, water and air, but for simplicity, this paper uses a two-phase formulation that neglects the effect of air. The volume fraction $\phi ^\nu$ of constituent $\nu$ determines what proportion of the mixture volume is occupied by that constituent. Using the constituent letters $\nu =w,g$ for water and grains, respectively, it follows that for a pure water phase, $\phi ^w=1$, while for a dry granular phase, the solids volume fraction $\phi ^g<1$ and $\phi ^w=0$. However, for a water saturated flow, the volume fractions sum to unity

(2.1)\begin{equation} \phi^g+\phi^w=1. \end{equation}

Mixture theory defines overlapping partial density $\varrho ^{\nu }$, partial velocity $\boldsymbol {u}^\nu$ and partial stress $\boldsymbol {\sigma }^\nu$ fields for each phase. Each distributed phase satisfies its individual conservation laws for mass,

(2.2)\begin{equation} \frac{\partial\varrho^\nu}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\varrho^\nu\boldsymbol{u}^\nu) = 0, \end{equation}

and for momentum,

(2.3)\begin{equation} \frac{\partial}{\partial t}(\varrho^\nu\boldsymbol{u}^\nu)+\boldsymbol{\nabla} \boldsymbol{\cdot}(\varrho^\nu\boldsymbol{u}^\nu\otimes\boldsymbol{u}^\nu) = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}^\nu + \varrho^\nu\boldsymbol{g} + \boldsymbol{\beta}^\nu, \end{equation}

where $t$ is time, $\boldsymbol {\nabla }$ is the gradient operator, $\boldsymbol {\cdot }$ is the dot product, $\otimes$ is the dyadic product and $\boldsymbol {g}$ the acceleration due to gravity. The interaction drag $\boldsymbol {\beta }^\nu$ is the force exerted on the $\nu$ phase by the other phase, and satisfies $\boldsymbol {\beta }^g+\boldsymbol {\beta }^w=0$ according to Newton's third law.

The partial density is related to its intrinsic counterpart (defined on a unit volume of that constituent) by a linear volume fraction scaling, while the partial and intrinsic velocity fields are identical

(2.4a,b)\begin{equation} \varrho^\nu=\phi^\nu\varrho^{\nu\star},\quad \boldsymbol{u}^\nu=\boldsymbol{u}^{\nu\star}, \end{equation}

where the superscript $\star$ represents an intrinsic quantity. The intrinsic pore fluid pressure in the water is $p^{w\star }$. Following de Boer & Ehlers (Reference de Boer and Ehlers1990), the partial stress in the grains is decomposed into

(2.5)\begin{equation} \boldsymbol{\sigma}^g={-}\phi^gp^{w\star}\textbf{1}-\boldsymbol{\sigma}^{e}, \end{equation}

where $\boldsymbol {\sigma }^{e}$ is the effective stress (which, by the convention used in soil mechanics, has the opposite sign to that of Cauchy stress) and $\textbf {1}$ is the unit tensor. The partial water stress is assumed to be

(2.6)\begin{equation} \boldsymbol{\sigma}^w={-}\phi^wp^{w\star}\textbf{1}+\boldsymbol{\tau}^w, \end{equation}

where the deviatoric stress $\boldsymbol {\tau }^w$ satisfies a Newtonian fluid law

(2.7)\begin{equation} \boldsymbol{\tau}^w=\phi^w\eta^w(\boldsymbol{\nabla}\boldsymbol{u}^w+ (\boldsymbol{\nabla}\boldsymbol{u}^w)^\text{T}), \end{equation}

in which $\eta ^w$ is the dynamic viscosity and $\textrm {T}$ is the transpose. This assumes that the interstitial fluid can be described as a Newtonian fluid with constant viscosity, which is a potential source of discrepancy. For instance, in Baumgarten & Kamrin's (Reference Baumgarten and Kamrin2019) two-phase fluid–sediment mixture model, Einstein's (Reference Einstein1906) linear solids volume correction is used to account for the increased viscosity of the fluid due to the suspended sediment. This idea can be extended to higher solids volume fractions. For instance, Boyer, Guazzelli & Pouliquen (Reference Boyer, Guazzelli and Pouliquen2011) showed how to relate the frictional behaviour of dense suspensions with classical viscous suspension rheology, so that their behaviour could be expressed either as a function of the solids volume fraction or the viscous number $I_\nu =J=\eta ^w\dot \gamma /p^g$, where $\dot \gamma$ is the shear rate and $p^g$ is the pressure due to grain contacts. This suspension rheology works well at low shear rates, but at higher shear rates, it has been found that there is a transition towards a dependence on the granular inertial number $I=\dot \gamma d/\sqrt {p^g/\varrho ^{g\star }}$ (Trulsson, Andreotti & Claudin Reference Trulsson, Andreotti and Claudin2012; Maurin, Chauchat & Frey Reference Maurin, Chauchat and Frey2016). While these suspension rheologies determine the composite response for a single effective medium, it is by no means clear how to partition the stress back between the individual phases.

Equations (2.6) and (2.7) are consistent with the mixture theory approach of Nunziato et al. (Reference Nunziato, Passman, Givler, Mactigue and Brady1986), although an ensemble averaging approach is also possible (Anderson & Jackson Reference Anderson and Jackson1967; Joseph & Lungdren Reference Joseph and Lungdren1990; Anderson, Sundaresan & Jackson Reference Anderson, Sundaresan and Jackson1995; Enwald, Pelrano & Almstedt Reference Enwald, Pelrano and Almstedt1996; Jackson Reference Jackson2000). The sum of the partial solid and fluid stresses, (2.5) and (2.6), determines the total stress

(2.8)\begin{equation} \boldsymbol{\sigma}^{total}={-}p^{w\star}\textbf{1}-\boldsymbol{\sigma}^{e}+\boldsymbol{\tau}^w, \end{equation}

which is consistent with the effective stress principle of soil mechanics (Terzaghi Reference Terzaghi1943). It also agrees with the experimental observation that the manometric pressure in the soil is the pressure as if the medium were a bulk fluid, unaffected by the presence of the solid constituent in the medium.

The interaction drags acting on the water and grain phases are defined by

(2.9)\begin{equation} \boldsymbol{\beta}^w={-}\boldsymbol{\beta}^g=p^{w\star}\boldsymbol{\nabla} \phi^w+C^d(\boldsymbol{u}^g-\boldsymbol{u}^w), \end{equation}

where the term $p^{w\star }\boldsymbol {\nabla }\phi ^w$ combines with the partial water pressure gradient $-\boldsymbol {\nabla }(\phi ^w p^{w\star })$ in the water momentum balance (2.3) to ensure that the fluid is driven by gradients in the intrinsic pore water pressure $-\phi ^w\boldsymbol {\nabla } p^{w\star }$, consistent with Darcy's law (Bear Reference Bear1972; Morland Reference Morland1992). The second term on the right-hand side of (2.9) is the drag due to the relative motion of the grains and water, with drag coefficient

(2.10)\begin{equation} C^d=\eta^w\frac{(\phi^w)^2}{k}, \quad\hbox{where the permeability}\ k=\frac{(\phi^w)^3d^2}{180(\phi^g)^2}, \end{equation}

is given by Carman's formula for the packing of spheres of diameter $d$. This agrees well with the sediment dynamics experiments of Goharzadeh, Khalili & Jørgensen (Reference Goharzadeh, Khalili and Jørgensen2005) and Ouriemi, Aussillous & Guazzelli (Reference Ouriemi, Aussillous and Guazzelli2009).

When the constitutive laws (2.5) and (2.6) and the interaction drags (2.9) are substituted into (2.3), the momentum balance equations for the grains and the fluid take the form

(2.11)$$\begin{gather} \frac{\partial}{\partial t}(\varrho^g\boldsymbol{u}^g)+ \boldsymbol{\nabla}\boldsymbol{\cdot}(\varrho^g\boldsymbol{u}^g\otimes\boldsymbol{u}^g) ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}^e - \phi^g\boldsymbol{\nabla} p^{w\star}+ \varrho^g\boldsymbol{g} + C^d(\boldsymbol{u}^w-\boldsymbol{u}^g), \end{gather}$$

(2.12)$$\begin{gather}\frac{\partial}{\partial t}(\varrho^w\boldsymbol{u}^w)+\boldsymbol{\nabla}\boldsymbol{\cdot} (\varrho^w\boldsymbol{u}^w\otimes\boldsymbol{u}^w)={-}\phi^w\boldsymbol{\nabla} p^{w\star}+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau}^w+\varrho^w \boldsymbol{g}+C^d(\boldsymbol{u}^g-\boldsymbol{u}^w), \end{gather}$$

assuming that the mixture is fully saturated. When there is no water ($\phi ^w=0$), the buoyancy force $-\phi ^g\boldsymbol {\nabla } p^{w\star }$ and the Darcy interaction force $C^d(\boldsymbol {u}^w-\boldsymbol {u}^g)$ are assumed to be zero, and the momentum balance (2.11) then reduces to that for dry grains. Conversely, in the case of pure water ($\phi ^w=1$), (2.12) naturally reduces to the Navier–Stokes equations. The momentum balances (2.11) and (2.12) are identical to those employed to investigate debris flows by Iverson & Denlinger (Reference Iverson and Denlinger2001), sediment flows by Ouriemi et al. (Reference Ouriemi, Aussillous and Guazzelli2009) and general fluid–sediment flows by Baumgarten & Kamrin (Reference Baumgarten and Kamrin2019), apart from the particular choices of the Darcy drag coefficients and effective fluid phase viscosity.

2.2. Boundary conditions

A Cartesian coordinate system $Oxz$ is introduced with the $x$-axis pointing down a slope inclined at an angle $\zeta$ to the horizontal and with the $z$-axis pointing upwards, as shown in figure 2. A third coordinate $y$, lying across the slope, can easily be added to model fully three-dimensional flows, but for ease of exposition, the theory presented here is purely two-dimensional. The debris flow is assumed to have velocity components $\boldsymbol {u}^\nu =(u^\nu, w^\nu )$ in the downslope and slope normal directions, respectively. The grain surface is defined by the function $F^g(z,x,t)=z-s^g(x,t)$, the water surface by $F^w(z,x,t)=z-s^w(x,t)$ and the basal surface by $F^b(z,x,t)=z-b(x)$. It follows that the upward pointing unit normal for each surface is $\boldsymbol {n}^\nu =\boldsymbol {\nabla } F^\nu /\vert \boldsymbol {\nabla } F^\nu \vert$, where $\nu =g, w, b$. In the undersaturated regime (figure 2b), the water surface is below the grain surface, i.e. $s^w< s^g$, whereas in the oversaturated regime (figure 2c), the water surface is above it, i.e. $s^g< s^w$.

Figure 2.

(a) Two-dimensional schematic diagram of a debris flow on a slope inclined at an angle $\zeta$ to the horizontal. Blue shading corresponds to water, while the grains occupy the region below the red free surface. Velocity shear results in the surface layer of (light) grains moving towards the front, while the (dark) grains near the base move backwards relative to the front. This creates a recirculating (green) frontal cell. The water free surface (blue) and the granular free surface (red) do not coincide. The debris flow therefore consists of a dry granular front, an undersaturated wet flow, an oversaturated wet flow and a pure watery tail. It is useful to resolve the vertical structure as shown in panels (b,c) for the undersaturated and oversaturated regimes. In panel (b), the granular free surface at $z=s^g(x,t)$ is above the water free surface at $z=s^w(x,t)$, while in the oversaturated regime, the situation is reversed. The base of the flow lies at $z=b(x)$ and the water and grain heights are $h^w=s^w-b$ and $h^g=s^g-b$. The internal interface height is $\textrm {min}(h^g, h^w)$ and the proportion of the water flow height that is occupied by grains is equal to $\mathcal {H}=\textrm {min}(h^g, h^w)/h^w$. Within this vertical structure, region (i) consists of dry grains, region (ii) is a mixture of grains and water, and region (iii) is pure water.

The grain phase is subject to kinematic conditions on the grain surface $F^g(x, t)=0$ and the basal surface $F^b(x)=0$. These conditions are given by

(2.13)$$\begin{gather} \frac{\partial F^g}{\partial t} + \boldsymbol{u}^g\boldsymbol{\cdot}\boldsymbol{\nabla} F^g = 0,\quad F^g(x, t)=0, \end{gather}$$

(2.14)$$\begin{gather}\frac{\partial F^b}{\partial t} + \boldsymbol{u}^g\boldsymbol{\cdot}\boldsymbol{\nabla} F^b = 0,\quad F^b(x)=0. \end{gather}$$

Similarly, the water phase is subject to kinematic conditions on the water surface $F^w(x, t)=0$ and the basal surface $F^b(x)=0$. They are given by

(2.15)$$\begin{gather} \frac{\partial F^w}{\partial t} + \boldsymbol{u}^w\boldsymbol{\cdot}\boldsymbol{\nabla} F^w = 0,\quad F^w(x, t)=0, \end{gather}$$

(2.16)$$\begin{gather}\frac{\partial F^b}{\partial t} + \boldsymbol{u}^w\boldsymbol{\cdot}\boldsymbol{\nabla} F^b = 0,\quad F^b(x)=0. \end{gather}$$

Grains can penetrate the water surface in the undersaturated regime (figure 2b) and water can penetrate the grain surface in the oversaturated regime (figure 2c). The material that crosses these surfaces can change its volume fraction abruptly, e.g. in the oversaturated regime, the volume fraction of water jumps from $\phi ^w=1-\phi ^g$ to $\phi ^w=1$ across the interface. To describe such singular interfaces, the mass jump conditions (Chadwick Reference Chadwick1976) are

(2.17)$$\begin{gather} \unicode{x27E6}\varrho^g(\boldsymbol{u}^g-\boldsymbol{u}^w)\boldsymbol{\cdot}\boldsymbol{n}^w \unicode{x27E7} = 0,\quad F^w(x, t)=0, \end{gather}$$

(2.18)$$\begin{gather}\unicode{x27E6}\varrho^w(\boldsymbol{u}^w-\boldsymbol{u}^g)\boldsymbol{\cdot}\boldsymbol{n}^g \unicode{x27E7} = 0,\quad F^g(x, t)=0, \end{gather}$$

where the jump brackets $\unicode{x27E6} ~\unicode{x27E7}$ indicate the difference between the enclosed quantity on the forward ($+$) and rearward ($-$) sides of the discontinuity. Equation (2.17) describes the conservation of grains across the water free surface in the undersaturated case, while (2.18) describes the conservation of water across the granular free surface in the oversaturated regime.

Grains and water satisfy traction-free conditions at their free surfaces

(2.19)$$\begin{gather} \boldsymbol{\sigma}^e\boldsymbol{n}^g=\boldsymbol{0},\quad F^g(x, t)=0, \end{gather}$$

(2.20)$$\begin{gather}\boldsymbol{\sigma}^w\boldsymbol{n}^w=\boldsymbol{0},\quad F^w(x, t)=0. \end{gather}$$

Iverson (Reference Iverson2003) used large-scale debris flow flume tests to show that the grain phase experiences a Coulomb-type friction. Additionally, the Chézy formula (Chaudhry Reference Chaudhry2008), an empirically derived expression to take into account turbulent friction arising from the bottom of the channel, can be employed to describe the basal friction experienced by the water phase. The solid bed friction and Chézy formula take the form

(2.21)$$\begin{gather} \boldsymbol{\sigma}^e\boldsymbol{n}^b-(\boldsymbol{n}^b\boldsymbol{\cdot} \boldsymbol{\sigma}^e\boldsymbol{n}^b)\boldsymbol{n}^b ={-}\frac{\boldsymbol{u}^g_b}{|{\boldsymbol{u}^g_b}|}\mu^b(\boldsymbol{n}^b \boldsymbol{\cdot}\boldsymbol{\sigma}^e\boldsymbol{n}^b),\quad F^b(x)=0, \end{gather}$$

(2.22)$$\begin{gather}\boldsymbol{\sigma}^w\boldsymbol{n}^b-(\boldsymbol{n}^b\boldsymbol{\cdot} \boldsymbol{\sigma}^w\boldsymbol{n}^b)\boldsymbol{n}^b =\varrho^{w\star}C^w|{\bar{\boldsymbol{u}}}^{w}|\bar{\boldsymbol{u}}^{w},\quad F^b(x)=0, \end{gather}$$

where $\mu ^b$ is the basal friction coefficient, $C^w$ is the Chézy drag coefficient (Hogg & Pritchard Reference Hogg and Pritchard2004) and $\bar {\boldsymbol {u}}^w$ is the depth-averaged water velocity. Detailed formulae for the coefficients $\mu ^b$ and $C^w$ will be discussed after the depth-integration process.

3.1. Scaling the governing equations

Debris flows typically have a flow depth $H$ that is much smaller than their length $L$. This makes it possible to derive a depth-averaged model, which does not depend on the normal coordinate $z$. Typical downstream flow speeds for both the grain and water phases are assumed to be of the order of the gravity wave speed $U=(gH)^{1/2}$, and the mass balances then imply that typical normal velocities are of magnitude $\varepsilon U$, where the aspect ratio $\varepsilon =H/L\ll 1$. Hydrostatic and lithostatic balances for the pore fluid pressure and the grain normal stresses suggest scalings of $\varrho ^{w\star }gH$ and $\varrho ^{g\star }gH$, respectively. A Coulomb friction law then suggests that the typical order of magnitude for the grain shear stress is $\mu \varrho ^{g\star }gH$, where $\mu$ is the friction coefficient. The scaling of the deviatoric water stress is determined by the Newtonian constitutive law (2.7). The governing equations and boundary conditions are therefore non-dimensionalised by the scalings

(3.1a–f)\begin{equation} \left.\begin{gathered} (x,z,s^w, s^g, b, d)=L(\hat{x}, varepsilon\hat{z}, \varepsilon\hat{s}^w, \varepsilon\hat{s}^g, \varepsilon\hat{b}, \varepsilon\hat{d}),\quad t=L/(gH)^{1/2}\hat{t}, \\ (u^\nu, w^\nu, \vert\boldsymbol{u}^\nu\vert)=(gH)^{1/2}(\hat{u}^{\nu}, \varepsilon\hat{w}^{\nu}, \vert\hat{\boldsymbol{u}}^\nu\vert),\quad C^d=\varrho^{w\star}(g/H)^{1/2}\hat C^d,\\ (\sigma^e_{xx}, \sigma^e_{zz},\sigma^e_{xz})=\varrho^{g\star}gH(\hat{\sigma}^e_{xx}, \hat{\sigma}^e_{zz}, \mu\hat{\sigma}^e_{xz}) ,\\(p^{w\star}, \tau_{xx}^w, \tau_{zz}^w, \tau_{xz}^w) =\varrho^{w\star}gH(\hat{p}^{w\star}, \varepsilon\hat{\tau}_{xx}^w, \varepsilon\tau_{zz}^w, \hat{\tau}_{xz}^w), \end{gathered}\right\} \end{equation}

where the hatted quantities are non-dimensional. These scalings imply that for the grains, the non-dimensional mass and momentum balance components are

(3.2)\begin{align} &\hspace{5pc}\frac{\partial\phi^g}{\partial\hat{t}}+\frac{\partial}{\partial \hat{x}}(\phi^g\hat{u}^g)+ \frac{\partial}{\partial \hat{z}}(\phi^g\hat{w}^g)=0, \end{align}

(3.3)\begin{align} & \varepsilon\left(\frac{\partial}{\partial\hat{t}}(\phi^g\hat{u}^g)+ \frac{\partial}{\partial \hat{x}}(\phi^g\hat{u}^g\hat{u}^g)+ \frac{\partial}{\partial \hat{z}}(\phi^g\hat{u}^g\hat{w}^g)\right) \nonumber\\ &\quad ={-}\varepsilon\gamma\phi^g\frac{\partial \hat{p}^{w\star}}{\partial \hat{x}}-\varepsilon \frac{\partial\hat{\sigma}^e_{xx}}{\partial\hat{x}} -\mu\frac{\partial\hat{\sigma}^e_{xz}}{\partial\hat{z}}+\phi^g\sin\zeta+ \gamma\hat{C}^d(\hat{u}^w-\hat{u}^g), \end{align}

(3.4)\begin{align} & \varepsilon^2\left(\frac{\partial}{\partial\hat{t}}(\phi^g\hat{w}^g)+ \frac{\partial}{\partial \hat{x}}(\phi^g\hat{u}^g\hat{w}^g)+ \frac{\partial}{\partial \hat{z}}(\phi^g\hat{w}^g\hat{w}^g)\right) \nonumber\\ &\quad ={-}\gamma\phi^g\frac{\partial \hat{p}^{w\star}}{\partial \hat{z}}-\varepsilon\mu \frac{\partial\hat{\sigma}^e_{xz}}{\partial\hat{x}} -\frac{\partial\hat{\sigma}^e_{zz}}{\partial\hat{z}} -\phi^g\cos\zeta+\varepsilon\gamma\hat{C}^d(\hat{w}^w-\hat{w}^g), \end{align}

where the density ratio

(3.5)\begin{equation} \gamma = \varrho^{w\star}/\varrho^{g\star}. \end{equation}

Similarly, the non-dimensional water mass balance and the downslope and normal components of the momentum balance are given by

(3.6)\begin{align} & \hspace{5pc}\frac{\partial\phi^w}{\partial\hat{t}}+\frac{\partial}{\partial \hat{x}}(\phi^w\hat{u}^w)+ \frac{\partial}{\partial \hat{z}}(\phi^w\hat{w}^w)=0, \end{align}

(3.7)\begin{align} & \varepsilon\left(\frac{\partial}{\partial\hat{t}}(\phi^w\hat{u}^w)+ \frac{\partial}{\partial \hat{x}}(\phi^w\hat{u}^w\hat{u}^w)+ \frac{\partial}{\partial \hat{z}}(\phi^w\hat{u}^w\hat{w}^w)\right) \nonumber\\ &\quad ={-}\varepsilon\phi^w\frac{\partial \hat{p}^{w\star}}{\partial \hat{x}}+\varepsilon^2 \frac{\partial\hat{\tau}^w_{xx}}{\partial\hat{x}}+\frac{\partial\hat{\tau}^w_{xz}}{\partial\hat{z}} +\phi^w\sin\zeta-\hat{C}^d(\hat{u}^w-\hat{u}^g), \end{align}

(3.8)\begin{align} & \varepsilon^2\left(\frac{\partial}{\partial\hat{t}}(\phi^w\hat{w}^w)+ \frac{\partial}{\partial \hat{x}}(\phi^w\hat{u}^w\hat{w}^w)+ \frac{\partial}{\partial \hat{z}}(\phi^w\hat{w}^w\hat{w}^w)\right) \nonumber\\ &\quad ={-}\phi^w\frac{\partial \hat{p}^{w\star}}{\partial \hat{z}}+\varepsilon \frac{\partial\hat{\tau}^w_{xz}}{\partial\hat{x}}+\varepsilon\frac{\partial\hat{\tau}^w_{zz}}{\partial\hat{z}} -\phi^w\cos\zeta-\varepsilon\hat{C}^d(\hat{w}^w-\hat{w}^g). \end{align}

3.2. Non-dimensional boundary conditions

The non-dimensional kinematic conditions (2.13) and (2.14) for the grains are

(3.9)$$\begin{gather} \frac{\partial \hat{s}^g}{\partial\hat{t}} + \hat{u}^g\frac{\partial\hat{s}^g}{\partial\hat{x}}-\hat{w}^g=0,\quad \hat{z}=\hat{s}^g(\hat{x},\hat{t}), \end{gather}$$

(3.10)$$\begin{gather}\frac{\partial \hat{b}}{\partial\hat{t}} + \hat{u}^g\frac{\partial\hat{b}}{\partial\hat{x}}-\hat{w}^g=0,\quad \hat{z}=\hat{b}(\hat{x}), \end{gather}$$

and the non-dimensional kinematic conditions (2.15) and (2.16) for the water are

(3.11)$$\begin{gather} \frac{\partial \hat{s}^w}{\partial\hat{t}} + \hat{u}^w\frac{\partial\hat{s}^w}{\partial\hat{x}}-\hat{w}^w=0,\quad \hat{z}=\hat{s}^w(\hat{x},\hat{t}), \end{gather}$$

(3.12)$$\begin{gather}\frac{\partial \hat{b}}{\partial\hat{t}} + \hat{u}^w\frac{\partial\hat{b}}{\partial\hat{x}}-\hat{w}^w=0,\quad \hat{z}=\hat{b}(\hat{x}). \end{gather}$$

To formulate the non-dimensional form of the mass-jump conditions (2.17) and (2.18), the kinematic conditions (2.15) and (2.13) are used to show that $\boldsymbol {u}^w\boldsymbol {\cdot }\boldsymbol {n}^w=-(\partial F^w/\partial t){\cdot} (1/\vert \boldsymbol {\nabla } F^w\vert )$ and $\boldsymbol {u}^g\boldsymbol {\cdot }\boldsymbol {n}^g=-(\partial F^g/\partial t)\cdot (1/\vert \boldsymbol {\nabla } F^g\vert )$. Substitution of these results implies that the non-dimensional jump conditions at the undersaturated water free surface and the oversaturated grain free surface are

(3.13)$$\begin{gather} \left[\!\!\left[ \phi^{g}\left(\frac{\partial\hat{s}^w}{\partial\hat{t}}+\hat{u}^{g} \frac{\partial\hat{s}^w}{\partial\hat{x}}-\hat{w}^{g}\right)\right]\!\!\right]=0,\quad \hat{z}=\hat{s}^w(\hat x,\hat t), \end{gather}$$

(3.14)$$\begin{gather}\left[\!\!\left[\phi^{w}\left(\frac{\partial\hat{s}^g}{\partial\hat{t}}+\hat{u}^{w} \frac{\partial\hat{s}^g}{\partial\hat{x}}-\hat{w}^{w}\right)\right]\!\!\right]=0,\quad \hat{z}=\hat{s}^g(\hat x,\hat t). \end{gather}$$

The non-dimensional downslope and normal components of the grain surface traction (2.19) are

(3.15)$$\begin{gather} -\varepsilon\hat{\sigma}^e_{xx}\frac{\partial\hat{s}^g}{\partial\hat{x}}+\mu\hat{\sigma}_{xz}^e=0,\quad \hat{z}=\hat{s}^g(\hat{x},\hat{t}), \end{gather}$$

(3.16)$$\begin{gather}-\varepsilon\mu\hat{\sigma}^e_{xz}\frac{\partial\hat{s}^g}{\partial\hat{x}}+\hat{\sigma}^e_{zz}=0,\quad \hat{z}=\hat{s}^g(\hat{x},\hat{t}), \end{gather}$$

and the non-dimensional components of the water surface traction (2.20) are

(3.17)$$\begin{gather} \varepsilon\phi^w\hat{p}^{w\star}\frac{\partial\hat{s}^w}{\partial\hat{x}}-\varepsilon^2\hat{\tau}^w_{xx} \frac{\partial\hat{s}^w}{\partial\hat{x}}+\hat{\tau}^w_{xz}=0,\quad \hat{z}=\hat{s}^w(\hat{x},\hat{t}), \end{gather}$$

(3.18)$$\begin{gather}-\phi^w\hat{p}^{w\star}-\varepsilon\hat{\tau}^w_{xz}\frac{\partial\hat{s}^w}{\partial\hat{x}}+ \varepsilon\hat{\tau}^w_{zz}=0,\quad \hat{z}=\hat{s}^w(\hat{x},\hat{t}). \end{gather}$$

The downslope and normal components of the basal granular traction (2.21) are

(3.19)$$\begin{gather} -\varepsilon\hat{\sigma}_{xx}^e\frac{\partial\hat{b}}{\partial\hat{x}} + \mu\hat{\sigma}^e_{xz} ={-}(\boldsymbol{n}^b\boldsymbol{\cdot}\hat{\boldsymbol{\sigma}}^e\boldsymbol{n}^b) \left(\frac{\hat{u}_b^g}{\vert\hat{\boldsymbol{u}}_b^g\vert}\mu^b|\boldsymbol{\nabla} F^b|+\varepsilon \frac{\partial\hat{b}}{\partial\hat{x}}\right),\quad \hat{z}=\hat{b}(\hat{x}), \end{gather}$$

(3.20)$$\begin{gather}-\varepsilon\mu\hat{\sigma}^e_{zx}\frac{\partial\hat{b}}{\partial\hat{x}} + \hat{\sigma}^e_{zz} ={-}(\boldsymbol{n}^b\boldsymbol{\cdot}\hat{\boldsymbol{\sigma}}^e\boldsymbol{n}^b) \left(\varepsilon\frac{\hat{w}^g_b}{|{\hat{\boldsymbol{u}}^g_b}|}\mu^b|\boldsymbol{\nabla} F^b|-1\right),\quad \hat{z}=\hat{b}(\hat{x}), \end{gather}$$

where $|\hat {\boldsymbol {u}}_b^g|=\sqrt {(\hat {u}_b^g)^2 + \varepsilon ^2(\hat {w}_b^g)^2}$ and $|\boldsymbol {\nabla } F^b|=\sqrt {1+\varepsilon ^2(\partial \hat b/\partial \hat x)^2}$. Similarly, the downslope and normal components of the water basal friction (2.22) are

(3.21)$$\begin{gather} -\varepsilon^2\hat{\tau}_{xx}^w\frac{\partial\hat{b}}{\partial\hat{x}}+\hat{\tau}_{xz}^w+\varepsilon \left(\phi^w\hat{p}^{w\star}_b+(\boldsymbol{n}^b\boldsymbol{\cdot} \hat{\boldsymbol{\sigma}}^w\boldsymbol{n}^b)\right) \frac{\partial\hat{b}}{\partial\hat{x}}=C^w\hat{\bar{u}}^w|{\hat{\bar{\boldsymbol{u}}}^w}\|\boldsymbol{\nabla} F^b|,\quad \hat{z}=\hat{b}(\hat{x}), \end{gather}$$

(3.22)$$\begin{gather}-\varepsilon\hat{\tau}_{zx}^w\frac{\partial\hat{b}}{\partial\hat{x}}+\varepsilon \hat{\tau}_{zz}^w-\left(\phi^w\hat{p}^{w\star}_b+(\boldsymbol{n}^b\boldsymbol{\cdot} \hat{\boldsymbol{\sigma}}^w\boldsymbol{n}^b)\right)=\varepsilon C^w\hat{\bar{w}}^w|{\hat{\bar{\boldsymbol{u}}}^w}\|\boldsymbol{\nabla} F^b|,\quad \hat{z}=\hat{b}(\hat{x}), \end{gather}$$

where $|\hat {\bar {\boldsymbol {u}}}^w|=\sqrt {(\hat {\bar {u}}^w)^2 + \varepsilon ^2(\hat {\bar {w}}^w)^2}$.

4.1. Depth-averaged mass balance equations

The mass balance equations (3.2) and (3.6) can be integrated through the flow depth using Leibniz’s rule (Abramowitz & Stegun Reference Abramowitz and Stegun1970) to exchange the order of differentiation and integration, and then simplified by using the kinematic conditions (3.9)–(3.12). In the oversaturated regime, the standard argument (Savage & Hutter Reference Savage and Hutter1989; Gray et al. Reference Gray, Wieland and Hutter1999) is complicated by the jump in the water volume fraction at the internal free surface of the grains. To overcome this, the integration of (3.6) is divided into two parts, which implies that the depth-integrated water mass balance is

(4.1)\begin{equation} \frac{\partial}{\partial\hat{t}}(\hat{h}^w\bar{\phi}^w)+ \frac{\partial}{\partial\hat{x}}(\hat{h}^w\overline{\phi^w{\hat{u}}^w}) +\left[\!\!\left[\phi^{w}\left(\frac{\partial\hat{s}^g}{\partial\hat{t}}+\hat{u}^{w} \frac{\partial\hat{s}^g}{\partial\hat{x}}-\hat{w}^{w}\right)\right]\!\!\right]=0, \end{equation}

where the jump bracket is evaluated at $\hat z=\hat s^g$, the water thickness $\hat {h}^w=\hat {s}^w-\hat {b}$, and the depth-averaged water concentration and depth-averaged water flux are

(4.2a,b)\begin{equation} \bar{\phi}^w=\frac{1}{\hat{h}^w}\int_{\hat{b}}^{\hat{s}^w}\phi^wd\hat{z},\quad \overline{\phi^w\hat{u}^w}=\frac{1}{\hat{h}^w}\int_{\hat{b}}^{\hat{s}^w}\phi^w{\hat{u}}^wd\hat{z}, \end{equation}

respectively. The jump bracketed term in (4.1) is zero by the water mass jump condition (3.14), so the depth-averaged water mass balance reduces to standard form, despite the jump in concentration. A similar argument carries through in the same way for the depth-averaged grain mass balance in the undersaturated regime. As a result, the depth-averaged grain and water mass balances are

(4.3)$$\begin{gather} \frac{\partial}{\partial\hat{t}}(\hat{h}^g\bar{\phi}^g)+ \frac{\partial}{\partial\hat{x}}(\hat{h}^g\overline{\phi^g{\hat{u}}^g})=0, \end{gather}$$

(4.4)$$\begin{gather}\frac{\partial}{\partial\hat{t}}(\hat{h}^w\bar{\phi}^w)+ \frac{\partial}{\partial\hat{x}}(\hat{h}^w\overline{\phi^w{\hat{u}}^w})=0, \end{gather}$$

where the granular flow thickness $\hat h^g=\hat s^g-\hat b$ and

(4.5a,b)\begin{equation} \bar{\phi}^g=\frac{1}{\hat{h}^g}\int_{\hat{b}}^{\hat{s}^g}\phi^gd\hat{z},\quad \overline{\phi^g\hat{u}^g}=\frac{1}{\hat{h}^g}\int_{\hat{b}}^{\hat{s}^g}\phi^g{\hat{u}}^gd\hat{z}. \end{equation}

To further simplify (4.3) and (4.4), the vertical distributions of the solid and water volume fractions must be prescribed. The simplest approach is to assume that $\phi ^g$ is constant, regardless of whether the flow is undersaturated or oversaturated, i.e.

(4.6)\begin{equation} \phi^g=\phi^c. \end{equation}

This is assumed in early debris flow models (Iverson Reference Iverson1997; Iverson & Denlinger Reference Iverson and Denlinger2001) and appears to be a good approximation in high solids volume fraction granular-dominated flows (Maurin et al. Reference Maurin, Chauchat and Frey2016; Chassagne et al. Reference Chassagne, Maurin, Chauchat, Gray and Frey2020). However, when the fluid can suspend grains, the solids volume fraction increases with depth (Egashira, Itoh & Takeuchi Reference Egashira, Itoh and Takeuchi2001; Maurin et al. Reference Maurin, Chauchat and Frey2016) and this will not be a good approximation. The fact that $\phi ^c$ is constant and uniform throughout the debris flow precludes pore pressure effects (Iverson & George Reference Iverson and George2014). There is therefore potential to improve the model (at this point) in future. The exclusion of pore pressure effects does, however, allow us to show that they are not necessarily needed to drive the forward motion of grains from an oversaturated tail towards the flow front, which is the traditional view in the debris flow community.

As sketched in figure 2(b,c) the volume fraction of water is assumed to have a vertical distribution that is dependent on whether the flow is undersaturated or oversaturated. In the undersaturated regime,

(4.7)\begin{equation} \phi^w=\left\{\begin{array}{ll} 0, & \hat z\in[\hat s^w,\hat s^g],\\ 1-\phi^c, & \hat z\in[\hat b,\hat s^w], \end{array}\right. \end{equation}

while in the oversaturated regime,

(4.8)\begin{equation} \phi^w=\left\{\begin{array}{ll} 1, & \hat z\in[\hat s^g,\hat s^w],\\ 1-\phi^c, & \hat z\in[\hat b,\hat s^g]. \end{array}\right. \end{equation}

It is useful to define volume-fraction-weighted depth-averaged velocities

(4.9a,b)\begin{equation} \hat{\bar{u}}^g=\frac{\displaystyle\int_{\hat{b}}^{\hat{s}^g}\phi^g\hat{u}^g\,\text{d}\hat{z}}{\displaystyle \int_{\hat{b}}^{\hat{s}^g}\phi^g\,\text{d}\hat{z}}=\frac{\overline{\phi^g{\hat{u}}^g}}{\bar{\phi}^g},\quad \hat{\bar{u}}^w=\frac{\displaystyle\int_{\hat{b}}^{\hat{s}^w}\phi^w\hat{u}^w\,\text{d}\hat{z}}{\displaystyle \int_{\hat{b}}^{\hat{s}^w}\phi^w\,\text{d}\hat{z}}=\frac{\overline{\phi^w{\hat{u}}^w}}{\bar{\phi}^w}. \end{equation}

Since the solids volume fraction is constant, (4.9a) reduces to the standard definition (Savage & Hutter Reference Savage and Hutter1989; Gray et al. Reference Gray, Wieland and Hutter1999; Pitman & Le Reference Pitman and Le2005; Gray & Edwards Reference Gray and Edwards2014). This is also true for the water in the undersaturated regime, but in the oversaturated regime, the internal concentration discontinuity (4.8) implies that the definition of the depth-averaged water velocity (4.9b) is different. Using the definitions (4.6)–(4.9a,b), it follows that the depth-averaged mass balance of the grains in both the undersaturated and oversaturated regimes is

(4.10)\begin{equation} \frac{\partial}{\partial\hat{t}}(\hat{h}^g\phi^c) + \frac{\partial}{\partial\hat{x}} (\hat{h}^g\phi^c\hat{\bar{u}}^g)=0, \end{equation}

while the depth-averaged mass balances for the water in the undersaturated and oversaturated regimes are

(4.11)$$\begin{gather} \frac{\partial}{\partial\hat{t}}(\hat{h}^w(1-\phi^c))+ \frac{\partial}{\partial\hat{x}}(\hat{h}^w(1-\phi^c)\hat{\bar{u}}^w)=0, \end{gather}$$

(4.12)$$\begin{gather}\frac{\partial}{\partial\hat{t}}(\hat{h}^w(1-\phi^c\hat{h}^g/\hat{h}^w))+ \frac{\partial}{\partial\hat{x}}(\hat{h}^w(1-\phi^c\hat{h}^g/\hat{h}^w)\hat{\bar{u}}^w)=0, \end{gather}$$

respectively. The factors $\phi ^c$ and $1-\phi ^c$ in (4.10) and (4.11) are constants that could be cancelled out. However, these factors are retained here to develop a theoretical framework that can switch seamlessly between regimes.

4.2. Normal components of the momentum equations

The normal components of the momentum balances play a crucial role in deriving normal stresses in depth-averaged models (Savage & Hutter Reference Savage and Hutter1989; Iverson & Denlinger Reference Iverson and Denlinger2001; Gray & Edwards Reference Gray and Edwards2014). To leading order in $\varepsilon$, the water normal momentum balance (3.8) implies that the fluid pressure is hydrostatic

(4.13)\begin{equation} \frac{\partial\hat{p}^{w\star}}{\partial \hat{z}}={-}\cos\zeta, \end{equation}

in both the mixed (ii) and pure water (iii) regions shown in figure 2(b,c). Excess pore pressure, which is important in large-scale debris flows containing fine particles (Iverson et al. Reference Iverson, Logan, Lahusen and Berti2010), is therefore neglected in this model. In the dry region (i), the water pressure and the Darcy drag vanish, and to leading order, the dry grain normal momentum balance (3.4) reduces to

(4.14)\begin{equation} \frac{\partial\hat{\sigma}_{zz}^e}{\partial \hat{z}} ={-}\phi^c\cos\zeta, \end{equation}

while in the mixed region (ii), the saturated grain normal stress gradient balances gravity and buoyancy

(4.15)\begin{equation} \frac{\partial\hat{\sigma}_{zz}^e}{\partial \hat{z}} ={-}\phi^c\cos\zeta-\gamma\phi^c\displaystyle\frac{\partial\hat{p}^{w\star}}{\partial\hat{z}}. \end{equation}

Equations (4.13)–(4.15), together with their respective boundary conditions, enable expressions for the pore fluid pressure and the grain normal stress to be derived, as presented below.

4.3. Undersaturated pore fluid pressure and grain normal stress

To leading order, the surface condition (3.18) reduces to

(4.16)\begin{equation} \hat{p}^{w\star}(\hat{s}^w)=0. \end{equation}

The integration of (4.13) subject to (4.16) implies that the water pressure is hydrostatic

(4.17)\begin{equation} \hat{p}^{w\star}=(\hat{s}^w-\hat{z})\cos\zeta,\quad \hat{z}\in [\hat b,\hat{s}^w]. \end{equation}

It follows that the pore fluid pressure on the base is

(4.18)\begin{equation} \hat{p}^{w\star}_b=\hat{h}^w\cos\zeta, \end{equation}

and that the depth-averaged water pressure is

(4.19)\begin{equation} \hat{\bar{p}}^{w\star}=\frac{1}{\hat{h}^w}\int_{\hat{b}}^{\hat{s}^w}\hat{p}^{w\star}\,\text{d}\hat{z}= \frac{1}{2}\hat{h}^w\cos\zeta. \end{equation}

To leading order, the grain surface condition (3.16) implies that the surface stress $\hat {\sigma }_{zz}^e$ vanishes. It follows that integrating (4.14) implies that the perpendicular grain normal stress in the dry region (i) is

(4.20)\begin{equation} \hat{\sigma}_{zz}^e=(\hat{s}^g-\hat{z})\phi^c\cos\zeta,\quad \hat z\in[\hat{s}^w,\hat{s}^g], \end{equation}

and hence that the perpendicular grain normal stress in the upper side of the water surface is $\hat {\sigma }_{zz}^{e+}=(\hat {s}^g-\hat {s}^w)\phi ^c\cos \zeta$. Appendix A shows that provided the granular velocity profile through the depth of the flow is continuous, then there is no jump in the normal effective stress, and hence

(4.21)\begin{equation} \hat{\sigma}_{zz}^{e-}=\hat{\sigma}_{zz}^{e+}=(\hat{s}^g-\hat{s}^w)\phi^c\cos\zeta,\quad \hat z=\hat s^w. \end{equation}

In the saturated regime (ii), substitution of (4.13) into (4.15) implies that

(4.22)\begin{equation} \frac{\partial\hat{\sigma}_{zz}^e}{\partial \hat{z}}={-}\phi^c(1-\gamma)\cos\zeta,\quad \hat{z}\in [\hat b,\hat{s}^w]. \end{equation}

Integration of (4.22) subject to the continuity condition (4.21) implies that the perpendicular grain normal stress in the saturated regime (ii) is

(4.23)\begin{equation} \hat{\sigma}_{zz}^e=(\hat{s}^g-\hat{z})\phi^c\cos\zeta-(\hat{s}^w-\hat{z})\gamma\phi^c\cos\zeta,\quad \hat{z}\in[\hat b,\hat{s}^w], \end{equation}

and hence on the base,

(4.24)\begin{equation} \hat{\sigma}_{zz}^e(b) = (\hat{h}^g-\gamma\hat{h}^w)\phi^c\cos\zeta,\quad \hat z=\hat b. \end{equation}

The depth-averaged perpendicular grain normal stress is

(4.25)\begin{equation} \hat{\bar{\sigma}}_{zz}^e=\frac{1}{\hat{h}^g}\left(\int_{\hat{b}}^{\hat{s}^w}\hat{\sigma}^e_{zz}d\hat{z}+ \int_{\hat{s}^w}^{\hat{s}^g}\hat{\sigma}^e_{zz}\,\text{d}\hat{z}\right)= \frac{1}{2}\hat{h}^g\phi^c\cos\zeta-\frac{1}{2}\gamma(\hat{h}^w/\hat{h}^g)^2\hat{h}^g\phi^c\cos\zeta, \end{equation}

which, provided the earth pressure coefficient is equal to unity (Savage & Hutter Reference Savage and Hutter1989; Gray et al. Reference Gray, Wieland and Hutter1999; Pouliquen Reference Pouliquen1999a), implies that the depth-averaged downslope grain normal stress is

(4.26)\begin{equation} \hat{\bar{\sigma}}_{xx}^e=\tfrac{1}{2}\hat{h}^g\phi^c\cos\zeta-\tfrac{1}{2} \gamma(\hat{h}^w/\hat{h}^g)^2\hat{h}^g\phi^c\cos\zeta. \end{equation}

It follows that to leading order, (3.20) implies that the basal normal stress

(4.27)\begin{equation} \boldsymbol{n}^b\boldsymbol{\cdot}\hat{\boldsymbol{\sigma}}^e\boldsymbol{n}^b=\hat{\sigma}^e_{zz}(b) = \hat{h}^g(1-\gamma\hat{h}^w/\hat{h}^g)\phi^c\cos\zeta,\quad \hat z=\hat b. \end{equation}

4.4. Oversaturated pore fluid pressure and grain normal stress

Since (4.13) is independent of $\phi ^w$, it can be integrated in the same way as the undersaturated case to show that the pressure is hydrostatic

(4.28)\begin{equation} \hat{p}^{w\star}=(\hat{s}^w-\hat{z})\cos\zeta. \end{equation}

Similarly, the basal water pressure and the depth-averaged water pressure are the same as the undersaturated case

(4.29a,b)\begin{equation} \hat{p}^{w\star}_b=\hat{h}^w\cos\zeta,\quad \hat{\bar{p}}^{w\star}=\tfrac{1}{2}\hat{h}^w\cos\zeta, \end{equation}

and the pore fluid pressure at the grain surface is

(4.30)\begin{equation} \hat{p}^{w\star}(\hat{s}^g)=(\hat{h}^w-\hat{h}^g)\cos\zeta,\quad \hat{z}=\hat{s}^g. \end{equation}

To leading order, the traction-free condition (3.16) reduces to $\hat {\sigma }_{zz}^e(\hat {s}^g)=0$. Substituting the gradient of the pore fluid pressure (4.13) into (4.15) and then integrating implies the perpendicular grain normal stress is

(4.31)\begin{equation} \hat{\sigma}_{zz}^e=(1-\gamma)(\hat{s}^g-\hat{z})\phi^c\cos\zeta,\quad \hat{z}\in[\hat b,\hat{s}^g], \end{equation}

and hence that the perpendicular basal normal stress is

(4.32)\begin{equation} \hat{\sigma}_{zz}^e(\hat{b}) = (1-\gamma)\hat{h}^g\phi^c\cos\zeta\quad \hat z=\hat b. \end{equation}

The depth-averaged perpendicular normal stress is then

(4.33)\begin{equation} \hat{\bar{\sigma}}_{zz}^e=\frac{1}{\hat{h}^g}\int_{\hat{b}}^{\hat{s}^g}\hat{\sigma}_{zz}^ed\hat{z}= \frac{1}{2}(1-\gamma)\hat{h}^g\phi^c\cos\zeta. \end{equation}

Given that the earth pressure coefficient is unity, the depth-averaged downslope normal stress is therefore equal to its counterpart in the vertical direction, i.e.

(4.34)\begin{equation} \hat{\bar{\sigma}}_{xx}^e=\tfrac{1}{2}(1-\gamma)\hat{h}^g\phi^c\cos\zeta. \end{equation}

To leading order, (3.20) shows the basal normal stress

(4.35)\begin{equation} \boldsymbol{n}^b\boldsymbol{\cdot}\hat{\boldsymbol{\sigma}}^e\boldsymbol{n}^b = \hat{\sigma}_{zz}^e(\hat{b})=(1-\gamma)\hat{h}^g\phi^c\cos\zeta. \end{equation}

4.5. Undersaturated depth-averaged downslope momentum balances

The downslope grain momentum balance (3.3) is now integrated from the bed $\hat {z}=\hat {b}$ to the grain surface $\hat {z}=\hat {s}^g$. The integration is divided into two parts, because the buoyancy and Darcy drag are active only below the water surface $\hat {z}=\hat {s}^w$. Leibniz’s integration rule is used to swap the order of differentiation and integration, and the surface and basal kinematic conditions, (3.9) and (3.10), and downslope surface traction (3.15) are used to simplify the resulting equations. The calculations are lengthy and are therefore shown in Appendix B. The resulting depth-averaged momentum balance for the grains is

(4.36)\begin{align} & \varepsilon\left(\frac{\partial}{\partial\hat{t}}(\hat{h}^g\overline{\phi^g\hat{u}^g}) +\frac{\partial}{\partial\hat{x}}(\hat{h}^g\overline{\phi^g(\hat{u}^g)^2}) +\frac{\partial}{\partial\hat{x}}(\hat{h}^g\hat{\bar{\sigma}}_{xx}^e) +\gamma\phi^c\frac{\partial}{\partial\hat{x}}(\hat{h}^w\hat{\bar{p}}^w)\right) \nonumber\\ &\quad =\hat{h}^g\phi^c\sin\zeta - \left(\varepsilon\hat{\sigma}_{xx}^e(\hat b) \frac{\partial\hat{b}}{\partial\hat{x}}-\mu\hat{\sigma}_{xz}^{e}(\hat b)\right) -\gamma\hat{C}^d\int_{\hat{b}}^{\hat{s}^w}(\hat{u}^g-\hat{u}^w)\,\text{d}\hat{z} -\varepsilon\gamma\phi^c\hat{p}^{w\star}_b\frac{\partial\hat{b}}{\partial\hat{x}}. \end{align}

Substituting (4.27) into (3.19) implies that the basal downslope traction

(4.37)\begin{align} \varepsilon\hat{\sigma}_{xx}^e(\hat b)\frac{\partial\hat{b}}{\partial\hat{x}}-\mu\hat{\sigma}_{xz}^{e}(\hat b) &= \frac{\hat{u}^g_b}{|\hat{\boldsymbol{u}}^g_b|}(1-\gamma\hat{h}^w/\hat{h}^g)\mu^b\hat{h}^g\phi^c\cos\zeta \nonumber\\ &\quad +\varepsilon(1-\gamma{\hat{h}^w}/\hat{h}^g)\hat{h}^g\phi^c\cos\zeta \frac{\partial\hat{b}}{\partial\hat{x}}+O(\varepsilon^2), \end{align}

where the factor $\hat {u}^g_b/|\hat {\boldsymbol {u}}^g_b|$ will be approximated as $\hat {\bar {u}}^g/|\hat {\bar {u}}^g|$. This approximation is employed in virtually all depth-averaged granular flow models (see e.g. Savage & Hutter Reference Savage and Hutter1989; Pitman & Le Reference Pitman and Le2005; Pudasaini Reference Pudasaini2012; Gray & Edwards Reference Gray and Edwards2014; Iverson & George Reference Iverson and George2014). To simplify the momentum transport terms, it is assumed that

(4.38)\begin{equation} \overline{\phi^g(\hat{u}^g)^2}=\frac{1}{\hat{h}^g}\int_{\hat{b}}^{\hat{s}^g} \phi^g(\hat{u}^g)^2\,\text{d}\hat{z}=\chi^g\bar\phi^g(\hat{\bar{u}}^g)^2, \end{equation}

where the shape factor $\chi ^g$ for specific velocity profiles will be considered in § 4.8. It is also useful to define the streamfunction for the grains

(4.39)\begin{equation} \hat\psi^g(\hat z)=\int_{\hat{b}}^{\hat{z}} \hat{u}^g(\hat z')\,\text{d}\hat z', \end{equation}

which allows the Darcy drag exerted on the grain phase to be expressed as

(4.40)\begin{equation} \gamma\hat{C}^d\int_{\hat{b}}^{\hat{s}^w}(\hat{u}^g-\hat{u}^w)\,\text{d}\hat{z}= \gamma\hat{C}^d(\hat{\psi}^g(\hat s^w)-\hat{h}^w\hat{\bar{u}}^w). \end{equation}

The gravitational force, the basal friction and the Darcy drag are all order unity quantities, but their sum is an order $\varepsilon$ quantity. To express this theoretically, the leading order balance $\hat {\mathcal {L}}^g$ is defined as

(4.41)\begin{equation} \varepsilon\hat{\mathcal{L}}^g = \underbrace{\hat{h}^g\phi^c\sin\zeta}_{Gravity} -\underbrace{\mu^b\frac{\hat{\bar{u}}^g}{|\hat{\bar{u}}^g|} \left(1-\gamma\frac{\hat{h}^w}{\hat{h}^g}\right)\hat{h}^g\phi^c\cos\zeta}_{Basal\ friction} -\underbrace{\gamma\hat{C}^d(\hat{\psi}^g(\hat s^w)-\hat{h}^w\hat{\bar{u}}^w)}_{Darcy\ drag}. \end{equation}

When (4.6), (4.9a), (4.37), (4.38) and (4.41) are substituted into (4.36), together with the basal pressure (4.18), the depth-averaged pore fluid pressure (4.19) and the depth-averaged normal stress (4.26), the leading-order depth-averaged grain momentum balance in the undersaturated regime becomes

(4.42)\begin{equation} \frac{\partial}{\partial\hat{t}}(\hat{h}^g\phi^c\hat{\bar{u}}^g) + \frac{\partial }{\partial \hat{x}}\left(\chi^g \hat{h}^g\phi^c(\hat{\bar{u}}^g)^2+ {\frac{1}{2}(\hat{h}^g)^2\phi^c\cos\zeta}-{\frac{1}{2}\gamma (\hat{h}^w)^2\phi^c\cos\zeta}\right)=\hat{S}^g, \end{equation}

where the source terms

(4.43)\begin{equation} \hat{S}^g=\hat{\mathcal{L}}^g -\underbrace{\gamma\hat{h}^w\phi^c\cos\zeta\frac{\partial\hat{h}^w}{\partial\hat{x}}}_{{Buoyancy}} -\underbrace{\hat{h}^g\phi^c\cos\zeta\frac{\partial\hat{b}}{\partial\hat{x}}}_{{Topographic\ gradient}}, \end{equation}

combine the leading order balance (4.41) with the lower order buoyancy and topographic gradient terms. Note that the gradient of $(1/2)\gamma (\hat {h}^w)^2\phi ^c\cos \zeta$ in (4.42) can be cancelled with the buoyancy term in (4.43). However, the equations are left in their present form to unify all the regimes in § 4.7.

Integrating the downslope water momentum balance (3.7) from the bed $\hat {z}=\hat {b}$ to the water surface $\hat {z}=\hat {s}^w$ is more straightforward, because the integration does not pass through a singular surface. Using Leibniz’s integration rule to exchange the order of differentiation and integration and using the kinematic conditions (3.11) and (3.12), the surface traction (3.17) and the Darcy drag (4.40), the depth-averaged downslope water momentum balance is

(4.44)\begin{align} & \varepsilon\left(\frac{\partial}{\partial\hat{t}}(\hat{h}^w(1-\phi^c)\hat{\bar{u}}^w) +\frac{\partial}{\partial\hat{x}}(\chi^w\hat{h}^w(1-\phi^c)(\hat{\bar{u}}^w)^2) +(1-\phi^c)\frac{\partial}{\partial\hat{x}}(\hat{h}^w\hat{\bar{p}}^{w\star})\right) \nonumber\\ &\quad -\varepsilon^2\frac{\partial}{\partial\hat{x}}(\hat{h}^w\hat{\bar{\tau}}^w_{xx}) =\hat{h}^w(1-\phi^c)\sin\zeta-\hat{C}^d(\hat{h}^w\hat{\bar{u}}^w-\hat{\psi}^g(\hat s^w)) -\varepsilon(1-\phi^c)\hat{p}^{w\star}_b\frac{\partial\hat{b}}{\partial\hat{x}} \nonumber\\ &\quad +\left(\varepsilon^2\hat{\tau}_{xx}^w(\hat{b}) \frac{\partial\hat{b}}{\partial \hat{x}}-\hat{\tau}_{xz}^w(\hat{b})\right), \end{align}

where the water shape factor $\chi ^w$ is defined as

(4.45)\begin{equation} \overline{\phi^w(\hat{u}^w)^2}=\frac{1}{\hat{h}^w} \int_{\hat{b}}^{\hat{s}^w}\phi^w(\hat{u}^w)^2\,\text{d}\hat{z}=\chi^w\bar\phi^w(\hat{\bar{u}}^w)^2. \end{equation}

The depth-averaged in-plane viscous stress $\hat {\bar {\tau }}^w_{xx}$ is smaller than the other forces, and is usually neglected (Savage & Hutter Reference Savage and Hutter1989; Pitman & Le Reference Pitman and Le2005). However, it plays a crucial role in some subtle situations, including (i) obtaining the correct cutoff frequency of roll waves (Forterre Reference Forterre2006; Gray & Edwards Reference Gray and Edwards2014), (ii) generating cross-stream velocity profiles (Iverson & Denlinger Reference Iverson and Denlinger2001; Baker, Johnson & Gray Reference Baker, Johnson and Gray2016b; Meng & Wang Reference Meng and Wang2018), (iii) regularising ill-posedness (Baker et al. Reference Baker, Johnson and Gray2016a) and (iv) forming critical phenomena like levees (Rocha, Johnson & Gray Reference Rocha, Johnson and Gray2019). An expression for the depth-averaged in-plane deviatoric water stresses in the undersaturated regime is derived in Appendix C. A similar depth-averaged in-plane viscous-like term, analogous to that derived by Gray & Edwards (Reference Gray and Edwards2014), could also be added to (4.43). However, both these terms are small, and are not important for reproducing the experimental observations of Davies (Reference Davies1988, Reference Davies1990) in § 5.

To leading order, (3.22) reduces to

(4.46)\begin{equation} \boldsymbol{n}^b\boldsymbol{\cdot}\hat{\boldsymbol{\sigma}}^w\boldsymbol{n}^b={-}\phi^w\hat{p}^{w\star}_b, \end{equation}

and together with (3.21), it shows

(4.47)\begin{equation} -\varepsilon^2\hat{\tau}_{xx}^w(\hat b)\frac{\partial\hat{b}}{\partial\hat{x}}+ \hat{\tau}_{xz}^w(\hat{b})=C^w\hat{\bar{u}}^w|{\hat{\bar{u}}^w}|+O(\varepsilon^2). \end{equation}

In an analogous way to the grain phase, the sum of the gravitational force, the Chézy drag and the Darcy drag is small, although they are all individually order unity quantities. To express this mathematically, the leading order balance $\hat {\mathcal {L}}^w$ is defined as

(4.48)\begin{equation} \varepsilon\hat{\mathcal{L}}^w= \underbrace{\hat{h}^w(1-\phi^c)\sin\zeta}_{Gravity} - \underbrace{C^w\hat{\bar{u}}^w\vert\hat{\bar{u}}^w\vert}_{Ch{\'e}zy\ drag} - \underbrace{\hat{C}^d(\hat{h}^w\hat{\bar{u}}^w-\hat{\psi}^g(\hat s^w)).}_{Darcy\ drag} \end{equation}

When (4.45) and (4.48) are substituted into (4.44), together with the basal pressure (4.18), the depth-averaged pore fluid pressure (4.19) and the depth-averaged shear stress (C3), the leading-order water momentum balance in the undersaturated regime becomes

(4.49)\begin{equation} \frac{\partial}{\partial\hat{t}}(\hat{h}^w(1-\phi^c)\hat{\bar{u}}^w) + \frac{\partial}{\partial\hat{x}} \left(\chi^w \hat{h}^w(1-\phi^c)(\hat{\bar{u}}^w)^2+{\frac{1}{2}(\hat{h}^w)^2\cos\zeta}\right)=\hat{S}^w, \end{equation}

where the source term

(4.50)\begin{equation} \hat{S}^w=\hat{\mathcal{L}}^w +\underbrace{\phi^c\hat{h}^w\cos\zeta\frac{\partial\hat{h}^w}{\partial\hat{x}}}_{Buoyancy} -\underbrace{\hat{h}^w(1-\phi^c)\cos\zeta\frac{\partial\hat{b}}{\partial\hat{x}}}_{Topographic\ gradient} +\underbrace{\frac{2\varepsilon}{Re}\frac{\partial}{\partial\hat{x}}\left(\hat{h}^w(1-\phi^c) \frac{\partial\hat{\bar{u}}^w}{\partial\hat{x}}\right),}_{Viscous\ term} \end{equation}

combines the leading order balance (4.48) with the lower order buoyancy, topography gradient and viscous terms.

4.6. Oversaturated depth-averaged downslope momentum balances

It is also straightforward to integrate the downslope granular momentum balance (3.3) from $\hat {z}=\hat {b}$ to the grain surface $\hat {z}=\hat {s}^g$, since there is no singular surface. By using Leibniz’s integration rule to exchange the order of differentiation and integration, together with the surface and basal kinematic conditions (3.9) and (3.10), the depth-averaged downslope granular momentum balance in the oversaturated region is

(4.51)\begin{align} & \varepsilon\left(\frac{\partial}{\partial\hat{t}}(\hat{h}^g\phi^c\hat{\bar{u}}^g)+ \frac{\partial}{\partial\hat{x}}(\chi^g\hat{h}^g\phi^c (\hat{\bar{u}}^g)^2)+\frac{\partial}{\partial\hat{x}}(\hat{h}^g\hat{\bar{\sigma}}_{xx}^e)+\gamma\phi^c \frac{\partial}{\partial\hat{x}}\left(\int_{\hat{b}}^{\hat{s}^g}\hat{p}^{w\star}d\hat{z}\right)\right) \nonumber\\ &\quad =\hat{h}^g\phi^c\sin\zeta - \left(\varepsilon\hat{\sigma}_{xx}^e(\hat b)\frac{\partial\hat{b}}{\partial\hat{x}}-\mu\hat{\sigma}_{xz}^{e}(\hat b)\right) -\gamma\hat{C}^d\int_{\hat{b}}^{\hat{s}^g}(\hat{u}^g-\hat{u}^w)\,\text{d}\hat{z} \nonumber\\ &\qquad +\varepsilon\gamma\phi^c\left(\hat{p}^{w\star}\frac{\partial\hat{z}}{\partial\hat{x}}\right)^{\hat{s}^g}_{\hat{b}}, \end{align}

where the basal granular traction is determined by the downslope component of the basal friction condition (3.19) and the basal normal stress (4.35), i.e.

(4.52)\begin{equation} \varepsilon\hat{\sigma}_{xx}^e\frac{\partial\hat{b}}{\partial\hat{x}}-\mu\hat{\sigma}_{xz}^e =\mu^b \frac{\hat{\bar{u}}^g}{|\hat{\bar{u}}^g|}(1-\gamma)\hat{h}^g\phi^c\cos\zeta +\varepsilon(1-\gamma)\hat{h}^g\phi^c\cos\zeta\frac{\partial\hat{b}}{\partial\hat{x}}+O(\varepsilon^2). \end{equation}

This time, it is useful to define the streamfunction for the water

(4.53)\begin{equation} \hat\psi^w(\hat z)=\int_{\hat{b}}^{\hat{z}} \hat{u}^w(\hat z')\,\text{d}\hat z', \end{equation}

which allows the Darcy drag to be expressed as

(4.54)\begin{equation} \gamma\hat{C}^d\int_{\hat{b}}^{\hat{s}^g}(\hat{u}^g-\hat{u}^w)\,\text{d}\hat{z}= \gamma\hat{C}^d(\hat{h}^g\hat{\bar{u}}^g-\hat{\psi}^w(\hat s^g)). \end{equation}

The difference of the gravitational force, the basal friction and the Darcy drag is of order $\varepsilon$. This leading order balance is expressed mathematically by defining the function

(4.55)\begin{equation} \varepsilon\hat{\mathcal{L}}^g = \underbrace{\hat{h}^g\phi^c\sin\zeta}_{Gravity} -\underbrace{\mu^b\frac{\hat{\bar{u}}^g}{|\hat{\bar{u}}^g|} (1-\gamma)\hat{h}^g\phi^c\cos\zeta}_{Basal\ friction} -\underbrace{\gamma\hat{C}^d(\hat{h}^g\hat{\bar{u}}^g-\hat{\psi}^w(\hat s^g))}_{Darcy\ drag}. \end{equation}

Substituting the hydrostatic water pressure (4.28) into the relevant terms in (4.51) implies that

(4.56)\begin{equation} -\frac{\partial}{\partial\hat{x}}\left(\int_{\hat{b}}^{\hat{s}^g}\hat{p}^{w\star}\,\text{d}\hat{z}\right) +\left(\hat{p}^{w\star}\frac{\partial\hat{z}}{\partial\hat{x}}\right)_{\hat{b}}^{\hat{s}^g} ={-}\hat{h}^g\cos\zeta\frac{\partial\hat{h}^w}{\partial\hat{x}}- \hat{h}^g\cos\zeta\frac{\partial\hat{b}}{\partial\hat{x}}. \end{equation}

Finally, substituting (4.54)–(4.56) together with the depth-averaged downslope normal stress (4.34) into (4.51) implies that the depth-averaged granular momentum balance is

(4.57)\begin{equation} \frac{\partial}{\partial\hat{t}}(\hat{h}^g\phi^c\hat{\bar{u}}^g) + \frac{\partial }{\partial \hat{x}}\left(\chi^g \hat{h}^g\phi^c(\hat{\bar{u}}^g)^2+ \frac{1}{2}(\hat{h}^g)^2(1-\gamma)\phi^c\cos\zeta\right)=\hat{S}^g, \end{equation}

where the source term $\hat {S}^g$ on the right-hand side is

(4.58)\begin{equation} \hat{S}^g=\hat{\mathcal{L}}^g-\underbrace{\gamma\hat{h}^g\phi^c\cos\zeta \frac{\partial\hat{h}^w}{\partial\hat{x}}}_{Buoyancy} -\underbrace{\hat{h}^g\phi^c\cos\zeta\frac{\partial\hat{b}}{\partial\hat{x}}.}_{Topographic\ gradient} \end{equation}

To leading order, the normal component of the Chézy formula (3.22) reduces to

(4.59)\begin{equation} \boldsymbol{n}^b\boldsymbol{\cdot}\hat{\boldsymbol{\sigma}}^w\boldsymbol{n}^b={-}\phi^w\hat{p}^{w\star}_b, \end{equation}

which implies the downslope component of the Chézy formula (3.21) becomes

(4.60)\begin{equation} -\varepsilon^2\hat{\tau}_{xx}^w(\hat b)\frac{\partial\hat{b}}{\partial\hat{x}}+\hat{\tau}_{xz}^w(\hat{b}) =C^w\hat{\bar{u}}^w|{\hat{\bar{u}}^w}|+O(\varepsilon^2). \end{equation}

In the oversaturated region, the water volume fraction is discontinuous across the grain free surface, so the integral of the downslope component of the water momentum balance (3.7) through the water depth is complicated. The detailed derivation can be found in Appendix D. It follows from this calculation that the depth-integrated water momentum balance in the oversaturated regime is

(4.61)\begin{align} & \varepsilon\left(\frac{\partial}{\partial\hat{t}}((\hat{h}^w-\phi^c\hat{h}^g)\hat{\bar{u}}^w) +\frac{\partial}{\partial\hat{x}}(\chi^w(\hat{h}^w-\phi^c\hat{h}^g)(\hat{\bar{u}}^w)^2) +\frac{\partial}{\partial\hat{x}}(\hat{h}^w\hat{\bar{p}}^{w\star}) -\varepsilon\frac{\partial}{\partial\hat{x}}(\hat{h}^w\hat{\bar{\tau}}^w_{xx})\right) \nonumber\\ &\quad =(\hat{h}^w-\phi^c\hat{h}^g)\sin\zeta -C^w\hat{\bar{u}}^w\vert\hat{\bar{u}}^w\vert -\hat{C}^d(\hat{\psi}^w(\hat s^g)-\hat{h}^g\hat{\bar{u}}^g) +\varepsilon\left(\hat{p}^{w\star}\frac{\partial\hat{z}}{\partial\hat{x}}\right)_{\hat{b}}^{\hat{s}^w} \nonumber\\ &\qquad +\varepsilon\phi^c\frac{\partial}{\partial\hat{x}} \left(\int_{\hat{b}}^{\hat{s}^g}\hat{p}^{w\star}\,\text{d}\hat{z}\right) -\varepsilon\phi^c\left(\hat{p}^{w\star}\frac{\partial\hat{z}}{\partial\hat{x}}\right)_{\hat{b}}^{\hat{s}^g}. \end{align}

This time, the leading order balance lies between the gravitational force, the Chézy drag and the Darcy drag. The sum of these terms is of order $\varepsilon$, which is expressed mathematically by defining the function

(4.62)\begin{equation} \varepsilon\hat{\mathcal{L}}^w= \underbrace{(\hat{h}^w-\phi^c\hat{h}^g)\sin\zeta}_{Gravity} - \underbrace{C^w\hat{\bar{u}}^w\vert\hat{\bar{u}}^w\vert}_{Ch{\'e}zy\ drag} - \underbrace{\hat{C}^d(\hat{\psi}^w(\hat s^g)-\hat{h}^g\hat{\bar{u}}^g).}_{Darcy\ drag} \end{equation}

Using (4.16), (4.29a,b), (4.56), (4.62) and the depth-averaged shear stress (C4), the leading-order depth-averaged oversaturated water momentum balance becomes

(4.63)\begin{equation} \frac{\partial}{\partial\hat{t}}((\hat{h}^w-\phi^c\hat{h}^g)\hat{\bar{u}}^w) + \frac{\partial }{\partial\hat{x}}\left(\chi^w(\hat{h}^w-\phi^c\hat{h}^g) (\hat{\bar{u}}^w)^2+{\frac{1}{2}(\hat{h}^w)^2\cos\zeta}\right)=\hat{S}^w, \end{equation}

where the source term

(4.64)\begin{equation} \hat{S}^w = \hat{\mathcal{L}}^w +\underbrace{\hat{h}^g\phi^c\cos\zeta\frac{\partial\hat{h}^w}{\partial\hat{x}}}_{Buoyancy} -\underbrace{(\hat{h}^w-\phi^c\hat{h}^g)\cos\zeta\frac{\partial\hat{b}}{\partial\hat{x}}}_{Topographic\ gradient} +\underbrace{\frac{2\varepsilon}{Re}\frac{\partial}{\partial\hat{x}}\left( (\hat{h}^w-\phi^c\hat{h}^g) \frac{\partial\hat{\bar{u}}^w}{\partial\hat{x}}\right).}_{Viscous\ term} \end{equation}

4.7. Unified dimensional system of depth-averaged equations

The leading-order system (4.10), (4.11), (4.42) and (4.49) in the undersaturated regime, and (4.10), (4.12), (4.57) and (4.63) in the oversaturated regime can be made dimensional by reversing the scalings

(4.65a–g)\begin{align} \left.\begin{gathered} x = L\hat{x},\quad (h^\nu, s^\nu, b)=H(\hat{h}^\nu, \hat s^\nu, \hat{b}),\quad t=(L/\sqrt{gH})\hat{t}, \quad u^\nu=\sqrt{gH}\hat{u}^{\nu}, \\ (S^\nu, \mathcal{L}^\nu) = (gH^2/L)(\hat{S}^\nu, \hat{\mathcal{L}}^\nu),\quad C^d=(\varrho^{w\star}\sqrt{g}/\sqrt{H})\hat C^d,\quad \psi^\nu=H\sqrt{gH}\hat\psi^\nu. \end{gathered}\right\} \end{align}

They can also be written in a unified form by introducing the variable

(4.66)\begin{equation} \mathcal{H}=\text{min}(h^w, h^g)/h^w, \end{equation}

which determines the proportion of the water height that is occupied by grains. As shown in figure 2(b,c), this implies that

(4.67)\begin{equation} \mathcal{H} =\begin{cases} 1, & h^g>h^w, \quad \text{(undersaturated)},\\ h^g/h^w, & h^g\le h^w, \quad \text{(oversaturated)}. \end{cases} \end{equation}

It follows that the unified system of dimensional equations can be expressed as

(4.68)$$\begin{gather} \frac{\partial}{\partial t}(h^g\phi^c) + \frac{\partial}{\partial x}(h^g\phi^c\bar{u}^g) = 0, \end{gather}$$

(4.69)$$\begin{gather}\frac{\partial }{\partial t}(h^w(1-\phi^c\mathcal{H})) +\frac{\partial}{\partial x}(h^w(1-\phi^c\mathcal{H})\bar{u}^w) = 0, \end{gather}$$

(4.70)$$\begin{gather}\frac{\partial}{\partial t}(h^g\phi^c\bar{u}^g) + \frac{\partial}{\partial x}\left(\chi^g h^g\phi^c(\bar{u}^g)^2 + \frac{1}{2}(h^g)^2\phi^c g\cos\zeta - \frac{1}{2}\gamma\phi^c(h^w\mathcal{H})^2 g\cos\zeta \right)=S^g, \end{gather}$$

(4.71)$$\begin{gather}\frac{\partial}{\partial t}(h^w(1-\phi^c\mathcal{H})\bar{u}^w) + \frac{\partial }{\partial x} \left(\chi^w h^w(1-\phi^c\mathcal{H})(\bar{u}^w)^2 + \frac{1}{2}(h^w)^2 g\cos\zeta \right)=S^w, \end{gather}$$

where the terms on the left-hand side of (4.70) and (4.71) represent the momentum transport and pressure gradients of each phase. The source terms on the right-hand side of (4.70) and (4.71) take the form

(4.72)\begin{align} S^g&=\mathcal{L}^g -\underbrace{\gamma\phi^c\mathcal{H} h^w g\cos\zeta\frac{\partial h^w}{\partial x}}_{Buoyancy} -\underbrace{h^g\phi^c g\cos\zeta\frac{\partial b}{\partial x},}_{Topography\ gradient} \end{align}

(4.73)\begin{align} S^w &= \mathcal{L}^w +\underbrace{\phi^c\mathcal{H} h^w g\cos\zeta\frac{\partial h^w}{\partial x}}_{Buoyancy} -\underbrace{(1-\phi^c\mathcal{H}) h^w g\cos\zeta\frac{\partial b}{\partial x}}_{Topography\ gradient} \nonumber\\ &\quad +\underbrace{2\nu^w\frac{\partial}{\partial x} \left(h^w(1-\phi^c\mathcal{H})\frac{\partial \bar{u}^w}{\partial x}\right),}_{Viscous\ term} \end{align}

where the kinematic viscosity of water

(4.74)\begin{equation} \nu^w=\eta^w/\varrho^{w\star}. \end{equation}

The dimensional leading order balances are

(4.75)$$\begin{gather} \mathcal{L}^g=\underbrace{h^g\phi^c g\sin\zeta}_{Gravity} -\underbrace{\frac{\bar{u}^g}{|\bar{u}^g|}\left(1-\gamma\mathcal{H} \frac{h^w}{h^g}\right)\mu^b h^g\phi^c g\cos\zeta}_{Basal\ friction} +\underbrace{\frac{C^d}{\varrho^{g\star}}\left(\psi^w(s^\nu)-\psi^g(s^\nu)\right)}_{Darcy\ drag}, \end{gather}$$

(4.76)$$\begin{gather}\mathcal{L}^w=\underbrace{h^w(1-\phi^c\mathcal{H})g\sin\zeta}_{Gravity} -\underbrace{C^w\bar{u}^w\vert\bar{u}^w\vert}_{Basal\ friction} -\underbrace{\frac{C^d}{\varrho^{w*}}\left(\psi^w(s^\nu)-\psi^g(s^\nu)\right),}_{Darcy\ drag} \end{gather}$$

where the internal interface height in the Darcy drag terms is defined as

(4.77)\begin{equation} s^\nu =\mathcal{H} h^w+b=\begin{cases} s^w, & h^g>h^w, \quad \text{(undersaturated)},\\ s^g, & h^g\le h^w, \quad \text{(oversaturated)}. \end{cases} \end{equation}

Equations (4.66)–(4.77) provide a convenient way of expressing the equations in all four regimes shown in figure 2, i.e. the dry granular front, the undersaturated region, the oversaturated region and the watery tail. In the undersaturated regime, the buoyancy terms in (4.72) and (4.73) can be written in conservative form. In this case, the whole system is in conservative form, which is advantageous for numerical methods (Kurganov & Tadmor Reference Kurganov and Tadmor2000) as well as for deriving jump conditions (Chadwick Reference Chadwick1976) at internal shocks. In the oversaturated regime, the buoyancy terms can not be written in conservative form, which is potentially a problem for numerical methods and if shocks develop. However, the presence of the viscous term in (4.73) will suppress the formation of water shocks, since the higher order gradient terms will smooth out any discontinuities. If granular shocks develop during oversaturated flows, it may also be necessary to add a viscous term to the depth-averaged granular momentum balance (see e.g. Gray & Edwards Reference Gray and Edwards2014).

4.8. Velocity profiles through the flow depth

To evaluate the shape factors in (4.70) and (4.71) and the Darcy drag terms in (4.75) and (4.76), it is necessary to make some assumptions about the velocity profiles through the depth of grains and the water. These assumptions must take account of the layered structure of the flow (4.6)–(4.8) in the undersaturated and oversaturated regimes (figure 2), and be consistent with the definitions of the depth-averaged velocities in (4.9a,b). This paper introduces shear into the model by defining the velocity for phase $\nu$ of the form

(4.78)\begin{equation} u^\nu=u_c^\nu(x,t)\left(\alpha^\nu + 2(1-\alpha^\nu)\left(\frac{z-b}{h^\nu}\right)\right), \end{equation}

where $u_c^\nu$ is independent of $z$ and the parameter $\alpha ^\nu \in [0,1]$ allows the velocity profile to vary from simple shear with no basal slip ($\alpha ^\nu =0$) to plug flow ($\alpha ^\nu =1$) (Gray & Thornton Reference Gray and Thornton2005; Gray & Kokelaar Reference Gray and Kokelaar2010; Johnson et al. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012). It would also be possible to use a power law or a Bagnold velocity profile (Turnbull, Bowman & McElwaine Reference Turnbull, Bowman and McElwaine2015; Baker et al. Reference Baker, Johnson and Gray2016b). The shear velocity profile (4.78) contrasts with widely used debris-flow models (e.g. Iverson & Denlinger Reference Iverson and Denlinger2001; Pitman & Le Reference Pitman and Le2005; Pelanti, Bouchut & Mangeney Reference Pelanti, Bouchut and Mangeney2008; Pudasaini Reference Pudasaini2012; Iverson & George Reference Iverson and George2014; Meng & Wang Reference Meng and Wang2016), which assume plug flow for both the grain and water phases.

This paper will show that velocity shear is sufficient to produce a dry granular front without the need to invoke ideas about particle segregation (Iverson Reference Iverson2003; Gray Reference Gray2018). Consider then the linear velocity profiles with depth given by (4.78). Since the volume fraction of grains is always equal to the constant $\phi ^c$ by (4.6), and the volume fraction of water is equal to the constant $1-\phi ^c$ in the undersaturated regime (4.7), it follows from (4.9a,b) that

(4.79a,b)\begin{equation} u_c^g=\bar{u}^g,\quad u_c^w=\bar{u}^w. \end{equation}

The corresponding granular and water shape factors, defined in (4.38) and (4.45), are

(4.80a,b)\begin{equation} \chi^g=\tfrac{4}{3}-\tfrac{2}{3}\alpha^g+\tfrac{1}{3}(\alpha^g)^2,\quad \chi^w=\tfrac{4}{3}-\tfrac{2}{3}\alpha^w+\tfrac{1}{3}(\alpha^w)^2, \end{equation}

respectively. For these cases, the shape factor $\chi ^\nu =1$ for plug flow and $4/3$ for simple shear. The situation for the oversaturated regime is more complex, because the volume fraction of water (4.8) is also height dependent. Substituting (4.8) and (4.78) into (4.9b) implies that

(4.81)\begin{equation} u_c^w=\lambda^w\bar{u}^w, \end{equation}

where

(4.82)\begin{equation} \lambda^w=\frac{(1-\phi^c\mathcal{H})}{1-\phi^c(\alpha^w\mathcal{H}+(1-\alpha^w)\mathcal{H}^2)}. \end{equation}

Note that when $\mathcal {H}=1$, $\lambda ^w=1$ and (4.81) reduces to the undersaturated case. Solving (4.45) for the shape factor and substituting the depth-averaged water velocity (4.9b) implies that

(4.83)\begin{equation} \chi^w=\frac{\overline{\phi^w(u^w)^2}\bar\phi^w}{(\overline{\phi^w {u}^w})^2}. \end{equation}

The oversaturated water shape factor is therefore

(4.84)\begin{equation} \chi^w=\frac{ (4\chi_0^w -2\chi_1^w\alpha^w +\chi_2^w(\alpha^w)^2)(1-\mathcal{H}\phi^c)} {3(1-\phi^c(\alpha^w\mathcal{H}+(1-\alpha^w)\mathcal{H}^2))^2}, \end{equation}

where the coefficients

(4.85a–c)\begin{equation} \left.\begin{gathered} \chi_0^w= 1-\mathcal{H}^3\phi^c,\quad \chi_1^w= 1+3\mathcal{H}^2\phi^c-4\mathcal{H}^3\phi^c,\\ \chi_2^w= 1-3\mathcal{H}\phi^c+6\mathcal{H}^2\phi^c-4\mathcal{H}^3\phi^c. \end{gathered}\right\} \end{equation}

Thus, (4.84) reduces to the undersaturated case (4.80b) when $\mathcal {H}=1$. For the velocity profiles (4.78), it follows from the definitions of the streamfunctions (4.39) and (4.53) that in the undersaturated regime,

(4.86)\begin{equation} \psi^w(s^w)-\psi^g(s^w)=h^w\bar{u}^w-\bar{u}^g\left(\alpha^g h^w+(1-\alpha^g)\frac{(h^w)^2}{h^g}\right). \end{equation}

Similarly, in the oversaturated regime at the internal grain surface,

(4.87)\begin{equation} \psi^w(s^g)-\psi^g(s^g)=\lambda^w\bar{u}^w\left(\alpha^w h^g+(1-\alpha^w)\frac{({h}^g)^2}{{h}^w}\right)-h^g\bar{u}^g, \end{equation}

where the factor $\lambda ^w$ enters due to the non-uniform concentration of water with $z$. This paper investigates the simplest possible shearing debris flow model, by assuming simple shear in the granular phase and plug flow in the water, i.e.

(4.88a,b)\begin{equation} \alpha^g=0,\quad \hbox{and}\quad \alpha^w=1. \end{equation}

With these assumptions, it follows from (4.86) and (4.87) that the difference of the streamfunctions, which arises in (4.75) and (4.76), can be expressed in the unified form

(4.89)\begin{equation} \psi^w(s^\nu)-\psi^g(s^\nu)= h^w\mathcal{H}\bar u^w-\frac{(h^w\mathcal{H})^2}{h^g}\bar u^g. \end{equation}

The assumption of plug flow in the water implies that $\chi ^w=1$, and for simplicity, it is further assumed that the grain shape factor $\chi ^g=1$, which is a common assumption in depth-averaged debris flow models (e.g. Iverson & Denlinger Reference Iverson and Denlinger2001; Pitman & Le Reference Pitman and Le2005; Pelanti et al. Reference Pelanti, Bouchut and Mangeney2008; Pudasaini Reference Pudasaini2012; Iverson & George Reference Iverson and George2014; Meng & Wang Reference Meng and Wang2016).

4.9. Closure of friction coefficients $\mu ^b$ and $C^w$

Considerable progress has been made in the field of the granular flows with the development of $\mu (I)$-rheology (GDR-MiDi 2004; Jop, Forterre & Pouliquen Reference Jop, Forterre and Pouliquen2006), where the friction coefficient $\mu$ is a monotonically increasing function of inertial number $I$. This functional dependence was determined from measurements of steady-uniform flows between the angles $\zeta _1$ and $\zeta _2$ on an inclined plane (Pouliquen Reference Pouliquen1999b; Pouliquen & Forterre Reference Pouliquen and Forterre2002). The ensuing empirical friction law is given by

(4.90)\begin{equation} \mu^b(Fr, h^g) = \mu_1 + \frac{\mu_2 - \mu_1}{1 + \displaystyle\frac{\beta h^g}{\mathscr{L}Fr}}, \end{equation}

where $\mu _1 = \tan \zeta _1$ and $\mu _2 = \tan \zeta _2$. The parameter $\beta$ is an empirical constant and $\mathscr {L}$ has the dimensions of a length and is dependent on the properties of the grains and on the bed roughness. The Froude number in (4.90) is defined as

(4.91)\begin{equation} Fr = \frac{|\bar{u}^g|}{\sqrt{h^g g\cos\zeta}}. \end{equation}

Maurin et al. (Reference Maurin, Chauchat and Frey2016) used a coupled fluid–discrete-element method to perform numerical simulations of turbulent bedload transport. Contrary to expectations, they found that for dense submerged granular flows, the effect of the interstitial fluid viscosity was negligible. Moreover, they found that the $\mu (I)$-rheology was able to collapse the simulated shear-to-normal stress ratio and the solids volume fraction over an unexpectedly wide range of inertial numbers. This lends weight to the use of the granular friction (4.90) to describe the basal friction of the grains. It enters the theory through the leading order granular source term defined in (4.75), which includes a factor $1-\gamma \mathcal {H}h^w/h^g$ to account for the reduced granular normal stress in the submerged grains.

The water bed friction coefficient is based on the Manning equation (Manning Reference Manning1891; Chaudhry Reference Chaudhry2008) for open channel flows, which implies that $C^w=g n^2/(h^w)^{1/3}$, where $n$ is the Manning coefficient (Chertock et al. Reference Chertock, Cui, Kurganov and Wu2015). This formula is modified to account for the presence of grains by defining it to be

(4.92)\begin{equation} C^w = (1-\phi^c\mathcal{H})\frac{g n^2}{(h^w)^{1/3}}. \end{equation}

This ensures that $C^w$ reduces to the classical form for shallow water flows in the absence of grains.

5.1. Description and summary of observations

The moving-bed apparatus consisted of a rectangular flume that was $2$ m in length and had a width of $50$ mm, with walls at the ends to retain all the material, see figure 3. The channel bed was made rough by using a corrugated nylon belt, driven by a toothed drive wheel connected by another toothed belt to a variable-speed electric motor and controller. The range of channel slopes used in the experiments was from $5^\circ$ to $19^\circ$ and the bed speeds ranged from 0.25 to $1.17$ m s$^{-1}$. Room temperature tap water was used and the plastic grains were $4$-mm-long cylinders cut from $4$-mm diameter polyvinyl chloride (PVC) rod with intrinsic density $\varrho ^{g\star }=1400$ kg m$^{-3}$. The maximum grain volume fraction was $0.56$. Approximately 10 % of the particles were painted white to act as tracers, and videos were used to track individual particle paths. In Davies’ (Reference Davies1988) experiments, the surge remains stationary, while the bed moves upward at a constant velocity $U$. The main experimental findings relevant to the current study are outlined below.

(i) Adding grains caused the flow depth to increase and a curved front appeared in the channel. As more grains were added, and the grain volume exceeded the threshold for uniform flow at a given speed and slope angle, the excess grains moved to the front and caused a thicker dry snout to develop. This dry bulbous flow head is similar in form to those observed in large scale debris flows.

(ii) As more grains were added, a substantial high-concentration stationary surge was formed at the lower end of the channel. This surge had a uniform depth body that extended from the curved front to an attenuated tail. The depth of this uniform body was not very sensitive to the grain volume, as shown in figure 9 of Davies (Reference Davies1990). The flows usually formed on a relatively gentle slope or when the bed speed was relatively low. An example is shown in online supplementary movie 2, at a slope angle of 15$^\circ$ and bed speed $U=0.2$ m s$^{-1}$.

(iii) When the bed speed was raised to a relatively high value, the stationary surge was curved rather than being of uniform depth (see e.g. online supplementary movie 3 at 15$^\circ$ inclination and bed speed $U=0.3$ m s$^{-1}$). The longitudinal profiles showed that a dry front is followed by an undersaturated and oversaturated body that degenerates down to pure water in the tail, as sketched in figure 2. This longitudinal profile is similar to those observed in the USGS debris-flow flume tests (Iverson et al. Reference Iverson, Logan, Lahusen and Berti2010; Johnson et al. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012) and in the small scale experimental curved channel of Scheidl et al. (Reference Scheidl, McArdell and Rickenmann2015).

(iv) Tracers show that grains near the free surface are sheared towards the front, are overrun there and transported back towards the tail near the base of the flow. Here they rise up to the surface again and are sheared forward to form a recirculating cell that conserves its mass.

5.2. Experimental parameters

Davies (Reference Davies1988) provided the values of the slope angle $\zeta$, the solids volume fraction $\phi ^c$, the density ratio $\gamma$, the particle diameter $d$ and the water viscosity $\eta ^w$. These experimental parameters are summarised in table 1. At the time, the $\mu (I)$-rheology had not been developed, so exact values of the parameters in the friction law (4.90) were not determined. However, Armanini et al. (Reference Armanini, Capart, Fraccarollo and Larcher2005) has conducted experiments using similar plastic cylinders, and Berzi & Jenkins (Reference Berzi and Jenkins2008) show that these cylinders are more frictional than glass beads (Pouliquen Reference Pouliquen1999b), but less frictional than sand. The values of $\zeta _1$ and $\zeta _2$ are therefore assumed to be slightly larger than those for glass beads presented in table I of Pouliquen (Reference Pouliquen1999b). The value of $\beta$ follows from Pouliquen (Reference Pouliquen1999b) and Gray & Edwards (Reference Gray and Edwards2014). The measurements of Forterre & Pouliquen (Reference Forterre and Pouliquen2003) show that the frictional length scale $\mathscr {L}$ depends on the grain diameter $d$, and that it is approximately equal to $1.65d$ for glass beads in dry granular flows. The parameters measured for dry flows are adopted here without modification. However, it should be noted that the basal friction in (4.75) includes a factor $1-\gamma \mathcal {H}h^w/h^g$ to account for the reduced granular normal stress when the grains are submerged in water. This is consistent with Maurin et al. (Reference Maurin, Chauchat and Frey2016) who measured quite large shear-to-normal stress ratios in their coupled fluid-discrete-element method simulations of bed load transport.

Table 1.

Physical parameters in the experiments and the computation.

The steady-uniform water flow at the tail in Davies’ (Reference Davies1988) experiments implies a balance between the gravitational force and the Chézy friction in (4.76),

(5.1)\begin{equation} h^w_0g\sin\zeta = gn^2(\bar{u}^w_0)^2/(h^w_0)^{{1}/{3}}, \end{equation}

where $h^w_0$ and $\bar {u}^w_0$ are the steady-uniform thickness and speed of the bed, respectively. It follows that Manning's coefficient is

(5.2)\begin{equation} n=\sqrt{\frac{(h^w_0)^{{4}/{3}}\sin\zeta}{(\bar{u}^w_0)^2}}. \end{equation}

For Davies’ (Reference Davies1988) experiment shown in the online supplementary material (movie 4), the slope inclination angle $\zeta =19^\circ$, the bed speed $\bar {u}^w_0=0.3$ m s$^{-1}$ and the steady-uniform water depth in the tail lies between one- and two-grain diameters. Assuming $h^w_0=5.5$ mm, Manning's coefficient is $n=0.0593$ s m$^{-(1/3)}$. All these additional parameters are also summarised in table 1. In particular, none of them are fitting parameters.

6.1. Case 1: formation of bulbous heads in the undersaturated regime

The online supplementary movie of Davies’ (Reference Davies1988, Reference Davies1990) experiment shows that when the bed speed is relatively small, a dry bulbous head develops that is connected to an undersaturated body of approximately constant depth. This configuration is equivalent to the case where the surge is steadily translating on an unmoved bed. It is therefore useful to transform the governing equations (4.68)–(4.71) into a frame of reference moving with speed $U$, by making the coordinate transformation

(6.1a,b)\begin{equation} \xi = x - Ut,\quad \tau = t. \end{equation}

The steady-state grain and water mass balances in the moving frame are therefore

(6.2)$$\begin{gather} \frac{\text{d}}{\text{d}\xi}(h^g\phi^c(\bar{u}^g-U))=0, \end{gather}$$

(6.3)$$\begin{gather}\frac{\text{d}}{\text{d}\xi}(h^w(1-\phi^c\mathcal{H})(\bar{u}^w-U))=0, \end{gather}$$

respectively. Integrating (6.2) and (6.3), subject to the condition that $h^g=0$ and $h^w=0$ at the front, implies that

(6.4a,b)\begin{equation} h^g\phi^c(\bar{u}^g-U)=0,\quad h^w(1-\phi^c\mathcal{H})(\bar{u}^w-U)=0. \end{equation}

These are satisfied either when $h^g$ or $h^w$ are identically zero, or when

(6.5)\begin{equation} \bar{u}^g=\bar{u}^w=U, \end{equation}

i.e. the depth-averaged velocities of the grains and the water are equal to the moving frame speed.

Since the depth-averaged water velocity is independent of $\xi$ and $b=0$ in Davies’ (Reference Davies1988, Reference Davies1990) experiments, the viscous term and the topography gradients vanish in the source terms (4.72)–(4.73). Assuming that the flow is undersaturated ($\mathcal {H}=1$) and that $\chi ^\nu =1$, the steady-state depth-averaged momentum balances (4.70)–(4.71) in the travelling frame reduce to

(6.6)\begin{align} & h^g\phi^c(\bar{u}^g-U)\frac{\text{d}\bar{u}^g}{\text{d}\xi}+h^g\phi^c g\cos\zeta \frac{\text{d} h^g}{\text{d}\xi}=h^g\phi^c g\sin\zeta \nonumber\\ &\quad - \mu^b h^g \phi^c g\cos\zeta (1-\gamma{h^w}/{h^g}) +\frac{180\eta^w h^w(\phi^c)^2}{\varrho^{g\star}(1-\phi^c)d^2}(\bar{u}^w-\bar{u}^g{h^w}/{h^g}), \end{align}

(6.7)\begin{align} & h^w(1-\phi^c)(\bar{u}^w-U)\frac{\text{d}\bar{u}^w}{\text{d}\xi} +h^w(1-\phi^c)g\cos\zeta \frac{\text{d} h^w}{\text{d}\xi}=h^w(1-\phi^c)g\sin\zeta \nonumber\\ &\quad -(1-\phi^c)g(n \bar{u}^w)^2/(h^w)^{{1}/{3}} -\frac{180\eta^w h^w(\phi^c)^2}{\varrho^{w\star}(1-\phi^c)d^2}(\bar{u}^w-\bar{u}^g{h^w}/{h^g}), \end{align}

where it is implicitly assumed that the sidewall friction in Davies’ (Reference Davies1988, Reference Davies1990) experiments is sufficiently small that this one-dimensional theory can be applied.

Using (6.5) to eliminate $\bar {u}^g$ and $\bar {u}^w$ in favour of $U$, it follows that (6.6) and (6.7) reduce to a pair of ordinary differential equations (ODEs) for the thicknesses $h^g$ and $h^w$

(6.8)\begin{gather} \frac{\text{d} h^g}{\text{d}\xi} =\tan\zeta - \mu^b(1-\gamma h^w/h^g) + \frac{180\eta^w\phi^c}{\varrho^{g\star}g\cos\zeta (1-\phi^c)d^2}\frac{h^w}{h^g}(1-h^w/h^g)U, \end{gather}

(6.9)\begin{gather}\frac{\text{d} h^w}{\text{d}\xi}=\tan\zeta-(nU)^2/((h^w)^{{4}/{3}}\cos\zeta) -\frac{180\eta^w(\phi^c)^2}{\varrho^{w\star}g\cos\zeta (1-\phi^c)^2d^2}(1-{h^w}/{h^g})U. \end{gather}

The grain thickness gradient in (6.8) is balanced by gravity, the buoyancy-reduced friction and the Darcy drag, while the water thickness gradient in (6.9) balances gravity, the Chézy drag and the Darcy drag.

In a dry snout, the position of the grain front $\xi ^g_f$ lies further downslope than the water front $\xi ^w_f$, and within this region, the water depth $h^w$ vanishes. As a result, (6.8) reduces to

(6.10)\begin{equation} \frac{\text{d} h^g}{\text{d}\xi} = \tan\zeta - \mu^b. \end{equation}

This equation was originally derived by Pouliquen (Reference Pouliquen1999a) to compute the shape of a dry granular front, and admits an exact solution (Gray & Ancey Reference Gray and Ancey2009) for the friction law (4.90). As $h^w\rightarrow 0$, the Chézy drag is singular in the water ODE (6.9). It follows that at the water front, the dominant balance is between the thickness gradient and Chézy drag, i.e.

(6.11)\begin{equation} \frac{\text{d} h^w}{\text{d}\xi}={-}\frac{(nU)^2}{(h^w)^{{4}/{3}}\cos\zeta}. \end{equation}

Integrating (6.11) subject to the condition that $h^w=0$ at $\xi = \xi ^w_f$, implies that

(6.12)\begin{equation} h^w=\left[\frac{7}{3}\frac{(nU)^2}{\cos\zeta}(\xi^w_f-\xi)\right]^{{3}/{7}}, \end{equation}

which determines the water free surface in the neighbourhood of the front.

Without loss of generality, the granular front is assumed to lie at $\xi ^g_f=0$, and (6.10) is integrated upstream using Matlab's ode45 initial value problem solver to determine $h^g$. This solution is valid until the water front is reached at $\xi =\xi ^w_f$, which is a free parameter in the problem. The approximate solution (6.12) is used to determine the water surface for small $h^w$, and then the ODEs (6.8) and (6.9) are integrated upstream to determine both $h^w$ and $h^g$. Figure 4 shows two solution profiles for a 15$^\circ$ slope and a bed speed $U=0.262$ m s$^{-1}$. In figure 4(a), the dry front is assumed to be 18 grain diameters long, and both the water and grain fronts are connected to the far field steady states by monotonically decreasing curves. Such solutions resemble those of monodisperse dry granular flows (Pouliquen Reference Pouliquen1999a), and agree with the experimental observation of Davies (Reference Davies1988) that for low speed flows on shallow inclines, the dry snout is connected to an upstream uniform-depth undersaturated flow. Although the depth-averaged velocities of grains and water are equal, the local velocities are not. In the lower, saturated part of the flow, the water flows downslope more quickly than the grains, and helps to drag the granular phase downslope. In the upper unsaturated part of the flow, velocity shear results in the grains (alone) being transported downslope faster than either the water or grains in the lower saturated part of the flow. This differs from the model of Berzi et al. (Reference Berzi, Jenkins and Larcher2010), in which the local grain and water velocities are assumed to be equal in the steady uniform flow far upstream of the front. This results in their solutions approaching an exactly saturated state $h^w=h^g$ upstream.

Figure 4.

Longitudinal profiles of the grain depth (brown lines) and water depth (blue lines) with $\zeta =15^\circ$ and $U=0.262$ m s$^{-1}$ for (a) $\xi ^g_f-\xi ^w_f=18d$ and (b) $\xi ^g_f-\xi ^w_f=36d$. The blue shaded regions represent water-saturated material and a few schematic grains are shown in the granular region to make it easier to identify. The flow direction is from left to right and the insets show a close-up view of the front.

When the dry front is doubled in length, to $\xi ^g_f-\xi ^w_f=36$ grain diameters (figure 4b), there is much more resistance from the dry snout, and the granular front has to have a steeper gradient to drive it downslope. As a result, both the granular and the water free surfaces develop a pronounced bulge just upstream of the bump, as shown in figure 4(b). In the dry snout, the component of gravitational acceleration is smaller than the basal friction, so that there is a negative thickness gradient. Upstream of the water front, the basal granular friction diminishes gradually, due to the buoyancy imparted by the interstitial fluid. This implies that the driving force gets close to the buoyancy-reduced grain friction, and exceeds it just upstream of the peak grain depth, implying the gradient of $h^g$ is positive there. In the wet region, the driving force, i.e. the downslope component of gravitational acceleration, balances the resistance due to Chézy drag and the Darcy interaction force for the water phase. Near the water front, the water thickness is sufficiently thin that the Chézy drag outweighs the gravity-driven force, and a steep water front forms with a negative thickness gradient. Note that far upstream, the two solutions shown in figure 4(a,b) approach the same steady uniform states.

6.2. Case 2: formation of bulbous heads and pure watery tails

When the bed speed $U$ is increased, the granular free surface becomes more curved, as shown in the supplementary online movies 3 and 4 of Davies' (Reference Davies1988, Reference Davies1990) experiments. As a result, the flow develops a dry snout, an undersaturated region that smoothly transitions into an oversaturated flow and then a watery tail, as shown schematically in figure 2. This is qualitatively different from the experiments and solutions described in § 6.1 for low bed speeds. In the travelling frame (6.1a,b), the steady-state mass balances (6.2) and (6.3), together with the grain and water free frontal condition, imply that either $h^g$ and $h^w$ are identically zero or the depth-averaged grain and water velocities are equal to the bed speed

(6.13)\begin{equation} \bar{u}^g=\bar{u}^w=U. \end{equation}

As in § 6.1, the solution starts at the front of the dry snout, which is assumed to lie at $\xi ^g_f=0$. The dry granular momentum balance (6.10) is then integrated back to the water front at $\xi ^w_f$ using Matlab's ode45 initial value solver. The approximate solution (6.12) at the water front is used to initiate the solution of the undersaturated ODEs (6.8) and (6.9). These are integrated upstream until the transition point is reached between the undersaturated and oversaturated regimes, where $h^g=h^w$. To advance the solution further, ODEs need to be derived for the oversaturated regime. Assuming that the shape factors $\chi ^\nu$ equal unity, the steady-state oversaturated depth-averaged momentum balances (4.70) and (4.71) in the travelling frame reduce to

(6.14)\begin{align} & h^g\phi^c(\bar{u}^g-U)\frac{\text{d}\bar{u}^g}{\text{d}\xi}+h^g(1-\gamma)\phi^c g\cos\zeta\frac{\text{d} h^g}{\text{d}\xi}+\gamma h^g\phi^c g\cos\zeta\frac{\text{d} h^w}{\text{d}\xi} =h^g\phi^c g\sin\zeta \nonumber\\ &\quad - \mu^b h^g\phi^c g\cos\zeta(1-\gamma)+ \frac{180\eta^w h^g(\phi^c)^2}{\varrho^{g\star}(1-\phi^c)d^2}(\bar{u}^w-\bar{u}^g), \end{align}

(6.15)\begin{align} & (h^w-h^g\phi^c)(\bar{u}^w-U)\frac{\text{d}\bar{u}^w}{\text{d}\xi} +(h^w-h^g\phi^c)g\cos\zeta\frac{\text{d} h^w}{\text{d}\xi}=(h^w-h^g\phi^c)g\sin\zeta \nonumber\\ &\quad -(h^w-h^g\phi^c)\frac{g(n\bar{u}^w)^2}{(h^w)^{{4}/{3}}}- \frac{180\eta^w h^g(\phi^c)^2}{\varrho^{w\star}(1-\phi^c)d^2}(\bar{u}^w-\bar{u}^g), \end{align}

respectively. Substituting the velocity condition (6.13) into (6.14) and (6.15) eliminates the acceleration and Darcy drag terms, even though locally, the grains, water and moving frame are all moving at different speeds. The resulting pair of equations can be solved for the gradients of the grain and water thicknesses

(6.16)$$\begin{gather} \frac{\text{d} h^g}{\text{d}\xi} =\tan\zeta - \mu^b + \frac{\gamma}{1-\gamma}\frac{(nU)^2}{(h^w)^{{4}/{3}}\cos\zeta}, \end{gather}$$

(6.17)$$\begin{gather}\frac{\text{d} h^w}{\text{d}\xi}=\tan\zeta-\frac{(nU)^2}{(h^w)^{{4}/{3}}\cos\zeta}. \end{gather}$$

The lubrication effect of the water phase on the grains is characterised by a buoyancy force that acts against the granular friction, and has the form of a Chézy drag. This is different from (6.8) in the undersaturated regime, where the lubrication effect is characterised by the buoyancy-reduced friction.

The ODEs (6.16) and (6.17) are integrated back from the transition point until $h^g$ reaches zero at $\xi ^g_t$. Since $h^g=0$ is a trivial solution of (6.14), the underlying system naturally degenerates to a single ODE (6.17) that can be solved for the shape of the pure water tail, upstream of $\xi ^g_t$. Figure 5 shows solutions for the free surfaces of the grains and the water at a slope angle $\zeta =15^\circ$, bed speed $U=0.39$ m s$^{-1}$ and for three values of the snout length $\xi ^g_f-\xi ^w_f$, which is a free parameter in the problem. Increasing the snout length results in an increase of the total volume of grains in the travelling wave. In the undersaturated regime, the water phase moves locally faster than the grains at the bottom of the flow, but slower than the dry surface grains, which are sheared forward to form the dry snout. In the oversaturated tail, the flow becomes sufficiently thin that the basal friction experienced by the water outweighs that acting on the grains. As a result, grains that reach the granular free surface are sheared forward to leave a watery upstream tail. This process can be seen in Davies’ (Reference Davies1988) experimental movies 3 and 4 in the online supplementary material. The solutions in figure 5 illustrate that it is possible to model the forward motion of the grains relative to the water, without the need to invoke the traditional view (in the debris flow community) that excess pore pressures are needed to drive the forward motion of grains in the oversaturated tail.

Figure 5.

Longitudinal profiles of the grain depth (brown lines) and water depth (blue lines) with $\zeta =15^\circ$ and $U=0.39$ m s$^{-1}$, for (a) $\xi ^g_f-\xi ^w_f=18d$, (b) $\xi ^g_f-\xi ^w_f=29.7485d$ and (c) $\xi ^g_f-\xi ^w_f=36d$. The blue shaded regions represent water-saturated material and a few schematic grains are shown in the granular region to make it easier to identify. The ‘$\circ$’ symbols in panel (b) represent the granular free surface measured by Davies (Reference Davies1990) at the same bed speed and inclination. The flow direction is from left to right.

For the solutions shown in figures 4 and 5, the snout length is the natural parameter to use when solving the ODEs, and the total volume of grains and water emerges as part of the solution. Conversely, in physical experiments, the total volume of grains and water is prescribed, and one then has to iterate on the snout length to select the solution with the correct water and grain volumes. Figure 5(b) shows a comparison of the computed grain free surface and that measured by Davies (Reference Davies1988) at the same bed speed and slope angle. The computed amplitude, length and slope of the granular free surface in the tail are in good agreement with the experiments. The main deviation from the experimental measurements lies in the frontal profile, which may be due to a collisional shear layer that forms near the moving bed in the experiments. Detailed profiles of the water free surface are not reported in Davies (Reference Davies1988, Reference Davies1990). However, Davies (Reference Davies1990) does say that the water free surface is approximately one-grain diameter below the granular free surface in the main part of the unsaturated body. Figure 5(b) shows that the water free surface lies within two grain diameters of the grain free surface for the majority of the undersaturated region. The theory therefore captures the key qualitative features of the flow as well as providing a good quantitative match to Davies’ (Reference Davies1990) granular free surface data without any fitting parameters.

6.3. Recirculation of the grains

In Davies’ (Reference Davies1988) experiments, approximately 10 % of the grains were painted white to reveal the particle paths. This can also be simulated by reconstructing the two-dimensional velocity field from the depth-averaged velocity and tracing out the streamlines. It follows from (4.78), (4.79a,b) and (4.88a,b) that the downstream velocity of the grains is given by the linear shear profile

(6.18)\begin{equation} u^g=2\bar u^g\left(\frac{z-b}{h^g}\right). \end{equation}

The granular mass balance (2.2) with the assumption of constant solids volume fraction (4.6) implies that in the moving frame, the grains satisfy bulk incompressibility

(6.19)\begin{equation} \frac{\partial u^g}{\partial \xi} + \frac{\partial w^g}{\partial z} = 0. \end{equation}

This can be integrated through the depth $z$, subject to the condition that $w^g=0$ on the bottom $z=b$, to give the normal velocity

(6.20)\begin{equation} w^g = \frac{\bar{u}^g}{(h^g)^2}\frac{\text{d} h^g}{\text{d}\xi}(z-b)^2. \end{equation}

Following Gray & Ancey (Reference Gray and Ancey2009), the particle paths in the moving frame

(6.21a,b)\begin{equation} \frac{\text{d}\xi^p}{\text{d} t}=u^g-U,\quad \frac{\text{d} z^p}{\text{d} t}=w^g \end{equation}

can be reconstructed by solving

(6.22)\begin{equation} \frac{\text{d} z^p}{\text{d}\xi}=\frac{w^g}{u^g-U}. \end{equation}

This can be simplified by introducing the streamfunction coordinate

(6.23)\begin{equation} \psi(\xi,z(\xi))= \int_b^{z(\xi)}u^g(\xi, z^\prime)-U\,\text{d} z^\prime. \end{equation}

Assuming that $\psi ^p(\xi )$ is the particle path $z=z^p(\xi )$, Leibniz’s rule can be used to swap the order of differentiation and integration in the derivative $\textrm {d}\psi ^p/\textrm {d}\xi$, and hence, using the bulk incompressibility (6.19), the no normal velocity condition $w^g(b)=0$ and (6.22), this implies

(6.24)\begin{equation} \frac{\text{d}\psi^p}{\text{d}\xi}=(u^g-U)\frac{\text{d} z^p}{\text{d}\xi}-w^g=0. \end{equation}

The streamfunction is therefore constant along a given particle path. An explicit expression for $\psi$ can be obtained by substituting the downstream velocity component (6.18) into (6.23) and recalling that the depth-averaged velocity $\bar u^g$ is equal to the bed speed $U$ in the travelling wave solution (6.5), to give

(6.25)\begin{equation} \psi =\frac{(z-b)^2}{h^g}U-(z-b)U. \end{equation}

It follows that the streamfunction $\psi$ equals zero on both the free surface $z=s^g$ and the base $z=b$. This implies that there is a closed streamline that runs along the grain free surface and base, which forms a recirculating cell. Grains on the free surface are sheared towards the flow front, where they are over-run and remain stationary on the bottom while the flow passes over them. As the tail propagates past, these basal grains are picked up again and sheared towards the front.

Within the $\psi =0$ contour, there are a nested set of recirculating closed particle paths, as shown in figure 6(a). At each $\xi$, the streamfunction reaches a local minimum $\psi _{min}=-h^gU/4$ at $z=z_{min}=b+h^g/2$. This is also the point where the downslope velocity (6.18) is equal to the frame speed $U$, i.e. it is the no-mean-flow line in the frame of the wave. Grains above the no-mean-flow line migrate towards the front and grains below advect backwards in the moving frame and are overtaken by the tail. An explicit relation for the particle paths can be obtained by solving the quadratic equation (6.25) for $z$ to give

(6.26)\begin{equation} z =\begin{cases} b+\frac{1}{2}h^g+\frac{1}{2}\sqrt{(h^g)^2+4\psi h^g/U}, & z\ge z_{min},\\ b+\frac{1}{2}h^g-\frac{1}{2}\sqrt{(h^g)^2+4\psi h^g/U}, & z< z_{min}.\\ \end{cases} \end{equation}

Figure 6.

(a) Nested (green) streamlines within the granular region for the same conditions as illustrated in figure 5(b). The granular free surface is illustrated in brown and the no-mean-flow line in the frame of the wave is shown with a dashed line. The water is shaded in blue. Each streamline has two intersections with the no-mean-flow line (as illustrated by the black markers on one loop). Grains that lie above the no-mean-flow line move forwards towards the front of the flow, while grains below it move backwards towards the tail. (b) Measured particle paths in Davies’ (Reference Davies1990) experiments over a one-second time interval (green arrowed lines) and the measured free surface height (circles). (c) Predicted particle paths using the flow field in panel (a). In panel (c), the original positions of grains with numbers (1) and (2) are ($-20.1326d$, $7.4855d$) and ($-26.9030d$, $7.7656d$), respectively, which are different from those in panel (b) due to the fact that their measured positions are beyond the predicted grain surface.

The predicted grain paths (figure 6c) are in good agreement with Davies’ (Reference Davies1990) measured particle paths over a one-second time interval (figure 6b). Both measurements and prediction show that grains near the surface are moving forward towards the front, while those near the base are transported backward (in the moving frame) and rise up again to prevent the surge from losing mass. This clearly indicates that strong shear develops in the flow. Some discrepancies also exist. This is largely because the velocity profile in Davies’ (Reference Davies1990) experiments has the shear concentrated at the bottom of the flow and it is sheared less strongly near the free surface.