Modelling snow failure with a fibre bundle model

Abstract

Dry-snow slab avalanches initiate from a failure in a weak snow layer below a cohesive slab. Snow is considered as a porous ice structure, and the strength distribution of the single elements of this structure, i.e. grains and bonds between grains, shows a high degree of disorder. On the bond or microstructural level, the failure process is believed to start if the fracturing of bonds between snow grains is not balanced by the formation of new bonds. We use a statistical fracture model – a fibre bundle model – to study the failure process in a weak snow layer. The model consists of fibres of various strengths representing single snow grains between two rigid plates which represent the slab above and the substratum below the weak layer. The fibres deform in a linear elastic manner and break instantly at their rupture strength. Broken fibres may sinter (re-bond) and regain strength after a finite sintering time. We show that the different characteristic times for breaking and sintering lead to the rate dependence of snow strength. This is, to our knowledge, the first statistical model to reproduce the ductile-to-brittle transition which snow exhibits with increasing strain rate. When the model is applied to simulate experimental stress–strain curves for different strain rates, the model and experimental results are in fair agreement.

Introduction

Dry-snow slab avalanches start with a failure in a weak snow layer below a cohesive slab (McClung, 1979; Schweizer and others, 2003). Fracture mechanical models of dry-snow slab avalanche release postulate an initial crack in the weak layer and study the circumstances under which this pre-existing crack will propagate in terms of crack energy (e.g. McClung, 1981; Schweizer, 1999; Heierli and others, 2008). However, how this initial crack originates cannot be determined within the theory of linear elastic fracture mechanics. In other words, the failure initiation problem related to spontaneous avalanche release is essentially unsolved.

We aim to model the mechanical behaviour, on the basis of heterogeneities within the weak layer, with a statistical failure model, namely a fibre bundle model. We are particularly interested in the ductile-to-brittle transition (DBT) which snow exhibits with increasing strain rate, and the formation of the initial crack (which may lead to crack growth and avalanche release). In the following, we use the terms ductile and brittle as described in Narita (1983) and commonly used in ice mechanics (Petrenko and Whitworth, 1999). Brittle fracture behaviour means that virtually no, or very little, plastic deformation occurs before fracture. Ductile behaviour implies large irreversible deformation before failure or no failure at all. Subsequently, we briefly introduce a few elements of the mechanical behaviour of snow such as the DBT, microstructure, sintering and disorder.

Field measurements of weak layer strength show large scatter in results (e.g. Jamieson and Johnston, 2001) because the mechanical properties of snow depend on microstructural characteristics which are variable both in space and time. This impairs the reproducibility of experiments performed at different times even at the same snowpack location. Therefore, various studies on shear failure of homogeneous snow samples under controlled laboratory conditions have been performed (e.g. Schweizer, 1998). Depending on snow type and temperature, these studies showed brittle behaviour of snow at deformation rates typically faster than 10⁻³ s⁻¹ and ductile behaviour at smaller deformation rates. Up to now, experiments on layered snow samples have only been performed by Fukuzawa and Narita (1993). For their experimental set-up they found the DBTat a strain rate of about 10⁻⁴ s⁻¹, an order of magnitude smaller than what Schweizer (1998) found for homogeneous snow samples under shear loading. Scapozza and others (2004) reported that the acoustic emission response of snow under compression shows a clear strain-rate dependency which is related to the DBT. Narita (1983) performed tensile experiments with homogeneous samples, and observed brittle fracture above strain rates of the order of 10⁻⁴ s⁻¹. The DBT increased with increasing temperature and with increasing snow density.

Density, which is an easily measured macroscopic parameter, is often used to characterize the bulk properties of snow. However, it is not a good predictor of mechanical strength (Shapiro and others, 1997): two snow samples with the same density but different types of microstructure may have strengths differing by as much as a factor of four (Keeler and Weeks, 1968). Bartelt and von Moos (2000) reported that snow properties change for a given density by as much as an order of magnitude due to differences in microstructure. It is therefore important to consider the microstructure of snow, especially when studying small-scale processes like failure initiation. The term microstructure refers to objects or geometries with characteristic length scales of about O(10⁻⁴ m) (form, arrangement, and size of grains and bonds), while macroscopic objects are of length scales of O(10⁻¹ m) and larger (the initial crack, snow samples used in the laboratory, avalanches).

During snow deformation (density: ρ ≤ 300 kg m⁻³), the microstructure (and hence the bulk mechanical properties) constantly changes because of rearrangement of single snow grains (Camponovo and Schweizer, 2001). In other words, continuous bond breaking must occur. Since snow in a natural snow cover exists close to its melting point, two snow grains which come into contact may easily bond to each other. The failure process preceding avalanche release is therefore believed to be related to two fundamental, but competing, processes at the microscale: bond fracturing and sintering (bond formation) (Schweizer, 1999).

Three-dimensional (3-D) images of snow microstructure (e.g. Schneebeli, 2004) clearly show that snow is a highly disordered material consisting of an ice matrix and open pores filled with air and water vapour. The degree of heterogeneity is supposed to affect failure initiation at various scales (e.g. Schweizer and others, 2008). Continuum (macroscopic) quantities like global stress, global strain and global strength are not evenly distributed over all micro-structural elements of the snow structure. We assume that failure starts where the weakest points meet the highest stress within the structure. Accordingly, the strength distribution is expected to play a major role in the failure process.

To model snow deformation and failure, both the microstructure and the disorder of snow need to be considered. There are several mechanical snow models that have included the microstructure in some way. These can be divided into three groups:

1. 1. Models that are essentially continuum models but include some parameterization of microstructure (e.g. Mahajan and Brown, 1993). This approach has been used to describe snow viscosity (Bartelt and von Moos, 2000) and to simulate the mass and energy balance of the snow cover (Lehning and others, 2002). The model by Gibson and Ashby (1997), where snow is considered as an open-cell foam of ice with two structural parameters (beam length and beam cross-section), also belongs to this group. Louchet (2001) followed their approach and considered healing using rate equations for bond rupture and bond formation.

2. 2. Models that try to reproduce the microstructure in simplified (or generalized) form, for example as an arrangement of beams (or spheres), possibly with some random variations of the local properties from one element to the next (Herrmann and Roux, 1990). Johnson (1998) used a dynamic finite-element computer program to study the rapid compaction of snow represented as an arrangement of randomly distributed spheres.

3. 3. With new experimental techniques, such as X-ray microtomography (μCT), the full 3-D representation of microstructure is used as input for a real microstructural (or specimen-specific) finite-element model (Schneebeli, 2004).

An approach that considers disorder, although on a macroscopic scale, was used by Fyffe and Zaiser (2004) who introduced time-dependent strength recovery and randomly varying shear strength into a model of snow slope failure.

Fibre bundle models (FBMs) (Daniels, 1945; Alava and others, 2006; Raischel, 2007) are statistical fracture models that can include a simple representation of the microstructure of porous, granular materials (see Kun and others, 2006, for a detailed description). They are especially helpful in describing materials with time-dependent failure effects (e.g. the fatigue failure of asphalt (Kun and others, 2007)). As such, the FBM technique offers the possibility of simulating the important ductile-to-brittle failure transition of snow, as well as competing microscale processes such as bond fracturing and sintering (Schweizer, 1999). FBMs belong to the second class of microstructural models mentioned above.

In this paper, we apply a fibre bundle to simulate snow deformation, damage, and failure of a weak snow layer. Our aim is to investigate the influence of microstructural parameters on the bulk mechanical response of weak snow layers under shear loading. In particular, we test the hypothesis that different timescales of fracturing and sintering of bonds explain the strain-rate dependence of the failure behaviour of snow (Schweizer, 1999). We use experimental results of snow under shear deformation (Schweizer, 1998) for comparison with our model output.

Methods

On an inclined slope, the gravitational force induces shear and normal deformation in the snowpack. Experimental studies and observations of natural buried weak layers show that the shear deformation is concentrated in the weak layer (Fukuzawa and Narita, 1993; Jamieson and Schweizer 2000). In our approach, the snow crystals (e.g. buried surface hoar) that form the weak layer correspond to fibres (Fig. 1). The fibres are located between two rigid plates which represent the slab above and the substratum below the weak layer. Since the layers above and below are substantially stronger and about an order of magnitude stiffer than the weak layer (Fukuzawa and Narita, 1993; Jamieson and Schweizer, 2000), the simple assumption of rigid plates is justified as a first approximation.

Fig. 1. Photo of a buried surface hoar layer, fractured on the left, intact on the right. Overlying the photo is a schematic drawing of the fibre bundle. Photo adapted from Jamieson and Schweizer (2000).

Each fibre i, where i = 1,…N, behaves in a perfectly elastic manner. In order to statistically model the spatial variability of strength which is presumed to be caused by variations in crystal orientation, size and bonding within the weak layer, the initial strengths σ _c, _i are taken from a Weibull probability distribution, where the density function is given by

(1)

with scale factor α and shape factor β. This is the most commonly used probability distribution in statistical fracture models (Herrmann and Roux, 1990; Chakrabarti and Benguigui, 1997) and has previously been used to describe the strength of snow (e.g. Sommerfeld, 1973; Kirchner and others, 2004). The effect of microstructure is incorporated into the model through this strength distribution.

Following the deformation-controlled shear experiments mentioned above, at each discrete time-step the upper plate is moved along the x axis by a constant amount Δx (Fig. 2). The constant global shear strain rate is therefore given by

(2)

where l ₀ is the initial fibre length (Fig. 2). As a first step we assume no vertical (z-direction) displacement of the upper plate. Furthermore, we treat the fibres as long, thin truss elements under uniaxial tension, i.e. shear deformation within a fibre is neglected.

Fig. 2. Schematic representation of model geometry.

Since the deformation of the (stiff) upper plate is imposed, the load on a single fibre is always given by the external deformation, no matter how many fibres are intact. The fibres break in order of increasing strength as determined by the Weibull distribution applied to each fibre, i.e. the weakest fibres break first. The elongation Δl of a fibre can be calculated via

(3)

where m is the number of time-steps the fibres have been intact. The force f_j (t) to deform a single fibre is

(4)

where σ_j (t) is the stress the fibre experiences, E is the Young’s modulus, a is the cross-section of the fibre, j = 1, …N _intact with N _intact the number of intact fibres, and ε_j (t) denotes the axial strain on the fibre.

If σ_j (t) acting on a fibre reaches its rupture strength, i.e. σ_j (t) = σ _c,j , the fibre breaks instantaneously. At the next discrete time-step the strength of the fibre σ _c,j is zero, and the fibre is considered broken. If more than half the fibres are broken, we consider the whole bundle as fractured. However, at each time-step there is the probability p _s that the broken fibre can start sintering (re-bonding). The sintering probability is proportional to the square of the number of broken fibres, p_s = p _max(N _broken/N)² since two fibre ends are necessary for sintering (second-order kinetics). We use periodic boundary conditions, i.e. the total number of fibres (broken and intact) remains constant and each broken fibre end always has an appropriate partner for sintering, no matter how far the upper plate has already moved. This seems realistic, since a weak snow layer’s height (in nature and also in the laboratory) is very small compared to its width and length. If a fibre is chosen to sinter at t_n , it is considered intact at the next time-step t_n +₁ = t_n + At and will then experience deformation and therefore also stress. Its final strength will only be reached after s time-steps, i.e. at t_n + t_s . The time evolution of fibre strength during sintering is given by

(5)

The new final fibre strength after sintering (t>t_s ), σ′_{c, i , final′} is again taken from the same probability distribution as the initial σ _c,i . So after sintering, the fibre strength can be smaller or larger than it was initially while the Young’s modulus stays the same.

The global force F(t) needed to perform the change of position of the upper plate (dynamics in the sense of accelerations are not considered) is calculated via

(6)

and the global stress σ _global(t) is then obtained via

(7)

where A is the area of the upper/lower plate. The FBM is governed by only three model parameters: the shape of the Weibull distribution, β, the maximum sintering probability, p _max, and the time it takes for a fibre to regain its full strength, t_s .

Model Results and Discussion

We first compare our model to the analytical solution which exists in the absence of sintering. The constitutive relation for a deformation-controlled FBM under tension is given by

(8)

where P(σ _c = Eε) is the cumulative probability distribution of fibre strength (Kun and others, 2006). If the sintering probability is set to zero and we simulate tensile deformation, i.e. Δl _fibre = Δx _bundle, our simulations agree with the analytical solution. The fact that we have shear deformation shifts the global stress–strain curve towards higher deformations (Fig. 3), because Δl < Δx (Equation (3)). Without sintering, the strain rate has no effect on the global stress–strain relation. As we assume stiff plates and imposed external deformation, there is no load sharing. A more detailed description of deformation- vs force-controlled fibre bundles is given by Kun and others (2006).

Fig. 3. Comparison of the analytical solution (tension) with FBM simulation results without sintering (tension: circles; shear: small solid dots) (N = 1000, β = 3, α = 2).

Compared to a FBM without sintering, our model produces stress–strain curves with a larger initial slope and a higher peak stress. This is due to a shift in fibre strengths to higher values. Although the newly assigned fibre strengths after sintering are taken from the same probability distribution as the strengths assigned initially, stronger fibres survive longer. In the natural snow cover, strengthening also occurs with time due to compaction and densification under its own weight which also leads to rearrangement, i.e. breaking of bonds and sintering, of grains.

Figure 4 shows the effect of different maximum sintering probabilities p _max and sintering times t_s on the mechanical behaviour for constant shape factor β = 1 and strain rate . Note that the model results are given in arbitrary units. At this strain rate, the bundle fractures only for low p _max. The slope of the stress–strain curve increases with increasing sintering probability. Brittle behaviour is favoured by either decreasing p _max or increasing t_s , or by increasing the strain rate (Fig. 5).

Fig. 4. FBM simulation results with varying sintering time t_s and maximum sintering probability p _max (N = 10 000, , β = 1, α = 2). The vertical arrows mark the point where the bundle fractures, while the horizontal arrows indicate that the bundle is still intact but the simulation was stopped (stress–strain curve would continue almost horizontally).

Fig. 5. FBM simulation results with varying shape factor β for the Weibull distributions of single-fibre strength: (a) β = 0.5; (b) β = 1; and (c) β = 3. For each shape factor, stress–strain curves for three different strain rates (10⁻², 7.2 × 10⁻⁵ and 5 × 10⁻⁶) are given (N = 100 000, α = 2, p _max = 0.0001, t_s = 450). The vertical arrows mark the point where the bundle fractures, while the horizontal arrows indicate that the bundle is still intact but the simulation was stopped (stress–strain curve would continue almost horizontally).

The snow microstructure is included in the model by the shape of the strength distribution. Since no experimental data yet exist, we have to assume the parameters. As Figure 5 suggests, they strongly influence the type of stress–strain curve. While the brittle fracture behaviour is not altered with increasing β, the ductile failure behaviour changes from creep (no failure) to strain softening. For ductile behaviour, the strength increases with the spread of the initial distribution. For β > 1 we find strain softening in the ductile case because the (always present) strengthening effect cannot compensate for the broken fibres which do not contribute to the global stress.

In case the gain of strength due to sintering is faster than the increasing internal stress which results from the external deformation, a sintering fibre survives, i.e. it sinters until it reaches its final strength. Since the breaking of a fibre occurs instantaneously (within one time-step) while the sintering takes s time-steps, slow deformation rates favour the sintering while the breaking dominates at fast deformation rates. This leads to a transition from ductile failure behaviour to brittle fracture with increasing strain rate. At high strain rates, the bundle breaks after little deformation. At slow deformation rates, the bundle does not break at all, as we use periodic boundary conditions. If open boundary conditions were imposed, the bundle would break when the deformation exceeded half the length, in the same way that a natural finite snow sample will break at some point.

To compare our model with experimental data, we assign physical units to the time of a fibre breaking (time-step) and the geometric dimensions of the model (Table 1). The parameters that are not directly accessible in an experiment have to be assumed. For the Young’s modulus and the tensile strength, we have taken typical values for ice (Petrenko and Whitworth, 1999). Density, height and area of the bundle are as in the experiment. We assume that 10% of the ice matrix contributes to the mechanical resistance of snow, in accordance with findings by Bartelt and von Moos (2000). The sintering time we use agrees with the range of values found by Szabo and Schneebeli (2007) in their study on sintering of ice.

Table 1. Overview of the model input parameters with typical values used

In Figure 6 we plot the experimental data from Schweizer (1998) together with our model results. Between the three modelled curves, only the strain rate is altered; all other parameters remain unchanged (in accordance with the experiment). At the slowest strain rate, the model results agree quite well with the experimental data. The higher residual stress reached in the simulation might be due to the periodic boundary conditions used for the model. For the intermediate strain rate, a higher fracture strain than in the experiment is found, but the fracture stress is roughly the same. At the highest strain rate, the model shows higher fracture strain and stress than found in the experiment. The convexity of the modelled curve is due to the geometric arrangement of the fibres. They are arranged vertically but the plates move horizontally. At small deformations, the exact orientation of the microstructural elements becomes increasingly important. This effect we cannot capture in our model containing vertical fibres only. However, the different failure behaviour which snow exhibits at different strain rates, i.e. brittle behaviour at high and ductile behaviour at low strain rates, is captured very well by our model.

Fig. 6. Comparison of FBM simulation with experimental results from Schweizer (1998).

Conclusions

We used a simple FBM to simulate the deformation and failure of a weak snow layer under shear loading. For the case of no sintering, the model can be solved analytically and our simulation results agreed with the analytic solution (Fig. 3). By incorporating sintering (i.e. re-bonding of fibres after fracture) into the model, we showed that the rate dependence of snow strength and the rate-dependent mechanical behaviour can be reproduced by introducing a sintering probability and two different characteristic times for sintering and for breaking of fibres (Fig. 5). Thus the competing effects of bond breaking and sintering of bonds (after rearranging) between snow grains might be sufficient to explain the strain-rate-dependent bulk behaviour of snow. This has so far only been hypothesized (Schweizer, 1999). Our results also suggest that incorporating sintering is crucial for realistically modelling the mechanical behaviour of snow, at least for timescales of the order of seconds or larger.

Despite the simplicity of the model, we found good qualitative agreement of the model output with the experimental data of deformation-controlled shearing of snow samples at different strain rates (Schweizer, 1998).

A more realistic modelling approach to describe the behaviour of snow will include dropping the assumption of stiff plates and introducing load sharing into the model. We could then study the influence of slab stiffness on the weak layer. Furthermore, metamorphism could be included in the model so that the strength distribution would explicitly change over time.

In the future, a more sophisticated, statistically based microstructural failure model might be coupled with a fracture mechanical model to simulate failure initiation, as well as fracture propagation in the snow slab-avalanche release process.