Snow Mechanics Near the Ductile‐Brittle Transition: Compressive Stick‐Slip and Snow Microquakes

The rate‐dependent mechanical response of snow is commonly believed to directly inherit the ductile‐to‐brittle transition (DBT) from the viscoplastic ice matrix. Recent work has however stressed the impact of microstructure evolution by rapid sintering during deformation. To understand these phenomena we have conducted deformation‐controlled compression experiments in an X‐ray tomography stage. By varying the strain rate over 3 orders of magnitude, we find an intermediate regime where the stress response changes from smooth to serrated behavior. This regime is accompanied by microstructure quakes associated with strain localization bands in the sample interior. To interpret the results we developed a minimal, scalar model with a rate‐dependent, elastoplastic constitutive law and healing which falls into the general class of rate‐and‐state models. The model correctly predicts the range of the instability as well as amplitude and frequency of the serrations and reveals a formal equivalence of the observed (compressive) stick‐slip with seismic faults.


Introduction
Snow is an archetype of a complex material with peculiar mechanical or rheological properties. The material behavior originates from the viscoplastic properties of ice as constituent material which is aggregated into a highly disordered, cohesive-granular microstructure. This combination causes the rich, rate-dependent behavior of snow (e.g., Narita, 1984;Schweizer, 1998).
The rate-dependent plasticity originates from the dynamics of dislocations in hexagonal ice. A scalar form of Glen's law ∼ . 1∕m (Schulson & Duval, 2009) is often used to describe the increase of viscoplastic flow resistance (and therefore of the average stress level) with increasing strain rate . . A well-known consequence is the ductile-to-brittle transition (DBT) (Kirchner et al., 2001) of snow and ice: Above a critical rate, the plasticity-induced buildup of internal stress concentrations, notably at grain boundaries, leads to local stress levels that may exceed critical stresses for crack nucleation. Thus, the material enters the brittle domain. In the absence of any other time-dependent process, snow at a constant temperature with a given microstructure would, accordingly, only show one qualitative change of its stress-strain characteristics as a function of strain rate, namely, the DBT.
However, recent experiments have provided ample evidence for the relevance of another time-dependent process on the DBT, namely sintering. While the potential impact of sintering on snow mechanics was discussed for a long time (Schweizer et al., 2003), the inclusion of healing mechanisms into models was accomplished only recently. Sintering can be conveniently included into fiber bundle models (Capelli et al., 2018;Reiweger et al., 2009) as a (time-dependent) healing of broken fibers. This extension of standard fiber bundle models explains the measured distribution of acoustic emission events during shear failure (Capelli et al., 2019). As discussed by Barraclough et al. (2017), time-dependent healing also explains the repeated propagation of compaction bands in snow compression. Thus, sintering not only modifies the parameters in It was first stressed by Montmollin (1982) that the DBT in snow is actually a double transition as a consequence of fast metamorphism. In shear experiments, three different regimes were classified as a function of strain rate. Between ductile and brittle behavior an intermediate regime of "cyclic brittle failures" was observed, in accordance with Yosida et al. (1958) who referred to this regime as "brittle of the first kind." The idea behind the intermediate regime is that the time scale of the externally imposed deformation can be either below, in between, or above the two intrinsic time scales of sintering and viscoplasticity. Sintering implies that the bond geometry can no longer be considered time independent and, accordingly, local stress concentrations may decrease over time. This stabilizing mechanism competes with the strain-induced reduction of strength (bond failure) that is associated with deformation in the brittle regime. While this competition is obvious and widely mentioned, its role for the emergence of "the brittle regime of the first kind" has never been detailed either experimentally or theoretically.
The competition between time-dependent stabilizing and driving-dependent destabilizing mechanisms is of broad interest in geosciences, and the resulting stick-slip phenomena have been extensively studied. The most prominent example is the interpretation of earthquakes as signatures of a dynamic instability of steady-state frictional sliding of tectonic plates (Scholz, 1998). The combination of the slip-induced transition to dynamic friction with a healing mechanism for regaining static friction during motion ("aging") causes a bifurcation in the governing equations that occur over a range of imposed slip velocities. The mechanisms underlying stick-slip are not limited to interface friction. In metal plasticity, the competition between external loading which causes dislocations to break free from obstacles, and the time-dependent process of dislocation pinning by diffusing solute atoms (McCormick, 1988) is known as the Portevin-Le Chatelier effect and produces "jerky flow" as an inherent feature of the mechanics.
Accordingly, it can be expected that equivalent dynamic instabilities can be revealed in snow mechanics and that a careful analysis of the rate-dependent deformation response of snow in view of aging phenomena can contribute to a deeper understanding of failure processes in geomaterials. It is the aim of the present letter to experimentally study serration effects in snow deformation and develop a minimal mechanical model for their explanation within the general context of stick-slip phenomena. We show that the "brittle regime of the first kind" (Yosida et al., 1958) near the DBT corresponds to an instability of the stationary solution of the nonlinear model equations and provide experimental evidence that the resulting stress serrations are accompanied by the nucleation of isolated, compressive strain bands in the sample interior. We establish a formal equivalence with rate-and-state friction models used for earthquakes which justifies an interpretation of the observed (compressive) stick-slip as microquakes in the ice matrix of snow.

Snow Sample Preparation
For sample preparation and microcompression experiments we follow the protocol established in . Snow samples were produced from snowmaker snow  where laboratory and water temperatures during production were set to T lab = −25 • C and T water = 20 • C. The new snow was then stored for 1 month under nearly isothermal conditions at T = −25 • C to promote moderate coarsening of the dendritic crystals. Specimens for mechanical testing were prepared at storage temperature by sieving the snow directly into the cylindrical sample holder (radius r 0 = 18 m) of the compression device . The filling level (nominal sample height) was 30 mm. The so obtained samples are characterized by an average density of 202 kg/m 3 and a specific surface area of 31 m 2 ∕kg. After sieving, the samples were transferred to a cold room at a different temperature of T = −15 • C where the samples were mounted in the compression stage for testing.

Microcompression Tests
Compression tests were conducted with the microcompression device (MCD) developed by . The MCD facilitates displacement-controlled, uniaxial compression experiments up to a peak stress of 13.4 kPa for piston velocities v in the range 0.1-200 μm s −1 with an accuracy of 50 Pa. Most experiments were carried out in rack mode, that is, without installing the MCD in the X-ray scanner. Piston velocities were used to compute externally imposed strain rates from nominal sample heights via . ext = v∕30 mm. An 0.2 6.7 · 10 −6 no 1 1 2 6.7 · 10 −5 no 2 1 10 3.3 · 10 −4 no 2 10 20 6.7 · 10 −4 yes 2 1 40 1.3 · 10 −3 yes 2 10 80 2.7 · 10 −3 no 2 10 160 5.3 · 10 −3 no 2 10 200 6.7 · 10 −3 no 2 10 overview of all tests with experimental conditions is given in Table 1. In order to resolve the relevant features of the deformation curves, the sampling rate was increased for higher strain rates.

X-ray Tomography
For two experiments (cf. Table 1), compression tests were conducted in the X-ray scanner (CT mode) to detect microstructural changes. Thereby, the MCD is mounted upside down in the scanner . Before and after the compression experiment, a CT scan was made from the entire sample at spatial resolution (voxel size) of 25 μm. The time required for each scan was ∼2 hr, while in rack mode the compression was started after a few minutes. Accordingly, in CT mode the snow undergoes enhanced sintering (cf. Supporting Information Text S1 for a comparison).

The 3-D Image Correlation
To characterize local displacements during the compression experiment we employed 3-D image correlation of CT images that were segmented by standard procedures (Hagenmuller et al., 2013). We implemented a multiresolution demons registration method (Thirion, 1998) in C++ using the Insight Toolkit (ITK) library (www.itk.org). The method is detailed in Text S2 in the supporting information. The registration was carried out on subimages with (946 × 946 × 1046) voxels that cover the full height of sample holder.

Stress-Strain Curves
As an overview over all experiments we show one example of the measured stress-strain curves for each strain rate in Figure 1a. Qualitatively, the deformation curves are smooth for low strain rates (dark blue) and recover the results from . At a (lower) critical strain rate we observe a transition to compressive stick-slip where the stress develops quasiperiodic serrations with large amplitudes and low  To quantify the serrations, we evaluated their relative amplitude as the mean height of a stress drop Δ normalized by the mean stress̄of an experiment. The serration frequency was determined from the mean time between two consecutive stress drops. A stress drop was thereby defined as any decreasing part of the stress strain curve after filtering noise close to the sampling frequency, which is an order of magnitude higher than any measured serration frequency. The results in Figure 1b demonstrate the systematic dependence of amplitude and frequency on imposed strain rate ext . For the highest strain rate the characteristic serrations disappear while random stress fluctuations remain.

Microquakes and Strain Localization
Additional insight was obtained from the two experiments (cf. Table 1) where microstructural changes during compression were monitored by CT.
Figure 2 a shows the (horizontally averaged) profiles of the density in the CT sample holder over the entire height before and after the compression. For one experiment (20 μm s −1 ) a high density region developed in the middle of the sample. For the other experiment (40 μm s −1 ), the densification is visible in the height range (0-15 mm), while the density near the piston remains at its initial level. These localized densification patterns are different from the homogeneous densification found at low velocities (0.2 μm s −1 ; cf. . In addition we conducted 3-D image correlation on the microstructure corresponding to the blue sample from Figure 2a. The derived strains are shown in Figure 2b and reveal a band of high compression that is consistent with the density peak in Figure 2a. Despite the uncertainties of the method (cf. Tet S2), the average of the local strains over the entire sample yields 0.115, which is in good agreement with the value of 0.111 that is computed from the piston displacement.

Elastoviscoplastic Model With Sintering
To interpret the experiments we generalize a scalar, mechanical Maxwell model to explicitly account for the evolution of intergranular bond area through competing processes of sintering and bond failure. To this end we consider purely axial deformation along one axis of a Cartesian coordinate system (cf. Barraclough et al., 2017) and assume a linear decomposition of the total strain rate vp into elastic and viscoplastic contributions. Elastic strain is related to stress via Hooke's law = E el with an elastic modulus E.  (6) and (8). We impose a fixed global strain rate . ext , such that the stress acting on the sample evolves according to where in general the tangent modulus E eff , which defines the maximum slope of the stress versus strain curve, is less than E due to machine effects and anelastic relaxation processes (cf. Text S3).
We describe the viscoplastic response by Glen's law in the following form where the rate exponent m ≈ 3 characterizes the viscoplasticity of ice and . 0 is a rate parameter which we set, without loss of generality, equal to . 0 = 1 s −1 . On the microstructure level, stress and viscoplastic strain localize in bonds connecting ice granules. Accordingly, in the expression for . vp , the macroscopic stress is divided by a stress parameter Σ which is proportional to the ratio between effective bond area and specimen cross section. Σ is therefore a monotonically increasing function of t s ∕ where t s is the mean sintering time, which serves as a microstructural state variable, and is a characteristic time for bond formation. In absence of irreversible deformation, the sintering time equals the time elapsed from the moment of snow deposition. During plastic flow, however, the sintering time is reset to zero once a bond is sheared and broken. The mean sintering time then becomes a dynamic variable with evolution equation where s is the characteristic macroscopic strain after which, on average, each bond is broken once (for derivation of equation (3), cf. Text S4).
Formally, equations (1) and (2) represent a nonlinear Maxwell model as employed by , for . ext = 0 and with fixed value of Σ, to fit stress relaxation in new snow samples prepared according to the same procedure as used here. To apply the same model to dynamic testing, it is necessary to specify the sintering function Σ. We determine this function by direct comparison with experimental data. To this end, we note that equations 1-(3) constitute a nonlinear, dynamical system for ( . vp , t s , ) which has the stationary solution The idea is now to compare equation (6) to the strain rate dependency of the flow stress in our deformation experiments which cover a wide range of strain rates . ext . To this end, we average over the serrated stress-strain curves to establish a rate-dependent average stress level. However, before this can be meaningfully done we need to remove the ascending trend of the stress-strain data which is manifest in Figure 1 and arises from the densification which occurs during the compression experiment. Densification can be accounted for by integrating the mass continuity equation for uniaxial compressive straining in the small strain limit, . = . , to obtain = 0 exp( ). Since the elastic modulus and the strength of snow increase with snow density according to E eff =Ẽ eff ( ∕ 0 ) p and =̃( ∕ 0 ) p , we introduce densification-corrected variables viã= where for simplicity we use a common density exponent p = 5 for all stress variables to enclose the high measured values for the elastic exponent (Gerling et al., 2017). Finally, we remove oscillations by averaging the resulting stress variable over the entire stress strain curve. The resultinḡ( . ext ) data are shown in Figure 3 (bottom, blue circles). These data are then divided by ( . ext ∕ . 0 ) 1∕3 to obtain the corresponding Σ values. The function Σ(t s ∕ ) is finally obtained using equation (5) to convert strain rate into mean aging time. It can be seen that Σ increases monotonically with increasing mean aging time (decreasing strain rate). For large strain rates (very short aging times), Σ approaches a constant, and for long aging times (low strain rates) it increases as a power law. In between we find a regime where Σ increases strongly with sintering time. The fit function yields an excellent representation of the experimental data.
The remaining model parameters are determined as follows: The effective modulus E eff is chosen to reproduce the upward slope of the serrations on the experimental stress-strain curves. The viscoplastic exponent m in Glen's law is chosen equal to m = 3, typical of creep deformation of ice. This choice separates the time dependency of deformation that is due to ice viscoplasticity from the time dependency of the snow microstructure contained in the sintering function Σ In summary, we obtain the following physical model parameters: The residual strength at low aging times Σ 0 = 4.8 kPa, the aging time-dependent strength parameter Σ 1 = 72.7 kPa, the long-term aging exponent n = 0.19, the characteristic strain for bond breaking s = 0.00243, the effective elastic modulus E eff = 400 kPa, and the viscoplastic exponent for ice matrix m = 3. The characteristic time is adjusted to correctly reproduce typical serration frequencies, yielding = 4 s.
We note that, for deformation after a sufficiently long sintering time t s , the model predicts that the initial yield stress Σ increases in proportion with t n s with n ≈ 0.19 which is in good agreement with independent sintering studies (Herwijnen & Miller, 2013;Podolskiy et al., 2014)

Numerical Solution
To simulate stress-strain curves we numerically integrate equations (1)-(3) with (8) and (7) and then restore the densification-related "trend" by setting → exp(p ), E eff → E eff exp(p ). This leads to the stress-strain curves depicted in Figure 4 (left). We find good qualitative agreement with the experimental observations. In particular, the increase of the oscillation frequency and the decrease of oscillation amplitude with increasing strain rate are well captured. As in the experiments, the serrations suddenly emerge with large amplitude at a lower critical strain rate . c,l (subcritical behavior) but disappear at an upper critical strain rate . c,u with an amplitude that decreases monotonically to zero, Δ̄∝ (

Discussion
By a thorough analysis of the strain rate dependence we have revealed the rich behavior of snow compression near the DBT with a sudden emergence of stress serrations with changing amplitude and frequency. The analysis sheds new light on the peculiar nature of the DBT in snow. Unlike ice, the DBT in snow cannot be considered as a single transition directly from ductile to brittle behavior but should be envisaged in the terminology of Yosida et al. (1958) as "brittle behavior of the first and second kind" which includes an intermediate strain rate regime of serrated deformation. Since a similar behavior has been observed in shear experiments (Montmollin, 1982), the existence of a "stick-slip regime" at intermediate strain rates may well be a generic feature of snow deformation.
We have derived a minimal, nonlinear Maxwell model with sintering that is able to reproduce the observations. The model reveals that the competition between the viscoplastic strain rate (embodied in Glen's law and its parameters) and the characteristic strain rate s ∕ of the bond renewal process leads to an transition region in the stress-strain rate characteristic (equation (6) and Text S5) where smooth deformations become unstable. Even though the a priori unknown parameters for the bonding state could be convincingly determined by fitting the model to experimental data (cf. Figure 3), it would be appealing to link them with established bonding metrics in snow (Hagenmuller et al., 2014).
"Compressive stick-slip" as studied here is one example of a wide class of phenomena that are discussed, in the context of materials mechanics, under the label of "strain rate softening instabilities" (Zaiser & Hähner, 1997a, 1997b, and in the tribological context, under the label "stick-slip instability." Our model has very close analogies to rate-and-state friction models of block spring systems that are widely used in earthquake modeling, as detailed in Table 2. In particular, the contact aging time as introduced by Dieterich (1979) in rate-and-state friction models fulfills exactly the same relaxation equation as equation (3) if the appropriate substitutions are made. The N-shaped stress-strain rate characteristics used in our model and depicted in Figure 3 corresponds to a "spinodal" friction law as introduced in the context of rate-and-state friction by Estrin and Brechet (1996): It exhibits two stable sliding branches at low and high velocities, separated by a potentially unstable regime of negative differential friction (friction force decreases with sliding velocity) at intermediate velocities. The corresponding nonlinear behavior has been studied by Putelat et al. (2010). In the regime of negative differential friction, steady compression may become unstable with respect to a Hopf bifurcation toward periodic orbits. For the "spinodal" friction law of Estrin and Brechet (1996), the Hopf bifurcation that leads from steady-state sliding to stick-slip is subcritical at the lower critical sliding velocity and supercritical at the upper critical sliding velocity (Putelat et al., 2010). Thus, at the lower critical velocity the periodic orbits emerge with finite amplitude whereas at the upper critical velocity stick-slip oscillations increase monotonically from zero. These analytical predictions match our experimental observations and simulation results: The regime of flow serrations matches the negative slope regime in the stress versus strain rate characteristics, and the serrations emerge with large amplitude at the lowest critical strain rate but disappear gradually at the upper critical strain rate (Figure 1). In the serration regime, an individual serration event corresponds to the nucleation of a single band of localized strain. Analysis of similar models in the context of earthquakes (Erickson et al., 2011) and metals (Lebyodkin et al., 2000) demonstrates a transition from localized bands at low strain rates in the serration regime to propagating bands at higher strain rates. The tomography analysis (Figure 2) indicates such a transition in the nature of the observed "snow microquakes" also in our experiments.
In view of localization as a spatiotemporal phenomenon, we note that the present scalar model can be reformulated in a tensorial manner and solved by the finite element method model to describe 3-D snow deformation under general boundary conditions. Similar features of stick-slip behavior should also emerge in 1-D or 2-D fiber bundle models if viscoplastic fiber effects were included to compete with the bond/fiber strength as dynamical variable (Capelli et al., 2018;Reiweger et al., 2009). And the same may be true when noninstantaneous, relaxational dynamics of bond strengths is implemented into discrete element models of snow deformation (Kabore & Peters, 2019;Mulak & Gaume, 2019). Applying such models to systematically investigate the mechanical behavior under general loading conditions in the vicinity of the DBT may be a fascinating task for future investigations.

Conclusions
We have experimentally and theoretically detailed the mechanisms underlying the ductile-to-brittle transition in snow as a double transition with an intermediate regime of stick-slip, due to microstructural healing, and in formal equivalence with other geomaterials. Notably, at high sampling rate even a barely constrained "squeeze test" of a human-made snowball shows compressive stick-slip prior to fracture (cf. Movie S1). So it will not come as a surprise if microquakes are prevalent but yet undetected events in natural snowpacks.