Temporal evolution of crack propagation propensity in snow in relation to slab and weak layer properties

If a weak snow layer below a cohesive slab is present in the snow cover, unstable snow conditions can prevail for days or even weeks. We monitored the temporal evolution of a weak layer of faceted crystals as well as the overlaying slab layers at the location of an automatic weather station in the Steintälli field site above Davos (Eastern Swiss Alps). We focussed on the crack propagation propensity and performed propagation saw tests (PSTs) on seven sampling days during a two-month period from early January to early March 2015. Based on video images taken during the tests we determined the 10 mechanical properties of the slab and the weak layer and compared them to the results derived from concurrently performed measurements of penetration resistance using the snow micro-penetrometer (SMP). The critical cut length, observed in PSTs, overall increased during the measurement period. The increase was not steady and the lowest values of critical cut length were observed around the middle of the measurement period. The relevant mechanical properties, the slab effective elastic modulus and the weak layer specific fracture, overall increased as well. However, the changes with time differed suggesting 15 that the critical cut length cannot be assessed by simply monitoring a single mechanical property such as slab load, slab modulus or weak layer specific fracture energy. Rather, crack propagation propensity is the result of a complex interplay between the mechanical properties of the slab and the weak layer. We then compared our field observations to newly developed metrics of snow instability related to either failure initiation or crack propagation propensity. The metrics were either derived from the SMP signal or calculated from simulated snow stratigraphy (SNOWPACK). They partially 20 reproduced the observed temporal evolution of critical cut length and instability test scores. Whereas our unique dataset of quantitative measures of snow instability provides new insights into the complex slab-weak layer interaction, it also showed some deficiencies of the modelled metrics of instability – calling for an improved representation of the mechanical properties.


Introduction 25
Dry-snow slab avalanche release is governed by failure processes within the layered snow cover. Whether a failure initiates and a resulting macroscopic crack will propagate, depends on the complex interaction between slab layers, the weak layer and to some extent the substratum, i.e. the layers below the weak layer. For example, the thickness and characteristics of the slab layers determine the magnitude of the stress at the depth of the weak layer due to a skier (Habermann et al., 2008;Monti et al., 2016;Thumlert and Jamieson, 2014), and also how much deformation energy can be released to drive crack propagation (Gaume et al., 2014b;Heierli et al., 2008;McClung, 1979). On the other hand, a weak layer is a prerequisite for failure initiation and offers a path for crack propagation. Sigrist and Schweizer (2007) first described the interaction of slab and weak layer properties for evaluating the critical cut length. By interpreting their results in a fracture mechanical framework they concluded that the energy that has to be exceeded to fracture a weak layer depends on the material properties 5 of the weak layer, whereas the energy that is available for crack propagation mainly depends on the material properties of the overlaying slab, and may also depend on the collapse height of the weak layer. Given the two most relevant processes in drysnow slab avalanche release, failure initiation and crack propagation (e.g., Schweizer et al., 2003a), van Herwijnen and  suggested a conceptual model on the effect of the slab properties, in particular slab depth, on these failure processes. With increasing slab depth conditions for failure initiation become less favourable whereas conditions for crack 10 propagation become more favourable.
Temporal changes in snow instability hence stem from changes in slab and weak layer properties -separately or in combination. For example, during a snowfall the probability of failure initiation in the weak layer increases due to the additional load. However, the additional load will also promote strengthening of the weak layer (e.g., Zeidler and Jamieson, 2006a) -though the strength increase may lag behind loading during a snowfall. Changes of weak layer strength have been 15 studied in detail (e.g., Chalmers and Jamieson, 2001;Gauthier et al., 2010;Jamieson and Johnston, 1999;Zeidler and Jamieson, 2006a, b) and more frequently than changes of slab properties. However, to the best of our knowledge, there are hardly any studies that investigated how temporal changes affect the complex interplay between slab and weak layer properties.
Repeated measurements of weak layer shear strength revealed how various types of weak layers gain strength over 20 time . Overall, an increase in strength was almost always observed in all studies while at the same time the load usually increased. Jamieson and Schweizer (2000) reported the shear strength over time for 19 buried surface hoar layers. Typical weak layer strength gain was on the order of 100 Pa d -1 during the initial weeks after burial (Schweizer et al., 1998). Occasionally measured decreasing strength with time was attributed to spatial variability within the study site or errors associated with measurement technique (Jamieson and Johnston, 1999). 25 Many of the above mentioned studies monitoring strength changes of weak layers focussed on relating stability trends with observed local or regional avalanche activity (e.g., Jamieson et al., 2007). However, observed avalanche activity was often not related to the stability index, i.e. the ratio of strength to stress (e.g., Föhn, 1987), calculated from the study plot measurements. While the shear strength of the weak layer is important for failure initiation, dry-snow slab avalanches release due to crack propagation which requires the release of deformation energy stored in the slab layers. This conceptual 30 mismatch has long been recognized. For example, Schweizer et al. (1998) pointed out that since the shear frame measurements will primarily provide information on the strength and strength changes of weak layers, the 'effective reactivity (propagation potential)' depending on the slab characteristics should be assessed by supplementary tests.
Nevertheless, temporal changes of slab characteristics were rarely monitored, apart of course, from the load.
The temporal evolution of shear strength has been modelled for persistent and non-persistent weak layers (e.g., Chalmers, 2001;Conway and Wilbour, 1999;Hayes et al., 2005;Lehning et al., 2004;Zeidler and Jamieson, 2006a, b) so that the evolution of the stability index can be monitored or even forecasted (Giraud, 1993;Vernay et al., 2015). Conway and Wilbour (1999) suggested a model for natural avalanches during storms by comparing the load to the strength by assuming that failure occurs at the base of the new snow layers; their strength solely depends on density which increases with 5 increasing overburden stress. Föhn and Hächler (1978) exclusively focused on slab properties as they studied how the slab layers settle during major snow storms. They proposed to follow the settlement (coefficient) over time to assess the probability of large natural dry-snow avalanches.
With the development of the propagation saw test (PST; Gauthier and Jamieson, 2006), now a well-established snow instability test exists providing a quantitative test result, the critical cut length. Furthermore, the PST allows 10 determining the relevant slab and weak layer properties Schweizer et al., 2011;Sigrist and Schweizer, 2007;van Herwijnen et al., 2016). Recently, Birkeland et al. (2014) repeatedly performed propagation saw tests on a layer of buried surface hoar; they focussed on conditions for fracture arrest.
A number of studies have focussed on the temporal evolution of spatial patterns on small uniform slopes -inter alia testing the hypothesis that variability should increase in the absence of major external forcing such as a snowfall (Birkeland 15 and Landry, 2002;Birkeland et al., 2004;Logan et al., 2007). Hendrikx et al. (2009) used the Extended Column Test (Simenhois and Birkeland, 2009) to investigate spatial variations of the propagation potential. They assessed the spatial variability of two sites each on two days and found increased spatial clustering on the second sampling day.
For clarification, we shortly define the following two terms: snow instability and crack propagation propensity. As Reuter et al. (2015) have pointed out both high failure initiation and high crack propagation propensity are required to 20 describe unstable snowpack conditions. More recently, Reuter et al. (2016b) suggested to complement the failure initiation and crack propagation criteria with a tensile criterion related to dynamic crack propagation. Hence, snow instability cannot be assessed with a single stability criterion as suggested in the past (e.g., Föhn, 1987) but only by a combination of indices related to the essential processes in dry-snow slab avalanche release. With regard to crack propagation propensity, we refer to this term to describe (1) in general, whether snowpack characteristics favour self-sustained crack propagation possibly 25 resulting in a snow slab avalanche, and (2) when interpreting propagation saw tests, whether the critical crack length is less than about one third of the column length and the crack propagates to the end of the column.
The aim of the present study is to repeatedly measure the slab and weak layer properties in a study plot to monitor their temporal evolution and to investigate their interaction in view of assessing snow instability. During the winter 2014-2015 we followed a layer of faceted crystals that was responsible for wide-spread avalanche activity in the region of Davos 30 (Eastern Swiss Alps) over the course of two months. We performed propagation saw tests, which we analysed based on the video recordings of the tests and compared the results to concurrently performed measurements of penetration resistance using the snow micro-penetrometer (SMP) (Schneebeli and Johnson, 1998). As we performed our measurements in a level study plot equipped with an automatic weather station, we also simulated the evolution of snow stratigraphy with the numerical snow cover model SNOWPACK. Finally, we compared our observations to metrics of instability derived from either the SMP signals or simulated with SNOWPACK. The acquired dataset provides a comprehensive time series of quantitative measures of snow instability; it allows insight into the complex interplay between slab and weak layer properties that jointly govern snow instability.

Methods 5
We followed the evolution of a weak layer of faceted crystals in the level study plot surrounding the automatic weather station WAN7 (2442 m a.s.l.) located in the Steintälli field site above Davos,eastern Swiss Alps (46.808° N,9.788° E).
Measurements were performed on eight days between 6 January and 3 March 2015, typically once a week during the two month study period (Table 1).

Weak layer formation 10
On 1 December 2014 the manually observed snow profile at the study plot Weissfluhjoch (2540 m a.s.l.) (located 3 km to the northeast of WAN7) showed a melt-freeze crust with 1 cm recently fallen snow on top. On 2-3 December 2014, a minor storm accumulated an additional 12 cm of new snow. During the following two weeks, the snow above the crust settled and transformed into a layer of faceted crystals due to near-surface faceting (Birkeland, 1998); this layer of faceted crystals was buried by a snowfall on 16 December 2014. The 2-3 cm thick melt-freeze crust was consistently found throughout the winter 15 below the layer of faceted crystals that formed the 5-8 cm thick weak layer. As of mid January 2015 the layer was classified as rounded facets (FCxr) with a grain size of 1-1.5 mm. Above the weak layer was another layer of faceted crystals overlain by well consolidated slab layers that had formed in late December 2014 and early January 2015 ( Figure 1). The weak layer was likely responsible for wide-spread avalanche activity in the region of Davos on 30-31 December 2014 when many natural dry-snow slab avalanches were observed. Since no snow profiles at fracture lines were recorded, we do not know the 20 depth and type of failure layer. However, the weak layer we monitored consistently failed in snow instability tests in early January 2015, and there were no other prominent weak layers within the snow cover. This observation suggests that it was also the primary failure layer of the late-December avalanches.

Field measurements
On each of the eight sampling days we observed a manual snow profile, including layer density, according to Fierz et al. 25 (2009). The detailed density profile is required for the analysis by particle tracking velocimetry (PTV) of the PSTs (see below). The manual snow profile served as reference for most other measurements and was completed with the two snow instability tests that can easily be performed in flat terrain: the compression test (CT) (Jamieson, 1999) and the extended column test (ECT) (Simenhois and Birkeland, 2009).
On each sampling day, except for 19 February 2015, we also conducted at least three propagation saw tests in the immediate vicinity of the snow profile. The tests were performed according to Greene et al. (2010), albeit with longer columns. Initially, column length was around 1.5 m, and it increased towards 2 m at the end of the study period as the weak layer became more deeply buried (slab thickness increased from 59 to 148 cm) Gaume et al., 2015). Using a 2-mm thick snow saw, we cut the layer of faceted crystals at its upper interface, where CT and ECT results indicated that 5 the failure occurred. Black markers (2.5 cm in diameter) were inserted into the snowpack; we filmed all tests with a video camera for subsequent analysis by particle tracking velocimetry . The video recording also allowed us to accurately determine the critical value of the crack length r c OBS , when the crack started to propagate.
Furthermore, we could also assess whether the tests were properly performed, e.g., whether the saw cut was within the weak layer. In some cases, we had to discard a test result since the cut was not performed consistently close to the layer interface; 10 this resulted in only two values of the critical cut length on 21 January and 3 February 2015.
In the PST experiments, cracks did not always propagate to the end of the column, but arrested with a slab fracture (Table 1). These propagation results are termed END and SF, respectively (e.g., Greene et al., 2010). In the following, while reporting the critical cut length when crack propagation initiates, we do not differentiate between these two propagation results since the critical crack length is independent of the subsequent phase of dynamic crack propagation. Whether or not a 15 running crack will arrest may depend on the tensile strength of the slab (Gaume et al., 2015;Schweizer et al., 2014).
However, the onset of crack propagation entirely depends on the balance between the energy available for fracture, i.e. the mechanical energy released due an incremental advance of the crack, and the specific fracture energy, i.e. the energy required for crack growth, or in other words the resistance to crack propagation.
Concurrently, we performed several SMP measurements (SMP version 2), at least three at the location of the 20 manual snow profile, and at least one at each of the PST locations; thus in total at least 6 measurements per sampling day.
The SMP measurements at the PST locations were conducted before isolating the columns, close to the end of the column where the saw cut was initiated. This procedure allowed comparing the SMP-derived properties with those from the PTV analysis.

PTV analysis 25
Using a particle tracking velocimetry algorithm (PTV; Crocker and Grier, 1996) the displacement field of the slab prior to crack propagation was calculated from the video image; it shows bending of the slab due to the saw cut (of length r) during the PST. Based on the displacement field of the slab, the effective elastic modulus of the slab E *PTV and the specific fracture energy of the weak layer w f PTV were determined (van Herwijnen and van Herwijnen et al., 2016). The approach is based on the work by Heierli et al. (2008) who suggested for the geometry of the PST (of unit width) an 30 expression for the total energy of the system V(r) as a function of crack length r which consists of the fracture energy V f (r) and the mechanical energy V m (r): The mechanical energy includes two terms, a fracture mechanical and a bending term: where w f is the specific fracture energy, D is the slab thickness, γ is a constant of about one, depending on Dundur's elastic mismatch parameter, ′ = /(1 − ν 2 ) is the plane strain elastic modulus of the slab, and ν the Poisson's ratio (assumed ν = 0.2). The load of the slab on the weak layer consists of the shear stress τ = -rgD sinθ and the normal stress with η = �4(1 + ν)/5 .

10
However, by comparing estimates obtained with the analytical expression for the mechanical energy provided by Heierli et al. (2008) with finite element (FE) simulations, van Herwijnen et al. (2016) recently showed that Eq. 2 underestimates the mechanical energy for realistic values of the ratio of crack length to slab thickness r/D. Therefore, they introduced a correction factor that accounts for the sensitivity to r/D and the slope angle θ to obtain an adjusted mechanical energy * ( ). 15 According to theorem of Clapeyron, for a linear-elastic system, the mechanical energy can be estimated from the gravitational potential energy. The gravitational potential energy was computed from the measured displacement field, assuming a layered slab and using the manually observed density profile, to determine the mechanical energy as a function of crack length. This measured mechanical energy was then fitted to the adjusted mechanical energy * ( ) to obtain the effective elastic modulus of the slab E *PTV , the fit parameter. It is defined as the modulus of a uniform slab of equal mean 20 density yielding the same displacement as the real slab. To determine the specific fracture energy w f PTV , the analytical expression for the adjusted mechanical energy with the best fit modulus is differentiated with regard to the crack length at the critical cut length r c OBS . In other words, the slope of * ( ) at the critical crack length corresponds to the specific fracture energy of the weak layer. For a more detailed description on deriving the effective elastic modulus and the specific fracture energy, the reader is referred to van Herwijnen et al. (2016). 25

SMP signal processing
We used the penetration resistance data acquired with the SMP to obtain detailed data on the layering of the snow cover to derive mechanical properties following the approach described in Reuter et al. (2015).
Based on the manually observed snow profile, layers were manually defined from the corresponding sections of the SMP signal, i.e. several slab layers, a weak layer and a basal layer. The shot-noise model by Löwe and van Herwijnen (2012) 30 was then applied to determine snow micro-structural parameters, namely the rupture force f, the deflection at rupture δ and the structural element size L. These three parameters were calculated over a moving window of 2.5 mm with 50% overlap and averaged over each layer. Furthermore, for each layer the snow density was derived according to Proksch et al. (2015), allowing to calculate the load on the weak layer. The micro-mechanical modulus for the slab layers and the strength of the weak layer σ WL were calculated according to Johnson and Schneebeli (1999). The weak layer specific fracture energy w f SMP 5 was calculated as the minimum of the integrated penetration resistance across each moving window within the weak layer (Reuter et al., 2013). Finally, the penetration depth PS was estimated by integrating the penetration resistance F from the snow surface to the depth PS where a threshold value of the absorbed energy was reached .
The effective elastic modulus of the slab E *SMP was determined by performing FE simulations of the experimental setup taking into account the slab layering (for details see Reuter et al., 2015). The FE model consisted of all the slab layers, 10 the weak layer and a basal layer, each with density, micro-mechanical modulus and thickness values derived from the SMP measurement. The mechanical strain energy V m FEM (r) was then calculated for a stratified slab bending over a crack of increasing length r. In order to recover an effective elastic modulus E *SMP , the analytical expression for the adjusted mechanical energy * ( )was fitted to the modelled values of mechanical energy V m FEM (r). Hence, we followed the same approach as for the PTV analysis, and we also used the newly developed correction factor to obtain the adjusted mechanical 15 energy * ( ).
The SMP-derived weak layer specific fracture energy and the effective elastic modulus were scaled by a linear factor of 2.8 and 2.5, respectively, to approximately match the corresponding values derived from the PTV analysis. Scaling of the specific fracture energy and the modulus was performed as there is no calibration yet of the microstructural mechanical properties that can be derived from the SMP signal. Whereas the microstructural properties derived from the 20 SMP are physically based, they cannot be expected to directly represent the corresponding macroscopic properties (Reuter et al., 2016a).
In addition, an alternative effective elastic modulus * SMP was derived, following the same approach with the same FE model as described above. However, rather than using the micro-mechanical modulus, for each layer the modulus was determined using the SMP-derived density and applying the parametrization provided by Scapozza (2004) Based on the mechanical parameters estimated from the SMP measurements, two metrics of point instability were derived, as suggested by Reuter et al. (2015).
The first metric is the failure initiation criterion S, a strength-over-stress criterion describing the propensity of the weak layer to fail in case of skier loading: 30 where σ WL is the SMP-derived micro-mechanical strength of the weak layer and ∆τ FEM the maximum additional shear stress at the depth of the weak layer due to skier loading. The maximum shear stress at the depth of the weak layer was modeled with the 2D linear elastic FE model originally presented by Habermann et al. (2008) to calculate the shear stress ∆τ FEM below a layered slab due to the weight of a skier.
The second metric is the crack propagation criterion r c SMP , an SMP-derived critical cut length. It is derived by numerically solving Eq. 7, which is obtained by finding the extremum of the total energy of the cracked system V(r) (Eq. 1) with respect to r and ensuring that it is a maximum: ( ) = + ( ) = 0, which yields (Schweizer et al., 2011): 5 with 0 = 3 2 4 2 , 1 = � + 3 2 � 2 + 3 2 + 2 , 2 = 2 + 9 2 + 3 2 2 , 3 = 3 2 , 4 = 3 2 . By inserting the effective elastic modulus of the slab * and weak layer specific fracture energy w f SMP into Eq. 7 the crack propagation criterion r c SMP was obtained. Hence this modelled critical cut length is independent of the critical cut 10 length r c OBS measured in the field. For a more detailed description on how to obtain the above mentioned mechanical properties as well as the SMP-derived metrics of point instability, the reader is referred to Reuter et al. (2015).

Snow cover modelling
We compared results from field measurements to output of the numerical snow cover model SNOWPACK (e.g., Lehning et al., 2004) driven by meteorological input from the automatic weather station WAN7. This weather station is located in the 15 study plot where we performed the field measurements. Meteorological input contained air temperature and relative humidity (Rotronic MP100H HygroClip, ventilated), wind speed and direction (YOUNG wind monitor), incoming short and long wave radiation (Campbell CNR1), and snow height (Campbell SR50). Data gaps shorter than one day were filled by linear interpolation. Gaps longer than one day were filled with data from the nearby AWS Weissfluhjoch (2540 m a.s.l.; 3 km to the northeast). Variables were filtered by introducing reasonable lower and upper limits, e.g. 5 and 100 % for 20 relative humidity. The modelling time step was 15 min after resampling the data from the AWS with a sampling rate of 10 min. The model was initiated on 1 October 2014, when no snow was present at the AWS and ran until the end of May.
Neumann boundary conditions for estimating the snow surface temperature and atmospheric stability corrections for estimating turbulent exchange were the preferred adjustments concerning the energy balance model (Stössel et al., 2010).
From the model output, we evaluated the skier stability index SK38 introduced by Jamieson and Johnston (1998). It is defined as the ratio of shear strength to shear stress: SK38 = +∆ with τ the shear stress due to the weight of the overlaying slab, ∆τ = 155/h the additional shear stress due to a skier with h the slab depth (Föhn, 1987;Monti et al., 2016), and s the shear strength as parameterized with density and grain type according to Jamieson and Johnston (2001).
In addition, we estimated the critical cut length r c SNP from the snow stratigraphy provided by SNOWPACK. We 5 used SNOWPACK model output to derive all the required mechanical properties of the snow layers. The critical cut length was then estimated based on the relation given by Gaume et al. (2014a, Eq. 5) (see Gaume et al., 2016 for a detailed derivation). Based on discrete element modelling they suggested the critical cut length, in the flat (for slope angle θ = 0°), to essentially depend on the plain strain elastic modulus of the slab ′ , slab load σ, and weak layer shear strength : with Λ = ( ′ WL / WL ) 1/2 , D the slab thickness, WL the weak layer thickness and WL its shear modulus. All the required parameters are provided by SNOWPACK, except for the elastic moduli ′ and WL . For the shear modulus of the weak layer we assumed a constant value of 0.5 MPa, based on laboratory experiments (Camponovo and Schweizer, 2001;. For the elastic modulus of the slab, we followed the same FE approach as described above for the SMP analysis to derive an effective elastic modulus taking into account slab layering rather than using a slab modulus based 15 on the average slab density. Hence, for each layer of the modelled snow stratigraphy an elastic modulus was calculated from density based on the relation provided by Scapozza (2004). With these properties (modulus, layer density and thickness) a FE simulation was performed to determine the effective elastic modulus E *SNP .

Avalanche activity
Study plot measurements are commonly used in operational forecasting to make assessments about the avalanche danger in 20 the surrounding terrain (e.g., Gauthier et al., 2010). We therefore compared the results obtained from the field measurements to the observed avalanche activity in the region of Davos. These observations include the number, type and size of avalanches recorded by personnel from the local ski areas, the avalanche warning service, SLF staff members, and others.
We then calculated the avalanche activity index per day as described by Schweizer et al. (2003b). The index is a weighted sum of the number of observed avalanches per day including natural as well as artificially triggered avalanches; the weights 25 depend on the avalanche size and are 0.01, 0.1, 1 and 10 for Canadian size classes 1 to 4, respectively (McClung and Schaerer, 2006).

Case studies
To explore the complex interaction between slab and weak layer properties on the critical cut length in a PST, we considered a few cases with exemplary temporal evolutions of slab load, slab elastic modulus and weak layer specific fracture energy. 30 To this end we numerically solved Eq. 7 to obtain the corresponding critical crack length, and hence its evolution over time.
These examples are not meant as a sensitivity study where one parameter is varied and the other held constant, but as an illustration of the interaction when most or all parameters change.

Propagation saw test results 5
On 6 January 2015, when we did the first propagation saw tests, cracks did not always fully propagate to the end of the column, but slab fractures were observed in three tests ( yielded a short cut length, 17 cm, the shortest value recorded during our sampling period. Subsequently, the cut lengths increased with time.

Load
On 6 January 2015, the weak layer was buried below a slab of 59 cm (total snow depth HS = 115 cm) and the initial load was about 1.4 kPa (Table 1, Figure 2b). The load did not change much during the following week but then continuously 20 increased due to snowfalls to almost 3 kPa by early February 2015. After a fair weather period in February with no new snow for more than two weeks, the snow depth increased again, and on the last sampling day (3 March 2015) the load was about 3.9 kPa.
The density derived from the SMP signal agreed well with the manually measured density (not shown) and accordingly the increase in load above the weak layer was well represented by the SMP measurements (Figure 2b). For the 25 numerical snow cover model SNOWPACK, values of load as calculated from average density and slab thickness were often slightly lower than those measured in the field. The underestimation is mainly due to the fact that the modelled snow depth was about 25 cm lower than measured at the location of the manual snow profiles.

Effective elastic modulus of the slab
The effective elastic modulus of the slab was derived from the bending of the slab during the propagation saw test via the PTV analysis E *PTV (Figure 3a) as well as from the SMP signal analysis using the FE model with either the micromechanical modulus or a modulus derived from SMP density, yielding * SMP and * SMP , respectively (Figure 3b). The effective elastic modulus of the slab obtained from the PTV analysis showed an overall increase from about 2.5 to 10 MPa 5 during the two-month sampling period. However, the increase was not steady; for example, on 13 February 2015 relatively low values between 2.4 and 6.1 MPa were obtained.
The SMP-derived effective elastic modulus * SMP initially did not change much with median values of approximately 3 MPa during the first 16 days. It then increased until the end of January to about 5 MPa. However, subsequently, it decreased and on the last measuring day low values of only about 1.8 MPa were derived. 10 The alternative SMP-derived modulus * SMP also increased from initially 8 MPa to about 22 MPa in early February.
Thereafter * SMP did not change much anymore with values between 19 and 26 MPa (median values per day).

Weak layer specific fracture energy
The PTV analysis suggests that the weak layer specific fracture energy w f PTV overall increased with time from about 0.6 J m -2 to about 2.2 J m -2 , mostly after the end of January 2015 (Figure 4a). The SMP analysis revealed a similar tendency of 15 increasing weak layer specific fracture energy with time. Overall, w f SMP increased from about 0.5 J m -2 to 1.4 J m -2 and most of the increase occurred towards mid February (Figure 4b).

SMP-derived metrics of instability
The failure initiation criterion S derived from SMP resistance data was initially about 300, indicating a transitional value for failure initiation propensity given the threshold reported by Reuter et al. (2015). They observed that a value of the initiation 20 criterion of about 230 divided between the cases with and without concurrently observed signs of instability. The index then increased to about 600 towards the end of January and further to about 1100 by the end of the sampling period ( Figure 5a).
The SMP-derived critical cut length r c SMP overall increased from about 40 cm to 70 cm (Figure 5b). Again the increase was not steady with a similar tendency with time as shown for the weak layer specific fracture energy w f SMP . At the three sampling days between end of January and mid February similar values of r c SMP were derived, namely about 55 cm. 25

SNOWPACK derived metrics of instability
The skier stability index SK38 (Figure 5c) for the investigated weak layer showed values between 1.1 and 1.35 in early January 2015; these values are in the range of transitional stability (1 to 1.5) according to Jamieson and Johnston (1998).
With the snowfalls at the end of January and early February ( Figure 1) the index then decreased towards 1.0 indicating increased triggering probability. Subsequently, during much of February, when there was no snowfall, the skier stability 30 index SK38 gradually increased towards 1.3 and decreased again to overall the lowest values of about 1 towards the end of the observation period.
The modelled critical cut length r c SNP (Figure 5d), on the other hand, steadily increased from about 20 cm in early January to 60 cm in early March 2015.

Avalanche activity 5
The highest avalanche activity in the region of Davos was observed around the end of the year 2014 (Figure 6), before our first measurements were performed. In January 2015, avalanches were occasionally observed, in particular on 18 January, a sunny Sunday after a snowfall. Towards late January and early February, there was a period of increased avalanche activity.
The last peak was on 3 March 2015, again after a major snowfall, but during this period avalanches did no longer run on the weak layer of facets we monitored; the observed fracture depths were much lower than the burial depth of the weak layer we 10 investigated. Overall, avalanche activity clearly decreased since early January 2015. However, after each significant snowfall avalanches were again observed (e.g. on 18, 28 and 31 January, Figure 6), i.e. the weak layer of facets was still reactive, a true persistent weak layer.
The overall decrease in avalanche activity until the end of February is in line with the observed results of the CTs and ECTs we performed concurrently on each sampling day. The scores increased from values just below 20 taps to values 15 of around 30 taps (red asterisks in Figure 6). On the first measuring day, no crack propagation was observed with the extended column test (Table 1). Subsequently, the fracture type in ECTs was always P (propagation across the entire column), except on 28 January when one ECT did not fracture (X). The increased avalanche activity towards the end of January and early February coincides with the lowest values of observed critical cut length.

Case studies 20
The observed temporal evolution of the critical cut length in our PSTs showed an unsteady increase with the lowest values in the middle of the measuring period and an increase thereafter. To explore whether this type of pattern is possible at all, we numerically solved Eq. 7 to find out how the critical cut length changes with time for various scenarios of the temporal evolution of the specific fracture energy of the weak layer w f , the load and the modulus of the slab E.
In the first simplified scenario, all input parameters were assumed to increase, corresponding to a situation when the 25 slab thickens and becomes stiffer, while at the same time the weak layer toughens due to the increased load (Figure 7a). In this scenario, the critical cut length did almost not change, showing a very slight increase (Figure 7b). The combination of increasing load and increasing slab modulus provided a bit more energy to drive the crack, but just as much about to compensate the increased specific fracture energy of the weak layer. This indicates that in this scenario the load had a larger effect than the modulus, since an increasing modulus reduces the amount of strain energy available due to less deformation. 30 In the second simplified scenario, the specific fracture energy of the weak layer increased as the load increased, while the slab modulus remained constant (Figure 7c). This scenario resulted in a decreasing critical cut length (Figure 7d). advance the crack due to the weak layer toughening.
In the third simplified scenario, the load remained constant, while the slab modulus and the weak layer specific fracture energy increased ( Figure 7e); this scenario corresponds to a fair weather period where the slab stiffens due to settlement and the weak layer toughens due to ongoing sintering with time. In this scenario, the critical cut length 5 substantially increased (Figures 7f). Due to the higher slab modulus less energy was released to drive crack propagation while at the same time the weak layer became tougher.
Finally, in our last scenario, we attempted to roughly mimic the temporal evolution of the critical cut length as observed in our PSTs. We assumed the load to increase approximately as observed, the modulus to triple and the specific fracture energy of the weak layer to increase as well, but not continuously (Figure 7g). With these assumptions, the critical 10 cut length first increased, had a local minimum approximately in the middle of the considered time period and finally increased again (Figures 7h).
Alternatively, the sensitivity of the critical cut length could be explored with Eq. 8. In fact, this would reveal the same trends for the critical cut length, if the temporal evolution of the shear strength followed the one assumed for the weak layer specific fracture energy. 15

Discussion
We repeatedly performed propagation saw tests in a level study plot to follow the temporal evolution of a weak layer of faceted crystals and of the overlaying slab layers by combining state-of-the-art measurement techniques and numerical snow cover simulations. Specifically, we used particle tracking velocimetry, the snow micro-penetrometer and the snow cover model SNOWPACK to estimate snow mechanical properties and derive snow instability criteria. 20

Observed critical cut length (PST)
While performing the propagation saw tests in the study plot, we initially observed a mixture of END and SF test results; mixed test results are typically not found in critical conditions. Only when the load had reached 2 kPa, all cracks fully propagated towards the end of the column. This observation suggests that the tensile strength of the slab was initially not large enough to allow crack propagation (Gaume et al., 2015), in particular at the more shallow locations. Other changes in 25 slab properties that might have favoured slab fractures, e.g. faceting, were not observed.
The initial absence of full crack propagation in our field tests contrasted with the high avalanche activity observed during that time. The contrasting observations may be due to the fact that we made our measurements in a wind-sheltered study plot that may not be very representative of wind-affected starting zones in the area. In typical lee-slope starting zones, the tensile strength of the slab might have been larger; moreover, the tensile stress might be lower on slopes since less 30 bending is expected prior to natural avalanche release than observed in PSTs in flat terrain (McClung and Borstad, 2012). In general, slab layers are more variable in terms of penetration resistance than weak layers reflecting the dynamic conditions of snowfall and wind during deposition of the slab . Since the properties of the slab layers are particularly important for crack propagation, variable crack propagation propensity has to be expected. The observed discrepancy highlights some of the limitations when extrapolating instability from flat field study sites.
As of 28 January 2015 all cracks in PSTs propagated to the very end of the column. This is in line with the results 5 of the other stability test performed: sudden fractures (SP/SC) in CTs and full crack propagation (P) in ECTs.
Overall the observed critical cut length doubled from about 25 cm to about 50 cm by early March. However, the increase was not steady. Despite the overall increase in critical cut length, the lowest values were observed at the end of January and in early February. This pattern of the temporal evolution of the critical cut length in PSTs was rather unexpected since weak layers typically gain strength with time (e.g., Jamieson and Schweizer, 2000). 10 Consistently low values of the critical cut length in PSTs, between 19 and 24 cm, were observed on 28 January 2015. On the following sampling day, the lowest value (17 cm) was observed, but also a rather high value of 38 cm. In general, low values or scores in snow instability tests are more trustworthy than high values. In the case of the propagation saw test, any measurement and observation errors increase the cut length. In our case, we were cutting at the top of a weak layer of faceted crystals, and the layer above was not much harder and as well consisted of faceted crystals. Hence, it was 15 relatively easy to move the saw out of the weak layer. We were able to identify such deviations on the video recordings, but only on the side facing the camera. The high values may therefore be due to imperfect sawing and show the difficulty in obtaining consistent PST results in weak layers which are not very well defined.
The median range, i.e. the difference between the highest and the lowest value, of the PST results on a given sampling day was 5.9 cm, so that the resulting uncertainty in the mean is about 2-3 cm. Reuter and Schweizer (2012) 20 reported a standard deviation of the critical cut length on days with surface warming of about 5 cm. Similar variations for PST results on a single day at a single location have been reported by Gauthier and Jamieson (2008). Therefore, we have no reason to dismiss these low values or attribute them to measurement errors. However, they may be related to spatial variations in weak layer and slab properties within the study plot. In fact, on 14 January 2015 we performed the PSTs at a location were the snow depth was below average compared to the rest of the study plot. On all other 25 days snow stratigraphy was very similar, exemplified by mostly similar SMP profiles (not shown). Similar snow stratigraphy in study plots has, for instance, been shown by Pielmeier and Schneebeli (2001). Previous snow instability studies performed in level study plots suggested that measurements in general are reliable and that the effect of spatial variations is relatively small. Jamieson (1995) reported a mean coefficient of variation of shear strength of about 15% for sets of measurements in study plots. Correspondingly, variations in stability indices derived from study plot measurements of load and shear strength, 30 two measurements that have comparable errors as the PST, were found to be indicative of avalanche activity (e.g., Jamieson et al., 2007). Though the PST results may be influenced by some small scale spatial variability of the snowpack in the study plot, we deem it unlikely that the observed minimum values are entirely the result of spatial variability, but indicate in fact a period of high propagation propensity. This interpretation is supported by the increased avalanche activity observed in late January and early February when many skier-triggered avalanches were observed (Figure 6).

Load, modulus and specific fracture energy
The temporal evolution of the parameters influencing the critical cut length, namely the load, the effective elastic slab modulus and the weak layer specific fracture energy, all exhibited similar, mostly increasing trends. The load obviously 5 increased (Fig. 2b) and the results derived from the manual density measurements and the SMP profiles agreed well -in accordance with recent comparison studies (Proksch et al., 2015;Proksch et al., 2016).
The PTV-derived effective elastic modulus also increased, but the increase was not steady with some low values on 13 February 2015 (Figure 3a). These low values were also observed with the SMP (Figure 3b). The SMP-derived effective elastic modulus * SMP , however, showed very low values on the very last sampling days -resulting in an overall decrease. 10 When estimating the effective elastic modulus from the SMP-derived density, the agreement in temporal evolution was better. Indeed, overall * SMP also increased. The observation that the SMP-derived modulus occasionally decreases compared to the previous measurement day, for example, on 3 February and 3 March 2015 can be explained by the fact that on these days each a snowfall period ended. With the addition of new snow on the top of the existing slab, the effective elastic modulus in general decreases, unless the old slab layers below substantially settle so that the penetration resistance 15 increases. In other words, the additional load due to the new snow leads to a thicker slab of lower average density, but also resulting in a lower effective elastic modulus.
In general, the weak layer specific fracture energy is expected to increase with time as sintering will increase bonding between crystals -unless the weak layer is subject to a large temperature gradient (e.g., van Herwijnen and Miller, 2013). Indeed, the PTV-derived weak layer specific fracture energy w f PTV and the microstructural related specific fracture 20 energy derived from the SMP signal w f SMP exhibited a similar overall increasing trend ( Figure 4). However, independent of the evaluation method, the increase occurred mostly after early February. This observation suggests that the weak layer toughening mainly occurred during the fair weather period without additional loading in February. A constant weak layer specific fracture energy in combination with additional loading by snowfall in January would have resulted in a decreasing critical cut length, in line with field observations. 25 With respect to the absolute values of the specific fracture energy, these should be considered as effective values since it is clear that they are too high compared to the specific fracture energy of ice (McClung, 2015). For the PTV-derived values, the discrepancy is most likely related to the fact that the observed bending in a PST includes non-elastic parts of deformation, thereby increasing the specific fracture energy to physically unrealistic values; for an in-depth discussion of this issue see van Herwijnen et al. (2016). 30 Compared to the weak layer specific fracture energy, the PTV-and SMP-derived values of the effective elastic modulus agreed less well. In particular, it is known that * SMP does not relate well to the PTV-derived modulus (Reuter et al., 2013). On the other hand, van Herwijnen et al. (2016) recently showed that the PTV-derived modulus fits relatively well with results from laboratory experiments in the same range of strain rates. The SMP-derived modulus * SMP includes the complex interaction between the cone and surrounding snow microstructure and does obviously not reflect simple elastic deformation, but breaking, jamming and other local effects occur (LeBaron et al., 2014;. Obviously, the SMP-derived modulus * SMP , which relies on the relatively robust density derivation, agreed somewhat better with the PTV-derived modulus. 5 Whereas the temporal evolution can well be compared, the absolute values cannot, since the SMP-derived values of * SMP and w f SMP were scaled to the corresponding PTV values to ease comparison. The scaling was performed using a larger dataset (unpublished) of side-by-side PST and SMP measurements. It is beyond the scope of this study to provide a quantitative comparison of the two methods; it will be the subject of a publication currently in preparation (Reuter et al., 2016a). To conclude the discussion on the PTV and SMP-derived values, below we present an error assessment -as far as 10 this is possible.
The errors associated with the parameters derived from the PTV analysis (i.e. the measurement uncertainty) were about 20% (or about 1 MPa) for the modulus and about 18% (or about 0.14 J/m 2 ) for the weak layer specific fracture energy.
These estimates are based on calculating these properties 1000 times for each experiment accounting for the uncertainty in the input parameters (uncertainty in the distance measurements in the field, density measurements, and especially the 15 location estimates of the dots in the PTV analysis). van  showed that the reproducibility of side-byside measurements for the slab modulus was good (even though they have greater uncertainty), whereas the reproducibility for the weak layer specific fracture energy was lower. The better reproducibility for the modulus might be related to the fact that the modulus is an integrated property over all slab layers and hence may be less sensitive to spatial variations of one layer. 20 The errors of the SMP-derived parameters are more difficult to assess. However, Proksch et al. (2015) showed that the SMP-derived density is a reliable measure. Their finding is supported by our measurements showing good reproducibility between SMP-derived and manually measured load (Figure 2b). However, in particular the derivation of the effective elastic modulus is more prone to errors. In general, the variability of SMP-derived microstructural parameters is rather high. Thus, any value which is indirectly derived from these microstructural parameters will exhibit large variations; 25 in particular values of the deflection at rupture δ are relatively unreliable (Löwe and van Herwijnen, 2012).

Metrics of instability
The SMP-derived metrics of point snow instability, the failure initiation criterion S and the crack propagation criterion r c SMP , were recently developed and validated . Here, we applied them for the first time to monitor the temporal evolution of instability. The failure initiation criterion S increased with time suggesting that initiating a 30 failure in the weak layer became increasingly difficult during the sampling period. This tendency is supported by the results of the two small column snow instability tests (Table 1, Figure 6). CT/ECT scores increased from around 20 to 30 taps. The absolute values of S were rather high, in the range that was considered as rather stable by Reuter et al. (2015). Again, this is in agreement with the rather high scores of the CTs and ECTs and is related to the fact that by the end of January 2015 the weak layer was buried below a thick slab of more than 1 m. Weak layers buried deeper than 1 m are not frequently triggered by skiers (e.g., Schweizer and Jamieson, 2007).
The crack propagation criterion r c SMP overall increased in agreement with the observations (Figure 2a). It showed a 5 similar evolution with time as the SMP-derived weak layer specific fracture energy, which is used to calculate r c SMP . A relative decrease towards the end of January was also observed but was less prominent than for r c OBS . Despite substantial scatter, until early February about half of the values were between 35 and 45 cm. Only after mid February r c SMP values increased to about 70 cm. Considering the threshold values indicating unstable conditions S < 234 and r c SMP < 41 cm suggested by Reuter et al. (2015), the two criteria predict unstable conditions in early January. In fact, signs of instability 10 were observed by the field team on 6 January and 28 January 2015; on the latter day the lowest values of critical cut length were observed in PSTs. However, at that day the failure initiation criterion was already high (S ≈ 590), since the slab thickness was >1 m.
Based on the simulated snow stratigraphy provided by SNOWPACK, we presented two corresponding metrics of instability. Overall, snow stratigraphy was well simulated, as exemplified by the comparison shown in Figure 1. Despite 15 large differences in vertical resolution, the simulated SNOWPACK profile, the manually observed profile and the SMP profile qualitatively agreed well.
The SNOWPACK-derived skier stability index SK38 varied between about 1.0 and 1.35 and did not show any trend with time. This is in contrast to the increasing scores obtained with CTs and ECTs. Whereas the skier stability index SK38 initially increased in agreement with the observation, it subsequently mainly reflected whether there was any loading by new 20 snow. After the end of January when the weak layer was deeply buried, the SK38 did no longer indicate skier triggering but yielded almost the same value as the natural stability index (not shown) since ∆τ, the additional shear stress due to a skier, became negligibly small. Accordingly, increasing load due to snowfall resulted in a decreasing SK38 as shown towards the end of February and in early March. In other words, the SK38 became dominated by the static shear stress. In contrast, the SMP-derived initiation criterion S does not include the static shear stress, and hence showed a different behaviour with time. 25 As Schweizer and Reuter (2015) pointed out, the dynamic load, rather than the static load due to the weight of the slab, is essential for initiating a failure due to the well-known deformation rate dependence of snow strength (e.g., . The modelled critical cut length r c SNP (Eq. 8) based on recent findings by Gaume et al. (2016) steadily increased over the two-month period. All three essential variables, namely the load, the slab modulus and the weak layer shear strength 30 overall increased with time. Whereas an increase of the load -all other variables remaining unchanged -would result in a decrease of the critical cut length, increasing slab modulus as well as increasing weak layer shear strength leads to increasing critical cut length. The latter two variables are directly related to snow density which in general increases with time -except for the average slab density which may temporarily decrease after a snowfall. In our case, with regard to the critical crack length r c SNP , the effect of increasing load was clearly over-compensated by the effects of increasing slab effective modulus and increasing weak layer shear strength.
The discrepancy with the observed critical cut length seems to be due to the weak layer shear strength which strongly increased with increasing density from initially 0.9 to 1.9 kPa in early March -though persistent weak layers are known to hardly settle despite increasing overburden pressure due to their anisotropic microstructure (e.g., Reiweger and 5 Schweizer, 2010;Walters and Adams, 2014). In fact, the SMP analysis of the weak layer strength (not shown) suggests that the increase was not as prominent as shown by SNOWPACK where the shear strength is parameterized based on the extensive shear frame measurements by Jamieson and Johnston (2001). Whereas Gaume et al. (2016) recently showed that the modelled critical cut length r c SNP of the weak layer we tested was lower than of the adjacent layers -suggesting that the modelled critical cut length can well discriminate between weak layers and other layers within a given snow stratigraphy, it 10 seems more challenging to predict changes over time since small changes of the contributing variables may decide whether the critical cut length increases or decreases.

Case studies
In general, considering the fracture mechanical approach (Anderson, 2005) reveals that in a first order approximation c ~f 2 ; the critical crack length decreases with increasing load, and increases with either increasing slab elastic modulus or 15 increasing specific fracture energy of the weak layer. To explore the complex interaction of these parameters described by Eq. 7, we presented four case studies. They show how the temporal evolution of the load and the modulus of the slab, and the specific fracture energy of the weak layer affects the changes of the critical cut length with time. The most influential parameter seems to be the load. Even with an increasing weak layer specific fracture energy, the increasing load caused the cut length to decrease (Figures 7c,d). This consistent decrease was however not observed in the field, where only towards the 20 end of January the critical crack length decreased and then clearly increased towards the end of the measurement period.
This suggests that the increase in slab modulus and/or specific fracture energy (over)-compensated the effect of the increasing load. The fourth case, more or less mimicking the observed changes with time of the three variables, shows that under certain conditions the cut length may increase or decrease with time, primarily due to subtle changes of slab properties. Exploring Eq. 8 relating slab modulus and load, and weak layer shear strength provided very similar results (not 25 shown) and confirms the findings of the case studies.

Conclusions
We monitored the temporal evolution of a weak layer of faceted crystals that was one of the critical weaknesses during the winter 2014-2015 in the region of Davos (Eastern Swiss Alps). We focused on the crack propagation propensity and performed propagation saw tests on seven sampling days during a two-month period from early January to early March 30 2015. Tests were completed with objective measurements, namely by resistance profiles acquired with the snow micro-penetrometer and particle tracking velocimetry based on video images of the PSTs. Our dataset represents the first comprehensive time series of quantitative measures of critical cut length and related mechanical properties of slab layers and weak layer. In addition, we compared our field observations to newly developed metrics of instability either derived from the SMP signal or from modelled snow stratigraphy.
The critical cut length, observed in the PSTs, showed an overall tendency to increase with time. However, the 5 lowest values were observed towards the end of January/early February. These low values were not expected and seem to be the result of the complex interaction of slab and weak layer properties -rather than of measurements errors or large spatial variations within the study plot.
The relevant mechanical properties, the effective elastic modulus of the slab and the weak layer specific fracture energy, overall increased, independent of the evaluation method. Only the SMP-derived effective elastic modulus * SMP 10 showed a different behaviour. However, the increase was not steady and these small changes likely affected the critical cut length as exemplified with the case studies. These findings suggest that it is not possible to assess the critical cut length, or in general crack propagation propensity, by simply monitoring a subset of these mechanical properties. One has to consider the complex interaction between the effective elastic modulus of the slab, the load due to the slab, and the weak layer specific fracture energy. Furthermore, these properties have to be determined with good accuracy since otherwise reliably modelling 15 the critical cut length is not possible.
We then applied traditional (such as the SK38) and newly-developed metrics of snow instability describing either the failure initiation or the crack propagation propensity. These metrics were derived from the SMP signal or calculated from simulated snow stratigraphy (SNOWPACK) and partially reproduced the observed temporal patterns. Whereas the SMPderived initiation criterion appropriately indicated that triggering probability overall decreased, the skier stability index 20 provided by SNOWPACK did not show this tendency, but indicated the period of increased avalanche activity towards the end of January. The predicted critical cut lengths, r c SMP and r c SNP , overall increased with time. Whereas the SMP-derived propagation criterion r c SMP partly reproduced the unsteady increase with some lower values towards the end of January, the SNOWPACK-derived critical cut length r c SNP did not show any of the observed variation in critical cut length, apart from the overall increase. 25 Whereas the PST combined with the PTV approach seems to provide the most reliable measure of propagation propensity and corresponding mechanical properties, the procedure is time consuming. However, this disadvantage can only be outweighed, if modelled metrics of instability become more reliable. This will require further validation studies -best performed by comprehensive field measurements in study plots equipped with an automatic weather station, and possibly an enhancement of the representation of mechanical properties in the model based on cold laboratory studies. 30 an Ambizione grant of the Swiss National Science Foundation (PZ00P2_161329). Tables   Table 1: Snowpack characteristics and snow instability test results on the eight measurements days. For the PST the total number of tests and the number of tests with crack propagation to the end of the column, for the CT the score and the fracture type (SP: sudden planar, SC: sudden collapse), and for the ECT the fracture type (X: no fracture, N: initiation, but no propagation, P: propagation to column end) and the score are given (Greene et al., 2010 19 Feb 2015 147 n/a n/a n/a n/a n/a n/a  (g,h), the situation during the sampling period is supposed to be roughly mimicked.