Dynamic anticrack propagation in snow

Continuum numerical modeling of dynamic crack propagation has been a great challenge over the past decade. This is particularly the case for anticracks in porous materials, as reported in sedimentary rocks, deep earthquakes, landslides, and snow avalanches, as material inter-penetration further complicates the problem. Here, on the basis of a new elastoplasticity model for porous cohesive materials and a large strain hybrid Eulerian–Lagrangian numerical method, we accurately reproduced the onset and propagation dynamics of anticracks observed in snow fracture experiments. The key ingredient consists of a modified strain-softening plastic flow rule that captures the complexity of porous materials under mixed-mode loading accounting for the interplay between cohesion loss and volumetric collapse. Our unified model represents a significant step forward as it simulates solid-fluid phase transitions in geomaterials which is of paramount importance to mitigate and forecast gravitational hazards.

The return mapping algorithm is the discrete equivalent to solving for a strain that satisfies the plastic flow rule in Equation 5 of the Methods section and that lies in the CamClay yield surface. In this section first we outline the method of Simo and Meschke [1] to derive the discrete equations from their continuous versions. This procedure starts by assuming there is no plastic flow and a return mapping algorithm is derived from the flow equations that shows how to project back to the yield surface if the assumption of no plastic flow is invalid.
Consider the evolution of b E from time t n to time t n+1 = t n + ∆t. We consider this evolution per particle, and thus it is useful to take a Lagrangian view. We outline the notation used in the Lagrangian view in the Methods section of the paper. Specifically useful here is the flow map φ : Ω 0 × [0, T ] → R d , and its relation to the deformation gradient F = ∂φ ∂X . Define the time t n configuration of the material as Ω t n = x = φ(X, t n ) for some X ∈ Ω 0 and defineφ : Intuitively,φ defines the deformation as if the time t n configuration Ω t n of the material is the reference configuration, rather than Ω 0 as in the standard Lagrangian view. This is some times called an updated Lagrangian view. While the deformation gradient F defines the deformation from the initial configuration (Ω 0 ) to the time t configuration (Ω t ), the JacobianF = ∂φ ∂x defines the deformation from the time t n configuration (Ω t n ) to the time t configuration (Ω t ), where t ≥ t n . Also these are related as F =FF n , or more precisely F(X, t) =F(φ(X, t n ), t)F(X, t n ) for all X ∈ Ω 0 .
Define b E * =F −1 b EF−T . Let us consider the difference between the evolution of b E * and b E in absence of plasticity at time t n < t < t n+1 . By the definition of no plastic flow. In this sense, we can see thatb E can be considered as the trial elastic state obtained without any plastic flow. If this does not satisfy the constraint, ∆γ and b E t n+1 must be defined to "project"b E to b E t n+1 . We use this process to define the projection.FF E is considered the trial elastic state, one obtained in the absence of plastic flow. Thus,b E =Fb EFT and we seek the solution of Equation 1 to define the projection to b E t n+1 . This can be done most easily by considering the singular value decomposition ofFF E .
If the singular value decomposition ofFF E is given byFF , then we may write (1) as Multiplying both sides of the previous equation by U T on the left and by U on the right, and taking log results in The model that we choose uses the Hencky-strain as a measure of deformation. By defininĝ tr := logΣ andˆ n+1 := log Σ, we may simplify and rearrange Equation (2) This is our discrete flow rule. In the return mapping algorithm, we want to solve forˆ n+1 satisfies Equation (4) subject to the constraint y(τ (ˆ n+1 )) ≤ 0.
Solving Equation (4) and (5) can be seen as a ray-ellipse intersection problem due to the ellipsoid shape of our CamClay yield surface.

Supplementary Note 2 -Shear behavior of the weak layer
We describe and discuss the mechanical behavior of the weak layer in shear by simulating an unconfined sheartest (simulation setup corresponds to experiment n • 2 with a wall horizontal velocity of 0.01 m.s −1 , Supplementary  Figure 1). After reaching the shear strength of the weak layer, significant softening is observed followed by weak layer collapse (almost 20% of weak layer thickness). Then, the weak layer has a pure frictional behavior with a residual shear strength. This mechanical behavior is very similar to what has been reported in laboratory experiments of shear failure of snow [2] and what has been assumed in interfacial shear models [3,4]. In particular, the so-called critical sliding displacement (displacement δ s associated to the softening zone) has the same order of magnitude as in [3] and [2,4] (δ s = 3.5 mm in [3], δ s = 2 mm in [2,4] and δ s ∼ 2 mm in our simulation).

Supplementary Note 3 -Influence of mesh resolution
We show the effect of mesh refinement on the amount of volumetric collapse of the weak layer. The collapse is controlled by the hardening rule (point (2*) to (3*) in Figure 1 of the paper). It characterizes the thickness of the anticrack (or compacting shear band) and can be evaluated by computing the change in volumetric plastic deformation as 1 − J P with J P = det(F P ) (0 < J P ≤ 1). The collapse height h c can be evaluated by computing h c = (1−J P )D wl . Here, we simulated the shear test presented in Supplementary Note 3 for different mesh resolutions and for different values of the parameters of our hardening rule. Supplementary Figure 2 shows that our model is consistent with mesh refinement as 1 − J P is almost not influenced by the mesh resolution.