Modeling of crack propagation in weak snowpack layers using the discrete element method

. Dry-snow slab avalanches are generally caused by a sequence of fracture processes including (1) failure initiation in a weak snow layer underlying a cohesive slab, (2) crack propagation within the weak layer and (3) tensile fracture through the slab which leads to its detachment. During the past decades, theoretical and experimental work has gradually led to a better understanding of the fracture process in snow involving the collapse of the structure in the weak layer during fracture. This now allows us to better model failure initiation and the onset of crack propagation, i.e., to estimate the critical length required for crack propagation. On the other hand, our understanding of dynamic crack propagation and fracture arrest propensity is still very limited. To shed more light on this issue, we performed numerical propagation saw test (PST) experiments applying the discrete element (DE) method and compared the numerical re-sults with ﬁeld measurements based on particle tracking. The goal is to investigate the inﬂuence of weak layer failure and the mechanical properties of the slab on crack propagation and fracture arrest propensity. Crack propagation speeds and distances before fracture arrest were derived from the DE simulations for different snowpack conﬁgurations and mechanical properties. Then, in order to compare the numerical and experimental results, the slab mechanical properties (Young’s modulus and strength) which are not measured in the ﬁeld were derived from density. The simulations nicely reproduced the process of crack propagation observed in ﬁeld PSTs. Finally, the mechanical processes at play were analyzed in depth which led to suggestions for minimum column length in ﬁeld PSTs.


Introduction
Dry-snow slab avalanches result from the failure of a weak snow layer underlying cohesive slab layers.The local damage in the weak layer develops into a crack which can expand if its size exceeds a critical length or if the load exceeds a critical value.Finally, crack propagation leads to the tensile fracture of the slab and ultimately, avalanche release (McClung, 1979;Schweizer et al., 2003).During the past decade, our understanding of the fracture process in snow has gradually evolved through the development of new theories as well as various field observations and experiments.The propagation saw test (PST), concurrently developed in Canada (van Herwijnen and Jamieson, 2005;Gauthier and Jamieson, 2006) and Switzerland (Sigrist and Schweizer, 2007), consists in isolating a snow column and initiating a crack of increasing length in the weak layer with a snow saw until the onset of rapid self-propagation of the crack.The PST allows observers to determine the critical crack length and evaluate crack propagation propensity.This field method has highlighted the importance of slab bending (due to the collapsible nature of weak snow layers) on crack propagation (e.g.van Herwijnen et al., 2010;van Herwijnen and Birkeland, 2014).On the other hand, theoretical and numerical models, based on fracture mechanics or strength of material approaches, were developed to investigate crack propagation and avalanche release (McClung, 1979;Chiaia et al., 2008;Heierli et al., 2008;Gaume et al., 2013Gaume et al., , 2014b)).While substantial progress has been made, application with regard to avalanche forecasting or hazard mapping is still hindered in Published by Copernicus Publications on behalf of the European Geosciences Union.© ASARC from Jamieson and Schweizer (2000).(b) Zoom on a surface hoar crystal © ASARC.
part by our lack of understanding of the dynamic phase of crack propagation.For instance, based on practitioners' experience, it is not uncommon to perform PST field measurements with widespread crack propagation in 1 day, while a few days later, with seemingly very few changes in snowpack properties, cracks will no longer propagate (Gauthier and Jamieson, 2008).Thus far, there is no clear theoretical framework to interpret such observations, and it is not clear how and which snowpack properties affect dynamic crack propagation.This limitation is due to the complex microstructure of snow and its highly porous character (Fig. 1) which are not taken into account in the continuous approaches previously mentioned.
In this paper, numerical experiments of the propagation saw test (PST) are performed by applying the discrete element (DE) method which allows us to mimic the high porosity of snow.The goal is to investigate the influence of weak layer failure and the mechanical properties of the slab on crack propagation.In the first section, field data as well as the proposed model are presented.Then, crack propagation speed and distance before fracture arrest are derived from the DE simulations using the same method as for the field experiments (particle tracking).In a parametric analysis, we show the influence of single system parameters on the crack propagation speed and distance.Finally, the interdependence of snowpack properties is accounted for in order to compare numerical and experimental results and the mechanical processes leading to fracture arrest are analyzed.

PST field data
Since the winter of 2004-2005, we collected data from 121 PST experiments at 46 different sites in Canada, USA and Switzerland (van Herwijnen and Jamieson, 2005;van Herwijnen and Heierli, 2009;van Herwijnen et al., 2010;Bair et al., 2012;van Herwijnen and Birkeland, 2014;Birkeland et al., 2014).At each site, we collected a manual snow pro- file and conducted one or several PSTs according to the procedure outlined in Greene et al. (2004).In many cases we used longer (than standard) columns to allow us to better investigate crack propagation.After column preparation, we inserted black plastic markers into the pit wall and used a digital camera on a tripod to make a video recording of the PST (Fig. 2).We used a particle tracking velocimetry (PTV) algorithm to analyze the motion of the markers and thus the displacement of the snow slab above the weak layer (Crocker and Grier, 1996).In this way, the position of the markers in each video frame can be determined with a mean accuracy of 0.1 mm.The displacement of a marker is then defined as the movement relative to its initial position, that is, the average position of the marker prior to movement.For propagating cracks, there is a delay between the vertical displacement of subsequent markers.A typical displacement-time evolution for a propagating crack is shown in Fig. 3 for four slab markers.As explained in van Herwijnen and Jamieson (2005), the time delay between the onset of movement between markers is proportional to the distance between the markers and was used to calculate the propagation speed c of the fracture (see also Appendix A).

Motivation and objectives
Discrete element (DE) modeling (Cundall and Strack, 1979) allows computing the motion of a large number of small grains by solving dynamic equations for each and defining a contact law between the grains.In addition, the DE method allows assessing mechanical quantities such as stress, displacement, deformation rate, porosity, etc., computed over representative elementary domains at each material point within the sample.Experimentally, this would be an impossible task.Hence, using DE, the mechanical and rheological behavior of the material can be explored locally, regardless of the spatial heterogeneities possibly displayed by the structure of the material and its mechanical quantities.This method can thus help to better understand physical processes at play in granular assemblies.The DE method has been widely used to study the flow of granular materials within industrial (e.g.Chaudhuri et al., 2006;Sarkar and Wassgren, 2010) or environmental applications such as avalanche dynamics (e.g.Rognon et al., 2008;Faug et al., 2009).However, to our knowledge, discrete elements have never been used to model crack propagation in layered systems or to describe slab avalanche release processes.The latter processes are generally modeled under a continuum mechanics framework, using methods such as finite elements (Podolskiy et al., 2013, and references therein).These methods can be used to assess the stability of a layered snow cover, i.e., determine the conditions of failure occurrence and the onset of crack propagation.However, they are less suited to study what occurs after failure, i.e., during the dynamic phase of crack propagation, due to a lack of relevant constitutive models for the WL, including softening which induces a loss of objectivity with respect to the mesh in dynamic problems.

The
The objective of the proposed approach is to use the discrete element method (DE) to study the dynamic phase of crack propagation in a weak snowpack layer below a cohesive slab.The DE method is adequate for our purpose because (1) no assumption needs to be made about where and how a crack forms and propagates, and (2) the model material is inherently discontinuous and well adapted to dynamic issues.We will show that this method allows us to capture all the main physical processes involved in the release of dry-snow slab avalanches, namely the complex mechanical behavior of the weak layer and the interplay between basal crack propagation, slab bending, and slab fracture.
However, an important preliminary issue to address concerns the scale of the considered model.In the weak layer, we intend to represent, through a simplified description, the particular collapsible and highly porous microstructure of the snow in order to be able to reproduce the main features of the failure envelope of this material.As will be shown, we achieve this by using triangular shapes of centimeter size.To account for the possible breakage of these elements, they consist of small cohesive grains of size r wl .In the slab, on the contrary, due to computational costs, it would be completely unrealistic to try to account for the complete microstructure of the snow at the scale of a real slope or even at the scale of a field test such as the PST.The slab is thus modeled as a continuous material with an elastic-brittle constitutive behavior.Yet, similar to what is classically done for concrete (Meguro and Hakuno, 1989;Kusano et al., 1992;Camborde et al., 2000;Hentz et al., 2004), the response of this layer to the dynamic propagation of failure in the WL is also computed with the DE method.In that case, however, the considered elements (grains of size r) have no physical meaning and should only be regarded as entities of discretization similar to the mesh size in FE models.The contact parameters need to be properly calibrated (through classical biaxial tests, for instance) in order to recover the correct macroscopic properties of the material.Other continuous methods, such as FE, could have been used to simulate the slab, but -apart from avoiding the non-trivial issue of coupling DE and FE regions -the use of DE is well suited to represent the large deformations involved in the bending of the slab and the spontaneous formation of the tensile crack.
To summarize, we contend that, unlike in other DE applications which are at the scale of the microstructure (e.g.Cundall (1989); Iwashita and Oda (2000) for frictional granular materials or Hagenmuller et al. (2015) for snow), the grains involved in the DE model developed in this study should not be regarded as snow grains, and that both r wl and r are only discretization scales whose choice will result from a compromise between resolution and computational cost as classically done to model concrete fracture (Hentz et al., 2004).We consider here a meter-scale model where the advantage of the DE method is its ability to mimic the poorly known mechanical response of the weak layer and to account for the different modes of failure displayed by snow (shear, compression, tension).The only microstructural scale directly accounted for is the size of the triangular elements in the weak layer, which are on the same order as the weak layer thickness, since it is a necessary ingredient for reproducing the particular mechanical behavior of this layer.

Software
The discrete element simulations were performed using the commercial software PFC2D (by Itasca), which implements the original soft-contact algorithm described in Cundall and Strack (1979).

Simulated system
The simulated system (see Fig. 4a) is two dimensional and is composed of a completely rigid (fixed) basal layer, a WL of thickness D wl and a slab of thickness D which were varied in the simulations.The slab is composed of grains of radius r = 0.01 m with a primitive square packing.The porosity of the slab is equal to 21 %.Hence the density of the slab ρ can be adjusted by changing the grain density ρ p, s (varied in the simulations).The WL is composed of grains of radius r wl = r/2 with a complex packing of collapsible triangular forms (see Fig. 4b) aimed at roughly representing the porous structure of persistent WLs such as surface hoar or depth hoar.The porosity of the WL is 70 % and the density of the WL grains is ρ p, wl = 400 kg m −3 , leading to a WL density ρ wl = 120 kg m −3 .The length of the system (column length) is L = 2 m and the slope angle is denoted ψ.As mentioned above, the numerical grains are not intended to represent the real snow grains which are obviously smaller and have a different density.Nevertheless, as will be shown, this set up allows to capture the main features observed in field PSTs.

Loading
The loading is applied by gravity and by advancing a "saw" (in red on Fig. 4a) at a constant velocity v saw = 2 m s −1 through the weak layer.This saw is composed by rigid walls and has approximately the same thickness as field saws h saw = 2 mm.The saw velocity was chosen relatively high to decrease the computational time, but lower than the lowest crack propagation speed observed in the field so as to correctly distinguish crack propagation from saw movement.

Contact law
The cohesive contact law used in the simulations is the PFC parallel bond model represented schematically in Fig. 5a.The cohesive part acts in parallel to the classical linear contact law (Radjai et al., 2011;Gaume et al., 2011).For the linear component (in grey in Fig. 5a), the normal force is the sum of a linear elastic and of a viscous contribution (springdashpot model), and the tangential force is linear elastic with a Coulombian friction threshold.The corresponding mechanical parameters, namely the normal and shear stiffness k n and k s (elasticity parameters), the restitution coefficient e (viscous parameter) and the friction coefficient µ are summarized in Table 1.The value of the normal stiffness k n was chosen in such a way that the normal interpenetration at contacts are kept small, i.e., to work in the quasi-rigid grain limit (da Cruz et al., 2005;Roux and Combe, 2002).Concerning the normal restitution coefficient e, we verified that the results presented below, and more generally all the macroscopic mechanical quantities obtained from the simulations, are actually independent of this parameter (in the range 0.1-0.9).This is due to the presence of the cohesive part of the contact law (see details below).Indeed, the restitution coefficient might have a strong influence for cases in which new contacts and collisions occur.In our case, the results are mostly driven by bond breaking which explains why e has no influence on the results.The cohesive component (in black in Fig. 5a) can be envisioned as a point of glue with constant normal and shear stiffness k n and k s acting at the contact point.This bond has a specified shear and tensile strength σ periphery are calculated via beam theory according to: where F n and F s are the bond normal and shear forces, M is the bending moment, r b the bond radius, A = π r 2 b the bond area and I = π r 4 b /4 its moment of inertia.If the tensile stress exceeds the bond tensile strength, the bond breaks and both the normal and shear contact forces are set to zero (Fig. 5b).If the shear stress exceeds the bond shear strength, the bond also breaks (Fig. 5c) but the contact forces are not altered, provided that the shear force does not exceed the friction limit, and provided that the normal force is compressive.The bond can also break if the bending moment exceeds σ b t I /r b (Fig. 5d).The ranges of parameters used for the bond model are summarized in Table 2.
As the two components are acting in parallel, the total stiffness of the bonded contact is equal to the sum of the contact stiffness k n and the bond stiffness k n : (3) Finally, in our case of a square grain assembly, the Young's modulus can be derived analytically from k n and k n according to: The contact stiffness k n was kept constant at 1 × 10 4 N m −1 and k n was varied between 1 × 10 3 and 1.5 × 10 6 N m −1 (for the slab) leading to realistic values of the Young's modulus E between 0.35 and 50 MPa.It was verified from biaxial tests that the macroscopic (bulk) Young's modulus of the slab effectively follows this relationship due to the specific (squared) structure of the slab.

Time step and elastic waves
The time step was computed classically as a function of the grain properties according to t = √ m/k ≈ r √ ρ/E where m, ρ and r are the smallest grain mass, density and radius and k and E the largest contact/bond stiffness and Young's modulus, respectively.The choice of this time step insures the stability of the algorithm and that the crack propagation speed is lower than the speed of the elastic waves in the sample.The order of magnitude of the speed of elastic waves in the sample is c e ≈ √ E/ρ.

Simulation protocol and illustration
First, the macroscopic properties of the slab have to be determined as a function of the microscopic properties of the bond.For the slab, bi-axial tests were carried out which allowed to validate that for a squared assembly, the macro-  scopic (bulk) Young's modulus depends on bond stiffness according to Eq. ( 4).
For the weak layer, similarly, simple loading tests were carried out to compute the macroscopic failure criterion (mixed-mode shear-compression) of the WL as a function of the bonds of WL grains (Gaume et al., 2014a).The failure envelope of the WL is strongly anisotropic as shown in Fig. 6.This failure envelope shows, for realistic values of the slope angle, a much lower strength in shear than in compression as well as a decrease of the shear strength with increasing normal stress.This type of behavior is similar to that reported in recent laboratory (Reiweger et al., 2015) and field (Chandel et al., 2014) experiments on persistent weak snow layers.Hence, the chosen WL structure allows to have different modes of failure (tension, shear, compression and mixed-mode) such as observed in real weak layers and thus has the essential characteristics to model the processes of slab avalanche release.
Then, PST simulations were performed.An illustration of a simulation result highlighting the displacement wave of the slab is shown in Fig. 7 and the associated vertical displacement y is represented in Fig. 8.The critical length is denoted a c and corresponds to the translation length of the saw required to obtain a self-propagating crack.
In order to determine the crack propagation speed, purely elastic simulations (infinite tensile and shear strength of the bonds between slab grains) were carried out.The propagation speed was computed using the same method as for field PSTs by analyzing the vertical displacement wave of the slab  (van Herwijnen and Jamieson, 2005).This procedure is presented in more detail in Appendix A.
The propagation distance was computed by taking into account the possible failure of the slab by setting finite values to the tensile and shear strength of the slab (σ t , σ s ).We define the propagation distance as the distance between the left wall of the system and the location where the tensile crack opens at the slab surface, as shown in Fig. 9.This measure of the propagation distance differs from that defined by Gauthier and Jamieson (2006) who defined it as the distance between the point of onset of crack propagation in the WL and the point of slab failure.However, we argue that the propagation distance, as we define it, is a more suitable measure since this is the one that influences the stresses in the slab and thus fracture arrest propensity.Similarly, for a real avalanche, the important distance for the bending induced stress and for the avalanche release size is the distance from the very first failure initiation point (whatever the nature of the initial trigger), to the location of the slab tensile failure, and not from the crack tip at the moment of the onset of crack propagation.
For the parametric analysis (Sect.3.2), we performed simulations for which only one system parameter was modified while the other parameters were kept constant.3. However to apply these results to slab avalanche release and in order to compare our results to field data (Sect.3.3), existing relations between snowpack properties were taken into account.Hence, simulations were performed for different slab densities with the Young's modulus varying according to an empirical exponential fit to the data reported by Scapozza ( 2004): and with a tensile strength varying according to a power-law fit to the data reported by Sigrist (2006): with ρ ice = 917 kg m −3 .

Displacement of the slab
The evolution of the vertical displacement y of the slab is represented in Fig. 8.In this figure together with the illustration of the displacement wave of the slab (Fig. 7), one can clearly observe the different processes acting before, during and after crack propagation.First, slab bending occurs prior to the onset of crack propagation and the dynamic propagation phase.These distinct phases (stable bending of the slab and crack propagation) are also clearly visible in the vertical displacement y, as shown in Fig. 8 for four different horizontal positions in the slab.Between 0 and 0.1 s nothing happens, then as the saw advances, the vertical displacement slowly increases.This phase corresponds to the bending of the under-cut part of the slab.Then, for t = 0.25 s approximately, the critical length a c was reached and the displacement increases abruptly, even beyond the saw, corresponding to the dynamic crack propagation phase.After t = 0.3 s, the slab made contact with the broken WL at the left-end of the slab for x = 0 m.After 0.32 s, the entire WL has collapsed leading to a constant vertical displacement of the slab approximately equal to y = 1.8 cm.This displacement is not perfectly equal to the WL thickness because of the grains remaining in the WL.The peak in the displacement around t = 0.38 s is an artefact associated with the movement of the saw after the crack has propagated which does not affect the results that we will present.

Crack propagation speed
For all the simulations carried out, the crack propagation speed varied between 5 and 60 m s −1 .Figure 10a shows that the crack propagation speed c strongly increases with the  3.
Young's modulus E of the slab, from less than 5 m s −1 for very soft slabs (E < 1 MPa) to 40 m s −1 for an almost rigid slab (E ≈ 50 MPa) where the increase levels off.The propagation speed c also strongly and linearly increases with the thickness of the slab D (Fig. 10b), from almost zero for a slab thickness lower than 10 cm to 60 m s −1 for a thickness of 80 cm.Similar to the thickness influence, the propagation speed increases almost linearly with the density of the slab ρ (Fig. 10c) and the slope angle ψ (Fig. 10d).The propagation speed seems not to be influenced by the thickness of the WL (Fig. 10b) as soon as the failure occurs under the same conditions (same critical length).Furthermore, the propagation speed decreases with increasing WL strength (Fig. 10e).This suggests that the crack propagation speed is mostly influenced by the failure conditions (load due to the slab and WL strength) rather than structural parameters such as the WL thickness.
Finally, the speed of the elastic waves in the slab (c e ) is always substantially higher (by about a factor of 10) than the crack propagation speed (Fig. 10).In addition, we can ob-    3.
serve that the speed of the elastic waves in the slab is not a good proxy to explain the trends in the crack propagation speed as shown for instance by the opposite trends with density (Fig. 10c).

Propagation distance
Figure 11 shows propagation distance as a function of different system parameters.Fig. 11a 1 , b 1 and c 1 show the increase of the propagation distance with increasing tensile strength of the slab σ t .This result was expected since a stronger slab requires a larger tensile stress in order to fracture and thus a larger propagation distance is required to obtain sufficient tensile stresses in the slab (induced by bending or by the shear component of the slab's weight additionally in case of ψ = 0).
The influence of the Young's modulus E of the slab is shown in Figs.11a 1 and a 2 .Overall, propagation distance decreases with increasing Young's modulus.Hence, the softer the slab, the lower the fracture arrest propensity.For a tensile strength of 2 kPa (Fig. 11a 2 ), the propagation distance l * sharply decreases from 2 m (column length) to an approximately constant value l * = 0.5 m for E ≈ 2 MPa.Also, Fig. 11a 1 shows that for higher Young's modulus, larger tensile strength values are required to obtain full propagation.The critical length a c for crack propagation was also represented in Fig. 11a  Then, the influence of WL thickness D wl is shown in Fig. 11b 1 and b 2 .The WLs have different thickness but the contact properties were calibrated to obtain the same critical length a c which is equal to 15 cm in this case.The propagation distance l * decreases with increasing WL thickness.For low values of the tensile strength of the slab, the propagation distance is small and almost independent of the WL thickness whereas an important decrease is observed for larger values of the tensile strength.In other words, thicker weak layers result in more slab bending so that slab failure becomes more likely due to high tensile (bending) stress.
Figure 11c 1 and c 2 show the influence of slope angle ψ on propagation distance.Similarly to WL thickness, slope angle seems to have no influence on the propagation distance for low values of the tensile strength.However, for larger values of σ t (σ t > 4 kPa), the propagation distance strongly increases with increasing slope angle ψ.Hence, fracture arrest propensity decreases with slope angle for large values of the tensile strength of the slab, typically higher than 4 kPa.This result is not trivial, since as the slope angle increases, there is a competition between the decreasing slab bending which results in a decrease of the tensile stresses in the slab and an increase of the tensile stresses due to the weight of the slab in the slope parallel direction.Hence, this result suggests that slab bending is the primary process influencing tensile failure of the slab (for homogeneous properties of the system).
Crack propagation distance slightly decreases with slab density as shown in Fig. 11d.For very low slab densities, the critical length a c is relatively high and thus the tensile failure across the slab occurs before the critical length is reached.Then, as the density of the slab increases, the critical length decreases and the propagation distance stabilizes around 0.4 m.
Whereas slab density ρ and slab thickness D have a similar influence on the stability of the system, on the crack propagation speed (Fig. 10b and c) and on the onset of crack propagation, as suggested by the decrease of the critical length a c with both ρ and D (Fig. 11d and e), their influence on fracture arrest propensity differs.Indeed, in contrast to the influence of slab density, the propagation distance strongly increases with increasing slab thickness (Fig. 11e).Hence, the thicker the slab, the lower the fracture arrest propensity.This result can be easily explained using beam theory (Timoshenko and Goodier, 1970) to express the maximum tensile stress in a bending slab which is inversely proportional to the slab thickness D (see Sect. 3.4 or Schweizer et al., 2014).
Finally, crack propagation distance decreases with WL strength (Fig. 11f) for low WL strength values for which the system is close to the overall failure (the critical length is close to zero).However, for higher values of the WL strength, the propagation distance is almost unaffected by this property of the WL.

Comparison with field data
The results of the previous parametric analysis should be interpreted with care since for snow, several of the system parameters are inter-related leading to more complex interactions.For instance, the result about the influence of Young's modulus on the propagation distance might seem contradictory to avalanche observations.Indeed, taken as it is, this result would imply that it is easier to trigger a tensile failure in stiff and thus hard snow than in soft snow.Consequently, hard slabs would result in smaller release areas than soft slabs which is clearly in contradiction with avalanche observations (van Herwijnen and Jamieson, 2007).Hence, even if the result behind Fig. 11a 2 is consistent, from a mechanical point of view, it cannot be directly applied to dry-snow slab avalanche release.To do so, one should take into account the relations between slab density ρ, Young's modulus E and tensile strength σ t according to Eqs. ( 3) and ( 4).Simulations were performed for slab densities ranging from 100 to 300 kg m −3 , corresponding to a Young's modulus E of the slab between 0.8 and 16 MPa [Eq.( 3)] and a tensile strength σ t between 1 and 16 kPa [Eq.( 4)].
In order to compare our numerical model to PST field data, we selected two simulation cases to show the overall trend of the propagation speed and distance with density, rather than simulating precisely each of the PSTs individually (which are prone to some variability) by using the available initial conditions from the field.
In the following, we distinguish two simulation cases: -Case #1 corresponds to simulations with a constant slab thickness D = 20 cm, slope angle ψ = 0 • and WL properties (σ WL c = 750 Pa); -Case #2 corresponds to a case with a slope angle ψ = 23 • which is the average slope angle of our field PSTs and a slab thickness D which is also a function of density according to field data (Table 4).In addition, we calibrated the strength of the WL bonds in order to have the same critical length for the different densities.This ensured we observed crack propagation and avoided the global and simultaneous failure of the entire WL.Indeed, as density increases, the critical length a c decreases and tends to zero (Fig. 11d) leading to the instability of the system without cutting the WL.
A similar choice was made by Gaume et al. (2015) who computed the tensile failure probability for two cases (constant depth or constant load) with similar trends in the results.

Displacement of the slab
Our numerical results (Fig. 8) obtained for a slab density ρ = 250 kg m −3 are in very good agreement with experimental results (Fig. 3) obtained for a similar density of ρ = 240 kg m −3 .Indeed, the same phases in the displacement curves, corresponding to slab bending and crack propagation, were observed in the measurements.Furthermore, the amount of slope normal displacement prior to crack propagation as well as the fracture time, defined as the time it takes for the slab to come into contact with the broken weak layer, were very similar.Finally, we would like to point out that the total slope normal displacement after crack propagation in our experimental results was not the same for all markers (Fig. 3), which has often been observed in previous studies (van Herwijnen et al., 2010;Bair et al., 2014), whereas it is approximately the same in the numerical simulations if no fracture arrest occurred (Fig. 8).This difference is presumably due to 3-D and edge effects such as wall friction at the right side of the column.

Crack propagation speed
The crack propagation speed c obtained in field PSTs and from the simulated PSTs is represented as a function of slab density in Fig. 12. Overall, the propagation speed obtained from field PSTs increases from 10 to 50 m s −1 as the density of the slab increased from 140 to 300 kg m −3 .The gray squares represent the cases with fracture arrest due to tensile fracture of the slab (SF) for which the crack propagation speed is not very accurate and generally lower than the velocity measured when the slab did not fracture (END: open squares for full propagation).
Overall, both simulation cases #1 and #2 reproduce the magnitude of the propagation speed c and the increasing trend with increasing slab density ρ.The case #2 model (relation between slab density, Young's modulus, thickness and slope angle) slightly overestimates the average propagation speed but provides good estimates for densities higher than 250 kg m −3 .
Furthermore, the simulations of case #2 were done for the same conditions of failure initiation, i.e., the strength of the WL bonds was calibrated in order to have the same critical length for the different densities.However, for the experiments, the critical length generally increases with increasing density due to the settlement which induces an increase of Young's modulus and a strengthening of the WL (Zeidler and Jamieson, 2006a, b;Szabo and Schneebeli, 2007;Podolskiy et al., 2014).In contrast, for case #1, a decrease in slab thickness and slope angle induces a decrease in the crack propagation speed (Fig. 10), explaining why the speeds for case #1 (ψ = 0 • and D = 20 cm) are lower.In addition, for case #1, the WL properties were kept constant, which together with the increase of the Young's modulus (less bending) with density resulted in an increase of the critical length with density.This is in agreement with field observations which might explain the better quantitative agreement with the experiments.
Finally, for a low slab density ρ = 100 kg m −3 (E = 0.83 MPa), the speed of the elastic waves in the slab c e is about 90 m s −1 , whereas the crack propagation speed is around 15 m s −1 .For a high density ρ = 300 kg m −3 (E = 16 MPa), c e ≈ 230 m s −1 , whereas the crack propagation speed is around 45 m s −1 .

Propagation distance
The proportion between the number of experiments for which fracture arrest was observed N SF and the total number of experiments N SF + N END decreases with increasing slab density ρ (Fig. 13a).This figure highlights the important decrease of fracture arrest propensity with slab density.For slab densities higher than 300 kg m −3 the number of experiments with slab fracture is very small (N SF < 20 %).
The crack propagation distance l * is represented as a function of slab density in Fig. 13b (SF) are represented.Overall, the propagation distance obtained from field PSTs increased with slab density and varies approximately from 0.4 to 2.1 m as the density increased from 140 to 300 kg m −3 .This trend is well reproduced by the discrete element simulations for both cases, but with a better qualitative agreement for case #1.For case #1 and for densities higher than about 300 kg m −3 , no fracture arrest is observed resulting in full propagation of the crack in the WL over the entire column length (END).For case #2, this transition occurs already for a density of about 200 kg m −3 .Besides, we would like to point out that field PSTs were not made systematically for the same column lengths.However, we checked numerically that, as soon as fracture arrest occurs within the column, the crack propagation distance is almost independent of the column length.For instance, if a propagation distance of 0.7 m is observed for a column length of 1.5 m, it will be the same for a column of 2.5 m.This is true as soon as the column length is higher than the length over which edge effect are observed (length typically lower than 1 m, Gaume et al., 2013).
The experiments and the simulations confirm that dense and hard snow slabs are more prone to wide-spread crack propagation than soft slabs.

Mechanical processes of fracture arrest
In order to better understand the underlying mechanical processes of fracture arrest in the slab, the normal stresses in the slab σ xx have to be compared with its tensile strength σ t .The normal stresses in the slab were computed for each grain from contact forces using the classic Love homogenization formula (Cambou and Jean, 2001).A tensile crack in the slab occurs when the maximum normal stress σ m xx exceeds the tensile strength.Hence, we analyzed the evolution of normal stresses in the slab during the process of crack propagation for case #1 with a slab density ρ = 250 kg m −3 leading to a tensile strength σ t = 10 kPa [Eq.( 4)].First, before the onset of crack propagation, an increase of tensile stress occurs at the top of the slab close to the crack tip of the WL (Fig. 14).The bottom of the slab is subjected to an increase in compression (σ xx < 0).This increase of tensile stress is due to the bending of the slab and increases with increasing crack length (Timoshenko and Goodier, 1970).
Then, once the critical length is reached, the crack becomes self-propagating.The crack length increase leads to an increase of the tensile stresses in the slab.Note that the maximum tensile stress σ m xx is always located at the top surface of the slab, not directly at the vertical of the crack tip but slightly shifted to the right above the undamaged weak layer (Fig. 14).At one point the maximum tensile stress meets the tensile strength of the slab (σ t = 10 kPa) which leads to the opening of a tensile crack and fracture arrest.This fracture arrest leads to the unloading of the slab where the stresses become close to zero everywhere, except at the position of the saw where some small local bending effects still occur.
In order to better understand why fracture arrest does not occur anymore for high densities, as shown in Fig. 13, we then analyzed the maximum tensile stress σ m xx as a function of slab density in the case of a purely elastic slab for which the Young's modulus was varied according to Eq. (3; case #1).The DE results were then compared to those predicted by the static beam theory.According to beam theory (Timoshenko and Goodier, 1970), the maximum theoretical tensile stress in a beam of length l, thickness D embedded on its right side and subjected to gravity, with an angle ψ (with regards to the horizontal) is equal to As we have seen before in Fig. 14, the tensile stress in the slab increases during propagation due to an increase of the crack length.However, this length is limited by the height of collapse h c of the WL (Fig. 15a).Indeed, once the left part of the slab comes into contact with the broken WL, the tensile stress reaches a maximum value and does not increase anymore (see also Fig. 7c).The length l 0 (already introduced by Heierli et al., 2008) required to come into contact with the broken WL can be obtained analytically by computing the vertical displacement of the slab according to beam theory and is equal to Hence, the overall maximum theoretical tensile stress σ m, th xx = σ th xx (l 0 ) is found by replacing l by l 0 in Eq. ( 7).For low values of slab density (ρ < 180 kg m −3 , Fig. 15b), the maximum tensile stress increases with increasing density and the DE model results are very well reproduced by beam theory.In addition, the maximum tensile stress is always higher than the tensile strength of the slab, leading to systematic fracture arrest for these low-density values.However, for higher densities, beam theory predicts a strong increase of the maximum tensile stress with density, so that the stress would always be higher than the tensile strength of the slab.This would lead to systematic fracture arrest for any value of the density which is in contrast to the results of both field and numerical PSTs.The DE model results, in particular, show that for a slab density higher than about a threshold density ρ = 180 kg m −3 , the maximum tensile stress starts to decrease with increasing density.Ultimately, for a density of approximately 280 kg m −3 , the maximum tensile stress becomes lower than the strength, leading to full propagation of the crack in the WL, in agreement with Fig. 13 (case #1).
This result highlights the limits of the static beam theory and thus the need to take into account dynamic effects when addressing fracture arrest propensity issues.Indeed, we suppose that the reason of this sudden decrease is due to the crack propagation speed which becomes higher as slab density increases and induces a loss of support in the slab where stresses do not have time to establish.In other words, the displacement of the slab due to gravity is too slow to establish a mechanical equilibrium between bending and gravity.For instance, if we assume that the crack would propagate at an infinite speed, then the tensile stresses in the slab would not increase after reaching the critical length.The maximum tensile stress in the slab would thus be the one obtained at the moment of the onset of crack propagation.Obviously, the propagation speed is not infinite but limits the establishment of the stresses in the slab.
Using the theoretical relationships for ρ < ρ The theoretical propagation distance l * th was represented in Fig. 15c for both zones (l * th = l * bt for ρ < ρ and l * th = l * dyn ρ > ρ ) as well as the characteristic distance l 0 .Again the beam theory reproduces the results for low densities well.For these low densities, the tensile failure in the slab occurs even before the onset of crack propagation due to the low value of the associated tensile strength.However, the important increase of the propagation distance for densities higher than ρ is not reproduced by beam theory.On the other hand, using the empirical relation (exponential fit of the maximum tensile stress in Fig. 15b), the strong increase of the propagation distance is well reproduced.After a certain value of the density ρ END = 280 kg m −3 , the propagation distance l * dyn becomes higher than l 0 which is technically not possible since the maximum tensile stress is obtained exactly for l = l 0 and cannot increase above l 0 .Hence the only solution is that no fracture arrest occurs for ρ > ρ END .In fact, for a density of 300 kg m −3 no fracture arrest was observed using the DE model (Fig. 13).Simulations were also performed with longer column lengths L up to 10 m which did not affect the full propagation.The corresponding maximum propagation distance for this case is about 2. agreement with field data for which the maximum propagation distance recorded was l * = 2.15 m (Fig. 13b).Obviously, the density ρ END which was 280 kg m −3 in our simulations, will vary depending on the geometry and material properties of the snowpack.For the cases presented in Fig. 15b and c, the Young's modulus was derived from density [Eq.( 3)], the slab thickness D was constant equal to 20 cm, the slope angle ψ = 0 • (case #1).This set of parameters resulted in a density ρ END = 280 kg m −3 .However, for a slope angle of 23 • and taking into account the dependence of slab thickness with slab density (Table 4), ρ END would be even less than 200 kg m −3 , as shown in Fig. 13b (case #2) since the transition between a regime of fracture arrest and full propagation is between 150 and 200 kg m −3 .Furthermore, as the propagation distance is also strongly influenced by the WL thickness D wl (Fig. 11b), we assume that ρ END increases with D wl as the maximum tensile stress in the slab increases with D wl .

Discussion
In this study, a numerical model based on the discrete element method was developed in order to perform numerical PST simulations and study the mechanical processes involved.Despite the apparent simplicity of the proposed DE model and of the structure of the simulated WL, we were able to quantitatively address the issue of the dynamic phase of crack propagation as well as fracture arrest propensity and to reproduce PST field data.
First, a parametric analysis was conducted to study the influence of snowpack properties on crack propagation speed and distance.It was shown that the propagation speed increases with increasing slab density ρ, slab Young's modulus E, slab thickness D and slope angle ψ.The propagation speed was almost not influenced by WL thickness.The increase of crack propagation speed with slab density is not compatible with the expression for crack propagation speed proposed by Heierli (2008) for which the speed decreases with increasing slab density ρ (for a constant value of the Young's modulus of the slab), as for a crack in a homogeneous material (Auld, 1973).However, this is obviously not the case here, since the crack propagates through the underlying WL.Therefore, the propagation speed is likely to decrease with increasing WL density (and thus to increase with increasing WL porosity) but to increase with slab density, as shown by our results.
In addition, it was shown that the tensile fracture of the slab always starts from the top surface of the slab.The propagation distance l * increases with increasing tensile strength of the slab σ t , slab thickness D and slope angle ψ.The latter result suggests that slab bending is the primary process influencing the tensile failure of the slab which corroborates the conclusions of van Herwijnen et al. (2010).In contrast, the propagation distance decreases with increasing slab Young's modulus E, slab density ρ and WL thickness D wl .These results are in agreement with Gaume et al. (2015)  that the tensile failure probability (fracture arrest propensity) decreased with increasing tensile strength σ t of the slab, increased with increasing Young's modulus E of the slab and decreased with increasing slab thickness D. This latter model is based on the finite element method and takes into account the weak layer strength heterogeneity, stress redistributions by elasticity of the slab, the possible tensile failure of the slab as well as dynamic effects.However, the weak layer was modeled in their approach as an interface and the bending of the slab induced by the collapse of the weak layer was not accounted for.Furthermore, by accounting for the relation between the mechanical properties of the snowpack, the increase of crack propagation speed and distance with increasing slab density was well reproduced.The slight overestimation of the propagation speed for low densities might be due to the fact that, to compute the propagation speed, the slab was considered as purely elastic and possible plastic effects in the slab that might induce energy dissipation were disregarded.
The in-depth analysis of the mechanical processes involved in fracture arrest showed that after a certain slab density value ρ , the evolution of the maximum tensile stress in the slab with slab density diverged from the static beam theory.This is due to dynamic effects during crack propagation that induce a loss of support of the slab where the stresses do not have time to establish.Ultimately, for a density ρ END , the maximum tensile stress in the slab decreases below its tensile strength leading to full propagation without fracture arrest.Consequently, for large densities, mechanical properties of the snowpack only marginally affect crack propagation distance.In that case, terrain characteristics and snowpack spatial variability will play a crucial role in the definition of the release area.
In addition, interestingly, in very few simulations both fracture arrest by tensile failure of the slab and full propagation was observed.In these cases, a portion of the WL on the right-side of the slab tensile crack was damaged over a sufficient length to exceed the critical length leading again to crack propagation.This process repeated itself until the end of the system leading to so-called "en-echelon" fractures (van Herwijnen, 2005;van Herwijnen et al., 2010;Gauthier and Jamieson, 2010).This is likely to happen for very unstable conditions (very low critical length) but for a slab of intermediate density, not too dense so fracture arrest can occur and not too loose so that the crack can still propagate.
Concerning the limitations of the model, we recall that the triangular shape of the WL structure is highly idealized and that more complex and more realistic geometries might have an influence on the presented results.In the future, the micro-structure of the WL could be derived from microtomographic images (Hagenmuller et al., 2014) in order to perform more realistic simulations.In addition, whereas we applied our model for cases for which the bending of the slab is important, our approach could still be used for cases with thinner weak layers and thus much lower height of bending.
Moreover, we would like to recall that the crack propagation speed was computed from the vertical displacement wave of the slab.However, for high values of the slope angle ψ, the collapse only constitutes a secondary process and the tangential displacement during propagation becomes higher than the vertical displacement.Typically, for ψ > 40 • , it is not possible to compute the propagation speed using the presented approach as the height of collapse becomes too low.An analysis of the tangential displacement revealed that the crack propagation speed on slopes where the shear component of the slab weight is very important (ψ > 40 • ) might be significantly higher than the propagation speed on gentle slopes.This analysis suggested propagation speeds up to 150 m s −1 , similar to those reported for real-scale avalanches by Hamre et al. (2014).However, they considered avalanches triggered artificially by explosives leading to even more complex interactions due to the propagation wave of the blast.
With regards to practical applications, the results of our study can help to choose the size of the column length in field PSTs.Indeed, we showed that the maximum length for which snowpack properties might affect the propagation distance is around 2 m, in agreement with the study of Bair et al. (2014).However, this result does not mean that all PSTs should be 2 m long.The chosen column length can be evaluated from slab thickness and density.As shown in Figs.11 and 13, slab Young's modulus and tensile strength which are related to slab density, as well as slab thickness strongly affect the propagation distance.Hence, for soft and relatively thin slabs, the standard column length of 1.2 m might be sufficient.However, for very strong and thick slabs, the column length should not be lower than 2 m in order to be able to still observe a possible arrest of the fracture due to slab tensile failure.If slab fracture is not observed in a PST for a column length of 2 m, fracture arrest is likely to be mainly driven by terrain and snowpack spatial variability and a 3Dterrain model including the snowpack might be required to evaluate where fracture arrest might occur (Veitinger et al., 2014).

Conclusions
We proposed a new approach to characterize the dynamic phase of crack propagation in weak snowpack layers as well as fracture arrest propensity by means of numerical PST simulations based on the discrete element method with elasticbrittle bonded grains.
This model allowed us to compute the crack propagation speed from slab vertical displacement as a function of snowpack properties.Furthermore, crack propagation distance was computed by taking into account the tensile strength of the slab.A parametric analysis provided the crack propagation speed and distance as a function of the different snowpack properties.We showed that the propagation speed increases with increasing Young's modulus of the slab, slab The Cryosphere, 9, [1915][1916][1917][1918][1919][1920][1921][1922][1923][1924][1925][1926][1927][1928][1929][1930][1931][1932]2015 www.the-cryosphere.net/9/1915/2015/depth, slab density and slope angle but decreases with increasing weak layer strength.The propagation distance decreases with increasing Young's modulus of the slab, slab density and weak layer thickness but increases with increasing slab tensile strength, slab depth, weak layer strength and slope angle.
Then, the existing relationship between slab thickness, Young's modulus and tensile strength with density was implemented.Accounting for this relationship, modeled propagation speed and distances were found in good agreement with those obtained from field measurements with the propagation saw test.In particular, for densities ranging from 100 to 300 kg m −3 , the propagation speed increased from approx.10 to 50 m s −1 and the propagation distance was found to increase from approximately 0.4 to 2 m (column length).Concerning the mechanical processes, the static beam theory predicts an increase of the maximum tensile stress with increasing density.However, we show that dynamic effects of crack propagation induce a loss of support of the slab which increases with increasing crack propagation speed and thus slab density.This produces a decrease of the maximum stress with density which ultimately becomes lower than the tensile strength of the slab for a critical density ρ END leading to the absence of slab tensile fracture and thus wide-spread crack propagation.According to our simulations, this critical density depends mostly on slab and WL thicknesses and slope angle.It decreases with slab depth and slope angle but increases with WL thickness.
For slab layers denser than ρ END , the slab tensile fracture in the field and thus the potential release area will mostly be controlled by topographical and morphological features of the path such as ridges, rocks, trees, terrain breaks, etc. but also by the spatial heterogeneity of the snow cover.In addition, we showed that the maximum propagation distance associated with the density ρ END was around 2 m, justifying why the column length of a propagation saw test should not be lower than 2 m for hard snow slabs, in order to be able to observe fracture arrest.This result is in agreement with the recent study of Bair et al. (2014) about PST edge effects.
In the future, an in-depth analysis of crack propagation speeds for large slope angles will be carried out in order to distinguish the speed associated with the collapse wave of the slab and the speed associated to its tangential displacement.Finally, different and more complex structures for the WL will also be implemented with the long-term objective to model the structure of the WL directly from segmented micro-tomographic images (Hagenmuller et al., 2013(Hagenmuller et al., , 2015)).

Figure 1 .
Figure 1.(a) Typical slab -weak layer configuration suitable for avalanche release.The weak layer is composed of surface hoar which is intact on the right and partially ruptured on the left.©ASARC fromJamieson and Schweizer (2000).(b) Zoom on a surface hoar crystal © ASARC.

Figure 2 .
Figure 2. Schematic drawing and picture of the propagation saw test (PST).The black dots are markers used for particle tracking in order to measure the displacement of the slab.The column length is denoted L. Adapted from van Herwijnen et al. (2010).

Figure 3 .
Figure 3. Temporal evolution of the measured vertical displacement y for a slab density of ρ = 240 kg m −3 obtained through PTV analysis of the marker's displacement from a field experiment.The different curves correspond to different horizontal positions in the slab, from the left-end (x = 0.1 m) to the right-end (x = 2.1 m).

Figure 4 .
Figure 4. (a) Simulated system of the propagation saw test (PST) composed of a slab, a weak layer and a rigid (fixed) substratum.The column is 2 m long.(b) Zoom on the weak layer structure.The blue lines represent the weak layer bonds.
Figure 5. (a) Schematic representation of the parallel bond contact model which is used.The bonded part is represented in black while the unbonded part is represented in grey.(b) Bond normal force F n as a function of normal interpenetration δ n between two grains.(c) Bond shear force ||F s || as a function of tangential interpenetration δ s between two grains.(d) Bond bending moment ||M|| as a function of bending rotation θ between two grains.

Figure 6 .
Figure 6.Failure envelope of the modeled WL obtained from mixed-mode shear-compression loading tests.The angle represented next to the data points is the slope angle, σ WL t and σ WL c

Figure 7 .
Figure 7. Snapshots of a PST numerical experiment.(a) Initial system t = 0.1 s, (b) onset of crack propagation t = 0.26 s, (c) dynamic propagation t = 0.28 s; (d) complete failure of the WL t = 0.45 s.

Figure 8 .
Figure 8. Temporal evolution of the modeled vertical displacement y of the slab for a slab density ρ = 250 kg m −3 .The different curves correspond to different horizontal positions of the slab, from the left-end (x = 0 m) to the right-end (x = 2 m).

Figure 9 .
Figure 9. Snapshot of a PST with fracture arrest due to tensile crack opening in the slab induced by slab bending.

Figure 10 .
Figure 10.Crack propagation speed c (continuous lines) and speed of elastic waves in the slab c e /10 (dotted lines) as a function of (a) the Young's modulus of the slab E, (b) slab thickness D and WL thickness D wl , (c) slab density ρ, (d) slope angle ψ and (e) WL compressive strength σ WL c .The parameters used for these figures are given in Table3.

Figure 12 .
Figure 12.(a) Boxplot of the propagation speed c versus slab density for all field experiments.The data were grouped by slab density classes of 50 kg m −3 .The red line represents the median value, the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points without considering the outliers, and outliers are plotted individually as a cross.(b) Crack propagation speed c vs. slab density ρ.The open squares correspond to field PSTs with full propagation (END) and the solid squares correspond to PSTs with fracture arrest (SF).The black line corresponds to the result of the DE model taking into account the relation between slab density, Young's modulus for a slope angle ψ = 0 • and a slab thickness D = 20 cm (case #1).The gray continuous line corresponds to the result of the DE model taking into account the relation between slab density, Young's modulus and thickness for a slope angle ψ = 23 • (case #2).The red lines represent the median values of the density classes (same as in a).The data consist of N = 121 PST experiments.
Figure 13.(a) Proportion between the number of experiments with slab fracture N SF and the total number of experiments N SF +N END for different classes of density.(b) Propagation distance l * vs. slab density ρ only for cases with fracture arrest (SF).The red line corresponds to case #1 and the gray line corresponds to case #2.The dashed line corresponds to the column length L = 2 m in the PST simulations.

Figure 14 .
Figure 14.Evolution of the normal stress σ xx in the slab during the process of crack propagation for a case #1 simulation and a density of 250 kg m −3 with a tensile strength of 10 kPa.The normal stress σ xx was calculated for each grain from contact forces according the Love homogenization formula and linearly interpolated between grains.The dashed lines represent the position of the crack tip in the WL.
(σ m xx = σ th xx ) and an empirical (exponential) fit to the data for ρ > ρ (σ m xx = σ dyn xx ), one can compute the theoretical propagation distance l * th by solving Figure 15.(a) Illustration of the characteristic length l 0 required for the slab to be in contact with the broken WL. h c represents the collapse height.(b) Maximum tensile stress σ m xx as a function of slab density for the discrete element model σ DEM xx , with a purely elastic slab and case #1, according to beam theory σ th xx (l 0 ) and from an empirical (exponential) fit σ dyn xx (l 0 ).The tensile strength of the slab σ t is also represented as a function of slab density ρ.(c) Propagation distance for the same cases as in (b) plus l 0 as a function of slab density.

Table 1 .
Mechanical parameters used in the simulations for the contact law.k n : normal contact stiffness; k s : tangential contact stiffness; µ: intergranular friction; e: normal restitution coefficient.

Table 2 .
Mechanical parameters used in the simulations for the cohesive law.k n : bond normal stiffness; k s : bond shear stiffness.σ t : macroscopic tensile strength; σ s : macroscopic shear strength.

Table 3 .
Table of the parameter values used for Figs. 10 and 11.The symbol "-" means that the associated parameter was varied.The symbol * means that the WL failure envelope was calibrated to obtain a constant critical length a c = 15 cm.

Table 4 .
Average slab thickness as a function of slab density for PST field data.