Influence of weak-layer heterogeneity on snow slab avalanche release: application to the evaluation of avalanche release depths

Abstract

The evaluation of avalanche release depths constitutes a great challenge for risk assessment in mountainous areas. This study focuses on slab avalanches, which generally result from the rupture of a weak layer underlying a cohesive slab. We use the finite-element code Cast3M to build a mechanical model of the slab/weak-layer system, taking into account two key ingredients for the description of avalanche release: weak-layer heterogeneity and stress redistribution via slab elasticity. The system is loaded by increasing the slope angle until rupture. We first examine the cases of one single and two interacting weak spots in the weak layer, in order to validate the model. We then study the case of heterogeneous weak layers represented through Gaussian distributions of the cohesion with a spherical spatial covariance. Several simulations for different realizations of weak-layer heterogeneity are carried out and the influence of slab depth and heterogeneity correlation length on avalanche release angle distributions is analyzed. We show, in particular, a heterogeneity smoothing effect caused by slab elasticity. Finally, this mechanically based probabilistic model is coupled with extreme snowfall distributions. A sensitivity analysis of the predicted distributions enables us to determine the values of mechanical parameters that provide the best fit to field data.

1. Introduction

The evaluation of avalanche release depths is of great importance for all applications related to hazard mapping or zoning. In particular, avalanche release depths represent a crucial input ingredient for dynamical runout models (Barbolini and others, 2000; Naaim and others, 2003) and are required for implementing statistical–dynamical simulations (Meunier and Ancey, 2004; Eckert and others, 2010). It has been shown in various studies (Hutter, 1996; Bartelt and others, 1999; Jamieson and others, 2008) that the outputs of these models in terms of runout distances and impact forces are strongly dependent on the release mass, as well as on other terms, such as friction, deposition and erosion.

It is commonly accepted that dry-snow slab avalanches are initiated by a shear failure in a weak snow layer (or at a weak interface) followed by tensile crown failure of the overlying slab (McClung, 1979; Schweizer and others, 2003). The shear failure is caused by a local loss of cohesion inside the weak layer, that may be due to (1) a localized surface loading, such as skiers or explosives (artificial release), (2) uniform loading due to a new snowfall (natural release) or (3) changes in the snowpack properties due to weather changes (natural release).

Most of the existing avalanche release models are based on the assumption first made by Palmer and Rice (1973) for an overconsolidated clay and then taken up for snow by McClung (1979), that a weak spot or shear band, i.e. a zone of zero shear strength, pre-exists inside the weak layer. Fracture mechanics is then applied to study the conditions for shear-band propagation. Based on an energy budget at the tip of the band, rapid propagation occurs when a critical length is reached after a phase of slow strain-softening. Schweizer (1999) produced a complete review of these fracture mechanics models and gave critical length values as a function of snow characteristics. More recently, Chiaia and others (2008) showed that a simple stress failure criterion from equilibrium equations could also be sufficient to predict shear-band propagation with good accuracy. Nevertheless, as already noted by Schweizer (1999), the very concept of a weak spot is questionable since, even if we can imagine how local weak zones could appear (e.g. around rocks, where the snow depth is reduced and thus the temperature gradient is increased), it is probably too simplistic, in general, to represent the complex heterogeneity of weak-layer mechanical properties.

Several studies (Conway and Abrahamson, 1988; Jamieson and Johnston, 2001; Birkeland and others, 2004; Kronholm, 2004; Schweizer and others, 2008; Bellaire and Schweizer, 2011) have shown that snow mechanical properties present considerable spatial variability. From field data, this variability is generally described as following Gaussian distributions with spatial correlations. Hence, the concept of weak spot may be replaced by mechanical models that take into account spatial stochastic processes to represent the heterogeneity. Such models would explain not only failure initiation in weak zones, but also fracture arrest in stronger zones (Schweizer, 1999). Recently, several studies have attempted to include this heterogeneity in mechanical models. These studies can be classified as follows, according to the numerical method used:

1. 1. Fiber bundle models (FBM) are simple statistical fracture models that are well adapted for representing spatially heterogeneous systems, including possible time-dependent effects, such as sintering. Using this framework, Reiweger and others (2009) model the weak snow layer as a discrete set of parallel brittle–elastic fibers. Spatial variability is accounted for by assigning each fiber an initial strength taken within a Weibull distribution. Despite the simplicity of the model, these authors are able to quantitatively reproduce laboratory shearing experiments on homogeneous snow samples. However, complex stress redistribution effects from elasticity of the slab cannot be taken into account in such models, which are therefore unable to reproduce full-scale avalanche release.

2. 2. Cellular-automata models (CAM) consist of a regular grid of cells, characterized by a state that can change over time and as a function of the state of neighboring cells. Fyffe and Zaiser (2004, 2007) and Kronholm and Birkeland (2005) applied such an approach to slab avalanche release with a heterogeneous shear strength of the weak layer represented by Weibull or Gaussian distributions. Using neighboring elements, their models account for stress redistribution between weak and strong regions. They also included mode II (shear) rupture of the weak layer with a strain-softening law, and mode I (tensile) rupture of the slab, as did Failletaz and others (2004). These studies demonstrated the influence of spatial variability characteristics (variance, nugget effect, etc.) on release depth distributions. Fyffe and Zaiser (2007) were also able to reproduce (under certain conditions) release depth distributions following power laws, as in field studies (Rosenthal and Elder, 2002; McClung, 2003; Failletaz and others, 2004). However, stress elastic redistribution effects are oversimplified in these models, and their applicability is limited to the case of a shear strength correlation length lower than the slab depth.

3. 3. Finite-element models (FEM) rely on the resolution of the complete mechanical equations of the problem. One of the main advantages of FEM compared to CAM is that they are able to capture large-scale stress redistribution effects due to elasticity. FEM has been successfully applied to modeling the mechanical response of sandwich specimens, including weak-snow layers with homogeneous properties (e.g. Bader and Salm, 1990; Stoffel, 2005; Mahajan and Joshi, 2008; Mahajan and others, 2010). However, only a few studies coupled FEM with a stochastic representation of the spatial variability. Recently, such a model was introduced by Griffiths and Fenton (2004) to study soil stability. To our knowledge, the same type of approach has not yet been undertaken for snow.

In addition, some recent studies (Johnson and others, 2004; Van Herwijnen and Heierli, 2009) relying on field data show that the shear failure of the weak layer tends to be systematically accompanied by a normal collapse. These authors argued that slope-normal and slope-parallel displacements occur simultaneously during release. Anticrack analytical models have been developed and have proved capable of reproducing these data (Heierli and Zaiser, 2008; Heierli and others, 2008). However, the issue of the influence of normal collapse on avalanche release is still a matter of debate. Some authors (Jamieson and Schweizer, 2000; Johnson and others, 2004) suggest that the simultaneous occurrence of weak-layer collapse and shear failure may facilitate fracture propagation due to bending effects. McClung (2011), however, showed that a model that does not account for slope-normal failure can reasonably reproduce critical length measurements obtained in field saw-cut tests. Hence, he argued that the slope-parallel propagation is very little influenced by the interaction between slope-normal displacement and stress. In the present study, the effect of normal collapse is not considered.

Here we use a finite-element method to build a mechanically based probabilistic model of the slab/weaklayer system, taking into account the two key ingredients: the redistribution effects by elasticity of the slab and the heterogeneity of weak-layer mechanical properties. The objective of this study is the evaluation of avalanche release depth statistical distributions, distributions that could be later coupled to propagation models for hazard mapping and zoning. In Sections 2 and 3, we present the system and validate the model using the well-known weak-spot case. This validation is also used to highlight a characteristic length of the system associated with the elastic redistribution of stresses. In Section 4, results obtained with a realistic spatial heterogeneity based on field data are presented. We analyze, in particular, how slab stability depends on slab depth and on the spatial correlation length of weak-layer properties. Finally, in Section 5, the obtained release depth distributions are combined with snowfall extreme value distributions and compared with field data of avalanche crown depths.

2. Mechanically Based Probabilistic Model

2.1. Objectives

As mentioned above, the objective of the model developed in this study is the evaluation of probability distributions of avalanche release depths, in particular in the context of absent or scarce data. More precisely, the aim is to produce a tool capable of predicting release depth distributions that are meaningful over relatively long timescales (typically several decades), and that could be used as input for hazard mapping procedures, such as statistical–dynamical approaches (Keylock and others, 1999; Eckert and others, 2008, 2010). Hence, the objective is not to develop a complete mechanical model of slab avalanche release accounting for all the complex processes at play. Both the geometry and the mechanical behavior of the system are greatly simplified, to reduce the number of poorly known parameters, while keeping the ingredients essential to describe the mechanics of slab release. Moreover, to be compared with field data, the predictions of this model then need to be coupled with a description of snowfall distributions, as shown by Gaume and others (2012). In the present paper, we focus only on the formulation and validation (numerical consistency) of the mechanical part of the model. In Section 5, however, the sensitivity of the predicted distributions (after coupling with snowfalls) to several mechanical parameters is presented.

2.2. Formulation of the model

The model is based on the finite-element code Cast3M (Verpeaux and others, 1988). The resolution procedure used (‘Pasapas’; Charras and Di Paola, 2011) enables the consideration of nonlinear models with an implicit integration scheme based on the weighted- residuals method. The momentum conservation equations, including inertial terms, are solved under the small deformations hypothesis:

(1)

(2)

with M the mass matrix, u the node displacement vector, σ the stress tensor, ε the deformation tensor, F the integrated force vector at nodes and D the damping matrix. In our study, the matrix D is taken as zero.

The system considered is a two-dimensional (plane stress conditions) uniform slope inclined at an angle θ, of length L = 50 m (Fig. 1). The x-axis is in the slope-parallel direction and the z-axis is orthogonal to the slope. The system consists of a slab of thickness h overlying a weak layer, modeled as an interface of zero thickness. The mesh is composed of 100 elements in the slope-parallel direction, x, and six elements in direction z. We used quadrilateral elements for the slab (QUA4: four nodes with two degrees of freedom per node) and joint elements for the weak-layer interface (JOI2: four nodes with two degrees of freedom per node). We checked that the mesh resolution is fine enough not to influence the results presented (Section 3.1).

Fig. 1. Geometry of the system: a weak-layer interface under a cohesive slab of depth h.

The boundary conditions applied to the slab are as follows: at the upper end of the slope (BC1) a shear stress, σ_xz = −ρg(z + h)sin θ, is applied in order to avoid bending of the slab linked to finite-size effects; at the lower end (BC2), a nil displacement in the slope-parallel direction, x, is imposed. The upper surface of the slab is free and the base is subjected to an interface law, i.e. a law relating shear stress to tangential displacement, which represents the weak layer.

2.3. Constitutive relationships

Snow is a very complex material whose mechanical behavior is still not fully understood. In the present model, only the ingredients necessary to produce realistic instability of the system, namely strain-softening of the weak layer and elasticity of the slab, are taken into account. Table 1 summarizes the values of the different mechanical parameters used in this study.

Table 1. Mechanical parameters used in this study and typical ranges of variation

Weak layer

Various studies (McClung, 1979; Föhn and others, 1998; McClung and Schweizer, 2006) have shown that weak snow layers behave as strain-softening (or quasi-brittle) materials. The softening is caused by the break of ice bridges at the microscopic scale. In existing mechanical models (McClung, 1979; Bažant and others, 2003; Fyffe and Zaiser, 2004, 2007; Mahajan and Joshi, 2008; Mahajan and others, 2010; Gaume and others, 2011), weak layers are generally characterized by a rupture displacement, u _p, and a critical softening displacement, δ. The pre-peak behavior is considered to be elastic, but stiffness values are very difficult to obtain, since these layers are generally very thin and unstable.

In the present study, the weak layer is modeled as a displacement-softening interface with a simple, linear piecewise relationship between shear stress, τ, and tangential displacement, u (Fig. 2). The value of the shear stress peak, τ _p, is governed by the Mohr–Coulomb criterion: τ _p = c + σ _n tan Φ, with c the weak-layer cohesion, σ _n the normal stress and Φ the friction angle. The friction angle is chosen as constant, Φ = 30° (De Montmollin, 1978; Van Herwijnen and Heierli, 2009), and the cohesion, c, is spatially heterogeneous, as described below. The tangential stiffness of the weak layer during the pre-peak phase is given by k _s = τ _p/u _p (Fig. 2). After the peak, the shear stress decreases (shear softening), until reaching a residual value, τ _r = σ _n tan Φ. This residual value corresponds to the situation where ice bridges are completely broken and only the friction between the slab and the underlying layer remains. Following McClung (1977), both the characteristic peak and softening displacement, u _p and δ, are taken equal to 2 mm (Bažant and others, 2003; Fyffe and Zaiser, 2004, 2007). Hence, the displacement to reach the residual stress u _r = u _p + δ = 4 mm. Note that, according to recent studies (McClung, 2009, 2011), the softening displacement, δ, under high loading rates tends to be smaller than the value assumed in this study, 0.1–0.3 mm. However, as shown below, the precise value of this parameter does not influence the outcomes of our model.

Fig. 2. Weak-layer constitutive law.

Slab

At high loading rates, such as those characteristic of slab avalanche release, laboratory experiments (Mellor, 1975; McClung, 1977; Narita, 1980; Navarre and others, 1992; Schweizer, 1998) and field measurements (Roch, 1966; De Montmollin, 1978; Jamieson and Johnston, 1990; Föhn and others, 1998; McClung and Schweizer, 2006) have shown that cohesive snow behaves as a brittle–elastic material. In the present study (as explained in Section 2.6) the possible tensile rupture of the slab is not directly modeled. Therefore, the mechanical behavior of the slab is modeled by an isotropic elastic law:

(3)

where E is the Young’s modulus and v the Poisson’s ratio. The following values have been used: E = 1 MPa, v = 0.2 and ρ = 250 kg m⁻³ for the slab density. In addition, in order to stabilize the computations, it has been necessary to slightly enrich the above constitutive law by the addition of a viscous term. All viscosity values within the range 10⁴– 10⁹ Pa s have been found to yield satisfactory results. We retain the value η = 10⁸ Pa s, which is in agreement with real snow viscosity measurements (Mellor, 1975; Camponovo and Schweizer, 2001; Schweizer, 1999).

2.4. Spatial heterogeneity

The spatial heterogeneity of the weak layer is modeled through a stochastic distribution of cohesion, c. Following Jamieson and Johnston (2001) and Kronholm and Birkeland (2005), we consider a Gaussian distribution of average 〈c〉 and standard deviation σ_c , with a spherical covariance function C(d):

(4)

where d is the distance between two points, ∊ is the spatial correlation length and the function is 1 if , 0 otherwise. The correlation length, ∊, represents the distance over which the cohesion values are significantly correlated. Note that in the present model, no nugget effect is considered (i.e. when d → 0). The effect of the nugget on avalanche size has been investigated by Kronholm and Birkeland (2005) using a CAM.

Figure 3 shows examples of cohesion field realizations with different values of the correlation length, ∊. These fields were generated using the turning bands method (Chilès and Delfiner, 1999) and we checked that, with the used mesh size, the obtained empirical covariance functions are in good agreement with the predictions of Eqn (4) for all values of ∊ investigated. Existing studies are not conclusive on the typical correlation length scale (Jamieson and Johnston, 2001; Schweizer and others, 2008; Bellaire and Schweizer, 2011) relevant for weak-snow layers. Schweizer and others (2008) recommended spacing snow pits at least 10 m apart in order to have independent results, thus suggesting that the correlation length, ∊, varies approximately within the 0.5–10 m range. In our study, ∊ was varied between 0.5 and 40 m (the lower limit, ∊ = 0.5 m, being imposed by our mesh size), but, due to finite-size effects, only the results with ∊ ≤ 10 m can be cross-compared (for ∊ > 10 m, the average cohesion, 〈c〉, begins to evolve with ∊). Lastly, the average cohesion, 〈c〉, was taken equal to 1 kPa and cohesion standard deviation, σ_c , to 0.3 kPa.

Fig. 3. Left: examples of the heterogeneity of the cohesion for different values of the correlation length, ∊. Right: comparison between the empirical normalized covariance function of the cohesion fields (computed from 100 independent realizations: circles) and the theoretical expression given by Eqn (4) (lines).

2.5. Loading

Gravity is the only applied external force and the system is loaded by progressively increasing the slope angle, θ, at constant slab depth until rupture. As will be shown, this loading procedure is equivalent to a progressive increase of the slab depth, h, at constant slope angle. The loading curve is shown in Figure 4. After an initial stage during which gravity is increased from 0 to 9.81 m s⁻², the loading is applied in two phases (Fig. 4).

Fig. 4. Applied loading curve. The time-step (t.s.) is 0.1 s. The blue curve represents a typical velocity evolution showing a marked precursor event.

First, the slope angle is increased from 0 to θ ₁ = Φ = 30° with a fast loading speed of 0.4° per time-step, since no failure can occur during this stage (τ < τ _p). The loading speed is then reduced to 0.04° per time-step until rupture occurs. We checked that, with this two-phase procedure, the chosen loading speed values do not influence the results presented. This simple loading is sufficient to study avalanche releases triggered by a progressive accumulation of snow. We emphasize that our model is not meant to account for the slow processes (snow metamorphism, viscous stress redistributions) active during the formation of the snowpack.

2.6. Rupture mechanism

Avalanche releases observed in our simulations are always induced by a local shear rupture inside the weak layer, which then propagates to both extremities of the slope. With reference to real avalanche releases, this corresponds to the case where the weak-layer heterogeneity is not sufficient to induce a tensile rupture within the slab, and slab rupture is thus triggered by morphological features, such as ridges, slope breaks, rocks, trees or other defects within the slab. We also performed simulations using a brittle–elastic constitutive law for the slab. These simulations showed that for a realistic tensile strength value (σ _t = 2 kPa; Mellor, 1975; Jamieson and Johnston, 1990; Sigrist, 2006), effectively no tensile failure occurred within the slab (the slab remained elastic) and the rupture mechanism was identical to that reported in this study. Note, however, that other sets of parameters may lead to tensile ruptures within the slab, which constitutes the subject of a further study (Gaume and others, 2011).

2.7. Avalanche release criterion

We consider that an avalanche occurs in our simulations when the following kinematic release criterion is met:

(5)

Hence, an avalanche is detected at time-step i when the velocity, , of any node of the weak layer is higher than N times the average velocity recorded over m previous time-steps, . Values of N = 10 and m = 10 were chosen in order to ensure that the criterion is not sensitive to small velocity variations triggered by potential precursor events (Fig. 4), which could lead to incorrect release angle values.

3. Mechanical Validation

In this section, the finite-element model presented above is validated against the classical case of release induced by a single weak spot. For this simple case, analytical solutions can be derived, following the approach of Chiaia and others (2008). The interaction between two weak spots is also considered, in order to illustrate the influence of an important characteristic length of the system that emerges from the analysis.

3.1. A single weak spot

Analytical solution

Here we follow the approach of Chiaia and others (2008), but considering a nonzero residual stress inside the weak spot due to friction (Fig. 2). Let us consider a weak spot of nil cohesion (c = 0), half-length a, inside a weak layer of homogeneous cohesion, c = 1 kPa, underlying a cohesive slab of depth h (Fig. 5). The equilibrium equation in the slope-parallel direction integrated over the slab depth is

(6)

with τ _g = ρgh sin θ the body weight shear stress, σ_uxx the normal stress in the slope-parallel direction and τ the shear stress in the weak layer. The shear stress, τ, is related to the tangential displacement, u, according to the interface constitutive law (Fig. 2). Two cases have to be distinguished: θ < Φ for which the shear stress, τ, depends on the tangential displacement, u, both inside and outside the weak spot; and θ > Φ for which the shear stress, τ, depends on the tangential displacement, u, only outside the weak spot.

Fig. 5. Geometry of the system with one weak spot of half-length a with nil cohesion.

Case 1: θ < Φ

For θ < Φ, the shear stress is τ(x) = τ _p u/ u _p outside the weak spot (|x| > a) and τ(x) = τ _ws = τ _r u=u _p inside the weak spot (|x| < a). Equation (6) and the linear elastic behavior of the slab lead to

(7)

with E′ = E/(1 − v ²) (plane stress hypothesis) and

(8)

(9)

Considering, in addition, the continuity of displacement and velocity at the interface between the weak layer and the weak spot, and the fact that the slope-parallel normal stress, σ_xx , vanishes far away from the weak spot, displacement and stress profiles can be determined by integrating Eqn (7).

Outside the weak spot (|x| > a):

(10)

Inside the weak spot (|x| < a):

(11)

with expressions for the constants r, r′, α and β given in Appendix A.

Case 2: θ > Φ

If θ > Φ, the shear stress inside the weak spot (|x| < a) meets the frictional criterion and thus no longer depends on displacement, u: T(x) = T _r = _17N tan 4). Equation (6) then becomes

(12)

for |x| ≤ a and Eqn (7) remains valid for |x| > a. Similarly to the previous case, the displacement and stress profiles can be determined.

Outside the weak spot (|x| > a):

(13)

Inside the weak spot (|x| < a):

(14)

with the expression for the constant r ₂ given in Appendix A.

We note that both the shear stress and the displacement present decreasing exponential profiles outside the weak spot (Eqns (10) and (13)). The characteristic length associated with these exponential decreases is the parameter Λ, which depends both on the slab and weak-layer characteristics (Eqn (9)). Far from the weak spot (|x| − a ≫ Λ), the shear stress tends to its body weight value, τ _g, and the displacement tends to u = u _p τ _g /τ_p (elastic behavior).

The shear band becomes unstable when the maximum stress, τ _max, at |x| = a reaches τ _p = c + σ_N tan Φ). Using Eqn (13), the theoretical critical stress τ _{g, s} for weak spot propagation can thus be expressed as

(15)

From this expression, the critical release angle, θ _r,s, can then be derived using

(16)

Comparison with simulations

As shown in Figure 6a and b, the overall agreement between theoretical predictions and FEM numerical results is very satisfactory, both for θ < Φ and for θ > Φ. In particular, the stress concentration at the weak-spot tip and the exponential decrease of stress outside the weak spot are very well reproduced. Similarly, the displacement profiles, which present a maximum at the center of the weak spot, are also well captured. Note, however, that the numerical model indicates the existence of slight variations with x of the normal stress, σ _n, which are not accounted for in the theoretical analysis.

Fig. 6. Displacement (u) and apparent friction coefficient (μ = τ/σ _n) profiles for a slab depth h = 1 m and a weak-spot half-length a = 4.5 m for (a) the case 0 < Φ and (b) the case φ < θ < θ _r, where θ _r is the release angle. Symbols: finite-element results (u and μ). Solid curves: analytical model (u _th and μ _th).

Figure 7 shows a comparison between the release angles obtained by the FEM calculations and those predicted by the stress rupture criterion (Eqns (15) and (16)). Here also, the agreement between the theory and numerical results is excellent for all values of weak-spot half-lengths. This agreement also holds for all tested values of slab depth, h. Globally, the results shown in Figures 6 and 7 constitute a validation of the various mechanical ingredients taken into account in our finite-element model. In particular, these results prove that the used mesh size is fine enough to account for cohesion heterogeneities with typical length scales (in this case, the weak-spot half-length) as small as 0.5 m.

Fig. 7. Evolution of the release angle, θ _r, as a function of the half-length, a, of the weak spot. Black squares: finite-element results Red dashed curve: theoretical stress rupture threshold.3.2. Two weak spots

3.2. Two weak spots

In order to illustrate the influence of the characteristic length, Λ, introduced above, we conducted simulations to investigate the interaction between two weak spots of length a separated by a distance d (Fig. 8). Different values of d were simulated and the effect of this parameter on the release angle was examined. (Note that for d = 0 the problem is the same as that presented previously with one weak spot of half-length a.)

Fig. 8. Geometry of the system with two weak spots of length a separated by a distance d.

Figure 9 shows the evolution of the release angle as a function of distance, d. Typical displacement profiles illustrating the behavior for different values of d/Λ are also shown. If the distance between the weak spots is high compared with A (typically for d/Λ≳10), then the displacement profiles generated by the weak spots do not interact. The release thus occurs for the same angle as for the case of only one weak spot of total length a. If the distance, d, between weak spots decreases to values of the same order as Λ (i.e. if 1 < d/Λ≲10), the displacement profile retains a bimodal shape but the two peaks progressively coalesce. As a consequence of this interaction, the release angle, θ _r, progressively decreases as d decreases. Empirically, the evolution of θ _r for d/Λ > 1 can be adjusted by an exponential function: θ _r = θ _∞(1 − γe^−d /^(kΛ)). As expected, the values of θ _∞ and γ depend on the slab depth, h, and weak-spot length, a (θ _∞ ≈ 39.1° and γ ≈ 0.17 in the present case), but the constant, k, is independent of these parameters (k ≈ 3). Finally, if the distance between weak spots is less than the characteristic length (d/Λ < 1), the release angle increases as d decreases and the displacement profile becomes unimodal. This indicates that, in this case, the slab does not ‘feel the effect’ of the cohesive zone between the two weak spots and only‘sees an equivalent weak spot’ of length L ≈ 2a + d.

Fig. 9. Evolution of the release angle, θ _r, with the inter-distance between the two weak spots, d (normalized by the characteristic length, A), for weak-spot lengths a = 4.5 m. The black dots are the results of the finite-element model. The red curve is an exponential adjustment for d/Λ > 1. Three different zones are distinguished: d < Λ, Λ < d≲10Λ and d≳10Λ, and the typical displacement profiles, recorded a few time-steps before release, in each of these zones are represented.

Hence, it appears that the interaction and progressive bridging between the two weak spots is primarily controlled by the characteristic length, A. Physically, this characteristic length represents the typical distance over which the stress redistribution induced by slab elasticity is felt. As illustrated in Figure 9, this stress redistribution actually amounts to smoothing out the effect of the structural heterogeneity of the weak layer as soon as the typical variations of this heterogeneity occur over distances less than A. Hence, A can be viewed as a characteristic smoothing length associated with slab elasticity. More generally, it should be noted that this parameter Λ appears to be the main length scale of the system. In particular, we checked that the softening length, δ, involved in the quasi-brittle weak-layer constitutive law, has essentially no influence on the results as long as it remains much smaller than Λ.

4. Results: Influence of Weak-Layer Heterogeneity

4.1. Simulation protocol

We now consider the case of a spatially heterogeneous weak layer, as described in Section 2.4. We conducted simulations for different values of the slab depth, h, varying between 0.25 and 4 m and different values of the correlation length, ∊, varying between 0.5 and 40 m. For each couple, (h, ∊), 100 simulations with different realizations of the heterogeneity were performed. In each of the simulations, the release angle, θ _r, was determined according to the release criterion given by Eqn (5). For reasons that are developed below, the results are primarily presented in terms of the release factor, F, defined as

(17)

4.2. Release angle and release factor distributions

Figure 10 shows the influence of slab depth, h, and correlation length, e, on the cumulative distributions of the release factor, F. Note that these distributions can also be interpreted in terms of the release angle, θ _r. First we observe that all the distributions obtained can be well represented by Gaussian laws, which can be interpreted as a consequence of the Gaussian nature of the cohesion heterogeneity. As shown in Figure 10a, the average and the variance of the release factor distributions decrease with the slab depth, h. In addition, the average appears to be approximately independent of the correlation length, ∊, while the variance increases with ∊ (Fig. 10b). These results will now be described in more detail.

Fig. 10. Cumulative distributions of release angle θ _r (top scale) and release factor F (bottom scale). (a) ∊ = 0.5 m and various values of slab depth, h; (b) h = 1 m and various values of ∊. Note that according to Eqn (17), the top scale represented is nonlinear in terms of release angle, θ _r.

4.3. Average release factor

In a homogeneous case, the release factor, F _h, is expected to decrease with h according to F _h = 〈c〉/ (ρgh). As shown in Figure 11a, the numerical results appear to closely follow this prediction. In detail, however, it can be noted that the average release factor is always slightly lower than F _h, the difference tending to vanish as the slab depth, h, increases. The same slight difference with F _h is seen in Figure 11b, where it appears to increase with increasing correlation length a for ∊ < 10 m. Recall that for ∊ > 10 m the results begin to be influenced by finite-size effects. These small discrepancies between the results and the theoretical homogeneous value, F _h, are due to the heterogeneity and the presence of local cohesion minima. However, globally, we can conclude that the average release factor (F) (and the average release angle〈θ _r〉) is almost unaffected by the weak-layer heterogeneity.

Fig. 11. Evolution of the average release factor, 〈F〉, and release factor variance, , as functions of (a) slab depth, h, for ∊ = 0.5 m and (b) heterogeneity correlation length, ∊, for h = 1 m. represents the release factor variance that would be obtained in the case of a completely rigid slab.

4.4. Variability and heterogeneity smoothing

Figure 11 a shows that the release factor variance decreases with slab depth, h, as a power law. The associated exponent is slightly smaller than −2 (∼ −2.16). In addition, this variance appears to be significantly smaller than the variance, , that would be observed if the stress field in the weak layer exactly followed the heterogeneity variations (the case of a completely rigid slab). This illustrates the smoothing of the heterogeneity due to the elastic redistribution of stresses in the slab.

Following Section 3, the stress redistribution effects induced by slab elasticity are characterized by the smoothing length, . Hence we can assume that the ratio , can be expressed only in terms of the ratio ∊/Λ. As shown in Figure 12a, all the data points corresponding to obtained from the simulations effectively collapse on a single master curve when plotted in terms of ∊/Λ. The following power-law expression provides a good fit to the data:

Fig. 12. Ratio between release factor variance, , and infinitely rigid slab variance , as a function of the ratio between correlation length, ∊, and elastic smoothing length, Λ, for all (h, ∊) couples (∼5000 simulations). (b) Evolution of the elastic smoothing length, Λ, with slab depth, h.

(18)

with δ = 5.45 × 10⁻². Note, however, that one can expect if ∊/Λ → ∞. Hence, the power law given by Eqn (18) is not expected to remain valid in this limit. In addition, Figure 12b shows that the elastic smoothing length, A, slightly increases with slab depth, h, in our simulations. This is consistent with the slightly less than −2 power-law exponent observed for the evolution of with h in Figure 11.

We can thus conclude that, at the limit of low values of correlation length, ∊, and/or high values of the smoothing length, Λ (and thus high values of slab depth, h), and the system behaves as in a homogeneous case. This also explains why the difference between the average release factor, 〈F〉, and the theoretical homogeneous value, F _h, decreases when h increases or ∊ decreases, as noted in Figure 11a and b. For large correlation lengths, since the effect of ∊ dominates the effect of h in Eqn (18), very thick slabs can be released even for moderate slope angles, which corresponds to the so-called knock-down effect (Kronholm and Schweizer, 2003; Schweizer and others, 2008).

5. Comparison With Field Data

5.1. La Plagne release depth data

La Plagne, in the French Alps, is one of the largest ski areas in the world, covering 100 km² with 225 km of ski tracks. Ski patrollers provided us with release depth data from 14 391 avalanches collected during the winters of 1998–2010. Since avalanche depths cannot always be directly measured, data come from a mix of visual estimates and precise measurements. From the complete database, 369 naturally released slab avalanches were extracted. These data have been analyzed in detail by Gaume and others (2012). Note that the same data were also used by Failletaz and others (2004), but with fewer avalanches (only three winters). These analyses showed that the release depth cumulative distribution at La Plagne seems to decrease as a power law for large slab depths (h > 0.7 m, corresponding to cumulative exceedance probability lower than ∼10%). Similar power-law trends have also been reported at other locations (Rosenthal and Elder, 2002; McClung, 2003). We note however that, due to the error associated with the data, the value of the power-law exponent is poorly constrained and strongly dependent on the cut-off considered for the power law. Typically, in agreement with McClung (2003), exponents in the [− 3; − 5] range provide a good fit to the data.

Our objective here is to examine whether our mechanical model is capable of reproducing these release depth data for both the core and the tail of the distribution. This comparison first requires computing the release depth distribution predicted by the model, which can be obtained in two phases: (1) the release factor distributions presented above have to be inverted to obtain release depth distributions for fixed angle values; (2) these release depth distributions must then be integrated over all slopes, since data mix avalanche paths of various slope angles. Lastly, to be compared with data, the release depth probability obtained from the mechanical model has to be combined with the local snowfall probability.

5.2. Release depth distributions obtained from mechanical model

Inversion of release factor distributions

We have shown that the distributions of the release factor, F, are normally distributed with mean 〈F〉 ≈〈c〉 / (ρgh) and variance , with f(∊/Λ) given by Eqn (18). In addition, since A varies only slightly in our results (see Fig. 12b), we assume it to be constant in the following. This approximation does not significantly influence the results to be presented, but it allows us to obtain analytical solutions. Hence, the variance, σ ² _F , can be written as

(19)

with f(∊) ≈ κ∊ ^2/3, and κ = δΛ^−2/3 assumed constant. Finally, the probability density of having a release factor, F, for a given slab depth, h, is given by

(20)

with C_μ = 〈c〉 = (ρg) and

The Mohr–Coulomb rupture criterion, which controls avalanche release in the simulations, can be written in terms of the release factor, F = sin θ _r − μcos θ _r, as

(21)

Hence, slab depth, h, and release factor, F, play similar roles in this criterion. Consequently, it can be shown that, if the probability density of release factor F for a given slab depth h is p(F|h) = g(F, h), then the probability density of h for a given value of F is p(h|F) = g(h, F). Equation (20) can thus be inverted to

(22)

A more detailed and rigorous demonstration of this inversion is provided in Appendix B.

Integration over all slopes

For the sake of simplicity, we chose, in this study, to consider a uniform probability distribution for the slope factor, F = sinθ − μcosθ, between F _min and F _max: p(F) = 1 / (F _max− F _min). Once again, this assumption enables us to obtain analytical expressions for the integrated release depth distribution, p _m(h):

(23)

From Eqn (22), we obtain

(24)

with

(25)

and

(26)

where we define U _min = (hF _min − C_μ ) / C_σ and _{U max} = (hF _max − C_μ ) / C_σ . In the following, without loss of generality, we assume F _min = 0 and F _max = 1.

5.3. Coupling mechanical and snowfall distributions

Gaume and others (2012) have shown that the global avalanche release depth probability, p _r(h), resulting from the coupling between the mechanical model presented above and snowfall distributions, can be related to p _m(h) as follows:

(27)

where p _sf(h _sf ≳ h) is the probability of having a snowfall whose thickness, h _sf, is higher than the depth, h, and C a normalization constant given by . This coupling relation expresses that the amount of snowfall represents a limiting factor weighting the mechanical probability density, p _m(h), derived from the stability criterion.

To define the snowfall distribution, p _sf(h _sf ≳ h), Gaume and others (2012) considered the 3 day snowfall annual maxima in La Plagne (MeteoFrance data: daily measurements from 1966) at the average altitude of 2200 m. These maxima follow a generalized extreme value (GEV) distribution so that

(28)

where μ _sf, σ _sf and ξ_sf are, respectively, the location, scale and form parameters. These parameters assume the values μ _sf = 0.98 m, σ _sf = 0.21 m and ξ _sf = 0.214.

5.4. Result of the coupling and sensitivity analysis

As shown by Gaume and others (2012) (see also Fig. 13), it is possible to find a set of mechanical parameters for which the coupled model described by Eqns (24), (27) and (28) provides a very good adjustment to La Plagne release depth data. The model effectively predicts a power-law behavior of the cumulative exceedence distribution for large slab depths, in good agreement with the empirical distribution, and also accounts well for the data corresponding to lower release depths. In spite of the various assumptions involved, this model is, to our knowledge, the first capable of reproducing release depth data with such good accuracy. To complement the results shown by Gaume and others (2012), we present below a detailed sensitivity analysis of the predicted distribution of the main mechanical parameters of the model.

Fig. 13. Slab release depth distributions predicted by the coupled model (Eqn (27)) for different values of average cohesion, 〈c〉, and comparison with field release depths from La Plagne. The numerical results have been obtained for a cohesion standard deviation σ_c = 0.3 kPa, a correlation length ∊ = 2 m, F _min = 0 and F _max = 1, the other parameters being the same as in Section 4. (a) Cumulative exceedance probability on logarithmic scale; (b) χ² of the model (χ² = ∑(p _data − p _model)² / p _model, where p _data and p _model are the cumulative exceedance probabilities derived from the data and from the model, respectively).

Figure 13 shows the comparison between data and the coupled model for different values of the average cohesion, (c). The goodness of fit (Fig. 13b) shows a pronounced minimum for 〈c〉 = 0.6 kPa, which indicates that the agreement between model and data strongly depends on this average cohesion value. In particular, the depth, h _m, below which no avalanche can occur is strongly dependent on 〈c〉. This depth can be approximated by (Gaume and others, 2012)

(29)

With the value 〈c〉 = 1 kPa retained in Sections 3 and 4, h _m is slightly overestimated compared to the data. For 〈c〉 = 0.6 kPa, a value still fully consistent with existing studies (Föhn and others, 1998; Jamieson and Johnston, 2001), an excellent agreement between the coupled model and the data is obtained.

The influence of the cohesion standard deviation, σ_c , is presented in Figure 14. First, it can be noted that σ_c plays a less significant role than 〈c〉 in the goodness of the fit. As shown in Figure 14b, values of σ_c in the 300–500 Pa range give similar results for constant values of the other parameters. In fact, the standard deviation, σ_c , mainly influences the curvature of the coupled cumulative exceedence distribution around the cut-off, h _m. Figure 14a also displays the distribution that would be obtained in the case of a completely rigid slab, which poorly reflects the data. Even if the tail of this distribution could be fitted to the data by tuning the other parameters, there would be no way to properly fit the core of the distribution. This highlights the major importance of the elasticity of the slab and of stress redistribution effects.

Fig. 14. Slab release depth distributions predicted by the coupled model (Eqn (27)) for different values of the standard deviation, σ_c , and comparison with field release depths from La Plagne. The numerical results have been obtained for a cohesion 〈c〉 = 0.6 kPa, a correlation length ∊ = 2 m, F _min = 0 and F _max = 1, the other parameters being the same as in Section 4. (a) Cumulative exceedance probability on logarithmic scale; (b) χ² of the model.

Finally, the influence of the correlation length, ∊, is shown in Figure 15. As with the standard deviation, σ_c , the correlation length, ∊, mainly modifies the curvature of the coupled cumulative exceedence distribution around h _m. Globally, the fit to the data remains good for correlation length values in the 0.5–15 m range.

Fig. 15. Slab release depth distributions predicted by the coupled model (Eqn (27)) for different values of the correlation length, ∊, and comparison with field release depths from La Plagne. The numerical results have been obtained for a cohesion 〈c〉 = 0.6 kPa, a cohesion standard deviation σ_c = 0.3 kPa, F _min = 0 and F _max = 1, the other parameters being the same as in Section 4. (a) Cumulative exceedance probability on logarithmic scale; (b) χ² of the model.

To conclude, in the range of realistic mechanical parameters for snow (for which stress redistribution effects play an important role), it turns out that the average cohesion, 〈c〉, has the most significant influence on slab avalanche release distributions predicted by the coupled model. Hence, provided GEV parameters are known, the adjustment of the model to the data essentially amounts to a one-parameter fit. Note also that the value of _{F max} which has been set to 1 in the previous results, in fact plays a role essentially similar to that of 〈c〉 (see Eqn (29)). Hence, changing the value of _{F max} would result in straightforward modifications of the best-fit value found for 〈c〉.

6. Conclusion and Perspectives

This paper investigates the influence of weak-layer heterogeneity on slab avalanche release using a finite-element model. A shear-softening interface underlying an elastic slab is modeled and the system is loaded by increasing the slope angle until failure and avalanche release. After validating the model for the case of a nil-cohesion weak spot, the effect of a heterogeneous weak-layer cohesion field was studied. The heterogeneity is represented through a Gaussian distribution, with a spherical covariance function characterized by a spatial correlation length. Release angle distributions were analyzed and a heterogeneity smoothing effect due to redistributions of stresses by elasticity of the slab was highlighted. This smoothing effect induces a reduction of the release angle variance compared to the case of a fully rigid slab. However, the average release angle is almost unaffected by this effect. The presented results show that the smoothing intensity critically depends on the ratio between the correlation length, ∊, and a characteristic elastic length of the system, Λ. Further work is required, however, to fully unravel the possible interplay between a and the cohesion variance, , on this smoothing effect. To be compared with field data, the obtained release angle distributions were inverted, yielding a release depth distribution integrated over all slopes. Coupling this mechanical distribution with the distribution of 3 day extreme snowfalls, we were able to reproduce with excellent accuracy field data from 369 natural slab avalanches. A detailed sensitivity analysis showed that this agreement is obtained with only one adjustable parameter, namely the average cohesion, 〈c〉. The mechanically based probabilistic model thus fulfills the objectives of the study, namely the evaluation of avalanche release depth distributions in any potential release zone, as soon as meteorological data are available. In the future, a straightforward extension to the three-dimensional case will be developed to predict distributions of avalanche release volumes. However, before being used in an operational context, additional tests on other datasets and in other locations need to be performed to further validate the model. Finally, let us recall that, with the parameters used in our model, the crown fracture always occurs at particular morphological features (e.g. ridges, rocks and trees) since the heterogeneity is not sufficient to directly trigger tensile rupture within the slab. Thus, another interesting perspective for future work is to study the release depth distributions obtained with different sets of parameters leading to other types of failure mechanisms.