Comparing measurements of snow mechanical properties relevant for slab avalanche release

Abstract

Snow properties relevant to the fracture processes involved in dry-snow slab avalanche release include weak layer specific fracture energy, slab elastic modulus and density. Various techniques exist to determine these snow mechanical properties, but it is presently unclear how values determined with different methods compare. In the laboratory, the 3-D microstructure of cm-sized snow samples is reconstructed by micro-computed tomography (μCT) so that density and elastic modulus can be computed. In the field, fracture energy and modulus are estimated based on particle tracking velocimetry (PTV) of the displacement field observed during propagation saw tests. Snow stratigraphy is measured with the snow micro-penetrometer (SMP) in either, field or laboratory. We compared SMP-derived properties to corresponding μCT- and PTV-derived values. Values of snow density related well to μCT results and so were SMP-derived elastic moduli related to PTV-derived values. By taking into account snow anisotropy a good relation between SMP- and μCT-derived moduli resulted suggesting the SMP-derived modulus characterizes the components of the modulus perpendicular to the axis of penetration. SMP- and PTV-derived values of fracture energy were correlated. The SMP can provide a bridge between scales and techniques, yet further improvements in signal interpretation are still needed.

1. INTRODUCTION

Dry-snow slab avalanche release involves a sequence of fracture processes, most importantly failure initiation and crack propagation, which eventually leads to the detachment of the snow slab (e.g. Schweizer and others, 2003). Following this idea, metrics related to these essential processes were suggested (Reuter and others, 2015a; Gaume and Reuter, 2017; Gaume and others, 2017) to determine whether a snow cover is prone to fail–rather than considering a single stability criterion (e.g. Föhn, 1987). By describing snow instability by failure initiation, crack propagation, it became apparent that the properties of the slab layers and the weak layer are equally important (e.g. van Herwijnen and Jamieson, 2007; Schweizer and Reuter, 2015). Recent modelling, as well as analytical studies, have confirmed this (e.g. Gaume and others, 2015; Monti and others, 2016; Reuter and Schweizer, 2018). Failure initiation is often modelled with a strength-of-material approach (e.g. Conway and Abrahamson, 1984) with the density of the slab and strength of the weak layer as the relevant parameters. For crack propagation, the critical cut length, as it would be measured in a propagation saw test (PST) (Gauthier and Jamieson, 2006; Sigrist and Schweizer, 2007), is considered the relevant measure. The critical crack length integrates slab and weak layer properties and can be modelled using the density and the elastic properties of the slab and the specific fracture energy of the weak layer, which is the resistance to crack propagation (Sigrist and Schweizer, 2007). It is, however, well-known that measuring snow mechanical properties is intricate and different techniques commonly yield different results, with obvious implications for predicting snow instability.

Density is measured in the field with a precision of ~5 % by weighing a snow sample of a given volume (Proksch and others, 2016). From micro-computed tomography (μCT) images the ice volume fraction is determined in a non-destructive way with high-spatial resolution (1 mm) providing the most detailed and accurate measurement to date (Proksch and others, 2016). Also, Proksch and others (2015) have shown that with snow micro-penetrometry the density of snow is obtained with good accuracy comparedwith 3-D reconstruction of μCT images (Schneebeli and Sokratov, 2004). More recently, Kaur and Satyawali (2017) derived snow density from SMP signals with a statistical model based on the number of peaks of maximum penetration resistance.

The elastic modulus is probably the most delicate parameter with respect to snow instability modelling since for typical densities of seasonal snow between 100 and 400 kg m⁻³ the modulus naturally spreads over at least two orders of magnitude (Mellor, 1975) (Table 1). The Young's modulus describes the relation between stress and strain during small reversible deformations. Elastic properties of snow were determined in shear (e.g. Schweizer, 1998), torsional shear (Camponovo and Schweizer, 2001), tension (e.g. Narita, 1980) and compression (e.g. Scapozza, 2004). Values of the elastic modulus reported from mechanical testing are of the order of ~1–100 MPa for typical densities of seasonal snow and snow temperatures between −5°C and −20°C (Schweizer, 1998; Camponovo and Schweizer, 2001; Scapozza, 2004).

Table 1. Selection of previous studies reporting the elastic modulus (E) of seasonal snow by mechanical testing, finite element simulations (FEM) of 3-D micro-computed tomography (μCT), particle tracking velocimetry (PTV) of propagation saw tests (PST) or acoustic wave propagation speed (AWS)

Strain rate, number of experiments (N), average density and snow temperature are listed if reported.

Traditional quasistatic loading experiments at low strain rates may not reveal the true elastic modulus, because the observed deformation contains non-elastic contributions (Mellor, 1975). In particular experiments with insufficient measurement resolution do not provide the true Young's modulus, as the initial elastic deformation cannot be resolved. Most appropriate values for Young's modulus are obtained during high-frequency dynamic experiments as suggested by Mellor (1975). Dynamic measurements with a cyclic loading device at 100 Hz were used by Sigrist and others (2006) who eventually calculated the energy release rate in mode II. Also, with snow micro-penetrometry sufficiently high strain rates can be reached and small strains are measured so that a micromechanical modulus can be derived (Marshall and Johnson, 2009). Sigrist and others (2006) found a good correlation of the SMP-derived micromechanical modulus with their dynamic measurements. They reported a linear relationship with a regression slope of 186.

A promising approach to circumvent the issues with mechanical testing is deriving the elastic modulus from measurements of the acoustic wave propagation speed in snow columns (Capelli and others, 2016). Comparing this method with microstructure-based finite element modelling Gerling and others (2017) found good agreement between values of the elastic modulus which likely reflects the Young's modulus of snow.

Microstructure-based modelling with finite elements allows deriving an elastic modulus of snow based on the elastic properties of ice. To this end, the stresses due to prescribed displacement are modelled from 3-D reconstructed μCT images (Schneebeli, 2004; Srivastava and others, 2010). Following the approach of Schneebeli (2004), Köchle and Schneebeli (2014) analyzed 32 samples containing natural snow layers to derive parameterizations of an isotropic elastic modulus for two density ranges below and above 250 kg m⁻³. By numerical homogenization to upscale from the microstructural image, Wautier and others (2015) determined the macroscopic Young's modulus for 31 samples and obtained values of up to an order of magnitude larger than Köchle and Schneebeli (2014). The differences partly relate to assuming orthotropic (rather than isotropic) material properties, but also to mechanical modelling assumptions (e.g. boundary conditions, smaller representative elementary volume) and differences in the μCT analysis (e.g. segmentation thresholds, measuring resolution). In this context, Srivastava and others (2016) introduced three fabric measures to account for anisotropy. In particular, in cases of faceted crystals and depth hoar, the prediction of the elastic modulus could be improved.

The specific fracture energy w _f is used to describe the resistance of weak snowpack layers to crack propagation and can be calculated based on fracture mechanical experiments or snow microstructure (Table 2). Kirchner and others (2002) were the first to report the critical fracture toughness of snow in tension and shear. From shear fracture experiments of natural layered snow samples, Sigrist (2006) and Sigrist and others (2006) derived values of the specific fracture energy. For a set of 27 weak layers containing facets, depth hoar and surface hoar they reported a mean fracture energy of w _f = (0.04 ± 0.02) J m⁻². LeBaron and Miller (2014) determined the fracture energy from 3D micro-tomography reconstructions of snow samples. For a snow sample of small rounded grains, they obtained w _f = 0.06 J m⁻². They calculated the critical energy release rate G _c ∝ 2 w _f along the minimum cut surface following the idea of Hagenmuller and others (2014a), who defined an area separating a sample at minimum energy cost. Later, LeBaron and Miller (2016) presented a discrete element modelling approach to derive the specific fracture energy and compared their values for three samples with shear laboratory experiments.

Table 2. Specific fracture energy (w _f) as reported in previous studies using different experimental methods: Finite element modelling (FEM), micro-computed tomography (μCT), propagation saw test (PST), particle tracking velocimetry (PTV) or snow micro-penetrometry (SMP). including failure mode (if reported), number of experiments (N), failure layer density (if not provided: average adjacent layer density in brackets)

Values obtained from field experiments range between 0.01 and 2 J m⁻² (Sigrist and Schweizer, 2007; Schweizer and others, 2011). McClung (2015) discussed the differences in view of strain rate dependencies during experiments and presented estimates of fracture energy from the shear fracture model of Palmer and Rice (1973) based on a large field dataset yielding values between 0.01 and 0.2 J m⁻². Most studies reporting values of the specific fracture energy are based on experiments with the PST, which for modelling purposes is the best-suited snow stability test due to its geometry and the well-defined loading state. Assuming linear elastic bending of the slab due to the saw cut in the weak layer the specific fracture energy can be derived from changes in strain energy. Modelling the strain energy with finite elements, which requires assuming an elastic modulus of the slab, Sigrist and Schweizer (2007) reported w _f = (0.07 ± 0.02) J m⁻² and Schweizer and others (2011) obtained w _f = (1.3 ± 0.8) J m⁻² for different sets of weak layers. The differences could partly be explained by the choice of elastic moduli used for the FE simulations. Also Gauthier and Jamieson (2010) reported values of the specific fracture energy from a dataset of 145 PSTs described in Gauthier and Jamieson (2008). They obtained the fracture energy with an analytical expression presented by Sigrist (2006) and reported a median of w _f ≈ 0.03 J m⁻². Requiring assumptions about the elastic modulus of the layered snow slab when determining the fracture energy, being a weak layer property, is the downside of using the finite element method to determine the strain energy. This is circumvented by the PTV technique where the displacement field in a PST is analyzed (van Herwijnen and others, 2010). The elastic modulus of the slab and the specific fracture energy of the weak layer are both determined simultaneously (van Herwijnen and others, 2016), however, at rather low strain rates likely involving contributions of non-elastic deformation. To model the critical crack length as would be observed in a PST Reuter and others (2015a) obtained the specific fracture energy from SMP measurements; their values of the critical crack length had an error of a few cm when they complemented an analytical expression with FE simulations to account for slab layering.

The snow micro-penetrometer (SMP) is a device recording profile of penetration resistance at a resolution of ~250 measurements per mm at a speed of 20 mm s⁻¹ (Schneebeli and Johnson, 1998), and allows measuring with one measurement all snow layers relevant for avalanche release modelling. Whereas the penetration resistance and its signal characteristics can be related to grain type and shape (e.g. Satyawali and others, 2009), models of micro-penetration (e.g. Johnson and Schneebeli, 1999) offer an alternative approach to conclude on the mechanical parameters of snow structure that produced the signal. Interpreting the measured signal of penetration resistance as a Poisson shot noise process, force and displacement-related parameters of individual snow particles are derived (Löwe and van Herwijnen, 2012), which allow calculating micromechanical snow properties (Johnson and Schneebeli, 1999).

Apart from the SMP snow mechanical properties relevant for modelling snow avalanche release can currently also be derived in the field by analyzing the displacement field by particle tracking as recorded during a PST experiment. Apart from preliminary comparisons (Reuter and others, 2013; Schweizer and others, 2016), it is presently unclear how well the values obtained with these two methods agree. Moreover, a comparison with an alternative method is lacking. Since the elastic modulus can be determined by microstructure-based finite element modelling with good accuracy (Gerling and others, 2017), the method can be seen as a reference measurement for this snow property.

Our aim is, therefore, to compare the two methods (SMP, PST) and relate them to data obtained from μCT. As the SMP is the only technique applicable in the field as well as in the laboratory, it has the potential to bridge these scales as well as the purely laboratory- and field-based measurement methods. We hence analyze two datasets, one linking PST with SMP measurements in the field (PST–SMP dataset), and the other one linking SMP with μCT measurements in the laboratory (SMP–μCT dataset). We provide comparisons between the properties of dry snow, namely density, elastic modulus and specific fracture energy.

2. DATA AND METHODS

2.1. Data

The PST–SMP field dataset consists of 83 side-by-side field measurements that were conducted in the Swiss Alps around Davos on 19 days between 1 March 2010 and 3 March 2015. On each of these days, SMP measurements and PST experiments were performed no further than 0.3 m from each other and adjacent to a snow pit at a distance of <5 m. The manual snow profile (CAA, 2007) included manual density measurements with a 100 cm³ density cutter. The SMP measurements were conducted along a contour line of the slope and snow columns were prepared for PTV analysis of PST experiments – each one close to an SMP measurement. The PTV data are part of a larger, recently published dataset (van Herwijnen and others, 2016). The SMP data were partly used earlier by Reuter and others (2015b). Here we show the overlapping data, i.e. the cases with SMP measurements and adjacent PST experiments, which allow deriving snow properties from every measurement. The snow properties derived from this dataset are for snow layers with a thickness of ~1 up to several centimetres.

The SMP–μCT laboratory dataset includes SMP measurements and μCT measurements that were previously presented by Riche and Schneebeli (2013) and Proksch and others (2015). The dataset consists of 27 natural snow samples and nature identical snow samples produced with the ‘snow maker’ (Schleef and others, 2014). From the centre of the snow samples, one μCT sample (≤0.5 cm³) was extracted and up to four SMP measurements were taken around the location of the μCT sample allowing direct comparison. A sketch of the original measurement setup is provided by Proksch and others (2015, Fig. 2a therein). The snow structural properties derived from these data correspond to cubes of several 0.1 cm³.

To compare the snow properties measured with different methods, we analyzed the correlations and report the coefficient of determination R ². In the case of linear regressions, we also provide the p-values of regression coefficients to describe the significance of the regression slope (assuming a level of significance p = 0.05). The measurement uncertainty of SMP-derived density was assessed with the RMSE in relation to the reference measurement. Moreover, we provide the mean absolute error (MAE), which weighs individual errors equally, and provides together with the RMSE a measure for the deviation among individual errors.

2.2. Miro-computed tomography and finite element modelling

μCT allows to reconstruct the full 3-D microstructure of snow (Schneebeli and Sokratov, 2004), which can be used in a subsequent finite element analysis to determine linear elastic properties (e.g. Schneebeli, 2004; Srivastava and others, 2010; Köchle and Schneebeli, 2014; Wautier and others, 2015). The μCT scans were performed with different nominal resolutions (voxel size) between (10 μm)³ for new snow samples and (18 μm)³ for depth hoar. The size of the scanned volume was selected manually to ensure representative elementary volumes, which ranged between (5.9 mm)³ for new snow and (7.1 mm)³ for depth hoar. The attenuation image (grey scale image) resulting from each scan was filtered using a Gaussian filter (σ = 1 voxel, kernel half-width = 2 voxel) following Kerbrat and others (2008) and then segmented into a binary image. The threshold for segmentation was constant for each sample and determined visually.

From the binary image, snow density (ρ _CT) was determined from the ice volume fraction and the density of ice. By standard means of a voxel-based finite element solution of the equations of elasticity for composites (Garboczi, 1998) the elastic stiffness tensor was computed. The same numerical method has been previously used in Srivastava and others (2016) or Gerling and others (2017). In the present analysis, we used a bulk modulus of 8.9 GPa and a shear modulus of 3.52 GPa, corresponding to polycrystalline ice at −16°C (Petrenko and Whitworth, 1999). The components of the elastic stiffness tensor E of the snow samples were then determined by averaging stresses and strains over the sample volume according to Garboczi (1998). We consider a Cartesian coordinate system with the two axes (1,  0,  0) and (0,  1,  0) spanning up the horizontal plane and the third axis perpendicular (0,  0,  1) and in the opposite direction of gravity. For the present analysis, we assumed snow to be a transversely isotropic material (Shertzer and others, 2011) and calculated E ₁₁ = E ₂₂, E ₃₃, ν ₁₂ = ν ₂₁, ν ₁₃ and ν ₃₁ from three FE simulations. In each simulation, a uniaxial strain of 0.001 was prescribed in one of the three coordinate directions. Solving a system of linear equations, the elastic moduli E ₁₁ = E ₂₂, E ₃₃ for a transversely isotropic material were recovered.

From the penetration resistance signal of the SMP, however, we do not obtain directional material properties. In fact, the deformation measured at the tip of the device during snow penetration is a superposition of the deformation of many individual elements (Johnson and Schneebeli, 1999). Hence, directional differences are obscured and with the current signal interpretation, only a single mean estimate of the modulus is retrieved from the SMP penetration resistance signal. In contrast, with FE simulations of μCT images, the components of the modulus can be computed and the directional differences of the material properties of snow become apparent. In the following, we show how to derive a modulus that relates to the deformation at the penetrometer tip. For comparison with SMP measurements, we will translate the transverse isotropic material properties obtained from the finite element analysis of μCT images into an isotropic-equivalent modulus. In other words, this modulus refers to an isotropic material showing the same strain vector as the transversely isotropic material and is obtained after setting equal the compliance forms of the stress-strain relations. Solutions only exist for particular stress states, two of which we describe in the following.

Under the assumption, the stress has a horizontal and a vertical component and the relation between them is a, i.e. σ = (aσ ₃₃,  0,  σ ₃₃,  0,  0,  0) in Voigt notation, an isotropic equivalent modulus E _iso can be derived for an orthotropic material:

(1)$$ E_{\rm iso} = {\displaystyle{E_{11}E_{33} (1 - a^2)} \over {E_{33} - a^{2} E_{11} + aE_{33}\nu_{13} - aE_{11}\nu_{31}}}. $$

In the case of the SMP, the ratio between the vertical and horizontal stresses a can be expressed by the SMP cone half angle (30°) and the friction coefficient between ice and steel (Johnson and Schneebeli, 1999) yielding 0.65 with a friction coefficient of 0.05. As we consider snow as a transversely isotropic material with the axis of symmetry along (0,  0,  1) we know, E ₃₃ν ₁₃ = E ₁₁ν ₃₁ and Eqn 1 simplifies to:

(2)$$ E_{\rm iso} = {\displaystyle{E_{11} E_{33} (1-a^2)} \over {E_{33} - a^2 E_{11}}}. $$

Alternatively, if the stress is confined to the two dimensions not including the axis of symmetry, i.e. σ = (bσ ₂₂, σ ₂₂,  0,  0,  0) in Voigt notation with b the ratio between the horizontal components, the isotropic equivalent modulus for a transversely isotropic material (i.e. E ₁₁ = E ₂₂) with the axis of symmetry along (0,  0,  1) equals the elastic modulus of the isotropic plane:

(3)$$ E_{\rm iso} = {\displaystyle{E_{11} E_{22} (1-b^2)} \over {E_{22} - b^2 E_{11}}} = E_{11}=E_{22}. $$

Both assumptions have been used to describe the stress state during penetration in the past. Johnson and Schneebeli (1999) described the stress exerted by the SMP tip with a vertical and one horizontal component. More recently, Ruiz and others (2017) used the penetration model by Bishop and others (1945) who described cone penetration with radial and circumferential stresses only. Hence, in this case, the stress is confined to the isotropic plane of a snow sample with transversely isotropic material symmetry. In the following, we will consider both cases as it is presently not clear, which better describes snow penetration with the SMP. We thus refer in the following to either Eqn 3, when the stress is confined to the isotropic plane or to Eqn 2, when the stress is not confined to the isotropic plane.

2.3. Snow micro-penetrometry

The SMP contains a high-resolution force sensor driven into the snow cover at constant speed which measures the penetration resistance of snow layers. By interpreting fluctuations of the penetration resistance as a Poisson shot noise process as proposed by Löwe and van Herwijnen (2012), the SMP signal can be described by three microstructural parameters, namely the rupture force f, the deflection at rupture δ and the structural element size L. The parameters were calculated over a moving window w of 2.5 mm with 50 % overlap and then averaged over layers which were manually identified. This setting was sufficient to capture even the thinnest layers in our dataset, i.e. the window size was small enough to resolve weak layers of several millimetres. Given the SMP measurement resolution, dealing with even thinner layers seems feasible, if the finer window and evaluation step sizes are selected. By visual inspection signals of lower quality were excluded following the procedure outlined in Pielmeier and Marshall (2009). Signal processing included correcting for penetration resistance offset determined from the noise level measured in the air signal prior to snow penetration. For the comparison with the μCT-derived snow properties, we averaged the SMP-derived microstructural properties, as up to four SMP measurements were available around one μCT sample. We used SMP devices of version 2, among others the device SMP 21.

2.3.1. Density

It has been shown that snow density is related to median penetration resistance $\widetilde {F}$ (Proksch and others, 2015; Kaur and Satyawali, 2017). Proksch and others (2015) modified the relation by Pielmeier (2003) by including the structural element length L to account for differences in snow microstructure. We used the formulation presented by Proksch and others (2015):

(4)$$ \rho_{{\rm SMP}} = a_1 + a_2 \ln(\widetilde{F}) + a_3 \ln(\widetilde{F}) L + a_4 L, $$

and their original coefficients a ₁ = 420.47, a ₂ = 102.47, a ₃ = −121150, and a ₄ = −169960 for $\widetilde {F}$ in N, L in m and ρ in kg m⁻³.

2.3.2. Elastic modulus

Based on foam theory, Johnson and Schneebeli (1999) suggested calculating an elastic modulus from the microstructural parameters f, δ and L according to

(5)$$ E_{{\rm macro}} = {\displaystyle{\,f} \over {\delta L}} {\displaystyle {\delta} \over {L}}. $$

where the first term is the microscale modulus and the second term (δ/L) accounts for the number of grains in contact with the SMP tip. Marshall and Johnson (2009) applied their formulation to penetration resistance signals of natural snow samples and obtained values for the micro-scale modulus in the range of 0.1–1 MPa i.e. about one order of magnitude lower than corresponding literature values such as provided in Mellor (1975). Also, in more recent comparisons SMP-derived microscale moduli were similarly low and not clearly related to those determined from acoustic wave propagation (Capelli and others, 2016; Gerling and others, 2017). These discrepancies were suggested to be partly due to a missing assumption for the particle shape in the SMP analysis (Marshall and Johnson, 2009).

In addition to Eqn 5, we suggest an alternative formulation for an effective modulus that can be derived from the microstructural parameters f, δ and L. To this end, based on dimensional arguments, i.e. considering a deformation energy per volume, the elastic modulus can alternatively be described as:

(6)$$ E_{{\rm SMP}} = a_5 {\displaystyle {\,f \delta} \over {L^3}} {\displaystyle {\delta} \over {L}}. $$

In other words, a single rupture event is assumed to release an energy proportional to the product fδ, which increases with increasing modulus. Scaling by the structural element volume (L ³) a microscale modulus is recovered. Again, the microscale modulus is multiplied by ${\displaystyle {\delta } \over {L}}$ to derive the macroscale modulus. The above formulation of the microscale modulus, i.e. the first part of Eqn 6 coincides with the formulation suggested by Johnson and Schneebeli (1999) for the macroscale strength. We will compare the two formulations Eqns 5 and 6 and determine the coefficient a ₅ in the Results section from comparisons with μCT measurements and literature data.

In the FE simulations of μCT samples, the linear elastic response of the material to a given strain is calculated, whereas the SMP-derived effective modulus may still contain contributions of non-elastic deformation – despite the relatively high strain rate. In fact, the strain rate during SMP penetration can be estimated to ${\dot \epsilon } \approx 100\,{\rm s}^{-1}$ assuming a typical deflection at rupture δ = 10⁻⁵ m, a typical structural element size L = 10⁻⁴ m and the penetration speed of the SMP (20 × 10⁻³ m s⁻¹). Hence, snow grains in front of the SMP tip are likely to fail after very little deformation, i.e. in a brittle manner. Nonetheless, from the SMP signal analysis, we do not obtain a truely elastic modulus, but rather an effective modulus. This is partly, because we do not account for geometric scaling due to indenter size (Huang and Lee, 2013). Moreover, we neglect possible snow compaction near the SMP tip (van Herwijnen, 2013; LeBaron and others, 2014) which leads to an overestimation of the deformation energy in the initial phase of deformation. Hence, with Eqn 6 we calculated a deformation energy fδ containing not only elastic contributions and eventually obtain values lower than the true elastic modulus. This said we refer to the SMP-derived effective modulus from Eqn 6 by E _SMP.

2.3.3. Slab modulus

To compare with the PTV-derived effective moduli, an elastic modulus of the slab consisting of several layers was calculated using a FE model with the SMP-derived layer properties as input (Reuter and Schweizer, 2018). The properties include snow layer thickness, density (from Eqn 4) and effective elastic modulus E _SMP (from Eqn 6). The geometry and loading conditions are similar to a PST and the model was previously used by van Herwijnen and others (2016). The base of the slab has rigid support on one side and is unsupported along a short crack on the other side. Increasing the length of this crack the change in mechanical energy was determined. Fitting pairs of crack length and mechanical energy with a beam equation (Heierli, 2008), we obtained the SMP-derived elastic modulus of the slab, which corresponds to the modulus of a uniform slab showing the same displacement as the layered slab. This modulus is termed SMP-derived elastic modulus of the slab or just SMP-derived slab modulus $E_{{\rm SMP}}^{{\rm slab}}$ (cf. Table 3).

Table 3. Overview of measured snow properties (rows) and measurement techniques (columns). Orange cells correspond to the SMP–μCT dataset, blue cells to the PST–SMP dataset

2.3.4. Specific fracture energy

The specific fracture energy (w _f,SMP) was determined by integrating the penetration force over the moving windows of width w and searching the window within the weak layer (WL) with the lowest value (Reuter and others, 2015a):

(7)$$ w_{{\rm f},{\rm SMP}} = {\rm min}_{\rm WL} \int_{-{w}/{2}}^{+{w}/{2}} F {\rm d}z. $$

In other words, the derivation of the fracture energy is based on the average penetration resistance, since the integral of the penetration resistance across a window of 2.5 mm can be expressed as the average of the measured SMP penetration resistance across this window times the window size. Under the assumption of elastic-brittle rupture events, the first cumulant of the shot noise process, i.e. the average of the measured penetration resistance, equals fδ/L ³. Hence, our alternative formulation to estimate the elastic modulus (Eqn 6) is as well based on the average SMP penetration resistance. This may seem paradoxical, but both mechanical properties are strongly related to material strength. As will be shown below, such SMP-derived estimates are useful to estimate the elastic modulus and the fracture energy of snow.

2.4. Propagation saw test and particle tracking

All PST experiments had similar dimensions (cross-slope width: 0.3 m, up-slope length: 1.2–2.2 m) and had slope normal, rather than vertical, column ends following recommendations and procedures of van Herwijnen and others (2016). A crack was introduced into the weak layer in up-slope direction with a 2 mm thick snow saw until a self-propagating crack started at a certain cut length, the critical cut length r _c. Numerous black markers were inserted in the snow above and below the weak layer and experiments were recorded with a video camera fixed on a tripod. Displacements were derived with PTV from all markers in the slab until the critical crack length was reached. The mechanical energy was calculated according to van Herwijnen and others (2016) from the product of snow density, thickness and average measured displacement of each row of makers. Snow layer densities and thicknesses were available from an adjacent manual snow profile. From the obtained pairs of mechanical energy and corresponding crack length the effective modulus of the snow slab ($E_{{\rm PTV}}^{{\rm slab}}$) was obtained by fitting an adjusted analytical expression for the mechanical strain energy provided by Heierli (2008).

Based on finite element simulations, van Herwijnen and others (2016) showed that the analytical energy formulation does not fully describe the influences of slope angle and the ratio of crack length to slab thickness. For typical values of slope angle and slab thickness encountered in field experiments, a correction factor is required to obtain an adjusted mechanical energy. We applied the corrections described by van Herwijnen and others (2016). Values of the mechanical energy measured in those field experiments likely overestimate deformation due to non-elastic contributions (van Herwijnen and others, 2016). The strain rate during PST experiments was estimated by van Herwijnen and others (2016) to between ${\dot \epsilon } = 6 \times 10^{-5}$ and 2 × 10⁻² s⁻¹ with a mean of 10⁻³ s⁻¹, i.e. in the range of the ductile-to-brittle transition for snow failure (Narita, 1980; Schweizer, 1998; Reiweger and Schweizer, 2010). The PTV-method does, therefore, not provide a truly elastic modulus, but rather an effective modulus, which we call PTV-derived effective modulus of the slab $E_{{\rm PTV}}^{{\rm slab}}$.

The specific fracture energy of the weak layer (w _f,PTV) was calculated as the first derivative of the mechanical strain energy with respect to the critical crack length, which includes the previously derived effective modulus. The specific fracture energy could thus only be derived if crack propagation occurred to the end of the column (43 out of 83 experiments).

3. RESULTS AND DISCUSSION

The SMP–μCT dataset included snow density and elastic modulus calculated from μCT images; the PST–SMP field dataset comprised the elastic modulus of the slab and the specific fracture energy derived from the PTV technique (Table 3). These values were compared to SMP-derived values, namely density, effective modulus and specific fracture energy that were available for both datasets. The following sections take the reader through the rows of Table 3 and contain each, first the results and then a discussion.

3.1. Density

Snow densities derived from the SMP signal (Eqn 4) were clearly related to those obtained from the μCT (R ² = 0.74, p _slope < 0.001; Fig. 1). The SMP-derived density was somewhat overestimated, in particular, for higher densities (≥200 kg m⁻³). The uncertainty of the SMP-derived density as indicated by error bars in Fig. 1 were determined from side-by-side measurements with the SMP in the same snow sample (up to four per sample). Averaging the standard deviation of those side-by-side measurements yielded a RMSE of 9 kg m⁻³. The deviation between the two measurement methods could not be explained by the structural element length (L) (Fig. 1).

Fig. 1. μCT- and SMP-derived snow density with least squares fit (solid line) and 1:1 line (dotted line) for 27 snow samples of the SMP-μCT dataset (full coloured circles). Error bars indicate standard deviation from side-by-side SMP measurements. Data cover a broad range of different snow types indicated by the structural element size L (colour bar). Dots show the values obtained with the model of Kaur and Satyawali (2017). Crosses show the values obtained after recalibrating the model of Proksch and others (2015) to the presented data.

Sources of the uncertainty of μCT-derived snow density are in part related to the segmentation process (the SMP-CT dataset from Riche and Schneebeli (2013) and Proksch and others (2015) depend on) to obtain binary images and the measurement resolution (Hagenmuller and others, 2016), which was chosen according to the snow structure. Proksch and others (2016) reported an agreement between volumetric measurements obtained with density cutters and μCT measurements of 5–9%. Assuming μCT as the more accurate measurement method and hence considering it as a reference, the accuracy of density prediction from SMP-derived parameters was estimated to RMSE = 49 kg m⁻³ and MAE = 41 kg m⁻³. The observed RMSE is somewhat higher than that reported in Proksch and others (2015) (30 kg m⁻³), although we used the same parameterization. Their fitting model is based on Alpine, Arctic and Antarctic data and dominated by the polar samples. Our data only include Alpine samples and are actually a subset of the data presented by Proksch and others (2015). Recalibrating the model of Proksch and others (2015) to our data subset, however, did not substantially change the correlation (R ² = 0.74). Applying the purely penetration resistance (F) based model of Kaur and Satyawali (2017) to our data, yielded a lower correlation (R ² = 0.63). Obviously, microstructural information is required to enhance the parametrization of snow density. Apart from the structural element size, other parameters may be useful. During snow sintering experiments, van Herwijnen and Miller (2013) found an increase of measured penetration resistance with time related to an increase in bond-to-grain ratio, while snow density remained essentially constant. To correctly reflect this, a microstructural property compensating the increase of the average penetration resistance F would be required – apart from L, which rather refers to the average distance between breaking microstructural elements.

The reported agreement between μCT- and SMP-derived density is based on the parameterization of Proksch and others (2015). As electronic components and software have changed in more recent versions of the SMP, data derived with different devices will differ. Hence, to reproduce the presented agreement a calibration specific to the SMP device is required to obtain the fitting coefficients for the parameterizations based on the microstructural parameters.

3.2. Elastic modulus

The isotropic equivalent elastic modulus derived from μCT images with the finite element method and Eqn 3 ranged between 5 and 258 MPa. The elastic modulus increased with increasing density (full circles in Fig. 2, R ² = 0.87).

Fig. 2. SMP-μCT dataset: Isotropic equivalent elastic modulus computed from μCT imaging (full circles, Eqn 3) and effective modulus computed from SMP signals (open circles, Eqn 6) versus μCT-derived density for 27 snow samples. Colour bar indicates structural element size L derived from SMP signals as in Figure 1. Also shown, empirical relations of Köchle and Schneebeli (2014) by dotted lines, of Sigrist and others (2006) by dash-dotted line, of Scapozza (2004) by dashed line and of van Herwijnen and others (2016) by full line.

The μCT-derived values of the elastic modulus were comparable to values modelled by Schneebeli (2004) for a sample of sifted snow with a density of 234 kg m⁻³ undergoing “temperature-gradient” metamorphism. During that process the modulus decreased from 226 to 62 MPa within 6 days. In Figure 2, the parameterizations covering the ranges 100–250 and 250–460 kg m⁻³ reported by Köchle and Schneebeli (2014) describe large parts of our μCT-derived values of the elastic modulus, except for low densities. The parameterization presented by Sigrist and others (2006) based on cyclic loading experiments at high strain rates of about ~2.7 × 10⁻²s⁻¹ are the experimental values closest to μCT-derived values. Scapozza (2004) calculated the modulus from the stress increase observed at small deformations (ε = 0.125%) in quasistatic uniaxial compression tests. The parameterization with density, he proposed yields lower values of the modulus and a stronger increase with density (Table 4). The differences compared with the FE results of the modulus are likely related to the rather low strain rates in his experiments ranging between 10⁻⁶ s⁻¹ and 10⁻³ s⁻¹. van Herwijnen and others (2016) derived an effective modulus by PTV of the (layered) snow slab during propagation saw tests. Those values represent among the presented data the ones obtained at the lowest strain rate of about ${\dot \epsilon } = 10^{-3}\,{\rm s}^{-1}$ similar to the ones by Scapozza (2004) and presumably for the highest strain.

Table 4. Fit parameters p ₁ and p ₂ for exponential fit models of the elastic modulus E and density ρ: E = p ₁ · exp (p ₂ · ρ) and coefficient of determination R ².

For the SMP-μCT dataset μCT-derived density and for the PST-SMP dataset values measured with a density cutter were used. Subscripts CT, SMP and PTV indicate μCT, SMP- and PTV-derived values, respectively. The superscript ‘slab’ refers to an effective modulus accounting for the slab layering. Also parameterizations from previous studies shown in Fig. 2 are listed: Scapozza (2004) (Sc), Sigrist (2006) (Si) and Köchle and Schneebeli (2014)(KS < 250, KS > 250), for densities below and above 250 kg m⁻³, respectively. The data of van Herwijnen and others (2016) (vH) were fitted with an exponential model.

We also present the effective modulus calculated from SMP measurements performed adjacent to the samples used for the μCT scans in Figure 2 (open circles). A coefficient a ₅ of 880 (Eqn 6) was determined to match the data of Scapozza (2004) (R ² = 0.75) and was used to graph the data. The SMP-derived values ranged from 0.18 to 41 MPa and increased with density. Also the scatter increased with density yielding R ² = 0.66. For a particular density, the scatter may partly be related to differences in structural element length L (compare blue and yellow open circles in Fig. 2). However, measurement uncertainties and possible temperature effects we did not measure cause the large part of the observed variations of modulus in the graph with density.

The formulation suggested by Johnson and Schneebeli (1999) (Eqn 5) yielded a lower coefficient of determination (R ² = 0.53), more scatter and a decreasing trend of the modulus with densities above 300 kg m⁻³ (not shown). With the presented, alternative formulation (Eqn 6, R ² = 0.66) the increase of the SMP-derived effective modulus with density was in the range of the parameterizations presented by Sigrist and others (2006), Scapozza (2004) and van Herwijnen and others (2016) (see Table 4). Note that the coefficient a ₅ for the SMP-derived effective modulus (Eqn 6) was estimated to fit the laboratory measurements of Scapozza (2004), since these data are numerous and cover a wide range of densities (Table 1). A calibration to the data of van Herwijnen and others (2016) is not straightforward due to different scales between the slab modulus obtained with the PTV method and the SMP-derived modulus for a single snow layer.

The SMP-derived values exhibited more scatter than FE model results of μCT images (cf. R ² in Table 4), probably due to measurement uncertainties. Moreover, for both, SMP- and μCT-derived values, the scatter was larger for densities above about 200 kg m⁻³. This was expected, as density alone cannot completely describe the complex snow structure (e.g. Shapiro and others, 1997; Schneebeli, 2004; Löwe and others, 2013). Our low-density snow rather consisted of decomposed and fragmented precipitation particles. At higher densities ≳ 200 kg m⁻³, on the contrary, our snow samples showed a more diverse snow microstructure and a much broader range with different mechanical behaviour (cf. Srivastava and others, 2016).

Figure 3 compares elastic moduli derived from μCT images and SMP signals of the SMP–μCT lab dataset. Assuming the stress was not in the isotropic plane (Eqn 2), the correlation between μCT- and SMP-derived elastic moduli was fair but significant (R ² = 0.52, p _slope < 0.01). Alternatively, assuming that the stress is confined to the isotropic plane (Eqn 3), as suggested by Ruiz and others (2017), the correlation was stronger and significant (R ² = 0.76, p _slope < 0.01). The SMP effective modulus was derived from Eqn 6 with a coefficient a ₅ = 5600, which we had determined from a linear regression with the isotropic equivalent modulus from μCT images assuming the stress to be confined to the isotropic plane. Employing the density parametrization reported by Köchle and Schneebeli (2014) with SMP-derived snow density yielded a similar R ² = 0.79 when graphed with μCT-derived modulus from Eqn 3, i.e. assuming the stress was in the isotropic plane. Whereas the μCT-derived elastic modulus refers to an average volume of ≈ (5 mm)³, the corresponding SMP-derived values were averaged on a signal of length >1 cm. Error bars show the standard deviation determined from repeated SMP measurements around the μCT sample. The error bars provide an estimate for the errors due to inhomogeneities in the snow sample that affect this comparison. Also, errors due to the procedure to derive the modulus from SMP measurements are included, such as errors in signal inversion or the applied formulation of the elastic modulus. The linear regression analysis indicated that stronger deviations were more often observed in samples with large snow structures (L > 0.5 mm).

Fig. 3. SMP-derived modulus is shown with μCT-derived elastic modulus under the assumption (a) that the stress is not confined to the isotropic plane, i.e. derived from Eqn 2 or (b) that the stress is confined to the isotropic plane, i.e. derived from Eqn 3. Colour bar indicates structural element size L derived from SMP signals. Error bars show standard deviation from repeated SMP measurements. Also shown by grey dots are elastic modulus derived from SMP-derived density with the relation of Köchle and Schneebeli (2014) versus μCT-derived elastic modulus. 1:1 line dashed. N = 27.

The uncertainty involved with μCT imaging depends on segmentation thresholds and measurement resolution, as pointed out by Hagenmuller and others (2016). Uncertainties of the SMP-derived elastic modulus are partly related to the structural parameters f, L and δ which are derived from the SMP signal. Directly computing the elastic modulus from the SMP parameters yielded a similar accuracy as computing the density from the SMP parameters and using a density parametrization of the elastic modulus (e.g. Köchle and Schneebeli, 2014). With the similar accuracy of either approach, it seems the accuracy is limited by the presently used interpretation of snow micro-penetrometry. Concerning the elastic modulus, Srivastava and others (2016) showed that with fabric-elasticity relations the derivation of the modulus based on μCT-measured snow structure could be enhanced. Moreover, SMP-derived signal characteristics such as correlation length and fractal dimension have the potential to describe snow structure (Satyawali and Schneebeli, 2010), which could be used in a statistical model to predict mechanical properties. So far, however, our interpretation of SMP penetration resistance allows us to derive a parameter set containing the rupture force f, the structural element length L and the deflection at rupture δ. Hence, to include further microstructural parameters into parameterizations, we need to advance the current interpretation of snow penetration or possibly derive additional parameters.

3.3. Slab modulus

From SMP measurements and PTV analysis, we obtained elastic moduli of the slab for the PST-SMP field dataset (Fig. 4). With the PTV method, some stratigraphic information is implicitly preserved, resulting in an effective slab modulus. To compare with SMP-derived moduli, we calculated the elastic modulus of the slab with a FE model, simulating the bending behaviour of a beam similar to the PST. PTV- and SMP-derived values of the slab modulus showed both an increase with snow density. PTV-derived values of the effective modulus showed more scatter than SMP-derived values (Table 4, cf. R²). From Figure 4 it becomes apparent that our PST-SMP dataset is a somewhat biased subset of the larger dataset of van Herwijnen and others (2016) who found a clearer relation with density (R² = 0.29).

Fig. 4. PST-SMP dataset: Effective slab modulus computed from SMP signals and FE modelling (open circles) and PTV analysis (full squares) versus SMP-derived snow density. Colours indicate different field sites. Grey open squares represent the entire dataset of van Herwijnen and others (2016).

With the PTV and the SMP method, we obtained lower values compared to the μCT-derived elastic modulus (cf. Fig. 2 and 4). As strain rates during PST experiments are relatively low, i.e. ${\dot \epsilon } \approx 10^{-3}\,{\rm s}^{-1}$ (van Herwijnen and others, 2016), we anticipated that effective moduli derived with the PTV analysis would be low, and certainly lower than the true elastic modulus. On the contrary, effective moduli derived from the SMP signals are determined at high strain rates (${\dot \epsilon } \approx 100\,{\rm s}^{-1}$), but may still include some contribution of non-elastic deformation due to possible snow compaction at the tip. In Figure 4 SMP-derived values coincide with the lower values of the PTV-derived moduli. This, however, is rather due to the chosen calibration (coefficient a ₅ in Eqn 6). Calibrating SMP-derived properties is necessary to account for geometric scaling because the microstructural properties determined from the SMP depend on the indenter, i.e. its size and shape (Huang and Lee, 2013). We decided to calibrate the SMP data to the laboratory measurements of Scapozza (2004) since they cover a broad range of densities, despite the rather low values of strain rate during his experiments.

Figure 5 shows the direct comparison between PTV- and SMP-derived effective moduli from adjacent field experiments. Relating PTV- and SMP-derived values yielded a rather weak correlation (R ² = 0.36, p _slope < 0.01). This relation does not include the origin (0,0). Due to different strain rates and possible time-dependent deformation during the PST, a constant term in the regression may be needed to describe this relationship.

Fig. 5. PTV-derived versus SMP-derived elastic moduli of the slab. Colours indicate field sites with same legend as in Figure 6. Black line represents linear regression. N = 83.

Fig. 6. Weak layer fracture energy derived with PTV analysis and from SMP signals including errorbars; colours indicate different field days and sites. One outlying PTV-derived value (w _f = 3 J m⁻²) excluded. N = 43.

With the PTV method an elastic modulus of the slab was determined for a bending snow column of ≤0.3 m³. From the SMP, we first calculated the elastic modulus (Eqn 6) and then computed a slab modulus using a FE model accounting for slab layering. This allowed comparing values derived from the two methods. After calibration the SMP-derived values of the slab modulus were in a similar range as the PTV-derived values. As shown in Figure 4 the subset of data we use is not representative of the entire dataset presented by van Herwijnen and others (2016). At low densities below $\lesssim \!\!200\,{\rm kg}\,{\rm m}^{-3}$ the data subset is biased towards higher moduli and at higher densities towards lower moduli. As only for this subset SMP measurements were available and in addition, PTV measurements naturally have an uncertainty of ~25%, the direct comparison with the limited data yielded a rather weak correlation. That said, uncertainties such as from snow density measurements likely play a less important role than the uncertainties from marker location tracking during PTV analysis.

3.4. Fracture energy

The specific fracture energy, the crucial weak layer property in view of crack propagation, describes the amount of energy needed to expand a crack over a given distance in the weak layer and form two new crack surfaces. The specific fracture energy derived from the SMP signal was related to the fracture energy derived from PST experiments with PTV (R ² = 0.40, p _slope < 0.001; Fig. 6). PTV-derived values were generally higher by a factor of 3.5 than SMP-derived values of w _f.

In contrast to PTV-derived moduli (Fig. 4 and 5), the fracture energy was only determined for those PST experiments where the crack propagated to the end of the column (43 of 83 cases). The error estimation was based on the repeated SMP measurements in the snow samples of the SMP–μCT dataset. From these samples, we estimated the variance of side-by-side penetration resistance measurements. As the derivation of the fracture energy is based on penetration resistance (Eqn 7), the uncertainty can thus be estimated by propagating the uncertainty of the measured penetration resistance. PTV-derived values were similar to those reported by Schweizer and others (2011) (w _f = 1.3 ± 0.8 J m⁻²), who modelled the specific fracture energy with finite element simulations with effective moduli from SMP microstructural parameters derived after Marshall and Johnson (2009). In general, PTV-derived values overestimated the specific fracture energy as they are higher than literature values for ice; for example Schulson and Duval (2009, p. 206) estimated w _f ≈ 0.32 J m⁻². As discussed in van Herwijnen and others (2016), the high values of fracture energy obtained with the PTV method suggest that a substantial fraction of the potential energy of the slab is dissipated during PSTs. Accounting for visco-plastic effects in the slab and in the weak layer is likely required to obtain more realistic values. Lower values of fracture energy w _f = (0.07 ± 0.02) J m⁻² were obtained by Sigrist and Schweizer (2007), who used a FE model with calibrated SMP-derived values of the effective modulus (Sigrist and others, 2006). Schweizer and others (2011) discussed the importance of the choice of the modulus for these methods. Some methods, however, offer a way to determine the fracture energy and are independent of the modulus – and the measurement uncertainties involved. LeBaron and Miller (2014) presented a method based on μCT reconstructions of snow samples and calculated the critical energy release rate along the minimum cut surface (Hagenmuller and others, 2014b). Another approach circumventing assumptions in particular for the modulus, was presented by McClung (2015). He demonstrated that the specific fracture energy may be obtained from field measurements based on the ‘shear band model’ originally presented by Palmer and Rice (1973). For his avalanche crown dataset he obtained values of w _f between 0.08 and 0.36 J m⁻² and similarly for a large dataset of 591 PSTs he reported values between 0.04 and 0.1 J m⁻². These two methods, namely the μCT-based and ‘shear band model’-based method, may allow calibrating SMP-derived values of the fracture energy independent of the modulus. So far, however, we did not calibrate SMP-derived values. Although the relationship between SMP- and PTV-derived values provides confidence in the SMP-derived property, calibrating with the PST–SMP dataset would shift the SMP-derived fracture energy towards unrealistically high values.

4. CONCLUSIONS

To advance snow instability modelling, accurate measurements of snow mechanical properties are required – in particular, snow density, elastic modulus and fracture energy. We, therefore, presented comparisons of these properties based on a PST–SMP and a SMP–μCT dataset. To facilitate comparisons between the measurements of the elastic modulus we employed two techniques: A slab modulus allowed comparing the SMP-derived elastic modulus with the PTV-derived effective slab modulus. For comparison with the SMP-derived elastic modulus, the anisotropic elastic modulus derived from μCT images was translated into an isotropic equivalent modulus.

Snow density was well reproduced in repeated SMP measurements. The data were in line with values derived from μCT image analysis for a range of seasonal snow types.

Regarding the elastic modulus, the presented alternative approach to derive the elastic modulus from SMP signals yielded a similar increase with density as values derived from μCT data, PTV data and previously published parametrizations. Transforming the μCT-derived elastic properties for a transversely isotropic material into an isotropic equivalent modulus allowed a direct comparison of SMP-derived and μCT-derived elastic moduli. This comparison showed good agreement of the elastic modulus derived from the SMP with the alternative approach in relation to μCT-derived values. The agreement of SMP- and μCT-derived elastic modulus was better, in particular, if the assumption was used that the stress due to snow penetration with the SMP is confined in the isotropic plane. Hence, our findings suggest that the SMP-derived elastic modulus reflects the elastic modulus perpendicular to the axis of penetration.

SMP-derived elastic modulus and fracture energy were related to PTV-derived values. Lacking a reliable method and with PTV-derived values of the fracture energy values too high compared with literature values of ice, assessing the SMP-derived data is difficult. Nonetheless, the specific fracture energy can likely be estimated from the SMP penetration signal, however, for calibration a dataset including an independent technique – possibly based on simulations of μCT-derived snow structure (LeBaron and Miller, 2014, 2016) – would be required.

Our results show that from the SMP signal the snow properties required to model failure initiation and crack propagation can be obtained. All relevant properties, in particular density, elastic modulus, specific fracture energy and strength (cf. Marshall and Johnson, 2009) were related to values obtained with other independent methods. However, the finding that the SMP signal contains information on the modulus components perpendicular to the direction of penetration shows that we not fully comprehend SMP–snow interaction during penetration. Better understanding the penetration process can potentially enhance signal analysis and allow calibrating SMP-derived properties to finally bridge different measurement scales.