Introduction
Reference McClungMcClung (1981), Reference SinghSingh (1980), and Salm (paper in preparation) have modelled shear-fracture propagation along thin, weak layers in snow; these theories assume the existence of initial shear-failure planes or deficit areas (Reference Conway and AbrahamsonConway, 1984) as a necessary condition for the start of fracture propagation. In this paper, we describe a model for the evolution of initial shear-failure planes at increasing load in the context of direct action avalanches. A simple model for ductile shear failure of low-density snow, combined with strength increase due to settling under increasing load and sintering depending on snow temperature, is used to predict fracture height of slab avalanches.
Interpretation of Field Measurements and Observations
Different types of field measurements and observations have been made at our institute to investigate natural slab releases. These include infrasonic acoustic emission measurements in potential release zones, stress and deformation measurements at the base of new snow layers, temperature profile measurements, and microwave radar stratigraphy determinations. The frequency of acoustic waves in snow originating from limited failure is of the order of O(υ) = c/(4πl), where c is shear-wave propagation velocity and l is characteristic linear dimension of the fracture. Measured frequencies correspond to typical extensions of fractures in the range 0.1-1 m occurring well before a slab is released (Reference Sommerfeld and GublerSommerfeld and Gubler, 1983). It can be concluded from these measurements that initial ductile failure planes have dimensions <l m2, and that the local failures develop independently from each other in potential slab-release zones. The observed increase in the rate of acoustic emissions before natural slab releases may by interpreted to mean that coalescence of different initial failure planes is a necessary condition, and causes the resulting failure plane to reach a size critical for the onset of fracture propagation. Shear-deformation measurements at the old-new snow interface show that mean shear strain-rates within a slab of 0.05 m thickness are of the order of 1 × 10−6 s−1 during heavy loading by precipitation. If this shear deformation is concentrated within an interface layer of a few millimeters, local strain-rates would increase to 1 × 10−5 s−1. Shear and normal stress measurements inhibiting strains at the old-new snow interface indicate that the range of influence of strain-inhibiting obstacles is small and decreases significantly with increasing snow temperature and therefore with decreasing snow viscosity. These findings support the idea that a local increase of stress or decrease of strength and viscosity may cause local failure. Because snow viscosity is proportional to some negative power of strain-rate (Reference NaritaNarita, 1983), high strain-rates may develop in the vicinity of strain-blocking distortions, such as trees and rocks, increasing the probability of local failures. Radar stratigraphy measurements (Reference GublerGubler, 1984) allowed for direct snow accumulation and settling measurements in release zones as well as estimates of the water equivalence of the released slab and determination of the sliding layer. Accumulation and settling rate have proved to be very important parameters in forecasting the occurrence of extreme avalanche cycles. Effective temperature conductivity in natural snow covers has been determined from periodical (1-24 h) temperature-profile measurements.
Failure Properties of Snow
Measurements of failure properties of snow (Reference NaritaNarita, 1983), as well as different kinds of field measurements and observations, support the idea that snow will initially fail in a ductile manner. Reference NaritaNarita (1983) has succeeded in demonstrating the start of snow failure as a function of increasing strain-rate. The factors affecting maximum stress, strength and strain-rates for low-density snow are summarized in Figure 1. According to Narita, with some extrapolations for lower snow densities and temperatures typical for snow layers involved in direct action avalanches, the mechanisms that determine tensile and approximate shear strength of snow as a function of strain, ε, and strain-rate, 
may be subdivided into four distinct categories:
Fig. 1. Basic dependence of maximum stress, strength (arbitrary units), and strain of snow on strain-rate (s−1)
Viscous behaviour without failure for 
.
Visco-plastic deformation with incomplete macroscopic failure at large strains for 
.
Ductile failure for 
.
Partly brittle failure for higher values of .
Corresponding strain to brittle failure is 0.1-0.3%, and increases with decreasing values for ductile failure from 1 to 10%. The ratio of maximum ductile strength to minimum brittle strengths is about 3.5:1; strengths are of the order of 5 × l02 to 5 × 103Pa. Reference McClungMcClung (1987) measured increasing strain softening with increasing shear strain (residual strength 50% for fine-grained snow and 25% for coarse-grained snow), basically confirming the findings of Narita. Strain-rates on inclined natural snow covers caused by the snow load itself hardly ever exceed 1 × 10−4s−1, even if concentrated in weak layers or by obstacles. Because this corresponds with shear-deformation measurements and model calculations, initial failure in connection with natural slab releases will be ductile. The existence of initial shear-failure planes is a necessary condition to start fracture propagation. Resulting stress concentrations at the tip of the crack, necessary to propagate the fracture, depend on the ratio of crack length to thickness of the low viscosity, the weak layer guiding the crack. The type of fracture propagation, slow ductile propagation at large strains before failure, and fast brittle-fracture propagation at small strains, depends on the amount of stress concentration available at the crack tip (Salm, paper in preparation).
Probable sequence of events for natural slab release
The following sequence of events is thought to be a plausible scenario for direct-action avalanches. It is based on measurements, observations, and experience in connection with artificial avalanche releases by skiers and explosives.
Development of layer interfaces of low cohesion by new snow on old snow, snow metamorphism, surface hoar, wind, melt water, etc.
Increasing stresses in weak layers caused by precipitation.
Viscous deformation partly concentrated in weak, low-viscosity layers.
Increase of strength because of settling and growth of bonds. Possible decrease of strength as the result of decrease in specific number of bonds caused by kinetic growth metamorphism.
Local concentration of stresses, strain-rates, and strains as the result of local in homogeneities in layering, or in the vicinity of obstacles, and in low viscosity layers.
Initial ductile failures locally increase stress and strain-rates.
Eventual coalescence of weakened zones to reach critical size for initiation of fracture propagation.
First propagation of ductile or brittle fractures initiated by ductile failures.
Stress repartitioning from weakened zones to more stable parts of slab may stop fracture propagation (infrasonic acoustic emission).
Coalescence of weakened zones by fracturing leads to unrestricted brittle-fracture propagation (primary fracture).
Model to Simulate Ductile Stability
From the facts, assumptions, and experiences mentioned above, it can be concluded that ductile failures of an identifiable minimum extension are a necessary condition for the onset of brittle unrestricted fracture propagation. The volume of the slab, in most cases some part of an interface layer, where a ductile failure of critical size occurs is called initial fracture volume. The model introduced below deals only with initial ductile failure as the first necessary step for direct action slab avalanche release. A microscopic model described earlier by Reference GublerGubler (1978a,b) is used to model initial failure in a weak layer as a function of increase of snow load by precipitation. The failure model is based mainly on observations and measurements on thin sections of snow by one of the authors and Reference NaritaNarita (1983).
Model assumptions
The model is based on the description of the strength of a bundle of fibres by Reference DanielsDaniels (1945). External stresses are distributed between all unbroken fibres and if the external stress is increased the weakest fibres break. The limiting strength of the bundle is reached if the external stress can no longer be redistributed between surviving fibres. Sequential breaking of links within a structure without storage of elastic energy at low strain-rate is a basic characteristic of ductile failure. Detailed assumptions of the model for snow are given inReference Gubler Gubler (1978a). To model low-density snow (ρ < 200 kg/m3) the fibres in the Daniels model are replaced by elementary units (EU) bridging the future failure surface. An EU consists of a chain of grains connecting grain clusters which are considered to be of low deformability and high strength, in contrast with the bonds connecting chain grains which are the weak and deformable parts of the EU. Stresses in ice bonds are significantly increased if the chains are bent. Observations by Reference NaritaNarita (1983) on thin sections made during the ductile failure process are in good agreement with this assumption for the development of ductile failure.
According to Reference GublerGubler (1978a), the distribution of strength. K(s), of the EU can be assumed to be logarithmic. For low-strength fine-grained snow, the density distribution of strength K(s) has its maximum at 1 × 10−5 − 10−4 N, and the specific number of EU per unit volume, Nv, is 1 × 10−11 − 1 × 10−10 m−3; the corresponding specific number of EU per unit area of the failure plane, Nf, is approximated as 
. For coarse-grained snow Nf is somewhat smaller for a similar density distribution of strength. These values correspond to a mean grain diameter of 30-100μm, a bond neck diameter of 10-30um, a bond ice strength of 1 × 106 Pa (Reference GublerGubler, 1982), and a two- to six-fold increase in maximum local stress in the bond ice due to chain-curvature radii equal to 0.5 mean grain diameters (Reference GublerGubler, 1976), compared with the mean stress in straight chains. The number of grains per EU is between 10 and 100, the mean chain length varies between 1 and 5 grain diameters.
Strength-determination procedures are given in Reference GublerGubler (1978b). Basically, strength within the future failure plane is determined by the surviving EU links as a function of stress and time. The strength distribution K(s), and also the total number of EU per unit area, Nf, varies as a function of time because of the increase of bond strength by sintering, variations of bond ice strength, and increasing density as a result of settling causing a proportional increase of Nv. The formation of new EU bridging the future failure plane as the result of shear strain, and the possibility of introducing a density distribution for EU loads instead of assuming equal loads (Reference GublerGubler, 1978a), have been ignored here, mainly because of the lack of confirmed values for the parameters involved. Knowing the strength distribution function K(s,t) for the EU as a function of time, the number, Ns, of surviving links per unit area with failure loads larger than S(t) (mean load per link increased to S(t)) can be determined as follows
The shear stress, ∑, transferred by the surviving links through the future failure plane is equal to
Ductile strength is reached if a redistribution of loads between the surviving EU links is no longer possible. This is equivalent to the condition ∑m ≤ shear stress in the failure plane, σ, with
To model the increase of strength of the EU as a function of both time and temperature, the following relationships have been determined (Reference KingeryKingery, 1969; Reference HobbsHobbs, 1974; Reference GublerGubler, 1982):
The temperature dependence for the growth rate of the bonds is approximated by
where T = snow temperature in K; 273K = 0°C;
R = 0.0038 kJ Mol−1 K.
The temperature dependence of ice strength is given by
Using Equations (4) and (5), EU strength, s(τ), at time τ, is found to be equal to
This model predicts an increase in strength proportional to a power of time, t, in correspondence with measurements obtained and the theory of snow sintering. For evaporation condensation as the limiting process for sintering Reference HobbsHobbs (1974) proposed n = 0.4, t0 is set = 1.
Snow density, ρ, determines the value of Nv. The relative increase of snow density in function of stress, σ, settling viscosity, η, and time, t, is
where G = weight per unit area of snow cover.
From measurements and theory of Claus (1978), Reference Feldt and BallardFeldt (1966), and Reference NaritaNarita (1983). η may be approximated by
where C = 7.06 × 10−11 Pa s, α = 0.02 m3 kg−1, E = 67.3 kJ Mol−1, R = 0.0083 kJ Mol−1 K−1, n = 0.737 for 
; for 
η is independent of 
.
Equations (6) and (8) indicate that snow temperature is one of the most important parameters in stability prediction; this IS in agreement with experimental observations. In this simplified model for ductile shear failure we determine stability as the ratio of shear strength to stress in a weak layer independent of viscosities of surrounding layers, although those viscosities are very important in modelling fracture propagation. In the model, calculation temperature in the critical layer is determined as a function of development of the snow surface temperature. The calculations are based on Equation (9), which describes the penetration of a sinusoidally varying surface temperature into a semi-infinite half-space, combined with an approximation of the observed situation of increasing snow depth and constant ground temperature.
where β = (ω/2m)½, ω = frequency of temperature wave, k = ω/νp, νp = (2mω>)½, m = coefficient of temperature diffusion = (3.8 × 10−7 m2 s−1).
These approximations have been checked using numerical solutions for a layered medium with time-dependent conductivities, densities, and layer thicknesses (new snow accumulation, settling). Details and an example are given in the next section.
For this preliminary modelling, location-independent shear stress along the future shear failure plane has been assumed. In reality, stress concentrations will occur at the boundaries of this initial failure plane, but so long as these initial failure planes are small, which acoustic emission measurements suggest they are, the above assumption is defensible. (Detailed calculations on this problem are currently in progress and will be published in a subsequent paper.) Stress concentrations will definitely affect ductile fracture propagation and the onset of brittle fracturing. In the model calculation, the governing equations were solved using a finite-difference method for time and depth. The results give density, temperature, viscosity, density distribution of EU strength, shear stress, ductile strength stability, and strain-rate in the critical layer as functions of time.
Numerical modelling of temperature in a layered snow-pack
We have restricted ourself to the treatment of layered snow-packs with translational symmetry parallel to the layers because this reduces the calculations to a one-dimensional problem. However, because of time-dependent layer thicknesses, the problem has moving boundaries.
For the example shown below the following boundary conditions and parameters have been chosen:
Boundary conditions: 0°C at ground-snow interface; sinusoidally varying surface temperature with a period of 1 d, amplitude 5°C; mean temperature −8°C.
Physical parameters: specific heat = c■ = 2090 J K −1 kg−1, heat conductivity = λ = (0.0007ρ + 0.026) J K−1 m−1 s−1, (10) ρi = density of layer i resulting from application of Equation (7) to each layer i, where Gi represents the weight of the layers above.
The one-dimensional heat equation for these parameters and boundary conditions was solved iteratively for time. The numerical methods used were explicit Euler and implicit Crank-Nicholson schemes generalized for moving boundaries. In each time-step, the new layer thicknesses were calculated by applying Equation (7). A more detailed description will be given in a subsequent paper but results for a typical problem are shown in Figures 2 and 3. Figure 2 shows the thicknesses of the layers due to settlement and accumulation and Figure 3 represents the temperature in the different layers as a function of time. The density of the old layer at time t = 0 was 300 kg m−3. We have assumed a constant snow fall rate of 30 mm during the 2 d precipitation period and in this period the density of the new snow was 70 kg m−3. Comparision of numerical solutions with the approximations made using Equation (6) shows good agreement, with only small deviations for very thin new snow layers.
Results of modelling
Model calculations have been made for precipitation periods of 48 h. Loading was chosen to be either constant at 30 Pa h−1 and 50 Pa h−1, or to vary daily between 0 Pa h−1 and 60-100 Pa h−1. The superimposed surface temperatures also had a period of 24 h with amplitudes of less than 5°C and means varying between −20°C and the melting point. The slope inclinations were fixed at 30° and 40°.
The example in Figure 4 for varying precipitation rate shows that critical stability in shear is not always well defined in respect of time. Because we model only initial failure and not slab failure, stability with respect to initial shear failure reaching a value of one does not necessarily mean that the slab will be released immediately. In most cases coalescence of many initial failure planes will be necessary to initiate unrestricted fracture propagation. Therefore, in this example, natural slab release may take place any time between 16 and 35 h after the beginning of the precipitation period. The main results are summarized in Figure 5. Approximate fracture height of the slab respectively critical shear stress in the failure plane for the time where stability crosses one are plotted as a function of slope angle and vertical stress rate (precipitation rate) or mean shear-stress rate. In addition, curves for different characteristic temperatures of the failure plane are given Because the temperature of the failure plane increases in all examples due to the zero degree ground-snow interface
Fig. 2. Thicknesses (m) of layers as a function of time (h), showing the settling of the new and old snow layers. The layer temperature partly responsible for the actual settling rate is given in Figure 3.
Fig. 3. Temperatures (°C) of the different layers as a function of time. Layer symbols correspond to those in Figure 2. Old snow layer at top of plot (0°C = temperature of ground-snow interface), snow surface with sinusoidally varying temperature at bottom of plot.
temperature after some variations during the starting phase at low new snow heights, the lowest temperature encountered by the future failure plane during the precipitation period has been chosen as a characteristic temperature. The temperatures at the time of slab fracture vary in a similar fashion but on a reduced scale, and from this we deduce that the characteristic temperature must be interpreted as an indicator of the varying effective temperature of the critical layer (the future failure plane), influencing local density increases and the strength distributions of the EU.
From the simulations performed so far, the following rules with respect to fracture height can be deduced:
1. Very small fracture heights result from very low initial layer strength with strong dependence on surface temperature variations on very steep slopes.
2. Small fracture heights result from low initial and sustained temperature of critical layer <-10°C, high snowfall intensity, high slope angle.
3. Small to medium fracture heights result from a temperature of critical layer very close to the melting point for all times.
4. Large fracture heights result from medium initial and sustained temperature of critical layer −5° to −8°C combined with moderate periodic snowfall intensity and lower slope angle.
Fig. 4. Simulation of initial stability (E), settling strain-rate (F), density (C), and temperature in critical layer as function of snow-accumulation rate (A) and surface temperature (B). Snow water equivalence of 30 kg m−2 is assumed below critical layer.
5. Very large fracture heights or often no avalanche releases, because stresses do not reach strength during a limited precipitation period, result from a temperature of critical layer > - 5°C but well below melting temperature, moderate snowfall intensity, and low slope inclination.
6. Fracture height increases at low layer temperatures with decreasing variance of the probability distribution, K(s), of strength of the EU, corresponding to an increase of homogeneity within the critical layer.
Fig. 5. Summary of results of stability simulation. Shear stress at initial shear failure plotted as function of slope angle and vertical stress rate (precipitation rate) for different snow temperatures. Temperatures (K) indicated are characteristic values for shear layer. L is linear stress increase in time; p is periodic stress increase in time. Equal total accumulation after 24 h for corresponding periodic and constant accumulation rates. Approximate shear-stress rates and fracture heights also given.
Resulting fracture heights and fracture timing are at least qualitatively in agreement with experience. Investigations by Föhn (personal communication) clearly show that the probability of occurrence of extreme avalanches with large fracture heights decreases with increasing settling rate of the snow cover at higher-than-mean settling rates for a given climatic region. Release zones of low angle of inclination are a necessary condition for occurrences of extreme avalanches. Increasing settling rates of snow (increase with rising snow temperature) at low to intermediate settling rates may either lead to increased fracture heights and slab sizes or may inhibit slab failures even at high load rates during a precipitation period of limited time. These aspects are clearly reproduced by the model which shows that fracture height increases with decreasing slope angle of the release zone with a preference of high precipitation rates on low inclinations rather than low precipitation rates on steeper slopes. Shallow slabs of originally cold snow on steep slopes are often triggered by short-wave radiation heating. This phenomenon is commonly experienced by ski patrollers and has been clearly shown by acoustic emission measurements. It also follows from model calculations, and is due mainly to the decrease in ice-bond strength with increasing temperature. In this context it has to be mentioned that morphological changes in near-surface layers caused by short- and long-wave balance at the snow surface may amplify the effect.
Final Remarks
The mechanical properties of newly fallen snow vary significantly during the first hours and days after fall as functions of temperature and stress. Modelling of stability with respect to initial shear failure using a semi-macroscopic model for failure combined with numerical determinations of snow settling and temperature can help to improve our knowledge of natural stability development. Modelling and measuring of mechanical parameters and temperatures of the different layers of the snow cover needs to be refined to make it possible to arrive at a substantial improvement in stability evaluation. As a next step, snow metamorphism, as a further cause of changes in mechanical properties, should be included in the model.
