Tijan Mede

Guillaume Chambon
email: guillaume.chambon@irstea.fr

François Nicot

Pascal Hagenmuller

Micromechanical investigation of snow failure under mixed-mode loading

Keywords: snow failure regimes, avalanche, second-order work, damage localization, force chain, collapse

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The concept of failure is somewhat elusive when it comes to geomaterials, but globally speaking, it could be defined as the state corresponding to an ultimate load, beyond which the system cannot sustain any larger loads [START_REF] Wan | Failure in Geomaterials: A Contemporary Treatise[END_REF].

Conventionally, the analysis of failure in geomaterials is based on the Mohr-Coulomb criterion which restricts the possible stress states to lie within a plastic limit surface, where unlimited strains can occur under constant stress through a plasticity flow rule [START_REF] Terzaghi | Theoretical soil mechanics[END_REF]. However, geomaterials that exhibit nonassociated plastic flow were shown to fail well inside the plastic flow limit under certain conditions [START_REF] Lade | Static instability and liquefaction of loose fine sandy slopes[END_REF].

A more general approach to material failure can be found in Hill's secondorder work criterion [START_REF] Hill | A general theory of uniqueness and stability in elastic-plastic solids[END_REF], which indicates the existence of an unstable mechanical state that, under suitable loading conditions, can lead to an effective failure. Recent studies have shown that the second-order work can succesfully detect developing material instabilities characterized by an outburst of kinetic energy [START_REF] Nicot | Bifurcation in granular materials: An attempt for a unified framework[END_REF][START_REF] Daouadji | Diffuse failure in geomaterials: Experiments, theory and modelling[END_REF][START_REF] Wautier | A micromechanical approach of suffusion based on a length scale analysis of the grain detachment and grain transport processes[END_REF], and is nowadays widely considered a fundamental building block of material failure theory [START_REF] Wan | Failure in Geomaterials: A Contemporary Treatise[END_REF].

While the second-order work criterion relates the onset of material instability to the evolution of macroscopic stress-strain state, a deeper understanding of failure can be sought by observing the deformational mechanisms at the microscopic scale. The overall behaviour of granular materials has been shown to rely on a limited number of grains, organized in linear columns referred to as force chains [START_REF] Drescher | Photoelastic verification of a mechanical model for the flow of a granular material[END_REF][START_REF] Radjai | Bimodal character of stress transmission in granular packings[END_REF][START_REF] Cambou | Micromechanics of granular materials[END_REF]. These force chains are transferring forces through the granular assembly [START_REF] Nicot | Force chain collapse as grain column buckling in granular materials[END_REF] and the macroscopic material stability is closely related to the stability of force chains. Tordesillas et al. [START_REF] Tordesillas | Buckling force chains in dense granular assemblies: physical and numerical experiments[END_REF] showed that macroscopic strain softening is associated to force chain buckling. Recently, it was also shown that an incremental stress loading leading to the vanishing of the macroscopic second-order work provokes the unjamming of force chains, which results in their bending accompanied by an increase in the kinetic energy [START_REF] Wautier | Micro-inertia origin of instabilities in granular materials[END_REF].

Cohesive granular materials are specific in their capacity to carry tensile forces, and to form looser assemblies of grains [START_REF] Valverde | Random loose packing of cohesive granular materials[END_REF][START_REF] Jiang | An investigation on loose cemented granular materials via DEM analyses[END_REF]. While loose cohesive assemblies exhibit, under growing confining pressure, important structural changes and an irreversible decrease of porosity, their mechanical behaviour has been shown to be predominantly influenced by the initial microstructure [START_REF] Gilabert | Computer simulation of model cohesive powders: influence of assembling procedure and contact laws on low consolidation states[END_REF][START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic loads[END_REF]. More specifically, contact density has been shown to have a decisive influence on the Young's modulus and compressive strength during the initial quasi-elastic phase, while the subsequent plastic consolidation phase was shown to be independent of the initial conditions [START_REF] Gaume | Scaling laws for the mechanics of loose and cohesive granular materials based on Baxter's sticky hard spheres[END_REF].

Snow is a somewhat specific geomaterial due to its extremely high porosity and a variety of existing microstructural arrangements [START_REF] Fierz | The International Classification for Seasonal Snow on the Ground, IHP-VI Technical Documents in hydrology[END_REF]. Despite the fact that snow failure is a key ingredient for numerical modelling of snow avalanche release, our understanding of related processes remains relatively limited [START_REF] Schweizer | Snow avalanche formation[END_REF]. Among the different types of snow avalanches, dry slab avalanches are known for their destructive force and complex release mechanism. Generally speaking, they are triggered by a failure in a weak layer of snow, underlying a cohesive slab on a slope [START_REF] Jamieson | Texture and strength changes of buried surface-hoar layers with implications for dry snow-slab avalanche release[END_REF]. This failure diminishes the support to the overlying slab, increasing the slab tensile forces, which can lead to a tensile crack in the slab and the release of a snow avalanche [START_REF] Mcclung | Mechanics of snow slab failure from a geotechnical perspective[END_REF]. Weak layer failure involved in avalanche release typically occurs under mixed-mode shear and normal loading, provided by the weight of the overlying slab (as well as potential additional perturbations such as skiers, animals etc.) [START_REF] Mcclung | Shear fracture precipitated by strain softening as a mechanism of dry slab avalanche release[END_REF], low temperatures (T < 0 • C) and relatively high strain rates ( ε > 10 -4 s -1 ) [START_REF] Schweizer | Snow avalanche formation[END_REF]. Under these conditions snow can be considered a highly porous brittle [START_REF] Narita | An experimental study on tensile fracture of snow, Contributions from the institute of Low Temperature[END_REF][START_REF] Schweizer | Laboratory experiments on shear failure of snow[END_REF] cohesive granular material [START_REF] Hagenmuller | Numerical simulation of microstructural damage and tensile strength of snow[END_REF][START_REF] Hagenmuller | Microstructure-based modeling of snow mechanics: a discrete element approach[END_REF]. Snow failure in the brittle regime has been investigated by several studies, and different models have been proposed to describe the plastic limit surface. Initially a Mohr-Coulomb criterion [START_REF] Mcclung | Direct simple shear tests on snow and their relation to slab avalanche formation[END_REF] has been employed and although this simple criterion can accurately model snow failure on steep slopes [START_REF] Gaume | Snow fracture in relation to slab avalanche release: critical state for the onset of crack propagation[END_REF], more recent studies indicated that snow features a closed failure envelope. Both Mohr-Coulomb criterion with a cap [START_REF] Reiweger | A new mixed-mode failure criterion for weak snowpack layers[END_REF] as well as elliptical failure envelopes [START_REF] Gaume | Dynamic anticrack propagation in snow[END_REF] have been proposed. In certain cases, snow failure has been observed to be coupled with a volumetric collapse [START_REF] Jamieson | Texture and strength changes of buried surface-hoar layers with implications for dry snow-slab avalanche release[END_REF][START_REF] Johnson | Seismic measurement of fracture speed in a weak snowpack layer[END_REF], somewhat resembling com-paction bands in porous rocks [START_REF] Mollema | Compaction bands: a structural analog for anti-mode I cracks in aeolian sandstone[END_REF]. Recently Gaume et al. [START_REF] Gaume | Dynamic anticrack propagation in snow[END_REF] have shown that the anticrack propagation is a key feature of dry slab avalanche release and stressed the importance of accurate description of snow mechanics at the material scale. Despite all the recent advances in snow mechanics, not much is known about the micromechanical aspects of snow failure [START_REF] Schweizer | Avalanche release[END_REF][START_REF] Schweizer | On recent advances in avalanche research[END_REF]. The main reason for this lack of knowledge is the low strength of snow and its sensitivity to environmental conditions, which render systematic experimental exploration difficult and complicate observation at the microscopic level. As an alternative to experiments, numerical approaches have been developed with the aim to reproduce the mechanical behaviour of snow, based on microstructural input data [START_REF] Hagenmuller | Numerical simulation of microstructural damage and tensile strength of snow[END_REF][START_REF] Hagenmuller | Microstructure-based modeling of snow mechanics: a discrete element approach[END_REF][START_REF] Gaume | Scaling laws for the mechanics of loose and cohesive granular materials based on Baxter's sticky hard spheres[END_REF]. These approaches offer the capacity to perform multiple loading simulations on the same sample, along with a perfect control over boundary conditions and access to microscopic descriptors. A recently developed discrete element (DEM) model of snow [START_REF] Mede | A medial axis based method for irregular grain shape representation in DEM simulations[END_REF] is thus used to obtain a deeper understanding of snow failure by studying the microscopic processes that lead to the volumetric collapse at the macroscopic scale. This paper can be viewed as a continuation of a recently published paper on the macroscopic response of snow to mixed-mode loading [START_REF] Mede | Snow failure modes under mixed loading[END_REF]. The aim of the present paper is to discuss these macroscopic responses from a microstructural perspective. Consequently, some aspects of previous work will be recalled and expanded. Section 2 presents the methodology. Section 3 recalls the main aspects of macroscopic failure, and introduces new results concerning the evolution of damage in the samples. The bulk of the original contribution of this paper is contained in section 4, where the micromechan-ical origins of snow failure are investigated. The paper is concluded by a discussion of results in section 5 and conclusions and perspectives in section 6.

Methodology

Discrete element simulations

Since the scope of this study is restricted to failure and immediate postfailure response of snow in the brittle regime, where deformation is dominated by inter-granular damage and grain rearrangements, individual snow grains can be modelled as rigid bodies [START_REF] Johnson | Identifying microstructural deformation mechanisms in snow using discrete-element modeling[END_REF][START_REF] Hagenmuller | Microstructure-based modeling of snow mechanics: a discrete element approach[END_REF]. This simplification sufficiently diminishes the number of degrees of freedom to render computational times manageable, and warrants the use of the discrete element method (DEM) [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. A specific DEM based approach has been developed by the authors of this paper in order to simulate the mechanical response of snow under mixed-mode loading [START_REF] Mede | A medial axis based method for irregular grain shape representation in DEM simulations[END_REF]. The method utilizes images of snow derived from X-ray microtomography as input information on snow microstructure. The grayscale attenuation images with a typical resolution of 10µm are binary segmented into pore space and a continuous ice matrix [START_REF] Hagenmuller | Energybased binary segmentation of snow microtomographic images[END_REF]. The binary image of the ice structure is then segmented into individual grains by detecting potential weak points in the ice matrix based on local geometrical criteria [START_REF] Hagenmuller | Snow as a granular material: Assessment of a new grain segmentation algorithm[END_REF]. In order to fully account for the snow microstructure in DEM, grain shapes are modelled by an original method, where the shape of a grain is captured by packing its volume with a set of overlapping spheres, distributed along the medial axis of the grain [START_REF] Mede | A medial axis based method for irregular grain shape representation in DEM simulations[END_REF] and clumped into an unbreakable entity. The method allows for grain shapes to be directly accounted for, while still harvesting the inherent efficiency of DEM in handling spherical contacts.

Even though the method enables a perfect reconstruction of the grain shape, an approximation rather than reconstruction of grain shapes was found to produce a good balance between the numerical cost and accuracy of the simulations [START_REF] Mede | A medial axis based method for irregular grain shape representation in DEM simulations[END_REF]. The approximation is produced by selecting a subset of all spheres needed for a complete grain reconstruction according to two parameters (D and R min ).

These approximated grains are spatially arranged according to the initial microstructure of the sample, as shown in Fig. 1. Grains themselves are assumed to be rigid and bonded to adjacent grains by elastic brittle cohesive bonds at the locations where grain contacts were identified by the segmentation algorithm. Since each grain consists of a multitude of spheres, the elastic brittle cohesive bonds have to be modelled between the constituting spheres of the two contacting grains (Fig. 1). Normal contact stiffness K N and adhesion A between each pair of spheres are set proportionally to the area of the grain contacting surface that they are occupying (S S ):

K N = S S E (1)
A = S S C (2)
and shear contact stiffness K S as a fraction of the normal contact stiffness:

K S = νK N , (3)
where E is the Young's modulus, C the tensile strength and ν the Poisson's ratio of the grain constituting material. Each initial cohesive sphere contact exhibits a linear relation between the normal displacement of spheres x N and the tensile force F N , where the maximal admissible tensile force is limited by the value of adhesion (Fig. 1):

F N = K N x N ≤ A. (4)
The shear force F S depends linearly on the shear displacement of spheres x S , where the maximal admissible tensile force is limited by the sum of adhesion and Mohr-Coulomb friction (Fig. 1):

F S = K S x S ≤ A + F N µ, ( 5
)
where µ is the friction coefficient. Once cohesion is broken in tension, it is no longer active in shear and vice-versa. As long as the spheres remain in contact after the cohesion has been broken, Mohr-Coulomb friction remains active in shear. Since the contact force between two grains is a sum of all sphere contact forces between the pair of grains, grain contacts exhibit qualitatively similar rheology as the sphere contacts (Fig. 1). Grain contact stiffness (in normal and shear respectively) is the sum of all sphere contact stiffnesses (in normal and shear respectively) and grain adhesion is the sum of all sphere contact adhesions. Thus the desired elastic brittle cohesive initial intergranular contact law is insured, where contact adhesion and stiffness depend linearly on the grain contact surface, independently of the choice of grain approximation parameters. It should also be added that new sphere contacts that arise due to grain rearrangement are modelled as non-cohesive elastic frictional. The contact stiffness in this case is calculated by applying:

K N = (l 1 l 2 /l 1 + l 2 )E (6)
K S = νK N , (7)
where l 1 and l 2 are the distances from the centres of the two contacting spheres to the contact point [START_REF] Šmilauer | Yade Reference Documentation[END_REF]. In this study, we considered a cubical sample of rounded grain (RG) snow with a side length of 5 mm and initial density of 250 kg m -3 , which was exposed to direct simple shear in order to imitate the naturally occurring mixed-mode loading of snow. A rigid boundary condition is applied to the top and bottom faces, while a periodic boundary condition is applied to all four side-faces (Fig. 1). The sample is loaded by applying a constant shearing velocity v = 1 cm s -1 and constant normal stress p to the top surface, while keeping the bottom surface fixed. Application of a constant shearing velocity, instead of a constant shear stress, enabled us to follow the postfailure response of the samples. The normal stress is applied first with a ramp function of appropriate duration to avoid inertial effects. Shearing is applied once kinetic energy due to the initial compression has dissipated. It has been verified that the applied shear rate is low enough to ensure a quasistatic shearing regime, i.e. that the stress-strain response is independent of v.

Series of simulations were performed under different values of normal stress

p, ranging from -5 kPa in tension to 10 kPa in compression. The simulations were stopped once the macroscopic shear strain reached the value γ = 0.05.

The simulations are performed with DEM solver YADE [START_REF] Šmilauer | Yade Reference Documentation[END_REF]. The initial cohesive sphere contacts are modelled with brittle elastic law, using a contact strength C = 10 6 Pa [START_REF] Haynes | Tensile strength of ice under triaxial stresses[END_REF] and Young's modulus E = 10 8 Pa. Intergranular sintering [START_REF] Colbeck | Sintering in a dry snow cover[END_REF][START_REF] Szabo | Subsecond sintering of ice[END_REF] is not accounted for in the scope of this paper. It has been verified that the utilized contact stiffness ensures the simulations are performed in the rigid grain limit [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF]. Contact friction coefficient µ = 0.2 is fixed according to typical ice values [START_REF] Hagenmuller | Microstructure-based modeling of snow mechanics: a discrete element approach[END_REF]. A Cundall's non viscous damping coefficient 0.02 is applied. Gravity is not taken into account as it has a negligible effect on the response of the samples compared to external loading. Lastly, the mass of each grain is derived from the original binary image. Due to volumetric errors induced by grain approximation [START_REF] Mede | A medial axis based method for irregular grain shape representation in DEM simulations[END_REF], this results in an effective density typically increased by 15% compared to that of ice.

Sensitivity analysis

Analogously to [START_REF] Mede | A medial axis based method for irregular grain shape representation in DEM simulations[END_REF], a sensitivity analysis of the simulated stress-strain response to the grain approximating parameters D and R min has been performed. Mechanical simulations were performed for nine combinations of grain approximating parameters (Fig. 2a). The combination A (D = 0.2, R min = 2) represents the control combination, offering a near-exact reconstruction of grains. This approximation leads to an error of the total ice volume of the reconstructed grains of the same order (≈ 2%) as the one due to X-ray image processing [START_REF] Hagenmuller | Sensitivity of snow density and specific surface area measured by microtomography to different image processing algorithms[END_REF]. The results (Fig. 2b) produced by parameter combinations B-I are compared against the reference response produced by parameter combination A. Mechanical error E M [START_REF] Mede | A medial axis based method for irregular grain shape representation in DEM simulations[END_REF] is introduced and defined as the normalized root mean square difference between the strain-stress response obtained in any given simulation and that obtained in the reference simulation. All simulations of this sensitivity analysis were performed for the value of normal stress p = 4 kPa and the shear stress responses were compared in the range of shear strains 0 < γ < 0.035. As expected, it is observed that while a coarse approximation results in low number of spheres and relatively high E M , vice-versa is true for a fine approximation (Fig. 2a). Furthermore, although E M is closely related to the number of approximating spheres (and consequently the numerical cost), a fast decay of E M with the number of spheres allows accurate mechanical results to be obtained at a relatively low computational cost. Parameter combination E (D=0.5, R min =5) was identified as the optimal one for this study.

Definitions and notations

A relative time scale is used throughout this paper where the time t is normalized by the time of failure t f ail . Such a time scale is employed as it enables comparisons between the different failure regimes (see section 3).

The definition of failure is based on the second-order work criterion. It has been shown that the vanishing second-order work is a necessary condition for the existence of material instabilities [START_REF] Nicot | Bifurcation in granular materials: An attempt for a unified framework[END_REF][START_REF] Wan | Diffuse instabilities with transition to localization in loose granular materials[END_REF][START_REF] Daouadji | Diffuse failure in geomaterials: Experiments, theory and modelling[END_REF]. In the present case, a normalized definition of the second-order work is utilized:

W N 2 (t) = ∂Ω δ f (t) • δ u(t) dS max t ∂Ω δ f (t) • δ u(t) dS , (8)
where δ u and δ f are incremental displacement and incremental force per surface unit along the sample boundary ∂Ω. Failure time t f ail can thus be defined by:

W N 2 (t f ail ) = 0 (9)
Force chains are identified by singling out the grains that comply to the following characteristics [START_REF] Peters | Characterization of force chains in granular material[END_REF]:

• A compressive principal stress higher than the mean principal stress in the observed granular medium. Principal stresses in a grain are calculated by accounting for the contact forces of all the grain constituting spheres according to Love-Weber formula [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF].

• The principal stress direction of force chain particles is aligned with the vectors connecting the mass centre of the observed grain with the mass centres of its neighbours in the force chain. An angle of no more that 45 • is allowed;

• Force chains are composed of at least three contacting grains.

In accordance with the utilized convention the principal stress angle of 0 • marks a pure compression, while 90 • refers to pure shear loading. Analogously a force chain angle of 0 • will refer to a vertical force chain, while 90 • will denote a horizontal force chain.

Macroscopic response

The response of the considered snow sample to shearing under varying values of normal stress p is presented by a series of snapshots in Fig. 3, while the evolution of vertical profiles of local damage within the samples for the same values of normal stress is shown in Fig. 4. Local damage (Fig. 3) is defined for each grain as the ratio of broken cohesive bonds with the neighbouring grains, whereas global damage is defined as the total ratio of broken cohesive bonds. The macroscopic loading and response of the sample are displayed in Fig. 5. Figure 5a summarizes the loading conditions: normal stress is first applied to the specimen, resulting in a principal stress angle close to 0 • . The fact that the principal stress angle is slightly negative during the normal loading phase presumably results from the fact that periodic boundary conditions in DEM, particularly for a material with such a low volume fraction, act like unconfined boundary conditions at low strains, as a certain deformation is needed to establish the mirror contacts at the periodic boundary. The macroscopic response curves in terms of shear stress, normal strain and global damage (Figs. 5b-d) during this initial phase are marked by dotted lines. Once the normal loading is applied and kinetic energy of the sample sufficiently low, a constant shear rate is applied, resulting in an increase of the principal stress angle (Fig. 5a). The macroscopic responses of the sample upon application of shearing are marked by solid curves in Figs.

5b-d.

Failure regimes

In general, three qualitatively different regimes of failure can be observed, depending on the value of applied normal stress:

Failure regime A: at low normal stresses (for the tested sample, p < 2 kPa), the normal loading phase does not produce any significant damage (Fig. 5d) or normal strain (Fig. 5c) of the sample. Subsequent application of shearing results in a quasi-elastic shear stress build-up (Fig. 5b) until a peak shear stress is reached, marking the sample failure. Normal strain and global damage at failure onset are low and local damage appears dispersed (Fig. 3). Failure is followed by strain softening, which gradually diminishes shear stress. The sample fails in a narrow horizontal region (Figs. 3 and4),

where the local damage damage is highly concentrated, while the rest of the sample remains largely intact, resulting in a low macroscopic normal strain and global damage throughout the simulation. It should be noted that a slight dilatancy is observed immediately after failure. This dilatancy might however be caused by numerical artifacts, as increase in sample height remains on the order of the resolution of our numerical samples. Failure regime B: at moderate normal stresses (for the tested sample, 2 ≤ p ≤ 8.5 kPa), the pre-failure response is qualitatively similar to that observed in failure regime A. However, as opposed to regime A, the peak shear force is followed by a volumetric collapse of the specimen (Fig. 5c) and rapidly propagating damage that eventually decomposes the specimen into a cohesionless assembly of grains (Fig. 3). The sample collapse appears to propagate from the bottom of the sample upwards (Fig. 4) and has been

shown to be a dynamical process, which initiates at the point of peak shear force and advances spontaneously once initiated [START_REF] Mede | Snow failure modes under mixed loading[END_REF]. In contrast to gradual strain softening and shear controlled post-failure response in regime A, regime B is marked by rapidly vanishing stresses and a time-driven response during the collapse phase. Once the sample porosity is sufficiently diminished, normal strain and global damage stabilize.

Failure regime C: at sufficiently high normal stresses (for the tested sample, p ≥ 9 kPa), the sample is observed to fail and collapse (Fig. 5c) during the application of the normal force, resulting in a rapid increase of damage (Fig. 5d). Similarly as in regime B, the collapse appears to propagate from the bottom of the sample upwards (Fig. 4) and decomposes the sample into a cohesionless assembly of grains (Fig. 3). The normalized second-order work in failure regimes A and B (Fig. 6) exhibits a roughly constant value during the normal loading phase, followed by a steep increase at the onset of shearing (t/t f ail ≈ 0.5). Peak value is reached at t/t f ail ≈ 0.6, after which the second-order work slowly diminishes until a fast decrease takes place just before failure (t/t f ail = 1). Normalized second-order work in failure regime C also exhibits a roughly constant value during the normal loading phase (t/t f ail < 0.99), although the plot becomes somewhat more noisy at t/t f ail > 0.5). In accordance with the chosen loading protocol, shearing is not started immediately after the concluded normal loading phase, due to increasing kinetic energy. The instability thus develops in the absence of shearing. Although the prescribed external force is kept constant, oscillations of the normal force are measured on the sample boundary and a negative second-order work is recorded as a result at t/t f ail = 1. It should also be noted that the different values of the normalized second-order work in the three failure regimes during the normal loading phase are a result of the normalization process, where the calculated value of the second-order work is divided by the highest recorded value before failure (Eq. 8).

Second-order work

Global damage and failure envelope

The level of global damage, necessary to trigger failure is extremely low through the entire range of simulated normal stresses (Fig. 7). Moreover, the global damage at failure appears almost completely independent of the applied normal stress. Evidently, failure within a small fraction of the load bearing contacts of the snow sample suffices for the entire sample to fail. Consistent failure regimes and qualitatively similar failure envelopes were observed for snow samples of varying snow types and densities, exposed to direct simple shear simulations [START_REF] Mede | Snow failure modes under mixed loading[END_REF]. The observed response can hence be regarded as a general snow characteristic, rather than a particularity of the studied snow sample. The reasons for the fundamentally different response of snow in the three identified failure regimes, particularly the occurrence of volumetric collapse in regimes B and C, while snow samples appear to main-tain a high normal load bearing capacity and structural integrity in regime A, are sought by observations at the microscopic level.

Microscopic investigation

DEM simulations provide access to all relevant microscopic descriptors, enabling the analysis of microstructural evolution and the onset of instability that leads to failure. The micromechanical analysis of instability is presented in two steps. In the first step, the distribution of all the contact forces within the snow sample is observed. In the second step, the analysis is focused on the distribution of contact forces within and around the load carrying fraction of the snow sample.

Contact force distribution

The distribution of particle contact forces during the course of the simulations is analysed for all three failure regimes. Let us recall that each grain contact consists of multiple approximating sphere contacts. To evaluate how close the different cohesive sphere contacts within the sample are to failing, we define contact charge C C as as the ratio between the sphere shear contact force F S and the maximal admissible contact force

F crit : C C = sgn(F N ) F S F crit = sgn(F N ) F S A + F N µ , ( 10
)
where the maximal admissible contact force is defined by the Mohr-Coulomb criterion (Eq. 5). Value C C = 0 refers to a non-loaded cohesive contact and is observed in regimes B and C, compared to regime A, the differences be-tween distributions of contact charges at the point of failure for all three failure regimes remain small, particularly in tension. The analysis of contact charges in the sample thus appears insufficient to explain the marked difference between different failure regimes of snow.

C C = ±1,

Force chain analysis

The ultimate strength of a granular material is controlled by the stability of the force chains [START_REF] Nicot | Force chain collapse as grain column buckling in granular materials[END_REF], and a developing material instability is generally reflected in changes of the force chain network. Force chain buckling has been shown to anticipate material failure [START_REF] Tordesillas | Buckling force chains in dense granular assemblies: physical and numerical experiments[END_REF], which is marked by rapid decrease in the force chain lifespans and an increase in the emerging/dying rate of force chains [START_REF] Wautier | Micro-inertia origin of instabilities in granular materials[END_REF]. Force chain network was thus identified in the loaded snow sample and its evolution was monitored through the course of the direct simple shear simulations in order to characterize the micromechanical origins of different failure regimes. An example of the spatial arrangement of force chains for the shearing simulation under normal stress p = 4 kPa is shown in Fig. 11. The force chains within the sample can be observed to align with the principal stress angle: they assume vertical orientations during the normal loading phase (Fig. 11a) and increasingly tilted orientations as the shear load is added (Figs. 11b and11c). Additionally, the fact that several identical force chains can be observed in Fig. 11b as well as Fig. 11c showcases the long force chain lifespans during the pre-failure phase. This is consistent with the fact that no grain rearrangement takes place in the sample, enabling force chains to persist in an unchanged configuration.

The average force chain lifespans during the course of the simulations for three different values of normal stress, corresponding to the three respective failure regimes, are analysed in Fig. 12. Figure 13 displays the average force chain angles for the same three simulations. While an alignment of the force chains with the principal stress angle was already noted in Fig. 11, this observation is underlined by the striking similarity between the plot of average force chain angles (Fig. 13) and the plot of principal stresses within the snow sample (Fig. 5a). Evidently the force chains align in the way that allows them to transmit the maximal load.

In failure regime A the normal loading phase (t/t f ail < 0.5) results in a monotonically increasing average force chain lifespan and a constant average force chain angle of -20 • . Although 0 • angle would be expected under normal loading, it should be recalled that the principal stress angle (Fig. 5a) during the normal loading phase is actually slightly negative (see section 3). Additionally, the highly porous rigid structure of snow limits the ability of force chains to adapt in terms of direction. Application of shearing (t/t f ail ≈ 0.5)

immediately results in a slight drop of the average force chain lifespan (Fig. 12), and a steep increase in force chain angle (Fig. 13). The average force chain angle showcases a steep increase at the onset of shearing, which slows down considerably at t/t f ail ≈ 0.6, consistent with the loading angle and the post-peak decrease of the second-order work (Fig. 6). Force chain lifespans then oscillate around a constant value and their average angle of inclination continues to slowly increase until failure (t/t f ail = 1). No sudden change in the force chain network is observed at failure -the average angle of inclination starts gradually diminishing while the average lifespan continues to oscillate around a roughly constant value.

In failure regime B, force chain development for t/t f ail < 1 appears very similar to that observed in regime A. The post-failure evolution of the force chain network is however significantly different. Consistent with the onset of the volumetric collapse, a major force chain breaking event is observed, resulting in the average force chain lifespan dropping to a value close to zero.

During the collapse, force chains align in the vertical direction due to inertial effects linked to the collapse.

In failure regime C average force chain lifespan exhibits a monotonic increase until a minor decrease is observed at t/t f ail ≈ 0.5. This point roughly coincides with the onset of oscillations of the second-order work (Fig. 6). The average force chain lifespan reaches a second peak at t/t f ail ≈ 0.95, when a local minimum of the normalized second-order work is recorded. Secondorder work at this point is still positive, but starts oscillating severely, along with a rapid decrease in of the average force chain lifespan, leading to the onset of collapse at t/t f ail = 1 and consequent grain rearrangement. The average force chain angle starts at -20 • and then gradually increases to 0 • during the normal loading phase. No abrupt change of the average force chain angle is observed at failure, consistently with a roughly constant principal stress angle in this phase of the simulation (Fig. 5a). This section focuses on contact forces within the force chains. Due to defining characteristics of force chains (see section 2.3), contact forces between individual grains of force chains are strictly compressive on the grain 29 scale. This is however not necessarily true for the sphere contacts, as can be observed in Fig. 14. While the majority of sphere contact forces in all three regimes are indeed compressive, a smaller fraction of contact forces appears to be tensile. Similarly as contact forces in the whole specimen (Fig. 9), the magnitudes of tensile as well as compressive sphere contact forces gradually increase in all three failure regimes as failure draws closer. However the number of high tensile sphere contacts within the force chains is distinctly higher in regimes B and C compared to regime A, especially in the short interval before the onset of failure (t/t f ail > 0.7). Since the sum of all sphere contact forces for each pair of contacting grains within force chains is compressive, the existence of tensile sphere contacts can be interpreted as force chain bending. More specifically, since the force chains were shown to align with the principal stress angle, this bending is likely related to buckling.

Force distribution within force chains

Volumetric collapse in regimes B and C thus appears to be anticipated by pronounced force chain buckling and failure can be expected when the force chains become unstable.

Force distribution around force chains

It is generally accepted that compressive confining forces insure the stability of force chains [START_REF] Nicot | Force chain collapse as grain column buckling in granular materials[END_REF] in dense cohesionless granular assemblies. The degree of lateral support has been shown to determine the extent to which force chains can resist buckling [START_REF] Tordesillas | On the modeling of confined buckling of force chains[END_REF]. Tensile lateral forces on the force chains, provided by the so-called weak network of grains, have been observed to play a similar stabilizing role in cohesive granular materials [START_REF] Radjai | Cohesive granular texture[END_REF][START_REF] Kadau | Contact dynamics simulations of compacting cohesive granular systems[END_REF].

In the case of snow, due to its highly porous nature, the lateral support on the force chains is exerted by a relatively low number of grain contacts.

Consequently, failure of any of these contacts, whether it be compressive or tensile, is expected to introduce a significant unbalanced force on the force chain, which might allow force chain buckling to develop. In order to analyse the stabilizing effect of the lateral forces on the force chains, contact charges

C F CS C
(Eq. 10) between the force chain grains and the surrounding grains were calculated. Analogue results were also recorded on other snow types, although normal stress values at the transitions between regimes A, B and C were observed to vary substantially depending on snow density and microstructure [START_REF] Mede | Snow failure modes under mixed loading[END_REF].

The range of normal stresses, where these transitions are observed, coincides with the relevant range of normal stresses typically involved in slab avalanche release (0.6 < p < 4 kPa).

The volumetric collapse of weak snow layers has been shown to have a key influence on the release of dry snow slab avalanches [START_REF] Gaume | Modeling of crack propagation in weak snowpack layers using the discrete element method[END_REF][START_REF] Gaume | Dynamic anticrack propagation in snow[END_REF]. Although this phenomena is well recorded [START_REF] Van Herwijnen | Measurement of the deformation field associated with fracture propagation in weak snowpack layers[END_REF][START_REF] Johnson | Seismic measurement of fracture speed in a weak snowpack layer[END_REF], the conditions under which it appears re-main unclear. The DEM model developed in the scope of this project, which successfully reproduces the collapse of snow samples, was already utilized to show how volumetric collapse at failure depends on the level of applied normal stress [START_REF] Mede | Snow failure modes under mixed loading[END_REF]. In this paper we go further and explore the microscopic mechanisms that lead to the volumetric collapse. Analysis of particle contact forces showed that compressive as well as tensile forces are present in the sample during the entire loading phase before the onset of failure. This is consistent with the work done by Hagenmuller et al. [START_REF] Hagenmuller | Snow as a granular material: Assessment of a new grain segmentation algorithm[END_REF][START_REF] Hagenmuller | Numerical simulation of microstructural damage and tensile strength of snow[END_REF], who recorded a large fraction of tensile inter-granular forces in FEM simulations on real snow microstructure. Gaume et al. [START_REF] Gaume | Scaling laws for the mechanics of loose and cohesive granular materials based on Baxter's sticky hard spheres[END_REF] on the other hand recorded a much lower fraction of tensile particle forces in their DEM simulations of snow compression. This discrepancy could be attributed to the fact that in the latter simulations, snow was modelled an assembly of spherical particles, whereas in the present work each snow grain is represented by a collection of spheres. Furthermore, while the distribution of tensile and compressive forces in regime A is rather symmetric, the ratio of compressive contact forces forces increases in regimes B and C, consistently with the higher value of applied normal stress. At the point of failure, a larger proportion of high contact charges in regimes B and C suggests that the number of contacts close to failing is increased compared to regime A, which can be considered a partial explanation of why a sudden collapse is observed at failure in these two regimes.

A more complete micromechanical explanation was sought by analysing the load carrying fraction of grains within the snow sample. The analysis of sphere contact forces within force chains showed that while the majority of sphere contact forces are compressive, a non-negligible amount of contact forces are tensile. Moreover, the ratio of high tensile forces in regimes B and C is decidedly higher than in regime A, particularly when the failure is near. Since the contact forces between the force chain members on the scale of the grain are strictly compressive, consistently with the utilized definition of a force chain, tensile sphere contact forces can be interpreted as force chain bending. As the force chains are aligned with the principal stress direction, it can be assumed that this bending is induced by buckling. We thus present evidence that the phenomena of force chain buckling, which has already been identified as the micromechanical process leading to material failure in the case of dense granular media [START_REF] Tordesillas | Buckling force chains in dense granular assemblies: physical and numerical experiments[END_REF] also represents an important mechanism of failure in the case of loose granular materials such as snow.

Furthermore, the analysis of contacts between the force chain members and the surrounding grains showed that at the point of failure in regimes B and C, a relatively high number of contacts are close to failing, which is not the case in regime A. Contacts between the force chain members and the surrounding grains insure the stability of force chains in snow in a similar way that has been observed in dense granular media [START_REF] Tordesillas | On the modeling of confined buckling of force chains[END_REF][START_REF] Nicot | Force chain collapse as grain column buckling in granular materials[END_REF][START_REF] Radjai | Cohesive granular texture[END_REF][START_REF] Kadau | Contact dynamics simulations of compacting cohesive granular systems[END_REF]. However, due to the extremely high porosity of snow, a very small number of grains supply the lateral support to the force chains and failure of a single supporting grain contact will introduce a significant unbalanced force on the force chain and increase the load on the remaining supporting contacts. This effect is even more pronounced in the case of tensile contacts, where the contact force is

Concluding remarks

Thanks to combined use of microstructure-based DEM modelling and state of art micromechanical analysis, macroscopic snow failure under mixedmode loading was related to microscopical physical mechanisms, leading to material instabilities. For sufficiently large values of applied normal stress, high contact charges in the vicinity of force chains were identified to trigger volumetric collapse of snow, by withdrawing support to the force chains and allowing force chain buckling to develop, leading to eventual failure of the load carrying fraction of the snow microstructure. These results offer new insights on the key processes to be accounted for in large-scale simulations of, e.g., snow failure or avalanche release [START_REF] Gaume | Dynamic anticrack propagation in snow[END_REF]. The complexity of the mechanical response of snow (localization, collapse, etc.) does however challenge the possibility of formulating a relevant macroscopic constitutive law for this material. Instead, the development of coupled multi-scale approaches might represent a fruitul avenue to account for this complexity in large-scale simulations (e.g., [START_REF] Desrues | From discrete to continuum modelling of boundary value problems in geomechanics: An integrated FEM-DEM approach[END_REF]).

Although small, the sample size considered in this study has been shown

to result in a representative mechanical response in the case of simulated large strain compression [START_REF] Hagenmuller | Microstructure-based modeling of snow mechanics: a discrete element approach[END_REF] on the same sample. Additionally, the size of the sample was proven sufficient to estimate the minimum cut density [START_REF] Hagenmuller | Characterization of the snow microstructural bonding system through the minimum cut density[END_REF],

which is related to the minimal damage needed to break the sample. We can thus argue that this sample size is also sufficient to monitor representative microstructural evolution under direct simple shear loading. Nevertheless, a proper study of REV size for the particular loading conditions under consideration, should be performed in the future. Additionally, in order to generalize the conclusions of this paper, the present micromechanical study would need to be applied to more snow types. the effects of the competition between damage and fast ice sintering processes [START_REF] Szabo | Subsecond sintering of ice[END_REF][START_REF] Mulak | Numerical investigation of the mixed-mode failure of snow[END_REF]

Figure 1 :

 1 Figure 1: Approximated snow sample, as modelled in DEM, along with boundary conditions and applied loading. Individual grains are marked with varying colours for clarity. Schematic zooms into the rheology of an individual grain contact (F * N and F * S are normal and shear grain contact force respectively, A * is the grain adhesion and µ is the friction coefficient), and into the rheology of an individual sphere contact (F N and F S are normal and shear sphere contact force respectively and A is the sphere contact adhesion) are displayed.

Figure 2 :

 2 Figure 2: a) Dependency of mechanical error E M on the number of approximating spheres for direct simple shear simulations performed with different approximating parameters combinations. b) Simulated direct simple shear responses for normal pressure p = 4 kPa and parameter combinations A, E and I. In both figures, A represents the reference response parameter combination and E the parameter combination that was identified as the optimal one.

Figure 3 :

 3 Figure 3: Snapshots of the distribution of local damage at different times and for three different values of normal stress: p = 1 kPa, p = 4 kPa and p = 9.5 kPa.

Figure 4 :

 4 Figure 4: Evolution of vertical profiles of damage during the course of simulations for three different values of normal stress: p = 1 kPa (failure regime A), p = 4 kPa (failure regime B) and p = 9.5 kPa (failure regime C). Each figure shows histograms of the average local damage computed at different times, and plotted against normalized vertical dimension of the sample (z/z max ).

Figure 5 :

 5 Figure 5: Macroscopic response of the snow sample for three different values of normal stress, corresponding to the three typical failure regimes. The dotted parts of curves denote the normal loading stage, and the grey vertical dotted line marks the onset of failure.

Figure 6 :

 6 Figure 6: Normalized second-order work evolution for three different values of normal stress, corresponding to the three typical failure regimes. For clarity, only the initial parts of second-order work curves are plotted -after the second order work becomes negative and the sample fails, the second-order work plot becomes extremely noisy. The dotted parts of curves denote the normal loading stage and the horizontal dotted grey line marks zero value of second-order work. A zoom is also provided into the vicinity of t/t f ail ≈ 1.

Figure 7 :

 7 Figure 7: Global damage at failure as a function of applied normal stress.

Figure 8 :

 8 Figure 8: Failure envelope of the tested snow sample. The dotted vertical line marks the threshold normal stress, above which a collapse (regimes B and C) of the sample is observed.

  Figure 9 displays histograms of normal contact forces between spheres at different times before the failure. The magnitude of normal contact forces steadily increases in all three regimes as the failure draws closer. Despite the fact that only compressive and shear loading is applied to the snow samples, a large amount of contact forces is observed to be tensile. While the distribution of tensile and compressive forces seems rather symmetric in failure regime A, the ratio of compressive to tensile forces increases with the applied normal stress which 22 produces a somewhat asymmetric distribution in regime B and a clearly asymmetric distribution in regime C.

Figure 9 :

 9 Figure 9: Time evolution of the distribution of normal contact forces between spheres during the simulations for three values of normal stress, corresponding to the three identified failure regimes. In the utilized convention, negative normal contact forces denote tensile forces.

  a contact load sufficient to break the cohesion. The sign of the normal contact force is used to distinguish between tensile (negative sign) and compressive contacts (positive sign).

Figure 10 :

 10 Figure 10: Distribution of contact charges C C at the point of failure (t/t crit = 1) for three values of normal stress corresponding to the three identified failure regimes.

Figure 10

 10 Figure 10 summarizes the distribution of contact charges at the point of failure for all three identified failure regimes. Again, it should be noted that even though the sample is loaded in compression and (in the case of regimes A and B) shear, close to half of the contacts at failure are loaded in tension. While a somewhat higher ratio of high contact charges (|C C | > 0.5)

Figure 11 :

 11 Figure 11: Snapshots of the simulations under normal stress p = 4 kPa with highlighted force chains.

Figure 12 :

 12 Figure 12: Average force chain life span evolution during the simulations for all three identified failure regimes. The dotted part of each curve denotes the normal loading stage and the vertical dashed line marks the onset of failure.

Figure 13 :

 13 Figure 13: Average inclination angle of force chains (for definitions see section 2.3) during the simulations for the three identified failure regimes. The dotted part of each curve denotes the normal loading stage and the vertical dashed line marks the onset of failure.

Figure 14 :

 14 Figure 14: Distribution of normal contact forces within the force chains at different times before failure for three values of normal stress, corresponding to the three identified failure regimes. In the utilized convention, negative normal contact forces denote tensile forces.

Figure 15 :

 15 Figure 15: Distribution of contact charges C F CS C

Figure 15

 15 Figure 15 displays the contact charges of sphere contacts between the force chain grains and the surrounding grains at the point of failure. Unlike what has been observed when analysing all the contacts (Fig. 10), the number of contacts around the force chains with a charge |C C | > 0.5 at the point of failure is relatively high in regimes B and C, while such high contact charges are not observed in regime A. In addition, the number of high contact charges in regimes B and C is particularly high in the case of tensile contacts.

  Force chain evolution was shown to reflect the failure regime of the snow sample -the force chains align with the principal stress angle in order to transmit the maximal load through the sample, while their lifespans clearly display the difference between failure in regimes A on one hand, and regimes B and C on the other hand. While a practically unchanged average force chain lifespan at the onset of failure in regime A is a sign of stable force chains, rapidly decreasing force chain lifespans at failure in regimes B and C showcase the force chain collapse. Consequently, the material maintains a relatively high normal load bearing capacity in regime A, associated to localized failure, whereas the failure of the load bearing fraction of the sample in regimes B and C triggers the complete volumetric collapse.

  completely withdrawn after the contact fails. If a large number of these supporting contacts are already close to failing, as observed just before the macroscopic failure in regimes B and C, a chain reaction of failing supporting contacts can be triggered, which withdraws support to the force chain and allows the force chain to buckle and eventually collapse. The average lifespans in regimes B and C are indeed observed to rapidly decrease at the point of failure, while a very slow gradual decrease of lifespans is observed after failure in regime A. The failure of the load carrying fraction of the sample in regimes B and C consequently leads to the volumetric collapse of the sample, while stable force chains result in a localized failure in regime A.

  Due to the relatively high deformation rates involved in the release of snow avalanches, which are well above the rate of transition from ductile to brittle behaviour in snow [ γ > 10 -4 s -1 , 58] [59], viscous effects and grain sintering were not taken into account in the present study. Expanding the model by applying a cohesive-viscous-elastic contact law and accounting for grain sintering could open interesting perspectives in terms of studying, e .g.,

Acknowledgements

We would like to thank Mohamed Naaim for fruitful discussions and support. This project is funded by Labex TEC21 (Investissements d'Avenir, grant agreement ANR-11-LABX-0030). CNRM / CEN and INRAE UR ETNA are part of Labex OSUG@2020 (Investisssments d'Avenir, grant agreement ANR-10-LABX-0056). We thank Jürg Schweizer and an anonymous reviewer for their constructive comments.