1. Introduction
It is commonly agreed that dry-slab avalanche release is triggered by shear failure of a thin weak layer underlying a thick cohesive slab (Reference SchweizerSchweizer, 1999). Several authors (e.g. Reference McClungMcClung, 1979; Reference Bader and SalmBader and Salm, 1990) have invoked fracture-mechanics concepts to understand the failure process. They consider the existence of an extended flaw (termed a shear band by McClung) which is, in fact, very similar to a Griffith crack in linear elastic fracture mechanics. As in fracture mechanics, the stability of the system depends on the extension and geometry of the critical flaw.
Recently, the concept of a critical crack has been criticized by several authors (e.g. Reference Arndt and NattermannArndt and Nattermann, 2001). It has been pointed out that, unless the nucleation problem is addressed and/or the size of pre-existing cracks can be determined experimentally, one simply replaces one unknown quantity (the failure stress) by another (the size of the critical crack). Reference Arndt and NattermannArndt and Nattermann (2001) also demonstrated that random variations in material strength may trigger athermal crack nucleation and may therefore play a crucial role in failure of bulk materials.
Large and apparently random variations in the shear strength of weak layers have been reported by Reference Conway and AbrahamsonConway and Abrahamson (1988). In their study they assess the implications for slow-slope stability by evaluating the probability of finding deficit zones where the strength of the weak layer falls below the acting (gravitational) stress. A problem of this approach is that the overlying slab may effect a stress redistribution such that a deficit zone is supported and “anchored” by surrounding stronger regions, and the situation may be further complicated by the presence and interactions of multiple deficit zones. In general, slope failure should therefore be considered as a collective process which involves both multiple local failure of the weak layer and concomitant stress redistribution to adjacent regions. To deal with such collective phenomena in slab avalanche release is the main rationale of the model which we formulate in the present paper.
2. Formulation of the Model
2.1 Stress redistribution in a thin slab
We consider the situation shown in Figure 1: a weak layer is running along the plane y = d′ in a snowpack of thickness d. For simplicity we assume the elastic constants (shear modulus G, Young’s modulus E) of the snow above and below the weak layer to be the same. Following the lines of Reference Palmer and RicePalmer and Rice (1973) and Reference McClungMcClung (1979), we consider a quasi-one-dimensional model: we assume the system to be homogeneous in the z direction, such that the displacement across the weak layer can be written as a function u(x) of the x coordinate only. The nature of the failure depends on the angle θ between the u and x directions: for θ = 0 we have mode II failure (shear failure along the weak layer propagating up-slope or downslope) and for θ = π/2 mode III failure (shear failure propagating sideways across slope). The extension of the slope is assumed much larger than any other characteristic length in the system such that we may neglect boundary effects. In generalization of previous models, we allow for random spatial variations of the shear strength of the weak layer as a function of the x coordinate.
Fig. 1. Snowpack geometry considered in the model.
Owing to the material inhomogeneity, J-integral methods as used by Reference Palmer and RicePalmer and Rice (1973) for analyzing shear band propagation do not work in the present case. Instead we have to explicitly consider the stress distribution which arises from a general distribution of slip u(x). This can be done using a dislocation representation of the internal stresses (Reference WeertmanWeertman, 1996; Reference Arndt and NattermannArndt and Nattermann, 2001). In particular, the internal shear stress acting at location x on the weak layer can be written as
Here τ(x) is the shear stress created at the point (x; d′) by a dislocation of unit strength located at (0; d′), and ▽xu(x) is the dislocation density associated with the inhomogeneous distribution u(x) of slip in the plane y = d′ of the weak layer.
To evaluate the internal stresses, we have to work out the stress field of a dislocation running in the z direction on the plane y = d′. The dislocation is contained in a layer with elastic moduli E, G sandwiched between media with E, G≈1for y < 0 (the bedrock) and E, G ≈ 0 for y > d (the open air). The stress field depends on the dislocation character, which in turn is governed by the angle θ between the u and x directions. For θ = π/2 (mode III failure) the crack dislocations have screw character, and we can use results given by Reference WangWang (1999) to write τ(x) as
We now make an important approximation: we assume that the snowpack is thin in the sense that variations in u occur on a characteristic length scale which significantly exceeds d and d′. We may then expand the slowly varying function 
in a Taylor series around 
Inserting into Equation (2) and retaining only the lowest-order non-vanishing term yields after some lengthy algebraic manipulations
where the stress redistribution factor
turns out to be proportional to the thickness of the slab above the weak layer but does not depend on the thickness of the underlying snowpack.
For other values of θ (mode II or mixed-mode fracture) the shear stress along the weak layer can be calculated using the stress fields of edge or mixed dislocations. The final result is again of the form of Equations (3) and (4). The only difference is that I multiplies with a scalar pre-factor of the order of unity which depends on the angle θ and on the value of Poisson’s ratio.
2.2. Equation for slab stability
The dependence of the shear strength of the weak layer on the shear displacement is characterized by a hardening– softening curve as shown schematically in Figure 2 (see, e.g., Reference McClungMcClung, 1979). The shear strength increases initially towards a peak value τm and then drops towards an asymptotic value τ∞. In general, the strength–displacement curve τs(u; x) may be position-dependent. In our simulations in section 4, we account for such position dependence in terms of random variations of the peak shear strength.
Fig. 2. Shear strength vs displacement across weak layer, and “equal-area” condition fulfilled by a shear band (schematically).
For the slab to be stable, the displacement field u(x) must fulfil the inequality
i.e. the locally acting (external and internal) stress must not exceed the local shear strength. Equation (5) is constitutive for our model. We consider rate-independent behavior, i.e. once condition (5) is violated for some x, the displacement u(x) at the unstable locations increases quasi-instantaneously until a new stable configuration is reached. If no such configuration exists, u increases indefinitely and the slab fails.
3. Failure of a Homogeneous Slab by Shear Band Propagation
We first consider the failure of a homogeneous slab (no strength variations) by formation and propagation of a shear band as investigated by Reference McClungMcClung (1979). A localized shear band is characterized by a displacement field u(x) which at x→–∞starts from a value u0 on the left stable branch of the 
curve, goes through a maximum u1 which without loss of generality we assume at x = 0, and then reverts to u0 for x→∞. The displacement field of a critical (marginally stable) shear band satisfies the equation
By analogy, this differential equation can be envisaged as describing the undamped motion of a particle of mass in I a potential 
It follows that the solution must satisfy the “energy conservation” criterion 
and hence u0 and u1 must fulfil the relation
Solutions can be found for any value of the external shear stress in the range τ1 ≤ τext ≤ τm. A geometrical illustration is given in Figure 2: areas I and II enclosed by the line τ = τext above and below the τs(u) curve must be equal.
Approximate analytical relations can be derived if the length of the shear band is large in comparison with that of the end region over which the strength of the weak layer drops to the residual value τ∞. Outside the end region the displacement profile is
Hence, 
The equal-area condition (7) can be written in the approximate form
where u* is the displacement where the descending branch of the stress–strain characteristics assumes the value 
Since a“long” shear band requires 
to be small, the integral on the righthand side is approximately equal to the area between the curve 
and the line 
We define the characteristic displacement ū
by
and, using Equations (8) and (9), we finally obtain the approximate relation
for a marginally stable shear band. A band of length L is stable up to the corresponding external stress. Slope failure then occurs if either the length of the band or the external stress is increased by an infinitesimal amount. The relation
(11) has been previously derived by Reference McClungMcClung (1979) from a quite different line of reasoning. The only difference is that here we consider mode III failure, whereas Reference McClungMcClung (1979) considers shear failure in mode II. This leads to a different combination of elastic moduli on the lefthand side of Equation (11) (McClung obtains 2E instead of 4G).
Formally, the above calculations are very similar to those carried out by Reference Zbib and AifantisZbib and Aifantis (1988) who used a second-order gradient approximation to the internal stress field to describe strain localization during shear band formation in a strain-softening material. Although we address a different problem, namely the “forward” extension of a shear band along the shear plane, the mathematical formalism is identical to that used by Zbib and Aifantis. The exercise of re-deriving relations obtained previously by other methods (see Reference Palmer and RicePalmer and Rice, 1973; Reference McClungMcClung, 1979) serves us mainly to demonstrate that the “gradient” approach to internal stresses can indeed be applied to the problem at hand. We now proceed to the case which is in the focus of the present work, i.e. failure of a slope where the weak layer exhibits random variations in strength.
4. Failure of a Slab with Random Strength Variations of the Weak Layer
4.1 Simulation method
When there are random spatial variations in the strength of the weak layer, no analytical solutions are available and one has to resort to computer simulations. To obtain reliable statistics, it is desirable to simulate huge ensembles of statistically equivalent systems. To improve computational efficiency, we use a lattice automaton technique where we allow the space coordinate to take only discrete values xi with constant spacing Δx. Accordingly, we replace uxx in Equation (5) by the discrete second-order gradient. Furthermorewe approximate the strength–displacement characteristics by a piecewise linear curve,
Randomness is introduced by considering the peak strength τm as a random function of space. We chose Δx to be equal to the characteristic correlation length ξ of the strength variations, such that we may envisage the peak strengths at the different “sites” xi as statistically independent random variables.
In the following, it is convenient to shift the zero of the stress axis to τ∞ and define non-dimensional space, displacement and stress variables via
In these variables, the stability condition for the “site” Xi reads
The peak strengths 
are assumed to be Weibull-distributed with the cumulative distribution function
This distribution is characterized by the parameters T and β. In the following, we use instead the mean value 〈TM〉 and relative variance σM/〈TM〉 of the peak strengths, which are related to the Weibull parameters by 
and 
where Г(x) denotes the gamma function. We have verified that the results remain qualitatively unchanged if other types of distribution (log-normal or box-shaped) are used.
A simulation is carried out as follows: Values of the average peak strength and peak strength variation are chosen, and random peak strengths are assigned to all sites according to the corresponding Weibull distribution. Then the external stress TEXT is increased from zero in small steps ΔT (this mimics slow loading as for instance by snowfall) until sites become unstable as the stress exceeds the local shear strength. The displacement at all unstable sites is increased by a small amount ΔU. Then, new internal stresses are computed for all sites and it is checked again where the sum of the external and internal stresses exceeds the local strength. The displacement at the now unstable sites is again increased, and so on. This is repeated until the system has reached a new stable configuration. Then the external stress is increased again and so on until the system has failed completely (Ui >1 for all sites). The corresponding critical stress is denoted by TC. The procedure is repeated for different values of ΔU and ΔT to ensure that the results do not depend on step size.
Typical values and ranges of the physical parameters of our system are compiled in Table 1. From these values, we find a typical range of the non-dimensional peak shear strength 1< TM <10. In the following we use, unless otherwise noted, the typical values TM = 5 and σM/〈TM〉 = 1 for the average peak shear strength and peak strength variation.
Table 1. Typical values and ranges of the physical parameters entering the model
4.2. Critical conditions for slope failure
Due to randomness of the local strength, the critical stress for failure of a slope of finite size is itself a random variable. The simulations indicate that the probability distribution of critical stresses is approximately Gaussian. The average failure stress decreases approximately in inverse proportion with the logarithm of system size XS.This is plausible since large systems (large slopes) have an enhanced probability to contain very weak configurations which may fail at low stress. However, the size dependence is not strong and, in view of the fact that the range of physical system sizes is restricted (XS = 1000 corresponds for ξ = 1m to a slope of 1km width), in the subsequent simulations we adopt a uniform system size of XS = 400.
The average failure stress increases linearly with the average peak strength of the weak layer (insert of Fig. 3). This figure also demonstrates that the average slope failure stress falls significantly below the average peak strength of the weak layer unless the strength distribution is very narrow. This effect is pronounced even at small values of σM/〈TM〉: already a relative scatter of the peak strengths of the order of 20% may decrease the slope failure stress by a factor of two. This indicates the importance of localized weak configurations for triggering failure. However, as we shall see in the next section, such configurations need not necessarily conform to the idea of point-like flaws or linearly extended shear bands.
Fig. 3. Dependence of slope failure stress on the average peak strength and peak-strength variation of the weak layer. Square data points: critical stresses for average peak strength 〈TM〉 = 5 and different values of the peak-strength variation σM/〈TM〉. Insert: dependence of failure stress on peak strength for two different values of the relative peak-strength variation. Each data point represents an average over an ensemble of 100 systems; the error bars represent the variance of the respective failure stress distributions.
4.3. Nature of the critical flaw
To investigate the nature of the critical flaw, we define a damage function D(x) which is unity at those sites where the strength of the weak layer has dropped to its residual value, and zero everywhere else. In our simulations we find that just before failure the distribution of damage is rather irregular (insert of Fig. 4). Instead of a single point-like flaw or an extended damaged area representing a shear band, one observes a “bar-code” pattern of irregular damage clusters which vaguely resembles a randomized Cantor set. To statistically characterize this pattern, we investigate the damage autocorrelation function C(X) = 〈D(X′)D′X + X′)〉. The average is evaluated over all values of X′ and, to obtain reliable statistics, we further average over an ensemble of many systems. Up to a characteristic damage correlation length XC ≈ 25 we observe a power-law decay C(X) α X-m, where m ≈ 0.46. This implies that, below XC, the damage distribution can be described as a random fractal set with fractal dimension D = (1 - m) = 0.54. Above XC there is a crossover to a Poissonian random pattern with fractal dimension 1. Investigation of systems of different sizes ranging between 100 ≤ XS ≤ 6400 indicates that the correlation length XC does not depend significantly on system size.
Fig. 4. Data points: two-point correlation function C(X) = 〈D(X′)D(X + X′)〉 of the damage pattern (average over final configurations of 1000 simulations); full line: power-law decay C(X) / X -0:46. Insert: distribution of damage at failure for one system.
4.4. Precursors to failure
In the presence of random strength variations, local failures may occur even before the system fails globally. In the simulations, we monitor such precursor activity by studying the displacement increments which occur after each stress increment. The insert in Figure 5 shows the increments in mean displacement which have been observed during a single simulation where the applied stress was increased in steps ΔT =0.02. One sees that the displacement increases in irregular bursts. Such bursts occur because the failure of one site may, through local stress redistribution, trigger the failure of others. (In other physical systems the term “avalanches” has become commonplace for this type of collective behavior. Here we prefer the term“bursts” for the obvious reason that we are talking about those precursor events which do not yet trigger an avalanche.) In real systems, such precursor events may be detectable through acoustic emission measurements and, hence, be used to gain experimental information about the internal stability of a slope.
Fig. 5. Average slope susceptibility (mean displacement increment per unit stress increment, average over 104 simulations) as a function of distance from critical stress. Insert: displacement increments in a single simulation with step size ΔT =0.02.
Figure 5 shows the susceptibility 〈ΔU〉/ΔT averaged over an ensemble of 104 systems. It is seen that the susceptibility increases as the critical stress is approached, but it does not exhibit any divergence towards failure. It is also instructive to calculate the fraction of failed sites. By analogy with equilibrium thermodynamics, this fraction may be considered an “order parameter” characterizing the failure process. This order parameter increases towards a finite value 51 as one approaches the critical stress from below, and then discontinuously jumps to unity. Hence, the behavior of the system (discontinuity of the order parameter, non-singular behavior of the susceptibility) is reminiscent of a first-order phase transition.
5. Discussion and Conclusions
In our model of slab failure the weak layer is modelled as a softening interface. Randomness is introduced by discretizing space and assigning random peak strength values to the different “sites”. Stress redistribution between strong and weak regions is expressed in terms of the (discrete) second-order gradient of the displacement, which makes the model very similar to fibre bundle models with local load sharing (Reference Harlow and PhoenixHarlowand Phoenix, 1991). As in such models, we observe a first-order-like transition at failure.
From our model, the following main conclusions are drawn:
Fluctuations in the peak strength of the weak layer have a strong knock-down effect on slope stability, as stability is governed by the probability of occurrence of localized weak configurations.
The critical flaw which leads to slope failure is neither a linearly extended shear band nor a point-like deficit. Rather it falls “in between” these two extremes: prior to slope failure, the distribution of failed “sites”can, up to a characteristic correlation length, be envisaged as a fractal set of dimension D ~ 0.5.
As the load increases, failure is preceded by small precursor events, which may be detectable by acoustic emission as proposed, for example, by Reference Sommerfeld and GublerSommerfeld and Gubler (1983). The precursor activity increases towards failure, but this increase is unspecific: since the susceptibility of the system does not diverge at the critical stress, it is difficult to pinpoint the failure stress from measurements of the precursor activity. (The situation may be compared with an attempt to determine the boiling point of a liquid from measurements of the density decrease prior to boiling.)
The present model relies on several strong idealizations which may affect the validity of our conclusions. (i) We describe stress redistribution in terms of a local load-sharing approximation. This is feasible as long as the correlation length ξ of strength variations is larger than the slab thickness. If (d - d′) > ξ, the stress redistribution effectively leads to an averaging of the fluctuations over a region of the order of (d - d′). Hence, for thick slabs our model tends to overestimate the impact of fluctuations. (ii) Snow is a rate-dependent material, and the assumption of rate-independent behavior should be relaxed in future investigations.
(iii) In its present form the model is quasi-one-dimensional. The assumption of homogeneity in one spatial direction is not very realistic, and, unfortunately, in statistical mechanics dimensionality matters. In the regime of large fluctuations, the behavior of the present one-dimensional model is governed by single very strong sites which may act as anchoring points for the whole slope. In two dimensions this is different since failure may bypass these sites, and ultimately the slope may avalanche, leaving the odd “super-strong” zone standing. A more detailed study of this issue is postponed to future investigations.
Acknowledgements
The author is grateful to M. Alava, J. Blackford, S. Zapperi and E. C. Aifantis for valuable discussions.
