On the compressive strength of weak snow layers of depth hoar

Abstract

We determined the compressive strength of weak layers of faceted crystals and depth hoar using artificially grown samples with a wide range of microstructural morphologies in a cold laboratory setup. Micro-computed tomography (µCT) imaging showed that the microstructures of the artificial samples were comparable to that of natural depth hoar. We performed compression experiments in a displacement controlled testing machine on 92 depth hoar samples with densities ranging from 150 kg m⁻³ to 350 kg m⁻³. The compressive strength spanned two orders of magnitude (1–150 kPa) at strain rates of about 10⁻³ s⁻¹ at $-5^{\circ}\mathrm{C}$ and followed a power law as a function of density. Several microstructural metrics such as the specific surface area, connectivity density and correlation lengths obtained from µCT measurements exhibited a statistically significant relationship with compressive strength. Analysis of the residuals of the power law fit showed that in addition to density, horizontal correlation lengths also correlated with strength. However, in this study, density remained the dominant predictor of the compressive strength of depth hoar.

1. Introduction

Slab avalanches initiate when a weak layer buried below cohesive slab layers fails (e.g. Reuter and others, 2015). Various types of weak layers with a wide range of microstructures are present in seasonal snowpacks, which originate from metamorphic processes, precipitation events or water vapor transport from the atmosphere. In the Alps, the most common persistent type of weak layers are faceted crystals and depth hoar (Schweizer and Jamieson, 2001). These crystals form as a result of nonequilibrium metamorphism of snow in the presence of a large temperature gradient. The resulting microstructure depends mainly on the temperature gradient, the average temperature, the snow density and time (e.g. Staron and others, 2014; Bouvet and others, 2023).

Slab avalanches are often triggered by the victims themselves. They are, therefore, responsible for the majority of avalanche fatalities, which makes them most relevant for avalanche forecasting (e.g. Schweizer and Jamieson, 2001). Since the processes involved in slab avalanche release are well known, deterministic models are expected to support avalanche forecasting in the future. These models require detailed knowledge of the mechanical properties of snow and in particular of weak layers. Snow properties strongly scale with density but can differ by orders of magnitude depending on the microstructure for identical snow densities (e.g. Mellor, 1975; Fukuzawa and Narita, 1993; Sundu and others, 2024b).

The mechanical properties of snow—particularly of rounded grains—have been studied since the 1940s. Snow is highly porous, has a complex microstructure and is a fast-changing material due to its high homologous temperature (>0.9) under natural conditions. It exhibits anisotropic, temperature and strain-rate dependent material behavior (e.g. Mellor, 1975). After settling on the ground, it immediately forms a porous matrix of sintered ice crystals and air (Kirchner and others, 2001). Consequently, the physical properties of snow are closely related to those of ice but strongly influenced by its microstructure. Due to its high porosity, snow can fail under (macroscopic) purely hydrostatic stress states in a volumetric, collapse-like manner. This failure can then propagate through the material like a mode I opening crack, which has been termed ‘anticrack’ (Heierli and others, 2008).

Slab avalanche release requires the failure of a weak layer, hence the material strength must be exceeded. Experimentally determining the strength of different weak layers is challenging due to the seasonal availability and the fragile nature of these layers. As shear strength of weak layers was found to be lower than the compressive strength, several field measurement campaigns aimed to determine this property for different weak layer types (e.g. Perla and others, 1982; Fukuzawa and Narita, 1993; Schweizer, 1998; Jamieson and Schweizer, 2000; Jamieson and Johnston, 2001). Only few studies looked into the compressive strength of snow (e.g. Jellinek, 1959; Kinosita, 1996; Fukue, 1979), but comprehensive data on the compressive strength of weak layers are missing. However, even in flat terrain, the macroscopic stress state in a weak layer is multiaxial due to the bending of the overlying cohesive slab (e.g. Heierli and others, 2008; Rosendahl and Weiß graeber, 2020; Adam and others, 2024). To model avalanche initiation, it is, therefore, necessary to know the full failure envelope from pure compression to pure shear stress states. Based on experimental results (e.g. McClung, 1977; Chiaia and Frigo, 2009; Nakamura and others, 2010), a Mohr–Coulomb criterion was proposed. This criterion describes material strength as a linear function of the ratio of shear to normal stress and is commonly used for geotechnical materials. As a result, the strength continuously increases for higher compressive stress states. However, since snow can fail in pure compression, the failure envelope must be closed. Accordingly, Chandel and others (2014) proposed an elliptical failure criterion based on shear experiments with dead weights. Reiweger and others (2015) performed multiaxial strength experiments on natural weak layers with a controlled loading rate and suggested a modified Mohr–Coulomb criterion by introducing a cap function. In recent numerical approaches based on FEM and DEM, elliptical failure envelopes were used (e.g. Chandel and others, 2015; Mede and others, 2018; Mulak and Gaume, 2019; Bobillier and others, 2020; Ritter and others, 2020; Singh and others, 2022).

Performing mechanical experiments on natural weak layers is challenging. If experiments are done in situ in the field, many parameters (e.g. the temperature) cannot be controlled and transporting the fragile samples to a cold laboratory is difficult and time-consuming (e.g. Lombardo and others, 2021). In both cases, the results are limited to the naturally available microstructures and densities, and the experiments are dependent on the seasonal availability of weak layers. Most of these challenges can be avoided by the artificial growth of weak layers, which is why artificial depth hoar weak layers have been used for experiments for many years (e.g. Fukuzawa and Narita, 1993; Reiweger and others, 2010; Reiweger and Schweizer, 2013; Capelli and others, 2018). A number of different setups have been described in the literature. The more recent publications use weak layers grown in situ between two dense snow slabs, providing a good interface for the subsequent testing procedure (e.g. Capelli and others, 2018).

Historically, quantifying snow microstructure was difficult and observations using magnifying glasses and crystal cards only characterized ‘grains’ representing the fragmented ice matrix and relied on the classification by an observer. Because density is relatively easy to measure, it has commonly been used to parameterize important physical quantities. However, despite its ease of measurement, it still leads to large uncertainties (Proksch and others, 2016). At the same time, the importance of microstructure has long been recognized (e.g. Good, 1987; Shapiro and others, 1997). Today, micro-computed tomography (µCT) imaging allows for a detailed analysis of snow microstructure (e.g. Coléou and others, 2001). Still there is a lack of a comprehensive understanding of the relationship between snow mechanics and snow microstructure. While Sundu and others (2024b) showed a first parameterization for the elastic modulus depending on snow density and microstructure, similar results on strength and fracture toughness are missing. Aside from the microstructural aspect, the experimental data on the strength of weak layers under multiaxial loading conditions are sparse, especially toward pure compression. Also it is still unclear, if artificially grown weak layers are a valid substitution for natural depth hoar in terms of their microstructure.

Therefore, our goal is to establish a link between the uniaxial compressive strength and microstructural information of depth hoar weak layers obtained with µCT imaging. To characterize a wide range of different microstructural morphologies, we used artificially produced specimens. We compared their microstructure with natural depth hoar to assess the suitability of artificial samples for conducting experiments. This study is considered as a proof of concept for subsequent experiments assessing the effect of microstructure on the strength under multiaxial loading conditions.

2. Methods

The workflow for sample preparation is visualized in Figure 1 and then described in detail below. First, we produced nature-identical snow or used natural snow from which we prepared a sandwich-like snow sample (sieved lower-density snow between denser slab layers). This sample was then exposed to a temperature gradient in a metamorphism box. From the resulting parent sample, five specimens were cut with a band saw for the experiments and one for the µCT scan. The workflow from snow production to actual testing took about 10 days for every parent sample.

Figure 1. Workflow from sample preparation to experimental testing and analysis of microstructure.

2.1. Sample preparation

Artificially grown depth hoar has been used in the past for studies of the mechanical behavior, thermal conductivity and settlement behavior (e.g. Fukuzawa and Narita, 1993; Reiweger and Schweizer, 2010; Staron and others, 2014; Capelli and others, 2018). To perform mechanical experiments, it is necessary to have the weak layer contained between dense and cohesive snow. This allowed us to handle the fragile weak layers and acts as an interface with the sample holders on the testing machine. We achieved this by sieving a layer of lower-density snow (20– $40\,\mathrm{mm}$ in thickness) onto a dense slab layer and covering it with a second slab. We made these 30– $50\,\mathrm{mm}$ thick slabs from nature-identical snow (Schleef and others, 2014) or natural snow consisting of partly decomposed precipitation particles (DF) or rounded grains (RG) according to Fierz and others (2009) which we compressed with 80 kPa between two $10\,\mathrm{mm}$ steel plates. To sieve the weak layer, we again used natural or nature-identical snow. Depending on the desired density of the weak layer, we varied the sieve mesh size (1.5– $4\,\mathrm{mm}$) and the crystal type of the snow (precipitation particles (PP), DF, RG). We found the compressed slabs needed to have a density of about 100 kg m⁻³ higher than the weak layer for the samples to fail in the weak layer. Overall, we produced 28 parent samples with weak layers of faceted crystals and depth hoar of different densities.

2.2. Snow metamorphism

To apply the temperature gradient to the prepared snow sandwich, we used custom build metamorphism boxes. The boxes consisted of $60\,\mathrm{mm}$ XPS foam panels within a wooden frame. Inside the insulation, a $10\,\mathrm{mm}$ thick aluminum plate with a 35 W heating foil attached allowed to control the bottom temperature of the snow sample using a calibrated DS18B20 temperature sensor (Maxim Integrated) and an Arduino-based PID controller. We used a $1.5\,\mathrm{mm}$ aluminum sheet on top of the sample to prevent sublimation and the formation of faceted crystals at the exposed surface due to fluctuations of the cold laboratory temperature (approximately $\pm1^{\circ}\mathrm{C}$). We applied gradients of 100–250 K m⁻¹ (over the whole sample, the relevant gradient in the weak layer is likely higher), which were determined by both the set temperature of the heating plate ( $-0.5^{\circ}\mathrm{C}$) and the cold lab temperature to which the top of the sample was exposed to ( $-10^{\circ}\mathrm{C}$ to $-20^{\circ}\mathrm{C}$). This allowed both the temperature gradient and the average temperature within the sample to be varied. To obtain different microstructural morphologies, we also varied the duration of metamorphism between 45 and 200 hours. From the resulting 28 parent samples ( $300\,\mathrm{mm}\,\times\,400\,\mathrm{mm}$), we cut out 5 specimens for the experimental tests ( $65\,\mathrm{mm} \times~140\,\mathrm{mm}$) and 1 specimen for the µCT scans (diameter $80\,\mathrm{mm}$) with a band saw. The edges of the parent sample (outer $50\,\mathrm{mm}$) were not used to avoid edge effects and the influence of the nonlinearity of the temperature gradient. The specimens were stored at $-5^{\circ}\mathrm{C}$ for at least 2 hours prior to the experiments to minimize temperature induced scatter in our results.

2.3. Experimental testing

For the mechanical experiments, we tested the compressive strength of five specimens per parent sample from which we calculated the average strength. The bottom of each specimen was frozen to a sample holder and the surface on the upper face was then leveled using a fine saw blade. Compression tests were performed using a uniaxial testing machine (Shimadzu Corporation) in combination with a 10 kN load cell with a sampling rate of 1 kHz (Fig. 2). The measurement accuracy is within ±1% of the indicated test force (at 1/500 to 1/1 load cell rating). The compression plate allows a rotation of $\pm 2.5^{\circ}$ around the vertical axis, thus allowing for a uniform load distribution, even if the sample surface is not perfectly aligned with the machine crosshead. Since snow properties are strain rate dependent, we aimed to have similar strain rates for all experiments. We set the displacement speed of the machine to 0.1– $0.4\,\mathrm{mm}~\mathrm{s}^{-1}$, which resulted in a theoretical strain rate (weak layer thickness in the range of 10– $40\,\mathrm{mm}$) of $10^{-2}\,\mathrm{s}^{-1}$, assuming that all the deformation is concentrated within the weak layer. The effective strain rate was estimated using digital image correlation (DIC) (see Section 2.4).

Figure 2. Setup for the mechanical tests showing the uniaxial testing machine, the machine crosshead with aluminum compression plate and the high speed camera. For reference, the compression plate is $140\,\mathrm{mm}$ wide.

For our analysis, we normalized the parameters with the properties of ice, which is commonly done for porous materials (Gibson and Ashby, 1997). To distinguish between the properties of the porous material and the bulk material, we introduce signifiers: ‘ ${\circ}$’ for the porous material (snow) and ‘ ${\bullet}$’ for the properties of the bulk material (ice). For strength of ice, we used the yield strength ${\sigma^{\bullet}_{\mathrm{y}}} \approx$ 2 MPa of polycrystalline ice (e.g. Haynes, 1978; Petrovic, 2003). We refer to polycrystalline ice because the properties of snow are governed by the bonds between crystals (Chandel and others, 2015). Snow grains can be seen as monocrystalline units (Riche and Schneebeli, 2013) with random lattice structures, so the bonds typically comprise a crystal lattice boundary. Since our experiments are in uniaxial compression, but the microscopic failure is likely to be due to bending moments in necks between grains (Ballard and McGaw, 1965), we use the yield strength of polycrystalline ice.

2.4. Video recordings and image correlation

We recorded each mechanical experiment with a Phantom Veo 710L high-speed camera (Vision Research - 1280 px × 720 px, which corresponds to approximately 9 px mm⁻¹) at 1 kHz to identify unsuccessful experiments (e.g. uneven loading, premature failures in upper or lower snow layers). It was relatively easy to identify failed experiments as can be seen in Figure 3. We also used a few of the recordings to perform DIC to estimate the effective strain within the weak layer. To this end, we used the tracking mode of the open source Digital Image Correlation Engine software (Turner and others, 2015). The surface of the snow samples was sprayed with black ink for better contrast. To evaluate the strain of the weak layer, we selected narrow regions of interest immediately below and above the weak layer and calculated the difference of the average displacement of the two regions. By evaluating the relative displacement between the two regions, we also eliminated noise due to micro-movements of the camera.

Figure 3. Examples of sample failure: (a) a sample that failed as intended in the weak layer and (b) a sample where a failure across the top layer occurred. The extent of the weak layer is highlighted in yellow. The samples are sprayed with black ink to improve the optical contrast.

This analysis showed that most of the displacement of the machine did not contribute to the strain within the weak layer (Fig. 4). Instead, a large proportion was absorbed at the sample–machine interface and by elastic and plastic deformation of the slab layers. We therefore analyzed six samples representing the extreme conditions of the effective strain rate (smallest vs largest weak layer thickness, lowest vs highest density, lowest vs highest displacement rate) using DIC. From this, we estimated the effective strain rate within the weak layer before failure to be in the range of $1 \times 10^{-3}$to $5 \times 10^{-3}\,\mathrm{s}^{-1}$.

Both the displacement and stress signals exhibit so-called pop-ins, which are abrupt changes observed in both the stress and strain signals. Figure 4 shows that for each pop-in the downward displacement of the machine crosshead increased, while all the other regions moved upward. This suggests that the pop-ins must originate from the interface between the machine crosshead and the specimen surface. We therefore assume these pop-ins are artifacts caused by the localized crushing of small irregularities at the sample–machine interface. Such events are not indicative of a failure within the weak layer and therefore had to be distinguished from a true failure of the weak layer, which was easily detectable in our experiments. Failure was indicated by a rapid decrease in stress after the peak and a relatively constant stress signal thereafter. In contrast, for pop-ins the stress increased again.

Figure 4. Exemplary displacement analysis using digital image correlation (DIC) of a sample during a compression experiment. (a) Shows the regions of interest tracked via DIC and (b) the corresponding average vertical displacements of the indicated regions during the compression test. The effective displacement within the weak layer is highlighted in yellow. (c) Shows the corresponding stress signal recorded by the testing machine.

2.5. Microstructure via $\mu$CT

Each artificial parent sample was characterized using µCT imaging (µCT 90, Scanco Medical), giving us a total of 28 scans. In order to avoid breaking the fragile weak layers while cutting them in preparation for the µCT scan, we had to use the largest sample holder (diameter 80 mm). This limited the effective resolution to a voxel size of (29.1 μm)³. For the scan, we used a source current of 145 μA and 55 kV. We chose an integration time of 300 ms and averaged the data twice. The reconstructed geometry was Gaussian filtered to reduce noise (width = 1.2 voxels, support = 2 voxels) and binary segmented (Hagenmuller and others, 2013). The µCT scans were used to (i) calculate average (bulk) properties of the weak layer (density, specific surface area (SSA), correlation lengths, anisotropy, connectivity density) and to (ii) perform a moving window analysis to obtain high resolution (about $3\,\mathrm{mm}$) vertical profiles of all samples (density, SSA).

For the bulk analysis (i), a cubic region of interest (ROI) containing only the weak layer was manually selected from the scan. This region had to be cubic for subsequent correlation length analysis and was, therefore, limited by either the height of the weak layer or the size of the sample holder. Typically, the ROI was larger than 1000 voxels ((29.1 mm)³). The minimum ROI size was 278 voxels ((8.09 mm)³), which we assume to be a sufficiently large representative volume in good agreement with the findings of Kaempfer and others (2005) and Coléou and others (2001). The average weak layer density $\rho^{\circ}_\mathrm{bulk}$ was calculated from the volume fraction of the segmented binary µCT image after triangulation (Proksch and others, 2016; Calonne and others, 2020b) and the density of ice ( $\rho^{\bullet} = 916.7\,\mathrm{m}^2\,\mathrm{kg}^{-1}$—Harvey (2017)). SSA $_\mathrm{bulk}$ was calculated from the triangulated surface of the binary image which was normalized by the ice volume (Schleef and Löwe, 2013). We calculated the weak layer vertical anisotropy α from the µCT exponential correlation lengths $p_\mathrm{ex}$ according to Calonne and others (2014):

(1)\begin{align} \alpha = \frac{p_\mathrm{ex,z}}{0.5 \left(p_\mathrm{ex,x} + p_\mathrm{ex,y}\right)}. \end{align}

Values greater than 1 for α indicate vertical structures and values smaller than 1 indicate horizontal structures. We also calculated the connectivity density (ConnD) of our CT scans. This parameter quantifies the degree of microstructural connectivity by normalizing a topological measure of the structure (based on the Euler characteristic) to the unit volume, as described by Odgaard and Gundersen (1993).

For the moving window analysis (ii), we chose a window size of 1000 × 1000 × 100 voxels (x, y, z) and a vertical overlap of 50% (Fig. 5a). Thus, the theoretical vertical resolution was 1.46 mm with a window height of $2.91\,\mathrm{mm}$. Density and SSA were calculated for each slice as described for the bulk properties above. The resulting profiles (an example is shown in Fig. 5b) generally showed that the properties were not uniform with height. Comparison of the profiles with the bulk measurements (solid vertical lines) showed that the bulk density was typically substantially higher than the minimum density. As strength is expected to be sensitive on snow density, we also used the minimum density within the weak layer $\rho^{\circ}_\mathrm{min}$ for our subsequent statistical analysis.

Figure 5. Comparison of the µCT bulk evaluation and moving window analysis. (a) Shows a vertical cross section through the µCT scan of an artificial weak layer sample and the size and overlap of the moving window. The dashed lines are the manually selected extent of the weak layer. (b) Shows a representative example comparison of the moving window analysis and the bulk measurements (representing weak layer averages, shown as vertical bars) for density (black squares) and specific surface area (SSA) (purple circles). For the moving window analysis, we used a window size of 1000 × 1000 voxels, a window height of 100 voxels and a vertical overlap of 50%. Therefore, each value represents the average of a $2.91\,\mathrm{mm}$ thick slice.

Geometries reconstructed from µCT measurements are subject to measurement uncertainties and limitations due to discretization. The measurement uncertainties arise from the measurement process itself. Then, the discretization of the image into voxels with a defined gray scale range affects the representativeness of the reconstruction results further. For snow, this is the main concern, as resolution limitations could lead to underrepresentation of fine structures in the binarized image and thus affect the calculated density and SSA. Kerbrat and others (2008) estimated the error of µCT-based SSA measurements to be 3% for SSA values in the range of 5– $64\,\mathrm{mm}^{-1}$ using a voxel size of 10–18 μm based on a comparison with gas adsorption values. Furthermore, Walter and others (2024) found that their CT-based SSA values were too low due to resolution limitations for SSA values larger than approximately $30\,\mathrm{mm}^{-1}$. For snow density and anisotropy, to the best of our knowledge, no objective estimate of the error has been reported in the literature. Freitag and others (2004) found that µCT density measurements of firn were in good agreement (deviations <1%) with results obtained using the gamma absorption method using a voxel size of 40 μm. Furthermore, Hagenmuller and others (2013) reported density deviations in the order of 5% resulting from the segmentation process.

The measurement uncertainties of the µCT mentioned above do not capture the variability within the parent sample and are, therefore, not sufficient to explain the scatter in our results. To obtain an estimate of the variability within a parent sample, we produced two additional samples, one with a low-density (LD) weak layer and one with a high-density (HD) weak layer. We took five samples distributed across each parent sample and analyzed them. The results are summarized in Table 1.

Table 1. Microstructural variability due to spatial variations for a low-density (LD) and a high-density (HD) weak layer. Values are calculated from five individual µCT scans distributed over the parent sample. Bold values indicate those used for error estimation

The primary source of uncertainty in density measurements was the spatial variability within a parent sample. As this variability strongly decreased at higher densities, we adopted the larger value of 11.9 kg m⁻³ as a fixed uncertainty instead of a percentage-based estimate. Similarly, for SSA and connectivity density, we observed a reduction in variability at higher densities and, therefore, used fixed uncertainty values of $0.635\,\mathrm{mm}^{-1}$ and $1.35\,\mathrm{mm}^{-3}$, respectively. Anisotropy remained similar across densities, leading us to use 0.78% as the estimated error.

These values are only applicable to weak layers. To also indicate the uncertainty in our time-series measurements (Fig. 6), we used the microstructure time-series data (see Sections 2.6 and 3.1) to calculate the average difference between the two samples for each layer (Table A1, Appendix). We assume that this difference between two identically produced samples is an upper limit for the spatial variability.

2.6. Microstructure time series

To better understand the depth hoar growth in our artificial samples and to quantify the repeatability of our setup, we performed daily µCT measurements on two identically produced parent samples that were placed in the metamorphism box for 1 week. To do this, we cut $100\,\mathrm{mm} \times 100\,\mathrm{mm}$ samples from the parent samples for scanning. The remaining holes were filled by sieving snow consisting of partly decomposed precipitation particles. Since the sieved patches probably influenced the global temperature gradient, the processes in an undisturbed metamorphism box might be slightly different. We used the daily µCT scans to show the evolution of the microstructural properties density, SSA and anisotropy over a week and to estimate the repeatability of our setup.

2.7. Comparison of microstructure datasets

To demonstrate the comparability of our artificially grown microstructural morphologies with natural depth hoar, we re-analyzed µCT scans of natural depth hoar from two large datasets. First, during the RHOSSA campaign (Calonne and others, 2020b), a dataset of µCT scans, snow micro pen (SMP) profiles and traditional profiles was collected over one winter season at the Weissfluhjoch, Davos. This dataset represents a typical alpine snowpack. The maximum temperature gradient within the snowpack was up to 50 K m⁻¹ and was reached in the early winter. We analyzed all µCT scans containing depth hoar and faceted crystals (about 30 individual scans—Calonne and others 2020a). Second, during the MOSAiC expedition (Macfarlane and others, 2023), an even larger dataset of microstructural properties of snow on sea ice was collected. It consists of SMP, µCT and density measurements as well as traditional snow profiles. This dataset represents rather extreme conditions for depth hoar formation (low snow heights and strong gradients, typically 50–150 K m⁻¹). From this dataset, we also analyzed about 30 scans from the period December 2019 to April 2020 (Macfarlane and others, 2022).

To create the artificial weak layers, we used higher temperature gradients than those typically found in alpine snowpacks to speed up the process and increase the number of samples. To asses the influence of the high temperature gradient, we also created an artificial sample with more realistic environmental conditions (temperature gradient over the whole sample: 30 K m⁻¹, average temperature in the weak layer $-13^{\circ}\mathrm{C}$ and a metamorphism duration of 44 days). This allowed us to assess the effect of the high temperature gradients used and to demonstrate the microstructural capabilities of our setup.

3. Results

3.1. Temporal evolution within the metamorphism box

The evolution of the microstructural properties density, SSA, anisotropy and connectivity density over the course of 1 week is shown in Figure 6 for two identical samples. Density increased only slightly in both the slabs and the weak layer. SSA decreased throughout the sample. For the vertical anisotropy, we observed the strongest increase within the first few days, after that a plateau seemed to be reached. The connectivity density decreased over time, especially in the slab layers. The values for the weak layer are one order of magnitude lower than for the slabs.

Figure 6. Temporal evolution of density, specific surface area and anisotropy in the metamorphism box. (a) Shows the different regions in a µCT scan, additionally the temperature gradient within the sample is qualitatively visualized. The graphs show time series of (b) density ( ${\rho_{\mathrm{bulk}}}$), (c) specific surface area ( $\mathrm{SSA_{bulk}}$), (d) anisotropy (α) and (e) connectivity density (ConnD). Measurements from the two different experiments are shown with different markers, the solid line represents the average of the two measurements. The error bars show the estimated error due to spatial variations.

3.2. Similarity to natural depth hoar

We compared the microstructure of our artificial samples with natural depth hoar data from the RHOSSA and MOSAiC campaigns. The SSA of our samples showed a clear decreasing trend with increasing density and the values were similar to natural samples above relative densities of ∼0.25 (Fig. 7a). The relation between bulk SSA and density for the artificial samples can be described with a linear regression (Eqn (2)). For the vertical anisotropy (Fig. 7b), we observed a strong increase for higher densities (Eqn (3)). The vertical anisotropy of the RHOSSA samples was consistently lower than that of our artificial samples. Half of the MOSAiC samples were comparably anisotropic, while the other half had lower vertical anisotropies. The connectivity density of our samples was comparable to the MOSAiC samples but systematically higher than most of the RHOSSA samples (Fig. 7c and Eqn (4)). The microstructure of the long-term artificial sample, which was grown under more realistic conditions, had a lower SSA and connectivity density than the other artificial samples, but comparable anisotropy.

(2)\begin{align} \mathrm{SSA_{bulk}} = -16\left( \frac{{\rho^{\circ}_{\mathrm{bulk}}}}{{\rho^{\bullet}}} \right) +14, \quad \mathrm{\textit{R}}^2 = 0.46 \end{align}

(3)\begin{align} \alpha = 1.6\left( \frac{{\rho^{\circ}_{\mathrm{bulk}}}}{{\rho^{\bullet}}} \right) +0.84, \quad \mathrm{\textit{R}}^2 = 0.71 \end{align}

(4)\begin{align} \mathrm{ConnD} = 19\left( \frac{{\rho^{\circ}_{\mathrm{bulk}}}}{{\rho^{\bullet}}} \right) -1.1, \quad \mathrm{\textit{R}}^2 = 0.75. \end{align}

Figure 7. Comparison of the microstructure of the artificial samples with natural depth hoar as a function of normalized density. (a) Shows $\mathrm{SSA}_{\mathrm{bulk}}$ of our artificial weak layers and the corresponding linear regression in comparison with samples from the RHOSSA and MOSAiC campaigns. (b and c) Show the same for anisotropy α and connectivity density (ConnD), respectively. All plots include data from one long-term artificial sample which has been grown under more realistic conditions, visualized with a cross.

3.3. Compressive strength

A total of 92 specimens from 24 parent samples were analyzed. Our results suggest a strong nonlinear relation of compressive strength and depth hoar density (Fig. 8). We used an orthogonal distance regression (ODR—Boggs and Rogers (1990)) to fit a power law using the bulk weak layer density:

(5)\begin{align} \frac{\sigma^{\circ}_\mathrm{c}}{\sigma^{\bullet}_\mathrm{y}} = 11 (\pm 7) \left(\frac{\rho^{\circ}_\mathrm{bulk}}{\rho^{\bullet}}\right)^{5.4 (\pm 0.5)}, \quad \chi^2_{\nu} = 2.9. \end{align}

Figure 8. Mean compressive strength per parent sample vs normalized density. The line represents the power law fit obtained via orthogonal distance regression (ODR). For comparison, we have included data from the literature on the compressive strength of faceted crystals and depth hoar. The dashed line shows the power law valid for the crushing strength of most cellular solids (Gibson and Ashby, 1997). The vertical error bars indicate the standard deviation of the individual compression tests for one artificial parent sample (five specimens).

The reduced chi-square parameter $\chi^2_{\nu}$ provides a measure of how well a model fits a set of data, normalized by the degrees of freedom. Values close to 1 indicate a good fit while values significantly higher suggest uncertainties in the data or inadequacies in the model (Andrae and others, 2010). For reference, we have also plotted the power law of Gibson and Ashby 1997 for cellular solids (dashed line Fig. 8). At high weak layer densities (bulk density >350 kg m⁻³), it was not possible to maintain the density difference of >100 kg m⁻³ between the weak layer and the slab layers, so many samples failed at the upper slab layer first. We, therefore, excluded these results from the power law fit.

3.4. Influence of microstructure

Various microstructural parameters were obtained from the µCT scans. We analyzed the statistical significance of each parameter with the experimentally determined mean compressive strength $\sigma_{\mathrm{c}}$ using a Spearman rank correlation. We chose this method because it is able to capture nonlinear correlations in a representative way and because it is less sensitive to outliers than other methods. Since the total number of observations was relatively small for this purpose, we performed a permutation test to obtain meaningful p-values (Phipson and Smyth, 2010). Excluding the samples which failed at the slab layers, we found that all SSA and density parameters as well as connectivity density correlated significantly with $\sigma_{\mathrm{c}}$ (Fig. 9). The anisotropy was statistically highly significant, although the individual exponential correlation lengths were not. The complete Spearman rank correlation matrix is shown in Figure A1, Appendix. We repeated the same analysis with the residuals obtained from the ODR power law fit (Eqn (5)) of compressive strength and bulk density to exclude the influence of density on other parameters. For the residuals, we found that only the horizontal correlation lengths $p_\mathrm{ex,x/y}$ correlated with the compressive strength (p = 0.04 and p = 0.03, respectively).

Figure 9. Results of the Spearman correlation. The orange p-values are calculated based on the correlation of weak layer parameters with compressive strength $\sigma_{\mathrm{c}}$ and the purple p-values are based on the correlation of weak layer parameters with residuals of the power law fit (Equation 5). Additionally, the Spearman correlation coefficients $r_\mathrm{s}$ are shown. The vertical dashed lines indicate levels of significance α.

4. Discussion

4.1. Temporal evolution within the metamorphism box

The µCT time series (Fig. 6) of two independent but identically setup metamorphism boxes allowed us to estimate the repeatability of our method and revealed the changes in microstructure due to the temperature gradient metamorphism. Both density and SSA of the two samples were very similar in the weak layer. However, anisotropy and connectivity density showed differences of up to 50%, which may be attributed to the spacial variability within the weak layer and to the perturbation of the linear temperature gradient due to the sample extraction for the µCT scan. The connectivity density of the weak layer was an order of magnitude lower than that of the slab layers. This is probably due to the unit volume normalization, which adds a strong density scaling to this parameter (see Fig. 7). In terms of microstructure, the time series shows a similar trend across all layers. Therefore, the main microstructural difference between the different regions of the sample was the density, highlighting the need to have higher densities in the slab layers compared to the weak layer to ensure that the fracture happens only in the weak layer during the mechanical experiments.

4.2. Similarity to natural depth hoar

Our setup allowed us to perform compressive strength experiments for a wide range of depth hoar densities (150–350 kg m⁻³), which would be difficult to achieve with natural snow. We used high temperature gradients (100–250 K m⁻¹ over the whole sample) to speed up the process. Gradients of this magnitude can occur during surface faceting due to diurnal temperature changes (Birkeland and others, 1998) and in sub-arctic (Sturm and Benson, 1997) and arctic (Macfarlane and others, 2023) snowpacks but are otherwise not commonly found in alpine snowpacks (e.g. Calonne and others, 2020b).

The analysis of the microstructure of our artificial samples (Fig. 7) showed a decreasing trend for SSA with increasing densities. This agrees well with findings using the gas adsorption method, for instance, the parameterization by Domine and others (2007) and the overall trend shown by Legagneux and others (2002) and µCT-based SSA evaluations (e.g. Schneebeli and Sokratov, 2004). For the vertical anisotropy, we found an increase with increasing snow density. This dependence is somewhat unexpected but is in line with the measurements of Calonne and others (2019) for faceted crystals and depth hoar and our analysis of the RHOSSA µCT scans (Fig. 7). Connectivity density also showed an increasing trend with density. This effect may be partly attributed to the closer proximity of grains at higher densities, as well as to the normalization by unit volume, which introduces a density-scaling behavior.

The comparison of the microstructural quantities SSA, anisotropy and connectivity density with natural samples (Fig. 7) suggests that our artificial depth hoar samples were comparable to natural ones. We found that the SSA of our artificial samples was systematically slightly higher than that of natural samples. As our more realistic long-term sample showed lower SSA values, we suggest that this offset is due to the short time of metamorphism of our samples. The anisotropy was similar to the majority of MOSAiC samples, while RHOSSA samples had systematically lower anisotropy. This could be due to the different temperature gradient, but also to the different weight of the overlying snowpack. The latter is supported by our long-term sample, which had a lower density and a greater vertical anisotropy than the RHOSSA campaign samples, although the SSA is comparable. The connectivity density of our artificial samples closely matched that of the MOSAiC samples, while most RHOSSA samples consistently showed lower values. In addition, our long-term sample showed reduced connectivity density, suggesting that these differences may be due to variations in metamorphism driven by different temperature gradients.

4.3. Compressive strength

Failure of weak layers in nature occurs if the multiaxial stress state exceeds the strength of the weak layer. In this study, we focused on the influence of snow microstructure on the uniaxial compressive strength and assessed only one dimension of the failure envelope.

Our strength values (Fig. 8) represented the mean of five mechanical tests per parent sample. This allowed us to obtain robust experimental results, although considerable standard deviations (represented by the error bars) of compressive strength (around 20%) were observed across the density range. This value is comparable to scatter reported in other experimental studies with snow (e.g. Jamieson and Johnston, 2001; Chandel and others, 2014; Reiweger and others, 2015) and indicates that the repeatability of our sample growth process was similar independent of the weak layer density. We attribute the observed scatter to variations in the microstructure of the parent sample (e.g. local differences in density, different notch effects, different cross sections of the load-bearing structures) and the test method (e.g. uneven load distribution due to a non-parallel sample interface and different contact points). It is also likely that the sample preparation process damaged some of the weak layer crystal structures at the edges. This would result in a systematic offset of the strength values, as the effective load bearing area may have been smaller than the sample dimensions. If we assume that only the very outer crystals were damaged (grain size on the order of 2– $3\,\mathrm{mm}$), this would correspond to an underestimation of the strength by less than 10%, which is a relatively small amount compared to the overall scatter in our data.

Although there are only few data available for comparison, our strength data agree well with most previously reported values for faceted crystals and depth hoar (see literature data in Fig. 8). Reiweger and Schweizer (2013) performed load-controlled experiments and reported compressive strengths of 2.5–6 kPa for different loading rates for a depth hoar weak layer with a density of 190 kg m⁻³. This is close to our data, though not directly comparable with our displacement-controlled test results in terms of loading rate. Chandel and others (2014) reported values between 0.3 kPa and 1.4 kPa for faceted crystals with densities between 100 and 200 kg m⁻³, which is also similar to our results. Strength data obtained using numerical modeling also shows comparable values. Mede and others (2018) found the compressive strength to be approximately 2 kPa for a sample consisting of faceted crystals and depth hoar with a density of 180 kg m⁻³ using discrete element modeling. The finite element results of Chandel and others (2015) showed lower compressive strength: approximately 6 kPa for faceted crystals with a density of 250 kg m⁻³ and 7.5 kPa for a density of 416 kg m⁻³, which suggests a very low dependence on density, not in line with our findings. Other numerical studies of multiaxial strength again show values comparable to our results (e.g. Mulak and Gaume, 2019; Bobillier and others, 2020; Ritter and others, 2020; Singh and others, 2022). However, most of these numerical studies used some sort of simplified snow microstructure and results cannot be compared directly.

As our data have uncertainties in both strength and density (due to inhomogeneities within the parent sample), we chose to use an ODR to fit the power law, which can account for measurement uncertainties in the fitting process. The $\chi^2_{\nu}$ value allows the goodness of the fit to be assessed, taking into consideration the degrees of freedom in the fitting process. The value obtained from our fit was $\chi^2_{\nu} = 2.9$, which indicates that our power laws fit the data well, but the errors in strength and density did not fully explain the observed scatter.

Mechanical experiments on weak layers are challenging even in the controlled environment of the cold lab. The samples are fragile and sample preparation (cutting, freezing to the bottom sample holder) can damage the weak layer. In addition, the top and bottom surfaces of a sample are typically not parallel, requiring manual post-processing upon sample preparation. It should be noted that even after post-processing, the load distribution during the compression test may not have been perfect due to uneven surface of the sample interface. Samples where tilting of the machine crosshead prior to failure was detected on the high speed videos were excluded from the analysis. Nevertheless, small tilting may have resulted in unwanted moments resulting in shear forces in the weak layer. It should also be noted that the effective strain rate in the weak layer is very different from the theoretical rate resulting from the machine speed. In our experiments, we found the strain rate in the weak layer to be on the order of $10^{-3}\,\mathrm{s}^{-1}$. This is close to the ductile-to-brittle transition (e.g. Löwe and others, 2020). The elastic behavior of the pop-ins in the displacement signal obtained by DIC (Fig. 4) suggests quasi-brittle material behavior.

4.4. Influence of microstructure

The large sample dimensions we used did not allow testing and scanning of the same sample. As a result, the strength and microstructure cannot be directly related using, for example, the finite element method. Instead, we extracted five specimens for experimental testing and one specimen for the µCT scan from one parent sample. We assume that the microstructure of the µCT scan is representative of the microstructure of the experimental test specimens, but small spatial variations are to be expected and have been discussed in Sections 2.5 and 4.1.

We found that minimum density, minimum SSA and maximum SSA from the µCT profiles added significant microstructural information at a significance level of 0.01, suggesting that the µCT profiles represented a valuable addition to the standard bulk analysis. Interestingly, while the correlation with anisotropy was statistically significant, the individual exponential correlation lengths were not. This suggests that shape and orientation of the microstructure influence the mechanical properties, which agree well with the findings of Walters and Adams (2014), Sundu and others (2024a) and Walters and Adams (2012). The thickness of the weak layer $h_{\mathrm{wl}}$ had no significant influence on the compressive strength. This seems plausible if we assume quasi-static stress distribution within the weak layer which would lead to failure at the weakest link, likely the region with the lowest density (Hagenmuller and others, 2014a).

As SSA, connectivity density and anisotropy were also correlated with density (see Fig. 7 and Figure A1, Appendix), we repeated the same analysis with the residuals obtained from the ODR power law fit of compressive strength and bulk density (Fig. 8 and Eqn (5)), thus excluding the direct influence of density. We found that only the horizontal correlation lengths $p_\mathrm{ex,x/y}$ correlated with the residuals of the power law fit (p = 0.04 and p = 0.03, respectively). This suggests that the extent of horizontal structures within the weak layer had some influence on the strength. The fact that anisotropy is not significant here could indicate that its effects are already accounted for through the correlation of density and anisotropy (see Eqn (3)), or that the driving mechanism is not related to orientation within the material. Instead, it may be that a specific horizontal length scale of the material plays an important role in the microscopic failure of the material. Since correlation lengths symmetrically capture the length scales of both phases (solid and void), their interpretation in this context remains ambiguous.

Overall, our results highlight that density was the primary factor influencing strength, with most of the strong correlations of other quantities with strength appearing to be artifacts of their dependence on density. However, an analysis of the residuals revealed that horizontal dimensions also play a role. This suggests that additional dependencies may exist but remain obscured by scatter in the data and correlations between parameters.

4.5. Power law fit and the relation to foams

Due to the similarity in terms of porosity, the principles of foamed engineering materials have often been used in the past to describe the mechanical behavior of snow (e.g. Schneebeli and others, 1999; Jamieson and Johnston, 2001; Blatny and others, 2023). Gibson and Ashby 1997 showed that the strength of many foam-like materials follows a power law as a function of density with an exponent of 1.5. Our results also show a power law relation of strength with density, but with a much larger exponent of 5.4 (Fig. 8). This discrepancy suggests that the foam model is not suitable for describing the material behavior of depth hoar under compression. Our exponents are also larger than those reported by Jamieson and Johnston (2001) for the shear strength of weak layers. This suggests that the compressive strength is more sensitive to variations in density than the shear strength, which may be due to the microstructure of snow.

Although snow is a bi-continuous material and has been considered a ‘foam of ice’ (Kirchner and others, 2001), its microstructure differs from other cellular solids in that it is composed of grains that form bonds with their immediate neighbors. These bonds are typically much smaller than the grains and therefore strongly influence the mechanical response of the material (Agrawal and Mittal, 1995). The inconsistency with other foamed materials becomes more apparent when the influence of density on the microstructure is assessed. Figure 10 compares the microstructure of an low density (182 kg m⁻³) and an high density (355 kg m⁻³) weak layer of depth hoar from our dataset. It can be see that the individual crystals have roughly similar shapes and sizes, the most obvious difference is the smaller distance between grains in the dense sample.

Figure 10. 3-D reconstruction (top) and 2-D section (bottom) of the microstructure of a low bulk density sample (left) and high density sample (right). The microstructural parameters of the left sample are: $\rho^{\circ}_\mathrm{bulk} = 182\,\mathrm{kg}\,\mathrm{m}^{-3}$, $\mathrm{SSA}_\mathrm{bulk} = 10.8\,\mathrm{mm}^{-1}$, $\mathrm{ConnD} = 1.51\,\mathrm{mm}^{-3}$, $p_\mathrm{ex,x} = 0.192\,\mathrm{mm}$, $p_\mathrm{ex,y} = 0.187\,\mathrm{mm}$, $p_\mathrm{ex,z} = 0.224\,\mathrm{mm}$, α = 1.18. The microstructural parameters of the right sample are: $\rho^{\circ}_\mathrm{bulk} = 355\,\mathrm{kg}\,\mathrm{m}^{-3}$, $\mathrm{SSA}_\mathrm{bulk} = 9.00\,\mathrm{mm}^{-1}$, $\mathrm{ConnD} = 5.93\,\mathrm{mm}^{-3}$, $p_\mathrm{ex,x} = 0.160\,\mathrm{mm}$, $p_\mathrm{ex,y} = 0.162\,\mathrm{mm}$, $p_\mathrm{ex,z} = 0.235\,\mathrm{mm}$, α = 1.46. Size: 300 voxels corresponding to $8.73\,\mathrm{mm}$.

This distance between grains defines an effective length scale $l_\mathrm{eff}$ for failure mechanism considerations. The hypothesis that the distance between grains influences compressive strength agrees well with the small, but statistically significant negative correlation coefficients we found for the horizontal correlation lengths $p_\mathrm{ex,x/y}$ from the ODR residuals (Fig. 9). At the same time, smaller distances between grains makes it more likely to have more bonds per grain. This density dependence of weak layer topology is to some extent reflected in the parameter connectivity density, which is a measure for the connected structures per unit volume (see Figs 6 and 10). The number and size of bonds add another scaling factor to the density-strength relation. An additional complication to the foam analogy is that the snow’s ice mass does not contribute equally to the strength, with a significant portion of the stress being carried by a relatively small fraction of the mass (Schneebeli and Sokratov, 2004; Hagenmuller and others, 2014b; Wiese and Schneebeli, 2017).

The foam model by Gibson and Ashby 1997 accounts for several failure mechanisms (e.g. bending and buckling), which dominate the response of cellular materials. However, we assume that the mechanisms controlling the macroscopic strength of the depth hoar layers we tested may differ. If this is the case, we expect to find different exponents for other load configurations and different microstructures, which may provide a new approach to microstructure-based weak layer strength assessment for practical applications such as snow instability modeling.

4.6. Conclusions

We performed 92 uniaxial compression tests on layered snow samples containing an artificially grown weak layer of depth hoar. Our results show that our method of artificially producing specimens with a weak layer of faceted crystals and depth hoar is well suited to produce nature-like specimens, suited for extensive experimental campaigns. These artificial specimens allowed us to determine the compressive strength of a wide range of microstructures and densities in a controlled environment. The compressive strength ranged from 1 to 150 kPa for weak layer densities between 150 and 350 kg m⁻³. The strong dependence on snow density was best described by a power law with an exponent of 5.4. In addition, several microstructural metrics such as SSA, connectivity density and correlation lengths obtained from µCT measurements were statistically significantly related to compressive strength. Since many of the metrics also depend on density, we repeated the analysis based on the residuals obtained from the power law fit of compressive strength and bulk density, thus excluding the direct influence of density. Here we found that the horizontal correlation lengths $p_\mathrm{ex,x/y}$ showed a weak negative, still significant correlation (p = 0.04 and p = 0.03, respectively) with the residuals of the power law fit. In conclusion, while compressive strength primarily depends on snow density, our findings suggest that the horizontal correlation lengths, a proxy for horizontal dimensions within the microstructure, may also play a role in determining the material’s strength, potentially influencing load transfer and failure mechanisms at a micro scale.

In the future, we will test other types of weak layers and evaluate the influence of microstructure on strength under multiaxial loading conditions in a laboratory setup. Field experiments using a portable, displacement-controlled testing device could also provide valuable insights, particularly into the temporal evolution of strength in a natural snowpack.