Modeling snow slab avalanches caused by weak layer failure – Part I: Slabs on compliant and collapsible weak layers

Dry-snow slab avalanche release is preceded by a fracture process within the snowpack. Recognizing weak layer collapse as an integral part of the fracture process is crucial and explains phenomena such as whumpf sounds and remote triggering of avalanches from low-angle terrain. In this two-part work we propose a novel closed-form analytical model for a snowpack under skier loading and a mixed-mode failure criterion for :: the nucleation of weak layer failure. In the first part of this two-part series we introduce a closed-form analytical model of a snowpack accounting for the 5 deformable layer. Despite the importance of persistent weak layers for slab avalanche release, no simple analytical model accounting for ::::::::: considering : weak layer deformations is available. The proposed model provides deformations of the snow slab, weak layer stresses and energy release rates of cracks within the weak layer. It generally applies to skier-loaded slopes as well as stability tests such as the propagation saw test. A validation with a numerical reference model shows very good agreement of the stress and energy release rates ::: rate results in several parametric studies including analyses of the bridging effect and slope 10 angle dependence. The proposed model is used to analyze 93 propagation saw tests. Computed weak layer fracture toughness values are physically meaningful and in excellent agreement with finite element analyses. In the second part of the series ::::::::::::::::::::::::::::: (Rosendahl and Weißgraeber, 2019) we make use of the present mechanical model to establish a novel failure criterion crack nucleation in weak layers. The code used for the analyses in both parts is publicly available ::::: under https://github.com/2phi/weac. 15


Introduction
Dry-snow slab avalanches can release when a persistent weak layer of, e.g., surface hoar or depth hoar breaks (see the wellknown image of a partially collapsed weak layer by Jamieson and Schweizer (2000) shown in Figure 1). Weak layer failure can be triggered by additional loads like a skier. If the conditions allow for crack propagation, ::: i.e., :: if ::: the :::::: energy :::::: release :::: rate :: of : a ::::::: growing ::::: crack ::::::: suffices, a triggered initial defect may extend across slopes and eventually cause the slab to fail and slide. 20 The earliest approaches to snowpack stability were so-called stability indices. They consider snowpack loading owing to the weight of the snow slab and owing to additional loading by a skier (Perla, 1977;Föhn, 1987). To account for snow stratification Figure 1. Weak layer of buried surface hoar. Left of the vertical slab fracture the weak layer has collapsed, whereas on the right hand side the porous weak layer is still intact. Reprinted from the Journal of Glaciology (Jamieson and Schweizer, 2000) with permission of the International Glaciological Society.
improved stability index models were proposed by Habermann et al. (2008) or Monti et al. (2015). However, these local models are insufficient to describe the stability of snowpacks across slopes (Bellaire and Schweizer, 2011). Many researchers suggested that these ::::::: stability indices are incomplete as they do not account for the propagation of the failure within the snowpack (van Herwijnen and Jamieson, 2007).
To evaluate stability indices or fracture mechanics criteria, a model of the stress distribution within the snowpack and 5 especially the weak layer is needed. Using the exact solution for a concentrated load on a homogeneous semi-infinite plate, Föhn (1987) proposed a model for the stress distribution below a skier. A first interface model with shear for snowpacks with weak layers was proposed by McClung (1979). Linear elastic finite element analyses have shown that the effect of layering can play an important role (Schweizer, 1993;Habermann et al., 2008). This was accounted for in the model of Monti et al. (2015) by using an equivalent homogenized snow layer with effective uniform Young's modulus and a modified total slab leading to smaller stresses in the weak layer. However, the effect of weak layer compliance in normal direction is unaccounted for in the aforementioned beam models (Heierli and Zaiser, 2008). The snowpack is assumed to rest on a rigid weak layer and only slab deformation in the region of a ::: the collapsed weak layer is modeled. Yet, it is known that the deformation of the weak layer is crucial for deformations and the local load transfer in the snowpack (Reiweger and Schweizer, 2010). (Chiaia et al., 2008;Gaume et al., 2013) :::::::::::::::: Chiaia et al. (2008) :: and :::::::::::::::::: Gaume et al. (2013) consider weak layer deformability only 5 in shear.
In the first part of this two-part contribution, a novel modeling approach for the description of weak layer failure is given.
It aims at providing a model that fully accounts for the weak layer's effect on deformations and load transfer and solutions to the mixed-mode energy release rates of cracks within the weak layer. The model is validated using finite element analyses and field experiments. In the second part we propose a new failure criterion , which ::: that physically links stability indices and 10 fracture mechanics models. Here, the necessary distinction of stress-strength and fracture mechanics approaches is highlighted and discussed.

Modeling approach
Deformation, stresses and consequently the energy release rate of cracks within the weak layer are controlled by loading and the complete stratigraphy of the snowpack. Deformations of the slab and in particular of the weak layer must be rendered 15 sufficiently accurate. The slab is loaded in local bending and stretching, : which we account for using beam and rod kinematics.
As in the analysis by Heierli et al. (2008) we model the slab as an elastic beam with bending and shear deformation. In order to account for weak layer deformations, the present model rests the beam on an elastic foundation of an infinite set of springs with compressive and shear stiffness commonly referred to as Winkler foundation. Base layers are assumed rigid. The model provides slab deformations, weak layer normal and shear stresses as well as the energy release rate of cracks. Skier penetration 20 is not considered. Like other models given in literature and discussed above we employ linear elasticity and neglect, e.g., viscous or plastic deformations. Consequential limitations are discussed in Section 4.

Governing equations
Consider the snowpack model on an inclined slope of angle ϕ depicted in Figure 2a. The plane strain beam representing the snow slab has an out-of-plane thickness b, height h and length l. The weak layer thickness is denoted t. The slab is assumed 25 to be homogeneous with Young's modulus E, Poisson's ratio ν and density ρ. The snow slab is loaded by its own weight prescribed as distributed loads q t = ρghb sin(ϕ) in tangential x-direction and q n = ρghb cos(ϕ) in normal z-direction where g is the acceleration of gravity. Additional loading (e.g., by the weight force of ::: the mass m of a skier, snowmobile or others) is represented by a concentrated force F = mgb/l o where l o is the effective out-of-plane length of the object such as the length of skis. F is split into normal and tangential concentrated forces F n = F cos(ϕ) and F t = F sin(ϕ), respectively. The weak layer of thickness t consists of an infinite set of springs with compressive stiffness k n = E weak b/t and shear stiffness k t = G weak b/t with the shear modulus G weak = E weak /(2(1 + ν)). Cracks within the weak layer are modeled by removing the support of the beam on the crack length a.
Timoshenko kinematics for the beam allow for shear deformation. Initially plane beam cross sections may rotate by an angle ψ (see Figure 2b) yet remain plane during deformation, which yields the general deformation kinematics of the beam where u is the horizontal displacement and w the vertical deflection of the beam middle plane, respectively. The index z is introduced to distinguish between midplane deformations u(x) and w(x) and the actual displacement fields u z (x, z) and Enforcing equilibrium of forces and moments and using the laws of elasticity allows for deriving a set of ordinary differential 10 equations (ODEs) with constant coefficients which ::: that describe the deformation of the snow slab (see Appendix A). The horizontal displacement u is obtained from where A = hb is the snow slab cross section. The vertical beam deflection w and the rotation of the beam cross section ψ are described by where κ = 5/6 is the shear correction factor for rectangular cross sections and I = bh 3 /12 the moment of inertia with respect to the y-axis.
The general solution of ODE (2) for the horizontal displacement is given by with the eigenvalue The free constants c 1 and c 2 must be determined from boundary conditions. The solution of the coupled ODEs (3) and (4) is of 5 exponential type as well. Depending on the material parameters, the eigenvalues of this solution may become real or complex.
When k n EI ≥ 4(κGA) 2 the eigenvalues are real numbers and the general solution is given by If k n EI < 4(κGA) 2 the eigenvalues are complex and read In this case the general solution of the normal deflection is given by Again, the constants c 3 to c 6 must be determined from boundary conditions.
For regions without elastic foundation, e.g., above a crack or when modeling propagation saw tests (PST), we obtain k t = k n = 0 and the differential equations simplify to The corresponding solutions are given by w(x) = c 9 + c 10 x + c 11 x 2 + c 12 x 3 + q n 24EI x 4 .
As before, the free constants must be determined from boundary conditions which are defined by the final assembly of the solution.

Assembling the solution
The present model readily applies to an inclined skier-loaded snowpack with or without crack in the weak layer as well as 5 to propagation saw tests. Therefore ::: For ::: this ::::::: purpose, considered snowpack configurations must be assembled from general solutions of beam segments with or without elastic foundation given above. The free constants are determined from boundary and transmission conditions.
To model a propagation saw test two beam segments are required as shown in Figure 3a. The left part of the snow slab rests on an elastic foundation representing the intact weak layer. The right part is a cantilever beam where cutting has removed the 10 weak layer support. The slab is loaded by its own weight q n and q t . Free left and right ends require vanishing section forces and moments. Using the expressions for the bending moment M , the lateral shear force Q and the normal force N yields 15 at free ends. At the transition between the unsupported and supported segments of displacements, cross section rotation, section forces and moments must be C 0 -continuous. This yields u l = u r , w l = w r , ϕ l = ϕ r , where the indices 'l' and 'r' denote quantities left and right of the discontinuity, respectively.

20
In skier loaded snowpacks (Figs. 3b and 3c) the skier point load adds discontinuities of normal and transverse shear forces and the solution is assembled from two or four regions. Skier loaded slopes extend far beyond the influence zone of local skier loading. Hence, we have chosen a length of 25 times the thickness of the slab (l = 25h) to avoid any edge effects. Of course, also shorter slabs can be modeled to account for the resulting edge effects on the mechanical response if necessary. Again, at either end of the assembly boundary conditions of vanishing section forces hold. At crack tips the continuity conditions Eq. (18) hold. At the concentrated skier load the transmission conditions must be modified according to The boundary and transmission conditions for the respective load case provide a linear system of equations with up to 24 unknown constants. The system can be solved easily using any mathematical toolbox. Closed-form solutions for u, w and ψ 5 can then be given in piecewise form.

Weak layer stresses
Since the weak layer is represented as an elastic Winkler foundation, also known as weak interface model (Lenci, 2001), it is not modeled as a complete continuum. Instead, it only transfers compressive and lateral shear loads. In this simplified continuum, the stress interaction is simplified and compression stress and shear stress :::::::::: compressive ::: and ::::: shear :::::: stresses : are uncoupled. Using 10 calculated slab displacements w in normal and u in tangential direction, : stresses are obtained from Weak interface approaches are well established in structural analysis (Krenk, 1992;Selvadurai, 1979), in particular when strong elastic contrasts are present. Stress solutions are generally of good quality. However, highly localized effects such as stress singularities cannot be captured. The effect of this limitation is discussed in part II of this work.

Fracture mechanical quantities
Fracture mechanics (Broberg, 1999;Anderson, 2017) is concerned with the behavior of cracks in continua. Propagation of existing cracks can be described in terms of energy considerations which allow for an assessment of the stability of cracks.
If the change in total potential energy dΠ with an infinitesimal crack advance da (denoted as differential energy release rate) equals the fracture toughness the crack will grow. This fundamental condition is called Griffith criterion. The fracture toughness describes the energy required for formation of a new crack surface of unit area. It comprises surface energy as well as dissipative energy terms. The latter render fracture processes irreversible. The fracture toughness is defined within the continuum framework. That is, it efficiently covers all local failure mechanisms on the micro-scale -even for complex and heterogeneous materials. 1

25
Within the framework of weak interface models the energy release rate corresponding to infinitesimal crack growth can be obtained from the local strain energy (cf. Krenk, 1992). For mode I and mode II contributions to the energy release rate hold. Here, w(a) and u(a) correspond to displacements at the crack tip. Again, this simplified framework typically provides good results but it does not capture certain effects such as the energy release rate of vanishing crack lengths (cf. discussion in section 4). 1 In linear elastic materials the Griffith criterion may be reformulated using stress intensity factors K which describe the asymptotic crack tip stress field.
Both, the critical energy release rate Gc and the critical stress intensity factor Kc are denoted fracture toughness. The units of these two quantities are different but can be converted using the material's stiffness :::: elastic :::::: modulus. Figure 4. Graphical representation of the mode I crack opening integral. An uncracked snowpack can be represented by removing a part of the weak layer and applying virtual stresses to both the slab and the removed weak layer segment such that they are deformed as in the original configuration (w0). Reducing the virtual stresses to zero quasi-statically increases the beam deformation (w b 1 ), relaxes the weak layer (w k 1 ) and finally yields the cracked configuration. The work done by the virtual stresses (integral of σ∆w on ∆a) corresponds to the total potential energy difference between uncracked and cracked configurations.
We distinguish different crack opening modes. Mode I loading is a crack opening mode normal to the crack faces. Strictly 5 speaking this comprises only symmetric deformations which typically does not hold for cracks along material discontinuities with different stiffnesses. Mode II and III are shear crack modes with tangential displacements of the two crack faces. In this work only modes I and II are considered. Following the concept of anticracks (Heierli and Zaiser, 2008) we extend the mode I definition to general deformation normal to the crack faces. This includes tearing deformations of the crack faces away from each other but also collapse where crack faces deforming towards each other. Since the micro-structural failure mechanisms 10 are of course very different for tearing and collapse the magnitude of the associated fracture toughnesses will also be different.
Hence, we must distinguish between mode I fracture toughness for tearing G + Ic and mode I fracture toughness for collapse G − Ic . Loading which causes any combination of the three crack opening modes are :: is called mixed-mode loading. Then mixedmode failure criteria must be used to account for the different contributions to the total energy release and the required energy to grow a crack under mixed-mode conditions. This will be subject in part II of this work.

15
In order to extend the scope of fracture mechanics, which is only applicable to infinitesimal crack growth da, Hashin (1996) considers finite differences in crack length ∆a. Then the so-called incremental energy release rate describes the difference in total potential energy ∆Π for a finite crack increment ∆a. The incremental energy release rate can also be obtained integrating the differential energy release rate over the finite crack advance ∆a The total differential energy release rate is the sum of contributions from different crack opening modes: Equivalently, the total incremental energy release rate can be split into contributions from crack opening (mode I) and crack sliding (mode II): The total potential energy difference ∆Π = ∆Π I +∆Π II between a cracked and an uncracked configuration can be calculated efficiently using the crack opening integral. The change in potential energy corresponds to the work done by stresses on crack flanks when reduced quasi-statically to zero as detailed in Figure 4. The mode I and II contributions amount to Using the stress to displacement relation (20) obtained in :::: from : the snowpack model and Eq. (23) yields The present model provides slab displacements, weak layer stresses and energy release rates for cracks within the weak layer as closed-form analytical expressions. In order to validate the model, stresses and energy release rates are compared against detailed finite element analyses (FEA) and existing models. Results of several parametric studies are shown in detail. Further, we compute the fracture toughness corresponding to critical cut lengths in propagation saw tests for a comprehensive set of 93 field experiments provided by Gaume et al. (2017) to investigate the capabilities of the present model.

Reference solution
Reference stress solutions and energy release rates are computed using the plane strain FEA model shown in Figure 5. Like the reference models used by Sigrist and Schweizer (2007) and Habermann et al. (2008) it considers an inclined snowpack  consisting of a homogeneous slab and a weak layer. The slab is loaded with a vertical volume gravity load. The weak layer is clamped at its bottom side. Cracks are introduced by removing all weak layer elements on the crack length a. The mesh of biquadratic 8-node elements with reduced integration is refined towards stress concentrations. Mesh convergence of the FEA 5 solutions has been controlled. The FEA total energy release rate of a crack of length ∆a is computed according to where the incremental energy release rate G(∆a) is determined using Eq. (23). In order to calculate the derivative of G(∆a), the incremental energy release rate G is evaluated for four different crack lengths closely around ∆a. The derivative G(∆a)/∂∆a is then obtained by differentiating the interpolating cubic spline of the four G values at ∆a.  Table 1Table 1.

5
Considering shear stress, Föhn provides an exact solution for transverse shear stresses in a homogeneous body. Theses ::::: These stresses are zero directly below the concentrated force and peak to the left and right of it as seen in Figure 6. Additional shear stresses arise from a horizontal displacement of the beam which also strains the weak layer. Lateral shear stresses originating from this effect peak below the load point. Superimposing both components yields the total weak layer shear stress obtained in FEAs. Note that lateral shear stresses owing to slope-parallel concentrated force components do not change their sign left 10 and right of the load point as shown in Figure 6. This is in contrast to transverse shear stresses owing to normal concentrated forceswhich : , :::::: which :: do ::::::: change their sign. However, as discussed in section 4 transverse weak layer shear stresses are not accounted for in the present model. It considers only lateral shear where its agreement with FEA results is better than Föhn's transverse shear solution.
The bridging effect of stiff slabs is studied in Figure 7. The results show that with thicker slabs the local loading of a skier is 15 transferred on a wider area of the weak layer. This dependence on the slab thickness is very pronounced as the bending stiffness of the slab increases in cubic dependence with the slab thickness 2 . Because the weight of the slab is proportional to the slab thickness, the stress level distant from the skier increases.  Table 1. layer allows for point load distribution over a larger area. In comparison to FEA data, a very good agreement of the present 5 model except is evident for either thickness. Very local deviations just below the point load are observed. The Föhn (1987) model is shown for comparison. As it does not account for layering, its normal stress distribution is independent of ::: the : weak layer thickness.
Let us consider the energy release rate solution. Figure 9 shows the incremental ::::::::: differential mode I and II energy release rates as a function of slope angle. Besides the implemented Timoshenko kinematics, the limit case of Euler-Bernoulli beam theory 10 is shown. In the numerical reference model only the total incremental energy release rate is evaluated. The comparison shows that for moderate slope angles the present model provides an excellent prediction of the energy release rate. While the mode I contribution of the energy release rate decreases for higher slope angles, the mode II contribution increases monotonously.
Using Euler-Bernoulli beam theory underestimates the mode I contribution and hence the total energy release rate. The mode II energy release rate is a result of tangential displacements and thus unaffected by the choice of beam kinematics.

15
In Figure 10 a PST experiment is studied. The mode I and mode II energy release rates are shown in dependence of the thickness of the slab above the weak layer. Two different slope inclinations are considered. Energy release rate results are normalized to account for the different orders of magnitude of both crack opening modes. For both angles the mode II contribution has a stronger dependence on the slab thickness than the mode I contribution as they are governed by different deformation mechanisms. Hence, the mode-mixity changes towards a more mode II dominated loading of the PST. No mode III contribution exists when PSTs are cut downsloping and no lateral loading occurs. other ::: than ::::::::::: E weak = 0.15 MPa ::: and ::::: t = 5 cm are chosen as given in Table 1. other ::: than ::::::::::::::: E weak = 0.15 MPa ::: and :::::: t = 5 cm : are chosen as given in Table 1.  Table 1.
20 Figure 14 correlates model predictions for the weak layer fracture toughness to data obtained using detailed FEAs and Eq. (28). The FE model features slope-normal faces such that the corresponding section force boundary conditions discussed above are applied to the beam model. Fracture toughnesses G − c are determined from critical cut lengths a c measured in 93 PST field experiments (Gaume et al., 2017). Measured cut lengths correspond to critical lengths required for crack propagation. No case in the data set showed maximum deflections exceeding the weak layer thickness which would indicate a base-touching 25 slab. According to the fundamental Griffith criterion of fracture mechanics, Eq. (21), the energy release rate corresponding to the critical cut length is the fracture toughness of the weak layer. We will denote it G − c to emphasize the difference between the fracture toughness of a collapsing weak layer G − c and a tearing fracture toughness G + c . The difference between the two will be discussed in section 4. Please note, that here only the total energy release rate is considered and no mode-mixity dependence is assumed. The 93 PST measurements allow for a model validation against the reference solution and we can obtain actual weak 30 layer fracture toughnesses. For comparison, : the closed-form analytical expression of Heierli (2008) given by Schweizer et al. (2011) is shown. In the notation of the present work this expression reads with w 1 to w 4 given in Appendix B.
Predictions of the present model are found within a narrow range around the one-to-one line indicating excellent agreement (R 2 = 0.99) with FEAs for this comprehensive set of real world parameters. However, the ::: The : Heierli solution significantly y = .
x + .  Gaume et al. (2017) are assumed as given in Table 1.  Table 1. underestimates the energy release rates as already seen in previous analysis of the influence of the weak layer thickness. The comparison to the FEA reference solution shows poor correlation (R 2 = 0.04). an increased mean relative error amounting to 12.7 :::: 19.3%. The relative error of the predicted 10 fracture toughness and the numerical reference model is between −27.9% and +45.0 ::::::: −28.3% ::: and :::::: +68.3%. 50% of the values have an error of less than 9.2 ::: 9.0%. The maximum value of the fracture toughness 3 calculated is 204.4 ::::: 124.0 Jm −2 with a relative error of −3.26 :::: 2.48% in the comparison to the reference model.
The previously given results of the energy release rate as given by the model show differential energy release rates G describing the growth of initial cracks introduced by the saw. In order to study crack initiation, : it is important to study ::::::: consider 15 the incremental energy release rate G of cracks of finite size, see Eq. (23). In Figure 15 the incremental energy release rate of finite cracks in the weak layer of a skier loaded snowpack are : is : shown. Besides the local load also the weight load of the slab is considered. As an example the effects of slab and weak layer thicknesses are studied. The incremental energy release rate is given as a function of the size of the finite cracks. Of course, the incremental energy release rate increases with crack length. It also increases with slab thickness and with increasing weak layer thickness.

Discussion
The presented closed-form analytical model of the snowpack contains two different levels of abstraction. The first level is the treatment of the snowpack as a linear-elastic continuum. This is common for models of skier-triggered avalanches (Schweizer and Camponovo, 2001; and agrees with the fact that the brittle failure process of dry-snow slab avalanches occurs within very short time scales (Narita, 1980). Further, (van Herwijnen et al., 2016) ::::::::::::::::::::::: van Herwijnen et al. (2016) 10 have shown a good agreement of displacements measured using particle tracking in field measurements with linear elastic models.
If it was :: In :::: order : to study long-term effects within a snowpack like temporal compaction or sintering/metamorphic processes such models are insufficient (Capelli et al., 2018;Mulak and Gaume, 2019). At first glance, the propagation saw test seems governed by the rather slow saw movement. However, the fracture process itself is associated to high strain rates (Reiweger 15 and Schweizer, 2010).
The second layer of abstraction is the simplification of heterogeneous media to continua and engineering structures like elastic foundations or deformable rods and beams. If chosen correctly, accurate representations of the deformation and stress fields within the continua can be recovered with such elements of structural analysis which are well established in civil and mechanical engineering Gross et al. (2014) ::::::::::::::: (Gross et al., 2014). Using appropriate simplifications, : numerical models can be 20 avoided and closed-form analytical solutions are obtained. Comparisons to finite element analyses (FEA) (see Figure 9) with continuum elements show that first order shear deformation theory is needed to account for the low shear stiffness of the slab.
Euler-Bernoulli beam theory does not suffice and the present work uses Timoshenko beam theory.
The studies shown in Figures 6 and 8 asses the quality of the stress solution. The present model shows very good agreement with FEA reference solutions with relative errors below 2-5% outside the stress localization domain just below the point 25 load. We exclude this domain from our discussion because considering a localized distributed load instead would resolve this typical problem of concentrated loads and weak interface models. Moreover, the effect of the point load limitation on the results of failure models is discussed in detail in part II. Stiffer weak layers yield smaller errors. For thicker weak layers stress concentrations below a skier load are less pronounced with larger areas of load transfer (cf. Figure 8). Hence, the peak stresses in the weak layer are reduced with increasing thickness of that layer. This result of the stress solution is in line with the FEA 30 reference solution but cannot be captured by the Föhn (1987) solution, : which makes use of the elastic half-space solution.
This model deficiency is discussed in detail by Gaume et al. (2017). However, the effect of reduced local load is in contrast with ::::::::: contradicts : observations that thicker weak layers are more likely to fracture and that avalanches are easier to trigger when the weak layer under the slab is thicker (Gaume et al., 2013). The failure analysis presented in part II of this work is able to correctly recover this effect by additionally considering the energy release rate as a requirement for failure. Figure 7 investigates the bridging effect. This effect describes the ability to distribute external loads on the slab on a wider area of the weak layer (Thumlert and Jamieson, 2014) as the slab thickness increases. The present model is able to capture the overall stress level due to the increased slab weight as well as the load distribution on a wider area in very good agreement with the numerical reference solution. The load transfer area increases by a factor of approximately five for a change of slab thickness from h = 10 cm to h = 80 cm. Although not shown, the classical solution by Föhn (1987) is capable of rendering this bridging effect correctly.  have introduced a bridging index as the product of slab hardness and thickness.
The results of the present analysis show that this is a good estimate but the local stresses decrease overproportionally with the slab thickness. This is due to the bending deformation which is controlled by the bending stiffness of the slab, which in turn is proportional to h 3 . The effect of the elastic contrast of the stiffnesses of weak layer and slab on this bridging behavior has been addressed by Monti et al. (2015). They proposed to scale the effective slab thickness with the cubic root of the elastic contrast of slab and weak layer. The bridging effect is very important and explains why thicker slabs can sustain higher loads before weak layer failure occurs (Schweizer and Camponovo, 2001). However, propagation of cracks is more likely below thicker slabs (van Herwijnen and Jamieson, 2007). This change of parameter effect is addressed in the failure model presented in the second part of this work.

15
The energy release rate of cracks within the weak layer is studied in detail in the validation studies shown in Figures 9 to 14. The effect of slope angle (Figure 9) on the energy release rates shows that both mode I and mode II contributions strongly depend on the slope angle. While the angle dependence of :: the : stress solution mainly stems from the split in normal and tangential force contributions, the angle dependence of the energy release rate also depends on the deformations of the slab above the cracked weak layer (cf. Figure 4). Therefore, the angle dependence of simplified models like McClung ( which is carefully set up to identify these important fracture mechanical quantities. The effect of the slab thickness on the energy release rate obtained for PST experiments is shown in Figure 10. As mode I and mode II are associated with very different failure modes of collapse and shear failure of the weak layer, the corresponding energy release rates can be of different magnitudes and at the same time both be relevant for the mixed-mode failure. This has 25 been studied in the comprehensive experimental work of Birkeland et al. (2019). They performed PST experiments in which the slab thickness was changed in order to assess the change of failure behavior at hand of the measured critical cut lengths.
The results of the present model show different ratios of mode I and II contributions when the slab thickness changes. This may allow for identifying a comprehensive mixed-mode fracture envelope using these extended PST experiments. The present model could also help to identify relevant geometric parameter ranges to be studied within experimental campaign ::::::::: campaigns.

30
As pointed out by Reiweger and Schweizer (2010) up to 90% of the total deformation of snowpacks is concentrated in the weak layer because of its considerably low stiffness. It is evident that in order to capture snowpack deformations adequately, a mechanical model must consider the weak layer. Increasing the thickness of the weak layer while keeping all other geometry and material properties constant increases its compliance. That is, Heierli's assumption of rigidity leads to increasing discrepancies as the weak layer thickness increases. This is reflected in Figure 11. The importance of considering the compliance of 35 the weak layer below the slab is also evident in Figure 14. The thickness of the weak layer (which influences its compliance linearly) has a pronounced effect on the energy release rate obtained for these cases. The model by Heierli and Zaiser (2008) :::::::::::::::::::: Heierli and Zaiser (2008)  (2016) also showed the insufficient approximation for longer cracks and resorted to use :::: using a numerical (FEM) correction factor. The energy release rate predicted by the present model, : which accounts for weak layer deformation : , is in good agreement with the numerical reference. Figure 9 indicates that with higher slope angles above 35°the agreement of present solution and reference solution would decrease by a few percent. However, the overall effect of the weak layer thickness (Figure 11) remains the same.
The analysis of 93 PST data points ( Figure 14) provides insight into the analysis of the energy release rate. The 93 data points 10 show a very wide span of the input parameters slab density (factor 4.6), slab thickness (factor 3.9), weak layer thickness (factor 44) and slope angle (from 0°to 44°). They provide a realistic overview on possible configurations along with corresponding measurements of the critical cut length. Applying the model to these data points now not only provides a comprehensive overview on the models capabilities by comparing it to the numerical reference model but also provides insight into the fracture toughness of a multitude of weak layers. The comparison to the numerical reference model shows that the present 15 model provides a good prediction of the energy release rate for most of the configurations. Only eight (8.6%of the data points :: 10 ::::::: (10.75%) data points deviate by more than 25 :: 30% from the numerical reference. Since the energy release rate has a quadratic dependence on (skier) loading : , this constitutes an error in a crack propagation analyses below 12 :: 14%. Because Heierli's model neglects weak layer deformations, it underestimates the fracture toughness significantly (see Figure 11) In the framework of continuum mechanics the PST must be considered a fracture mechanical experiment aimed at identifying the fracture toughness. The fracture toughness is the energy required to form a new surface on the idealized plane of fracture.
This energy comprises the dissipative processes at the micro-scale. These processes are much different between crack growth 30 under tension and under compression with a collapse of the weak layer. In the latter case local dissipative damage processes are much more pronounced, leading to a significantly higher value of the (effective) fracture toughness in compression G − c than in tension G + c . This has been in observed in lab experiments of cracks on glass foams leading to critical energy release rates two orders of magnitude higher in compression than in tension Heierli et al. (2012) :::::::::::::::: (Heierli et al., 2012) 4 . For such porous materials with thin micro-structure : , local stability failure also contributes significantly to the local damage processes :::::: process. Sigrist (2006) and McClung (2007) show that typical values of the fracture toughness in tension are between 10 −1 Jm −2 and 10 0 Jm −2 , which provides the same relation of tensile and compression fracture toughness as in the experiments by Heierli et al. (2012). The tensile fracture toughness of ice (G + c ) is at the order of 10 0 Jm −2 (Dempsey, 1991). This is, of course, larger than the fracture toughness of snow in tension. However, the fracture toughness for the present case of mode I in compression may lie well above this value as the corresponding failure process of collapse dissipates significantly more energy than a simple tensile bond breaking. Compared to temperature and time-driven transformation processes of snow, PSTs represent rather fast experiments with little or no impact of viscous effect ::::: effects. Hence, material properties determined from PSTs may be used for analyses of skier-triggered snowpack instability, : which is associated to fast loading and high strain rates.
As pointed out by Gaume et al. (2017), fracture parameters obtained from experimental results are always linked to the failure mode that is considered in the model employed to evaluate the data. Hence, fracture energies obtained under the assumption of a rigid weak layer (Heierli et al., 2010) or in pure shear models (McClung, 1979) will always be different from models that take into account weak deformation and mixed-mode fracture. If the failure mode is modeled insufficiently, the obtained failure parameters are not true material parameters and hence not consistent.
15 Figure 15 shows the result of incremental energy release rate of cracks of finite size (as opposed to infinitesimal crack growth). These effects are very similar to those the previously discussed case of the differential energy release rate. When differential energy release rates can be calculated as a function of the crack length, incremental energy release rates can always be obtained from integration. The integration is highly efficient using closed-form analytical solutions and established (numerical) integration schemes. However, computational effort can be reduced even further employing the crack opening 20 integral. While integrating the differential energy release rate requires its evaluation for several crack lengths, the crack opening integral only needs data from the uncracked configuration and the configuration with the final crack length. It computes the difference in total potential energy between an uncracked state and a configuration with finite crack from the work done by stresses as the crack opens. Similar to differential energy release rates, incremental energy release rates increase with increasing crack length ( Figure 15). More importantly, they also increase with both weak layer thickness and slab thickness. Maximum 25 deflections of the slabs considered in Figure 15 are smaller than half the weak layer thickness even for the longest cracks.
At first glance, Figures 7 and 15 demonstrate a contradiction: Thicker slabs distribute loads more evenly over the weak layer, reduce local stress peaks and should be more difficult to trigger. However, they release more energy favoring crack propagation.
A failure criterion using both stress and incremental energy release rates to asses the nucleation of finite cracks within the weak layer resolving this contradiction will be introduced in part II of this work.

30
The present model does not account for layering of the slab above the weak layer. It is well known that the layering can have a significant effect (e.g., Schweizer, 1993;van Herwijnen and Jamieson, 2007;Habermann et al., 2008;Reuter et al., 2015).
It influences the load distribution, the fracture process of the weak layer and also the bridging effect. Monti et al. (2015) have proposed a stability index based on the Föhn solution that accounts for layering by using adapted values of the thickness and the stiffness of the slab. In the present work, rigid base layers are assumed. However, as shown by Jones et al. (2006) and Habermann et al. (2008) substratum elasticity is not negligible. The elastic foundation should not only account for weak layer compliance but also for deformable base layers. Future research is required to make use of these considerations in a refined version of the present model to combine the capabilities.
Penetration depths leading to forces acting not on top of the slab but well within the slab are not considered at this point.
To understand loading scenarios where the penetration depth is high, i.e. :::: e.g., because of a soft top layer of snow or a highly localized load (recreationist by foot or a snowmobile) further modeling and experimental effort is needed.
In this work, interaction of shear and compression is limited. Slab extension causes only tangential and slab deflection causes only normal weak layer deformations. This modeling strategy can be pictured as weak layer springs attached to the mid-surface 10 of the beam. With this simplification the governing system of equations is uncoupled and can be solved independently for tangential and normal displacements, respectively. Of course, the weak layer does not interact with the mid-surface of the slab but with its bottom side. Tangential displacements on the slab bottom side are caused by both slab extension and the rotation of beam cross sections, : which in turn depends on the slab deflection. Hence, the rotation of beam cross sections couples tangential and normal displacements. This will increase the solution effort but simultaneously provide improved accuracy, in particular 15 concerning weak layer shear stresses and mode II energy release rates (see, e.g., Figure 6).

Conclusions
By considering a deformable weak layer the present work provides a simple but comprehensive closed-form analytical model for snowpack deformations, weak layer stresses and energy release rates of cracks within weak layers: 1. The model applies to skier loaded slopes as well as PST experiments. 20 2. Providing closed-form solutions, : the present analysis framework is highly efficient and evaluates in real-time.
3. Comparisons with numerical reference solutions indicate good agreement of weak layer stress and energy release rates. 4. The model renders physical effects such as the bridging effect of thick slabs, the influence of slope inclination and most importantly the impact of weak layer compliance. Evaluating a comprehensive set of 93 PSTs fracture toughnesses of a multitude of weak layers are analyzed and discussed. 25 5. Providing weak layer stress and energy release rates of cracks within the weak layer, : the present framework allows for evaluating comprehensive criteria for skier-triggered crack nucleation and crack propagation in the weak layer which will be discussed in part II of this two part work.

Appendix A: Derivation of governing equations
The set of differential equations governing the present problem can be derived from the equilibrium of forces and moments of a beam element and elasticity laws for the bending moment M , the transverse shear force Q and the normal force N . The latter read (Timoshenko and Goodier, 1951) where A = hb is the snow slab cross section, κ = 5/6 is the shear correction factor for rectangular cross sections and I = bh 3 /12 the moment of inertia with respect to the y-axis. In the limit κGA → ∞ we obtain the classical Euler-Bernoulli beam theory. Equilibrium of forces and moments of a beam element yields where k t u and k n w are distributed reaction forces acting on the beam, : which originate from the elastic foundation. Therefore, they depend on the normal and tangential displacements along the the top of the weak layer and its normal and lateral shear 15 stiffnesses k n and k t .
Preliminary studies showed that two-parameter foundation models like Pasternak or Vlazow foundations (Selvadurai, 1979), which consider additional stiffnesses such as a transverse foundation shear stiffness, do not have a significant effect on the outcome of the analysis.
The ODE governing horizontal displacements, i.e., the deformations of an elastically bedded rod is obtained by inserting the 20 derivative of (A1) into (A4). The following rearrangements are necessary to obtain differential equations of the deflection w and the curvature ψ of a Timoshenko beam on an elastic foundation: Plugging the derivative of (A2) into (A5) yields κGA (w + ψ ) = k n w − q n .
Adding 0 = EIw −EIw to the balance of moments (A6), differentiating twice and plugging the result into the first derivative of (A6) yields Now we use the third derivative of the elasticity law of shear deformations (A2) to substitute the first term of the right hand side and obtain Q = EI κGA Q + EIw .
Substituting Q using (A5) and Q using the second derivative of (A5) yields the ODE of the beam deflection: q n = EIw − EIk n κGA w + k n w.
In order to obtain the ODE of the rotation of normals to the beam mid-surfce we insert (A2) and the derivative of (A3) into 5 (A6) which yields EIψ = κGA(w + ψ).
Code availability. The analysis code of both the modeling framework in part I and the mixed-mode failure criterion based on this framework are available under https://github.com/2phi/weac.
Author contributions. PLR and PW designed, conceived and developed the model. PW drafted the paper and PLR carried out the numerical reference analysis. Both authors contributed equally to the the presented model and the final version of the paper.