2. 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
where
and
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);
,
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;
,
are the local width of the cross-section and local inclination of the channel axis;
g is the gravity acceleration,
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;
is the so-called hydraulic radius, and
r is the curvature radius of the channel axis. The values of
,
,
q are specified as follows. We follow Grigorian [
4] and Grigorian and Ostroumov [
5] for the friction force
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
the normal pressure on the sliding surface.
The quantities
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
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
and
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,
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
.
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
and
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
in general. These considerations also lead to the following expression for the force contribution
([
1], Note 5):
where the condition
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
,
, and
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
is the (constant) coordinate of the avalanche tail.
The boundary conditions at the avalanche front are evident:
where
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 (), 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
initial conditions
with
, and the boundary conditions
Here, again,
t and
S are the time and the space coordinate along the avalanche path (see
Figure 3),
is the coordinate of the avalanche front, etc.;
is the local depth of the avalanche;
,
are the initial distributions of
and
, 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
,
,
,
,
k,
,
, is assumed known. The function
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
and
—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 (
,
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 , , 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.
3. 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
The second step of the simplification procedure introduces the integrated form of the mass conservation law instead of its differential form (
16). Integration of
over the coordinate
S leads to the relation
where
,
, and
.
is a parameter characterizing the“fullness” of the graph of
. 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
,
,
so that changes of
can be neglected within each segment
. This allows us to integrate (
22) over
, where we get the relation
where
(we consider only the case
,
). 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
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
instead of
h and calculate
l by the relation
substituting
S for
because
S-coordinates in the dynamically most significant zone near the cross-section with
are close to
. In [
6], values in the range
–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.
4. Comparison of Computations with the Full and Simplified Models
The calculations with the full model were performed in a simple geometry with
. On the interval
, different initial distributions
,
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
,
,
,
,
kPa ([
1], Note 8).
was prescribed as triangular or parabolic forms with different values of
; for
, a uniform distribution was used (see
Figure 4).
Figure 5a presents the distribution of
as a function of
S in different stages of the motion.
Figure 5b plots the evolution of
,
H,
w, and
against
, where
is the
S-coordinate of the point where
. 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 , and H, w and become smooth and slowly varying functions of t with , for all the variants calculated.
Note that the asymptotic behavior of the solution emerges after the avalanche body has spread to the point where , so that the Coulomb friction law governs the entire avalanche body. Before this stage, significant changes in and 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
,
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
. 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
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.
5. Improvement of the Simplified Model
Including in the momentum equation the collisional term
and using all other simplifying considerations introduced in the previous section, we get the relations
Introducing the definitions
we can rewrite (
27) in the form
The solution of (
32) on the segments
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
, we get
The graph of
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
, we will be interested only in solutions of (
38) satisfying the conditions
On the
y-
S plane, the hyperbolas
([
1], Note 11) are the isoclines of (
38), i.e., the lines along which
On the hyperbolas
we have
, so all the solutions have extrema there (more precisley, these are maxima, for if
then
and at
). As a result, we can visualize the integral curves of Equation (
38) as shown in
Figure 7. Depending on the initial data, when
an integral curve tends to the hyperbola (
39),
, either monotonically or after some initial rise to a maximum value of
y. It can easily be shown that
on the curve
and that, at any given point of (
40) with coordinate
, a solution of (
38) touches some hyperbola
depending on the initial data. For
, this solution becomes “captured” between the curve (
40) and the hyperbola
. 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
. This value, thus, can be recommended for practical use over a wide range of parameters.
6. 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),
, and include the hydrodynamic resistance term,
, 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
.
Consider now the evaluation of the simplest model for several hypotheses about the bed shear stress .
6.1. 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
,
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
is known). So, the numerical evaluation of these expressions reduces to rather simple calculations. Note that if the slope is steep enough, i.e.,
, 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
and the second horizontal (
). Using (
47) and assuming
, we have ([
1], Note 14)
With
, where
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
, and substituting
, we have ([
1], Note 15)
In the second segment, the avalanche decelerates and completely stops at some distance
. So, for this distance, we have (substituting
,
in (
51))
Note that this runout distance does not depend at all on the avalanche size.
Using (
50) and (
52) and the definition of
, we obtain the full runout distance
as ([
1], Note 16)
If , 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
) ([
1], Note 17)
instead of (
49) and (
50).
On the horizontal path segment, we find
So, if we set
here, we get
instead of (
52).
Finally, we have for the full runout distance ([
1], Notes 16 and 18)
When
, Equation (
54)–(
57) reduce to (
49)–(
53), as they should.
6.2. 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
beyond which the friction law
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
, we have
and
For the second (horizontal) part, we have
and at
, where
, we obtain
The full distance
is given by the formula ([
1], Notes 16 and 19)
An important feature of the relations obtained in this case is their pronounced dependence on the avalanche size H. The dependence of the runout distance L on H is particularly interesting: The bigger (deeper) the avalanche is, the longer is its runout. In this respect, big avalanches governed by the limited-stress friction law with differ significantly from small avalanches governed by the classic Coulomb law.
This conclusion is, of course, valid not only for the simplest model considered here but also for the “exact” ones as it is a direct consequence of the physical difference between these two friction laws. The restriction of the bed shear stress to the limiting value leads to extremely high mobility of the snow flow on the mountain slope and the horizontal plane, with a dramatic rise of runout distance with increasing avalanche size.
It is also interesting to see how these results are modified when hydraulic resistance, represented by the term
, is included in Equation (
44). In this case, integration over the initial path segment gives
and in the second segment, one obtains ([
1], Note 20)
For the runout distance
, (
67) and
lead to
Taking the limit
in (
68), one retrieves (
63).
The strong dependence of
on avalanche size
H given by (
68) is obvious, and all the conclusions relating to the properties of big avalanches made above in the simplest case with
remain valid. However, the hydraulic resistance term significantly modifies the relationships (
63) and (
64) between runout distance and avalanche size
H: They are no longer linear since Equation (
68) contains logarithmic and exponential terms. As a result, the runout
rises much more slowly with
H than is predicted by the linear relations (
63) and (
64).
Consider a numerical example with the following parameter values ([
1], Note 8):
;
kPa;
;
;
;
,m/s
;
m;
m;
. The Coulomb friction law (
53) gives
,
m. In the case of the limited-stress friction law, the relation (
63) gives
,
m. In this case, the lower limit on the flow depth for the shear stress to be limited by
is
m. Thus, the value
m used in the calculations is only slightly above
. This is why
for this friction law is only about twice the value obtained with the Coulomb law. A value
kPa seems more adequate, in which case we get
m and
with
m. One can see that, in this case, the runout
is one order of magnitude longer than for the Coulomb law.
Consider now influence of hydraulic resistance on these numbers. Using the relation (
56) for the Coulomb law, we have ([
1], Note 21)
From (
68), we obtain (for
kPa)
Reducing
to 5 kPa results in
. 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
([
1], Note 22). The problem of how to adequately estimate the hydraulic part of the flow resistance remains open and needs additional research.
7. 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
,
,
kPa,
kPa,
,
,
, and
, 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
Figure 9 in terms of the distributions of
h and
u at the instance when the front reaches the point
m. Corresponding plots for the other simulations can be found in the
Supplementary Materials Document 2, Figures S2.10–S2.14.
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
, the end values of the range lead to avalanche dynamics of the second, “weak” type, while the intermediate value
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.
8. 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 (, , , , 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 , 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.