A grain-size driven transition in the deformation mechanism in slow snow compression

Compression of snow at low strain rates governs the natural densification of snow and firn on Earth. Different deformation mechanisms (e.g., grain boundary sliding and dislocation creep) are widely employed in models and are mainly believed to change depending on relative density or snow type. To explore this picture, we conducted compression experiments with a nominally constant strain rate ̇𝜀 ≈ 10 −6 s −1 and systematically varied snow type, density, and specific surface area. We used a micro-compression stage to enable X-ray micro-tomography imaging and microstructure characterization before and after the experiment. Using repeated load-relaxation cycles, we eliminate changes in the microstructure that occur during the first loading, allowing us to probe the viscoplastic behavior of the intact ice matrix. To evaluate the effective mechanical behavior of snow, we derived an implicit, analytical solution of a non-linear scalar model that characterizes loading and relaxation in terms of the stress exponent 𝑛 in Glen’s law. Our results evidence a sharp transition of 𝑛 as a function of the geometrical grain size, and rule out that this transition is caused by density or snow type. The derived values of 𝑛 ≈ 1 . 9 for small grains to 𝑛 ≈ 4 . 4 for large grains are consistent with a transition from grain boundary sliding to dislocation creep at a transition grain size of around 0.5-0.6 mm. Our results shed new light on the old discussion about deformation mechanisms in snow and firn.


Introduction
Slow creep of snow or low-density firn is a key process for understanding densification of seasonal (alpine) or perennial (polar) snowpacks.Accurate predictions of densification rates have important applications.In paleoclimatology, firn densification models are required to interpret climate signals of entrapped gases (e.g.[1]).In natural hazards, snow densification determines the mechanical stability of layers for forecasting avalanches (e.g.[2]).Yet a clear picture of contributing or dominant deformation mechanisms underlying slow creep in snow and firn is still lacking.
A characteristic, though not unique and therefore widely debated, signature of a deformation mechanism is the stress exponent characterizing the (power-law) creep rate as a function of stress.For icy materials, this relation is commonly termed Glen's law [3] and widely used to characterize viscoplastic rheologies of ices in lab experiments and in glacier flow [4].Glen's law of the form ε =   relates the strain rate ε to the stress  where the value of the stress exponent  varies depending on the deformation mechanism.
For snow and low-density firn, it is widely accepted that the relevant deformation mechanism is grain boundary sliding (GBS).This was originally put forward by Alley [5] based on the theory from Raj and Ashby [6] for diffusion-accommodated GBS.Accordingly, the stress exponent attributed to GBS is  = 1.With increasing density, deformation is argued to transition to a dislocation creep (DC) regime with  ≈ 3 − 4 [3,5].The transition from GBS to DC is believed to occur at a relative density of around ≈0.6, corresponding to a random close packing [7].There, particle rearrangements essentially cease to exist, although it was pointed out [8,9] that particle rearrangements may be still effective at even higher densities.As a consequence of [5] nowadays physical models of snow and firn densification [10][11][12] often include transitions from  ≈ 1 to  ≈ 3 − 4 at the snow-firn transition density.However, experimental evidence is still lacking that a transition from GBS to DC (or at least a transition in exponents) actually exists in snow [13].
A considerable number of results implicitly support or contradict the picture from Alley [5].On the supporting side, the GBS-based model could be successfully optimized to match firn density data [13].
In triaxial experiments on Alpine snow, a crossover of  ≈ 1.8 to  ≈ 3.5 as a function of density was found [14].Creep experiments complemented by microstructural analysis consistently found values between  ≈ 4.1 and  ≈ 4.6 for firn densities well above the snow-firn transition [15].A value of the stress exponent of  ≈ 2 https://doi.org/10.1016/j.actamat.2023.119359Received 15 May 2023; Received in revised form 15 September 2023; Accepted 18 September 2023 was found [16] for low-density snow in compression.On the contradictory side, the crossover in Scapozza and Bartelt [14] occurs at a relative density of ≈0.3, well below the close packing density.The well-known increase in the effective viscosity of depth hoar [17,18] occurs without any considerable changes in density.Values of  ≈ 3 were used for low-density snow and yield consistent results in three-dimensional numerical homogenization [19] or one-dimensional models [20].And lastly, it was shown [21] that, as an alternative to GBS, anisotropic (mono-crystalline) material behavior with  = 3 also leads to a reasonable agreement with experiments.
While snow and firn densification are certainly affected by concurrent microstructural processes (sintering, metamorphism) in the evolving ice matrix [22], snow, as porous ice, must inherit viscoplastic properties from the ice matrix [23].The stress exponent in polycrystalline ice has been widely and controversially investigated since [24].From laboratory experiments conducted on polycrystalline ice at a stress range between 0.01 to 10 MPa, it is well-known that the value of  generally changes with stress level [25][26][27], with  ≈ 2 at low stresses and  ≈ 3 − 4 at high stresses.Observations of natural ice flow [28][29][30][31][32] also consistently yield stress exponents in the range of 1 to 5. For polycrystalline ice, the transition from GBS to DC was emphasized by Peltier et al. [33] and Goldsby and Kohlstedt [26] involving the concept of superplasticity [34].Accordingly, the DC relation in ice is characterized by a stress exponent of  ≈ 4 and is independent of grain size, while GBS (dislocation-accommodated) [26] is characterized by  ≈ 2 and is strongly dependent on grain size [26,35].This grain size dependence facilitates the reconciliation of lab and field data [36] in explaining the range of naturally measured exponents as an apparent outcome of the underlying composite rheology.On the other hand, superplastic flow from [26] was also critically discussed [37,38] due to potential inconsistencies with fabric evolution in natural ice.
Overall we are left with a widespread controversy about the existence of and drivers for transitions in the stress exponent in snow and ice, and implications for the underlying deformation mechanisms.For snow, systematic attempts at characterizing  are virtually absent in the literature.The main difficulty is that access to small-strain, viscoplastic properties of the ice matrix in snow compression experiments are usually concealed by rheological variability associated with changes in snow microstructure.Compression experiments can be conveniently conducted nowadays in micro-compression stages combined with micro-computed tomography to link the mechanical response to the microstructure [15,16,18,20,[39][40][41][42][43][44].However, all studies thus far have focused on the evolving microstructure during initial loading.This prevents disentangling the effects of micro-failure and metamorphic processes from viscoplasticity in the ice matrix.According to [45,46] this problem can be avoided if a sample has undergone preliminary loading before testing.
Thus, complementary insight can be gained by conducting microcompression experiments with repeated load and relaxation cycles [45].This allows us to access the viscoplastic response of the intact ice matrix.Our paper aims to make a substantial contribution to the controversial discussion of stress exponents in snow and ice through systematic experiments on the viscoplastic properties of the ice matrix in snow.We explore the influence of the (intact) porous microstructure at low strain rates and small strains.We achieve this through consecutive loading-relaxation experiments to eliminate structural transients (due to microfailure, for example) combined with X-ray micro-tomography (micro-CT) imaging.We derive a novel solution of a nonlinear transient model involving generic power law creep.This allows us to estimate effective ice mechanical parameters from loading and relaxation cycles within a single framework.Due to a careful selection of snow samples of different types (rounded, faceted, depth hoar, and melt forms), relative densities (0.18-0.43) and geometrical grain sizes (0.2 mm-1 mm), we reveal for the first time a clear transition of the stress exponent in snow from  ≈ 1.9 to  ≈ 4.4 that is driven by (geometrical) grain size.
This finding contradicts the previously hypothesized dependence of the value of  on density or snow type.
The manuscript is organized as follows.In Section 2, we describe the snow samples and micro-tomography setup, and detail the parameter estimation through the transient model and data from repeated loading relaxation experiments.The mechanical data and micro-structural results are shown in Section 3. Our discussion in Section 4 will attempt to reconcile the variety of existing results discussed above.

Snow samples and sample preparation
For our experiments, we chose natural Alpine snow samples that were collected in the field near Davos, Switzerland in the form of large snow blocks during the 2021/2022 and 2022/2023 winter season and stored at −24 • C in a cold laboratory.Since the blocks were not collected at the same time, storage times in the cold laboratory prior to mechanical testing varied from days to months, whereas the protocol for the tests (see below) was exactly the same.Samples were selected to include a variety of combinations of density, specific surface area and snow type.In total 11 experiments were conducted on samples with microstructures as shown in Fig. 1.The relative density (or volume fraction of ice)  ranged from 0.18 to 0.43.As a geometrical grain size metric, not to be confused with the crystallographic grain size [47], we used the optical equivalent diameter  opt = 6∕SSA which is unambiguously defined through the inverse of the surface area per ice volume, SSA, [48].Values range from 0.28 mm to 0.97 mm (see Table 1).The microstructures of the snow samples included rounded grains (RG), faceted crystals (FC), depth hoar (DH), and melt forms (MF) which were visually classified from the 3D tomography images following the international classification for seasonal snow on the ground [48].Since the four MF microstructures exhibit clear differences (see Fig. 1h-k), we included the sub-types (MF(RG) and MF(DH)) [48].The unusual combination of MF snow types with rather small grains reflects seasonal particularities at the sampling site where the snowpack was temporarily isothermal during a warm period in the early season.For visual distinction throughout the paper, we chose marker or highlight colors in figures (such as in Fig. 1) that allow to discern the different snow types (RG, FC, DH, and MF) by their standard color defined in the classification [48].
The preparation of each specimen for the compression experiments involved the following steps.Before the experiment, a rectangular specimen of snow was extracted from the stored block and leveled to create a smooth top and bottom surface, and then carefully cut into a cylindrical shape using a band saw.To obtain homogeneous samples, we aimed at small specimen heights to limit the vertical variability of the natural snow in the sample.The resulting sample geometries (approximately 20 mm in height and 36 mm in diameter) have therefore a considerably smaller height-to-diameter ratio as the firn tests from [15].The sawed samples were then carefully slid into the cylindrical holder with 36 mm diameter.The sample holder is made of polyoxymethylene (POM) to minimize the friction in contact with the snow sample.The true contact area between snow sample and sample holder is limited by the sawing accuracy, which is difficult to assess and therefore constitutes an uncertainty in the boundary condition of our experiments.
After preparation of the specimen, each sample was stored at −14 • C for around 1 h for thermal equilibration.The sample holder was then mounted upside down in the micro-compression device (MCD) and installed in the micro-CT scanner where they resided for the entire test. .See text and Table 1 for the description of the snow samples.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Cyclic micro-compression experiments
Displacement-controlled, uniaxial multi-cycle micro-compression experiments were conducted on snow samples using the MCD developed by Schleef et al. [16].All experiments were conducted at a temperature of −14 • C. Each specimen underwent four cycles (I, II, III, IV) of loading and relaxation.During the loading phase, the sample was compressed with the lowest possible piston velocity of 0.1 μm s −1 until the maximum force of approximately 13 N on the force sensor was reached.The motor was then stopped and the stress relaxation was recorded.We held the sample at a constant displacement for roughly 4 h for the first I and II relaxation cycles and 8 h for the III and IV relaxation cycles.An example of the data and the cycles is shown in Fig. 2 for illustration.From the experiments, we obtain the nominal stress (force per piston area) and the nominal strain (ratio of sample height at time t to initial height) as a function of time.The nominal (imposed) strain rate is obtained from the piston velocity and the nominal sample height according to ε = ∕20 mm leading to a nominal strain rate of about 5 × 10 −6 s −1 .Note that the initial sample heights varied leading to an average strain rate of 6.05×10 −6 s −1 for all samples.
The microstructure of the snow samples was characterized by microcomputed tomography using a SCANCO CT80 scanner.The scanner is equipped with a temperature control to balances the heat production of the X-ray tube.A temperature increase of ≈1 • C during the scan cannot be avoided, though.To minimize temperature fluctuations, we have conducted X-ray scanning only before and after the mechanical test for monitoring the microstructural changes.To this end, the same  The scans took ≈2 h to complete.For the evaluation, cubic images with 412 × 412 × 412 voxels were selected and processed using established segmentation procedures [49].Snow microstructural properties (ice volume fraction and geometrical grain size before and after the experiment) were obtained from the binary images by calculating the derivative of the two-point correlation functions at the origin.Table 1 gives a summary of all experiments with snow type, volume fraction of ice , optical grain size  opt , and nominal total strain  in each cycle.

Non-linear Maxwell model
In order to analyze the transient loading and relaxation stages of the experiments in view of an effective, non-linear, stress-strain rate relation (Glen's law), we derive a novel analytical solution of a scalar Maxwell model involving the physical parameters of the problem.
To this end, we employ the common superposition of total strain rate ε into an elastic εel and viscoplastic part εvp according to For the viscoplastic strain, the widely used form of Glen's law εvp =   is inconvenient for parameter estimation since the physical units of the prefactor depend on the exponent .Therefore we cast Glen's law in the form In this form, the characteristic stress scale  * can be chosen arbitrarily for non-dimensionalization, while the time scale  vp can be interpreted as the compactive viscosity.The exponent  reflects the non-linearity in Glen's law and should be indicative of the deformation mechanism.The elastic strain rate in Eq. ( 1) is given by εel = σ∕, with E being Young's modulus.Introducing a normalized stress by  = ∕ * and defining the dimensionless parameter  = ∕ * , the strain superposition leads to the ordinary differential equation that applies to both the loading and relaxation stages.

Relaxation phase
For the stress relaxation phases, we have ε = 0, and the ordinary differential equation (ODE) can be integrated, yielding the analytical solution previously used in Schleef et al. [16] as: with a relaxation time  =  vp ∕.

Loading phase
In the loading phase, we have ε ≠ 0, and a closed-form expression for () cannot be found.However, the inverse function () is analytically accessible and sufficient for estimating the experimental data's physical parameters in monotonic regions.Separation of variables of Eq. (3) leads to with the dimensionless parameter  ∶= 1∕( ε vp ).This equation can be integrated from initial time  0 to final time  with corresponding stress values  0 to  and written in closed form according to Here (, , ) denotes Lerch's transcendent [50] which is numerically accessible via lerchphi in the python sympy library.Thus Eq. ( 6) provides us with a closed-form expression for the inverse function () with ( 0 ) =  0 which is sufficient for the estimation of the physical parameters.
By computing the relaxation time scale  from the loading fit parameters  and , the consistency of the parameters obtained from independent loading and relaxation fits can be verified by

Parameter estimation
The first cycle of the experiment is subject to large strains due to structural relaxations.On average, the first cycle contributes to 75.65% of the total strain of the sample (see Fig. 2).Therefore disregarding the first cycle allows us to access the mechanical behavior of the ice matrix without microstructural changes in the II, III, and IV cycles.Only for these cycles we expect our model to be applicable.For parameter estimation, we separate the stress curves into loading and relaxation curves (see Fig. 2).For each cycle considered, the loading phase begins at the minimum stress (represented by the red dots) and ends at the maximum stress (represented by the green stars).Green stars mark the start of the relaxation phase which continues until the red dot.After separation, the loading and relaxation curves were fit with respective fit functions, as described in Section 2.3.1.
In principle, the loading fit will involve three free parameters (, , and ), and the relaxation fit two parameters,  and .Since the loading phase only involves very slight deviations from a straight line, an unconstrained optimization involving all three parameters is poorly constrain.To this end, we prescribe the parameter  in the loading function Eq. ( 6) with the value obtained from the fit of the subsequent relaxation curve.In the loading function, care has to be taken for the argument  since (, , ) has a branch point at  = 1.To reside on a unique (continuous) branch we choose  * as the maximum stress and enforce the dimensionless parameter  < 1 as a parameter bound in the least squares optimization.We include the full range of the stress curve for the loading fit function, while for the relaxation we always restricted the fit to the first 2 h of the relaxation phase.
The uncertainties of the fit parameters were characterized by standard errors that were calculated from the covariance matrix of the fit provided by the minimize method from scipy in Python.The uncertainty of derived parameters is calculated through error propagation.In addition, the cost function (sum of squared residuals) was evaluated in the neighborhood of the optimum to visualize differences in the minimization problem.

Integrity of the microstructure during compression
An important assumption of our experimental approach and the subsequent analysis is that the ice matrix largely remains intact during compression.This is visually confirmed by Fig. 1, where the same  slice in the scanned volume before and after the experiment is shown for each sample.Corresponding microstructural features can be easily recognized and the displacement of the grains in vertical direction is attributed to the strain in the first compression cycle.The similarity of the microstructures before and after the test is further consolidated by quantitative analysis of the 3D images.Fig. 3 shows the comparison of the microstructural parameters  and  opt for all samples before and after the experiment, indicating hardly any changes (agreement with  2 = 0.961 and  2 = 0.999 respectively).This implies that structural changes occurring during the first cycle were consistently restricted to the top and bottom boundary of the sample container where the microstructure is weakened through the preparation procedure [45] and therefore prone to failure upon first loading.The largely unchanged microstructure in the center of the sample (cf.Fig. 2) justifies an interpretation of our model (fit) parameters as the viscoplastic response of the intact ice matrix.

Estimation of mechanical model parameters
Fig. 4 shows the normalized stress  as a function of the rescaled time  ε for the three evaluated loading and relaxation cycles (II,III,IV) together with Eqs. ( 4) and ( 6) fitted to the data.For clarity, only one example per snow type is shown here.Overall the visual agreement of the model with the data is very good for both loading and relaxation.The contour plots of the cost function (cf.Fig. 5) reveal the differences between the fitting for loading and relaxation: In loading, the model Eq. ( 6) must essentially constrain the slope of a straight line (cf.Fig. 4(f,i)) emanating from a given start point ( 0 ,  0 ) with two parameters.The contour plots confirm this very weakly constrained situation which leads to generally higher uncertainties in the loading phase.This is different when the stress attains the non-linear region (cf.sample C0010107 in Fig. 4(a), cf.Fig. 5(a)) Applying the fitting from Section 3 consistently to all 11 samples and all three cycles (II, III, IV) we obtained the set of mechanical relaxation parameters ,  and loading parameters ,  vp for further analysis.

Consistency of loading and relaxation time scales
Before turning to the main result, we demonstrate the consistency of the fit parameters in Fig. 6.First, we compared the parameter  obtained from the loading fit (as described in Eq. ( 7)) with the estimate  from the relaxation fit (as described in Section 2.3.1).On the one hand, we observe a fair agreement between these independent estimates (cf.Fig. 6a with  2 = 0.79).This correlation confirms that the stress evolution for subsequent loading and relaxation phases can be consistently described within the same model (see Eq. ( 3).On the other hand, we also find high correlations between corresponding estimates of  from different cycles ( 2 = 0.79 for the loading phase in Fig. 6b) and  2 = 0.96 for the relaxation phase in Fig. 6c).The correlation between consecutive cycles tends to increase with cycle number indicating a ''convergence'' of the response.This is confirmed by the total strains in the last two compression cycles, which are virtually identical (cf.Table 1).

Linking stress exponent to the snow microstructure
Finally we examine the relationship between the stress exponent  (Fig. 7a-b) and the time scales  and  vp (Fig. 8a-d), and the microstructural properties  and  opt of the samples.
For the stress exponent, we find essentially two clearly separated groups of values,  ≈ 1.9 and  ≈ 4.4 (Fig. 7a-b).Neither density nor snow type are able to explain these two groups.In contrast, a systematic behavior of the exponent is revealed when plotted as a function of grain size (Fig. 7b).The stress exponent exhibits a clear transition from one group to the other within the narrow range 0.5 <  opt < 0.6.This transition is independently visible for each cycle, and the estimated difference in  clearly exceeds the uncertainties.The sample C0011067 is thereby located in the transition region with  ≈ 3.4 mm.The considerably larger error bars for sample C0010102 in all cycles suggest to regard this sample as an outlier.Neglecting this outlier, the time scale  shows neither a clear dependence on density  nor on grain size  opt (Fig. 8a and b).In contrast, the time scale  vp increases with increasing  but lacks a clear trend with  opt (Fig. 8c and d).This density dependence of the time scale  vp is expected since it is essentially proportional to the effective viscosity.In contrast, the time scale  is essentially the ratio of the effective viscosity and an effective elastic modulus.Both are known to increase with density in snow which should cancel, or significantly reduce, the density dependence of .
The considerable differences in error bars in Fig. 8 are explained by the fact that  is obtained from relaxation while  vp is obtained from loading where the uncertainties are higher due to the reasons explained above.

A grain size driven transition in Glen's law for snow
Motivated by old work [45], we made a novel attempt to probe the viscoplastic response of the ice in snow using state-of-the-art instrumentation (micro-compression stage and tomography imaging).By repeated loading and relaxation cycles, we could overcome mutually influencing viscoplastic and micro-failure processes in the first load cycle (reflected by the large strains (cf.Table 1)) and probe the intact microstructure (Fig. 3) in the interior through slow, small-strain (on average <1%) compression relevant for natural densification processes.Effective mechanical parameters for this process were obtained through a novel implicit solution of a non-linear, scalar model based on Glen's law.This led to an excellent visual agreement with the data (Fig. 4).Various consistency arguments were given for the comparison of the derived parameters between cycles (Fig. 6b-c), between loading and relaxation (Fig. 6a), and their dependence on microstructure (Figs. 7  and 8).With this approach, we found for the first time a marked transition in the effective exponent in Glen's law for snow driven by geometrical grain size (Fig. 7b).This finding is robust, given the uncertainties of the method.This result will be argued to be consistent with results in ice mechanics and suited to reconcile several old results in snow mechanics.

Interpretation as a transition from GBS to DC
The apparent transition in  as a function of grain size must not be confused with an explicit dependence of the exponent on grain size.It rather confirms the interpretation of a composite rheology [34,36,51,52].When the different deformation mechanisms (GBS, DC) are active simultaneously, the total (viscoplastic) strain εvp = εGBS + εDC is an additive superposition of the respective strain rates.The individual contributions have different stress exponents, which are reported to be  ≈ 1.8 for GBS and  ≈ 4 for DC in ice [35].The GBS rate constant depends on grain size while the DC rate does not [35], which changes the relative contribution of each mechanism as a function of grain size.Suppose this situation is interpreted within a single, non-linear stressstrain rate relation, as usually done by using Glen's (single power) law as in Eq. ( 2).This leads to an apparent (effective exponent) that transitions from the GBS value to the DC value [26,36].This is precisely what we also observe here for snow.
The transition between GBS and DC in ice depends on stress and strain rate (and thus temperature) [36,52] such that the transition grain size decreases with increasing stress.The nominal volume-averaged stress (around 0.01 MPa), the nominal volume-averaged strain rate (6.05 × 10 −6 s −1 ), and the temperature (−14 • C) are constant in our experiments.Deviations of true stresses and strain rates from the mean at the micro-scale in the porous structure are evident, though.Under these conditions, we consistently fall into the low-stress regime and found that the transition is in the 0.5-0.6 mm range.In snow, no comparison data exists.In ice, an estimate of 0.2-0.3mm for the transition grain size was given [36] for higher stresses, which appears to be qualitatively consistent with our result.For this comparison, we stress that the relevant grain size in ice deformation is the size of the crystallographic grains, which must not be confused with the geometrical grain size [47] (optical equivalent diameter) used here due to the lack of crystallographic orientation measurements.In a first approximation, both grain sizes are of the same order of magnitude as in snow as suggested by texture images from snow [53] and firn [54].

Reconciling strain rate peculiarities in snow
Our experiments characterize the transient flow stress covering time scales of hours to days; thus, an extrapolation to longer times cannot be directly justified.Despite this caveat, that similarly applies to all lab experiments, our results strikingly reconcile observations beyond the well known super-linear increase of the effective viscosity with density [55,56].This is consistent with the increase of  vp from Fig. 8c since the effective viscosity is proportional to  vp by virtue of Eq. ( 2).
A striking feature in snow is the creation of depth hoar (DH) under temperature gradients, which exhibits strongly decreased strain rates compared to RG [17], even if temperature, stress, and relative density are approximately the same [57].DH is said to have ''almost always'' a specific surface area of 10 m 2 /kg, which translates to an optical grain size of 0.65 mm.This is directly above but close to the transition regime observed here.More specifically, [18] observed a considerable decrease of the strain rate in temperature gradient metamorphism, despite the small changes in SSA values between 13.2 mm −1 and 11.7 mm −1 .The conversion to optical grain size yields the range 0.49 <  opt < 0.56, i.e. located exactly in the transition region.This provides an explanation for the sensitivity in the measured strain rates in [18].
Attributing differences in snow mechanical behavior to the categorical variable ''snow type'' is still widely adopted, though unphysical.In [58], differences in the effective exponent of  = 4.1 for faceted compared to  = 3.2 for rounded grains at −15 • were also attributed to snow type.However, grain sizes of 0.2-0.3mm were only characterized visually by a handheld microscope.An explanation of the differences based on optical grain size seems likely.In addition, our results (Fig. 7b) comprise two samples, C0011059 and C0011067, of the same type (MF/DH), with clearly different exponents, ruling out a spurious dependence on snow type.
We are not the first to show that the simple adoption of Glen's law with an exponent  = 3 is incompatible with slow snow compression.In [21], creep experiments were conducted with two RG snow samples with different optical diameters of  opt = 0.28 mm, and  opt = 0.41 mm, thus well in the regime attributed here to GBS.The reported disagreement of the experiments with simulations using Glen's (isotropic) law with  = 3, and the improvement by incorporating anisotropy, should therefore be revisited in view of an alternative explanation involving GBS.In contrast, the exponent of  ≈ 2 found in [16] (optical grain size  opt = 0.18 mm) are consistent with the present results.
In [14], variations of the effective stress exponent were attributed to density.Our data rule this out.First, in contrast to [14], our high-density samples have the lowest exponents, and second, the two exponent groups clearly overlap in density (Fig. 7).This was facilitated by including samples with untypical combinations of high density and small optical diameter.Since optical diameter and density are essentially correlated in nature (the denser, the older, the coarser), a transition in the deformation mechanism that is physically driven by grain size may be easily misinterpreted as a transition driven by density.

Implications on firn modeling
Correlations between optical grain size and density may also cause mutually concealing impacts in firn densification since the introduction of GBS by [5].
X-ray tomography of shallow firn cores from different polar sites show that the optical grain size spans a wide range of values between 0.4 mm and 1.6 mm [59].For all these sites, the relative density was below the close packing density of 0.6.If these sites were to be modeled with common firn densification models [1,12,60,61] the same deformation mechanism (i.e. the same ) would be employed, in contradiction to our findings.Despite leading to consistent results, [13] emphasized the missing evidence for the transition from GBS to DC.Our results seem to provide such evidence, however, with slightly different ingredients than commonly adopted.First, we expect the transition to be driven by grain size driven and not by density.Thus, macroscopic constitutive models that do not involve a composite flow law for firn, where the transition would be automatically taken care of [36], may miss independent grain size and density controls.Second, our stress exponents are more in favor of dislocationaccommodated GBS discussed by [26], rather than the linear behavior from the diffusion-accommodated GBS [6] employed by [5].
In contrast, our experiments are in line with unconfined creep experiments with X-ray tomography characterization conducted in [15].By analyzing the relationship between the stress at the minimum strain rate for total strains between 0.7% and 7.1%, the authors consistently found stress exponents between 4.1 and 4.6, for Greenland firn samples extracted from a depth of 50-80 m with ice volume fractions between 0.82 and 0.93 and optical diameters between 2.4 mm and 4.6 mm.These optical diameters lie well above the grain size transition found in the present work.Therefore, exponent values and interpretation in terms of dislocation creep as the dominant deformation mechanism given in [15] is in agreement with our results.However, our findings suggest to attribute this dominant deformation mechanism to the grain size rather than the density.For high-density firn, repeated loadingrelaxation experiments, similar to ours, have also been conducted by [62].Parameter estimation from that data using the transient model developed here would be straightforward and insightful to further link different density and grain size regimes.

Limitations and perspectives
Our experiments and analysis can only provide indirect evidence on deformation mechanisms from the macroscopic mechanical response.The macroscopic behavior should inherit several dependencies from the micro-scale constitutive law (of ice) [23], which are modulated by geometrical influences from the microstructure.A similarity in exponents does not necessarily imply that GBS is the dominant mechanism [63] and direct grain boundary analysis [64] or surface observation methods [65] are not possible in the interior of porous snow.It is thus necessary to further constrain additional parameter dependencies (e.g.activation energies through temperature variations) to tighten the link to ice mechanics and involved deformation mechanisms.Another candidate parameter is the GBS grain size exponent [34].With additional measurements in the small grain size regime, it seems possible to explicitly resolve a grain size dependent prefactor and estimate a grain size exponent for consolidating the present findings.This would require macroscopic mechanical models for a composite (GBS-DC) rheology, which have been similarly employed for ice [36] or materials like Tie6Ale4V [65].However, a formulation of the macroscopic response using mixture rules solely based on volume fractions [65] is certainly not applicable to the complex, porous microstructure of snow.
For developing advanced models it seems crucial though to improve the characterization of the experimental boundary conditions.For the present analysis within the simple, scalar model this is less critical, since boundary conditions do not enter explicitly.This will be different for a comparison to three-dimensional mechanical models.A promising approach to bypass the poorly constrained experimental situation was suggested by [66] where digital volume correlation on tomography images before and after the compression can be exploited to supply the relevant boundary conditions for a model.

Conclusions
We have shown how the viscoplastic behavior of the (intact) ice matrix in snow can be experimentally accessed in deformation-controlled compression experiments without being concealed by predominant micro-failure processes in the porous material.Through the analysis within a minimal model for the effective, non-linear mechanical response, we demonstrated that snow undergoes a transition in Glen's law with  ≈ 1.9 to  ≈ 4.4.Due to the careful selection of snow samples comprising various combinations of density, grain size, and snow type, we have shown that this transition is driven by the geometrical (optical) grain size and not by relative density or snow type.By analogy with corresponding results in ice, we interpret the transition as a transition in deformation mechanisms, from grain boundary sliding (fine-grained snow) to dislocation creep (coarse-grained snow).Despite working with our setup's slowest possible strain rates, we acknowledge the remaining gap in stress and time scales between our experiments and natural conditions.But our results clearly suggest that the geometrical grain size must be taken into account when formulating macroscopic densification rates of snow and firn.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. Example of the MCD experiment showing stress and strain as a function of time during four consecutive loading-relaxation cycles.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Comparison of initial and final relative density  (left) and geometrical grain size  opt (right)from the CT images before and after the experiment.

Fig. 4 .
Fig. 4. Normalized stress  as a function of rescaled time ε data for the II (column 1), III (column 2), and IV (column 3) loading and relaxation cycles for one sample each (rows) of RG, FC, DH, and MF together with the model fits, represented by white line.

Fig. 5 .
Fig. 5. Contour plots of the cost function in (A) loading and (B) relaxation for the data shown in Fig. 4 as a function of the free parameters.White dots represent the optimal parameters, rows correspond to samples, and columns to cycles (II, III, IV).

Fig. 6 .
Fig. 6.(a) Comparison of  from loading fit with  estimated from relaxation for cycles II, III, and IV.(b) Comparison of  from loading cycle III with  from loading cycles II and IV.(c) Comparison of  (s) from relaxation cycle III with  relaxation cycles II and IV.The marker sizes indicate cycle number and increase from II to IV.

Fig. 7 .
Fig. 7.The estimated stress exponent  from the relaxation fit as a function of volume fraction  (a) and geometrical grain size  opt (b) for all loading cycles.The black dashed line represents the grain boundary sliding value  = 1.8, and the magenta dashed line represents the dislocation creep value  = 4 from [26].The marker sizes indicate cycle number and increase from II to IV.The error bars are the standard errors computed in the optimization.

Fig. 8 .
Fig. 8. Link between mechanical time scales and microstructure: Dependence of the relaxation time scale  on relative density (a) and geometrical grain size (b) and dependence of the loading time scale  vp on both parameters (c,d).The marker sizes indicate cycle number and increase from II to IV.The error bars are the standard errors computed in the optimization.

Table 1
Overview of the experiments/samples with microstructural properties ( and  opt , snow type) and nominal strains for each cycles (I-IV).
opt (mm)  opt (mm)  (%)  (%)  (%)  (cylindrical region of interest (height 10.3 mm, full diameter) in the middle of the sample holder was scanned with a voxel size of 25 μm.