2.1. Eroding Bingham fluid

In the absence of entrainment, a two-dimensional, stationary, gravity-driven simple shear flow of a Bingham fluid exhibits a plug layer that extends from the surface of the flow at $z=h$ to a height $z=r$ above the bed:

The shear stress profile is linear,

the shear rate diminishes linearly from the bed to the height $r$ ,

and the velocity profile is parabolic in the shear layer if there is no sliding at the bed,

With entrainment, the velocity profile will remain constant in time in our idealized situation, but deviate from (2.8). The shear stress being zero at the surface $z=h$ , there will still be a plug layer. As there is no acceleration, the plug layer depth is unchanged from the expression (2.5). Upon applying (2.4) and (1.3) for a Bingham fluid to (2.2), we obtain

Temporarily assuming the entrainment speed $w_{e}$ to be given and applying the boundary condition $\dot{{\it\gamma}}(r)=0$ , we find the solution

From this we obtain the following velocity and shear-stress profiles if we assume no slip at the bed:

By expanding the exponential functions to second order, one easily verifies that (2.10)–(2.12) reduce to (2.6)–(2.8) in the limit $w_{e}\rightarrow 0$ .

These formulae suggest the following scaling for lengths, velocities and stresses, with dimensionless quantities designated by a tilde:

note that $\tilde{r}=1-\tilde{{\it\tau}}_{y}$ . Then the velocity and stress profiles and the momentum form factor (or Boussinesq coefficient), defined by

are given by the expressions listed in table 1.

Table 1.

Vertical profiles of shear stress and longitudinal velocity in steady-state flows of Bingham fluids with and without entrainment. All quantities are non-dimensionalized according to (2.13). In the limit of vanishing entrainment, i.e. $\tilde{w}_{e}\rightarrow 0$ , all profiles with entrainment reduce to the corresponding profiles without entrainment. The Boussinesq form factor $f_{B}$ given by (2.14) is evaluated for $\tilde{r}=0.5$ , and for $\tilde{{\it\tau}}_{c}=0.7$ in the case of entrainment. The value of $\tilde{w}_{e}$ indicated in brackets corresponds to the physically consistent entrainment rate.

It remains to implement the stress boundary condition $\hat{{\it\sigma}}_{xz}(0)=\hat{{\it\tau}}_{c}$ or $\tilde{{\it\sigma}}_{xz}(0)=\tilde{{\it\tau}}_{c}$ . Using the pertinent expression from table 1 at $\tilde{z}=0$ yields the transcendental equation

for the dimensionless erosion rate $\tilde{w}_{e}$ in terms of the dimensionless yield strength $\tilde{{\it\tau}}_{y}$ and the dimensionless shear strength of the snow cover, $\tilde{{\it\tau}}_{c}$ . Note the physical requirement $0\leqslant \tilde{{\it\tau}}_{y}<\tilde{{\it\tau}}_{c}<1$ . Equation (2.15) can be recast as

where

Here, $v$ takes values in the range $-\infty <v<-1$ so that $v\text{e}^{v}\in D:= ]-1/\text{e},0[$ . The general solution can be expressed in terms of Lambert’s non-elementary $W$ function, which is defined as the solution $x(y)$ of the equation $x\text{e}^{x}=y$ . In the possible range of $v$ , $W$ is two-valued, with $W_{0}(-1/\text{e})=W_{-1}(-1/\text{e})=-1$ and $W_{0}(D)=\, ]-1,0[$ , $W_{-1}(D)= ]-\infty ,-1[$ . The equation for the dimensionless erosion speed now reads $w=W(v\text{e}^{v})$ . As $v<w$ , $v$ corresponds to the lower branch $W_{-1}$ and $w=W_{0}(v\text{e}^{v})$ . The erosion speed thus becomes

the dimensionless quantity $\tilde{w}_{e}$ is plotted as a function of $\tilde{{\it\tau}}_{y}$ and $\tilde{{\it\tau}}_{c}$ in figure 3(a). As figure 3(b) shows, favourable conditions for entrainment lead to low flow velocities.

Figure 3.

Quasi-stationary flow of a Bingham fluid with non-dimensional yield strength $\tilde{{\it\tau}}_{y}$ over an erodible bed with dimensionless shear strength $\tilde{{\it\tau}}_{c}$ : (a) dimensionless entrainment rate $\tilde{w}_{e}$ and (b) ratio of surface velocities of entraining and non-entraining flows. Quasi-stationary eroding flow is not possible for $\tilde{{\it\tau}}_{y}>\tilde{{\it\tau}}_{c}$ .

One may obtain a first-order approximation and circumvent the evaluation of the nonlinear equation (2.15) by using the parabolic shape of the velocity profiles of non-entraining flows even with entrainment or deposition. To this end, apply the boundary condition for the bed shear stress to the integrated form of (2.2) to get $w_{e}u(h)=hg_{x}-\hat{{\it\tau}}_{c}$ , and to (1.3). Together with the relation $u(h)=u(r)=2r\dot{{\it\gamma}}(0)$ for a parabolic velocity profile, one arrives at

The Boussinesq form factor (which appears in the momentum balance) changes only by a small amount as long as dynamically sustainable entrainment rates are specified, see the last column of table 1. However, the relative error in the entrainment rate and velocity increases rapidly with increasing entrainment rate, i.e. with decreasing ratio $({\it\tau}_{c}-{\it\tau}_{y})/(hg_{x})$ . In the case $\tilde{{\it\tau}}_{c}=0.7$ , $\tilde{{\it\tau}}_{y}=0.5$ (implying $\tilde{r}=0.5$ ), the approximation yields $\tilde{w}_{e}=6.0$ and ${\tilde{u}}(1)=0.1$ , whereas the self-consistent solution is $\tilde{w}_{e}=4.46$ and ${\tilde{u}}(1)=0.134$ .

Comparison of the predicted entrainment rates with experimental data is difficult because erosion rates have rarely been measured, and then mostly in dry-snow avalanches (Gauer & Issler Reference Gauer and Issler2004; Sovilla Reference Sovilla2004) for which a Bingham fluid is a poor approximation. For an order-of-magnitude check, table 2 lists entrainment speeds and rates that would result from a typical slow, humid granular snow avalanche on a 30° slope, with a flow depth of 1 m, a density of $200\ \text{kg}\ \text{m}^{-3}$ and a surface velocity of $20\ \text{m}\ \text{s}^{-1}$ . The values of ${\it\tau}_{c}$ were chosen in a realistic range while ${\it\tau}_{y}$ and ${\it\nu}$ result from the requirement $u(h)=20\ \text{m}\ \text{s}^{-1}$ .

Table 2.

Examples of entrainment rates obtained for different combinations of snow-cover strength ( ${\it\tau}_{c}$ ), yield strength ( ${\it\tau}_{y}$ ) and Bingham viscosity ( ${\it\nu}$ ) in the flow. Idealized stationary flow with density $200\ \text{kg}\ \text{m}^{-3}$ , flow height 1 m and surface velocity $20\ \text{m}\ \text{s}^{-1}$ on a 30° slope.

To conclude this section and illustrate the similarities and differences between deposition and entrainment, we modify our model, supplying material at the top at surface velocity and removing it at the bottom after it has come to rest. Deposition makes the non-dimensionalized shear stress profile (cf. figure 4 a) convex instead of concave in the shear layer; the value at the bed interface exceeds the gravitational traction. Similarly, the velocity normalized by its value at the surface is larger in the shear layer than in the flow without deposition, opposite to the effect of entrainment (see figure 4 c). In contrast to the eroding flow, the material properties do not impose a boundary condition for the bottom shear stress in the depositing flow. The rate at which mass is supplied to the flow at its surface speed determines the flow velocity according to (2.11), with negative $w_{d}$ instead of positive $w_{e}$ . We stress, however, that determination of the deposition rate in more realistic situations requires explicit modelling of the flow dynamics and consolidation processes at the bed–flow interface, which cannot be carried out strictly within the framework of depth-averaged equations.

Figure 4.

Non-dimensionalized profiles of (a) shear stress and (b) velocity for quasi-stationary flows of Newtonian and Bingham fluids on an inclined plane, with and without entrainment or deposition. The parameter values are $\tilde{r}=0.5$ for the Bingham fluid, $\tilde{{\it\tau}}_{c}=0.7$ for the entraining flows, and $\tilde{{\it\sigma}}_{xz}(0)=1.5$ for the depositing flows. (c) Velocity profiles for non-entraining, entraining and depositing Newtonian and Bingham fluids, scaled with their surface velocity. The same parameter values are used as in (a) and (b). (d) Dependence of the entrainment rate on the bed shear strength ${\it\tau}_{c}$ for various values of the yield strength ${\it\tau}_{y}$ of the Bingham fluid. All quantities are non-dimensionalized according to (2.13), with ${\it\sigma}_{g}\equiv {\it\rho}hg_{x}$ .

2.2. Eroding frictional–collisional fluid

Next we apply the general momentum-balance formula (2.2) for quasi-stationary entraining flows in a simple-shear configuration to the FC granular fluid described by (1.4). We obtain a nonlinear ordinary differential equation for the shear rate if we again assume the erosion speed $w_{e}$ to be known:

First, we divide by $2K\dot{{\it\gamma}}(z)$ and separate variables, then integrate one side over $z$ from $h$ to $z$ , the other side over $\dot{{\it\gamma}}$ from $\dot{{\it\gamma}}(h)=0$ to $\dot{{\it\gamma}}(z)$ to obtain

This becomes dimensionless,

with $\tilde{w}_{e}(1)=0$ at the flow surface, if we define

The solution of this transcendental equation cannot be expressed in terms of elementary functions, but with the help of Lambert’s $W$ function one obtains

after a few manipulations on (2.23). In the interval $[-1/\text{e},0[$ , where Lambert’s $W$ function is two-valued, the upper branch $W_{0}$ has to be selected. The exponent in the exponential function in the argument of $W_{0}$ ranges from $-\infty$ for small $\tilde{z}$ and large $q$ to $-1$ at $\tilde{z}=1$ , thus the argument of $W_{0}$ lies in the range $[-1/\text{e},0[$ , with the function value in the range $[-1,0[$ so that $\tilde{w}_{e}(\tilde{z})\in [0,1[$ .

Next we want to apply the boundary condition ${\it\sigma}_{xz}(0)={\it\tau}_{c}$ to determine $w_{e}$ from $\tilde{w}_{e}(0)$ , but the solution $\tilde{w}_{e}(0)=W_{0}(-\text{e}^{-1-q})$ given above is of no help because both $q$ and $\tilde{w}_{e}$ contain the unknown $w_{e}$ . However, inserting the expression for $\tilde{w}_{e}(0)$ from (2.24) into the definition of $q$ allows us to eliminate one of the instances of $\tilde{w}_{e}$ by transforming (2.23) at $\tilde{z}=0$ into

with

For lack of a better name, we will call $p$ the entrainment stress parameter as it is the ratio of the shear stress available for entrainment and the non-static shear stress; $p$ and all factors in $\tilde{w}_{c}$ except $w_{e}$ are given by the general set-up of the system ( $g_{x}$ , $g_{z}$ , $h$ ) and the rheology of the flow ( ${\it\mu}$ , $K$ ) or the bed ( $\hat{{\it\tau}}_{c}$ ). Thus the entrainment speed can be obtained once the solution to (2.26) has been found numerically.

Several important points can be read off (2.26) and (2.27) immediately. First, for a stationary flow to be possible, the down-slope gravitational force must be larger than the Coulomb friction, thus $g_{x}+{\it\mu}g_{z}>0$ is required (recall that $g_{z}<0$ ). Second, erosion is possible only if the down-slope gravitational bed shear stress is larger than the snow-cover strength, $g_{x}h>\hat{{\it\tau}}_{c}$ . This constrains $p$ to positive values. Third, if the Coulombic part of the shear stress alone were larger than the snow-cover strength, the snow cover would fail at once to a finite depth even if the flow is at rest. Thus we have to demand $\hat{{\it\tau}}_{c}+{\it\mu}g_{z}h>0$ , which in turn leads to $p<1$ . Incidentally, one arrives at the same conclusion by expanding (2.26) around $\tilde{w}_{c}=0$ . Fourth, if (2.26) is to be real, $\tilde{w}_{c}<1$ follows immediately. Fifth, for $p$ very close to 0, i.e. barely above the erosion threshold, we expand to first order in $p$ and to third order in $\tilde{w}_{c}$ and find that $\tilde{w}_{c}\approx 3p/2$ . Sixth, when $p$ approaches 1, $\tilde{w}_{c}$ is also driven towards 1, $p\tilde{w}_{c}^{2}+\tilde{w}_{c}\longrightarrow 1+p$ . Expanding (2.26) in ${\it\delta}=1-p$ and ${\it\epsilon}=1-\tilde{w}_{c}$ , one finds ${\it\epsilon}\approx \text{e}^{-1/(2{\it\delta})}$ .

When we plot the numerical solution of (2.26) as a function of $p$ , we find that it is surprisingly well approximated by the simple function $\sin ({\rm\pi}/2p)$ (figure 5). This appears to be a coincidence, and the agreement is not perfect, but the relative difference between the approximate and exact solutions is less than 6 % for small values of $p$ and diminishes for larger $p$ . An approximation of this quality may be useful in numerical flow models as it avoids solving the nonlinear equation millions of times: for a given value of $p$ one may use

Figure 5.

Dimensionless entrainment speed of a FC fluid in quasi-stationary flow as a function of the dimensionless parameter $p$ , (2.27). It is closely approximated by the function $\sin ({\rm\pi}p/2)$ .

The next steps are to find the shear stress and velocity profiles. The former is given by

Recall that $w_{e}$ is chosen such that $\dot{{\it\gamma}}(0)=\dot{{\it\gamma}}_{c}$ , thus $\hat{{\it\sigma}}_{xz}(0)=\hat{{\it\tau}}_{c}$ . Since $W_{0}(-1/\text{e})=-1$ , one confirms $\hat{{\it\sigma}}_{xz}(h)=0$ . To find the velocity profile, we have to evaluate the integral of (2.25), with an exponential function as the argument of $W_{0}$ . To this end, we substitute $w=W_{0}(-\text{exp}(-1-q+q\tilde{z}))$ , equivalent to $w\text{e}^{w}=-\text{exp}(-1-q+q\tilde{z})$ with $w$ running from $w_{1}(\tilde{z}):=W_{0}(-\text{exp}(-1-q+q\tilde{z}))$ to $w_{2}:=W_{0}(-\text{exp}(-1-q))$ . After switching the direction of integration, the integral becomes elementary:

One immediately confirms $u(0)=0$ . Ascertaining $u^{\prime }(h)=0$ directly from (2.30) is less straightforward, but can be accomplished by expanding the left- and right-hand sides of $\text{exp}(-1-q+q\tilde{z})\equiv x=w\text{e}^{w}$ around $\tilde{z}=1$ (or $x=-1/\text{e}$ ) and $w=-1$ , respectively, and using this in the total differential of (2.30). Because of $w_{1}(1)=W_{0}(-1/\text{e})=-1$ and (2.24) and (2.25), one obtains the surface velocity from (2.30) as

in agreement with what momentum conservation dictates. The same technique allows evaluation of the depth-averaged flow velocity, ${\bar{u}}$ , and the form factor (2.14), but the resulting expressions are lengthy and not particularly illuminating.

Table 3 gives examples of erosion rates and velocities as predicted by this theory. The slope angle ${\it\theta}$ , flow height $h$ , and rheological parameters ${\it\mu}$ and $K$ are typical of medium-size to large snow avalanches. Assuming an avalanche density of $250\ \text{kg}\ \text{m}^{-3}$ , ${\it\tau}_{c}$ values between 0.6 and 0.9 kPa result, which can also be considered realistic. Note that the calculated erosion speeds are the rates at which the avalanche flow depth increases (if the eroded snow is not redistributed). If the snow-cover density is $125\ \text{kg}\ \text{m}^{-3}$ , say, the snow cover is eroded twice as fast as indicated in the table. Cases 4, 6 and 7 are closest to the situation at the measurement location at the Vallée de la Sionne test site. With the chosen parameter values, the velocities are much lower than observed near the fluidized front, but may approximately match the velocity in the dense core of the avalanches in case 4. However, in the events we have profiling radar data of, the supply of erodible snow was already exhausted by the head of the avalanche and ${\it\tau}_{c}$ of the old snow was too large for erosion to occur.

Table 3.

Erosion speed ( $w_{e}$ ), mean and surface velocity ( ${\bar{u}}$ , $u_{s}$ ), and form factor ( $f_{B}$ ) for selected flow configurations of a quasi-stationary FC fluid on an inclined plane. Rows 6–8 ( ${\bar{u}}_{0}$ , $u_{s,0}$ ) represent non-eroding flows with the same parameters except $\hat{{\it\tau}}_{c}>hg_{x}$ .

For case 4 of table 3, figure 6 compares the resulting velocity profile (normalized by its value at the flow surface) to the well-known Bagnold profile $1-(1-\tilde{z})^{3/2}$ that would result if there were no entrainment. The difference is qualitatively similar to the one observed for the Bingham rheology. If the gravitational traction $hg_{x}$ exceeds the shear strength $\hat{{\it\tau}}_{c}$ only slightly, the mean and surface velocity are close to the ones obtained in a non-eroding flow with the same parameters (cases 6 and 7). If the difference between $hg_{x}$ and $\hat{{\it\tau}}_{c}$ becomes substantial, the erosion speed increases and the flow velocity decreases. The velocity profiles then approach straight lines except near the surface, where they curve very sharply to satisfy the condition $u^{\prime }(h)=0$ . This is also reflected in the form factor, which is $\frac{5}{4}$ for a Bagnold profile and $\frac{4}{3}$ for a linear profile. Cases 1–3 and 6, 7 have $f_{B}$ close to $\frac{4}{3}$ and $\frac{5}{4}$ , respectively. Since $\dot{{\it\gamma}}(\tilde{z})\leqslant \dot{{\it\gamma}}(0)=\dot{{\it\gamma}}_{c}$ , the mean velocity is limited to ${\bar{u}}_{max}=h\dot{{\it\gamma}}_{c}$ . As for the Bingham fluid, the shear stress profile is curved, starting from $\hat{{\it\tau}}_{c}$ at $z=0$ – in contrast to the linear profile in a non-entraining flow.

Figure 6.

Quasi-stationary flow of an eroding FC fluid: profiles of (a) the velocity normalized by its value at the surface and (b) the shear stress normalized by the gravitational traction $g_{x}h$ for parameter values of case 4 in table 3. For comparison, the Bagnold velocity profile and the shear stress in a non-entraining stationary flow (gravitational traction) are also shown. Note that the surface velocity of the entraining flow is only 60 % of the surface velocity of the non-entraining flow in this case.

3.1. General dynamical considerations

Gravity mass flows in nature start from rest, go through successive acceleration and deceleration phases due to the non-uniform topography, and finally come to rest again. Thus they are never in a truly (quasi-)stationary phase. Also, in reality no skimming occurs; instead, the flow depth grows, and eventually so does the velocity even on an infinitely long inclined plane. Therefore we need to address the question of how the exact solutions for quasi-stationary flows developed in § 2 can be adapted to non-stationary flow conditions. At best, one may hope to find exact solutions in very simple topographies, e.g. an infinitely long inclined plane (Eglit & Yakubenko Reference Eglit and Yakubenko2014). In cases of practical interest, i.e. with complex topography, some simplifying assumptions have to be made. However, they should be consistent with the physics and the assumptions used in deriving depth-averaged models. The uniqueness of the solutions achieved in the quasi-stationary case will be lost due to the approximations being non-unique.

In order to concentrate on the essentials, here we consider an infinitesimal column of variable mass and flow depth $h(t)$ moving along a path $Z(X)$ , parameterized by its arclength $s$ , at depth-averaged speed ${\bar{u}}(t)\equiv {\dot{s}}(t)$ . This allows us to drop terms with partial derivatives in the tangential plane, but we will account for the internal velocity profile. The flow depth may change due to erosion at the speed $w_{e}$ from the bed. In order to show the differences with the quasi-stationary system studied in § 2, we also allow skimming from the flow surface at the speed $w_{s}$ , but set $w_{s}=0$ afterwards. The equation of motion involves gravitational driving, the specific bed shear stress $\hat{{\it\sigma}}_{xz}\leqslant \hat{{\it\tau}}_{c}$ , and the decelerating effect of entraining mass from rest and losing mass moving at the speed $u(h)>{\bar{u}}$ :

Here, $\tilde{r}$ is the sheared fraction of the flow depth, see § 2.1. In the quasi-stationary case ( $qs$ ) treated previously, erosion at the bed was balanced by skimming at the surface, $w_{e}^{qs}=w_{s}^{qs}$ , thus $h=\text{const.}$ from (3.1). In contrast, the non-stationary case ( $ns$ ) is characterized by $w_{s}^{ns}=0$ . If entrainment occurs, the equation of motion becomes

In the stationary case, the problem consisted of finding $w_{e}^{qs}$ and $u(h)$ for a given value of their product. The non-stationary case raises the additional problem of apportioning the excess specific shear stress $hg_{x}-\hat{{\it\tau}}_{c}$ between entrainment ( $w_{e}^{ns}{\bar{u}}$ ) and acceleration ( $h\dot{{\bar{u}}}$ ).

Simply setting $w_{e}^{ns}({\bar{u}},{\it\theta},h,\hat{{\it\tau}}_{c})\approx w_{e}^{qs}({\it\theta},h,\hat{{\it\tau}}_{c})$ in (3.4) is a poor choice because the entrainment rate is determined by the stress gradient at the interface between bed and flow, which depends strongly on the flow velocity. In the following, we will discuss two possible extensions based on different approximations, one for the Bingham fluid and another for the FC fluid. While they cannot claim the same degree of rigour as the results obtained so far, they can be applied to real-world problems and eventually compared to numerical approaches that resolve the bed-normal dimension (Crosta et al. Reference Crosta, Imposimato and Roddeman2009b ; Issler & Pastor Pérez Reference Issler and Pastor Pérez2011; Eglit & Yakubenko Reference Eglit and Yakubenko2014).

For completeness, we state the full set of equations for a quasi-three-dimensional depth-averaged flow model (with $w_{s}=0$ ). The fields $h(\boldsymbol{x},t)$ and $\boldsymbol{u}(\boldsymbol{x},t)\equiv (u_{x}(x,y,t))$ , $u_{y}((x,y,t))$ replace $h(t)$ and ${\dot{s}}(t)$ . Neglecting curvature-dependent terms for simplicity, we obtain the mass balance and the momentum balance in conservative form,

The indices $i$ , $j$ , $k$ , $l$ take the values $x$ , $y$ , and summation over repeated indices is understood. In (3.6), the stresses are averaged over the flow depth except for the bed shear stress $\hat{{\it\sigma}}_{iz}^{(b)}$ , and the momentum flux term can be expressed in terms of ${\bar{u}}_{i}$ and the (unknown) form factor $f_{B}$ as $h\overline{u_{i}u_{j}}=f_{B}h{\bar{u}}^{2}$ . The term proportional to $w_{e,s}$ in (3.2), quantifying the deceleration of the flow due to the acceleration of the entrained mass from rest to the flow velocity, does not appear explicitly in (3.6); it is recovered, however, if one transforms the momentum balance equation into the equation of motion by substituting the mass balance equation into it. Also note that (2.2) should be extended to

(with summation over the indices $i,j\in \{1,2\}$ ) when normal-stress gradients occur.

3.2. Approximate solution for the Bingham fluid

A central ingredient in the present theory of erosion and entrainment is the assumed brittle behaviour of the bed material, which locks the bed shear stress at the value of the bed shear strength, ${\it\tau}_{c}$ . For a given rate-dependent rheology, this determines the shear rate at the bed–flow interface, $\dot{{\it\gamma}}_{c}$ , independent of the depth-averaged flow velocity, ${\bar{u}}$ . The derivative of the shear rate, however, is linked to the flow velocity, and at the same time it is a determining factor in the equation for the erosion rate, see (2.2). (Note that Issler & Pastor Pérez Reference Issler and Pastor Pérez2011 derived this equation for a non-stationary flow.) There is a one-to-one relation between $\dot{{\it\gamma}}^{\prime }(0)$ and ${\bar{u}}$ in the quasi-stationary case, but not in the non-stationary case.

In order to obtain an approximate erosion rate for non-stationary flows, we therefore need to make a reasonable assumption concerning the derivative of the shear rate at the bed–flow interface. To this end, at any given time we adopt the velocity profile of a quasi-stationary eroding ‘proxy’ flow with the same flow depth, shear-layer depth and critical shear rate. However, we use a different viscosity ${\it\nu}^{\prime }$ for the depth-averaged velocities to match between the original non-stationary flow and the proxy flow. In calculating the erosion speed, we use the original value of ${\it\nu}$ . Conceptually, this is analogous to the assumption, routinely made in depth-averaged flow models without entrainment, that the velocity profile is identical to that of an equilibrium flow with the same velocity.

The approximating profile function $v(z,t)$ has to satisfy four conditions: (i) no slip at the interface, $v(0,t)=0$ , (ii) vanishing shear in the plug layer, $\partial _{z}v(z,t)=0$ for $z\geqslant r$ , (iii) critical shear stress for erosion at the interface, $\partial _{z}v(0,t)=\dot{{\it\gamma}}_{c}$ , and (iv) matching depth-averaged velocities between the non-stationary flow and quasi-stationary proxy flow, $\bar{v}(t)={\bar{u}}(t)$ . The first and second conditions are automatically fulfilled by the solution (2.11) if $r$ or, equivalently, $\tilde{r}\equiv r/h$ match between the two flows. Applying the scaling of (2.13) with the viscosity ${\it\nu}^{\prime }$ , we obtain the depth-averaged velocity as

(We write ${\it\chi}$ instead of $\tilde{w}_{e}$ to emphasize that this equation does not give us the entrainment rate of the non-stationary flow.) Now we implement the condition (iii) stated above and find $1-\text{e}^{-\tilde{r}{\it\chi}}={\it\chi}{\it\nu}^{\prime }\dot{{\it\gamma}}_{c}/(hg_{x})$ , giving us ${\it\nu}^{\prime }$ in terms of ${\it\chi}$ . By substituting this in (3.8) and rearranging terms, we obtain the following equation for  ${\it\chi}$ :

Again, this is a transcendental equation – this time depending on the dimensionless parameter ${\it\alpha}\equiv {\bar{u}}/(h\dot{{\it\gamma}}_{c})$ – and needs to be solved numerically. We can substitute ${\it\chi}({\it\alpha})$ into the expression for $\dot{{\it\gamma}}^{\prime }(0)$ ,

and obtain the erosion rate by evaluating (2.2) at $z=0$ :

If the depth-averaged velocity is equal to that of the quasi-stationary flow, ${\bar{u}}={\bar{u}}_{qs}$ and ${\it\alpha}={\it\alpha}_{qs}$ , then ${\it\nu}^{\prime }={\it\nu}$ in (3.8) and ${\it\chi}({\it\alpha}_{qs})=\tilde{w}_{e}^{qs}$ . Recalling the scaling of $w_{e}$ and that ${\it\nu}\dot{{\it\gamma}}_{c}/(hg_{x})=\tilde{{\it\tau}}_{c}-\tilde{{\it\tau}}_{y}$ , one verifies that (3.11) recovers the erosion speed of the quasi-stationary flow, given by (2.15). This is exemplified in figure 7, where the corresponding solutions of the quasi-stationary equations are marked by symbols.

Figure 7.

The dimensionless erosion speed $w_{e}h/{\it\nu}$ of a Bingham fluid as a function of the dimensionless depth-averaged velocity ${\it\alpha}={\bar{u}}/(h\dot{{\it\gamma}}_{c})$ for two combinations of the dimensionless bed shear strength, $\tilde{{\it\tau}}_{c}=\hat{{\it\tau}}_{c}/(hg_{x})$ , and the relative shear-layer depth, $\tilde{r}=r/h$ . The points corresponding to the solutions of the corresponding quasi-stationary systems are marked with symbols. The curves end at the maximum value of ${\it\alpha}$ that is compatible with the chosen value of $\tilde{r}$ .

The maximum function appears in (3.11) because we demand the erosion rate to be positive for the reasons discussed in § 1.3. However, this condition is not equivalent to the condition ${\it\chi}>0$ , which separates velocity profiles that are straighter ( ${\it\chi}>0$ ) or more rounded ( ${\it\chi}<0$ ) than a parabolic profile, respectively. The threshold velocity can be determined from (3.11) and (3.8): some algebra on (3.11) shows the threshold value of $\tilde{w}_{e}$ to be $\tilde{w}_{e}^{thr}=-(\tilde{{\it\tau}}_{c}+\tilde{r}-1)^{-1}$ , which can be substituted into (3.9) to obtain  ${\bar{u}}_{thr}$ .

Figure 7 shows plots of the erosion speed (3.11), scaled with ${\it\nu}/h$ , versus the depth-averaged velocity (3.8), scaled with $h\dot{{\it\gamma}}_{c}$ , for two combinations of the shear-layer depth $r$ divided by $h$ and the bed shear strength divided by $hg_{x}$ . The graphs end at ${\bar{u}}/(h\dot{{\it\gamma}}_{c})=\tilde{r}(1-\tilde{r}/2)$ , which is the limit velocity for flows with convex velocity profiles, as we have assumed here. The corresponding limit value of the erosion speed is $w_{e}h/{\it\nu}=(\tilde{{\it\tau}}_{c}+\tilde{r}-1)^{-1}$ .

3.3. Approximate solution for the frictional–collisional fluid

For the FC rheology, we explore a different approach. Rather than matching profile functions of a quasi-stationary eroding flow to the depth-averaged velocity and critical shear rate of the actual non-stationary flow, we approximate the profile function directly. The numerical simulations of accelerating flows by Issler & Pastor Pérez (Reference Issler and Pastor Pérez2011) as well as the results of § 2.2 show that the profiles approach the straight line $u_{\infty }(z)=\dot{{\it\gamma}}_{c}z$ as the erosive power (and the velocity in the case of non-stationary flows) increases. The profiles become increasingly curved near $\tilde{z}=1$ . The two-parameter family of functions

captures this behaviour. $U$ , the surface velocity of the flow, and the exponent ${\it\beta}$ depend on variables of the actual non-stationary flow.

As in the case of the Bingham fluid, the first and second of the four conditions (i) $v(0,t)=0$ , (ii) $\partial _{z}v(h,t)=0$ , (iii) $\partial _{z}v(0,t)=\dot{{\it\gamma}}_{c}$ , and (iv) $\bar{v}(t)={\bar{u}}(t)$ are automatically fulfilled by the functional form (3.12) if ${\it\beta}>1$ . One easily calculates $\partial _{z}v|_{z=0}={\it\beta}U/h$ and ${\bar{u}}={\it\beta}/({\it\beta}+1)\,U$ . Substituting the latter relation into the former and equating $\partial _{z}v$ to $\dot{{\it\gamma}}_{c}$ , we obtain

The admissible range of ${\it\beta}$ leads to a restriction for the depth-averaged velocity,

the limit corresponding to a straight-line velocity profile. Values of ${\it\beta}$ larger than the Bagnoldian value of $\frac{3}{2}$ may arise, e.g. during start-up of the flow: the velocity profile evolves progressively from a plug-like shape with very high values of curvature near the bottom towards nearly linear profiles with the curvature concentrated at the top. Similarly, when a flow enters a steeper segment of the path, ${\it\beta}$ will initially increase rapidly, possibly beyond $\frac{3}{2}$ , and then diminish again as the flow adjusts to the new conditions.

The approximate erosion rate of the non-stationary flow follows from (2.2):

The scale for the entrainment speed is set by

but the actual entrainment speed is always smaller by a finite amount because ${\bar{u}}$ is limited. It grows as a nonlinear function of ${\bar{u}}$ from zero at the threshold velocity

to its maximum value $(g_{x}+{\it\mu}g_{z})/\dot{{\it\gamma}}_{c}$ at ${\bar{u}}=h\dot{{\it\gamma}}_{c}/2$ . Here ${\bar{u}}_{thr}$ may differ significantly from the value $\frac{2}{5}h\dot{{\it\gamma}}_{c}$ obtained by using the non-eroding equilibrium velocity profile. The exponent ${\it\beta}$ has been eliminated from (3.15), but can be recovered from (3.13) for calculating the form factor $f_{B}=({\it\beta}+1)/({\it\beta}+1/2)$ . Erosion is possible for $hg_{x}>\hat{{\it\tau}}_{c}-5K\dot{{\it\gamma}}_{c}^{2}$ if the velocity is high enough, but if the slope is too gentle, the flow is decelerated and erosion will stop eventually.

Figure 8 shows the non-dimensional approximate erosion speed $w_{e}/w_{e\infty }$ as a function of the non-dimensional depth-averaged velocity ${\bar{u}}/(h\dot{{\it\gamma}}_{c})$ for three different values of the non-dimensionalized threshold velocity $2{\bar{u}}_{thr}/(h\dot{{\it\gamma}}_{c})$ . Note that the curves do not pass precisely through the points corresponding to solutions of the quasi-stationary problem because the ansatz functions (3.12) do not capture the quasi-stationary profile exactly, but the agreement is satisfactory.

Figure 8.

Approximate velocity dependence of the erosion speed in the case of the FC rheology, for three different values of the threshold velocity. All quantities are scaled by $w_{e,\infty }$ and $h\dot{{\it\gamma}}_{c}$ , respectively. The symbols mark the quasi-stationary erosion rates of the corresponding example cases listed in table 3.

We may gain some insight into the properties of this formula by considering the asymptotic behaviour of the model consisting of (3.1), (3.2) and (3.15) in the case ${\it\mu}=0$ . (With ${\it\mu}>0$ , the model eventually breaks down as the flow depth increases, see § 4.2.) If there is a solution for which the erosion rate tends to a constant value, the flow depth tends to linear growth in time, $h\sim w_{e}t$ , and becomes very large. As a consequence, the bed resistance term in (3.2) becomes negligible and the equation of motion simplifies to

solved by ${\bar{u}}(t)\simeq g_{x}t/2$ , i.e. constant acceleration at half the value of the component of the gravitational acceleration along the slope (Issler & Pastor Pérez Reference Issler and Pastor Pérez2011). Eglit & Yakubenko (Reference Eglit and Yakubenko2014) showed that this asymptotic solution has a linear velocity profile and that the flow particles accelerate at the rate $g_{x}$ , i.e. at twice the rate of the centre of mass. The asymptotic velocity corresponds to the limit velocity for which the model ceases to be applicable (the limit itself growing linearly with the flow depth). Inserting the limits of (3.16) and (3.17), one sees that $w_{e}\rightarrow w_{e,\infty }$ as $h\rightarrow \infty$ and

With the values from table 1 in Issler & Pastor Pérez (Reference Issler and Pastor Pérez2011) and $g=10.0\ \text{m}\ \text{s}^{-2}$ , one finds $w_{e,\infty }=0.161\ \text{m}\ \text{s}^{-1}$ , which agrees with the result from the numerical simulation (figure 3e in Issler & Pastor Pérez Reference Issler and Pastor Pérez2011). Furthermore, the fact that the erosion speed tends to $w_{e,\infty }$ implies that the velocity profile becomes linear asymptotically.

We can compare the present model in detail with the published depth-resolved simulation by solving the mass-point equations (3.1) and (3.2) for the same parameter values, i.e. ${\it\mu}=0$ , $K=0.005\ \text{m}^{2}$ , initial depth (normal to the slope) $h_{0}=1.05\ \text{m}$ , ${\it\rho}_{f}=300\ \text{kg}\ \text{m}^{-3}$ , ${\it\tau}_{c}=1450\ \text{kPa}$ , and constant inclination angle ${\it\theta}=30^{\circ }$ . The actual model implementation differs somewhat from those equations (see appendix B for details), but this is of no consequence in the present context.

Figure 9(a) shows erosion to start at a lower velocity in the depth-averaged model than in the depth-resolved one. For this reason, the velocity grows more slowly in the initial phase than in the depth-resolved simulation. Asymptotically, however, the acceleration of the centre of mass tends to $(g\sin {\it\theta})/2$ in both models. The erosion speeds also tend to the same value, but the approach to the asymptotic value differs strongly between the models (figure 9 b): $w_{e}$ is a monotonically increasing function of time in the depth-resolved model. In the approximate analytic erosion model, it rapidly jumps to a value about twice as large as $w_{e,\infty }$ after the threshold speed is crossed, then drops approximately 20 % below $w_{e,\infty }$ , reaches another local maximum and finally approaches $w_{e,\infty }$ from above. Accordingly, the flow depth is approximately 25 % larger after a drop of 10 km (figure 9 a). This discrepancy between the models is closely related to the shape of the velocity profile, as figure 9(b) reveals: the depth-resolved model computes a profile close to the equilibrium Bagnold shape with $f_{B}$ slowly growing from $\frac{5}{4}$ to $\frac{4}{3}$ . In contrast, the depth-averaged flow has a blunter shape than the Bagnold profile $(\,f_{B}=\frac{5}{4})$ during the entire simulation, but will cross-over to an increasingly linear profile if given enough time. Why the approximate model chooses blunt profiles and how it can be modified to agree more closely with the depth-resolved model remains to be investigated.

Figure 9.

Avalanche flow on a 30° slope as simulated with a slab model for the FC rheology and the erosion formula (3.15), compared to a depth-resolved simulation presented by Issler & Pastor Pérez (Reference Issler and Pastor Pérez2011). See the text for the parameter values. (a) Velocity (left scale) and flow depth (right scale). (b) Erosion speed (relative to flow density, left scale) and shape factor $f_{B}$ (right scale) from the depth-averaged and depth-resolved models. Low values of $f_{B}$ correspond to blunt, plug-like velocity profiles.

4.1. Implications of the proposed erosion and entrainment theory

Perhaps the most fundamental result of the present work is the proof of concept that erosion and entrainment in depth-averaged models of GMFs can be formulated in a mathematically self-consistent and physically meaningful way. This is not a priori obvious because depth-averaging discards information about the state of the system at and near the bed–flow interface, and it is precisely there that the erosion rate is determined.

Two additional hypotheses helped to make the problem tractable: first, the mechanical problem becomes much simpler if we assume the bed to be perfectly brittle, with a characteristic shear strength ${\it\tau}_{c}$ . We argued that this is a reasonable approximation for snow. The central role of this hypothesis in the present approach suggests testing its applicability in other types of cohesive substrates overflowed by GMFs, like soil with a substantial content of silt and clay. In this context, excess pore pressure generated by the flow-induced shearing of the soil may be an important factor (Iverson Reference Iverson2012) that necessitates a different approach.

Secondly, the brittleness assumption suggests adopting an analogue of Owen’s (Reference Owen1964) hypothesis, stipulating that the interfacial shear stress is locked at the critical value ${\it\tau}_{c}$ whenever there is erosion. The validity of this hypothesis cannot be proved or refuted in the framework of depth-averaged models, but it might be possible to study this question through discrete-element simulations if the frictional granular bed can be modified so as to exhibit a finite yield strength at vanishing normal load (Rognon et al. Reference Rognon, Roux, Naaim and Chevoir2008).

Despite its unrealistic features, the toy model introduced in Issler & Jóhannesson (Reference Issler and Jóhannesson2006) and developed further in § 2 proved to be invaluable because it constrains the system to the point where it becomes rigorously solvable. It is an open question at present what other rheological models can be solved in this way, and it should be worthwhile to explore this direction further. A particularly interesting candidate rheology is the Jop–Forterre–Pouliquen model for granular materials (Jop, Forterre & Pouliquen Reference Jop, Forterre and Pouliquen2006).

The explicit analytic solutions (even though in terms of non-elementary functions) of the quasi-stationary problem were very helpful both in developing a numerical approach (Issler & Pastor Pérez Reference Issler and Pastor Pérez2011) and to find analytical extensions to non-stationary flows. Sections 3.2 and 3.3 illustrate that this process is not uniquely determined. It should be illuminating to approximately parameterize the numerical solutions for non-stationary flow of a Bingham fluid presented by Eglit & Yakubenko (Reference Eglit and Yakubenko2014) in a similar way as was done in § 3.3, or to adapt quasi-stationary eroding solutions of the FC model to the non-stationary case, as we did for the Bingham fluid in § 3.2.

Extensive comparison of the predictions of this theory to measurements is not at present possible. On the one hand, there do not seem to be any laboratory experiments studying entrainment from a brittle cohesive bed into a flow describable by the rheologies studied in this paper. On the other hand, full-scale measurements to date (mainly on snow avalanches) have not measured all the parameters that are needed to determine the flow rheology and the bed properties. The preliminary back-calculation (appendix B) of a small avalanche whose deposits were studied in great detail exemplifies this: at first sight, the model underestimates the erosion rate greatly. However, if our interpretation of the event (appendix A) is correct, most of the eroded mass was not really entrained into the flow and the simplistic slab model may in fact estimate the entrainment rate fairly well. However, this cannot be ascertained because the bed shear strength, which plays a central role in our model, was not measured and we had to guess it.

4.2. Validity limits of the theory

At several points in this paper, limitations of the applicability of the presented approach became apparent from physical considerations or mathematical analysis. These limitations can be traced back to the assumptions that erosion is gradual and that there is a critical shear strength of the snow pack. Together, they limit the shear rate to values below the critical value at the interface, $\dot{{\it\gamma}}_{c}$ . Since the shear rate decreases with increasing distance from the bed, the depth-averaged flow velocity is strictly limited.

Further restrictions arise because the shear stress exerted by an avalanche at rest must be less than ${\it\tau}_{c}$ in the model. For the Bingham fluid, this implies ${\it\tau}_{y}<{\it\tau}_{c}$ , but does not restrict the flow depth. The FC rheology contains a Coulomb friction term, which limits the flow depth to $h<{\it\tau}_{c}/({\it\mu}g\cos {\it\theta})$ . An eroding avalanche on an inclined plane will eventually exceed this limit, but the theory ceases to apply even before this point is reached: as $h$ grows, $\dot{{\it\gamma}}_{c}$ decreases towards 0, and so ${\bar{u}}_{max}=h\dot{{\it\gamma}}_{c}/2$ will at some point fall below ${\bar{u}}$ , which is increasing.

We qualitatively discuss what probably happens in an FC avalanche beyond this point; detailed study and mathematical modelling must be deferred to future work, however. Due to the growth of $h$ , $\dot{{\it\gamma}}_{c}$ decreases, but this is felt only near the bottom of the flow while the upper layers continue to accelerate. The decrease of $\dot{{\it\gamma}}_{c}$ directly leads to a pronounced increase of $w_{e}$ . The velocity profile is originally convex throughout the entire flow depth, but this cannot be maintained: it remains convex near the top, but becomes concave near the bottom. The family of profile functions (3.12) cannot accommodate this shape, hence the breakdown of the approximate model. It remains to be seen whether approximate profile functions can be constructed that capture this transition while still leading to tractable equations.

However, even if one succeeds in this regard, the system will reach the point where the Coulomb friction alone overwhelms the shear strength of the bed. The condition $\dot{{\it\gamma}}|_{z=0}=\dot{{\it\gamma}}_{c}$ cannot be maintained as $\dot{{\it\gamma}}|_{z=0}$ would become imaginary. The only possible reaction of the system seems to be catastrophic failure of the bed over a finite depth ${\rm\Delta}b$ down to a layer where the strength matches the load; this is the erosion mechanism termed ripping in § 1.1. One can probably consider this slab fully disintegrated if it is perfectly brittle as assumed hitherto, and treat it as an FC fluid. The shear stress at $z=0$ , which no longer is constrained at the value ${\it\tau}_{c}(0)$ , will progressively accelerate the eroded mass, the velocity profile being concave near $z=0$ initially. Over time, the shear stress at the new interface $z=-{\rm\Delta}b$ may grow to the value ${\it\tau}_{c}(-{\rm\Delta}b)$ as the velocity of the eroded slab and the shear rate increase, and then erosion by scour or by ripping resumes.

If future work confirms these conjectures, the failure of the model at ${\bar{u}}=h\dot{{\it\gamma}}_{c}/2$ is reminiscent of a second-order phase transition that is not captured by the model: all physical quantities like the velocity or the erosion rate evolve continuously, but the character of the profile of the shear rate changes discontinuously from monotonic to non-monotonic. In contrast, the breakdown at $\dot{{\it\gamma}}_{c}\rightarrow 0$ has the characteristics of a first-order phase transition where the erosion rate diverges and then drops to zero while the flow depth increases by a step. Thus the approximate non-stationary FC erosion model needs to be amended to (i) incorporate profiles with an inflection point and (ii) accommodate discontinuous changes in the flow depth.

4.3. Steps towards practical applications

In hazard mapping, the central task is to predict the run-out distance and other properties of (future) mass-flow events with long return periods. Practical application of the presented erosion models in such a context will therefore, among other parameters, require knowledge of the likely spatial distribution of brittle shear strength during the first few days after snow storms, when most avalanches are released. One has to expect that these values depend significantly on the snow climate, the altitude, the exposure to radiation and wind, and the return period. These dependences need to be investigated through extensive field measurements. For example, snow around the fjords in western Norway typically is warm and humid near sea level, but under exceptional conditions it may be cold and have low cohesion along the entire avalanche path. Under such conditions, avalanches have run more than a kilometre further than usual (Lied et al. Reference Lied, Instanes, Domaas, Harbitz and Hestnes1998). Also, the shear strength of snow increases significantly with depth in the snow pack. The complex relationships governing this effect need to be distilled into simple formulae so that the numerical model can compute the instantaneous shear strength ${\it\tau}_{c}(x,y,b(t))$ as erosion progresses.

Our preliminary attempt at back-calculating a real avalanche event with the proposed erosion formula for FC flows (appendix B) exemplifies that the inclusion of entrainment renders the system highly nonlinear and sensitive to the initial and boundary conditions. We therefore anticipate that – even if we had the perfect model available – it will take a substantial effort to establish practical procedures that allow predictive and reliable simulations to be carried out.

As mentioned in § 1.1, there are other important erosion mechanisms besides scour. The back-calculation of the small avalanche event in appendix B indicates that ripping may affect much larger masses than scour even though it had little influence on the flow dynamics in that case. Meaningful validation and application of the model will probably only be possible when all erosion modes are included in a more comprehensive model. That entrainment model should then be implemented, not in a slab model, but in a fully fledged quasi-three-dimensional code.

Finally, the present study shows how erosion and entrainment are inextricably linked to the flow rheology. This conclusion has immediate consequences for experimental studies of entrainment: experiments have to be devised such that the rheology of the flow can be determined simultaneously with the erosion rate and information on the erosion mechanism gained. This is a great challenge, but if progress can be made, it will give answers simultaneously to the two main problems in GMF dynamics.