Comments on “On a Continuum Model for Avalanche Flow and Its Simplified Variants” by S. S. Grigorian and A. V. Ostroumov

Abstract

This note first summarizes the history of the manuscript “On a Continuum Model for Avalanche Flow and Its Simplified Variants” by Grigorian and Ostroumov—published in this Special Issue—since the early 1990s and explains the guiding principles in editing it for publication. The changes are then detailed and some explanatory notes given for the benefit of readers who are not familiar with the early Russian work on snow avalanche dynamics. Finally, the editor’s personal views as to why he still considers this paper of relevance for avalanche dynamics research today are presented in brief essays on key aspects of the paper, namely the role of simple and complex models in avalanche research and mitigation work, the status and possible applications of Grigorian’s stress-limited friction law, and non-monotonicity of the dynamics of the Grigorian–Ostroumov model in the friction coefficient. A comparison of the erosion model proposed by those authors with two other models suggests to enhance it with an additional equation for the balance of tangential momentum across the shock front. A preliminary analysis indicates that continuous scouring entrainment is possible only in a restricted parameter range and that there is a second erosion regime with delayed entrainment.

1. On the History of the Paper by Grigorian and Ostroumov

The paper “On a Continuum Model for Avalanche Flow and Its Simplified Variants” by Grigorian and Ostroumov [1] has a long, tortuous history, the knowledge of which is a prerequisite for assessing its value and for understanding why it is published in this Special Issue of Geosciences about a quarter century after it was written and almost half a century after important parts of the work had been completed.

Starting in the early 1960s, a vigorous research program on snow avalanche dynamics was started in the Soviet Union, with Samvel S. Grigorian (1930–2015) and Margarita E. Eglit at Lomonossov Moscow State University (MSU) as the key figures. They applied modern concepts of fluid mechanics and hydraulics to the problem, thus going far beyond Voellmy’s work [2]. While the latter was designed for simple engineering applications using a slide ruler, the effort at MSU included both powerful mathematical analysis and numerical simulation on digital computers, applying state-of-the-art numerical techniques. Already in 1966, a consistent quasi-two-dimensional (depth-integrated) hydraulic model with frontal snow-cover entrainment was formulated, implemented as a computer code and its results compared to avalanche experiments in the Caucasus [3]—more than two decades before the influential work of Savage and Hutter [4].

A fair part of the snow avalanche work by Eglit and her co-workers has been made accessible to western researchers through translations [5,6] and papers written in English [7,8,9,10,11], but only three papers by Grigorian in English translation are known to me [12,13,14], with [14] not explicitly focusing on snow avalanches. Snow avalanches represent but a small part of Grigorian’s research interests—they ranged from the mechanics and fracture of rocks and soil, surging glaciers, seismology, explosion dynamics and missile penetration over practical problems arising in deep-well drilling, bio-mechanical questions connected to blood circulation and its regulation, sports mechanics, to the theory of ball lightning, typhoons and tornadoes, tectonic volcanism, celestial mechanics and cosmology [15].

After the end of the Soviet Union, scientists in Russia and the other successor states were in dire straights, with many research institutions—notably defense-related ones—dismantled or salaries grotesquely lagging behind the galloping inflation. Recognizing the potential threats from this situation, the European Union financially supported INTAS (International Association for the Promotion of Cooperation with Scientists from the Independent States of the Former Soviet Union; see for example, [16]), which promotes joint research projects between research institutions in the FSU (former Soviet Union) and in Europe and provides funding for the partners from the FSU.

Having long-standing relations with leading Soviet snow and avalanche researchers, Dr. Bruno Salm from the Swiss Federal Institute for Snow and Avalanche research (SLF) obtained an INTAS grant for supporting Grigorian’s work on avalanche dynamics. Salm mentioned the project occasionally during lunch converations at SLF in the mid-1990s. The final report arrived by telefax in 1994 or 1995, but for unknown reasons was not circulated. At the occasion of SLF’s 60 year jubilee symposium in Davos in the fall of 1996, Grigorian was one of the invited speakers. When the overhead projector failed right before his talk and could not be repaired on the spot, Grigorian closed his eyes while facing the audience, put his finger tips against his temples, concentrated for two minutes, and then delivered an improvised, yet impressive talk without any slides or other graphical support. However, as the reader will immediately see by skimming the paper [1], it is not possible to convey any of the (important) mathematical details about the models or the results obtained with them in this way.

Before or shortly after the symposium, Grigorian with co-author Alexander V. Ostroumov, who had contributed to the original work in the 1970s and again in the 1990s, submitted a manuscript for the planned proceedings volume. For various reasons, putting together the proceedings volume at SLF took about three years. At that point, I was able to obtain a photocopy of Grigorian and Ostroumov’s paper. Even though finally completed, the proceedings volume has never been published, and it seems that the files with the electronically submitted papers are now lost.

In the context of the present Special Issue of Geosciences, I pondered the question whether Grigorian would be interested in finally publishing this manuscript, preferably in an updated form. From Margarita E. Eglit I learned, however, that Grigorian had passed away in 2015 and that Ostroumov had given her an old manuscript on this topic in electronic form, which turned out to be the very same paper except for the figures, which were not included in the electronic version.

To the best of my knowledge, Grigorian’s manuscript was not subjected to peer review or substantial editing at SLF between 1996 and 1999. Publishing it unedited two decades later would neither do justice to this work nor comply with the standards of a peer-reviewed journal. Since Grigorian no longer can approve or reject changes to the manuscript, I decided to (i) make the original manuscript accessible as Supplementary Materials Document (SMD) 1, (ii) to edit the text correcting some minor errors and emphasizing grammatical correctness, standard terminology and ease of reading over strict adherence to the original wording, (iii) to write figure captions and insert selected figures in the text, relegating the others to SMD 2 for the sake of succinctness, and (iv) in the present Comment paper to explain both my editorial choices and a number of technical points that may not immediately be understandable for readers unfamiliar with the early work of the MSU school. Ostroumov confirmed that he had no objections against publication and carefully checked the edited version of the paper as well as the Supplementary Materials. Four short essays (Section 3, Section 4, Section 5 and Section 6) discuss points of potential relevance for today’s research. I hope that this, together with Reference [17], will give western researchers a fuller appreciation of the pioneering work on avalanche dynamics from the Soviet era and further the progress in this field.

2. Notes on Editor’s Changes and on Technical Points in the Paper by Grigorian and Ostroumov

Note 0

The title was changed following the suggestion of an anonymous reviewer.

Note 1

The original manuscript has an abstract of considerable length, whereas the original Section 1 starts medias in res. To adhere to the guidelines for papers in Geosciences, the original abstract was made into the Introduction (Section 1) and a new, shorter abstract was written by the editor.

Note 2

Equation (3) can be derived from the modified dry-friction law stated below in Equation (4) by assuming that the normal pressure $p_n$ decreases linearly with the height. Then the two channel sidewalls combined contribute a force per unit length of $2\mu \rho a{h}^2/2$ in the case ${\tau }_1<{\tau }_{*}$. Dividing by the density $\rho$, the flow depth h and the channel width L to obtain an acceleration, one verifies the upper part of Equation (3). In order to obtain the formula for the case ${\tau }_1>{\tau }_{*}$, one first needs to find the height $z_{*}$ where the internal shear stress equals ${\tau }_{*}$ from the equation

$$\mu \rho a(h-z_{*})={\tau }_{*},$$

and then add the contributions from the channel bottom, ${\tau }_{*}L$, the channel walls up to $z_{*}$, that is, $2{\tau }_{*}{z}_{*}$, and the channel walls above $z_{*}$, which amount to $2{\tau }_{*}(h-z_{*})/2$. Again, the sum is to be divided by $\rho L$ to obtain the frictional deceleration.

Note 3

See Section 4 for the editor’s comments on Grigorian’s modified friction law. An obvious misprint in the lower line ($\mu p_n\le {\tau }_{*}$) was corrected to $\mu p_n\ge {\tau }_{*}$. The optical appearance of the conditions for the two alternatives was made consistent between Equations (3) and (4) without changing their mathematical content.

Note 4

The jump conditions (5) and (6) can be derived by considering the mass and momentum balances in infinitesimally thin boxes enclosing a short piece of the erosion front (see Figure 1). These boxes are assumed to move with the erosion front. The mass and momentum changes inside each of these boxes during a time interval $\Delta t$ must vanish—if this were not the case, the density and/or the speed of the boxes would diverge because the box volume and thus the mass and momentum content are infinitesimally small.

Let $\mathrm{d}l$ be the length of the box measured along the interface and $\omega$ and v the interface-normal velocities of the interface and the moving snow, respectively, right above the interface (measured relative to the fixed Earth coordinate system). Then the snow mass per unit width flowing from the snow cover into the box in a time interval $\mathrm{d}t$ is

$$\mathrm{d}{m}_{\mathrm{in}}={\rho }_0\omega \mathrm{d}l\mathrm{d}t.$$

On the avalanche side, the mass

$$\mathrm{d}{m}_{\mathrm{out}}={\rho }_1(\omega -v)\mathrm{d}l\mathrm{d}t$$

leaves the box. Equating these two mass increments and dividing by $\mathrm{d}l\mathrm{d}t$ immediately gives equality of the fluxes according to Equation (5).

In deriving the momentum jump condition, care is required with regard to the signs of momentum fluxes and the choice of coordinate system. Here we choose to describe the process in a coordinate system that moves with the interface and whose x-axis is perpendicular to the interface, pointing from the snowcover into the avalanche. Then the momentum flux (i.e., momentum per unit area and unit time) into the box is given by

$$J_{\mathrm{in}}={\rho }_0{\omega }^2+p_{*}.$$

$p_{*}$ is the compressive strength of the snowcover that is mobilized at the interface and is a non-negative quantity. It leads to momentum entering the box in the $+x$-direction and is therefore counted positive. The corresponding momentum flux on the avalanche side is

$$J_{\mathrm{out}}={\rho }_1{(\omega -v)}^2+p(h,u).$$

The pressure from the avalanche acts in the $-x$-direction, reducing the momentum inside the box. This corresponds to a positive momentum flux out of the box; the pressure $p(h,u)$ being a positive quantity, it is added to the mass-related momentum flux ${\rho }_1{(\omega -v)}^2$. Equating $J_{\mathrm{in}}$ to $J_{\mathrm{out}}$ and using ${\rho }_1(\omega -v)={\rho }_0\omega$ from the jump condition for the mass, one arrives at (6). Equation (7) then follows by simple algebra. For further comments on the assumptions behind this erosion model, see (Reference [9], [A.1.2.1]) and Section 5 of this paper.

Note 5

The correct formulation of the momentum balance equation in the presence of entrainment has been a source of confusion and heated debates for decades. The authors formulate the problem such that the eroded snow enters the avalanche with speed 0. Moreover, they write an equation of motion rather than a momentum balance equation and therefore—correctly—include the pseudo-force $-qu/F$, which arises because the eroded snow must be accelerated to the mean speed of the avalanche, u. See the beginning of Section 5.4 of this paper for details.

Note 6

The original computer code used in Reference [1] is lost, but A. V. Ostroumov [18] described its key features, which are summarized below. The approach anticipates many features of the numerical technique now known as the Material Point Method (MPM) in that two separate grids, one moving and one fixed, are used for numerically solving Equations (1)–(10) in Reference [1]. The equations for $F(S,t)$ (cross-sectional area of the flow), $u(S,t)$ (depth-averaged velocity), $w\left(t\right)$ (front velocity and $S_f\left(t\right)$ (front position) are solved on a moving grid G${}_1$ with a constant, user-selected number of nodes, $N_1$. The uppermost node $i=1$ is held fixed at the fracture line of the avalanche while the last node $i=N_1$ moves with the avalanche front. To avoid loss of precision due to large gradients of the flow depth $h\left(S\right)$ near the front of the avalanche, the cell length $\Delta S_1$ is kept constant at a user-selected value between the nodes $i=N_1-n_1$ and $i=N_1$ near the front. In contrast, the cell lengths between the rear nodes $i=1,\dots ,N_1-n_1$ grow uniformly in the process of avalanche propagation.

The variables $\delta (S,t)$ (erodible snow depth), $\alpha (S,t)$ (inclination of the snow surface relative to the terrain) and $q(S,t)$ (volumetric entrainment rate integrated over the avalanche width) are calculated on a spatially fixed grid G${}_2$ with uniform, user-selected spacing $\Delta S_2$. At each time step, the currently active portion of G${}_2$ is determined as follows: The rearmost active node of G${}_2$, $j_e$, is chosen such that $\delta (S_j,t)=0$ for $j<j_e$ and $\delta (S_j,t)>0$ for $j\ge j_e$, that is, behind node $j_e$ the snow cover is already completely eroded. The frontmost node, $j_f$, is set at $j_f=\lceil S_f\left(t\right)/\Delta S_2\rceil +5$; thus the number of active G${}_2$ nodes is $N_2\left(t\right)=j_f\left(t\right)-j_e\left(t\right)+1$.

Simulations are advanced by one time step from t to $t_1=t+\Delta t$ by first calculating the new values $\delta (S_j,t_1)$, $\alpha (S_j,t_1)$ and $q(S_j,t_1)$ on G${}_2$ by applying the formulas (5)–(9) to the known values $F(S_j,t)$ and $u(S_j,t)$, which are obtained by interpolating from G${}_1$ to G${}_2$. After interpolating the so obtained values from G${}_2$ to G${}_1$, $\delta (S_i,t_1)$, $\alpha (S_i,t_1)$ and $q(S_i,t_1)$ can be used to compute the new values $F(S_i,t_1)$, $u(S_i,t_1)$, $w(t-1)$, $S_f\left(t_1\right)$ on G${}_1$.

For calculating each of the different cases listed in (Reference [1], Table 1), suitable values of $N_1$, $n_1$, $\Delta S_1$ and $\Delta S_2$ are selected to obtain sufficiently precise results. The variable time step $\Delta t$ is determined from the stability condition of the numerical scheme.

Note 7

In the original manuscript, the authors erroneously state that they omitted the last two terms of Equation (1) in Reference [1], that is, $-f_1-f_2$. In reality, they drop the longitudinal pressure gradient $-(L/2F){\partial }_S(a{F}^2/L^2)$ and the entrainment deceleration $-f_2$, but retain the resistance term $-f_1$.

Similar assumptions have been used in applications to various types of gravity mass flows. If the momentum balance equations are integrated over the entire (instantaneous) length of the flow, the pressure-gradient term integrates to zero. Grigorian and Ostroumov do not explicitly carry out such an integration and therefore have to neglect the pressure-gradient term explicitly. If one assumes a constant shape (but variable length and height) of the longitudinal section of the flow, the local velocities and flow depths can be expressed in terms of the center-of-mass velocity (or, alternatively, the front velocity), the maximum instantaneous height over the entire body of the flow, and the assumed shape functions. With this, important local variables like the bed shear stress or the entrainment rate can be computed and integrated. In an early model of powder-snow avalanches, a half-ellipse was assumed as the shape function [19], whereas uniform height leads to so-called “box models” (e.g., Reference [20]).

Grigorian and Ostroumov integrate the mass balance equation over the avalanche length, but do not carry out this step for the momentum balance (or equation of motion). Instead, they argue that only the frontal region is of importance and, effectively, only write down the equation of motion of the front. The combined assumptions of no entrainment and constant longitudinal profile shape allow expressing the frontal flow depth simply as $H=H_0{l}_0/l$.

Note 8

In the original manuscript, the gravitational metric system (GMS) is used. This means that densities and stresses have units of kg s${}^2$ m${}^{-4}$ and kg m${}^{-2}$, respectively, that is, they are densities and stresses divided by the standard gravitational acceleration, $g\approx 9.80665\mathrm{m}{\mathrm{s}}^{-2}$. In the edited text, SI units have been used consistently together with the approximation $g\approx 10\mathrm{m}{\mathrm{s}}^{-2}$, but GMS units are used in some of the original figures.

Note 9

The Tables 3.2 and 4.1 in the original manuscript were combined into the single Table 1 to make the text more concise and to facilitate direct comparison between the results of the full model and different approximations.

Note 10

In the argument of the first exponential function of the second term under the root sign, a minus sign was omitted in the original manuscript. Note also that this formula assumes tacitly that the value of D is constant within each segment. If D changes at some point because $\tau$ crosses the value ${\tau }_{*}$, the segment needs to be subdivided at that point.

Note 11

The original manuscript reads $\gamma S=A/B$, which clearly is a misprint.

Note 12

The point-mass model can be derived from the full model by integrating over the length of the avalanche body. The gradient of the hydrostatic pressure thereby integrates to the difference of its values at the front and tail of the avalanche, which individually are 0. this reflects the fact that this is an internal force and does not influence the motion of the center of mass. Note, however, that in this derivation $sin\psi \left(S\right)$ and $cos\psi \left(s\right)$ should be understood as their mean values along the avalanche body at this instant.

Note 13

Subsection titles added by the editor.

Note 14

In the original manuscript, $\psi \left(\xi \right)$ is used instead of ${\psi }_1$. Since $\psi$ is declared constant, including an explicit $\xi$-dependence may be confusing. To emphasize that the value of $\psi$ to use is the one from the first (inclined) segment, the subscript ${}_1$ was added. Moreover, the condition $u_0\ge 0$ was added because the formula is not valid in the (admittedly unrealistic) case where the avalanche has an uphill initial velocity.

Note 15

The original manuscript reads u instead of $u^2$.

Note 16

The full run-out distance is measured along the path and not along its projection onto the horizontal plane, as would be more customary in hazard-mapping applications.

Note 17

The formula given in the original manuscript is

$$u=u_0\sqrt{1-\frac{gH(sin{\psi }_1-\mu cos{\psi }_1)}{k{u}_0^2}\left[exp\left(-\frac{2k}{H}(S-S_0)\right)-1\right]}.$$

This is correct only in the case $u_0=0$. If $u_0\ne 0$ and $S\to \infty$, it gives

$$u⟶\sqrt{u_0^2+{\displaystyle \frac{gH}{k}}(sin{\psi }_1-\mu cos{\psi }_1)},$$

instead of tending to the correct asymptotic value

$$u_{\infty }=\sqrt{{\displaystyle \frac{gH}{k}}(sin{\psi }_1-\mu cos{\psi }_1)}.$$

Solving the differential equation, one finds that a factor $exp[-(2k/h)(S-S_0)]$ was omitted from the first term under the root. Moreover, if $u_0<0$, the original equation gives $u\left(S\right)<0$ all the while the avalanche is moving in the $+S$-direction.

After Equation (5.16) in the original manuscript (corresponding to (57) in the edited text), the authors stated that the run-out distance does not depend on the size of the avalanche, even though both the original and corrected equation show an approximately linear dependence on the flow depth H. This statement has therefore been deleted.

Note 18

A small misprint ($L_2$ instead of L on the left-hand side) has been corrected.

Note 19

The original manuscript erroneously subtracts $l_1$ from the right-hand side of Equation (64). That would give the same result as for $l_2$ in (63).

Note 20

A spurious factor of H in the denominator of the expression for $u_{**}$ has been removed.

Note 21

The original manuscript states $l_2/l_1\approx 0.125/8=0.0156$ in the first case with Coulomb dry friction and $l_2/l_1\approx 0.25/8=0.0313$ for the case of stress-limited dry friction.

Note 22

Compared to friction parameters typically used in the simulation of medium-size to large avalanches with Voellmy-type models (see, e.g., Reference [21]), $k=0.02$ (corresponding to $\xi =500{\mathrm{m}\mathrm{s}}^{-2}$) appears rather large. With a more typical value $k=0.005$ (or $\xi \approx 2000{\mathrm{m}\mathrm{s}}^{-2}$), $l_2/l_1\approx ...$ results with Coulomb-type dry friction and $l_2/l_1\approx ...$ for the stress-limited friction law.

Note 23

A. V. Ostroumov confirmed that the authors observed this behavior only in simulations with both entrainment and the stress-limited friction law [18].

Note 24

This table was compiled by the editor on the basis of the information contained in Figures 8 and 9 and S2.10–S2.14.

3. Grigorian and Ostroumov’s Work in the Context of Present-Day Avalanche Dynamics Research—Is There Still a Use for Simple Avalanche Models in the 21st Century?

The rapid increase of computational power at the users’ fingertips has made it possible to numerically simulate complex natural processes in an ever increasing degree of detail. During the past decade, codes simulating depth-averaged snow avalanche flow on arbitrary surfaces have become the de facto standard in hazard mapping. With the widespread availability of digital terrain models with a grid resolution of 10 m or better and user-friendly geographical information systems, quasi-3D avalanche simulations can now be set up and carried out more rapidly than calculations with simpler models along a profile line, which needs to be carefully defined by an expert. When it comes to avalanche dynamics, the question has been raised, however, whether our knowledge of the underlying physical processes, the material properties and the initial conditions is sufficient to justify the use of complicated models [22]. Moreover, complicated non-linear models may exhibit unexpected solutions or chaotic behavior, that is, infinitesimal differences in the initial conditions may lead to qualitatively different solutions.

This ongoing debate has not always been very fruitful because one approach is usually put against the other. In contrast, the results of this paper suggest that combining simple and more comprehensive models may often increase insight. The simple models often allow one to understand the effect of one process or one specific parameter by inspecting the analytic solution of the problem. They also make it easier and quicker to scan a multi-dimensional parameter space. Grigorian and Ostroumov show, however, in Section 5 that large deviations from the solution provided by the full model can arise if simplification is carried only a little too far. In the present case, the evolution of the maximum flow depth near the front was found to be crucial because it strongly influences the retarding force per unit mass. Thus the simplest (mass-point) model is useful for comparing the general behavior of the solutions for different friction laws, but hardly for obtaining approximate solutions to the full equations. In contrast, the simplified (box) model captures the effect of decreasing maximum flow depth quite well, but at the expense of two dynamical variables and a substantially more involved analytic treatment.

The analysis presented in Reference [1] applies to models of dense-flow avalanches with the Voellmy friction law and variants thereof. Much of the approach to simplification should be applicable to other rheological assumptions or to models of powder-snow or mixed avalanches.

My personal experience with students at the master’s level in the geosciences suggests that the present paper might also serve as teaching material in self-study or courses on natural hazards. These students often have a good background in geomorphology, geology, geotechnics or geomatics, but most of them lack a solid understanding of the physical and mathematical foundations underlying the modeling of the dynamics of gravity mass flows. Modern software packages make it easy to run simulations, but they provide little assistance in the critical interpretation of the simulation results. A student who has worked through Section 2 of Reference [1] will easily understand the assumptions and simplifications upon which modern two-dimensional depth-averaged models are built. After assimilating the rest of the paper, their mathematical skills will have increased significantly (even though hardly anything beyond the level of first-year or at most second-year mathematics in a physics or engineering curriculum is required). Additionally, they will have gained insight into the key features of the solutions to the simplified equations in an archetypical situation in hazard mapping, namely an avalanche path with a hockeystick profile. Such a foundation in the mechanics of these phenomena will contribute greatly to the quality of their hazard-mapping and hazard-mitigation work in their future professional activities.

4. A Brief Discussion of Grigorian’s Stress-Limited Friction Law

Like the majority of snow avalanche models developed and used in the West over the past six decades, the Soviet work on avalanche dynamics—including large parts of the present paper—mostly uses Voellmy’s [2] friction law, which can be expressed as

$${\tau }_b=\mu {\sigma }_n+\frac{g}{\xi }{u}^2\equiv \mu {\sigma }_n+k{u}^2.$$

(1)

The bed shear stress, ${\tau }_b$, is written as the sum of a velocity-independent Coulomb-type contribution proportional to the normal stress on the bed, ${\sigma }_n$, and a Chézy-type term proportional to $u^2$, the square of the velocity. The dimensionless coefficients $\mu$ and k can be considered as the tangent of the bed friction angle and a drag coefficient, respectively. The Voellmy friction law has several properties that are widely believed to at least qualitatively agree with observations [2]: (i) Avalanches may stop on sufficiently gentle slopes, so there must be some velocity-independent component of friction. (ii) Avalanches may have a tendency to reach a terminal velocity on a long inclined plane, which requires the bed shear stress to increase with velocity (see, however, References [23,24] for recent inferences based on statistical analysis of aggregated measurements). (iii) Larger avalanches (with larger flow depth) attain higher velocities and longer run-out [24]. Besides this, Voellmy’s friction law has important shortcomings and can only be considered a heuristic model of a granular flow.

One should expect the Voellmy friction law to apply equally well to rock avalanches and landslides. These phenomena exhibit an intriguing dependence on the size: The run-out ratio R is defined as the total drop height divided by the horizontal run-out distance from the escarpment or fracture line to the distal end of the deposit and is essentially equal to the effective friction coefficient ${\mu }_{\mathrm{eff}}$. Observations of giant landslides show that $R\propto V^{-a}$ with $a\sim \frac{1}{6}$ and that R can attain values that are smaller than any plausible value of ${\mu }_{\mathrm{eff}}$ or $\mu$. With the Voellmy friction law, however, the inequality

$$R\approx {\mu }_{\mathrm{eff}}>\mu$$

(2)

must necessarily hold so that it cannot reproduce the scaling behavior of large landslides, unless one chooses $\mu$ according to the size of the landslide.

In a very short paper [14], Grigorian proposed a simple solution to this conundrum. He asserted that the Coulomb friction law (or more precisely, the Coulomb failure criterion), which stipulates proportionality of shear stresses and normal loads, cannot be valid for normal stresses beyond some material-dependent limit. He therefore postulated the existence of an upper limit ${\tau }_{*}$ of the shear strength that a given material can mobilize and modified the Coulomb criterion accordingly (Equation (3) of Reference [1]; see Section 5.5 for a discussion of the relation of ${\tau }_{*}$ to the snow strength $p_{*}$ in the erosion model):

$${\tau }_b=min\left(\mu {\sigma }_n,{\tau }_{*}\right)+k{u}^2.$$

(3)

No limitation is imposed on the Chézy-type contribution to the shear stress. With Grigorian’s friction law (GFL), the run-out ratio will decrease as $h^{-1}$ for large flow depths h, or as $V^{-1/3}$ if one supposes that the length and width scale linearly with h. Grigorian considered this stress-limited friction law a break-through with many diverse engineering applications [15] and also applied it to snow avalanches [1]. It has, however, received very little attention in the West—the author is not aware of any publication of non-Soviet researchers citing it.

Grigorian does not discuss the physical mechanisms that would limit the shear stress in a geophysical medium irrespective of the applied normal stress, but some insight can be gained by discussing possible scenarios. Let us first consider highly contractive and completely saturated soils or quick clays. They are examples of materials whose intricate texture will be destroyed if the shear stress exceeds some threshold. Large amounts of pore water are liberated thereby and lead to complete liquefaction. The residual yield strength can be orders of magnitude smaller than the peak strength during initial loading. Since the effective stress vanishes in this case and the whole load is borne by the (essentially incompressible) pore water, the shear stress will depend on the shear rate, but there is no obvious reason for the yield strength to grow with the normal stress. The yield strength of sensitive clay is almost independent of the normal stress applied in a controlled experiment; it depends, however, about linearly on the overburden it was exposed to in its natural environment. Not having stress-history dependence built in, the GFL law is too simple to fully capture such materials, but they support his notion of an upper limit of the frictional shear stress, independent of the normal stress.

Another important scenario, especially in connection with mega-landslides, concerns a granular material at the interface under such high load that the granules are crushed. Crushing will lead to denser packing and a larger average number of grain–grain contacts because small debris can fill the voids between larger particles. In this way, the load will be distributed over a larger number of force chains and the compressive or shear stress acting on a single particle will eventually decrease below the strength of the material. At that stage, one expects, however, the Coulomb friction law to become applicable again, perhaps with a somewhat lower friction angle because the particles are more rounded. The comminution of the granular material is expected to have a significant effect on the collisional shear stresses: According to the kinetic theory of granular flows, the collisional stresses are proportional to the square of the particle size. This would mean that comminution will allow large flows to attain much higher speeds than small ones, but it is questionable whether this mechanism can explain the observed volume dependence of the run-out angle of mega-landslides.

One may further consider flows in which the normal and shear stresses are large enough to melt a thin shear layer at the interface between the flow and the substrate. Prime candidates for such a mechanism are glaciers, wet-snow avalanches or very large landslides. Such melting will not occur along the entire bed–flow interface at once, but patch-wise; the size of the melted patches will grow with the normal and shear stresses. Accordingly, one should expect the average value of the dry-friction coefficient $\mu$ to decrease with the normal stress. The GFL may well capture the qualitative features of the bed shear stress evolution in such situations.

The GFL appears to be less suited for describing the surprisingly low friction in the (fluidized) head of dry-snow avalanches that is primarily observed under conditions with a substantial layer of cold, light new snow. In such avalanches, the largest normal basal stresses tend to occur in the non-fluidized, dense body of the flow. While Reference [1] does not consider density changes in avalanche flow, full-scale experiments at Mt. Elbrus (see Reference [25], Figure 60 on page 195) for a succinct summary of the experimental findings), in which Grigorian was involved, indicated the existence of an intermediate, fluidized flow regime. Grigorian and Ostroumov may well have been aware of the challenges offered by the experimental results when they wrote the last paragraph of Reference [1].

These considerations suggest that the GFL should be considered a heuristic description of diverse and complex mechanisms that limit the growth of the shear stress in some systems when they are subjected to normal loads that are large compared to some system-specific, intrinsic scale. As the examples discussed above suggest, detailed physical models of these mechanisms may be complicated, difficult to validate in detail and unsuitable for practical applications. In such situations, the GFL offers a simple and elegant alternative, especially if one can estimate the order-of-magnitude of ${\tau }_{*}$ from a simplified formulation of the relevant process. For example, an estimate of ${\tau }_{*}$ for rock avalanche events exhibiting basal melting and frictionite formation would allow comparing the predictions of Grigorian’s model to those of a model that incorporates the melting process explicitly [26].

To conclude this discussion of the GFL, we ask whether it could be used to absorb the strong volume dependence in the standard calibration of the friction coefficients $\mu$ and k of Voellmy-type models for dry-snow avalanches. Figure 2 shows the volume dependence of $\mu$ according to Reference [21] for avalanches with a return period of 300 years. The different curves correspond to selected combinations of three altitude zones and four classes of terrain characteristics. Experience confirms that these values give plausible run-out distances in the majority of cases. In the range $5000{\mathrm{m}}^3<V<60,000{\mathrm{m}}^3$

$$\mu (V,\dots )\propto V^{-0.25}$$

(4)

captures the volume dependence of the calibrated values reasonably well. The GFL implies instead

$$\mu \propto h^{-1}\mathrm{if}{\mu }_0\rho hgcos\theta >{\tau }_{*}.$$

(5)

${\mu }_0$ is the dry-friction coefficient for $\left|\tau \right|<{\tau }_{*}$. To compare the calibration of the Voellmy model to the GFL, one has to assume a (statistical) relation between the avalanche volume and the fracture depth/flow depth. Postulating $h\propto V^{0.25}$, GFL matches the behavior of Equation (4) in the range $5000{\mathrm{m}}^3<V<60,000{\mathrm{m}}^3$ (Figure 2). However, to avoid unrealistically low values for very large avalanches, h must be assumed to scale with an exponent much closer to 0 (Figure 2, yellow curve).

In which range is ${\tau }_{*}$ expected to lie? To match the calibration [21] of $\mu$, ${\tau }_{*}$ should typically be reached in avalanches with a volume of 5000–10,000 m${}^3$, corresponding to fracture or flow depths $h\sim 0.5$ m. With the density in the range 200–300 $\mathrm{kg}{\mathrm{m}}^{-3}$, one estimates ${\tau }_{*}=1$–1.5 kPa, which is encouragingly similar to the shear strength of a new-snow layer. To test the practical applicability of the GFL, one needs to compare its overall performance in a large number of avalanche cases to that of a Voellmy-type model with the standard calibration. The test cases should span a wide range of avalanche sizes, climatic conditions and terrain characteristics.

5. Some Remarks on the Grigorian–Ostroumov Erosion Model

From today’s perspective, the most relevant aspect of Reference [1] is how erosion and entrainment are modeled. The Grigorian–Ostroumov formula for the erosion rate was applied in References [27,28], but is not implemented in any code used in practice today. Below we discuss some aspects of this remarkable formula and try to assess its applicability in dynamical models of snow avalanches.

5.1. How Does the Grigorian–Ostroumov Erosion Formula Relate to Other Erosion and Entrainment Models?

Already the first depth-averaged avalanche model by the Moscow State University group [3]—called the Eglit–Grigorian–Yakimov erosion model (EGYEM) below—described entrainment, in this case plowing entrainment at the front, formulated as a jump condition at the moving boundary of the computational domain. This entrainment mechanism can often be observed visually in slowly moving wet-snow avalanches but is absent in fast-moving dry-snow avalanches, which tend to have relatively low density at the front. Such avalanches erode the snow cover more gradually along the bottom of the flow. The Grigorian–Ostroumov erosion model (GOEM) was developed to address this situation [1]. A simpler, yet conceptually related model for basal entrainment in depth-averaged flow models that assume a uniform velocity profile and sliding at the bottom was formulated in References [29,30] and is termed tangential-jump entrainment model (TJEM) in the following. A brief comparison of the three models will help understanding the particular features of the GOEM.

These three models conceptualize the avalanche–snow cover interface as a moving non-material shock front where the eroded snow abruptly changes its velocity and usually also its density. Accordingly, jump conditions for mass and momentum (and total energy) apply at the shock. The compressive or shear strength of the snow cover enters the jump conditions for the momentum and together with the dynamical variables—flow depth, velocity, density—determines the entrainment rate uniquely; it can, in principle, be measured independent of an avalanche event, thus there are no empirical parameters that need to be fixed a posteriori. This is a fundamental difference from the vast majority of entrainment models and crucial for practical applications because no further parameter is introduced that needs to be fitted from scarce or high-uncertainty data.

The concept of the bed–flow interface being a shock applies to situations where the velocity and possibly the density change considerably in the process. One expects this to be the case if the bed material can be characterized as brittle so that erosion is a sudden process rather than a relatively slow deformation over an extended period. If the eroded material is strongly dilatant or contractive, a density shock is to be expected. These conditions are met in dense, strongly eroding gravity mass movements like snow avalanches, debris flows and pyroclastic flows that flow over beds that are granular or fracture into a granular material. They are also fulfilled in dilute flows if the erosion rate is strong, for example, in powder-snow avalanches, turbidity currents and nuées ardentes. The shock concept applies neither to slow flows over a highly plastic, cohesive bed, in which it is hard to define an interface between substrate and flow, nor to flows that only sporadically manage to erode a particle from a very strong bed, for example, a torrent flowing over bedrock.

Figure 3 shows the control volumes across which the jump conditions are evaluated. The EGYEM evaluates the jump conditions of longitudinal momentum across a vertically oriented control volume at the front end of the (expanding) computational domain. In contrast, the TJEM considers the jump in tangential momentum (determined by the longitudinal velocity and the shear stresses); it neglects the variable inclination of the erosion interface relative to the terrain surface. The GOEM accounts for the inclination of the interface relative to the terrain and evaluates the jump conditions for the momentum components normal to the interface.

5.2. The Eglit–Grigorian–Yakimov Model for Frontal Entrainment

The jump conditions in the EGYEM are essentially the same as in a hydraulic jump or bore in hydraulics. The only difference is that the pressure from the downstream side is not hydrostatic fluid pressure, but the resistance of the (solid) snow cover. A derivation of the shock velocity $\omega$ in terms of the depth $\delta$, density ${\rho }_0$, and compressive strength ${\sigma }_{*}$ of the snow cover as well as the frontal flow depth $h_f$, density ${\rho }_1$, and near-front velocity $u_f$ can be found in, for example, Reference [17], the result being

$$\omega =\frac{h_f^2{\rho }_{1}gcos\psi }{2{\rho }_0{u}_f\delta }-\frac{{\sigma }_{*}}{{\rho }_0{u}_f}.$$

(6)

The EGYEM has a number of notable features: (i) The differential mass balance equation and the equation of motion do not contain entrainment terms—all entrainment occurs at the boundary of the computational domain and manifests itself as a boundary condition. (ii) Rearranging Equation (6) as $(h_f^2/2){\rho }_{1}gcos\psi ={\sigma }_{*}\delta +{\rho }_0\omega u_f\delta$ shows directly how the hydrostatic pressure exerted by the avalanche head is used partly for overcoming the resistance ${\sigma }_{*}\delta$ of the entrainable snow cover, partly for accelerating the eroded snow cover mass ${\rho }_0\omega \delta$ to the internal velocity $u_f$ of the avalanche. (iii) The density change from ${\rho }_0$ in the snow cover to ${\rho }_1$ in the avalanche is assumed given, but it could, in principle, be calculated from an appropriate constitutive law. (iv) The model idealizes the front as an infinitesimally thin shock, but in reality it is extended in the flow direction and the snow motion inside it is complicated, with comminution of the snow cover, compression or dilation from ${\rho }_0$ to ${\rho }_1$ and lifting of the center-of-mass from the height $\delta /2$ to $h_f/2$. As in hydraulic jumps, some of the pressure work done by the avalanche is dissipated across the shock. The shape of the front is largely determined by the effective friction angle, implying that the control volume in reality extends over a length of the order of the flow depth.

Equation (6) assumes the hydrostatic pressure in the avalanche front to be sufficient to erode and entrain the entire erodible snow cover. In (Reference [17], Equation (5)) three alternative modifications of the model in the case of Equation (6) predicting $\omega <0$ are suggested, but no unique solution is singled out. However, the mass jump equation ${\rho }_0\omega \delta ={\rho }_1(\omega -u_f)h_f$ implies the strong constraint $\omega /u_f=1/\left(1-\frac{{\rho }_0\delta }{{\rho }_1{h}_f}\right)$, that is, $\omega >u_f$; this problem deserves further study.

5.3. Main Features of the Tangential-Jump Entrainment Model

The TJEM neglects the contribution of normal pressure to the erosion process; instead, the snow cover is assumed to fracture at a critical value of the shear stress, $\left|\mathit{\tau }\right|\ge {\tau }_c$. The difference between the shear stress exerted by the avalanche and ${\tau }_c$ is available for accelerating the eroded mass from rest to the avalanche velocity $\mathit{u}$. Thus the tangential momentum jump condition leads to the entrainment rate

$$q=-{\rho }_0\omega =\Theta \left(\right|\mathit{\tau }|-{\tau }_c)\frac{\left|\mathit{\tau }\right|-{\tau }_c}{\left|\mathit{u}\right|}.$$

(7)

Here, $\Theta \left(x\right)=1$ if $x>0$ and 0 otherwise. When erosion occurs, $\omega <0$.

Although disregarded in Reference [29], the jump conditions for mass and bed-normal momentum may also be evaluated. The mass jump condition, ${\rho }_0\omega ={\rho }_1(\omega -w)$, implies that the flow body must have a small negative (positive) velocity w normal to the terrain if the substrate is contractive (dilatant). The jump condition for bed-normal momentum states that the eroded mass is accelerated from 0 to w by the jump in normal pressure across the interface. This effect has little practical significance, however.

At first sight, the TJEM appears to achieve closure without any adjustable parameters and is the unique consistent solution under the stated conditions. However, there is a critical hidden assumption: The bed shear stress ${\mathit{\tau }}_b$ over an eroding bed as a function of flow velocity $\mathit{u}$ and depth h differs in general from that over a non-erodible bed under the same flow conditions. If one considers that the interface is, in reality, a thin layer across which the eroded particles are accelerated, it is evident that the shear rate or sliding velocity at the interface between this accelerating layer and the quasi-stationary flow is less than it is at the bottom of a non-eroding flow. A simple friction law like Voellmy’s, which depends only on the depth-averaged velocity $\overline{\mathit{u}}$ and h does not contain enough information for finding a consistent solution [30]. If one uses the standard friction law for non-erodible beds in Equation (7), one will likely overestimate the erosion rate.

5.4. Particular Features of the Grigorian–Ostroumov Erosion Model

The GOEM is based on a control volume parallel to the instantaneous interface between the avalanche and the snow cover, which is inclined by the angle $\alpha (S,t)$ relative to the terrain. Like the EGYEM, the GOEM evaluates the jump conditions for mass and the interface-normal momentum component across the control volume, which in this case read

$$\begin{array}{cc}\hfill {\rho }_0\omega & ={\rho }_1(\omega -v),\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\rho }_0·0·\omega -p_{*}& ={\rho }_{1}v(\omega -v)-p_{\alpha }.\hfill \end{array}$$

(9)

$\omega$ and v are the interface-normal components of the shock propagation velocity and flow velocity, respectively. $-p_{\alpha }<0$ is the interface-normal stress exerted by the avalanche, $p_{*}$ the (compressive) strength of the snow cover. As in the EGYEM and TJEM, the avalanche density, ${\rho }_1$, is assumed known rather than computed from a constitutive equation. Equation (9) is valid only if $p_{\alpha }>p_{*}$; if not, $\omega =0$. This leads to

$$\omega =\Theta (p_{\alpha }-p_{*})\sqrt{\frac{p_{\alpha }-p_{*}}{{\rho }_0(1-{\rho }_0/{\rho }_1)}}.$$

(10)

Due to the factor $1-{\rho }_0/{\rho }_1$, $\omega$ would become imaginary if ${\rho }_1<{\rho }_0$. In contrast to the EGYEM and TJEM, the GOEM is thus not applicable to the fluidized and suspended parts of mixed avalanches.

Surprisingly, the Formula (10) is independent of the avalanche velocity u, and the erosion speed $\omega$ grows with the square root of the excess pressure, $p_{\alpha }-p_{*}$, whereas both the EGYEM and TJEM lead to an inverse dependence on the flow velocity and a linear dependence on the relevant excess stress. Technically, this comes about because v is eliminated in (9) with the help of (8). The same procedure could be applied to the EGYEM, leading to

$$\omega =\Theta (h_f{\overline{p}}_f-{\sigma }_{*}\delta )\frac{h_f{\overline{p}}_f-{\sigma }_{*}\delta }{{\rho }_0\left(1-\frac{{\rho }_0}{{\rho }_1}\frac{\delta }{h_f}\right)\delta }$$

(11)

instead of Equation (6), with ${\overline{p}}_f\equiv (h_f/2)gcos\psi$. Conversely, in deriving the GOEM one can choose to substitute ${\rho }_0\omega$ for ${\rho }_1(\omega -v)$ in Equation (9) to obtain

$$\omega =\Theta (p_{\alpha }-p_{*})\frac{p_{\alpha }-p_{*}}{{\rho }_{0}v},$$

(12)

showing that the erosion or entrainment speed has the same dependence on the excess stress and the interface-normal component of the flow speed in all three models.

In their review of erosion and entrainment models, Eglit and Demidov ([9], A.1.2.1) criticized the expression that Grigorian and Ostroumov use for the pressure at the interface between the avalanche and the snow cover,

$$p=\rho ahcos\alpha +C\rho {\overline{u}}^2.$$

(13)

$\overline{u}$ is the depth-averaged velocity, h the flow depth, and $\rho$ the mean avalanche density, which may differ from the density at the interface, ${\rho }_1$; $a\equiv gcos\psi +\varkappa {\overline{u}}^2$, where $\varkappa$ is the curvature of the path profile. There is no explanation given in Reference [1] why the second term, which resembles dynamic pressure or turbulent friction drag, should be present; neither is the origin and physical significance of the parameter C discussed. Eglit and Demidov argue that only the so-called static pressure should enter the momentum balance equation from which the jump condition is derived; the stagnation or dynamic pressure is already accounted for by the advective terms. This may explain why the entrainment speeds simulated in Reference [1] are an order of magnitude larger than observed in Nature—whenever erosion occurs in the simulations, a snowpack of 1 m is typically eroded over a distance of 5 m or less.

The simplest way to amend the GOEM would be to drop the velocity-dependent term in (10), simultaneously disposing of the empirical parameter C. This leads to smaller, more realistic entrainment rates with significantly lower values of $p_{*}=\mathcal{O}\left(1\mathrm{kPa}\right)$ rather than the values $\mathcal{O}\left(5\mathrm{kPa}\right)$ used in (Reference [1], Table 3). However, a more thorough analysis is warranted because the GOEM does not enforce the jump condition for the momentum tangential to the control volume, effectively treating an oblique shock as a normal shock. This work remains to be done, but a few pointers can be given here. First consider that, even though the GOEM only involves the normal stress on the interface between the avalanche and the snow cover, there is also a shear stress, ${\mathit{\tau }}_{\alpha }$. Accordingly, the maximum shear stress (inclined relative to the interface) in the new-snow layer is

$${\tau }_{\max }=\sqrt{\frac{1}{4}{\left(\frac{1}{1-\nu }\right)}^2{p}_{\alpha }^2+{\parallel {\mathit{\tau }}_{\alpha }\parallel }^2}\equiv \sqrt{b^2{p}_{\alpha }^2+{\tau }_{\alpha }^2},$$

(14)

where $\nu \approx 0.2$–0.35 is Poisson’s ratio for snow [31] and $b\equiv 1/(2-2\nu )$ is in the range $1/2\le b\le 1$, typically $b\approx 0.6$–0.8. With the typical range ${\tau }_{\alpha }\equiv \parallel {\mathit{\tau }}_{\alpha }\parallel \sim (0.2-0.6)p_{\alpha }$, the two terms in (14) are of similar magnitude, but there may be a significant velocity dependence in ${\tau }_{\alpha }$. The threshold condition for erosion in (10) thus should be $\Theta ({\tau }_{\max }-{\tau }_c)$, with ${\tau }_{\max }$ from Equation (14) and ${\tau }_c$ the shear strength of the snow cover. Next, examine $p_{\alpha }$ and ${\tau }_{\alpha }$; they refer to the interface, which is inclined at the angle $\alpha \equiv arctan\left({\partial }_S\delta (S,t)\right)$ to the terrain, relative to which the flow model is formulated. The tensor transformation laws lead to

$$\begin{array}{cc}\hfill p_{\alpha }& =p_b+{\tau }_{b}sin\left(2\alpha \right),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {\tau }_{\alpha }& ={\tau }_{b}cos\left(2\alpha \right).\hfill \end{array}$$

(16)

The hydrostatic pressure $p_b\equiv {\rho }_{1}ah$ is assumed isotropic. If ${\tau }_b$, the bed shear stress tangential to the terrain, is calculated from the Voellmy friction law, Equations (13) and (15) behave qualitatively similar.

The next task is to relate $p_{*}$ to the shear strength ${\tau }_c$ so that the shock propagation velocity $\omega$ can be computed. The GOEM disregards the shear stresses at the interface so that Equation (14) simplifies to ${\tau }_{\max }=p_{\alpha }/(2-2\nu )$ and one obtains $p_{*}=2(1-\nu ){\tau }_c$ in terms of the shear strength of the snow cover. However, the problem is more involved with shear stresses because there is no justification for assuming that the ratio between $p_{\alpha }$ and ${\tau }_{\alpha }$ remains unchanged across the shock. At this point, one needs to use the constraint from the jump condition for the interface-parallel momentum component:

$${\rho }_0·0·\omega -{\tau }_{*}={\rho }_{1}w(\omega -v)-{\tau }_{\alpha }.$$

(17)

If one assumes a uniform velocity profile and neglects the small component of $\mathit{\omega }$ parallel to the terrain, $v\approx \overline{u}sin\alpha$ and $w\approx \overline{u}cos\alpha$, implying that the jump in the shear stress accelerates the eroded mass to the mean avalanche speed as in the TJEM. ${\tau }_{*}$ (not to be confused with ${\tau }_{*}$ in Grigorian’s modified friction law, see Section 4) indicates the critical shear stress along the interface; it is related to $p_{*}$ and ${\tau }_c$ by

$${\tau }_c^2=b^2{p}_{*}^2+{\tau }_{*}^2.$$

(18)

The mass jump condition (8) allows to transform (17) into

$$\omega =\Theta ({\tau }_{\max }-{\tau }_{*})\frac{{\tau }_{\alpha }-{\tau }_{*}}{{\rho }_{0}w}.$$

(19)

The shock speed $\omega$ has to be the same, whether calculated from Equation (12) or (19), requiring

$$p_{\alpha }-p_{*}=({\tau }_{\alpha }-{\tau }_{*})tan\alpha .$$

(20)

$p_{\alpha }$ and ${\tau }_{\alpha }$ are given by Equations (15) and (16) and related to the maximum shear stress by (14). This leads to a quadratic equation with the solution

$$p_{*}=\frac{p_{\alpha }-{\tau }_{\alpha }tan\alpha }{1+b^2{tan}^2\alpha }\pm \frac{\sqrt{{\tau }_c^2(1+b^2{tan}^2\alpha )-{(p_{\alpha }-{\tau }_{\alpha }tan\alpha )}^2(b^4{tan}^4\alpha +2b^2{tan}^2\alpha )}}{1+b^2{tan}^2\alpha }.$$

(21)

To find the approximate conditions for the existence of a physical solution, assume the erosion angle $\alpha$ to be small—measurements at Vallée de la Sionne point towards typical values of 0.01 or less. Then, $p_{\alpha }\approx p_b$ and ${\tau }_{\alpha }\approx {\tau }_b$. Furthermore, the ratio ${\mu }_{\mathrm{eff}}\equiv {\tau }_b/p_b$ is typically in the range 0.3–1. A physical solution to Equation (20) thus only exists if ${\tau }_{*}\gtrsim 0$, $p_{*}\lesssim {\tau }_c/b$, ${\tau }_{\max }>{\tau }_c$, and $p_{\alpha }\approx {\tau }_c/b+{\mu }_{\mathrm{eff}}{p}_{b}tan\alpha$. Hence, a rough estimate for the admissible range of the hydrostatic pressure is

$$\frac{{\tau }_c}{b}\frac{1}{\sqrt{1+{({\mu }_{\mathrm{eff}}/b)}^2}}<p_b<\frac{{\tau }_c}{b}\frac{1}{1+{\mu }_{\mathrm{eff}}tan\alpha }.$$

(22)

The lower limit of $p_b$ is typically (0.9–1.5) ${\tau }_c$, the upper limit (1.2–1.7) ${\tau }_c$.

If $p_b$ falls below the lower limit, entrainment ceases. What happens if $p_b$ exceeds the upper limit? On the one hand, the conditions for erosion, that is, the propagation of the shock, are still fulfilled and $\omega$ obtained from Equation (12) increases with $p_b$. On the other hand, ${\tau }_b$ is too small to satisfy Equation (19). This means that the eroded mass cannot be accelerated to the full avalanche speed at once and forms a bottom layer that is dragged along by the main flow above.

5.5. Can the Grigorian–Ostroumov Erosion Formula Be Used in Practice?

The importance of modeling entrainment has been acknowledged among avalanche dynamics researchers for half a century, but it will probably still take years until flow models with entrainment (and deposition) become a standard in avalanche hazard mapping procedures. There are both a fundamental modeling problem and a number of practical issues to be overcome.

On the practical side, the following challenges persist: (i) Data suitable for testing entrainment models is still scarce. (ii) Entrainment models can be validated only in conjunction with a specific avalanche rheology or bed friction law. (iii) Many authors argue that entrainment models just add to the list of poorly constrained initial conditions and friction parameters (see, e.g., Reference [22]). First, published data [27,32] already allows for at least qualitative validation. The situation can be further improved by pooling observations and data from different sites to create a sound basis for quantitatively validating dynamical models. The second problem can be approached by implementing the EGYEM, GOEM and TJEM in flow models based on different friction laws or rheologies and comparing their performance with and without entrainment for selected test cases. Thirdly, a practical method needs to be developed for prescribing the erodible snow depth and the shear strength as functions of the return period, regional precipitation data, winter temperature, slope orientation, and the dominant wind directions. This task appears feasible where distributed data of temperature, snow depth and precipitation are available; however, additional measurements of snow shear strength at different altitudes and under different topographical conditions might also be required.

The first theoretical problem is to understand how the snow-strength parameter $p_{*}$ in the GOEM, the limit shear strength ${\tau }_{*}$ in the GFL, ${\sigma }_{*}$ in the EGYEM and ${\tau }_c$ in the TJEM are related to each other and to measurable properties of the snow cover. At this point, only a tentative answer can be offered. First, note that all four process models for friction or erosion originally refer to the top layer of the substrate (the new-snow layer in the case of snow avalanches). However, a weak layer underneath the top layer or poor bonding between the layers will greatly influence the behavior of the system; with regard to the erosion models, see the discussion below. In the GFL and TJEM, shear stresses are transmitted across a competent top layer to the weak layer, which will fail before the top layer. Thus, the shear strength of the weakest layer is decisive. In the original GOEM, the propagation of the primary shock front is determined by the strength $p_{*}$ of the top layer, but the normal stress $\min (p_{\alpha },p_{*})$ is transmitted across the top layer and may induce failure of a weak layer or interface underneath. Describing this situation is, however, beyond the scope of the GOEM. Perhaps surprisingly, the weak-layer shear strength does not contribute to the momentum jump condition in the EGYEM if the shock is infinitesimally thin. This implies that the strength of the top layer is the decisive parameter in the EGYEM whereas the weak layer determines the erosion depth. In a refined version of the EGYEM with a finite-width transition zone, the weak-layer strength will contribute to some degree.

A second point to note is that the parameters of all four process models are failure criteria for a solid under given normal and shear stresses. All models implicitly assume that there are confining lateral stresses (only in the spanwise direction for the EGYEM). This implies that both the normal and the shear stress have to be taken into account to determine the maximum shear stress, as shown in Equation (14). The failure criterion is ${\tau }_{\max }$ exceeding the shear strength of the snow cover. The GOEM, TJEM and GFL neglect one of the two contributions to ${\tau }_{\max }$, but this can easily be improved. The shear strength of snow can, in principle, be determined as half the difference between the major and minor normal stress at failure in a triaxial compression test—a standard procedure in soil mechanics. However, with snow these tests have to be conducted at high compression rates to enter the brittle-failure regime.

The second theoretical question in entrainment modeling concerns the layered structure of natural snowpacks: Profiling-radar measurements at the test site Vallée de la Sionne in 1999 gave evidence for continuous, scouring-type erosion and entrainment ([33], Figure 3) as described by the GOEM or TJEM. However, later measurements at the same locations frequently show step-wise (“ripping” [34]) erosion [27,28] where the erosion front at the radar location quasi-instantaneously jumps, presumably either because the layer is weak or because a weak layer underneath fails. In general, one has to expect that both erosion modes coexist in the same event. We conjecture that scouring occurs when the new-snow shear strength increases with depth and the layer bonds well to the underlying snowpack. With step-wise erosion, the eroded slab is too massive to be entrained immediately. It is compressed, comminuted and gradually accelerated while the erosion front advances with the avalanche front.

A simple approach to this problem is to disregard the weak layer and to use a model for scouring erosion, hoping that the effects of delayed entrainment are small on the scale of the entire flow. Alternatively, one may envisage an extended model that treats the new-snow layer as dynamical. Where the shear stress at the interface to the old snow exceeds the weak-layer shear strength ${\tau }_{w,c}$, fracture is assumed and ${\tau }_{w,c}$ is replaced by a suitable friction law. The slab is accelerated by the difference in the shear stresses from the avalanche and the substrate and by the longitudinal pressure gradient. The mass-exchange term between the two layers should describe instantaneous entrainment of a slab segment when the relative velocity difference is nearly zero. The model should also account for possible scouring entrainment from the moving slab into the avalanche. Transforming this concept (Figure 4) into a consistent set of equations and a working model is left for future work.

The time appears right for a new attempt to solve the problem of entrainment and deposition—Sovilla et al. [35] found inclusion of entrainment to improve the modeling results when compared to observations and to substantially reduce the unphysical scatter of friction parameter values needed to reproduce the observed run-out distances. Recently, Rauter and Köhler [36] demonstrated that two rather dissimilar avalanches released in the same path on the same day could be modeled with the same parameter set, provided that the avalanche model accounts for both entrainment and deposition. The GOEM (with the modifications discussed above and preferably an extension for dynamically calculating the avalanche density) provides a physically consistent description of continuous, scouring erosion of a homogeneous snowpack without empirical parameters. Being one of the few serious contenders in this field, it clearly deserves further study.

6. Do Avalanches Behave Non-Monotonically in the Friction Coefficients?

Probably the most surprising finding in Reference [1] is the non-monotonic behavior with regard to variations of the dry-friction parameter $\mu$ under otherwise identical conditions: The avalanche “ignited” only in a limited range of $\mu$, both with larger or smaller friction coefficients entrainment soon ceased and the flow stopped. Such behavior was observed only in the complete model and only if both GFL and the GOEM were activated [18]. If such behavior were to generally occur in models with non-smooth friction law and a threshold for entrainment, the practical usefulness of advanced, complex models would be put in question. One way to investigate this problem is to reimplement the “complete” model as a numerical code to confirm that the observed non-monotonous behavior is not caused by numerical artifacts and to explore the phase diagram of the model.

Grigorian and Ostroumov did not observe non-monotonic behavior in the “simplified” models, but those models do not include entrainment, which seems to be an essential prerequisite for this phenomenon to occur. It may therefore be instructive to study the “simplified” model enhanced by an entrainment term derived from the GOEM because it has only two degrees of freedom—the maximum of the instantaneous flow depth, H, and the front velocity U—in contrast to the infinite number of degrees of freedom contained in the fields $h(S,t)$ and $u(S,t)$ of the “complete” model. One can thus study such an “enhanced simplified” model with the techniques developed in the mathematical theory of dynamical systems. In principle, this can be done for any given path profile described by $\psi \left(S_f\right)$, but in practice one will select a few elementary types like an inclined plane or the hockeystick profile with slope angle ${\psi }_1>{\psi }_0\ge 0$.

Such a study remains yet to be done, but we briefly derive the set of equations here. The mass balance, initially expressed as an integral over space and time,

$$F\left(t\right)=F_0+{\int }_0^t{\int }_{S_0}^{S_f\left(t^{\prime }\right)}{w}_e\left(p(S,t)\right)\mathrm{d}S\mathrm{d}{t}^{\prime },$$

(23)

transforms into a spatial evolution equation for the longitudinal section area $F\left(S_f\right)=\varkappa (S_f-S_0)H\left(S_f\right)$ with $F\left(S_{f,0}\right)\equiv F_0=\varkappa (S_{f,0}-S_0)H_0$ if one uses the front position $S_f$ as the independent variable instead of the time and sets $\mathrm{d}X/\mathrm{d}t=U\mathrm{d}X/\mathrm{d}{S}_f\equiv U\left(S_f\right)X^{\prime }\left(S_f\right)$:

$$H^{\prime }=\frac{1}{S_f-S_{f,0}}\left(\frac{1}{\varkappa U}{\int }_{S_{f,0}}^{S_f}{w}_e\left(p(S,t)\right)\mathrm{d}S-H\right).$$

(24)

The equation of motion reads

$$U{U}^{\prime }=gsin\psi \left(S_f\right)-\mathrm{sgn}\left(U\right)f(U,H,\psi \left(S_f\right))-w_e\left(p\left(S_f\right)\right)\frac{U}{H}.$$

(25)

The appearance of an integral in (24) but not in (25) is a direct consequence of the assumptions that (i) the shape of the longitudinal section of the avalanche remains self-similar and (ii) only the front is relevant for the dynamics of the flow. To evaluate Equation (24), some further simplifying assumptions are required. Grigorian and Ostroumov mention that the shape of the avalanche is well described by the parabola

$$h(S,t)/H\left(t\right)={(S/S_f\left(t\right))}^2,$$

(26)

which leads to the value $\varkappa =1/3$ for the shape parameter, close to the recommended value 0.36. If there were no entrainment, self-similarity of the shape through time would imply

$$u(S,t)/U\left(t\right)=S/S_f\left(t\right).$$

(27)

In the presence of entrainment (which diminishes with distance from the front), the velocity must be smaller than the value given by Equation (27) so that part of the eroded snow is transported towards the tail to maintain the parabolic shape. However, for the purpose of investigating the qualitative properties of this dynamical system, it seems justifiable to use Equation (27) as an approximation. The integral over the erosion rate can then be evaluated in terms of H and U, the resulting expression depending on the one choice for the pressure $p\left(H\right(t),U(t\left)\right)$, as discussed in Section 5.4.

The two Equations (24) and (25) determine the trajectories of the system in the phase space $(H,U)$. Of particular interest are its fixed points, that is, pairs of values $(H_{*},U_{*})$ for which $H^{\prime }(H_{*},U_{*})=U^{\prime }(H_{*},U_{*})=0$, and critical lines with $H,U\to \infty$. Fixed points $(H_{*},0)$—corresponding to avalanches that stop—can exist only if $tan\psi \left(S_f\right)<\mu$ over a suitable interval in $S_f$. In contrast to the classic examples of dynamical systems found in textbooks, the Grigorian–Ostroumov model on a hockeystick profile with slope angles ${\psi }_1>{\psi }_2>0$ has a continuum of such fixed points, with their basins of attraction degenerated to a single line.

In Reference [37], a rather general asymptotic solution to a simplified model for an entraining flow on an infinite inclined plane was found; it corresponds to a critical line tending to a straight line parallel to $U=cH$. A characteristic property of that model is that the bed shear stress is limited to the shear strength of the snow (the TJEM is a limiting case with the depth of the shear layer shrunk to 0). It remains to be seen whether such a solution also exists for the Grigorian–Ostroumov model.

Once the structure of the phase space with its fixed points, critical lines and basins of attraction is determined for given values of $\mu$, ${\tau }_{*}$, k, and $p_{*}$, one may proceed to study the dependence of this structure on $\mu$, with the other three parameters held fixed. Figure 5 visualizes the dependence of the phase-space structure for different values of $\mu$ in a strongly simplified way. The non-monotonicity found in Reference [1] would manifest itself in the way that some starting point $(H_0,U_0)$ in the phase space lies within the basin of attraction of some fixed point $(H_{*},0)$ for $\mu <{\mu }_1$ and $\mu >{\mu }_2$, but belongs to the basin of attraction of another fixed point $(H_{*}>0,U_{*}>0)$ or tends towards a critical line with $H,U\to \infty$ if ${\mu }_1<\mu <{\mu }_2$.

Almost half a century after the main work by Grigorian and Ostroumov summarized in Reference [1] and close to a quarter century since the writing of that paper, our insight into the flow behavior of snow avalanches has been dramatically deepened through new experiments, and a large number of new numerical models have been developed and are being applied in practical work. Yet, the span of topics discussed in the present Comments paper, the open questions raised by the Grigorian–Ostroumov model and the new avenues for research that a fresh reading of their paper suggests are proof enough that this work can inspire the work of today’s generation of avalanche researchers. The same holds for early and more recent work summarized in the review paper by Eglit et al. [17], which extends the scope of Reference [1].