On a Continuum Model for Avalanche Flow and Its Simplified Variants

Mathematical models of different degrees of complexity, describing the motion of a snow avalanche along a path with given center line and spatially varying width, are formulated and compared. The most complete model integrates the balance equations for mass and momentum over the cross-section and achieves closure through an entrainment function based on shock theory and a modified Voellmy bed friction law where the Coulombic contribution to the bed shear stress is limited by the shear strength of the snow cover. A simplified model results from integrating these balance equations over the (time-dependent) length of the flow and postulating weak similarity of the evolving avalanche shape. On path segments of constant inclination, it can be solved for the flow depth and speed of the front in closed form in terms of the imaginary error function. Finally, the very simplest model assumes constant flow height and length. On an inclined plane, the evolution of flow depth and velocity predicted by the simplified model are close to those from the full model without entrainment and with corresponding parameters, but the simplest model with constant flow depth predicts much higher velocity values. If the friction coefficient is varied in the full model with entrainment, there can be non-monotonous behavior due to the non-linear interplay between entrainment and the limitation on the Coulomb friction.

The present paper was originally written in 1996 for the proceedings of a symposium held in Davos, Switzerland, in late 1996, but that volume-and thus this paper-were never published due to circumstances recounted in Section 1 of the Comment paper [1] written by the editor. In the editor's opinion, it still contains a host of valuable material that merits publication, but the manuscript needed to be edited to conform to the standards of a refereed journal in 2019. Since the lead author died in 2015, he could not do this himself or comment on the (substantial) changes made by the editor. Therefore, the original text and figures are reproduced in Section 1 of Comment paper [1]. In the original manuscript, all figures are contained in an appendix, without captions. The most essential ones were inserted into the main text and given captions; Section 2 of Comment paper [1] contains the remaining figures. Besides the history of the paper, Ref [1] contains notes that relate to references in the present paper and explain why certain changes were made or expand on points that may be difficult to understand for readers not familiar with the earlier work of the pioneers of avalanche research in the Soviet era. Two particular aspects of this paper-the value of analytic solutions and the modified friction law-are also briefly discussed in [1].

Introduction
The mathematical modeling of snow avalanche movement is a basic task when estimating and predicting possible hazard parameters for given climatic and geographical conditions ( [1], Note 1). There are three main topics arising in the quantitative description of snow avalanches: (1) the physics and mechanics of the snow cover on mountain slopes and its evolution in time, potentially culminating in the release of a snow avalanche; (2) the motion of snow masses on mountain slopes; (3) the impact of moving snow on obstacles and structures. The mathematical modeling of these processes requires application of different approaches and methods based on suitably simplified descriptions of these complicated phenomena. A comprehensive survey of essential progress in this field and an analysis of necessary future research was given years ago in [2]. Since that time, many new results and publications have appeared, but there is still need for more efficient and practically applicable methods in this field.
The topic of the present work is to present the state-of-the-art in avalanche dynamics by discussing mathematical models for evaluating the main parameters of snow avalanche motion at different degrees of simplification with regard to physics content and mathematical formulation. Such a suite of models is necessary and useful: In situations where one does not have detailed information about the slope and the snow cover parameters, etc., one needs an instrument for quick, but rough estimations. In contrast, a detailed quantitative description is desirable when this information is available. In all cases, constructing a mathematical model of moving snow masses implies significant simplification by introducing averages of dynamical quantities and effective values of physical and mechanical properties of flowing snow and the underlying bed. Of course, the models need to be tested-and improved if necessary-by comparing the simulation results to data collected from real avalanche events.
The problems of snow cover evolution and snow flow impact are not touched upon here. The mathematical model of snow avalanches presented here was developed at the Department of Mechanics of Natural Processes at the Institute of Mechanics of Moscow State University.

The Full Model
The mathematical problem of modeling the release and movement of snow masses on a mountain slope is rather complicated and ill-posed due to our lack of knowledge about the failure processes and gradual fragmentation of snow in the early phase of an avalanche event, as well as corresponding open questions about the mechanical interactions between the snow fragments and particles and the underlying soil or rock bed. To a first approximation, which is ordinarily sufficient for practical needs, we can consider a rough model that averages the main dynamical variables, in the flow direction over a scale much larger than the typical bed roughness and over the full local cross-section of the flowing snow mass.
Such a procedure leads to so-called hydraulic models, which are widely used in engineering hydraulics and give adequate results for practical purposes. So, it is reasonable to hope that the hydraulic approach will also be appropriate when modeling the dynamics of snow avalanches.
Following this approach, we shall use the "full" mathematical model [3] described by the differential equations ∂F ∂t where sign(u) = In deriving the system of Equations (1) and (2), we have assumed that the snow masses move in a channel with rectangular local cross-section that varies with the longitudinal coordinate along the channel axis. The structure of (1) also applies if the geometry of the channel cross-section varies but is prescribed as a function of the longitudinal coordinate. Here, t, S are the time and space coordinate along the avalanche trajectory on the bed (see Figure 1); u(S, t), F(S, t) are the unknown values of the S-component of snow velocity averaged over the full cross-section of the flow and the area of this cross-section, respectively; L(S), ψ(S) are the local width of the cross-section and local inclination of the channel axis; g is the gravity acceleration, a(S, u) is the projection of the full acceleration (gravitational plus centrifugal) of the snow "particle" onto the direction normal to the bed. Recall that u is averaged over the cross-section and along the S-coordinate with averaging scale larger than the depth of the flow. k is the hydraulic resistance coefficient; R = FL/(L 2 + 2F) is the so-called hydraulic radius, and r is the curvature radius of the channel axis. The values of f 1 , f 2 , q are specified as follows. We follow Grigorian [4] and Grigorian and Ostroumov [5] for the friction force f 1 on the bed surface and adopt the expression ( [1], Note 2) where µ is the Coulomb friction law coefficient, ρ is the density of the flowing snow, h is the local flow depth, and τ * is the upper limit of the friction shear stress that appears in the modified friction law of Grigorian [4], expressed as ( [1], Note 3) with p n the normal pressure on the sliding surface. The quantities f 2 and q relate to entrainment of the (non-moving) snow cover at the frontal surface of moving snow. This process is difficult to model mathematically, and our approach is as follows. We suppose that the entrainment rate, represented by q, depends on the load p exerted on the snow cover by the moving snow. The kinematics of the erosion process are shown in Figure 2, where δ 0 is the thickness of the (erodible) snow cover, δ is its thickness in the zone where that snow is gradually being entrained, and ω is the propagation velocity of the front where new snow is being eroded and mixed into the moving snow mass. Using the momentum and mass conservation laws, we can write ( [1], Note 4) where ρ 0 and ρ 1 are the densities of undisturbed new snow and eroded snow; ω and v are the interface-normal velocity components of the erosion front and of avalanching snow particles, p * is the strength of the snow cover, and p is the full pressure (hydrostatic plus dynamic) exerted on the bed (new snow) surface. Furthermore, we shall consider the case ρ 1 > ρ 0 . Figure 2. Entrainment near the front of the avalanche. The interface between the snow cover (with original depth δ 0 (S) and instantaneous depth δ(S, t)) and the flowing avalanche is regarded as a shock front that propagates into the snow cover with normal velocity ω.
The relations (5) and (6) lead to with the parameter σ characterizing the degree of snow compaction during the erosion process. The pressure p is expressed by the relation where h(S, t) and δ(S, t) are the local depth of moving snow and the local thickness of the erodible snow cover, respectively; C is an empirical constant. Using (7), we get the expression for the entrainment rate q. The average density ρ is assumed to differ from ρ 1 in general. These considerations also lead to the following expression for the force contribution f 2 ([1], Note 5): where the condition δ = 0 means that the new-snow layer is completely eroded and entrained at a certain distance behind the avalanche front. An alternative expression for q can be derived by connecting the snow erosion rate to the shear forces, leading to an empirical relation of the type This alternative formula can easily be implemented in the mathematical model, but we will not discuss this in detail here.
Finally, we have the set of relations (1)-(3), (8)-(10), giving the closed mathematical formulation of the model. Before they can be used to solve a given problem, it remains to formulate the initial and boundary conditions.
The initial conditions are where the functions u 0 (S), F 0 (S), and δ 0 (S) are to be specified for each case under consideration. The boundary conditions describe the motion of the front and the "tail" part of the snow flow. Even though the moving snow mass not only has a sharp frontal boundary but also a tail moving down the mountain side, we suppose that the latter boundary has limited importance for the dynamics of the main snow mass. Therefore, we adopt a simpler scheme with the tail of the avalanche at rest. In other words, we will set the tail boundary conditions as where S 0 is the (constant) coordinate of the avalanche tail.
The boundary conditions at the avalanche front are evident: where S f (t) is the coordinate of the front-a function of t to be determined in the course of solving the problem. This completes the full mathematical formulation of the initial-boundary value problem of snow avalanche dynamics. Now we consider a special case of the problem where the channel is infinitely wide (L → ∞), making the motion effectively two-dimensional, and the slope is devoid of erodible snow. Such a simplification of the problem allows to obtain results in a form suitable for analysis and parametric comparisons, which serve as guideposts in more complicated cases.
In this case, we have the differential equations ∂h ∂t with S ∈ [S 0 , S f (0)], and the boundary conditions Here, again, t and S are the time and the space coordinate along the avalanche path (see Figure 3), S f (t) is the coordinate of the avalanche front, etc.; h(S, t) is the local depth of the avalanche; u 0 (S), h 0 (S) are the initial distributions of u(S, t) and h(S, t), respectively; f is the retarding acceleration given by the relation where again µ is the Coulomb friction coefficient, ρ is the density of flowing snow, and τ * is the upper limit of the frictional stress τ on the bed surface, as specified in the modified friction law proposed in [4]. In the mathematical formulation of the problem (15)-(20), the set of parameters ψ(S), u 0 (S), h 0 (S), µ, k, ρ, τ * , is assumed known. The function ψ(S) describes the geometry of the slope and is known if the topography of the slope is available. The other quantities in this parameter set have to be estimated from the mechanical properties of snow, which depend on climatic and geographical factors and-in the case of u 0 (S) and h 0 (S)-on the probable and admissible initial distributions of u and h. The problem formulated in terms of the full model can usually only be solved numerically on a computer ([1], Note 6). This is recommended for "precise" calculations when detailed information about the input data (ψ(S), k, etc.) is available and the mathematical predictions are to be calculated as exactly as possible. The formulation (15)-(20) provides the necessary tools for such work.
However, in many real cases, there is insufficient information available on ψ(S), µ, and other input quantities. Moreover, one often needs a simpler mathematical instrument for estimates-either an explicit formula or a simple computational procedure that can be applied without a powerful computer. For such purposes, we can construct different classes of simplified models, and those derivations are the topic of the following parts of this work.

A Simplified Version of the Full Model
A possible simplification procedure, considered initially in [6], is as follows. On the right-hand side of (1), we neglect the third and last term and rewrite the resulting equation in a Lagrangian coordinate system as where ξ is the Lagrangian coordinate. This transforms (1) into 1 2 The second step of the simplification procedure introduces the integrated form of the mass conservation law instead of its differential form (16). Integration of h(S, t) over the coordinate S leads to the relation κ is a parameter characterizing the"fullness" of the graph of h(S, t). Now we make the crucial and plausible assumption that κ is constant ([1], Note 7). Of course, κ changes over time but, as will be demonstrated by comparing numerical results obtained from the full and simplified models, the deviations of κ from a constant value are not dramatic; therefore, this hypothesis looks acceptable for estimations. To evaluate the model (21)-(23) with simple calculation tools, one can divide the avalanche path into several segments, introducing the dividing points S 0 , S 1 , S 2 , . . . , S N so that changes of ψ(S) can be neglected within each segment S i ≤ S ≤ S i+1 . This allows us to integrate (22) over [S i , S i+1 ], where we get the relation (we consider only the case u > 0, sign(u) = 1). To calculate the velocity u from (24), we need information about h to evaluate (25). As calculations with the full model demonstrate, the important unknowns of avalanche motion-velocity, depth, and runout distance-are mainly determined by the flow properties in the zone near the coordinate with the maximum value of h. The full-model calculations show that this region is closely attached to the avalanche front and that the Lagrangian coordinate of the phase with h(S(ξ)) = max S (h(S, t)) changes little over time. For this reason, we can further simplify the calculations by supposing that this coordinate remains unchanged in time and coincides with the Lagrangian coordinate of the avalanche front. All these considerations allow us to use the value H = max(h) instead of h and calculate l by the relation substituting S for S f because S-coordinates in the dynamically most significant zone near the cross-section with H ≈ max(h) are close to S f . In [6], values in the range κ ≈ 0.3-0.4 are recommended in cases where the geometry of the slope changes smoothly with S. Of course, to assess the accuracy of the simplified model presented here, one needs to compare parallel calculations for the same problem using the full and the simplified model. This was done, and the results are presented in the next section.

Comparison of Computations with the Full and Simplified Models
The calculations with the full model were performed in a simple geometry with ψ = const. On the interval [S 0 = 0, S f (0) = 100 m], different initial distributions h 0 (S), u 0 (S) were prescribed, and different combinations of the values of the parameters µ, k, ρ, τ * , and ψ were tested. Table 1 details the parameter values for the different variants. The basic variant was ψ = 30 • , µ = 0.5, k = 0.02, ρ = 500 kg m −3 , τ * = 10 kPa ([1], Note 8). h 0 (S) was prescribed as triangular or parabolic forms with different values of H 0 ; for u 0 (0), a uniform distribution was used (see Figure 4). Figure 5a presents the distribution of h/H as a function of S in different stages of the motion. Figure 5b plots the evolution of κ, H, w, and S * /S f against S f , where S * is the S-coordinate of the point where h = max(h) = H. All these data relate to the basic variant. Figures S2.1-S2.9 in SMD 2 illustrate analogous results for other variants.  The common feature of the solutions for all the variants is their fast approach towards an asymptotic form, independent of the initial distribution. The asymptotic distribution of h can be described approximately by the relation h/H = (S/S f ) 2 , and H, w and S * /S f become smooth and slowly varying functions of t with S * /S f → 1, κ → 0.36 for all the variants calculated.
Note that the asymptotic behavior of the solution emerges after the avalanche body has spread to the point where H < h * ≡ τ * /(ρ a µ), so that the Coulomb friction law governs the entire avalanche body. Before this stage, significant changes in κ and S * /S f are visible, and the character of these changes depends essentially on the initial data.
We can conclude that the asymptotic simplification of the flow structure with κ = const., S * /S f = const. ≈ 1 will occur in cases with slowly changing slope geometry, and for this very case the simplified model is acceptable for simple estimations. In Table 2, the results of calculations with the full and simplified model are presented in terms of w and H as functions of S f . One can see that the simplified theory gives acceptable accuracy for H and overestimates the values of the front velocity w. This last effect is mainly due to our neglecting the resistance term −ku 2 /h in the equation of motion of the simplified model and indicates the necessity of improving the simplified theory. Such an improvement will be presented in the following section.  Table 2. The results of calculations with the full model (Section 2), the simplified model of Section 3 (formulas (23) and (24)), and the improved simplified model described in Section 5 (formulas (36), (39) and (40)) are presented in terms of the flow depth H and the front velocity w as functions of the front position S f . For the simplified models, the value of the shape factor was always chosen as κ = 0.36 ([1], Note 9).

Improvement of the Simplified Model
Including in the momentum equation the collisional term −ku 2 /h and using all other simplifying considerations introduced in the previous section, we get the relations 1 2 Introducing the definitions we can rewrite (27) in the form The solution of (32) on the segments S ∈ [S i , S i+1 ] with initial condition is given by the relation ([1], Note 10) When calculating the coefficients A, B, and D, one should use the average of ψ over the interval Defining the function and setting ξ ≡ S √ B, we get The graph of Z(ξ) is shown in Figure 6. Consider now the qualitative features of solutions of (32). Using a new variable y, we have, instead of (32), the equation and as u 2 = y − D/B ≥ 0, we will be interested only in solutions of (38) satisfying the conditions y ≥ D/B ≥ 0. On the y-S plane, the hyperbolas yS = γ = const. ( [1], Note 11) are the isoclines of (38), i.e., the lines along which y = const. On the hyperbolas we have y = 0, so all the solutions have extrema there (more precisley, these are maxima, for if γ < A/B then y > 0 and at γ > A/B y < 0). As a result, we can visualize the integral curves of Equation (38) as shown in Figure 7. Depending on the initial data, when S → ∞ an integral curve tends to the hyperbola (39), yS = A/B, either monotonically or after some initial rise to a maximum value of y. It can easily be shown that y = 0 on the curve and that, at any given point of (40) with coordinate S = S * , a solution of (38) touches some hyperbola yS = γ * depending on the initial data. For S > S * , this solution becomes "captured" between the curve (40) and the hyperbola yS = γ * . This leads to simple estimates for large values of S: On the basis of this improved model, a series of calculations have been performed with the same parameter values as for the non-improved model, and the results are shown in Table 2, illustrating the performance of the improved model compared to the "exact" results of the full model. The correspondence between these two sets of results is very good over a wide range of the key parameters of the problem. At large values of S, even the simplest relation (39) delivers estimations of high accuracy. Note that all these calculations were performed for κ = 0.36. This value, thus, can be recommended for practical use over a wide range of parameters.

The Simplest Model for Rough Estimations
In the literature, one can find many articles that use a very simple scheme for snow avalanche dynamics. In this scheme, the avalanche is considered as a material point, neglecting the spatial extent of the flowing snow mass and its changes in time. Here, we will consider the problem in this case and establish what kind of assumptions are needed to derive such a model from the "exact" one. The main differential equations of the material point model are as follows: where m is the total mass, which is assumed to be constant, r is the resistance force between the moving body and the underlying base surface, and other denotations are as before. We can write the relation where ρ is the snow density (assumed constant), H, l, b are the height, length, and width of the avalanche body, respectively. For the resistance r, we have where τ is the average frictional stress on the bed surface. So, we reduce (42)-(44) to and it is obvious that, to derive Equation (45) from the corresponding momentum Equation (1) of the full model, we must neglect the term representing the gradient of normal force acting at the cross-section of the avalanche body ([1], Note 12), (ρg/2F)∂F/∂S, and include the hydrodynamic resistance term, −kρu 2 /R, in the expression for τ, if desired. Up to this point, there is no difference between the simplified model of Section 3 and that under consideration here. Such a difference appears when we look at the next simplification in this simplest point-mass model, which assumes the flow depth H to be constant. As calculations with the full or simplified model reveal, H varies significantly in the course of the avalanche propagation along the mountain slope. One may, nevertheless, hope that it is admissible for rough estimations to use H = const. Consider now the evaluation of the simplest model for several hypotheses about the bed shear stress τ.

Solutions for the Coulomb and Voellmy Friction Laws
In point-mass model calculations ( [1], Note 13), often the Coulomb friction law is used, with τ expressed as in which the hydrodynamic term (quadratic in u) is neglected. In this case, the problem is solved by quadrature: where S 0 , u 0 are the initial position and velocity of the avalanche body.
For a given slope geometry, the right-hand sides of (47) and (48) are known functions of S (as ψ(S) is known). So, the numerical evaluation of these expressions reduces to rather simple calculations. Note that if the slope is steep enough, i.e., sin ψ − µ cos ψ > 0, the avalanche motion will accelerate in time. In the opposite case, the motion decelerates and even stops on the slope if the inclined part of the trajectory is long enough.
To obtain explicit formulas for analyzing the avalanche behavior, we consider an avalanche path consisting of two segments, the first inclined at an angle ψ 1 = const. and the second horizontal (ψ 2 = 0). Using (47) and assuming u 0 ≥ 0, we have ( [1], Note 14) With l 1 ≡ S 1 − S 0 , where S 1 is the coordinate of the end of the first slope segment, we get the initial velocity in the second segment as Using (47) in this segment with (50), instead of u 0 , and substituting ψ 2 = 0, we have ([1], Note 15) In the second segment, the avalanche decelerates and completely stops at some distance S 2 . So, for this distance, we have (substituting u = 0, S = S 2 in (51)) Note that this runout distance l 2 does not depend at all on the avalanche size. Using (50) and (52) and the definition of l 1 , we obtain the full runout distance L = S 2 − S 0 as ([1], Note 16) If sin ψ 1 − µ cos ψ 1 < 0, the avalanche may stop on the slope; in the opposite case, it runs out on the horizontal path segment.
We can now easily extend this simplest model with a resistance term that is quadratic in the velocity. In the inclined segment of the path, we have (assuming u 0 ≥ 0) ( [1], Note 17) instead of (49) and (50). On the horizontal path segment, we find So, if we set u = 0 here, we get instead of (52). Finally, we have for the full runout distance ( [1], Notes 16 and 18) When k → 0, Equation (54)-(57) reduce to (49)-(53), as they should.

Solution with the Stress-Limited Friction Law
Consider now the case of very big avalanches for which the snow depth H almost everywhere in the avalanche body exceeds the critical value H * = τ * /(µρ cos ψ) beyond which the friction law τ = τ * = const. is applicable instead of the Coulomb law ([1], Note 13). In this case, after integrating (44), we have the following relations in the initial segment of the avalanche path: At the mountain foot, at S − S 0 = S 1 − S 0 = l 1 , we have and For the second (horizontal) part, we have and at S − S 1 = S 2 − S 1 = l 2 , where u = 0, we obtain The full distance L = l 1 + l 2 is given by the formula ( [1], Notes 16 and 19) From (68), we obtain (for τ * = 10 kPa) Reducing τ * to 5 kPa results in l 2 /l 1 ≈ 0.92/8. These numbers show that the influence of the quadratic resistance term is strong enough to reduce the runout distance significantly. In both cases, the reduction of the runout distance due to the quadratic resistance seems to be highly overestimated with k = 0.02 ( [1], Note 22). The problem of how to adequately estimate the hydraulic part of the flow resistance remains open and needs additional research.

Influence of Snow Entrainment on Avalanche Dynamics
We now return to numerical examination of the full model to see how entrainment of new snow at the front of a propagating avalanche influences its dynamics. We have carried out a series of computations with the full model for a simplified setting with a new-snow cover of constant thickness of 1 m (measured perpendicular to the surface) on an inclined plane. Starting from the basic parameter set µ = 0.25, k = 0.1, τ * = 4 kPa, p * = 4 kPa, C = 1, σ = 1, ρ = 300 kg m −3 , and ρ 0 = 300 kg m −3 , we systematically varied the parameters one by one, as shown in Table 3. Except for simulation number 18 ( Figure S2.14) with a parabolic, 40 m long and up to 4 m high head followed by a 0.2 m thin tail, the avalanches started as 100 m long and 1 m deep slabs. Five selected simulations are presented in Figure 8 and 9 in terms of the distributions of h and u at the instance when the front reaches the point S = 500 m. Corresponding plots for the other simulations can be found in Section 2 of Comment paper [1], Figure S2.10-S2.14.  Table 3). (b) Same as (a), except for the flow density increased from 300 to 500 kg m −3 (case 2 in Table 3).
Depending on the combination of the values of the governing parameters, the avalanche behaves according to one of two quite different regimes: One is characterized by intense new-snow entrainment in the frontal zone and the formation of a huge avalanche head. In the other regime, the avalanche is unable to maintain new-snow erosion even if erosion had occurred in the starting zone, resulting in a "weak" avalanche with a spreading body and decreasing velocity. It is very important to note that the governing parameters (µ, k, etc.) determine in a non-monotonic fashion which of these two flow regimes will be realized. As an example, Figure 9 illustrates how, for 0.4 ≤ µ ≤ 0.55, the end values of the range lead to avalanche dynamics of the second, "weak" type, while the intermediate value µ = 0.42 is associated with dynamics of the first, "strong" type! ( [1], Note 23) This seemingly strange dependence of the mathematical solution of the avalanche dynamics problem is, of course, a consequence of the strong non-linearity of the modified friction law of Grigorian [4]. If the avalanche flow depth is sufficiently small, so that τ < τ * everywhere and snow accumulation in the avalanche body due to entrainment is weak, then the condition τ > τ * will never be fulfilled and the avalanche dynamics remains "weak" at all times. If, in contrast, such mass accumulation does occur in some zone (e.g., due to sharp changes in the cross-sectional area of the channel, locally large values of new-snow depth, etc.) and the condition τ > τ * is met, the avalanche behavior changes in accordance with the friction law, the avalanche "ignites" and becomes highly mobile. As a result, the output parameters characterizing the avalanche motion (mean velocity, thickness, runout distance, etc.) depend on the initial parameters in a stochastic and unpredictable way. From a mathematical point of view, this feature of avalanche dynamics is the most interesting, and it is very important for our general understanding of the avalanche problem. Table 3. Parameter values of numerical simulations with the full model for investigating the influence of entrainment on avalanche dynamics ( [1], Notes 8 and 24). Figure numbers preceded by 'S2' refer to Section 2 of Comment paper [1].

Conclusions
On the one hand, the analytical and computational investigations presented here show that the mathematical problem of quantitative description of snow avalanche dynamics is highly nonlinear and that simple explicit solutions are difficult to obtain. On the other hand, there is a need for simple explicit relations in practical applications. Unfortunately, this need cannot be satisfied to the desired degree for only rough estimations are possible on the basis of very "finite" simplifications made with more or less exact full models. The corresponding procedures are presented and illustrated by numerical calculations in this paper.
One of the most important parameters for practical applications is the maximum runout distance of avalanches at given geographical and climatic conditions. The maximum runout distance corresponds to the most powerful and massive avalanches, in which case one can expect the details of the full model not to be significant. A sufficiently accurate result can in this case presumably be obtained from the simplified model or even from the simplest model described above.
A more difficult problem is quantitative description of the evolving and subsiding avalanche motion, which depends on uncertain information relating to natural data (ψ(S), µ, ρ, τ * , etc.), and mathematical modeling of the process of snow entrainment into the avalanche body in the frontal zone. The computational simulation of the avalanche dynamics made with the full model shows that avalanche behavior is stochastic in some sense. Namely, depending on slight differences in the initial data, the bed geometry, etc., mass velocities, depth distribution and other parameters of motion exhibit finite differences and pulsating behavior as the avalanche propagates down the mountain slope. This phenomenon is associated with the dual nature of the assumed bed resistance law: the transitions between the Coulomb law and the τ * -law are connected to changes in the flow-depth distribution H(S), which in turn depends non-linearly on the resistance law. This stochasticity of avalanche dynamics prevents a deterministic description of the full problem; instead, only the statistical characteristics of the main avalanche parameters, as calculated by the (seemingly deterministic) full model, should be evaluated and used in practice.
Nevertheless, there is a program for improving the "exact" full model. It comprises tests of various schemes of snow entrainment at the avalanche front and more adequate modeling of the hydraulic or viscous parts of the resistance, as well as development of a complete software package that integrates all the components of numerical flow simulation and parameter variation that are needed to solve real-world problems. This work remains yet to be done.