2.1. Summary of model assumptions

Paper [IV] starts from the long-established depth-averaged balance equations of mass and momentum. For the earliest formulation in snow-avalanche dynamics, see (Eglit, Reference Eglit1968; Plam and others, Reference Plam, Radok and Taylor1984). Parker and others (Reference Parker, Fukushima and Pantin1986) give a detailed derivation in a two-layer situation similar to the one discussed in Section 4. In a coordinate system following the terrain, these equations read as follows if terms due to non-orthogonality and curvature (Gray and others, Reference Gray, Wieland and Hutter1999; Gray, Reference Gray2001; Bouchut and Westdickenberg, Reference Bouchut and Westdickenberg2004) are neglected:

(2) $$\partial_t h + {\bf \nabla} \cdot \lpar h \hskip 0.9pt\bar{{\hskip -0.9pt{\bi u}}}\rpar = Q,$$

(3) $$\eqalign{\partial_t \lpar h\hskip 0.9pt{\bar{\hskip -0.9pt{\bi u}}}\rpar + {\bf \nabla} \cdot \lpar h\hskip 0.9pt {\bar{\hskip -0.9pt{\bi u}}}\hskip 0.9pt {\bar{\hskip -0.9pt{\bi u}}}\rpar = & -\!{\bf \nabla} \left({1\over 2} g_{z} h^{2} \right) + h {\bi g} \cr & - {{\bar{\hskip -0.9pt{\bi u}}} \over \vert \hskip 0.9pt{\bar{\hskip -0.9pt{\bi u}}}\vert} \left[\mu h g_z + {g\over \xi} \hskip 0.9pt{\bar{\hskip -0.9pt{\bi u}}}^2 \right]\comma }$$

where $\bar {{\bi u}} = (u,\,v)^{\rm T}$ , ∇ = (∂ _x , ∂ _y ) and g  = (g _x , g _y )^T denote the depth-averaged mean flow velocity, the two-dimensional (2-D) gradient operator, and the slope-parallel components of the gravitational acceleration vector, respectively. Finally, h is the flow depth, and Q is the volumetric entrainment rate, which we will not discuss here. Incidentally, Eqn (3) assumes $\overline {\hskip -0.9pt{\bi uu}} \approx \hskip 0.9pt{\bar {\hskip -0.9pt{\bi u}}}\hskip 0.9pt{\bar {\hskip -0.9pt{\bi u}}}$ , i.e., a uniform velocity profile. This is a poor approximation if RKE is important, as the velocity profiles presented in [I] show.

The basic RKE model departs from standard Voellmy-type models such as RAMMS in that the friction parameters μ and g/ξ in Eqn (3) are postulated to depend on R _K , the RKE density due to fluctuations of the particle velocities about their mean values, as

(4) $$ \mu\lpar R_K\rpar = \mu_0 {\rm e}^{-R_K/R_0} \comma \quad {g\over \xi\lpar R_K\rpar } = {g\over \xi_0} {\rm e}^{-R_K/R_0} . $$

With R ₀, a new parameter enters the model that needs to be determined empirically. A further balance equation describing advection, production and dissipation of RKE complements the balance equations shown above:

(5) $$ \partial_t\lpar R_K h\rpar + {\bf \nabla} \cdot \lpar R_K h \hskip 0.9pt\bar{{\hskip -0.9pt{\bi u}}}\rpar = \alpha {\bi S}_{{\rm b}} \cdot\hskip 0.9pt \bar{{\hskip -0.9pt{\bi u}}} - \beta R_K h \comma $$

where S _b is the bed shear stress given by the second line in (3), and 0 < α < 1, β > 0 are two parameters. Crucial ingredients of the model are the production and dissipation terms for R _K as well as the postulated dependence of the Voellmy friction parameters (Eqn (4)).

While papers [I]–[IV] refer to the kinetic theory of granular flows, they depart from one of its well-established results without mentioning or justifying their choice: both for dilute and dense granular flows (and also for turbulent fluid flows), the dissipation rate is found to grow as R _K ^3/2 (Jenkins and Richman, Reference Jenkins and Richman1985; Jenkins and Berzi, Reference Jenkins and Berzi2010). This is known as Haff's law and is a basic result of kinetic theory (Haff, Reference Haff1983). If dissipation is assumed to grow only linearly with R _K as in Eqn (5), much higher equilibrium values of R _K result for a given production rate. Equation (5) contains another strong assumption for which no justification is given, namely that a fixed fraction of the shear dissipation rate is converted into RKE.

2.2. Velocity dependence of the effective friction coefficient in the RKE model

From Eqns (4) and (5), one can deduce the speed of steady, uniform flow for a given slope angle θ and flow depth h, and thus the effective friction law: the left-hand side of Eqn (5) vanishes in steady, uniform flow so that

(6) $$ R_K^\infty = {\alpha\over \beta} {{\bi S}_{{\rm b}} \cdot \hskip 0.9pt\bar{{\hskip -0.9pt{\bi u}}}\over h} = {\alpha\rho\over \beta h} \left[\mu_0 g_z h + {g\over \xi_0}\hskip 0.9pt \bar{{\hskip -0.9pt{\bi u}}}^2 \right]\vert \hskip 0.9pt\bar{{\hskip -0.9pt{\bi u}}}\vert {\rm e}^{-R_K^\infty / R_0} . $$

For simplicity, let us define $U \equiv \vert\bar {{\bi u}}\vert$ , the non-dimensional RKE r ≡ R _K ^∞/R ₀, the velocity scale

(7) $$ U_0 \equiv {\beta R_0\over \alpha \mu_0 \rho g_z} \comma $$

and the effective Voellmy friction coefficient,

(8) $$ \mu_{\rm eff}^{\rm V} \equiv c_{\rm V}\mu_0 \equiv \left(1 + {g U^2 \over \xi_0 \mu_0 g_z h} \right)\mu_0. $$

This brings Eqn (6) into the form

$$ r \lpar U\comma h\rpar = c_{\rm V} {U\over U_0} \, {\rm e}^{-r} \comma $$

which is solved by

$$ r\lpar U\comma h\rpar = W_0\lpar c_{\rm V} U / U_0\rpar . $$

W₀ is the upper branch, defined on (−e⁻¹, ∞), of Lambert's W function, which is the solution to x = W(x) exp (W(x)). As x → ∞, also W(x) → ∞, but for x > e, W(x) < ln x. Applying this to Eqn (4) and carrying out a few algebraic manipulations involving the defining equation of W(x), we can rewrite the effective friction coefficient of the RKE model in steady, uniform flow as

(9) $$ \mu_{\rm eff}^{\rm RKE}\lpar U\comma h\rpar = \mu_0 {U_0\over U}\, W_0\lpar c_V U / U_0\rpar . $$

2.3. Comparison with kinetic theory

It is interesting to compare this heuristic bed-friction law to one derived using methods of the kinetic theory for collisional grain flows. We consider an essentially passive slab that is suspended and transported on a relatively thin region of intensely sheared grains at its base (Jenkins and Askari, Reference Jenkins and Askari1999). This may represent a slab avalanche in a very early phase, when the slab is only about to disintegrate and is gliding on the thin weak layer whose collapse caused the avalanche to release. We do not propose this model as a replacement for any other model, but have chosen it because it describes the same flow configuration as the Voellmy model (Salm, Reference Salm1993), namely a deformable but essentially passive heap riding on a thin, intensely sheared basal layer. Yet, the Jenkins–Askari model shows distinctly different behaviour from both the original Voellmy model or the basic RKE model of [IV].

Expressions for the shear stress, normal stress, energy flux and rate of collisional dissipation at the base result from detailed consideration of the transfer of momentum and energy in particle collisions with a bumpy, rigid boundary (Richman, Reference Richman1988). Their values are obtained as solutions of boundary value problems for the RKE and average particle velocity in the thin region of intense shearing at the base. The resulting expression for the dynamic, or rate-dependent part, of the ratio of the shear stress, S, and normal stress, P, at the base, when added to the rate-independent part, μ ₀, gives the effective friction coefficient (Jenkins and Askari, Reference Jenkins and Askari1999):

(10) $$\eqalign{\mu_{\rm eff}^{\rm JA} & = \mu_0 \cr & \quad - \alpha \left({4g_z h\over U^2}\right)^{1/2} + \left[\alpha^2{4g_z h\over U^2} + {12+\pi\over 5\pi}\lpar 1-e\rpar \right]^{1/2}.}$$

Here, α measures the difference between slip working and collisional dissipation at the boundaries; it is a function of the boundary roughness and the coefficient of restitution, e _w, in a collision between a flow sphere and the boundary. e is the coefficient of restitution in a collision between two particles in the flow. Note that the model is restricted to low and moderate velocities so that the slab on top of the thin shear layer does not become strongly agitated.

Figure 1 compares the velocity dependence of the effective friction coefficient of the original Voellmy friction law (in two different calibrations), the kinetic theory and the RKE-enhanced Voellmy model. Note that these curves are valid only for the specific choice of parameters. For all friction laws, we chose a slope angle of 30° and a flow depth of 1 m. For the Voellmy model, we set μ = 0.15 and g/ξ = 0.0049 in the usual calibration, but μ = 0.4 and g/ξ = 0.00018 in the calibration proposed by Gauer (Reference Gauer2014). For the Jenkins–Askari model, we selected μ ₀ = 0.4, volume fraction at the base ν ₀ = 0.45, and restitution coefficients e = 0.85 inside the shear layer and e _w  = 0.80 at its boundaries. Assuming a rough boundary of smaller particles, we set α = 0.463. For the RKE model, μ ₀ = 0.4, g/ξ ₀ = 0.017, α = 0.1, β = 0.5 s⁻¹ and R ₀/ρ = 15 m² s⁻² were assumed. With the traditional calibration of the Voellmy model, the resistance is dominated by the velocity-dependent term already at speeds of 15–20 m s⁻¹ so that the model generally underpredicts avalanche velocities substantially. This defect is mitigated in Gauer's calibration – at the expense of an almost linear dependence of μ on the mean slope angle of the path and a linear dependence of ξ on the drop height (Gauer, Reference Gauer2014). The Jenkins–Askari model exhibits rapid growth of μ _eff at low velocity, but moderate growth at high velocity, thus allowing the avalanche to reach high velocity even on moderately steep slopes. Finally, the RKE model shows completely different, non-monotonic behaviour, with μ _eff decreasing with increasing velocity above some threshold speed. We caution, however, that Fig. 1 takes the models beyond their range of applicability. Within the realistic range 0–60 m s⁻¹, a constant value μ ≈ 0.42 and g/ξ = 0 – corresponding to the Savage–Hutter model (Savage and Hutter, Reference Savage and Hutter1989, Reference Savage and Hutter1991) – would give fairly similar behaviour.

Fig. 1.

Comparison of the velocity dependence of the effective friction coefficient for an avalanche with flow depth 1 m on a 30° slope. (i) Voellmy model with traditional calibration (Christen and others, Reference Christen, Kowalski and Bartelt2010), μ =  0.15, g/ξ = 0.0049. (ii) Voellmy model with the calibration suggested by Gauer (Reference Gauer2014), for an average slope angle β ≈ 30°: μ ₀ = 0.4, g/ξ = 0.00018. (iii) Jenkins and Askari (Reference Jenkins and Askari1999) model with μ ₀ = 0.4, ν ₀ = 0.45, e = 0.85, e _w  = 0.8 and α = 0.463. (iv) RKE model [IV] with μ ₀ = 0.4, g/ξ ₀ = 0.017, α = 0.1, β = 0.5 s⁻¹ and R₀/ρ = 15 m² s⁻².

The Jenkins–Askari model predicts that stationary flows are only possible in a limited range of slope angles (22°–38° for the choice of parameters in Fig. 1). The original Voellmy law (1) clearly has a lower bound $\theta \geq \arctan \mu $ , but no upper bound because μ _eff grows as the square of the velocity. A stationary flow can thus be attained in arbitrarily steep terrain. From Figure 1, one can infer the following scenario in the RKE model: an avalanche starting in sufficiently steep terrain will quickly attain high velocity and low effective friction. If the slope angle decreases in such a way that always tanθ(x) ≳ μ _eff ^RKE(U(x)), the avalanche will continuously accelerate while the slope angle tends to 0. Of course, slopes on Earth are too short for this to happen, but the theoretical possibility illustrates that Eqn (4) is an extremely strong assumption.

It has been known for a long time that the shear stress in rapidly sheared granular materials increases with the square of the shear rate if the volume is held constant. If the experiment is carried out at constant normal stress (as in free-surface chute flows), the shear stress increases also, but much more slowly due to the material expanding. For example, experiments on a wide range of slope angles (Holyoake and McElwaine, Reference Holyoake and McElwaine2012) show that the effective friction always increases with velocity for a fixed flow rate. Kinetic theory successfully predicts this behaviour.

It is conceivable that one sometimes can obtain satisfactory simulations of observed events if completely different mechanisms have a similar friction-reducing effect, e.g., lubrication by a thin water layer, progressive comminution of snow particles or excess pore pressure. However, if the objective is to construct a physically founded model of snow avalanche motion, it appears more promising to adopt the key results from the theory of granular flows and to add such non-granular effects in a more specific manner.

2.4. Comparison with experiments

Paper [IV] points out that chute experiments with snow (Platzer and others, Reference Platzer, Bartelt and Kern2007) show significantly smaller ratios of basal shear stress S to basal normal stress N in the head of the flows than at the snout and in the tail. The authors attribute this difference to different levels of RKE and argue that, therefore, a separate balance equation for RKE is needed. However, this conclusion does not follow: the entire flow mass started from rest simultaneously, so its flow depth, velocity and RKE should evolve in the same way unless some other process differentiates between head and tail. Longitudinal normal stresses reinforce gravity at the snout and counteract it in the tail, leading to a significant decrease of the velocity from snout to tail (Platzer and others, Reference Platzer, Bartelt and Kern2007). Thus, a shear-rate-dependent rheology may also be able to capture the observed variation in the ratio of shear stress to normal stress.

There is, however, a finding from the chute experiments (Platzer and others, Reference Platzer, Bartelt and Kern2007) that deserves closer consideration: independent of whether the flow is stationary or not, S/N is the effective bed friction coefficient, μ _eff. The snow chute experiments indicate that μ _eff is higher in the slower tail than in the faster head, while the snout, which moves at the same speed as the head, also has a high value of μ _eff. This appears to be in stark contrast to the findings of the granular–flow experiments and the granular rheology, in which μ _eff increases monotonically with the velocity. It should be illuminating to reanalyse the snow-chute experiments in terms of granular rheology in order to confirm or refute the discrepancy with the granular experiments. One may then try to relate the variation of μ _eff to, e.g., the dominant particle size at different locations in the flow.

2.5. Is a dynamical equation for RKE necessary?

Another interesting question is whether the extra balance equation for RKE is really necessary. Certainly, it is needed in a complete theory of granular flows because RKE (or granular temperature) may be produced at one location, transported by advection and diffusion, and finally dissipated somewhere else. Bartelt and co-workers emphasise that RKE is produced at the bed interface and diffuses into the avalanche body. Much of the RKE is certainly produced near the bed, where the shear rate is highest, but in highly agitated flows with a Bagnold-type velocity profile, significant RKE production occurs also inside the flow so that bed-normal diffusion of RKE need not play a dominant role in the balance equation. In fact, kinetic theory shows that boundary effects decay exponentially, with a decay length that is only a few particles long for dense inelastic flows. This means the equilibrium granular temperature is slaved to the local shear rate with an extremely rapid relaxation time. Moreover, the flow model is depth-averaged and bed-normal RKE diffusion is not directly visible in the model; this makes it much less compelling to use an extra equation for RKE unless the dissipation coefficient β is very small, which hardly can be the case in a dense granular flow with rather inelastic collisions. Typical values of the inverse RKE decay constant 1/β ≳ 1 s used in [IV] are short compared with the macroscopic timescale of avalanche flow, i.e., the RKE is almost always close to its instantaneous equilibrium value. In addition, it is inconsistent to include a process with such rapid relaxation while excluding processes with much slower relaxation.

The RKE model, using the Voellmy friction law as a basis, inherits from it the general disadvantage that it only specifies the basal shear stress. If the shear stress inside the flowing material or the normal stresses are needed, then additional assumptions have to be made. Bartelt and others (Reference Bartelt, Buser and Platzer2006) assume the Voellmy model for the basal shear stress, a viscous-frictional model inside the flow and an exponential decrease of the RKE with distance from the bed in order to approximate velocity profiles measured at the test site, Vallée de la Sionne. In particular, these rheological assumptions determine the bed-normal stresses inside the flow. This approach was, however, abandoned in subsequent papers.

3.1. The approach of [V]: analogy with an ideal gas.

In the text following (V.6), Buser and Bartelt postulate that the equation of state of a granular snow avalanche is essentially equivalent to that of an ideal gas, i.e., they assume

(11) $$ N^{\lpar b\rpar } h = \gamma R_K h $$

in their notation. Next, they take a total time derivative of the equation of motion (V.4) to arrive at (V.5) and then do the same with the equation of state shown above. However, in doing so, they only keep $\gamma\hskip 0.9pt \dot{\hskip -0.9pt R}_K h$ on the right-hand side and omit $\gamma R_K {\dot h}$ :

(V.6) $$ {{\rm d} \lpar N^{(b)} h\rpar \over {\rm d} t} = \gamma \hskip 0.9pt \dot{\hskip -0.9pt R}_K h. $$

With ${\hskip2pt}\dot {{\hskip-2pt}N}_g^{(b)} \equiv 0$ and ${\hskip2pt}\dot {{\hskip-2pt}N}_K^{(b)} \equiv {\dot {{\hskip-2pt}N}}^{(b)}$ , this mathematical error leads to the third-order equation adopted in [V],

(V.7) $${\ddot k} + \left( g_z + {\ddot k}\right) \displaystyle{{\dot k}\over {k}} = \displaystyle{{\gamma \hskip 0.9pt \dot{\hskip -0.9pt R}_K} \over {m}}.$$

Indeed, if one used Eqn (11) in (V.4) and set N _K ^(b) = N ^(b) + mg _z (note that g _z  < 0), one would immediately arrive at the second-order equation

(12) $$ {\ddot k} = \gamma {R_K\over m} + g_z . $$

As the avalanche expands perpendicularly to the bed, the RKE density R _K tends to decrease both because the energy is distributed over a larger volume and energy is expended in working against gravity.

A more detailed analysis of the analogy between thermodynamics and the kinetic theory of granular materials suggests to replace the constant γ by a function of the particle volume fraction, f(ν). This function depends on the details of the particular approximation one chooses. Using ν = ν ₀ k ₀/k, we can express f(ν) in terms of k as $\tilde{\hskip -0.9pt f}(k) = f(\nu _0 k_0 / k)$ and replace γ with $\tilde{\hskip -0.9pt f}(k)$ in Eqn (12).

3.2. The approach of [VI]: configuration energy

Paper [VI] arrives at Eqn (VI.25), which is essentially equivalent to (V.7), without invoking analogy with an ideal gas or making an explicit assumption for the constitutive equation; hence, we need to discuss this approach separately. As mentioned above, an important difference between [V] and [VI] is the use of energy balances throughout [VI]. In particular, the notion of configuration energy density (CED), R _V , is introduced. In [VI], Section 3, it is defined as the gravitational potential energy per unit volume, averaged over the depth of a bed-normal column, relative to the configuration with maximum random packing at the same location. Buser and Bartelt introduce the production rate of CED, $\hskip 0.9pt \dot{\hskip -0.9pt P}_V$ , and postulate it to be a fixed fraction γ of the net production rate of the sum of RKE and CED.

The critical step of their derivation is described between Eqns (VI.20) and (VI.21), which we quote here (note that the authors switched notation from h to V _Φ, but we will use h in what follows):

The total work done per unit time by the normal pressure at the bottom of the avalanche N […] must be in balance with the total working of the particle interactions per unit volume. We have termed this change in potential energy as the configurational energy production $\hskip 0.9pt \dot{\hskip -0.9pt P}_V$ . Therefore, the total change in the volume is

(VI.21) $$ \partial_t \lpar N V_\Phi\rpar = \hskip 0.9pt \dot{\hskip -0.9pt P}_V V_\Phi. $$

Equation (VI.21) above is identical with (V.6) if one identifies the coefficients γ in [V] and [VI] with each other, extends R _K in (V.6) to R = R _K  + R _V  = R _K /(1 − γ), and replaces the partial time derivative ∂ _t with the advective derivative D _t  ≡ ∂_t  +u⋅∇ to account for the translational motion of the control volume. We find this text passage somewhat ambiguous, but understand it as making the following interrelated statements: (i) At the mesoscopic level, particle interactions in the avalanche do mechanical work, which manifests itself as (and is quantitatively equal to) the mechanical work done by the basal pressure at the macroscopic level. (ii) The rate of mechanical work done by the pressure due to particle collisions must be equal to the rate of change of the configuration energy. (iii) The work rate of the bed pressure is D _t (Nh). The paragraphs below will analyse these three statements in detail. Note that we will henceforth write P _V instead of $\dot {P}_V$ used by Bartelt and co-workers. Their notation suggests that the production rate is the time derivative of some other quantity, which, however, is never introduced and would have no other meaning than the total of the RKE produced in an advected unit volume – in particular, it is not the RKE.

3.3. Pressure, work rate and configuration energy

The first issue to note is that, contrary to the statement (i) above, the pressure at the bottom does not do mechanical work because the bed-normal velocity w(0) vanishes at that boundary. At the top surface, w(h) ≠ 0, but the pressure (relative to atmospheric pressure) vanishes. The mechanical work is being done inside the mass, where w(z) ≠ 0 and N(z) ≠ 0. This is not merely a semantic point because it leads to extra coefficients in the expression for the work rate that are missing on the left-hand side of Eqn (VI.21). We will detail this after mentioning the remaining issues.

Second, statement (ii) neglects the change of kinetic energy associated with bed-normal motion that unsteady expansion necessarily induces. A third issue is intertwined with this one: the text as well as Eqn (VI.21) set the change in configuration potential energy equal to the work done by the total pressure N = N _g  + N _K (statements (ii) and (iii) combined). However, only N _g (z) contributes to the change of CED, whereas the gradient of dispersive pressure, ∂_z N _K , accelerates the avalanching mass in the z-direction and changes the corresponding contribution to the kinetic energy density ${\cal K}_z \equiv \overline {\rho w^2}/(2\,h)$ .

Fourth, Eqn (VI.21) stipulates that the work rate of the pressure is D _t (Nh); no further explanation for this assertion is given. The expression D _t (Nh) contains the three terms $\dot {N}_K h$ , $N_g {\dot h}$ and $N_K {\dot h}$ ( $\hskip 1pt \dot{\hskip -1pt N}_g h$ is zero if there is no net mass flux into or out of the control volume). If Eqn (VI.21) were true, pressure would do mechanical work whenever it increases, even if h is held constant. A simple example is heating of a gas in a rigid container: the gas does not do mechanical work in this process, but its capacity to do so increases.

In Appendix B, we will compare the model [VI] to the general balance equations for mass, momentum and fluctuation energy in a depth-averaged flow model. It may be instructive, however, first to illustrate the issues mentioned above in a simple, quasi-1-D setting, disregarding variations in the x- and y-directions. We assume, however, that shearing motion in these directions produces RKE. The control volume under consideration and the stresses and body forces in the z-direction are schematically represented in Fig. 2. We first consider the momentum balance and then turn attention to the energy balance.

Fig. 2.

Schematic representation of an infinitesimally thin column of an avalanche on a plane inclined at an angle θ and an infinitesimal control volume within that column. Only the forces relevant for the bed-normal motion are indicated.

The momentum balance for a thin slope-parallel slice of thickness dz at height z is given by

(13) $$ \partial_t \lpar \rho w\rpar + \partial_z \lpar \rho w^2\rpar = \rho g_z - \partial_z N. $$

Integrating this equation over z from 0 to ∞, using the boundary conditions ρ(∞) = w(0) = N(∞) = 0 and introducing the mass hold-up $m \equiv \int \rho (z) {\rm d} z$ as well as the centre-of-mass velocity $\bar {w} \equiv \int \rho w {\rm d} z$ , one obtains

(14) $$ m {\dot{\bar w}} = mg_z + N^{(b)}\comma $$

where N ^(b) is the value of N at the bed. If one assumes the density to be uniform across the flow depth but variable in time, $w(z,\,t) = 2 \bar {w}(t) z/h$ for kinematic reasons. Similarly, $ \hskip 0.9pt\dot{\hskip -0.9pt w}(z,t) \equiv (\partial_t + w(z ,t ) \partial_z)w(z , t)= 2 {\dot {\bar w}}(t) z/h$ in the Lagrangean sense. The dynamics, governed by Eqn (13), then requires the pressure to vary with ζ ≡ z/h as

(15) $$ N(z\comma t) = m \left[-\lpar 1-\zeta\rpar g_z + \lpar 1 - \zeta^2\rpar {\dot{\bar w}}\lpar t\rpar \right]. $$

The kinetic and configuration energies are given by

$$ \bar{{\cal K}}h + R_V h \equiv h \int_0^1 {\rho\over 2} w^2 {\rm d}\zeta - \int_0^1 \rho g_z h\zeta {\rm d}\zeta + m k_0 g_z . $$

If the density profile is uniform, these energies become

(16) $$ \bar{{\cal K}}h = {4\over 3} \times {1\over 2} m \hskip 0.9pt\bar{\hskip -0.9pt w}^2 \quad {\rm and} \quad R_V h = -m \left({h\over 2} - k_0\right)g_z. $$

We can calculate the depth-integrated work rate $ \hskip 0.9pt\dot{\hskip -0.9pt W}$ from the well-known expression for the work rate per unit volume, $\hskip 1pt\dot{\hskip -2.4pt{\cal W}}$ , of the internal (Cauchy) stresses in a continuous medium:

(17) $$\hskip 1pt\dot{\hskip -2.4pt{\cal W}} = \sigma_{ij} D_{ij} \equiv \sigma_{ij} {1\over 2}\lpar \partial_j u_i + \partial_i u_j\rpar.$$

In our case and with the assumption that the density profile is uniform, the strain rate has only one non-zero component, given by $D_{zz} = \partial _z w = 2\bar {w}/h$ ; the stress tensor has the corresponding component σ _zz (z, t) = − N(z, t). A straightforward calculation of $ \hskip 0.9pt\dot{\hskip -0.9pt W}$ using Eqn (15) then yields

(18) $$ \hskip 0.9pt\dot{\hskip -0.9pt W} = {4 \over 3} m{\bar w}{\dot{\bar w}} - m {\bar w} g_z . $$

Comparison with Eqn (16) immediately shows $ \hskip 0.9pt\dot{\hskip -0.9pt W} = {\rm D}_t[({\bar {\cal K}} + R_V) h]$ , as it should be.

Now, if D _t (N ^(b) h) indeed is the work rate of pressure, directly evaluating it must give the same result as in Eqn (18). From Eqn (15) and the assumed uniform density profile, we obtain

(19) $$ {\rm D}_t \left[N^{\lpar b\rpar } h \right] = 2\, m \left[{\rm D}_t\left( k{\ddot k}\right) - \bar{w}g_z \right]. $$

The right-hand sides of Eqns (18) and (19) differ significantly. Equation (18) is derived directly from the principles of continuum mechanics, with the only additional assumption of w(z, t) being linear in z. Equation (19) follows from the postulate (VI.26) combined with the definition (VI.3), the balance Eqn (VI.8), the postulate (VI.12) and assumed linearity of w(z, t). This leaves two alternatives: either D _t (N ^(b) h) is not a correct expression for the depth-integrated work rate of the pressure, or Eqn (VI.21) must be considered an implicit constitutive assumption for the granular pressure. The first alternative has far-reaching consequences: (VI.21) must be abandoned and the basis for the mathematical development in the rest of Sections 5 and 6 of [VI] is invalidated. Instead, one has to adopt Eqn (18) and state a suitable constitutive equation for the granular pressure as a function of the flow variables, e.g., N = N(R _K ).

Now consider the second alternative: in this case, D _t (N ^(b) h) is not the total work rate of the granular pressure. In order for D _t (N ^(b) h) to be equal to D _t (R _V h), the equation ${\rm D}_t(N^{(b)}h) = {\dot W} \,-\, {\rm D}_t({\bar {\cal K}}h)$ must hold. In the case of linear w(z, t) with h ≡ 2k and $\bar {w} \equiv {\dot k}$ , this leads to

$$ {\rm D}_t\lpar N^{\lpar b \rpar} h \rpar \equiv 2{\rm D}_t\lpar N^{\lpar b\rpar} k \rpar = -m {\dot k} g_z. $$

If m and g_z are constant along the avalanche path, then we can perform the time integration (more precisely, the integration along the characteristic line of the control volume) easily and obtain

$$ N^{\lpar b\rpar }\lpar t\rpar = -{1\over 2} m g_z + {\rm cst.} $$

If the avalanche is at rest (and its depth is 2k ₀), N ^(b) must equal the weight, thus cst. = − (1/2)mg_z . However, this leaves no room for a dynamical evolution of the flow depth, and we conclude that the second alternative is not viable.

In our opinion, the most immediate solution of this dilemma in the energy formalism is to explicitly state a constitutive equation for N ^(b), to express the work rate of the granular pressure in terms of N ^(b) as

(20) $$\eqalign{ \hskip 0.9pt\dot{\hskip -2pt W} &= 2{\bar w} \int_0^1 \left[N_g^{\lpar b\rpar } \lpar 1-\zeta\rpar + \lpar N^{\lpar b\rpar } + m g_z\rpar \lpar 1-\zeta^2\rpar \right]{\rm d}\zeta \cr &= \left[-m g_z + {4\over 3} \lpar N^{\lpar b\rpar } + m g_z\rpar \right]{\bar w} }$$

and to set it equal to the rate of change of kinetic and configuration energy. Using Eqn (16), dividing by $\bar {w}$ and rearranging terms, we arrive at

(21) $$ m{\dot {\bar w}} = N^{\lpar b\rpar } + m g_z . $$

Not surprisingly, this is the same as Eqn (12) if one makes the constitutive assumption N ^(b) = γR _K . We note that such modification of the extended RKE model also requires modifying the RKE balance equation to properly account for the conversion of RKE to kinetic and configuration energy:

(22) $$\eqalign{\partial_t\lpar R_K h\rpar & + {\bf \nabla} \cdot \lpar R_K h \hskip 0.9pt{\bar{\hskip -0.9pt{\bi u}}} \rpar \cr &= \alpha \hskip 0.9pt\dot{\hskip -0.9pt W}_f - \beta_K R_K h - \left({4\over 3}\gamma R_K + {1\over 3}m g_z\right)w.}$$

This will be discussed further when comparing the extended RKE model with the general form of the balance equations for mass, momentum and RKE in Appendix B.

3.4. First-order equation for the flow depth implied by the model assumptions

Our fifth remark is that the modelling assumptions put forth in the text of [VI] together with the energy partitioning postulate (VI.12) imply a simple, first-order evolution equation for h. As the cited text from [VI] states, Eqn (VI.21) is to describe the (advected) rate of change of CED, ${\rm D}_t R_V \equiv \partial _t R_V +\hskip 0.9pt {\bar {\hskip -0.9pt{\bi u}}}\,\cdot \,{\bf \nabla} R_V$ . Its right-hand side is the production rate of CED, P _V , integrated over the flow depth. Buser and Bartelt postulate the following equations:

(VI.12) $$ P_V = \gamma P,$$

(VI.8) $$P h = \alpha \hskip 0.9pt\dot{\hskip -0.9pt W}_f^{xy} - \beta_K R_K h,$$

(VI.6) $$ \hskip 0.9pt\dot{\hskip -0.9pt W}_f^{xy} = {\bi S}_{\rm b} \cdot \hskip 0.9pt{\bar{\hskip -0.9pt{\bi u}}}_\parallel.$$

In an avalanche starting at rest, R _K (0) = R _V (0) = 0. Due to(VI.12), P _K  = (1 − γ)P so that the advected rates of change of R _K and R _V are proportional. This immediately leads to

(23) $$ \lpar R_V h\rpar \vert _{{\bi x}\comma t} = {\gamma\over 1-\gamma} \lpar R_K h\rpar \vert _{{\bi x}\comma t} = \gamma \lpar Rh\rpar \vert _{{\bi x}\comma t}. $$

The system of the first four equations defined by (VI.30), (VI.37)–(VI.39), is closed by the assumptions (VI.34)–(VI.36) for S _b together with Eqn (23) and specified values for α, β _K and γ. Moreover, throughout [VI] (and [VII]), the authors assume in addition that the density profile is uniform. Then the depth-integrated change rate of the CED can be expressed straightforwardly in terms of ${\dot h}$ and the mass hold-up: the integral of the potential energy over z relative to the reference configuration is

(24) $$ R_V h = {\bar \rho} g_z h \lpar k-k_0\rpar = m g_z \left({h\over 2}-k_0 \right). $$

We will disregard terrain curvature and entrainment for the sake of simplicity. In the Lagrangian point-of-view, D _t m = 0. From this, we obtain directly

(25) $$ {\rm D}_t h = {{\rm 2}\gamma\over m g_z} \left[\alpha {\bi S}_{{\rm b}} \cdot{\bar {\bi u}} - \beta_K R_K h \right]. $$

This first-order differential equation describes relaxation of h to a (quasi-)steady-state value governed by the RKE. This enslavement of R _V and h is a direct consequence of the very strong assumption (VI.12). We emphasise that this first-order evolution equation for h is implied by the model assumptions stated in [VI] and that the additional three equations in (VI.30), (VI.37)–(VI.39) at best are superfluous. If one desires a model with more complicated dynamics in the z-direction, one has to replace the assumption (VI.12) by a weaker one that does not enslave R _V h.

Closer examination reveals that Eqn (25) does not strictly conserve energy. According to the assumptions made in [VI], part of the frictional work is directly converted to heat, the rest to Rh = (R _K  + R _V )h, where R _V is potential energy and R _K is fluctuation energy, i.e., 〈w′²〉 = 0 for a suitably defined time or ensemble average 〈.〉. However, if the density changes, D _t h ≠ 0 and $\overline {\rho w^2} \gt 0$ , thus there is a contribution from the slope-normal expansion to the total kinetic energy that is not accounted for in the modified model. As mentioned earlier, it might be more natural to include this non-random part of the kinetic energy in the CED; if one does so, energy is conserved. However, Eqn (24) is then no longer valid and a separate evolution equation for the flow depth must be constructed. As explained above, Eqns (VI.5)–(VI.7) are not suitable for this.

As in the case of the basic RKE model, an analysis of the timescales associated with different processes – in this case, relaxation of the depth-averaged velocity, the RKE, the velocity profile and the density to their steady-state values – is required for a consistent approximation of avalanche flow. We cannot pursue this question further here, but the tight coupling between h, $\hskip 0.9pt\bar {{\hskip -0.9pt{\bi u}}}$ and R _K revealed by Eqn (15), the rapid relaxation of R _K and our experience from studying alternative snow avalanche models all suggest that algebraic equations for R _K and h instead of differential ones would produce a simpler, yet more consistent and equally accurate model.

3.5. Mathematical formulation of the flow model

A final remark concerns [VI], Section 6, where all model equations are reformulated for implementation in a numerical code. The procedure for doing so is well known and correctly applied for the conserved quantities mass, momentum and RKE in Eqns (VI.30)–(VI.34). Section VI.5 then states that the three first-order evolution Eqns (VI.27)–(VI.29) for k, w and N _K (obtained from the erroneous third-order Eqn (VI.17)) are extended to include advection, which indeed is necessary. However, Eqns (VI.37)–(VI.39) present these equations in conservation form, which does not follow from the advected form for non-conserved quantities. For example, the difference between the conservative and advective extensions of (VI.27) is

$$ \partial_t k + {\bf \nabla}\cdot\lpar k\hskip 0.9pt{\bar {\hskip -0.9pt{\bi u}}}\rpar - {\rm D}_t k = {\bf \nabla}\cdot\hskip 0.9pt{\bar {\hskip -0.9pt{\bi u}}}. $$

Clearly, the (2-D) divergence ${\bf \nabla}\,\cdot \,\bar {{\bi u}}$ does not vanish identically.

4.1. General considerations

The notion of powder-snow avalanche is somewhat fuzzy in the literature. To avoid ambiguity, we will use the term ‘mixed snow avalanche’ (MSA) for flows that simultaneously feature three different flow regimes, namely dense flow (DF), light (or intermediate-density or fluidised) flow (LF) and suspension flow (SF) (see Sovilla and others (Reference Sovilla, McElwaine and Louge2015) for an in depth discussion of these regions). The original RKE-extension of the Voellmy model (Section 2) with constant density is, in principle, applicable only to the DF, but is in practice used to model both the DF and LF regimes. The variable-density model (Section 3) attempts to explicitly model transitions between the DF and LF regimes in a single-layer model. In [VII], the SF regime is added to the model through a second layer.

Based on today's knowledge from observations and measurements, the LF regime may be attained in small avalanches, but is typically more strongly developed in larger avalanches. The SF regime will not be attained unless a considerable part of the avalanche has reached the LF regime. Due to their elevated velocity, parts of the avalanche in the LF regime will reach farther than the parts in the DF regime, and the part in the SF regime may travel yet farther (by several kilometres in extreme cases).

There is probably a smooth transition between flow regimes. This would favour a mathematical description in terms of a multi-phase model (air and snow particles of different sizes), where the density and the stress tensor depend on the volumetric concentrations of the different particle size classes and where air turbulence plays an important role at low particle concentration. However, such a model would have to be formulated as a genuine 3-D model and would at present be poorly suited for practical applications. The different deposit characteristics (Issler and others, Reference Issler, Gauer, Schaer and Keller1996) of the DF and LF regimes as well as some measurements with Frequency-Modulated Continuous-Wave (FMCW) radar suggest that the three flow regimes nevertheless may often be fairly distinct, with large density gradients at their boundaries. This opens the way for models with several layers corresponding to different flow regimes and with depth averaging applied to each layer. Virtually all models proposed so far follow this path. We note in passing that recent detailed measurements in large MSAs at the Vallée de la Sionne test site in Switzerland (Sovilla and others, Reference Sovilla, McElwaine and Louge2015; Köhler and others, Reference Köhler, McElwaine, Sovilla, Ash and Brennan2016) suggest a more complex picture in which sudden, intermittent bursts of rather large and high-density volumes of snow particles play an important role – perhaps not unlike horse-shoe vortices detaching from the bottom surface in turbulent flows. If confirmed, these measurements may question the traditional approach of modelling avalanches in terms of continuum models with slowly varying, depth-averaged fields.

Considered in isolation, the SF is a turbulent suspension of small snow grains in air and a sub-type of particulate gravity currents (Simpson, Reference Simpson1987). The volume concentration of the grains is usually very low (≪0.1) so that grain interactions are not important, though due to the high density of the grains they carry most of the momentum. The excess density of the mixture over the air is referred to as the buoyancy of the current. The density is not constant as on the upper surface ambient air is entrained by turbulence and on the lower surface snow can be lost due to particle settling or gained due to entrainment. All models must therefore have at least three equations for momentum, air mass and snow mass or – equivalently and more commonly – buoyancy and volume. In its initial and final stages, the SF is in the Boussinesq regime where the average density is almost the same as the ambient air density, but it is far in the non-Boussinesq regime when fully developed.

There is a large body of experimental, theoretical and numerical work on density and particulate gravity currents in a variety of idealised settings. It is an important question which of these results remain valid in the case of MSAs and must be taken into account in the modelling of the SF regime. Density and turbidity currents in the laboratory are typically produced from a dilute initial suspension and run over a relatively smooth, non-erodible bed. In contrast, the SF in an MSA forms at the front and on top of the highly agitated LF, and mass is exchanged between the LF and the SF at a high rate. The particle concentration in the SF layer typically being less than 0.01 and the particle settling velocity <1 m s⁻¹, the mechanism of particle suspension inside the layer has to be the same as in other dilute particulate gravity currents and the processes at the upper surface that govern entrainment of ambient air have to follow the general laws observed in high-Reynolds number jet flows, plumes or thermals.

A number of theoretical analyses determine the front velocity of gravity currents by treating them as inviscid, energy-conserving flows without explicitly taking into account turbulence (Benjamin, Reference Benjamin1968; McElwaine, Reference McElwaine2005; Nokes and others, Reference Nokes, Davidson, Stepien, Veale and Oliver2008). The applicability of the results from this approach appears a priori questionable, but recent numerical studies nevertheless appear to confirm it at least for flows without a dense undercurrent (Konopliv and others, Reference Konopliv, Llewellyn Smith, McElwaine and Meiburg2016). In all these situations, one finds that the bed shear stress is negligible compared with the effect of ambient-fluid entrainment along the upper surface. It remains to be seen, however, whether this also holds true in MSAs, where there presumably is a pronounced vertical density gradient across the depth of the SF layer and where the surface of the DF/LF layer beneath can be strongly agitated. Both these flow properties increase the shear stress at the lower interface of the SF layer.

4.2. Basic modelling assumptions in [VII]

The approach proposed in [VII] is a two-layer formulation: a lower layer of intermediate to high (variable) density, consisting of large snow balls, fine snow grains and air for the DF and LF regimes, and an upper layer of low (variable) density containing air and fine snow grains in the SF regime. The air in the SF layer is treated as incompressible and the layer depth and mass per unit footprint area are used instead of air and snow mass. The relative motion between snow grains and air is neglected so that a single momentum balance equation is sufficient for the SF layer. Perhaps surprisingly, the balance of turbulent energy is not considered here, even though the concept of RKE is borrowed from the theory of turbulence, and turbulence is instrumental in maintaining the snow grains in suspension (Parker and others, Reference Parker, Fukushima and Pantin1986; Fukushima and Parker, Reference Fukushima and Parker1990).

We need not discuss the left-hand sides of the balance equations further because they have standard form. Paper [VII] postulates the following source terms for the conservation equations for mass, x− and y-momenta and volume of the suspension layer:

(VII.35) $$ {\bi G}_\Pi = \left(\!\matrix{\hskip 0.9pt\dot{\hskip -0.9pt M}_{\Phi \rightarrow \Pi} + \hskip 0.9pt\dot{\hskip -0.9pt M}_{\Lambda \rightarrow \Pi} \cr \hskip 0.9pt\dot{\hskip -0.9pt M}_{\Phi \rightarrow \Pi} u_\Phi - S_{\Pi x} \cr \hskip 0.9pt\dot{\hskip -0.9pt M}_{\Phi \rightarrow \Pi} v_\Phi - S_{\Pi y} \cr {\dot V}_{\Phi \rightarrow \Pi} + \hskip 0.9pt\dot{\hskip -0.9pt V}_{\Lambda \rightarrow \Pi} }\!\right). $$

$\dot{\hskip -0.9pt V}_{\Phi \rightarrow \Pi }$ and $\dot{\hskip -0.9pt V}_{\Lambda \rightarrow \Pi }$ are the volume fluxes from the dense core (Φ) and the ambient air (Λ) to the suspension layer (Π), and ${\hskip .5pt \dot{\hskip -0pt M}}_{\Phi \rightarrow \Pi }$ and $\dot{\hskip -0.9pt M}_{\Lambda \rightarrow \Pi }$ are the associated mass fluxes. S _Π denotes the surficial shear stress on the suspension layer. In what follows, we will change the notation from $\dot{\hskip -0.9pt M}_{\Phi \rightarrow \Pi }$ , $\hskip 0.9pt\dot{\hskip -0.9pt M}_{\Lambda \rightarrow \Pi }$ , $\hskip 3pt \dot{\hskip -0.9pt V}_{\Phi \rightarrow \Pi }$ , $\dot{\hskip -0.9pt V}_{\Lambda \rightarrow \Pi }$ to $Q_{\Phi \rightarrow \Pi }$ , $Q_{\Lambda \rightarrow \Pi }$ , $W_{\Phi \rightarrow \Pi }$ , $W_{\Lambda \rightarrow \Pi }$ because these quantities are not advective derivatives of M _Φ, M _Π, h _Φ and h _Π, as the dot notation would imply.

4.3. The role of gravity

A particular feature of the MSA model of [VII] – immediately apparent from the second and third components of G _Π in (VII.35) – is that gravity is neglected in the dynamics of the SF layer, both as the driving force and as the cause of sedimentation. This is different to every other model of gravity currents and contradicts a large body of well-documented research showing that gravity and air entrainment at the top surface are the dominant terms in the momentum balance of density currents (Hopfinger, Reference Hopfinger1983; Meiburg and others, Reference Meiburg, McElwaine, Kneller and Fernando2012). Such an approximation might be justified when describing the motion of jets of almost particle-free air ejected from the avalanche, but this would be devoid of practical relevance. Neglecting sedimentation is less grave unless one is interested in the late run-out phase of the PSA, where the flow rarely does damage.

The missing processes are easily introduced into the source term (VII.35):

(26) $$ {\bi G}_\Pi' = \left(\!\!\! \matrix{ Q_{\Phi\rightarrow\Pi} - Q_{\Pi\rightarrow\Phi} + Q_{\Lambda\rightarrow\Pi} \cr \lpar M_\Pi-h_\Pi\rho_a\rpar {\bi g} + \lpar Q_{\Phi\rightarrow\Pi}-Q_{\Pi\rightarrow\Phi}\rpar {\bi u}_i - {\bi S}_\Pi \cr W_{\Phi\rightarrow\Pi} - W_{\Pi\rightarrow\Phi} + W_{\Lambda\rightarrow\Pi} }\!\!\! \right)\comma $$

where ρ _a is the density of the ambient air. The buoyancy term (the left-most term in the middle row of Eqn (26)) is present in all earlier PSA models we are aware of. The grain-borne shear stress (middle term in the second row) represents the momentum flux from one layer to the other due to the mass flux. It is the product of the net entrainment rate (entrainment minus sedimentation) and the (slope-parallel) velocity at the interface, u _i . The latter needs to be modelled as a function of $\bar {{\bi u}}_\Phi $ , $\bar {{\bi u}}_\Pi $ and the densities in the two layers, but we will not go further into this question. The relation between mass and volume loss of the PSA cloud due to settling of snow grains is $Q_{\Pi \rightarrow \Phi } = \rho _i W_{\Pi \rightarrow \Phi }$ , assuming that all particles in the PSA cloud are snow grains with the density of ice, ρ _i. A candidate model for the sedimentation rate is the one used by Parker and others (Reference Parker, Fukushima and Pantin1986),

(27) $$ Q_{\Pi\rightarrow\Phi} \approx c_{{\rm b}} \left({M_\Pi\over h_\Pi} - \rho_{{\rm a}}\right)w_{{\rm s}} \cos\theta \comma $$

where M _Π/h _Π is the depth-averaged PSA density, w _s the average settling velocity of the snow particles, and θ the local slope angle. c _b, the ratio of bottom snow concentration to depth-averaged concentration, needs to be assumed, the most plausible values being in the range 3–10.

4.4. Air entrainment and drag on the suspension layer

Paper [VII] does not specify the air entrainment rate $Q_{\Lambda \rightarrow \Pi }$ appearing in (VII.35) and in the first row of Eqn (26). This entrainment rate has been measured repeatedly in inclined plumes and particulate gravity currents since the pioneering experiments by Ellison and Turner (Reference Ellison and Turner1959). It is well understood by now that it is governed by the Richardson number (Turner, Reference Turner1973). (Turnbull and others (Reference Turnbull, McElwaine and Ancey2007) provide a summary of these ideas applied to avalanches.) The Richardson number is the ratio of the potential energy to the kinetic energy of a parcel of fluid. For an entire layer, one defines the bulk Richardson number,

(28) $$ {\rm Ri} = {\lpar \rho-\rho_a\rpar g h \cos\theta\over \rho_a u^2}\comma $$

where θ is the angle between the slope normal and vertical. This implies that the entrainment rate depends both on the slope inclination, as shown experimentally, e.g., by Ellison and Turner (Reference Ellison and Turner1959) and Beghin and Olagne (Reference Beghin and Olagne1991), and on the density difference.

Based on laboratory experiments, Turner (Reference Turner1986) proposed the following formula for the entrainment coefficient:

(29) $$ E\lpar {\rm Ri}\rpar = \left\{\matrix{{} \displaystyle {0.08 - 0.1\, {\rm Ri}\over 1 + 5\, {\rm Ri}} \comma & \quad {\rm Ri} \lt 0.8 \comma \cr 0\comma & \quad {\rm Ri} \geq 0.8.} \right. $$

Ancey (Reference Ancey2004) more recently fitted unpublished data of Beghin by the function

(30) $$\hskip -38pt E({\rm Ri}) = \left\{\matrix{ {\hskip -4pt}{\rm e}^{-\lambda {\rm Ri}^2}\comma & \quad {\rm Ri} \leq 1\comma \cr {\rm e}^{-\lambda}/{\rm Ri}\comma & \quad {\rm Ri} \gt 1\comma } \right. $$

where λ = 1.6. This provides an even better closure for the volume and air-mass balances

(31) $$ W_{\Lambda\rightarrow\Pi} = E\lpar {\rm Ri}\rpar \vert \hskip 0.9pt{\bar{\hskip -0.9pt{\bi u}}}_\Pi\vert \quad{\rm and}\quad Q_{\Lambda\rightarrow\Pi} = \rho_a W_{\Lambda\rightarrow\Pi}. $$

Bartelt and co-workers assume the sum of the shear stresses on the upper and lower interfaces of the SF layer to be proportional to the square of the cloud velocity, u _Π ², and the cloud density, ρ _Π, writing

(VII.37) $$ {\bi S}_\Pi = -{\hskip 0.9pt{\bar {\hskip -0.9pt{\bi u}}}_\Pi\over \vert \hskip 0.9pt{\bar {\hskip -0.9pt{\bi u}}}_\Pi\vert } {g\over \xi_\Pi} \rho_\Pi\hskip 0.9pt {\bar{\hskip -0.9pt{\bi u}}}_\Pi^2.$$

The drag coefficient ξ _Π is considered a constant to be selected by the user. The authors state that the dominant contribution to the drag is from air entrainment at the upper boundary.

There are two issues with these assumptions: (i) If indeed air entrainment is the dominant contribution to the retarding forces on the SF layer, S _Π ≈ 0 would result because entrainment of ambient air at rest does not remove momentum from the SF layer, but distributes it over an increasing mass and thus decelerates the flow. To see this, multiply the mass-balance equation

$$ \partial_t\lpar h{\bar \rho}\rpar + {\bf \nabla}\cdot\lpar h{\bar \rho}\hskip 0.9pt{\bar {\hskip -0.9pt{\bi u}}}\rpar = Q $$

by $\bar {{\hskip -0.9pt{\bi u}}}$ and subtract it from the momentum-balance equation

$$ \partial_t\lpar h{\bar \rho}{\bar {\bi u}}\rpar + {\bf \nabla}\cdot\lpar h{\bar \rho}{\bar {\bi u}}{\bar {\bi u}}\rpar = \lpar {\bar \rho}-\rho_a\rpar h {\bi g} - {\bf \nabla} \lpar h {\bar p}\rpar . $$

This produces the equation of motion

$$ h{\bar \rho} {\rm D}_t {\bar {\bi u}} = \lpar {\bar \rho}-\rho_a\rpar h {\bi g} - {\bf \nabla}\lpar h{\bar p}\rpar - Q{\bar {\bi u}}\comma $$

which correctly features the decelerating force due to the acceleration of the ingested mass if this mass originally is at rest. (ii) Drag at the lower boundary should depend, not on the SF layer velocity, but on the difference between SF layer and DF/LF layer velocities. It is thus seen that the model presented in [VII] effectively assumes the interfacial shear stress between the DF/LF layer and the SF layer to be given by Eqn (VII.37), whereas the entrainment function $W_{\Lambda \rightarrow \Pi }$ is not specified in the paper. Replacing $\bar {{\bi u}}_\Pi ^2$ by $(\bar {{\bi u}}_\Pi \,-\, \bar {{\bi u}}_\Phi )^2$ in (VII.37) and specifying $W_{\Lambda \rightarrow \Pi }$ and $Q_{\Lambda \rightarrow \Pi }$ according to Eqn (31) with either (29) of (30) would then provide a physically consistent, but rather crude closure, provided the opposite term + S _Π is added to the momentum-balance equation of the DF/LF layer: the modified Eqn (VII.37) does not take into account that the interfacial shear stress will depend strongly on the densities in both layers.

The authors correctly note that air entrainment at the upper interface is the main source of resistance for the powder-snow cloud. However, there are also pressure forces acting, though the distinction between pressure drag and turbulent entrainment is not as straightforward as it first appears. For pressure to have a net retarding force there must be a separation region behind an object and a turbulent wake, otherwise we have D'Alembert's paradox that there is no drag. Thus, the momentum transfer from the pressure drag goes into a momentum deficit in the turbulent wake. This turbulent mixing region contains snow and if we regard it as part of the SF layer, then correctly accounting for the momentum balance here means that we do not need to directly consider the pressure drag. A more detailed discussion and quantitative analysis from a direct numerical simulation is contained in (Konopliv and others, Reference Konopliv, Llewellyn Smith, McElwaine and Meiburg2016).

This is correct if we consider the momentum balance over a large region containing the front, but if we wish for a more detailed model this can be improved. The pressure distribution along the SF layer surface is obtained by solving an elliptic potential flow problem in the ambient fluid. At the nose of the avalanche there will be a positive (with respect to the background pressure) stagnation pressure, but moving back along the top surface the pressure decreases linearly and may become negative (McElwaine, Reference McElwaine2005). Perhaps the best method would be to solve this potential flow problem using a boundary element method (Nokes and others, Reference Nokes, Davidson, Stepien, Veale and Oliver2008). This approach has not yet been directly applied to geophysical flows (see (De Blasio and others, Reference De Blasio2004) for a simplistic approximation). In addition, the shallowness approximation is violated since the front angle will be approximately 60 degrees (McElwaine, Reference McElwaine2005). Instead, the more usual approach is to apply a front condition by imposing a constant Froude number, thus setting the front velocity as a function of the front height.

In addition, added-mass effects (which give additional drag when an avalanche is accelerating) may also be relevant for powder-snow clouds of low density (Turnbull and others, Reference Turnbull, McElwaine and Ancey2007).

4.5. Formation of the suspension layer

The authors of [VII] interpret the cleft-and-lobe structure of MSA fronts as firm evidence for what they term blow-out. This proposed mechanism for SF layer formation can be summarised as a consequence of a hypothetical ‘breather’ mode in the motion of the dense core: the latter would expand periodically (normal to the slope) due to dispersive pressure and then contract again. During the expansion, air would be ingested and mixed with fine snow dust; when the large particles fall down again, the nascent powder-snow cloud would be expelled from the core as plumes (or rather puffs in conventional terminology), forming the observed lobes and clefts.

Paper [VII] specifies the boundary fluxes for the suspension layer as (in our notation)

(VII.24) $$W_{\Lambda\rightarrow\Phi} = 2 w_\Phi \Theta(w_\Phi),$$

(VII.25) $$ Q_{\Lambda\rightarrow\Phi} = \rho_{\Pi_0} W_{\Phi\rightarrow\Pi} \Theta(w_\Phi),$$

(VII.26) $$ W_{\Phi\rightarrow\Pi} = 2 w_\Phi \Theta(-w_\Phi),$$

(VII.27) $$ Q_{\Phi\rightarrow\Pi} = \rho_{\Pi_0} W_{\Phi\rightarrow\Pi} \Theta(-w_\Phi),$$

where the Heaviside distribution Θ(x) is 1 for x < 0 and 0 otherwise. Precisely speaking, (VII.24) and (VII.25) refer to the avalanche core, and as such (VII.24) is merely a kinematic statement. Together they describe the first half of one cycle in the ‘breather’ mode that is invoked here. The second half of the cycle transfers mass and momentum from the core to the SF layer and is described by (VII.26) and (VII.27).

As stated above and by (VII.24), the core entrains air (without ice dust) at the rate $W_{\Lambda \rightarrow \Phi }$ . In contrast, (VII.25) states that the core entrains mass from the ambient air at the rate $Q_{\Lambda \rightarrow \Phi } = \rho _{\Pi _0} W_{\Lambda \rightarrow \Phi }$ , i.e., air laden with ice dust. This is likely a typographical error, and ρ _Π0 should be replaced by ρ _a. Even so, $Q_{\Lambda \rightarrow \Phi }$ does not appear in (VII.34) and (VII.35), while the terms $W_{\Lambda \rightarrow \Pi }$ and $Q_{\Lambda \rightarrow \Pi }$ appearing in (VII.35) are not explicitly defined in [VII]. Thus, the first, second, third and fifth components in Eqn (VII.34) should be amended to read

(32) $$ {\bi G}_\Phi = \left(\!\! \matrix{ Q_{\Sigma\rightarrow\Phi} - Q_{\Phi\rightarrow\Pi} + Q_{\Lambda\rightarrow\Phi} \cr G_x - S_{\Phi x} - \lpar Q_{\Phi\rightarrow\Pi} - Q_{\Pi\rightarrow\Phi}\rpar u_{ix} + S_{\Pi x} \cr G_y - S_{\Phi y} - \lpar Q_{\Phi\rightarrow\Pi} - Q_{\Pi\rightarrow\Phi}\rpar u_{iy} + S_{\Pi y} \cr \cdots \cr - h_\Phi {\bf \nabla}\cdot{\bi u}_\Phi + 2 w_\Phi \cr \cdots }\!\! \right). $$

The first term of the fifth component corrects for the conservation form of the left-hand side of the balance Eqn (VII.28); see Section 3 for further necessary modifications of (VII.32) and (VII.34).

Besides these technical issues, the proposed SF formation mechanism presents two conceptual problems. First, lobes and clefts cannot be considered evidence for this mechanism. They form inevitably as instabilities in all types of highly turbulent gravity flows, with or without suspended particles, on completely flat and smooth surfaces (Simpson, Reference Simpson1987). Even if there is no dense underflow and the proposed ‘blow-out’ mechanism cannot be operative, the same type of structure is observed. Engulfing of ambient air by large eddies (similar in size to the flow depth) is the main entrainment mechanism in turbulent gravity currents and may indeed lead to the impression that dense jets are ejected from the flow. In MSAs, the density gradients are expected to be large except in the late stage; this will accentuate the impression of jets shooting out. To make their proposal more than speculation, the authors ought to show that the conventional mechanisms fail to reproduce central aspects of MSA behaviour or, ideally, perform measurements at the interface of the core and the suspension layer.

Another crucial aspect of the ‘blow-out’ mechanism is not addressed in [VI] or [VII], but in [V]: the core must undergo some kind of intermittency or oscillatory behaviour with a rather large amplitude to eject large volumes of snow–air mixture. Such behaviour indeed arose in the model described in [V] since the second-order equation for the vertical velocity derived in that paper has no damping term. However, as we show in Section 3, that equation lacks a proper physical basis.

Finally, we note that the ‘blow-out’ mechanism as described in [VII] requires that there is no suspension layer above the core during the first (expansion) half of the cycle; otherwise, the core could not entrain ‘pure’ ambient air, as is stated before (VII.24) and (VII.25). This means that the puff of snow–air mixture that is ejected during core contraction must immediately move away from the core segment that ejected it before the next expansion begins. This may work for the front, but core segments behind the nose will typically have a puff above them when they expand, and thus ‘swallow’ the suspension layer rather than ingesting ambient air. There is no indication in [VII] of how the model prevents the core from repeatedly incorporating what it has just expelled.