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 (Reference Keylock, McClung and MagníssonKeylock and others, 1999; Reference Eckert, Parent, Naaim and RichardEckert and others, 2008, Reference Eckert, Naaim and Parent2010). 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 Reference Gaume, Chambon, Eckert and NaaimGaume 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 (Reference Verpeaux, Charras, Millard, Fouet, Ladevèze and OhayonVerpeaux and others, 1988). The resolution procedure used (‘Pasapas’; Reference Charras and Di PaolaCharras 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 (Reference McClungMcClung, 1979; Reference Föhn, Camponovo and KrüsiFöhn and others, 1998; Reference McClung and SchweizerMcClung 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 (Reference McClungMcClung, 1979; Reference Bažant, Zi and McClungBažant and others, 2003; Reference Fyffe and ZaiserFyffe and Zaiser, 2004, Reference Fyffe and Zaiser2007; Reference Mahajan and JoshiMahajan and Joshi, 2008; Reference Mahajan, Kalakuntla and ChandelMahajan and others, 2010; Reference Gaume, Chambon, Naaim, Eckert, Bonelli, Dascalu and NicotGaume 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° (Reference De MontmollinDe Montmollin, 1978; Reference Van Herwijnen and HeierliVan 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 Reference McClungMcClung (1977), both the characteristic peak and softening displacement, u _p and δ, are taken equal to 2 mm (Reference Bažant, Zi and McClungBažant and others, 2003; Reference Fyffe and ZaiserFyffe and Zaiser, 2004, Reference Fyffe and Zaiser2007). Hence, the displacement to reach the residual stress u _r = u _p + δ = 4 mm. Note that, according to recent studies (Reference McClungMcClung, 2009, Reference McClung2011), 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 (Reference MellorMellor, 1975; Reference McClungMcClung, 1977; Reference NaritaNarita, 1980; Reference Navarre, Taillefer, Danielou, Brugnot, Brun, Gubler, Norem and ReynaudNavarre and others, 1992; Reference SchweizerSchweizer, 1998) and field measurements (Reference RochRoch, 1966; Reference De MontmollinDe Montmollin, 1978; Reference Jamieson and JohnstonJamieson and Johnston, 1990; Reference Föhn, Camponovo and KrüsiFöhn and others, 1998; Reference McClung and SchweizerMcClung 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 (Reference MellorMellor, 1975; Reference Camponovo and SchweizerCamponovo and Schweizer, 2001; Reference SchweizerSchweizer, 1999).

2.4. Spatial heterogeneity

The spatial heterogeneity of the weak layer is modeled through a stochastic distribution of cohesion, c. Following Reference Jamieson and JohnstonJamieson and Johnston (2001) and Reference Kronholm and BirkelandKronholm 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 Reference Kronholm and BirkelandKronholm 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 (Reference Chilès and DelfinerChilè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 (Reference Jamieson and JohnstonJamieson and Johnston, 2001; Reference Schweizer, Kronholm, Jamieson and BirkelandSchweizer and others, 2008; Reference Bellaire and SchweizerBellaire and Schweizer, 2011) relevant for weak-snow layers. Reference Schweizer, Kronholm, Jamieson and BirkelandSchweizer 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; Reference MellorMellor, 1975; Reference Jamieson and JohnstonJamieson and Johnston, 1990; Reference SigristSigrist, 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 (Reference Gaume, Chambon, Naaim, Eckert, Bonelli, Dascalu and NicotGaume 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.1. A single weak spot

Analytical solution

Here we follow the approach of Reference Chiaia, Cornetti and FrigoChiaia 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.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 (Reference Kronholm and SchweizerKronholm and Schweizer, 2003; Reference Schweizer, Kronholm, Jamieson and BirkelandSchweizer and others, 2008).

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 Reference Gaume, Chambon, Eckert and NaaimGaume and others (2012). Note that the same data were also used by Reference Faillettaz, Louchet and GrassoFailletaz 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 (Reference Rosenthal, Elder and StevensRosenthal and Elder, 2002; Reference McClungMcClung, 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 Reference McClungMcClung (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

Reference Gaume, Chambon, Eckert and NaaimGaume 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), Reference Gaume, Chambon, Eckert and NaaimGaume 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 Reference Gaume, Chambon, Eckert and NaaimGaume 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 Reference Gaume, Chambon, Eckert and NaaimGaume 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 (Reference Gaume, Chambon, Eckert and NaaimGaume 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 (Reference Föhn, Camponovo and KrüsiFöhn and others, 1998; Reference Jamieson and JohnstonJamieson 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〉.