Evaluating snow weak-layer failure parameters through inverse finite element modelling of shaking-platform experiments

Snowpack weak layers may fail due to excess stresses of various natures, caused by snowfall, skiers, explosions or strong ground motion due to earthquakes, and lead to snow avalanches. This research presents a numerical model describing the failure of “sandwich” snow samples subjected to shaking. The finite element model treats weak layers as interfaces with variable mechanical parameters. This approach is validated by reproducing cyclic loading snow fracture experiments. The model evaluation revealed that the Mohr–Coulomb failure criterion, governed by cohesion and friction angle, was adequate to describe the experiments. The model showed the complex, non-homogeneous stress evolution within the snow samples and especially the importance of tension on fracture initiation at the edges of the weak layer, caused by dynamic stresses due to shaking. Accordingly, a simplified analytical solution, ignoring the inhomogeneity of tangential and normal stresses along the failure plane, may incorrectly estimate the shear strength of the weak layers. The values for “best fit” cohesion and friction angle were ≈ 1.6 kPa and 22.5–60. These may constitute valuable first approximations in mechanical models used for avalanche forecasting.


Introduction
Dry snow avalanche release mechanics present a key research question.Various mechanical models that have been used to address the dry snow slab avalanche release problem focused on weak layer failure: e.g.crack models inspired by the over-consolidated clay theory (McClung, 1979), cellularautomata models (Fyffe and Zaiser, 2004), fibre-bundle model (Reiweger et al., 2009), physical-statistical models (Chiaia and Frigo, 2009), multiple finite element method, FEM (Stoffel, 2005;Podolskiy et al., 2013), and analytical and empirical models (Zeidler and Jamieson, 2006).Recent studies, based on FEM with interfacial constitutive laws for weak layers, have shown that one of the key uncertainties in avalanche forecasting, spatial heterogeneity of weak layers, can be treated by statistical methods and that its importance is reduced for greater snow slab depths (Gaume, 2012;Gaume et al., 2012Gaume et al., , 2013)).Moreover, merging of FEM with terrestrial laser scanning input data (e.g.Teufelsbauer, 2009Teufelsbauer, , 2011) ) and the growth of computer performance promise that this decade will see the possibility of precise estimation in terms of statistical distributions of potentially unstable snow masses for feeding into models of avalanche dynamics (Naaim et al., 2003).Accordingly, further investigation concerning the weak layer mechanical behaviour and constitutive law, along with their implementation in FEM, are certainly needed for better quantitative understanding of the avalanche formation process.
For studying dry snow slab avalanches, various approaches have emerged and have been employed in FEM models to represent a snow weak layer under a cohesive slab; for detailed review refer to Podolskiy et al. (2013).Previous studies were mainly designed to investigate: (1) the stress state of a snow slab on a slope (Smith et al., 1972;Curtis and Smith, 1974;Smith and Curtis, 1975;McClung, 1979), (2) snow deformation (Lang and Sommerfeld, 1977), (3) skier loading (Schweizer, 1993;Wilson et al., 1999;Jones et al., 2006;Habermann et al., 2008;Mahajan et al., 2010), (4) weak layer heterogeneity, super weak zone length and stress concentration, as well as avalanche release slope Published by Copernicus Publications on behalf of the European Geosciences Union.
Previous FEM studies may roughly be classified into three principally different numerical approaches in terms of representation of weak layers (or potential failure surfaces): (1) a thin isotropic (or anisotropic) continuum (Bader and Salm, 1990;Jamieson and Johnston, 2001;Miller et al., 2011), (2) an interface with zero thickness, which may be vertically "collapsible" or not (McClung, 1979;Stoffel and Bartelt, 2003;Gaume et al., 2012) or (3) a combination of the first two methods as a thin collapsible/non-collapsible layer with interfaces at the bottom and the top of it (Mahajan and Senthil, 2004;Mahajan and Joshi, 2008;Mahajan et al., 2010).The above-mentioned approaches are chosen based on the objectives of the study, and at the same time it may be noted that there is no universal, generally accepted framework for treatment of the "slab-failure surface" system.Due to computational difficulties related to the size of avalanche release zone, it generally appears preferable to represent the weak layer by an interface, since its thickness is significantly smaller than the total snow height.Furthermore, by referring to volumetric layers, the FEM mesh size in the weak layer would have to be smaller than the size of the crystals and thus may put the validity of the continuous approach into question.More importantly, it is known from fracture line studies that poor bonding between layers may be a more significant cause of avalanching than low strength within weak layers (McClung and Schaerer, 2006).Accordingly, the idea of treating weak layers as interfaces appears attractive for large-scale applications and will be explored further in this paper.

Objectives and scope of the study
The aims of the present work are two-fold.First, we revisit the mechanical behaviour of weak layers through FEM simulations of previous experiments on failure of layered snow by Podolskiy et al. (2010b).Second, we show that the well known Mohr-Coulomb failure criterion with cohesion, which includes normal pressure effects and tensile strength (one of the most common approaches in mechanics of granular materials), may be used as a first approximation to reproduce these dynamic experiments.Accordingly, this paper reports on an evaluation of the performance of this failure criterion, as well as an evaluation of associated parameters (cohesion and angle of internal friction), through a detailed numerical-experimental cross-comparison.
The idea of describing the failure of snow according to the Mohr-Coulomb failure criterion emerged since the pioneering studies of Haefeli (1963), Roch (1966), andMellor (1975).According to Mellor's review (Mellor, 1975) cohesion can be associated with time-dependent inter-crystalline bonding (sintering), while the angle of internal friction can be imagined as initial or residual strength of snow with broken bonds.Recently, this criterion was used, for example, for modelling snow erosion by flowing avalanches (Louge et al., 2011), for predicting critical inertial loads for failure of weak layers in seismically active regions (Matsushita et al., 2013;Pérez-Guillén et al., 2013), or for analysing the packing of snow against sensor surfaces caused by wet avalanche (Baroudi et al., 2011).
One important prediction of Mohr-Coulomb failure criterion is the dependence of shear strength on the normal stress imposed to the sample.Mellor (1975) suggested that this criterion could be problematic for snow due to changes of the material state under pressure.Since then, numerous experimental studies investigated the effects of normal load on shear strength of snow and snow weak layers, mainly through shear frame or shear vane tests.Results showing an influence of normal stress on various snow types were reported by Roch (1966), DeMontmollin (1978, 1982), Mc-Clung (1977), Perla and Beck (1983) and Navarre et al. (1992).Jamieson and Johnston (1998) reported similar influence on non-persistent weak layers, but found no significant effect on persistent weak layers, thus proposing φ = 0 • .Recently Matsushita et al. (2012) conducted tests with artificial precipitation snow to investigate temporal variation of the shear strength and concluded that the influence of normal load on the strength was more significant than temperature.Overall, most of these studies investigated the influence of normal pressure using shear-frames.Results obtained with alternative methods, like shaking-platform tests (Nakamura et al., 2010;Podolskiy et al., 2010b), may provide valuable new insights on these issues of normal stress influence and applicability of Mohr-Coulomb criterion, and thus still remain to be analysed in that context.

Shaking-platform experiments
This paper is based on a series of experiments using the shaking platform as described by Nakamura et al. (2010) and Podolskiy et al. (2010b).The procedure can be briefly summarised as follows: snow samples were frozen to the platform and loaded via inertia through horizontal oscillations with a constant amplitude of 1.65 cm and with growing frequency (see Sect. 4.1.4for details).The frequency increase caused gains in velocity, acceleration, and thus stresses within the samples; at the point when the stress exceeded the strength of snow, the sample failed.High-speed video records, accelerometers and measurement of the fractured mass revealed the instant of failure and the corresponding peak acceleration (in the range of 2.23-6.36g) (Nakamura et al., 2010;Podolskiy et al., 2010b).Originally, this dynamic experimental approach was developed for studying the shear strength properties of snow and their relationship to vibrations (Abe andNakamura, 2000, 2005;Nakamura et al., 2000aNakamura et al., , b, 2010)).These previously reported tests were performed on homogeneous blocks of snow.Podolskiy et al. (2010b) introduced a weak layer into the blocks and the possibility to incline the platform at 0, 25 and 35 • ; these two points make the study more relevant to dry snow slab avalanche release.Nevertheless, free surfaces on five sides of the sample and the probable occurrence of stress heterogeneities in the snow block due to edge effects restrict the possibility of simple stress assessments and relating the experimental results to a real snowpack at slope scales.For example, an attempt by Nakamura et al. (2010) to calculate dependence of shear strength on presumably constant overburden pressure produced surprisingly high values of internal friction angle (73.4-83.1 • ) with zero cohesion, thus exemplifying the importance of understanding normal stress variations in the experiments for reliable interpretations.

Nat
The experiments, considered in this study, were performed in a cold laboratory (with an ambient air temperature of −10 • C) on artificial "sandwich" snow samples constituted of two blocks of snow, with a weak layer made of low density snow placed approximately at mid height.The samples were prepared by sieving artificial precipitation snow over a cohesive slab with a density around 234 kg m −3 , covering it with another slab, leaving for 74 h of sintering, and later vertically cutting the resulting structure into smaller blocks.The resulting weak layer density was around 100 kg m −3 , and its thickness was around 1-2 cm.If we attempt to identify the closest type of natural weak layer to the artificially created horizons in the middle of snow samples, it would be a nonpersistent precipitation layer, made of low density, partly decomposed dendrite crystals, or DFdc according to classification by Fierz et al. (2009).The length, width and height of specimens were 0.3, 0.2, and 0.2-0.45m, respectively.The masses overlaying the weak layers ranged between 1.3 and 4.6 kg.This difference in mass was created by varying the height of the upper block.For the inclined tests, the geometry of the sample side cuts was always kept vertical.
From the different types of tests performed by Podolskiy et al. (2010b), we selected here only the weak layer fracture tests made with horizontal single-degree-of-freedom oscil-lations (at the same time we recall that a sample can have various inclinations: 0, 25 or 35 • ).
In total, 19 individual tests with varying properties were modelled.Most relevant parameters and results of experiments are indicated in Table 1; for more details refer to Podolskiy et al. (2010b).

Some further experimental conditions relevant to construction of the model
Four specific points, relevant to the construction of the numerical model, need to be highlighted.First, in the majority of experiments weak layer fractures were observed at the lower interface (between the weak layer and the lower block).No significant vertical collapse within the weak layer could be recognised during tests (based on video quality we could only restrict the maximum possible collapse as less than 1 mm) (Podolskiy, 2010).Moreover, with the particular inertial loading considered in this study, we do not expect vertical collapse to play a major role in the failure process.Hence, for the purpose of simplification, the possibility of vertical collapse was not considered in the modelling and the weak layer was represented as a non-collapsible interface.Second, analysis of video records shows no noticeable horizontal strains in the two snow blocks surrounding the weak layer; with the available video quality, the upper bound for strain is less than 0.33 % (Podolskiy, 2010).This means that the whole blocks can be regarded, as a first approximation, as rigid bodies.Such assumption amounts to considering that most of the possible deformation is concentrated within the weak layer, in agreement with previous studies (Reiweger and Schweizer, 2010;van Herwijnen and Jamieson, 2005;van Herwijnen et al., 2010).For the purpose of reducing computing costs, we will thus omit the lower block in the model and only consider a system made of an upper block and an interface (Fig. 1).Furthermore, in agreement with this discussion, it will be shown (Sect.5.3) that the elastic properties of the upper block do not influence failure properties.
Third, we note that the size of the samples used in the experiments is much smaller than the critical crack length required for failure self-propagation (McClung, 2011;Gaume et al., 2013).In other words, weak layer failure in the experiments is driven only by the applied loading and not by stress redistributions, which remain negligible at the scale considered.Global sample failure occurs when the inertial stresses induced by the oscillations reach the failure criterion in the whole weak layer.Therefore, there is no need for considering the post-failure behaviour of this layer or the progressive accumulation of damage during the successive loading cycles.In this sense, these experiments appear particularly well-suited to focusing on the weak layer failure criterion, independently of post-failure propagation phenomena.
Lastly, since the experiments had high rates of loading (failure occurred within a second; strain rates were higher than 10 −3 s −1 ), we do not consider viscous behaviour of  snow.High loading rates guarantee a brittle range for all observed fractures.Such high loading rates are representative of stress variations involved in natural avalanche releases, loading produced by skiers/snowmobilers, explosive air blasts, as well as strong ground motion due to earthquakes or mine blasting (Podolskiy et al., 2010a).

Model geometry
The upper block is represented by a rectangle or parallelogram (for inclined tests), which is 0.3 m long and 0.14-0.36m high.A 1 cm × 1 cm quadrilateral element shape with four nodes is used for the mesh (QUA4); there are about 14-36 elements in the vertical dimension (depending on the sample height) and 30 in the horizontal dimension.We noted that sensitivity tests with twice as many elements produced similar, but much more computationally costly results.
The weak layer is treated as an interface, modelled by joint elements with four nodes (JOI2) and zero thickness, i.e. an element is created between two segments of two points (Fig. 1c).There are 30 joint elements (each 1 cm long).The "lower" part of the joint ( 1 A -2 B ; Fig. 1c) is fixed to the bottom boundary, meaning that there are no vertical or horizontal displacements relative to the boundary.Displacements on the upper part ( 4 A-3 B) are not restricted, thus allowing free deformation.

Constitutive laws of the block
The upper block is considered as a uniform and isotropic elastic material similar to many slab models presented in literature (Mahajan et al., 2010;Heierli et al., 2008;Bazant et al., 2003;Borstad and McClung, 2011).Accordingly, its behaviour is controlled by Young's modulus, E, and Poisson ratio, υ (Table 2).We use Young's modulus values vary-ing with density after Mellor (1975), in the range of 1.2-1.5 MPa.We follow the study by Teufelsbauer (2011) and select a Poisson's ratio to be equal to 0.04 (for temperature −10 • C).Also we note that since the problem deals with dynamics and vibration, non-physical viscosity of the block, η, is introduced into the damping matrix of the model for numerical stability reasons.Sensitivity tests to Young's modulus, E, Poisson ratio, υ, and viscosity, η, showed that they have a negligible influence on failure properties (see Sect. 5.3).

Constitutive laws of the interface
The interface is governed by a Mohr-Coulomb failure criterion, with an angle of internal friction, φ, and a cohesion, c: where τ is shear stress and σ is normal stress (Fig. 1d).
The cohesion is related to the tensile strength, σ st , as follows (Fig. 1d): Accordingly, we may refer in the following text to both tensile strength (σ st ) and cohesion (c), depending on the context.We note that considering an interface law that includes a tensile strength is crucial to reproduce our experimental tests, since these may involve significant tension stresses.Additionally to failure criterion, for joint elements we also specify values of shear and normal stiffness, K s and K n (Table 2).To the best of our knowledge, there are hardly any experimental data for weak layer elastic properties (Föhn et al., 1998;Jamieson and Schweizer, 2000).After conducting sensitivity tests for different couplings of K s and K n (within the range 10 5 -10 8 N m −3 ) for a full set of experiments, the shear and normal interface stiffnesses were set to 10 8 N m −3 .We found negligible effects of K s and K n on failure, which will be discussed later in Sect.5.3.

Cyclic displacements, inertial loadings and gravity
Before initiating inertial loading for a particular set of parameters, we subject the system to its actual weight.Gravity is applied gradually at a rate of 2.45 g s −1 during 4 s until reaching its 100 % value, in order to avoid numerical instabilities.Furthermore, the Poisson ratio of the block, υ, is set to 0 during this initial phase in order to obtain homogeneous normal stresses within the sample, i.e. without any stress concentrations at the edges.In the next procedural step, the value υ = 0.04 is introduced and the system is allowed to stabilise during 0.4 s.
Next, we reproduce horizontal shaking of the platform by imposing displacements onto the lower boundary.Consistently with the experimental conditions, the system base is subjected to the following cyclic displacement: where 0.0165 is the displacement amplitude in metres.The angular frequency, ω, increases linearly as a function of time (after initial sample preparation): where the coefficient k ω varies between 0.44 and 1.43 s −2 (see below).This angular frequency increase introduces the gradual growth of velocities and accelerations, and thus stresses, with every oscillation (Fig. 2).Since a sample's failure always occurs at an instant when acceleration reaches a peak (caused by a change of the platform's direction of movement), and since the corresponding peak acceleration is known from the experimental measurements, we individually adjusted the coefficient k ω for each test in order to recover the right value of peak acceleration at the instant of failure.An example of k ω adjustment for one test is provided in Fig. 3a and b.Values of k ω obtained for all tests are listed in Table 1.
Here, it is also appropriate to recall the simplified analytical evaluation of the shear force evolution (τ ex ) previously used, in the horizontal case, to estimate weak layer shear strength during experiments (Nakamura et al., 2010;Podolskiy et al., 2010b): where m f is the mass of the upper block, a(t) is block acceleration (second derivative of s(t) with respect to time), and A is the area of the failure plane.This analytical approximation corresponds to a purely static model, since it does not account for dynamic stress inhomogeneities caused by inertia and geometry.Since m f = h s A ρ s , where h s is the height of the block and ρ s is its density, Eq. ( 5) can be rewritten for further comparisons as: In the inclined case, gravity effect should also be taken into account so that: where α is the inclination of the boundary.

Definition of sample failure
We define sample failure as the first instant when all nodes of the interface, N, satisfy the Mohr-Coulomb failure criterion: where N f is a number of failed nodes.This instant is denoted t m .We emphasise here that when the stresses decrease during the loading cycles, the interface nodes recover their initial elastic behaviour and strength, and that no progressive accumulation of damage is considered (see Sect. 3.2).

Search for optimal failure parameters
Systematic numerical simulations were run to find values of cohesion c and friction angle φ that minimise the time difference between instants of failure predicted by the model (t m ) and those measured in the experiments (t e ).By adjusting these two degrees of freedom c and φ, we investigated the ability of the assumed weak layer failure threshold to predict correct failure time values.The adjustment was performed on all tests, involving a variety of experimental conditions (different masses and inclination angles), since all used samples were expected to be characterised by similar weak layer properties.Accordingly, we defined a numerical optimisation search procedure based on the following constrained singleobjective cost function, C FEM : for a number of simulated tests n.We consider a parameter space where cohesion, c, is limited to a range between 0.5 and 2.8 kPa (Föhn et al., 1998), and the internal friction coefficient range is to 0.18-3.73,or 10-75 • (e.g.Keeler and Weeks, 1968;Nakamura et al., 2010); more detailed discussion is provided in Sect.6.In order to reduce computational costs, instead of covering the c-φ parameter space by all possible discrete combi-nations and for the whole ensemble of 19 tests, we followed cost function gradients manually by selecting a smaller representative sample of experiments.This reduced "calibration" sample consisted of 5 (or 9) individual tests selected with different inclinations, masses and sizes to avoid possible biases.The results obtained with the "calibration" sample would then be verified on the complete "validation" sample (as will be explained in Sect.5.2).

Mechanical behaviour of samples and failure
Figure 4 provides examples of stress fields within the blocks caused by motion and the geometry of the system.Two principal observations can be made for all inclinations.First, as the block changes its direction of movement and thus experiences high accelerations, we observe the expected emergence of maximum shear stress (see instant t 2 at Fig. 4).These stresses then decay as the block moves backward and passes through the central position of its trajectory (t 3 ).At the opposite side of the oscillation (t 4 ), shear stresses re-peak with an opposite sign (Fig. 4).Second, we see that at the critical points (t 2 and t 4 ), normal stress remains quasi-constant in the middle of the block, but shows important variations of opposite signs at the edges.Due to the inertia of the mass, one side will experience an increase in normal stress, while the other a decrease.With higher accelerations, these decreasing normal stresses progressively turn into tension.As the block leaves the point t 2 and reaches the opposite critical point (t 4 ), signs of normal pressure reverse.
Similarly, in the interface the imposed oscillations produce shear stresses with changing directions and strong oscillations of normal stress at the edges of the joint (Figs. 5 and 6).Tensile stresses appearing at the edges clearly illustrate that tensile strength of the weak layer needs to be taken into account for realistic representation of tests (Fig. 5). Figure 6 also shows the differences between the analytical (Eqs.limitations of the analytical approach for samples of limited length.For the inclined tests (25 and 35 • ), the differences between the analytically and FE-derived shear stresses in the middle of the interface are slightly smaller (Fig. 6b and c).However, edge effects are more significant for these inclined tests (Figs. 5 and 6).
Figure 7 shows the growth of the number of nodes, N f , that have reached failure criterion with time.As expected, as the block passes the critical point and reverses its direction, the stresses start dropping so that the number of failed nodes N f progressively diminishes and is always zero at the central position of the trajectory.The next peak is then larger than the previous one, because of the progressive acceleration increase.Accordingly, we observe progressive enlargement of the failure zone, purely driven by the external loading, during the successive oscillations (Fig. 7).The time difference between the instant of "total failure" (t m ) and experimental failure (t e ) is also indicated in Fig. 7.

Mohr-Coulomb parameter optimisation
The failure time delay t m − t e obtained for all tests and for different pairs (c, φ) is indicated in Fig. 8.For the considered range of parameters, experimental time-to-failure, t e , is reproduced within 20 % for the majority of tests with only a few outliers.The figure shows that if the modelled joint has a cohesion that is too high, failure will be delayed compared to t e ; on the contrary, if it is too low, failure will occur earlier.
In general, the responses of all individual tests to changes in c and φ appear similar, thus justifying the choice of a smaller sample for adjusting these parameters.
More specifically, Fig. 9 shows cost function, C FEM , sensitivity to a selection of a different number of tests and illustrates that earlier introduced sub-sampling (Sect.4.2) is reasonable for the optimal parameter search.For particular variations in parameters the sample's C FEM5 (where subscript 5 indicates the number of tests considered) responds similarly to C FEM computed for the complete population of tests.In the following, two validations samples will be considered, "#! !respectively (blue corresponds to the middle of the interface; red to the lower edge; green to the upper corresponding to C FEM5 (tests 27, 30, 33, 35, 41) and C FEM9 (idem with also tests 23,26,32,39).
The variations of C FEM5 with cohesion and friction angle, c and φ, are shown in Fig. 10 (and in Table 3).The figure shows all tested combinations of c and φ together with some sensitivity tests.The most important outcome of the parameter optimisation (Fig. 10) is the lack of a clear global minimum in C FEM5 .Instead we can observe a valley, which is narrow in cohesion c, but wide in φ, and which is characterised by very close values of C FEM5 (this is more clearly seen in the colour contours based on a cubic interpolation).Accordingly, simulation results appear more sensitive to the cohesion than to the angle of friction.
Following this finding that simulations with different pairs of c-φ resulted in comparable values of cost function, a similar study was performed on C FEM9 (Fig. 10, Table 3).However, even with additional tests a minimum in φ did not become evident.We could nevertheless identify the three pairs of c-φ corresponding to the lowest computed cost function values: 1.57 kPa-30 • , 1.57 kPa-35 • , 1.6 kPa-30 • with C FEM9 = 0.365, 0.373 and 0.385 s, respectively.
Three "validation" sample simulations (including the remaining tests: 25,31,37,40,42,43) were run with these parameters resulting from the optimisation of the calibration sample.Excluding test 25, which behaved as an outlier for the three cases, the "validation" samples produced similarly low C FEM values to those that were obtained with the "calibration" sample (Table 3; C FEM5v = 0.406, 0.377 and 0.394 s, respectively).For example, for simulations with c = 1.57kPa and φ = 35 • , the time difference between modelled and observed failures correspond, on average, to 5 % of the total duration of each individual test.

Sensitivity analysis
In this section we briefly describe the sensitivity tests that were performed in order to confirm that none of the results provided in this study were affected by other parameters of the model.These tests were performed during different stages of the model development and testing, therefore here we just summarise the main conclusions.The ranges of values used for this sensitivity analysis are specified in Table 2.
Sensitivity tests with a higher Young's modulus, E, of the block have shown negligible increase in the magnitude of stresses within the joint (about 1-2 %), and negligible effects on computed time-to-failure (t m ).Numerical experiments (s6y; for nine tests) with the same cohesion, c, and angle of internal friction φ as in the s6 simulations (Table 3), but with a Young's modulus twice or thrice as high produced similar C FEM9 values (0.383 and 0.380 s compared to 0.373 s in the s6 simulations; see Fig. 10 with details also shown in Table 3).(27, 30, 33, (27, 30, 33, tests (25, 31, 37, 35, 41) 35, 41, 23, 40, 42, 43)  Sensitivity calculations with respect to Poisson's ratio, υ, in the block showed that a value of 0.23 instead of 0.04 produces slightly higher normal stresses within the joint (1.1 % at the largest), and thus may delay the timing of total failure but for only one time step at the most.
Variations in the numerical viscosity of the block, η, by two or four orders of magnitude (from 10 4 Pa s up to 10 6 or 10 8 Pa s or down to 10 2 Pa s) have a negligible effect on failure time.The baseline value of 10 4 Pa s was found to be optimal for the overall stability of the model (no artificial high-frequency stress oscillations).
Thus, in short, none of the parameters tested in this sensitivity analysis have effects on the computed failure time comparable to the impact of the failure criterion parameters c and φ.

Interpretation of observed behaviour
The main objective of the study was to investigate the applicability of the Mohr-Coulomb failure criterion, which is one of the most common failure criteria in mechanics of granular material, in relatively complex experiments performed on sandwich snow samples.A first important result is that FEM modelling appeared necessary to capture the stress inhomogeneities arising within the sample that were disregarded by previous analytical analyses (Podolskiy et al., 2010b).Our results also show that even with a simple set of model assumptions, it was possible to reproduce correct failure times for very different experimental cases (i.e. with various inclinations, masses and sizes).Importantly, it appears that the normal stress dependence of the failure criterion is an important ingredient that should be taken into account.In particular, criteria involving only a cohesion cannot reproduce as well the considered set of experiments.We also recall that the occurrence of normal stress oscillations impose a requirement of the interface to have tensile strength, σ st .This means that the weak layer cannot be described by a purely cohesive form of the Mohr-Coulomb failure criterion.
To illustrate the meaning of the cost function C FEM results indicated in Fig. 10, all corresponding Mohr-Coulomb failure envelopes are represented in Fig. 11.In Fig. 11 we have used green shading and red lines to highlight results for which the C FEM is lower than 0.5 s (for both sample sizes, i.e. with 5 or 9 tests).It clearly appears that these optimal simulations provide strong constraints on the value of cohesion c, which lies in the range 1.6-1.8kPa.These cohesion values derived through our inverse simulations fall well within the range of measurements reported for weak layers composed of precipitation particles or interfaces (Föhn et al., 1998).
Unlike for cohesion, the simulations do not provide strong constraints on the values of friction angle φ (Figs. 10 and 11).Thus the overall behaviour of the observed failures in the considered experiments appears mostly controlled by cohesion, c, while friction angle φ plays only a secondary role (in the range 20 to 60 • ).It is probable that this behaviour is partly due to a limited range of sample heights and inclinations, and thus to insufficient variations of the normal stresses between the different experiments.Slight variability between the tests may further enhance the poor localisation of the minimum in φ.
An additional effect may also contribute to relatively poor resolution in friction angle provided by our numerical optimisation.For a fixed value of cohesion, c, the friction angle, φ, controls both the value of the tensile strength, σ st , and the linear "strengthening" of the interface with higher normal stress (e.g.Fig. 11).Hence, with a higher (resp.lower) friction angle, the tensile strength of the interface becomes smaller (resp.higher) than the cohesion, and at the same time, the compressive part of the criterion has a steeper (resp.lower) inclination and requires higher (resp.lower) shear stress for failure.Due to stress inhomogeneity along the interface the above described dual effects are always superimposed in simulations.Thus, for an instance of high φ, if some edge nodes easily "failed" in tension at a given oscillation, the rest of the nodes will be stronger in compression.This could explain the comparable times of model failure obtained for some tests computed with fixed cohesion, but different values of φ (e.g.Fig. 8).In this respect, we can expect the friction angle to play a stronger role for the inclined tests, in which the tensile component of stress is higher.As shown in Fig. 12, this suggestion is supported by computations of the cost function performed only on inclined tests (26,27,32,33).Values of C FEM display more pronounced variations with φ in this case, with a clear minimum for φ = 30-35 • .This range may be considered as the optimal friction angles resulting from our simulations.Experimental data on the angle of internal friction of snow are very scarce.Previously published values of φ vary strongly depending on the literature source (Fig. 13a and b).Approximately 30 • is commonly reported (Schweizer et al., 2004;Gaume et al., 2012), but the value may range from 5.7 to 57.7 • in experimental data or fracture line analysis (Keeler and Weeks, 1968;Mellor, 1975;McClung and Schaerer, 2006;Jamieson et al., 2001).This wide range probably indicates that further clarification and distinction between different snow types will be necessary.Nevertheless, the value around 30-35 • derived from our analysis appears consistent with these previous experimental data, and we ar- gue that it represents the most physically realistic value for weak layers of the type and density used in our experiments.
Finally, we note that while our results stress the importance of accounting for a normal stress dependence of snow failure criterion, the linear shape assumed with Mohr-Coulomb criterion is just an approximation at this stage.Refinements of the criterion by testing more complex envelope shapes (e.g.Haefeli, 1963) or including a closure in compression remain open issues for future work.

Conclusions
This paper presents a FEM study to simulate snow weaklayer failure under cyclic acceleration loading and to analyse the performance of the Mohr-Coulomb failure criterion.The model is tested by comparison with previous cold-laboratory results for shaking-platform experiments (Podolskiy et al., 2010b).An ensemble of individual experiments is simulated and analysed for overall sensitivity to the adjustment of the constitutive parameters.Based on more than 500 simulations, we found that the Mohr-Coulomb failure criterion for the weak layer is sufficient and adequate for the analysis of the experiments.The best couplings of cohesion and friction angle, c and φ, were found to be [1.6 kPa, 22.5-60 • ].The wide range of φ highlights the fact that the reproduction of experiments is largely controlled by the value of cohesion and has relatively low sensitivity to friction angle (within the limit shown above).Based on values of the cost function for a limited sample of inclined tests (Fig. 12) and on previous experimental evidence, we could suggest that φ around 30-35 • is the most optimal value, which may be further clarified with follow-up studies.In addition, the requirements to consider effects of normal stress on failure and to include the tensile strength of the interface were evident, meaning that a purely cohesive form of the Mohr-Coulomb failure criterion is not applicable.The tensile strength could be limited to a range between 0.9 and 3.8 kPa (Table 3), which is comparable to previously reported results (Gaume, 2012).
The FE results were also compared with the previously used analytical solution (Nakamura et al., 2010;Podolskiy et al., 2010b), which was found to be inadequate for estimating shear stresses along the failure plane; in particular, for cases with an inclination of the platform.Shear stresses produced during the inclined tests (25 or 35 • ) were found to be highly non-homogeneous and thus poorly represented by the analytical approach.Accordingly, the interpretation of experiments through the previously used analytical (or "static") solution is limited.Finally, we are aware that our model with the weak layer representation employed here is only one of many possible approaches, which could have been used to fit the data, and that we confronted the method against only one type of weak layer (composed from precipitation particles) used in previous experiments.Nevertheless, the reasonable results described in this paper suggest that our approach may be further verified and developed (for instance, for non-linear shapes of the failure criterion) and may be also applied to other types of loadings and weak layers.Such work, along with computationally expensive comparison against other failure criteria, could constitute follow-up studies.

Figure 1 .
Figure 1.(a) 2-D geometry of the discussed experiments and (b) an example of corresponding geometry in the finite element model; (c) schematic of the joint element; (d) Mohr-Coulomb failure criterion.

Figure 3 .
Figure 3. Examples of fitting angular frequency by adjusting k ω : (a) k ω = 0.33 s −2 (in black) and 1.43 s −2 (in blue), (b) the same zoomed.The markers indicate an example of observed peak acceleration reached at observed failure time.

Figure 4 .
Figure4provides examples of stress fields within the blocks caused by motion and the geometry of the system.Two principal observations can be made for all inclinations.First, as the block changes its direction of movement and thus experiences high accelerations, we observe the expected emergence of maximum shear stress (see instant t 2 at Fig.4).These stresses then decay as the block moves backward and passes through the central position of its trajectory (t 3 ).At the opposite side of the oscillation (t 4 ), shear stresses re-peak with an opposite sign (Fig.4).Second, we see that at the critical points (t 2 and t 4 ), normal stress remains quasi-constant in the middle of the block, but shows important variations of opposite signs at the edges.Due to the inertia of the mass, one side will experience an increase in normal stress, while the other a decrease.With higher accelerations, these decreasing normal stresses progressively turn into tension.As the block leaves the point t 2 and reaches the opposite critical point (t 4 ), signs of normal pressure reverse.Similarly, in the interface the imposed oscillations produce shear stresses with changing directions and strong oscillations of normal stress at the edges of the joint (Figs.5 and 6).Tensile stresses appearing at the edges clearly illustrate that tensile strength of the weak layer needs to be taken into account for realistic representation of tests (Fig.5).Figure6also shows the differences between the analytical (Eqs.8 and 9) and FEM solutions for shear stresses.In general, the FEM gives larger shear stresses, by about 20 % in the middle of the horizontally inclined joint, thus clearly indicating the

Figure 5 .
Figure 5. Example of evolution of normal stresses in the middle and at edges of the interface (blue corresponds to the middle of the interface; red -to the left edge; green -to the right edge).(a) Horizontal test (Test 23); (b) and (c) -inclined tests (25 and 35 • ; Tests 33 and 27).
For example, three sensitivity simulation sets performed for all 19 tests with the same c = 1.6 kPa and φ = 45 • , but vary-"#!!F ig. 8 Example of delays between observed and modeled failures (t m t e ) for different tests as a function $%%&! of adjustment parameters (

Figure 10 .
Figure 10.Effects of c and φ adjustments on time delay between modelled and experimental failures (C FEM , or RMSE; shown for a sample of five tests by empty circles, for a sample of nine tests by crosses, and for Young's modulus sensitivity tests, s6y and s6yy, by pentagrams).Colour contours are based on cubic interpolation for generalisation of results (C FEM5 ).

Figure 11 .
Figure 11.Illustration of all tested pairs of c and φ as parameters of the Mohr-Coulomb failure criterion (blue dashed lines); red lines (with green shading) indicate the most successful simulations (i.e. when both C FEM , for the representative sample of 5 or 9 tests, are ≤ 0.5 s).

Figure 12 .
Figure 12.Effect of the angle of friction, φ, on C FEM for simulations with the same cohesion 1.57 kPa (shown for a sample of 5 tests by blue empty circles, for a sample of 9 tests by red crosses, for a sample of 4 inclined tests by black diamonds).

Figure 13 .
Figure 13.(a, b) Values of the angle of friction obtained from different studies.Y axis in a corresponds to tan φ (it is equal to c/σ st and shown as a blue curve); it visualises the ratio between cohesion and tensile strength).

Table 1 .
Podolskiy et al. (2010b)or validation of the model, afterPodolskiy et al. (2010b)and prescribed modelling parameters for each test.

Table 2 .
Properties of FEM model (values in square brackets correspond to sensitivity tests).