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 Gubler (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 Narita (1983).

### Model assumptions

The model is based on the description of the strength of a bundle of fibres by Daniels (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 in Gubler (1978a). To model low-density snow (ρ < 200 kg/m^{3}) 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 Narita (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 Gubler (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, N_{v}, is 1 × 10^{−11} − 1 × 10^{−10} m^{−3}; the corresponding specific number of EU per unit area of the failure plane, N_{f}, is approximated as . For coarse-grained snow N_{f} 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 × 10^{6} Pa (Gubler, 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 (Gubler, 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 Gubler (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, N_{f}, 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 N_{v}. 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 (Gubler, 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, N_{s}, 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 (Kingery, 1969;
Hobbs, 1974;
Gubler, 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 Hobbs (1974) proposed n = 0.4, t_{0} is set = 1.

Snow density, ρ, determines the value of N_{v}. 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), Feldt (1966), and Narita (1983). η may be approximated by

where C = 7.06 × 10^{−11} Pa s, α = 0.02 m^{3} 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} m^{2} 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 G_{i} 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.