# The influence of snow microstructure on the compressive mechanical properties of weak snowpack layers

Jakob Schöttner, Berit Zeller-Plumhoff, Pascal Hagenmuller, Philipp

Weißgraeber, Philipp L Rosendahl, Henning Löwe, Jürg Schweizer, Alec van Herwijnen

► To cite this version:

Jakob Schöttner, Berit Zeller-Plumhoff, Pascal Hagenmuller, Philipp Weißgraeber, Philipp L Rosendahl, et al. The influence of snow microstructure on the compressive mechanical properties of weak snowpack layers. Acta Materialia, 2025, 302, pp.121657. 10.1016/j.actamat.2025.121657. hal-05451252

HAL Id: hal-05451252

https://hal.science/hal-05451252v1

Submitted on 9 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Full length article

# The influence of snow microstructure on the compressive mechanical properties of weak snowpack layers

Jakob Schöttner<sup>a,∗</sup>, Berit Zeller-Plumhoff<sup>b,c</sup>, Pascal Hagenmuller<sup>d</sup>,
Philipp Weißgraeber<sup>e</sup>, Philipp L. Rosendahl<sup>f</sup>, Henning Löwe<sup>a</sup>, Jürg Schweizer<sup>a</sup>,
Alec van Herwijnen<sup>a</sup>

<figure><img src="image_3.png" /></figure>

Check for
updates

## ARTICLE INFO

Dataset link: https://www.doi.org/10.16904/e
nvidat.648

Keywords:

Snow avalanches

Snow mechanics

Snow microstructure

Snow strength

Morphology

Scaling behavior

Stiffness

Strength

Machine learning

## ABSTRACT

To find a link between microstructure and the mechanical properties of weak snow layers, we acquired a comprehensive dataset combining mechanical tests for compressive strength and stiffness on 278 weak layer samples with high-resolution 3D data from 70 co-located micro-computed tomography (µCT) scans. The samples encompass three main categories, often associated with dry-snow slab avalanche release: faceted crystals and depth hoar, decomposing particles and rounded grains, and buried surface hoar. The results for stiffness and compressive strength highlight different scaling behaviors with densities for different grain types, but also fundamentally different relations between stiffness and strength for persistent and non-persistent grain types. Using a neural network, we identified key clusters of microstructural features that are most predictive of mechanical behavior. These clusters were determined based on the intercorrelation of features and reflect different geometric and topological aspects of snow microstructure. Stiffness was primarily influenced by features describing the alignment and efficiency of load-bearing paths, such as anisotropy and tortuosity. In contrast, compressive strength was more sensitive to the local geometry of the ice-air interface, though it also depended on the load path characteristics. By distinguishing between these controlling factors, the model predicted the mechanical properties with an accuracy comparable to the experimental variance. While our results confirm that the mechanical properties primarily depend on density, they also demonstrate that snow microstructure shapes the specific relation between density and mechanical properties. This relationship differs markedly between persistent weak layers and those composed of decomposing particles and rounded grains.

## 1. Introduction

The mechanical properties of snow have been studied since the 1940s, primarily to understand and predict snow avalanches. In this context, layers of low strength, called weak layers, are critical, as they are a key determinant for the formation of dry-snow slab avalanches [1]. These layers typically consist of large grains and exhibit a highly porous structure [2–5]. Weak layers are generally categorized as persistent or non-persistent, based on their microstructural evolution. Persistent types — such as faceted crystals (FC), depth hoar (DH; Fig. 1a), and buried surface hoar (SH; Fig. 1c) — can remain mechanically weak over extended periods due to slow rates of strengthening through processes such as sintering and densification [6,7]. Non-persistent layers,

including precipitation particles (PP), decomposing fragments (DF), or rounded grains (RG; Fig. 1b), tend to stabilize more rapidly. Improved avalanche prediction thus requires robust mechanical characterization of snow, and particularly of weak layers [e.g., 8,9].

The mechanical behavior of snow strongly depends on temperature and loading rate, resulting in completely different failure behaviors, ranging from ductile to brittle [e.g., 10,11]. Because of its high homogeneous temperature, snow's microstructure is constantly evolving; new bonds form rapidly when particles touch, and broken bonds can reform through sintering [e.g., 12]. Like for other porous materials, the density of snow serves as a good first-order predictor for physical properties. Indeed, not only does its microstructure show similarities, but in terms

<footnote><sup>a</sup> WSL Institute for Snow and Avalanche Research SLF, Flüielastrasse 11, Davos Dorf, 7260, Switzerland</footnote>

<footnote><sup>b</sup> Chair of Data-Driven Analysis and Design of Materials, University of Rostock, Albert-Einstein-Str. 2, Rostock, 18059, Germany</footnote>

<footnote><sup>c</sup> Institute of Metallic Biomaterials, Helmholtz-Center Hereon, Max-Planck-Str. 1, Geesthacht, 21502, Germany</footnote>

<footnote><sup>d</sup> Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, 1441 rue de la Piscine, Grenoble, 38400, France</footnote>

<footnote><sup>e</sup> Chair of Lightweight Design, University of Rostock, Albert-Einstein-Str. 2, Rostock, 18059, Germany</footnote>

<footnote><sup>f</sup> Institute for Structural Mechanics and Design, TU Darmstadt, Franziska-Braun-Str. 3, Darmstadt, 64287, Germany</footnote>

<footnote>* Corresponding author.
E-mail address: jakob.schoettner@slf.ch (J. Schöttner).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Exemplary microstructures obtained from µCT images for each tested category of grain types: (a) faceted crystals and depth hoar, (b) decomposing particles and rounded grains, and (c) buried surface hoar.</figcaption></figure>

of (compressive) strength and Young's modulus, snow can be found among most other porous materials and cellular solids (see Fig. 2).

Density-based parameterizations have historically struggled to capture the extensive variability in snow's mechanical properties [e.g., 13–16], a challenge compounded by the logistical difficulty of accurately sampling density in thin or heterogeneous snow layers. Compilations of snow mechanical data show that for similar densities, there is a large scatter in measured properties. This variability is attributed to a combination of factors, including intrinsic differences in microstructure as well as extrinsic differences in experimental conditions, such as temperature, strain rate, and sample preparation. This discrepancy is especially pronounced for weak layers [17], as most research focuses on rounded grains and thus has limited applicability to the distinct crystal morphologies that constitute these layers [e.g., 16,18].

While some recent work incorporated structural anisotropy into stiffness parameterizations [19], most parameterizations still rely on density only, and comprehensive datasets remain sparse [16,20–22].

## 1.1. Theoretical framework

Since snow is effectively a bicontinuous material, research in its microstructural analysis has been closely linked to developments in related fields — such as foams, bone, aerogels, and nanoporous gold — where techniques such as micro-computed tomography (µCT) scanning have paved the way for significant advancements.

A foundational approach to theoretically explain the mechanical behavior of porous solids is the Gibson–Ashby theory [34], which relates mechanical properties to volume fraction via idealized beam-based unit cells. For example, the compressive strength $\sigma^\circ$ of the porous material scales with that of the bulk phase ($\sigma^*$) as:

$$ \frac{\sigma^\circ}{\sigma^*} = C_1 \left( \frac{\rho^\circ}{\rho^*} \right)^{1.5} = C_1 \phi^{1.5} \quad (1) $$

where $\rho^\circ$ and $\rho^*$ are the densities of the porous and bulk material, $C_1$ is a dimensionless constant, and $\phi$ is the volume fraction. The scaling exponent in this simplified model reflects the dominant deformation mode: stretching-dominated structures yield nearly linear scaling between strength and stiffness (exponent of 1.0 in the Ashby plot - Fig. 2), while bending-dominated structures exhibit a characteristic relation where the strength exponent is about 0.75 times that of stiffness [35]. In practice, this framework is limited to explaining scaling exponents up to a value of 2 [34], necessitating its extension for more complex materials where higher exponents are often observed.

Since then, the scaling exponent has been associated with a broader range of drivers, including microstructural complexity [e.g., 36], the effective load-bearing volume fraction [e.g., 37], and the dominant deformation modes at smaller scales [e.g., 38,39]. This has focused research on the morphology of the porous material, aiming to explain its observed mechanical behavior through the geometric and topological properties of the microstructure. For example, the theoretical work by Huber [40] showed that mechanical properties can be described as a function combining geometry and topology.

In the literature, topology plays a crucial role in bridging mechanics and microstructure [41]. It can be described on a local scale, for example, through the coordination number $N_3$ of a node, or on a global scale using measures such as the genus $g$ — a measure for the number of “handles” or “loops” in the structure. While the local-scale parameter is easy to determine for a granular or cellular structure, it is difficult to obtain a meaningful parameter for materials that lack a clear unit cell definition, such as bicontinuous materials. Here, only a global perspective remains. The primary issue with topological descriptors such as the genus density ($g_v$) is that, as a dimensional property, it is not scale-invariant. Its value, therefore, changes if a given structure is scaled uniformly, even though the underlying topology is unchanged. To enable meaningful comparisons between materials with different microstructures and different scales, the genus must be scaled by a representative volume of the structure. This yields a dimensionless quantity known as the scaled genus density ($g_v^*$), which conceptually represents a mean coordination number and allows for direct comparison across different materials and scales [40,42,43].

## 1.2. Empirical background

Exploring the connection between the microstructure and mechanical properties of snow presents an empirical challenge, as many easily obtained microstructural parameters are highly inter-correlated and were found to have poor predictive power for the mechanical properties [e.g., 44]. For this reason, research in this area has recently been dominated by numerical studies. While experimental and analytical studies on this topic exist [e.g., 15,45,46], numerical simulations have the advantage that microstructural parameters can be well-controlled.

Although topology is a key predictor of mechanics in many porous materials [41], its application to snow presents unique challenges. Traditionally, and in many numerical models, snow is assumed to be a granular material. In reality, this assumption holds only for certain snow types, while many others — particularly complex morphologies
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Relating compressive strength and Young's modulus of snow — specifically in the brittle failure regime relevant for slab avalanche release to a wide range of technical and natural materials. Comparison data from the CES EduPack 2023, Grant a Design Limited, Cambridge, UK, 2023. The data for snow, firn and ice are taken from [16,18,23–28], the data for bone from [29,30], and the data for nanoporous gold from [31–33].</figcaption></figure>

such as depth hoar — exhibit a bicontinuous structure where individual grains are not clearly defined. This makes the segmentation of complex morphologies such as depth hoar very difficult [6,47–49]. As a result, µCT analyses often rely on bulk descriptors, which represent statistically averaged information over a large volume and cannot resolve granular-level features such as bond dimensions.

An exception to this is the minimum cut density approach by Hagenmuller et al. [50], which identifies the smallest surface whose removal disconnects two opposite faces of the structure [51]. Although it is not a formal topological invariant and does not resolve individual grains, the minimum cut density provides a topology-sensitive measure that reflects local bond-scale properties. However, the predictive power of such minimum cross-sectional area approaches for microstructures with large heterogeneity is debated [52] and it has not yet been systematically tested over a wide range of snow microstructures.

Due to these limitations in extracting topological information from snow samples, numerical modeling has played a key role in advancing our understanding of the connection between snow microstructure and mechanics. Stochastically generated, particle-based analogues offer unambiguous topological frameworks and allow independent tuning of volume fraction and coordination number. Löwe and Picard [53] derived sticky-hard-sphere (SHS) model parameters — ice volume fraction, sphere diameter, and stickiness — by fitting the µ CT-based two-point correlation function to the model's analytical form [54]. These were later linked to coordination number, connecting microstructure to contact mechanics [55]. The SHS framework was then used to construct discrete element models, where Gaume et al. [56] and Ritter et al. [57] showed that strength and stiffness scale with contact density $v_c = N_3\phi$, consistent with the geometric-topological scaling proposed by Huber [40].

More recently, Blatny et al. [58] extended numerical methods to bicontinuous structures, which better represent the true microstructure of snow. They used Gaussian random fields (GRFs) to generate statistically bicontinuous microstructures, and simulated their response under various loading conditions with the Material Point Method (MPM). These simulations produced scaling exponents of 3.22 and 3.86 for

the bulk modulus and compressive strength, respectively, and showed sensitivity to internal friction.

### 1.3. *Research gap*

Traditional density-based parameterizations have long been used to characterize the mechanical properties of weak snow layers. However, these models neglect critical geometric and topological features that govern load-bearing behavior and failure, leading to large uncertainties. While numerical studies have offered valuable insights, progress has been limited by the lack of experimental datasets linking the mechanical response to the three-dimensional microstructure.

In this work, we address this gap by introducing a dataset that pairs compressive mechanical tests on weak snow layers — both artificially grown and natural — with quantitative microstructural descriptors derived from µCT images. We adapt morphological metrics developed for other porous materials to characterize snow geometry and topology and statistically evaluate their relationship to Young's modulus and compressive strength. This analysis aims to identify the structural parameters that most strongly influence the mechanical and failure behavior of snow. This could provide a foundation for improved snow stability prediction and more accurate mechanical models in future snow cover simulations.

## 2. Methods

Fig. 3 provides an overview of our methodology. Weak layer samples were either artificially produced or collected from the field, scanned using µCT, and subjected to uniaxial compression tests. High-speed video recordings enabled the extraction of stress–strain curves and strain rates, from which the elastic modulus was estimated. Subsequently, microstructural features were quantified and statistically analyzed. A neural network was trained and used to determine the feature importance, identifying which features were most strongly associated with the mechanical behavior.
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Schematic workflow of the study, showing each method used, the principal results of every stage, and how those results feed into the subsequent steps of the analysis.</figcaption></figure>

## 2.1. Weak layer samples

To perform mechanical experiments, the weak layers need to be contained between dense snow layers. This allows for cutting and transporting the samples without damaging the weak layer and provides an interface for the mechanical tests. This sandwich-like configuration also helps to localize the failure within the weak layer and is therefore required for experimental testing of any type of weak layer. To quantify the effect of microstructure on the mechanical properties for different morphologies, we used 278 weak layer samples of three different types: faceted crystals and depth hoar (FC&DH - 114 samples), decomposing particles and rounded grains (DF&RG - 132), and buried surface hoar (53). Examples of the corresponding microstructures are shown in Fig. 1.

The starting point for our analysis was a “parent sample” — a large block of snow with an internal weak layer that was either created in the laboratory or collected from the natural snowpack. Each parent sample provided the source material for 5–8 rectangular samples (65 mm ×140 mm with variable thickness) intended for mechanical experiments and one specimen for µCT imaging.

For FC&DH and DF&RG, it would be very time-consuming to find natural samples with different microstructures in a sandwich-like configuration for mechanical testing. We therefore produced these weak layers artificially in the cold laboratory. In the present study, we used the same dataset on FC&DH as in [44] and complemented it by performing additional tests with high-density samples. To overcome the inherent data scatter in snow experiments, a sufficiently large number of samples is necessary. We achieved this by using elevated temperature gradients, which sped up production at the cost of creating more extreme microstructures. The methods to grow weak layers of faceted crystals and depth hoar were described in detail in [44]. Weak layers of DF&RG were produced by sieving different snow types in a sandwich-like configuration. These samples were then stored for between one

week and 6 months at temperatures between -2.5 °C and -20 °C to allow the grains to sinter.

Surface hoar, on the other hand, was difficult to produce in a laboratory environment. While we were able to grow surface hoar in a setup similar to [59], covering it with another layer to obtain a sandwich-like stratigraphy usually damaged the surface hoar severely. We therefore used natural weak layers of buried surface hoar for this study, which we harvested from the snowpack along a creek where surface hoar frequently forms. We sampled two different buried surface hoar layers at separate times in January 2024. One of the layers was located close to the snow surface and was covered by only a very thin slab. To achieve the desired slab thickness for our experiments, we sieved an additional layer of rounded grains on top.

## 2.2. Experimental testing

For the mechanical experiments, we followed the procedure previously described in [44]. Compression tests were conducted on a uniaxial testing machine with a 10 kN load cell (Shimadzu, Japan - sampling at 1 kHz and force accuracy within ±1% in the range of 20 N to 10 kN). This load cell was therefore well-suited for the experiments, as the observed sample failure loads (4.6 N to 10.1 kN, mean value of 415 N) fell mostly within its specified accuracy range. The compression plate can rotate ±2.5° to ensure a uniform load even if the sample surface is slightly misaligned. To account for the strain rate dependency of snow and to ensure that the experiments remained within the brittle regime, we maintained displacement rates between 0.1 mm s⁻¹ and 0.4 mm s⁻¹, yielding a theoretical strain rate of approximately 10⁻² s⁻¹ (for weak layer thicknesses of 10–40 mm, assuming the deformation is concentrated within the weak layer). The actual strain rate is analyzed later. As reported in our previous study [44], the signals of the compression experiments showed artifacts (pop-ins), due to localized failures at the sample-machine interface.
<figure><img src="image_7.png" /><figcaption>Fig. 4. Determining the lag λ between the independent stress and strain signals with the help of cross-correlation: (a) shows both signals in their reference time frame, (b) shows the absolute increments of the stress signal (|Δσᵢ|) and the strain signal (|Δεᵢ|), which were used for the cross-correlation, (c) shows the cross-correlation function R<sub>σε</sub>(λ) according to Eq. (2), the global maximum defines the lag of the two signals.</figcaption></figure>

We recorded all experiments with a Phantom Veo 710L high-speed camera (AMETEK Vision Research, USA - 1280 px × 720 px) at 1,000 frames per second. These videos allow us to assess the quality of the tests (e.g., tilting, failure outside the weak layer), but also to calculate the deformation within the sample using Digital Image Correlation (DIC). The camera was positioned both horizontally and vertically perpendicular to the sample to minimize perspective distortion. Lens distortions were negligible compared to other sources of error and the overall experimental scatter and were therefore not compensated in the analysis.

For the image correlation, we used the tracking mode of the open-source Digital Image Correlation Engine software [60]. The analysis is described in detail in [44]. To evaluate the strain in the weak layer, we selected narrow regions of interest (ROI) immediately below and above the weak layer and calculated the difference in the average displacement of the two regions. To calculate the strain, we used the vertical distance between the center points of the two ROIs. The strain signal showed the same pop-ins as observed in the stress signal, which helped in synchronizing both signals.

Due to limitations in our experimental setup, the strain and stress signals did not have a common time coordinate (see Fig. 4a). We therefore determined their lag (λ) by cross-correlating the increments (Δσᵢ = σᵢ - σᵢ⁻¹ and analogously for Δεᵢ) of the two signals (Fig. 4b):

$$R_{\sigma\epsilon}(\lambda) = \sum_{t=t_{\text{start}}}^{t_{\text{end}}} |\Delta\sigma_t| |\Delta\epsilon_{t+\lambda}| \quad (2)$$

where R<sub>σε</sub>(λ) is the cross-correlation function of the stress and strain increments at a given lag λ. A positive λ shifts the strain signal forward in time relative to stress. No normalization was applied because only the peak position is required for the alignment. The limits t<sub>start</sub> and t<sub>end</sub> were chosen manually to focus the analysis on the relevant timing of the failure, as we also recorded additional post-failure stress data. Because pop-ins occurred quasi-periodically [61], R<sub>σε</sub>(λ) often exhibited multiple peaks (Fig. 4c); we therefore accepted only those cases where a visual verification of the cross correlation success was possible. To this end, we looked at the moment of failure of the original signals, at which a correct alignment of both signals leads to opposing trends. The procedure yielded synchronized stress-strain curves for 241 of the 278 experiments.

We then used the resulting stress-strain plots to estimate the effective elastic modulus of our weak layer samples. To obtain more robust results, we also used the slopes of the pop-ins for this analysis (shown for one example in Fig. 5). From the 241 samples where we were able to correlate the stress and strain signal, we identified 987 regions showing

a consistent stress-strain response (pre-failure behavior and slopes of pop-ins). We identified each of these regions by their local minimum and maximum. The first region analyzed was always the slope leading to the ultimate strength of the samples. In addition, we selected regions earlier if pop-ins were present.

We estimated an effective elastic modulus, E<sub>eff</sub>, for each of the 987 stress-strain segments with an incremental window-growing regression. Starting with the first 10% of the stress-strain data points, we fit a regression line and computed the coefficient of determination, R². If R² ≥ 0.9, the slope was accepted as E<sub>eff</sub>. Otherwise, the window was enlarged in steps of one data point until the criterion was met or the window encompassed 50% of the segment, beyond which the test was discarded. This universal approach, based on the R² value, allowed us to reproducibly identify a sufficient linear region for most samples, while standard methods (e.g., the 0.2% method [62]) did not reliably work due to the limitations of our data in terms of resolution and scatter. We also calculated a strain rate ε<sub>DIC</sub> based on the final regression window.

For a more generalized comparison, our experimental results were normalized against reference mechanical properties of the bulk material. Specifically, the yield strength of polycrystalline ice (σ* ≈ 2 MPa) [18] served as the normalization factor for strength, while the Young's modulus of polycrystalline ice (E* ≈ 10 GPa) [63] was used for elasticity.

## 2.3. Microstructural analysis

We used µCT imaging (µCT 90, Scanco Medical, Switzerland - nominal resolution < 8 µm at 10% of the modulation transfer function) to quantify the microstructure of every parent sample (70 scans in total). As scan parameters, we chose a source current of 145 µA at 55 kV, an integration time of 300 ms for 1500 projections, and we averaged the data twice per projection. The reconstructed geometry was Gaussian filtered to reduce noise (width = 1.2 voxels, support = 2 voxels) and binary segmented [64]. We used three different sample holders (50 mm, 70 mm and 90 mm) resulting in voxel sizes of (16.2 µm)³, (22.8 µm)³, and (29.1 µm)³. The selection of different sample holders represented a necessary trade-off between maximizing imaging resolution and preserving the structural integrity of the fragile weak layer. Lower resolutions were generally only used for microstructures with a specific surface area (SSA<sub>ice</sub>) below a value of 20 mm⁻¹, to avoid the loss of structural detail [e.g., 65].

The µCT scans were utilized for a comprehensive analysis of the weak layer microstructure, encompassing the following:
<figure><img src="image_8.png" /><figcaption>Fig. 5. Example illustrating the determination of the elastic modulus from linear segments of the stress-strain curve. Distinct regions, identified by local maxima and minima, are color-coded together with their corresponding regression lines. The regression analysis started with the initial 10% of data points per region, incrementally including additional points until the criterion R² ≥ 0.9 was satisfied.</figcaption></figure>

(i) To calculate average (bulk) properties of the weak layer, we manually selected a cubic region of interest (ROI) containing only the weak layer. The size of the ROI was limited by either the thickness of the weak layer or the size of the sample holder. The minimum ROI size was 278 voxels ((8.09 mm)³), which we assumed to be a sufficiently large representative volume, similar to [17,66,67]. Typically, the ROI was much larger (> 1000 voxels). This ROI was then used to calculate a wide range of parameters: Volume fraction, representing the proportion of ice and air in the sample, was calculated based on the triangulated surface mesh ($\phi_{tri}$) and based on the voxel fraction ($\phi_{vox}$). The triangulated mesh was generated using a marching cubes algorithm [68]. For the calculation of the specific surface area (SSA) — a measure describing the surface-to-volume ratio of the structure — we again used the triangulated surface mesh. The total surface area was normalized by either the volume of ice phase (SSA$_{ice}$) or the volume of the pore space (SSA$_{pore}$). Furthermore, several microstructural parameters were calculated from the binary 3D data: chord length distributions ($l_c$), which describe the statistical lengths of straight lines within a single phase; connectivity density ($C_v$), a measure of the structure's connectiveness calculated from the Euler characteristic; and the exponential correlation length ($\xi$), which represents the characteristic size of the structural features.

(ii) To resolve heterogeneity within the weak layer, we performed a moving window analysis on each sample to generate high-resolution vertical profiles. Here, we chose a window size of 1000 × 1000 × 100 voxels (x, y, z) and a z-overlap of 50%. Volume fraction and SSA were calculated for each slice as described above. We used the profiles to find the minimum values of volume fraction and SSA (min ($\phi_{tri}$), min (SSA$_{ice}$)) of each weak layer sample [44].

(iii) We calculated further morphology metrics using additional software or custom python code:

First, we used the open-source software ImageJ for additional topology and geometric evaluations. We calculated Euler characteristics ($\chi$) and local thickness distributions ($l$), where local thickness is defined as the diameter of the largest sphere that fits within the structure and contains a given point with the BoneJ plugin for Fiji (ImageJ [69–71]).

Second, we used a plugin developed by Roque and Costa [72] for ImageJ [73] to calculate tortuosity ($\tau_{ice}$) of the ice structure (a measure for the complexity or sinuosity of paths through porous structures). This plugin uses a geodesic reconstruction algorithm, which gives robust estimates for $\tau_{ice}$ without the need to use skeletonization [74]. The algorithm calculates the geodesic distance $l_g$ of each pixel of the

microstructure relative to a starting plane and compares this against the Euclidean distance $l_e$. $\tau_{ice}$ is then estimated from the slope of the best fit line of $l_g$ as a function of $l_e$. The plugin allows selecting the number of neighboring pixels (6, 18, or 26) for geodesic distance calculation; we consistently chose 6 neighbors.

Third, we quantified the curvature of the ice-air interface of the microstructure. To this end, we used the triangulated surface meshes of all samples, which we analyzed using the VTK library [75]. The mean curvature $H$ was computed as a weighted average of dihedral angles along the edges connected to each vertex, capturing how much a vertex deviates from its neighbors. The Gaussian curvature $K$ was calculated using the angle defect method, which measures deviation from flatness by summing the angles of all adjacent faces at a vertex and normalizing the difference from $2\pi$ by the associated area. From $H$ and $K$, we derived the principal curvatures $\kappa_1 > \kappa_2$:

$$ \kappa_1 = H + \sqrt{H^2 - K}, \quad \kappa_2 = H - \sqrt{H^2 - K} \qquad (3) $$

These are the extrema of normal curvature in orthogonal directions. The principal curvatures enable us to classify points on the surface mesh into three distinct categories: concave regions ($\gamma_{conc}$: $\kappa_1 < 0$, $\kappa_2 < 0$), convex regions ($\gamma_{conv}$: $\kappa_1 > 0$, $\kappa_2 > 0$), and saddle-shaped regions ($\gamma_{saddle}$: $\kappa_1 > 0$, $\kappa_2 < 0$). Using these definitions, we classified all points on the surface mesh and computed the fractions corresponding to each category. As we did not employ adaptive mesh refinement, we assumed that the element-based fractions provide a representative distribution throughout the sample.

Fourth, we used the algorithm of Hagenmuller et al. [50] to find the minimum cut plane. This method uses a network-based heuristic to identify the minimal cross-sectional ice area that, if removed, would split the snow sample into two disconnected parts. The resulting minimum cut plane allows for the calculation of the cut area ($M_{c,A}$) and the number of bonds ($M_{c,N}$) that need to be removed to achieve this separation.

Building upon the previously summarized microstructural parameters, which are direct measurements of the microstructure, we can derive higher-level descriptors. These parameters often result from combining or normalizing basic measurements to capture invariant aspects of the microstructure. Such descriptors can provide a more comprehensive understanding of material properties. Grain size is a common microstructure parameter, recorded in every traditional snow profile. Although defining a discrete grain size is conceptually paradoxical for a bicontinuous medium such as snow, the parameter represents a
characteristic length scale of the material and has thus implications for remote sensing and mechanical considerations. Determining the grain size exactly for all types of snow microstructures from µCT data is non-trivial; therefore, it is usually simplified as optical grain size $d_{opt}$, which assumes spherical grains. In theory it can be calculated both from SSA$_{ice}$ or the mean chord length $\bar{l}_{c,ice}$ [42,76]:

$$d_{opt} = \frac{6}{\text{SSA}_{\text{ice}}} = \frac{3\bar{l}_{c,ice}}{2} \qquad (4)$$

The genus *g* is a dimensionless topological descriptor that counts the number of through-going tunnels (holes or loops) in a material's microstructure [77]. In three-dimensional porous media, it is typically derived from the Euler characteristic $\chi$. For a fully connected structure that contains no isolated pores (= bicontinuous structure), the genus per unit volume $g_v$ equals the connectivity density $C_v$. These parameters describe the topology of the microstructure, but they are dimensional and not scale-invariant and therefore not suitable to compare topology across different microstructures. To obtain a true, scale-invariant metric for topology, it is customary to normalize $g_v(\equiv C_v)$ by an appropriate characteristic volume of the microstructure [e.g., 42,43], yielding the scaled genus density ($g_v^*$) [40,78–80]:

$$g_v^*[\ell] \propto C_v (\ell_{\text{rep}})^3 \qquad (5)$$

In the context of snow microstructure, there are several length scales $\ell_{\text{rep}}$ that would be suited to calculate $g_v^*$, such as the correlation and chord lengths, ligament- and pore-size distributions, the inverse specific surface area [81], and the inverse principal curvature [82], which we will evaluate in the results section of this study. One more high-level metric we used was anisotropy. Snow anisotropy is commonly defined by the relation of vertical correlation length to horizontal correlation lengths according to the general formulation [83]:

$$\alpha = \frac{\ell^z}{0.5(\ell^x + \ell^y)} = \frac{\ell^z}{\ell^{xy}} \qquad (6)$$

In addition to this metric, we similarly calculated anisotropy based on the minimum cut area, minimum cut bond number, chord lengths, and tortuosity. All the parameters used in this analysis are summarized in Table 1.

## 2.4. Data-driven analysis

To identify the microstructural parameters that govern macroscopic stiffness and strength across all considered grain types, we derived feature importance scores from a neural network. The workflow consisted of four steps: (i) data standardization, (ii) supervised model training with hyperparameter tuning using a grid search, (iii) correlation-based clustering to handle multicollinearity, and (iv) permutation-importance scoring within each cluster. The workflow was implemented in a Jupyter Notebook with scikit-learn [84] and scipy [85].

(i) Data and targets: 44 geometric and topological descriptors extracted from µCT were z-standardized (mean 0, variance 1; [86]) and used to predict the parent-sample-averaged elastic modulus and compressive strength. The target variables were log-transformed to reduce the disproportionate influence of the largest values.

(ii) Model training: We chose a fully connected multilayer perceptron (MLP) because it can approximate arbitrary continuous mappings yet remains interpretable with post-hoc tools [87]. Architectures comprising one to three hidden layers, each containing 10, 20, 50 or 100 neurons and using identity, tanh or relu activations, were explored. Hyperparameters were selected by a grid search embedded in a 10 × 10 repeated k-fold cross-validation [88]; the best configuration was then refitted to 80% of the data and evaluated on the remaining 20% hold-out set. As a benchmark metric during training, we used the Mean Absolute Percentage Error (MAPE), which quantifies the average magnitude of prediction errors relative to the actual values [89]. MAPE

<caption>**Table 1**</caption>

List of microstructural descriptors with abbreviations, full names, and data sources. Lucid Xamflow (XF) is our standard µCT evaluation software.

<table><thead><tr><th>Abbr.</th><th>Full name</th><th>Source</th></tr></thead><tbody><tr><td>Φ<sub>tri</sub></td><td>triangulated bulk volume fraction</td><td>XF</td></tr><tr><td>Φ<sub>vox</sub></td><td>voxel-based bulk volume fraction</td><td>XF</td></tr><tr><td>min (Φ<sub>tri</sub>)</td><td>minimal profile-based volume fraction</td><td>XF</td></tr><tr><td>SSA<sub>ice</sub></td><td>Bulk SSA norm. by ice volume</td><td>XF</td></tr><tr><td>min (SSA<sub>ice</sub>)</td><td>minimal SSA profile</td><td>XF</td></tr><tr><td>SSA<sub>pore</sub></td><td>bulk SSA norm. by pore volume</td><td>XF</td></tr><tr><td>ξ<sup>z</sup></td><td>vertical exponential correlation length</td><td>XF</td></tr><tr><td>ξ<sup>xy</sup></td><td>mean horizontal exponential correlation length</td><td>XF</td></tr><tr><td>α<sub>z</sub></td><td>vertical anisotropy based on ξ</td><td>Eq. (6)</td></tr><tr><td>l̄<sub>ice</sub></td><td>mean chord length ice</td><td>Krol and Löwe [42]</td></tr><tr><td>l̄<sup>z</sup><sub>ice</sub></td><td>mean vertical chord length ice</td><td>Krol and Löwe [42]</td></tr><tr><td>l̄<sup>xy</sup><sub>ice</sub></td><td>mean horizontal chord length ice</td><td>Krol and Löwe [42]</td></tr><tr><td>α<sub>c,ice</sub></td><td>vertical anisotropy, based on l̄<sub>c,ice</sub></td><td>Eq. (6)</td></tr><tr><td>l̄<sub>c,ice</sub></td><td>mean chord length pore</td><td>Krol and Löwe [42]</td></tr><tr><td>l̄<sup>z</sup><sub>c,ice</sub></td><td>mean vertical chord length pore</td><td>Krol and Löwe [42]</td></tr><tr><td>l̄<sup>xy</sup><sub>c,ice</sub></td><td>mean horizontal chord length pore</td><td>Krol and Löwe [42]</td></tr><tr><td>α<sub>pore</sub></td><td>vertical anisotropy based on l̄<sub>pore</sub></td><td>Eq. (6)</td></tr><tr><td>C<sub>v</sub></td><td>connectivity density</td><td>XF,[70]</td></tr><tr><td>ξ<sup>SSA</sup></td><td>scaled genus density calculated using SSA<sub>ice</sub></td><td>Eq. (5)</td></tr><tr><td>l̄<sub>ice</sub></td><td>mean local ice thickness</td><td>Domander et al. [70]</td></tr><tr><td>std (t<sub>ice</sub>)</td><td>standard deviation of t<sub>ice</sub></td><td>Domander et al. [70]</td></tr><tr><td>max (t<sub>ice</sub>)</td><td>maximal t<sub>ice</sub></td><td>Domander et al. [70]</td></tr><tr><td>l̄<sub>pore</sub></td><td>mean local pore thickness</td><td>Domander et al. [70]</td></tr><tr><td>std (t<sub>pore</sub>)</td><td>standard deviation of t<sub>pore</sub></td><td>Domander et al. [70]</td></tr><tr><td>max (t<sub>pore</sub>)</td><td>maximal t<sub>pore</sub></td><td>Domander et al. [70]</td></tr><tr><td>χ</td><td>corrected Euler characteristics of the fully-connected structure</td><td>Domander et al. [70]</td></tr><tr><td>H</td><td>mean interface curvature</td><td>Sullivan and Kaszynski [75]</td></tr><tr><td>K</td><td>Gaussian interface curvature</td><td>Sullivan and Kaszynski [75]</td></tr><tr><td>κ̄<sub>1</sub></td><td>mean first principal curvature</td><td>Eq. (3)</td></tr><tr><td>std (κ<sub>1</sub>)</td><td>standard deviation of κ<sub>1</sub></td><td>Eq. (3)</td></tr><tr><td>κ̄<sub>2</sub></td><td>mean second principal curvature</td><td>Eq. (3)</td></tr><tr><td>std (κ<sub>2</sub>)</td><td>standard deviation of κ<sub>2</sub></td><td>Eq. (3)</td></tr><tr><td>γ<sub>conv</sub></td><td>convex interfacial fraction</td><td>Eq. (3)</td></tr><tr><td>γ<sub>conc</sub></td><td>concave interfacial fraction</td><td>Eq. (3)</td></tr><tr><td>γ<sub>saddle</sub></td><td>saddle interfacial fraction</td><td>Eq. (3)</td></tr><tr><td>τ<sup>z</sup></td><td>vertical tortuosity</td><td>Roque and Costa [72]</td></tr><tr><td>τ̄<sup>xy</sup></td><td>mean horizontal tortuosity</td><td>Roque and Costa [72]</td></tr><tr><td>α<sub>z</sub></td><td>vertical anisotropy based on τ</td><td>Eq. (6)</td></tr><tr><td>M<sup>z</sup><sub>c,A</sub></td><td>vertical MinCut bond area per m<sup>2</sup></td><td>Hagenmuller et al. [50]</td></tr><tr><td>M<sup>xy</sup><sub>c,A</sub></td><td>mean horizontal MinCut bond area per m<sup>2</sup></td><td>Hagenmuller et al. [50]</td></tr><tr><td>α<sub>M<sup>z</sup><sub>c,A</sub></sub></td><td>vertical anisotropy based on M<sup>z</sup><sub>c,A</sub></td><td>Eq. (6)</td></tr><tr><td>M<sup>z</sup><sub>c,N</sub></td><td>vertical MinCut bond number per m<sup>2</sup></td><td>Hagenmuller et al. [50]</td></tr><tr><td>M<sup>xy</sup><sub>c,N</sub></td><td>mean horizontal MinCut bond number per m<sup>2</sup></td><td>Hagenmuller et al. [50]</td></tr><tr><td>α<sub>M<sup>xy</sup><sub>c,N</sub></sub></td><td>vertical anisotropy based on M<sup>xy</sup><sub>c,N</sub></td><td>Eq. (6)</td></tr></tbody></table>

was chosen because it is less sensitive to outliers than other common error metrics.

(iii) Handling multicollinearity: To avoid inflating importance scores, we grouped highly correlated predictors before interpretation. Absolute Spearman rank correlations (|r<sub>s</sub>|) were converted to distances $d = 1 - |r_s|$ and hierarchically clustered with average linkage [90]. Cutting the dendrogram at $d = 0.25$ yielded clusters whose members share $|r_s| \ge 0.75$.

(iv) Permutation importance and ranking: We estimated the contribution of every descriptor with permutation importance [91] (cross-validation with 50 splits and 1000 permutations per fold on the hold-out set). Importance values of variables belonging to the same cluster were averaged, producing a ranked list of microstructural motifs —
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Mean strain rate per parent sample from our experiments plotted in temporal order, obtained through Digital Image Correlation (DIC). The color corresponds to the weak layer category FC&DH (yellow), DF&RG (blue) and SH (purple). The vertical error bars represent the observed standard deviations per parent sample. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

rather than single, collinear descriptors — that most strongly influence
the predicted mechanical response.

3. Results

We conducted displacement-controlled compression tests on 278 weak layer samples, including faceted crystals and depth hoar (114 samples, of which 92 were already published in [44]), precipitation particles and rounded grains (132 samples), and buried surface hoar (53 samples). In parallel, we performed 70 µCT scans to obtain detailed microstructural information for the samples. The combined dataset was then used to find the microstructural features associated with the mechanical response. For clarity, we present the results in three sections: experimental results, microstructural results, and statistical results. These are discussed in the same order.

3.1. *Experimental results*

To verify the experimental conditions, we first address **strain rate**, a parameter to which snow's mechanical properties are highly sensitive. Schöttner et al. [44] found a large discrepancy between the displacement rate applied by the testing machine and the actual strain rate observed within the weak layer. Using DIC on all our samples, we found that the mean strain rate during our experiments was relatively consistent across all sample types, averaging 2.1(±0.9) ×10⁻³ s⁻¹ (Fig. 6).

The effective elastic behavior derived from the DIC analysis (Fig. 7a) revealed distinct relations on volume fraction for different types of weak layers. Notably, faceted crystals and depth hoar (FC&DH) exhibited a stiffer effective elastic behavior compared to decomposing particles and rounded grains (DF&RG) at comparable densities. Our experimental results are best described using a power law relation with the bulk weak layer volume fraction, which we fitted using an orthogonal distance regression (ODR - [92]). The largest power law exponent, 5.1, was observed for faceted crystals and depth hoar. For decomposing and fragmented particles and rounded grains, the exponent was slightly lower (4.6), and their curve consistently was below that of depth hoar, although the difference is within the margin of error. In contrast, buried surface hoar exhibited a different behavior, characterized by a shallower slope with a power law exponent of 1.7. For the effective elastic modulus, Fig. 7a includes the data from both the initial pop-ins and the slope leading to the ultimate strength. If we only considered the elastic regions of the slope leading to ultimate strength (see Fig.

A.1), we found similar results, but with more scatter in the data. The quality of the fits was quantified using the reduced chi-squared value, $\chi_{\text{red}}^2$. The SH layers yielded a $\chi_{\text{red}}^2$ of 0.2 — substantially below the ideal value of 1 — indicating that the experimental uncertainties are relatively large compared to the well-defined linear trend of the mean values. In contrast, FC&DH and DF&RG layers had $\chi_{\text{red}}^2$ values of 0.8 and 2.2, respectively, suggesting a more balanced correspondence between the orthogonal distance regression fits and the data scatter.

To validate our approach of calculating the effective elastic modulus from the slopes of the pop-ins, we investigated whether this introduced any systematic errors. The example in Fig. 5 showed an increasing trend in the elastic modulus for each pop-in. To determine whether this was a systematic effect, we examined how the elastic modulus evolved with the number of pop-ins prior to failure. Specifically, we normalized each pop-in's elastic modulus and strain rate by the mean value of all pop-ins in a given sample and then averaged these normalized values across all samples. As shown in Fig. A.9, no clear trend emerged for the elastic modulus when considering all samples. The strain rate, however, increased for each pop-in. We therefore conclude that our approach — using the slopes of pop-ins to obtain additional data points and thereby reduce scatter — is valid.

The results for the **ultimate compressive strength** (Fig. 7b) also exhibited a distinct scaling for different crystal types. Faceted crystals and depth hoar again exhibited the largest exponent, characterized by a power law exponent of 5.7. However, in contrast to the trends seen for the elastic modulus, the power law exponents for rounded grains and buried surface hoar were similar, with exponents of 2.8 and 2.1, respectively. The lines of best fit for these two crystal types were closely aligned, differing by only a small vertical offset. The quality of the fits for the ultimate strength data was reasonable, with $\chi_{\text{red}}^2$ values ranging from 1.4 to 3.0. These values suggest that the error bars may be slightly underestimated, as the observed scatter in the data is not fully captured by the current error bars.

Further insight is gained from the material's **scaling behavior**. Rather than independent data points, the results regarding stiffness and strength reveal fundamental differences in the mechanical behavior of each grain category. This is best illustrated by plotting normalized strength versus normalized elastic modulus, thereby removing the influence of volume fraction. As shown in Fig. 8, two distinct behaviors emerge: the persistent grain types FC&DH and SH showed steep slopes of 1.1 and 1.2, respectively, while the DF&RG was characterized by a much shallower slope of 0.65.

3.2. *Microstructure results*

We extracted 44 microstructural descriptors for each of the 70 µCT scans. Among these, the minimum cut bond area — representing the minimal area that, if removed, separates the sample in the vertical direction — emerged as the single best predictor of mechanical properties. Remarkably, this parameter enables a collapse of the mechanical response across different grain types onto a single curve (Fig. 9). This is likely because, while both strength and minimum cut bond area exhibit distinct dependencies on volume fraction for each grain type (see Fig. A.8, their respective variations appear to compensate for each other when strength is related to minimum cut bond area, resulting in this unified behavior. Regarding the strength scaling, Fig. 9b) reveals two distinct behaviors: a plateau in strength for $M_{c,A}^z$ values below approximately 1%, followed by a power-law relationship at higher values.

While other microstructural metrics also showed interesting relations, we do not present them individually here. Many of these parameters were strongly inter-correlated [44], making it challenging to isolate their individual contributions to the mechanical behavior. For completeness, the most relevant of these additional descriptors are included in the Appendix. To facilitate comparison, we express each of them as a function of volume fraction: SSA (Fig. A.2), correlation
<figure><img src="image_10.png" /><figcaption>Fig. 7. (a) Elastic modulus and (b) compressive strength as a function of volume fraction and density. The different colors correspond to the three groups of grain types. The values are normalized with the Young's modulus (10 GPa) and the tensile strength (2 MPa) of polycrystalline ice. The vertical error bars represent the standard deviation of the experimental results, the horizontal error bars are an estimate for the spatial variability of volume fraction within a parent sample [44].</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Normalized strength vs normalized elastic modulus shown for our experimental data. The values are normalized with the elastic modulus and tensile strength of polycrystalline ice (see Section 2.2). Additionally we show the power laws of Gaume et al. [56], Blatny et al. [58]. For the latter we transformed the bulk modulus power law assuming a poison ratio of 0.26.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Elastic modulus (Panel a) and compressive strength (Panel b) vs the vertical minimum cut area $M_{c,A}^z$. The vertical error bars represent the standard deviation of the experimental results. The Mean Absolute Percentage Error (MAPE) of the presented fits is 53.2% for Young's modulus and 63.6% for strength.</figcaption></figure>

lengths and anisotropy (Fig. A.3), interface curvature (Fig. A.4) and
curvature classification (Fig. A.5), local thickness (Fig. A.6), tortuosity
(Fig. A.7), and minimum cut density (Fig. A.8). The interested reader
can find more details in the supplemental materials.

3.3. *The role of topology*

To investigate the influence of topology on snow microstructure and
mechanical behavior, it is essential to identify a scale-invariant topolog-
ical descriptor, namely the scaled genus density $g_v^*$. A typical approach
in the literature is to adopt the inverse specific surface area (1/SSA$_{ice}$
— see Eq. (4)) as the characteristic length scale [e.g., 81,93]. Using
this, we computed $g_{v,SSA}^*$ according to Eq. (5), which yielded an almost
perfectly monotonic trend with volume fraction for the persistent grain
types (Fig. 10). The DF&RG samples, however, were well off this curve.
While this is promising in separating the different crystal types, it did
not allow for a convincing data collapse of our experimental mechanical
data as a function of volume fraction and topology. This remained true
even when the $g_v^*$ was calculated using other microstructural length
metrics from our dataset, including correlation- and chord lengths,
ligament- and pore-size distributions, the inverse specific surface area,
and the inverse principal curvature [82].

To further investigate the relationship between topology and mi-
crostructure, we approached the problem from a different perspective.
Rather than seeking a length scale that would allow the scaled genus
density to explain the mechanical results, we instead aimed to identify a
length scale that collapses the scaled genus density $g_v^*$ across all grain
types onto a single master curve. Such a collapse would indicate the
length scale most representative of the internal structure — a proxy
for the size of a hypothetical structural unit cell. To this end, we re-
evaluated all available microstructural length scales in our dataset. For
each candidate, we computed the corresponding scaled genus density
and assessed the quality of the resulting collapse using the coefficient of
determination ($R^2$) from a global log-log power law fit across the grain-
type groups. As shown in Fig. 11a, the best collapse was achieved by
the mean chord length of ice in the vertical direction $\bar{l}_{c,ice}^z$, followed by
the mean overall chord length of ice $\bar{l}_{c,ice}$ and the mean local thickness

$\bar{l}_{ice}$. The corresponding master curve obtained using the best-performing
metric is shown in Fig. 11b.

3.4. *Data-driven exploration*

We used a neural network to identify the features that most strongly
control the mechanical behavior of weak layers and to predict mechani-
cal properties from microstructural descriptors. For the elastic modulus,
the search for the hyperparameters settled on a relatively small network
— three hidden layers of 10 neurons each with tanh activation — yet
cross-validation (CV) still gave a relatively high MAPE (≈ 71 (±62) %),
even though the held-out test error was lower (≈ 60%). For peak stress,
a slightly smaller architecture — two hidden layers of 10 neurons
each with relu activation — proved optimal, yielding a CV MAPE of
≈ 65 (±35) % and a test MAPE of ≈ 40%.

We grouped features into a small number of non-overlapping clus-
ters based on their pairwise Spearman correlation, using a distance
threshold of 0.25 (i.e. $|r_s| \ge 0.75$ — see (Fig. 12a). For each cluster,
the mean permutation importance across all its features was computed
and used as the cluster-level importance.

The resulting cluster importance shows that the elastic modulus
(Fig. 12b) was primarily influenced by Cluster 4, which was defined
by anisotropies calculated from the minimum cut density, along with
vertical tortuosity. In contrast, compressive strength (Fig. 12c) was
associated with a broader range of clusters exhibiting relatively high
permutation importance. The most influential was Cluster 6, which
included various interfacial curvature parameters, followed by vertical
tortuosity. Cluster 1 and the saddle fraction also ranked relatively high
in importance.

Our main analysis identified global relationships across the dataset,
which can obscure trends within distinct subsets like different grain
types. To investigate these specific behaviors, we performed a separate
analysis for each grain type. The results, presented in Fig. A.11, show
that each subset yields different clustering patterns and variable impor-
tance rankings. Especially the different clustering makes it difficult to
draw further conclusions. The combination of these divergent results
and the subjectivity of the underlying grain classification reinforces
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Scaled genus density, $g_{SSA}^* = C_v \times SSA_{ice}^{-3}$, obtained by normalizing the connectivity density with the cube of the inverse specific-surface area of ice. The horizontal error bars are an estimate for the spatial variability of volume fraction within a parent sample [44].</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Finding a representative length scale of topology in snow. Panel (a) shows $R^2$ ranking of length-scale candidates for collapsing $g_v^*$ of the different weak layer types on a single curve. Panel (b) shows the best collapse using the mean vertical chord length of the ice phase $\bar{l}_{c,ice}^z$. The horizontal error bars are an estimate for the spatial variability of volume fraction within a parent sample [44].</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Resulting feature importance based on the clustered features. Panel (a) shows a dendrogram of the distance (1-correlation) of all features, the clusters based on a distance threshold of 0.25. These clusters were then used to calculate an aggregated feature importance for each cluster for (b) elastic modulus and (c) compressive strength.</figcaption></figure>

our focus on identifying universal relationships that are independent of grain type.

We also tested other state-of-the-art statistical methods to understand the link between microstructure and mechanics, such as principal component analysis (PCA - Fig. A.12) and polynomial regression. However, these methods performed poorly, likely due to the strong multicollinearity present in our dataset.

## 4. Discussion

We performed 278 compression tests on weak layers to explore the link between microstructure and mechanical behavior. The results show that the mechanical properties relate differently to volume fraction depending on the grain type category. The observed scatter in mechanical properties spans roughly one order of magnitude and can be partly attributed to differences between grain type categories. Within a single category, the variability is smaller and consistent with the scatter reported in previous mechanical studies on snow [e.g., 14,94].

### 4.1. *Experimental results*

For the **elastic modulus in compression**, FC&DH and DF&RG layers exhibited volume fraction power law exponents of 5.1 and 4.6, respectively, whereas SH layers showed a much lower exponent of 1.7. The quality of the fits was reasonable with $\chi_{\text{red}}^2$ values in the range of 0.2 to 2.2.

A comparison with literature data shows consistent trends. Fig. 13 compares our measured elastic moduli with predictions from Gerling et al. [20] (Panel a) and Sundu et al. [19] (Panel b). In both cases, values are systematically offset from the 1:1 line, which is expected for snow given the well-documented discrepancies between the observed

relation between stress and strain, and the theoretical Young's moduli [e.g.94]. These discrepancies likely arise from non-ideal behaviors such as minor viscoplasticity during loading. The resulting apparent stiffness is termed the “effective elastic modulus”, $E_{\text{eff}}$, a necessary simplification of snow's time dependent behavior that enables the application of non-time-dependent models to snow mechanics [e.g., 9].

In particular, [20] reported scaling exponents between 4.0 and 4.6 for rounded-grain snow using a FEM-based approach combined with µCT scans, and acoustic measurements. When comparing their FEM-based power law to our data, we observed not only the previously mentioned offset in elastic modulus but also an additional discrepancy between FC&DH and DF&RG layers of approximately a factor of 2. Also, because the power law exponents differ, the power law fits are not parallel to the 1:1 line (Fig. 13a). Layers with buried SH showed a much larger exponent, indicating that the density-based parameterization overestimated the effect of density on the elastic modulus.

In contrast, employing the parameterization by [19] resulted in a better alignment. As shown in Fig. 13b, the fitting lines for FC&DH and DF&RG now coincide and run nearly parallel to the 1:1 line, albeit with an offset of roughly a factor of 8. Although this parameterization corrects the slope discrepancy observed with the density-based parameterization — highlighting the role of anisotropy — scatter within each grain type category (reflected by the $R^2$ values) remains unchanged. As with the density-based parameterization, the SH data points remain closely aligned with the 1:1 line. A plausible explanation for this observation is, that SH weak layers often incorporate smaller grains from the subsequent snowfall (see Fig. 1). These grains may increase the overall volume fraction, even though only the large SH crystals carry the load. Therefore, the bulk volume fraction may be an inaccurate input parameter for the parameterizations.

For the **ultimate strength in compression**, DF&RG layers exhibited a volume fraction lower power law exponent (2.8) compared to
<figure><img src="image_16.png" /><figcaption>Fig. 13. Comparison of our experimentally determined elastic moduli with the two most recent parameterizations in the literature. Panel (a) shows the comparison with the density-based parameterization of [20], panel (b) shows the comparison with the volume fraction- and anisotropy-based parameterization of [19]. The error bars represent the standard deviation of the experimental results.</figcaption></figure>

FC&DH (5.7), which is slightly higher than the previously reported value (5.4; [44]) due to the inclusion of additional high-density data. The relation for SH layers had a similar exponent of 2.3. The quality of the fits was reasonable, with $\chi^2_{red}$ values ranging from 1.4 to 3.0, meaning that the error estimate did not fully capture the scatter in the data and potential natural variability.

A detailed comparison for FC&DH with literature data was provided in [44]. In brief, our strength data agreed well with most previously reported values. For example, Reiweger and Schweizer [95] reported compressive strengths of 2.5–6 kPa for depth hoar layers with a density of 190 kg m⁻³, which is close to our results despite differing loading conditions. Chandel et al. [96] found strengths of 0.3–1.4 kPa for faceted crystals at lower densities (100–200 kg m⁻³), also consistent with our data. Numerical studies (e.g., Chandel et al. [96], Mede et al. [97]) reported strengths ranging from 2 to 7.5 kPa, though they found a weaker density dependence than observed in our experiments. Overall, despite methodological differences, the literature supports our findings for FC&DH layers.

Many studies reported compressive strength data for rounded grains, but variations in strain rate, temperature, and experimental setup limit direct comparability. Wang et al. [98] measured strengths around 100 kPa at densities of 350 kg m⁻³ and loading rates similar to ours. This is substantially higher than our DF&RG values (~30 kPa from our fit), but closer to our FC&DH values. This likely reflects differences in sample preparation and testing conditions. Kinosita [24] reported a density-scaling exponent of 3, consistent with our DF&RG results. Similarly, Smith [23] found strengths of ~70 kPa at 400 kg m⁻³, in good agreement with our fit.

For the compressive strength of SH, two experimental studies exist in the literature. Chandel et al. [99] found the compressive strength of artificially grown SH to be in the range of 1.0 to 2.3 kPa (with a density estimated at approximately 40 kg m⁻³). Similar values (although no pure compression tests were performed) were reported by Reiweger

et al. [100] for weak layers with densities between 60 and 70 kg m⁻³. The values of both studies are close to our results.

Two studies directly comparing different weak layer types offer context for the overall trends in our results. Hagenmuller et al. [101] found that microstructures formed under temperature-gradient metabolism were less resistant to compression than those formed under isothermal conditions, even at equal densities. In contrast, our data suggest an overall opposite trend. However, at low volume fractions (below 0.2), we observed similar strength values for both microstructure types. The power law behavior further indicates that, at such low volume fractions, isothermal samples (DF&RG) may even be stronger than temperature-gradient samples (FC&DH). Jamieson and Johnston [14], though focused on shear experiments, reported a comparable trend: a larger power law exponent for persistent grain types than for non-persistent ones. However, their exponents were substantially lower, pointing to a weaker density dependence in shear compared to compression.

The strain rates measured via DIC agreed with the values estimated by [44]. Fig. 6 shows that the strain rates of all of our samples were comparable with most samples exhibiting strain rates above 10⁻³ s⁻¹, which was one order of magnitude lower than expected based on the displacement rate of the testing machine. This highlights the need for detailed strain measurements of the weak layer. Our findings suggest that strain rate estimates from experimental studies in the literature, when not based on DIC measurements, are likely biased.

The average strain rate across all samples was 2.1×10⁻³ s⁻¹, suggesting the experiments were within the brittle failure regime but close to the brittle-to-ductile transition [61]. Notably, the FC&DH samples showed slightly lower strain rates than this average, a behavior we attribute to the faceting of the slab layers. Although the overall conditions were brittle, some minor plastic (ductile) effects were observed, such as irreversible strain offsets in the stress-strain curves (Fig. 5). We also observed a change in the strain rate with every pop-in. Fig.
A.9 shows that the measured strain rate increased with every pop-in by about 10%. This may be due to plastic processes at the sample-machine interface, leading to a more direct transmission of deformation towards the weak layer.

The presence of these ductile components suggests that faster loading rates would be needed to obtain purely brittle failure. However, the limited temporal resolution of our testing equipment means that higher rates would drastically reduce the stress–strain signal’s resolution, potentially leading to increased data scatter.

We also investigated the strain upon failure to search for differences in the behavior of the different types of weak layers. Fig. A.10 shows the strain in the weak layer when the ultimate strength was reached. We found an average value of 0.2%, and no clear differences for different weak layer types. We also identified two samples with very high strain values at failure, which we double-checked for irregularities, but we did not find a clear cause. The observed strain at failure agreed with the value reported by Chandel et al. [99] for SH.

## 4.2. Data-driven exploration

Identifying the microstructural features that govern mechanical behavior proved challenging, as many of them are highly correlated. This high degree of multicollinearity rendered simple statistical methods ineffective. A principal component analysis (Fig. A.12), for example, showed that while six components explained 95% of the variance, their loadings were broadly distributed, making physical interpretation impossible. To overcome this, we employed a machine learning approach and addressed the multicollinearity by clustering features into correlated groups. This method improved both model robustness and the interpretability of feature importance, enabling our model to reveal that distinct structural aspects govern the different mechanical properties, summarized in Fig. 12.

A methodological point to consider is that our analysis was performed on the dataset as a whole. Consequently, the results reflect the average correlation structure across all data points, aiming to find links between mechanics and microstructure fundamental to snow. This, however, can mask relations and can differ from the feature relationships that exist within distinct subsets (i.e. grain types) of the data. The results must therefore be interpreted with care. For instance, while volume fraction is a known predictor of stiffness [19], it ranked low in our global analysis. This discrepancy arises because the strong relationship between volume fraction and stiffness is grain type-dependent and becomes obscured when all data are combined. This is also true for other features which showed a grain-dependent behavior such as minimum cut density (Fig. 9), tortuosity (Fig. A.7), curvature (Fig. A.5), and local thickness (Fig. A.6). For practical applications, we suggest focusing on the most dominant clusters identified in the main analysis, while also considering potential grain type-specific behavior that may not be fully captured by global correlation metrics.

The **elastic modulus** in our model was primarily governed by Cluster 4, which comprised anisotropies derived from the minimum cut density. The second most influential feature was vertical tortuosity. Both anisotropy and tortuosity are physically plausible predictors, as they are representative for the alignment and efficiency of load-bearing paths in the material. The significance of these parameters for elastic modulus is well-supported in the literature [19,102].

For **compressive strength**, our analysis identified Cluster 6 — encompassing various interfacial curvature metrics — and vertical tortuosity as the most critical features. Interestingly, several other clusters also exhibited relatively high feature importance, indicating that compressive strength is influenced by a broader range of microstructural characteristics. From a geometrical standpoint, the importance of Cluster 6 may stem from its role as a proxy for the snow’s bonding network. Interfacial curvatures, which vary significantly between grain types, likely scale with geometric features of the grain and bond structure [103]. We suspect that while curvatures may not directly dictate

mechanical behavior, they are likely correlated with more fundamental topological and geometric aspects of the microstructure. Likewise, vertical tortuosity characterizes the linearity of load paths through the material; more tortuous paths may induce unfavorable bending and stress concentrations, ultimately reducing the overall strength. Collectively, these findings highlight the dominant role of the geometry of the bonding system in determining the compressive strength of weak snow layers.

For both stiffness and strength, it is noteworthy that classical **topological metrics**, such as the Euler characteristic or scaled genus density, were not among the most predictive features. This finding is unexpected, as such parameters are known to be highly important in other porous materials [41]. We propose two possible interpretations for why geometric features related to the bonding system proved more predictive than topological ones in this analysis. First, our topology metrics may not adequately capture the relevant aspects of structural connectivity. This may be due to an unrepresentative choice of length scale in the scaled genus density, limitations in the size of the considered microstructure volume to be representative for topology [40], or the presence of artifacts that distort topological measurements. For instance, in SH samples, we often observed many tiny particles between the large surface hoar crystals, which influenced the topological metrics strongly but are likely irrelevant to mechanical performance. Second, the results may highlight that the specific geometric configuration of the bonding network — such as bond length, cross-sectional thickness, and spatial orientation — provides critical information beyond a purely topological perspective. The challenge for deeper analysis lies in identifying and quantifying these precise local features, especially given the tension between viewing snow as a granular structure held together by bonds versus a more realistic bicontinuous material.

Some literature suggests that **pore characteristics** play a key role in explaining mechanical behavior. For instance, Roberts and Garboci [104] concluded that at high volume fractions, the shape of pores influences the elastic modulus of the material. Initially, our analysis did not seem to support this, as general pore-specific descriptors (such as SSA$_{pore}$, $\bar{I}_{pore}$, and $I_{c,pore}$) failed to rank among our top predictors. We reconcile this apparent discrepancy by proposing that the controlling factor is not pore shape per se, but rather the specific structural features that control load distribution and stress concentration. Indeed, the top predictors identified by our analysis were representations of these mechanically important aspects of the microstructure.

To further explore the grain type-dependent relationship between mechanics and microstructure, we repeated our analysis for each category individually (Fig. A.11). A particularly interesting insight from this is the variation in the number of correlated feature clusters, which we interpret as an indicator of structural similarity across different volume fractions; a low number of clusters suggests a common “building principle” for a given grain type. This is best illustrated by comparing two examples. For FC&DH samples, many microstructural parameters were highly correlated and formed one large cluster. This indicates that temperature-gradient-driven metamorphism causes structural similarity over a wide range of densities. In contrast, the correlations for SH were relatively weak, resulting in many loosely connected clusters. This could imply that SH crystals lack a consistent building principle, or that this result is an artifact of the measurement scale, with weak correlations caused by either data noise (e.g., from fine particles) or an unrepresentative sample volume. While a detailed discussion of all per-type findings of this additional analysis is beyond the scope of this paper, these results reveal a fundamental link between the genesis of a snow type and its final microstructural characteristics.

In summary, our findings show that stiffness is critically controlled by the shape of the load-bearing paths, while compressive strength is additionally influenced by the local geometry and topology of the bonding system. To better understand the role of the bonding system, future work should move beyond indirect bulk parameters and focus on a detailed analysis of deformation mechanisms and mechanically relevant length scales at the microstructural level.
### 4.3. Implications for analytical and numerical models

The high power law exponents found in our data (up to 5.7) for volume fraction scaling make it clear that simple **analytical models** based on bulk volume fraction are insufficient. This suggests that mechanical behavior is governed by a suite of complex factors — such as microstructural complexity, the effective load-bearing structure, and small-scale deformation modes. This view is supported by recent findings showing that scaling exponents are directly influenced by factors like internal friction and local geometric effects, which may alter load paths and deformation mechanisms [e.g., 39,58].

We argue that bulk volume fraction is merely a proxy for a cascade of changes in more fundamental properties like topology, local geometry, and stress concentrations, which are the true drivers of mechanical response. This issue is exacerbated in snow science by the practical difficulty of sampling thin, heterogeneous weak layers, making bulk volume fraction a particularly problematic parameter. To better isolate the true drivers, we turn to a volume fraction-independent analysis, focusing on the mechanical parameters as proposed by Köchle and Schneebeli [17].

Directly analyzing the stiffness-strength relationship (Fig. 8) allowed us to identify the dominant deformation mechanisms, which were found to govern the slopes found in the plot, thereby bypassing the need to consider volume fraction. A slope of 0.5 suggests failure driven by crack propagation [105], whereas in porous materials and cellular solids, values distinguish between bending-dominated (0.75) and stretching-dominated (1.0) structures [34,35]. Exponents greater than 1.0 are attributed to stability-driven failure, such as the elastic buckling of slender elements [106,107]. In this framework, the scaling exponent is also directly linked to the topology of the material. Specifically, a higher degree of connectivity in the microstructure promotes a shift from bending-dominated to stretching-dominated behavior [35].

Based on this theoretical framework, our results indicate two distinct mechanical regimes: DF&RG are bending-dominated (slope ≈ 0.6), while persistent grain types (FC&DH and SH) are stretching-dominated (slope ≈ 1.1–1.2). While it is currently unclear if this distinction reflects a true microstructural property or a macroscopic proxy, a Finite Element analysis would be the necessary next step to validate this framework at the microscopic level. This finding is practically significant. The distinct, grain-type-dependent power laws shown in Fig. 8 suggest that the underlying microstructural features per grain type coevolve in a predictable way. This may allow us to bypass the use of volume fraction as a proxy for mechanical properties altogether. If one property — such as strength, which is often easier to measure in the field [108] — is known, the other could be inferred, provided the grain type is also identified.

From a microstructural standpoint, the stretching-dominated behavior of persistent grain types implies they possess a higher degree of connectivity [35] (i.e., a higher average coordination number), a conclusion our topological data supports (Fig. 10). This finding is counter-intuitive, as these weak layers are conventionally described as “poorly bonded”. This apparent paradox can be resolved by the frame of reference: while persistent grains may have fewer connections per unit volume, their larger structural length scale means each grain still has a higher average number of connections than the smaller particles in DF&RG. Should this prove to be true, the fragile nature of persistent weak layers does not stem from their topology. Instead, their role in slab avalanche release would then likely originate from other factors, such as a higher susceptibility to shear loading, a lower fracture toughness, or their characteristic slow rate of densification within the snowpack.

Comparing these findings to recent **numerical studies** on snow shows a coherent picture. While the relation between stiffness and strength has not been a focus, we found that the volume fraction power laws for stiffness and strength of [56,58] result in similar scalings in the stiffness-strength plot (Fig. 8).

Gaume et al. [56] used DEM simulations based on sticky hard-sphere packings and showed that stiffness and strength scale with cohesive contact density, a product of coordination number and volume fraction. Their exponent ratio for stiffness and strength (~0.6) closely matches our findings for DF&RG. Interestingly, results based on different coordination numbers still show the same exponent ratio, thus effectively canceling out the effect of topology in the theoretical framework of [35].

Blatny et al. [58] employed Gaussian random fields (GRF) combined with the Material Point Method (MPM) to simulate porous materials and observed high scaling exponents for both strength and stiffness. Notably, the exponent ratio between stiffness and strength was around 1.2, closely matching the ratio we found for persistent grain types in our samples. Here, identifying the role of topology is not trivial, because GRF models vary volume fraction and topology together. However, [109] showed that GRF morphologies exhibit a nearly linear increase in scaled genus density with volume fraction. Therefore topology should play a role in the resulting higher exponents compared to Gaume et al. [56].

### 4.4. Limitations of this work

Our study has several limitations that should be considered when interpreting the results. First, part of our dataset is based on artificially grown weak layers. Although previous work has shown that some of the microstructural parameters are similar to the parameters of natural samples [44], it remains uncertain whether the volume fraction dependence we observed fully captures the behavior of weak layers formed under different natural conditions. Especially, the high temperature gradients used to produce a sufficient amount of samples in a given time might hinder a direct applicability to natural snow.

Second, our experimental results are subject to typical uncertainties inherent in snow science, including sensitivity to experimental conditions such as temperature, loading rate, sample preparation, and handling.

A third limitation lies in our methodical approach. It presumes that the analyzed µCT scans represent the sample's bulk mechanical properties, despite some expected spatial variability. By focusing on bulk parameters, our method assumes a homogenized microstructure is sufficient for predicting mechanical behavior, potentially masking the effects of local features such as internal layering. At the same time, the homogenized properties may not be predictive of mechanical behavior. A notable example is the occurrence of small rounded grains embedded within large surface hoar crystals, which can distort µCT-based measurements and reduce their interpretability.

Fourth, the concept of a Representative Volume Element (RVE) is difficult to apply to thin weak layers. For layers composed of large crystals such as surface hoar, which can span the entire sample height, it is questionable whether the selected volume for our analysis can be truly representative. This conceptual challenge was compounded by practical limitations, as we had to balance the competing demands of sample holder diameter, image resolution, and a sufficiently large REV.

Finally, the data-driven analysis has several limitations. Our dataset is relatively small by machine learning standards and contains considerable scatter. These factors, combined with the high correlation between input features, make the analysis challenging. Although our results appear robust, the identified correlations between microstructure and mechanics should be viewed with caution, as they do not necessarily imply causation or reflect the actual underlying mechanisms.

## 5. Conclusions and outlook

This study presents a comprehensive dataset linking the mechanical behavior of weak snow layers with detailed microstructural descriptors derived from µCT imaging. By combining experimental uniaxial compression tests, digital image correlation, and extensive morphological
analysis, we explored the microstructure and mechanisms that govern the compressive strength and stiffness of snow weak layers. Our main findings are:

*   **Distinct volume fraction scaling per grain type:** Different grain types exhibited unique power law relationships between volume fraction and mechanical properties. Faceted crystals and depth hoar showed the largest power law exponents (5.1 for the elastic modulus and 5.7 for compressive strength), followed by decomposing particles and rounded grains (4.6 and 2.8, respectively), while buried surface hoar exhibited a distinctly weaker volume fraction dependence (1.7 and 2.1, respectively).

*   **Different microstructural drivers for stiffness and strength:** Our neural network identified distinct sets of microstructural parameters governing the mechanical response. Stiffness was primarily influenced by features describing the alignment and efficiency of load-bearing paths, such as anisotropy and tortuosity. In contrast, compressive strength was more sensitive to the local geometry of the ice-air interface, though it also depended on the load path characteristics.

*   **Limited predictive power of topology metrics:** Classical topological metrics, such as the Euler characteristic and scaled genus density, did not rank among the most important features in our analysis. This suggests that for snow, either these specific descriptors are inadequate or that other factors — such as the geometry of the bonding network — play a more dominant role.

*   **Volume fraction-independent scaling regimes:** By plotting strength versus stiffness, we identified a volume fraction-independent relationship that reveals two distinct mechanical regimes. Persistent grain types (FC&DH and SH) exhibited a steep scaling slope (~1.1), characteristic of stretching-dominated structures. In contrast, decomposing particles and rounded grains (DF&RG) showed a shallower slope (~0.65), characteristic of bending-dominated structures. In the context of analytical models such as [34,35] this implies that persistent grains have a higher average coordination number — a counter-intuitive implication for persistent weak layers.

While our study provides valuable insights, several open challenges remain in linking snow microstructure to its mechanical properties. The primary blind spot remains the influence of combined compressive and shear loading, as weak layer failure attributed to slab avalanche release happens under multiaxial loading conditions. Future experiments are needed to complement the picture.

Clarifying the precise role of topology in snow mechanics remains a key challenge. Future work should explore other topological descriptors and their appropriate length scales. Investigating whether certain snow types exhibit structural similarity, for instance, could simplify structure-function relationships and guide the development of synthetic microstructure models. Advanced imaging and machine learning methods offer a path forward, not only to better quantify grain-based metrics such as coordination number, but also to test the validity of the grain-centric model itself. This is crucial, as viewing snow as a bicontinuous medium would render such discrete descriptors physically questionable.

One key question yet to be answered is the exact origin of the high volume fraction scaling exponents. Numerical simulations and Finite

Element studies could help to explore the mechanisms at work at the microstructure level.

A complete mechanical framework for snow must also address its behavior at the volume fraction extremes. While we have investigated a broad range of densities, it remains uncertain over which volume fraction ranges the existing formulations are valid, as they fail to converge to the known mechanical properties of polycrystalline ice as volume fraction approaches unity. A general formulation must also account for the material's percolation limit, which likely differs between grain types. Certain weak layer types, such as buried surface hoar, can persist at very low densities, suggesting a lower percolation threshold that must be considered to develop a universally valid model.

Finally, to facilitate cross-study comparisons and model validation in snow science, we advocate for the standardization of morphological descriptors, using a consistent set of physically meaningful and computable parameters such as a selected set of those listed in Table 1.

### CRediT authorship contribution statement

**Jakob Schöttner:** Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Supervision, Validation, Visualization, Writing – original draft. **Berit Zeller-Plumhoff:** Conceptualization, Formal analysis, Methodology, Software, Writing – original draft. **Pascal Hagenmuller:** Data curation, Investigation, Software, Writing – original draft. **Philipp Weißgraeber:** Conceptualization, Writing – original draft. **Philipp L. Rosendahl:** Conceptualization, Funding acquisition, Writing – original draft. **Henning Löwe:** Conceptualization, Methodology, Resources, Software. **Jürg Schweizer:** Conceptualization, Resources, Supervision, Writing – original draft. **Alec van Herwijnen:** Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – original draft.

### 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.

### Acknowledgments

We thank M. Walet, V. Adam, and S. Kraus for their assistance with sample preparation, data analysis, and for many helpful discussions. We are grateful to M. Jaggi and Lucid Concepts AG for their support in preparing and evaluating the µCT scans. We also extend our thanks to J. Gaume, L. Blatny, and M. Hedvard for insightful discussions on the role of morphology in analytical and numerical models. Special thanks are extended to W. L. Roque for providing the geodesic reconstruction plugin for ImageJ, which was instrumental in calculating tortuosity; to T. Theile for testing his grain segmentation algorithm on some of our µCT scans; and to M. Hach for his assistance with the data for Fig. 2.

We used deepl.com/write and gemini.google.com to enhance the language and readability of this manuscript.

This research has been supported by the Swiss National Science Foundation (SNF) under Grant No. 201071 and German Research Foundation (DFG) under Grant No. 460195514.

## Appendix A

See Figs. A.1–A.12.
<figure><img src="image_17.png" /><figcaption>**Fig. A.1.** Evaluation of the elastic modulus based solely on the slope of the stress–strain curve leading up to the ultimate strength. The vertical error bars represent the standard deviation of the experimental results, the horizontal error bars are an estimate for the spatial variability of volume fraction within a parent sample [44].</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. A.2.** Specific surface area (SSA) vs. volume fraction $\phi$ for our weak layer samples. (a) shows the normalization by ice volume SSA$_{ice}$, (b) the normalization by pore volume SSA$_{pore}$. The different colors correspond to the different weak layer types. The horizontal error bars are an estimate for the spatial variability of volume fraction within a parent sample [44]. We observed a decreasing trend of SSA$_{ice}$ with increasing volume fraction for all types of weak layers, except for the two high-density FC&DH samples, where SSA$_{ice}$ increased again. In contrast, SSA$_{pore}$ shows an increasing trend with volume fraction for all weak layer types. Notably, DF&RG exhibited higher values and greater scatter than the other types.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. A.3.** Correlation lengths in different directions and the resulting anisotropy. The correlation lengths in the horizontal (x and y) directions are nearly identical, so we show their average ($\bar{\epsilon}^{xy}$), whereas those in the vertical ($\bar{\epsilon}^z$) direction show different values but follow similar trends. The horizontal error bars are an estimate for the spatial variability of volume fraction within a parent sample [44]. Overall, DF&RG exhibited the lowest correlation lengths, with a slight increase as volume fraction rises, while FC&DH displays the opposite trend. For SH, the horizontal correlation lengths increased with volume fraction, but the vertical lengths remained relatively constant. Anisotropy was quantified by the ratio of the z-direction correlation length to the average horizontal correlation length, as defined in Eq. (6). The plot shows clear differences in anisotropy among the weak layer types. DF&RG typically had values near or below 1, indicating a predominance of horizontal structures, with a slight increase as volume fraction increases. In contrast, FC&DH shows strong vertical anisotropy, with most values exceeding 1.2 and increasing with volume fraction. SH, on the other hand, exhibited values both above and below 1 and showed a clear decrease with increasing volume fraction.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. A.4.** Example comparison of the principal curvatures of the three different weak layer types: (a) faceted crystals and depth hoar, (b) decomposing particles and rounded grains, (c) buried surface hoar. We normalized the curvatures with the inverse SSA$_{ice}$ to allow comparison with other cellular materials. Both FC&DH and SH exhibited similar distributions, with many points clustering near $\kappa_1 = \kappa_2 = 0$, indicating planes and areas with small curvatures. In contrast, FC&RG displayed a more uniform distribution—especially in the region $\kappa_1 > 0, \kappa_2 < 0$, which corresponds to saddle shaped structures. In all three samples, distinct lines appear in the plots, representing points where the relationship between the principal curvatures remains constant. This indicates the presence of size-invariant microstructural features.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. A.5.** Classification of the snow-air interface into (a) concave, (b) convex and (c) saddle shaped regions based on the principal curvatures. The horizontal error bars are an estimate for the spatial variability of volume fraction within a parent sample [44]. Both FC&DH and SH exhibited similar distributions across all three plots. The concave and saddle fractions increased linearly with volume fraction, while the convex fraction decreased. DF&RG showed similar trends, yet displayed a clear offset relative to the other types of weak layers for the concave and convex fractions. In contrast, the saddle fraction for DF&RG exhibited considerably more scatter than observed in the other two types.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. A.6.** Local thickness defined as the maximum-sized sphere that can be inscribed within a structure (ice or pore) analysis for weak layer structures (the error bars represent the standard deviation). (a) shows the mean ice thickness versus volume fraction, (b) shows the mean pore thickness versus volume fraction, and (c) shows the mean ice thickness plotted against mean pore thickness. Although the mean ice thickness (Panel a) for all structures lies between 0.1 and 0.4 mm, an increasing trend with volume fraction is observed for DF&RG. For DF&RG, the mean pore thickness (Panel b) remained constant, whereas a decreasing trend was apparent for FC&DH and SH. Plotting the mean ice thickness against the mean pore thickness (Panel c) revealed clear separations among the different weak layer types.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. A.7. Tortuosity obtained via geodesic reconstruction considering 6 neighboring voxels, averaged for the two horizontal directions (x and y - Panel a) and the vertical direction (z - Panel b). The colors correspond to the different weak layers. The horizontal error bars are an estimate for the spatial variability of volume fraction within a parent sample [44]. Tortuosity decreased with increasing volume fraction in all three directions. Notably, DF&RG remained nearly constant in the horizontal directions, while it showed a distinct decrease in the vertical direction, similar to FC&DH. Moreover, the absolute values for DF&RG were approximately 1.1 times larger than those for FC&DH. SH exhibited more scatter, likely because only a few grains determine the outcome of the algorithm.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. A.8. Minimum cut density for different types of weak layers, represented by various colors, with volume fraction. The horizontal error bars are an estimate for the spatial variability of volume fraction within a parent sample [44]. Panel (a) shows the percentage of the minimum cut area in the vertical direction whose removal separates the sample into two parts. The scaling of this parameter with volume fraction was best described by power laws for all three types of weak layers. Specifically, FC&DH exhibited the strongest scaling with an exponent of 3.4, DF&RG follows with an exponent of 3.0, and SH showed an exponent of 1.6. Panel (b) displays the corresponding bond density (number of bonds per mm²). For FC&DH and SH, the bond density followed a behavior similar to that of the minimum cut area, whereas for DF&RG no clear trend was observed; instead, the bond density remained relatively constant with increasing volume fraction and exhibited strong scatter. Panel (c) presents the anisotropy, calculated from the minimum cut areas in different directions according to Eq. (6). It followed similar trends to those based on correlation lengths (Fig. A.3), albeit with generally higher overall values.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Fig. A.9.** Mean elastic modulus vs. mean strain rate for the last four pop-ins and the final slope leading to failure for all tested samples. To enable comparison between samples, values were normalized by each sample's average elastic modulus and strain rate. The error bars represent the standard deviations. The plot shows that, on average, the Young's modulus remained constant across pop-ins, while the strain rate increased with each successive pop-in.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. A.10.** Strain at failure for all our samples. The mean value is 0.21%, the error bars show the experimental standard deviation.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Fig. A.11.** Feature clustering and feature importance analyzed for each grain type category individually (FC&DH - Panel (a), DF&RG - Panel (b), and SH - Panel (c)).</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Fig. A.12. Principal component analysis (PCA) of our microstructural dataset. Panel (a) shows that only six principal components are needed to explain 95% of the variance. The loading matrix (Panel b) revealed that no single parameter dominates; instead, most parameters contributed similarly. To quantify feature importance, we calculated a score defined as the sum of loadings normalized by the proportion of variance explained by each principal component (Panel c). Here we see that most parameters have a relatively high feature importance score, with SSA leading the top of the list.</figcaption></figure>
## Data availability

The dataset is available on EnviDat: https://www.doi.org/10.1690/4/eenvidat.648.

## References

[1] J. Schweizer, B. Jamieson, M. Schneebeli, Snow avalanche formation, Rev. Geophys. 41 (4) (2003) http://dx.doi.org/10.1029/2002RG000123.

[2] A. van Herwijnen, B. Jamieson, Snowpack properties associated with fracture initiation and propagation resulting in skier-triggered dry snow slab avalanches, Cold Reg. Sci. & Technol. 50 (1-3) (2007) 13-22, http://dx.doi.org/10.1016/j.coldregions.2007.02.004.

[3] J. Schweizer, J.B. Jamieson, Snowpack properties for snow profile analysis, Cold Reg. Sci. & Technol. 37 (3) (2003) 233-241, http://dx.doi.org/10.1016/S0165-232X(03)00067-3.

[4] J. Schweizer, J. Jamieson, Snow cover properties for skier triggering of avalanches, Cold Reg. Sci. & Technol. 33 (2-3) (2001) 207-221, http://dx.doi.org/10.1016/S0165-232X(01)00039-8.

[5] O. Abe, Shear strength and angle of repose of snow layers including graupel, Ann. Glaciol. 38 (2004) 305-308, http://dx.doi.org/10.3189/172756404781815149.

[6] K. Fourteau, J. Freitag, M. Malinen, H. Löwe, Microstructure-based simulations of the viscous densification of snow and firn, Cryosphere 18 (6) (2024) 2831-2846, http://dx.doi.org/10.5194/tc-18-2831-2024.

[7] N. Calonne, B. Richter, H. Löwe, C. Cetti, J. ter Schure, A. van Herwijnen, C. Fierz, M. Jaggi, M. Schneebeli, The RHOSSA campaign: multi-resolution monitoring of the seasonal evolution of the structure and mechanical stability of an alpine snowpack, Cryosphere 14 (6) (2020) 1829-1848, http://dx.doi.org/10.5194/tc-14-1829-2020.

[8] V. Adam, B. Bergfeld, P. Weißgräber, A. van Herwijnen, P.L. Rosendahl, Nat. Fracture toughness of mixed-mode anticracks in highly porous materials, Nat. Commun. 15 (1) (2024) 7379, http://dx.doi.org/10.1038/s41467-024-51491-7.

[9] P.L. Rosendahl, P. Weißgräber, Modeling snow slab avalanches caused by weak-layer failure - Part 1: Slabs on compliant and collapsible weak layers, Cryosphere 14 (1) (2020) 115-130, http://dx.doi.org/10.5194/tc-14-115-2020.

[10] H. Narita, Mechanical behaviour and structure of snow under uniaxial tensile stress, J. Glaciol. 26 (94) (1980) 275-282, http://dx.doi.org/10.3189/s0022143000010819.

[11] J. Schweizer, Laboratory experiments on shear failure of snow, Ann. Glaciol. 26 (1998) 97-102, http://dx.doi.org/10.3189/1998AoG26-1-97-102.

[12] A. van Herwijnen, D.A. Miller, Experimental and numerical investigation of the sintering rate of snow, J. Glaciol. 59 (214) (2013) 269-274, http://dx.doi.org/10.3189/2013JoG12J094.

[13] F. Domine, J. Bock, S. Morin, G. Giraud, Linking the effective thermal conductivity of snow to its shear strength and density, J. Geophys. Res. 116 (F4) (2011) http://dx.doi.org/10.1029/2011JF002000.

[14] B. Jamieson, C.D. Johnston, Evaluation of the shear frame test for weak snowpack layers, Ann. Glaciol. 32 (2001) 59-69, http://dx.doi.org/10.3189/1727564017819472.

[15] L. Shapiro, J. Johnson, M. Sturm, G. Blaisdell, Snow Mechanics: Review of the State of Knowledge and Applications, Technical Report, US Army Corps of Engineers Cold Regions Research and Engineering Laborator, 1997, p. 40.

[16] M. Mellor, A review of basic snow mechanics, in: Symposium at Grindelwald 1974 - Snow Mechanics, IAHS Publ., 114, 1975, pp. 251-291.

[17] B. Köchle, M. Schneebeli, Three-dimensional microstructure and numerical calculation of elastic properties of alpine snow with a focus on weak layers, J. Glaciol. 60 (222) (2014) 705-713, http://dx.doi.org/10.3189/2014JoG13J220.

[18] J.J. Petrovic, Mechanical properties of ice and snow, J. Mater. Sci. 38 (1) (2003) 1-6, http://dx.doi.org/10.1023/A:1021134128038.

[19] K. Sundu, J. Freitag, K. Fourteau, H. Löwe, A microstructure-based parameterization of the effective anisotropic elasticity tensor of snow, firn, and bubbly ice, Cryosphere 18 (4) (2024) 1579-1596, http://dx.doi.org/10.5194/tc-18-1579-2024.

[20] B. Gerling, H. Löwe, A. van Herwijnen, Measuring the elastic modulus of snow, Geophys. Res. Lett. 44 (21) (2017) http://dx.doi.org/10.1002/2017GL075110.

[21] H.O.K. Kirchner, H. Peterlik, G. Michot, Size independence of the strength of snow, Phys. Rev. E 69 (1) (2004) 011306, http://dx.doi.org/10.1103/PhysRevE.69.011306.

[22] J. Schweizer, G. Michot, H.O. Kirchner, On the fracture toughness of snow, Ann. Glaciol. 38 (2004) 1-8, http://dx.doi.org/10.3189/172756404781814906.

[23] J.L. Smith, The elastic constants, strength and density of Greenland snow as determined from measurements of sonic wave velocity, in: US Army Cold Regions Research & Engineering Laboratory, vol. 167, US Army Cold Regions Research & Engineering Laboratory, 1965.

[24] S. Kinosita, Compression of snow at constant speed, in: H. Oura (Ed.), Physics of Snow and Ice - Proceedings International Conference on Low Temperature Science, Sapporo, Japan, 14-19 August 1996, Hokkaido University, Sapporo, Japan, 1967, pp. 911-927.

[25] C. Scapozza, Entwicklung eines dichte- und temperaturabhängigen stoffgesetzes zur beschreibung des visko-elastischen verhaltens von schnee (Ph.D. thesis), Eidgenössische Technische Hochschule Zürich, Zürich, 2004.

[26] A. Capelli, J.C. Kapil, I. Reiweger, D. Or, J. Schweizer, Speed and attenuation of acoustic waves in snow: Laboratory experiments and modeling with Biot's theory, Cold Reg. Sci. & Technol. 125 (2016) 1-11, http://dx.doi.org/10.1016/j.coldregions.2016.01.004.

[27] Y. Zhang, Z. Qian, S. Lv, W. Huang, J. Ren, Z. Fang, X. Chen, Experimental investigation of uniaxial compressive strength of distilled water ice at different growth temperatures, Water 14 (24) (2022) 4079, http://dx.doi.org/10.3390/w14244079.

[28] L.W. Gold, Engineering properties of fresh-water ice, J. Glaciol. 19 (81) (1977) 197-212, http://dx.doi.org/10.3189/S0022143000215608.

[29] C. Ma, T. Du, X. Niu, Y. Fan, Biomechanics and mechanobiology of the bone matrix, Bone Res. 10 (1) (2022) 59, http://dx.doi.org/10.1038/s41413-022-00223-y.

[30] E.F. Morgan, G.U. Unnikrisnan, A.I. Hussein, Bone mechanical properties in healthy and diseased states, Annu. Rev. Biomed. Eng. 20 (1) (2018) 119-143, http://dx.doi.org/10.1146/annurev-bioeng-062117-121139.

[31] J. Li, J. Li, Y. Chen, J. Chen, Strengthening modulus and softening strength of nanoporous gold in multiaxial tension: insights from molecular dynamics, Nanomaterials 12 (24) (2022) 4381, http://dx.doi.org/10.3390/nano12244381.

[32] L-Z. Liu, H.-J. Jin, Scaling equation for the elastic modulus of nanoporous gold with "fixed" network connectivity, Appl. Phys. Lett. 110 (21) (2017) http://dx.doi.org/10.1063/1.4984108.

[33] T.J. Balk, C. Eberl, Y. Sun, K.J. Hemker, D.S. Gianola, Tensile and compressive microspecimen testing of bulk nanoporous gold, JOM 61 (12) (2009) 26-31, http://dx.doi.org/10.1007/s11837-009-0176-6.

[34] L.J. Gibson, M.F. Ashby, Cellular Solids, Cambridge University Press, 1997, http://dx.doi.org/10.1017/CBO9781139878326.

[35] M. Ashby, The properties of foams and lattices, Philos. Trans. R. Soc. A: Math., Phys. Eng. Sci. 364 (1838) (2006) 15-30, http://dx.doi.org/10.1098/rsta.2005.1678.

[36] G. Bruno, A.M. Efremov, A.N. Levandovskyi, B. Clausen, Connecting the macro- and microstrain responses in technical porous ceramics: modeling and experimental validations, J. Mater. Sci. 46 (1) (2011) 161-173, http://dx.doi.org/10.1007/s10853-010-4899-0.

[37] L-Z. Liu, X.-L. Ye, H.-J. Jin, Interpreting anomalous low-strength and low-stiffness of nanoporous gold: Quantification of network connectivity, Acta Mater. 118 (2016) 77-87, http://dx.doi.org/10.1016/j.actamat.2016.07.033.

[38] H. Zhong, T. Song, C. Li, R. Das, J. Gu, M. Qian, The Gibson-Ashby model for additively manufactured metal lattice materials: Its theoretical basis, limitations and new insights from remedies, Curr. Opin. Solid State Mater. Sci. 27 (3) (2023) 101081, http://dx.doi.org/10.1016/j.cossms.2023.101081.

[39] H. Zhong, R. Das, J. Gu, M. Qian, Low-density, high-strength metal mechanical metamaterials beyond the Gibson-Ashby model, Mater. Today 68 (2023) 96-107, http://dx.doi.org/10.1016/j.matod.2023.07.018.

[40] N. Huber, Connections between topology and macroscopic mechanical properties of three-dimensional open-pore materials, Front. Mater. 5 (2018) 69, http://dx.doi.org/10.3389/fmatrs.2018.00069.

[41] C. Richert, N. Huber, A review of experimentally informed micromechanical modeling of nanoporous metals: From structural descriptors to predictive structure-property relationships, Materials 13 (15) (2020) 3307, http://dx.doi.org/10.3390/ma13153307.

[42] Q. Krol, H. Löwe, Relating optical and microwave grain metrics of snow: the relevance of grain shape, Cryosphere 10 (6) (2016) 2847-2863, http://dx.doi.org/10.5194/tc-10-2847-2016.

[43] R. Inoue, T. Aoki, S. Fujita, S. Tsutaki, H. Motoyama, F. Nakazawa, K. Kawamura, Spatial variation in the specific surface area of surface snow measured along the traverse route from the coast to Dome Fuji, Antarctica, during austral summer, Cryosphere 18 (8) (2024) 3513-3531, http://dx.doi.org/10.5194/tc-18-3513-2024.

[44] J. Schöttner, M. Walet, P. Rosendahl, P. Weissgräber, V. Adam, B. Walter, F. Rheinschmidt, H. Löwe, J. Schweizer, A. van Herwijnen, On the compressive strength of weak snow layers of depth hoar, J. Glaciol. 71 (2025) e54, http://dx.doi.org/10.1017/jog.2025.16.

[45] K.C. Agrawal, R.K. Mittal, Influence of microstructure on mechanical properties of snow, Def. Sci. J. 45 (2) (1995) 93-105.

[46] P.R. Kry, The Relationship between the Visco-Elastic and Structural Properties of fine-grained Snow, J. Glaciol. 14 (72) (1975).

[47] P. Hagenmuller, G. Chambon, F. Flin, S. Morin, M. Naaim, Snow as a granular material: assessment of a new grain segmentation algorithm, Granul. Matter 16 (4) (2014) 421-432, http://dx.doi.org/10.1007/s10035-014-0503-7.

[48] T. Theile, M. Schneebeli, Algorithm to decompose three-dimensional complex structures at the necks: tested on snow structures, IET Image Process. 5 (2) (2011) 132-140, http://dx.doi.org/10.1049/iet-ipr.2009.0410.

[49] C. Willibald, H. Löwe, T. Theile, J. Dual, M. Schneebeli, Angle of repose experiments with snow: role of grain shape and cohesion, J. Glaciol. 66 (258) (2020) 658-666, http://dx.doi.org/10.1017/jog.2020.36.
[50] P. Hagenmuller, N. Calonne, G. Chambon, F. Flin, C. Geindreau, M. Naaim, Characterization of the snow microstructural bonding system through the minimum cut density, Cold Reg. Sci. & Technol. 108 (2014) 72-79, http://dx.doi.org/10.1016/j.coldregions.2014.09.002.

[51] R.W. Rice, Evaluation and extension of physical property-porosity models based on minimum solid area, J. Mater. Sci. 31 (1) (1996) 102-118, http://dx.doi.org/10.1007/BF00355133/METRICS.

[52] R.W. Rice, Comparison of physical property-porosity behaviour with minimum solid area models, J. Mater. Sci. 31 (6) (1996) 1509-1528, http://dx.doi.org/10.1007/BF00357860.

[53] H. Löwe, G. Picard, Microwave scattering hard of snow in MEMLS and DMRT-ML revisited: the relevance of sticky hard spheres and tomography-based estimates of stickiness, Cryosphere 9 (6) (2015) 2101-2117, http://dx.doi.org/10.5194/tc-9-2101-2015.

[54] R.J. Baxter, Percus-Yevick equation for hard spheres with surface adhesion, J. Chem. Phys. 49 (6) (1968) 2770-2774, http://dx.doi.org/10.1063/1.1670482.

[55] S. Pilia, F. Baroni, A. Lapini, S. Paloscia, S. Pettinato, E. Santi, P. Pampaloni, M. Valt, F. Monti, On the relationship between stickiness in DMRT theory and physical parameters of snowpack: Theoretical formulation and experimental validation with SNOWPACK: Snow model and X-Band SAR data, IEEE Trans. Geosci. Remote Sens. 60 (2022) 1-14, http://dx.doi.org/10.1109/TGRS.2022.3201669.

[56] J. Gaume, H. Löwe, S. Tan, L. Tsang, Scaling laws for the mechanics of loose and cohesive granular materials based on Baxter's sticky hard spheres, Phys. Rev. E 96 (3) (2017) http://dx.doi.org/10.1103/PhysRevE.96.032914.

[57] J. Ritter, H. Löwe, J. Gaume, Microstructural controls of anticrack nucleation in highly porous brittle solids, Sci. Rep. 10 (1) (2020) 12383, http://dx.doi.org/10.1038/s41598-020-67926-2.

[58] L. Blatny, H. Löwe, S. Wang, J. Gaume, Computational micromechanics of porous brittle solids, Comput. Geotech. 140 (2021) 104284, http://dx.doi.org/10.1016/j.compgeo.2021.104284.

[59] T. Ozeki, M. Tsuda, Y. Yashiro, K. Fujita, S. Adachi, Development of artificial surface hoar production system using a circuit wind tunnel and formation of various crystal types, Cold Reg. Sci. & Technol. 169 (2020) http://dx.doi.org/10.1016/j.coldregions.2019.102889.

[60] D. Turner, P. Crozier, P. Reu, Digital Image Correlation Engine, Technical Report, Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA ..., National Technology & Engineering Solutions of Sandia, LLC (NTESS), 2015.

[61] H. Löwe, M. Zaiser, S. Mösinger, S. Schleef, Snow mechanics near the ductile-brittle transition: Compressive stick-slip and snow microquakes, Geophys. Res. Lett. 47 (4) (2020) http://dx.doi.org/10.1029/2019GL085491.

[62] ASTM International, ASTM E8/E8M-22: Test Methods for Tension Testing of Metallic Materials, Technical Report, ASTM International, West Conshohocken, PA, 2022, http://dx.doi.org/10.1520/E0008(_)E0008M-13A.

[63] S. Nanthikesan, S. Shyam Sunder, Anisotropic elasticity of polycrystalline ice, Cold Reg. Sci. & Technol. 22 (2) (1994) 149-169, http://dx.doi.org/10.1016/0165-232X(94)90026-4.

[64] P. Hagenmuller, G. Chambon, B. Lesaffre, F. Flin, M. Naaim, Energy-based binary segmentation of snow microtomographic images, J. Glaciol. 59 (217) (2013) 859-873, http://dx.doi.org/10.3189/2013JoG13J035.

[65] B. Walter, H. Weigel, S. Wahl, H. Löwe, Wind tunnel experiments to quantify the effect of aeolian snow transport on the surface snow microstructure, Cryosphere Discuss. 2023 (2023) 1-27, http://dx.doi.org/10.5194/tc-2023-112.

[66] T.U. Kaempfer, M. Schneebeli, S.A. Sokratov, A microstructural approach to model heat transfer in snow, Geophys. Res. Lett. 32 (21) (2005) 1-5, http://dx.doi.org/10.1029/2005GL023873.

[67] C. Coléou, B. Lesaffre, J.B. Brzoska, W. Ludwig, E. Boller, Three-dimensional snow images by X-ray microtomography, Ann. Glaciol. 32 (2001) 75-81, http://dx.doi.org/10.3189/172756401781819418.

[68] W.E. Lorensen, H.E. Cline, Marching cubes, in: Seminal Graphics, ACM, New York, NY, USA, 1998, pp. 347-353, http://dx.doi.org/10.1145/280811.281026.

[69] S. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid, J.-Y. Tinevez, D.J. White, V. Hartenstein, K. Eliceiri, P. Tomancak, A. Cardona, Fiji: an open-source platform for biological-image analysis, Nature Methods 9 (7) (2012) 676-682, http://dx.doi.org/10.1038/nmeth.2019.

[70] R. Domander, A.A. Felder, M. Doube, BoneJ2 - refactoring established research software, Wellcome Open Res. 6 (2021) 1-21, http://dx.doi.org/10.12688/WELLCOMEOPENRES.16619.2.

[71] M. Doube, M.M. Kłosowski, I. Arganda-Carreras, F.P. Cordelières, R.P. Dougherty, J.S. Jackson, B. Schmid, J.R. Hutchinson, S.J. Shefeline, Bone: Free and extensible bone image analysis in ImageJ, Bone 47 (6) (2010) 1076-1079, http://dx.doi.org/10.1016/j.bone.2010.08.023.

[72] W.L. Roque, R.R. Costa, A plugin for computing the pore/grain network tortuosity of a porous medium from 2D/3D MicroCT image, Appl. Comput. Geosci. 5 (2020) 100019, http://dx.doi.org/10.1016/j.acags.2020.100019.

[73] C.A. Schneider, W.S. Rasband, K.W. Eliceiri, NIH Image to ImageJ: 25 years of image analysis, Nature Methods 9 (7) (2012) 671-675, http://dx.doi.org/10.1038/nmeth.2089.

[74] C.J. Gommes, A.J. Bons, S. Blacher, J.H. Dunsmuir, A.H. Tsou, Practical methods for measuring the tortuosity of porous materials from binary or gray-tone tomographic reconstructions, AIChE J. 55 (8) (2009) 2000-2012, http://dx.doi.org/10.1002/AIC.11812.

[75] C. Sullivan, A. Kaszynski, PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK), J. Open Source Softw. 4 (37) (2019) 1450, http://dx.doi.org/10.21105/joss.01450.

[76] M.Q. Edens, R.L. Brown, Changes in microstructure of snow under large deformations, J. Glaciol. 37 (126) (1991) 193-202, http://dx.doi.org/10.3189/S0022143000007206.

[77] K. Michelsen, H. De Raedt, J. De Hosson, Aspects of mathematical morphology, in: Advances in Imaging and Electron Physics, vol. 125, 2003, pp. 119-194, http://dx.doi.org/10.1016/S1076-5670(02)80016-7.

[78] K. Hu, M. Ziehmer, K. Wang, E.T. Lilleodden, Nanoporous gold: 3D structural analyses of representative volumes and their implications on scaling relations of mechanical behaviour, Phil. Mag. 96 (32-34) (2016) 3322-3335, http://dx.doi.org/10.1080/14786435.2016.1222087.

[79] R. Mendoza, K. Thornton, I. Savin, P. Voorhees, The evolution of interfacial topology during coarsening, Acta Mater. 54 (3) (2006) 743-750, http://dx.doi.org/10.1016/j.actamat.2005.10.010.

[80] R.T. DeHoff, E.H. Aigeltinger, K.R. Craig, Experimental determination of the topological properties of three-dimensional microstructures, J. Microsc. 95 (1) (1972) 69-91, http://dx.doi.org/10.1111/j.1365-2818.1972.tb03712.x.

[81] Y. Kwon, K. Thornton, P. Voorhees, Morphology and topology in coarsening of domains via non-conserved and conserved dynamics, Phil. Mag. 90 (2010) 317-335, http://dx.doi.org/10.1080/1478643093260701.

[82] M. Ziehmer, K. Hu, K. Wang, E.T. Lilleodden, A principle curvatures analysis of the isothermal evolution of nanoporous gold: Quantifying the characteristic length-scales, Acta Mater. 120 (2016) 24-31, http://dx.doi.org/10.1016/j.actamat.2016.08.028.

[83] N. Calonne, F. Flin, C. Geindreau, B. Lesaffre, S. Rolland du Roscoat, Study of a temperature gradient metamorphism of snow from 3-D images: time evolution of microstructures, physical properties and their associated anisotropy, Cryosphere 8 (6) (2014) 2255-2274, http://dx.doi.org/10.5194/tc-8-2255-2014.

[84] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournaapeau, M. Brucher, M. Perrot, E. Duchesnay, Scikit-learn: machine learning in Python, J. Mach. Learn. Res. 12 (2011) 2825-2830.

[85] P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournaapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S.J. van der Walt, M. Brett, J. Wilson, K.J. Millman, N. Mayorov, A.R.J. Nelson, E. Jones, R. Kern, E. Larson, C.J. Carey, I. Polat, Y. Feng, E.W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E.A. Quintero, C.R. Harris, A.M. Archibald, A.H. Ribeiro, F. Pedregosa, P. van Mulbregt, A. Vijaykumar, A.P. Bardelli, A. Rothberg, A. Hilboll, A. Kloeckner, A. Scopatz, A. Lee, A. Rokem, C.N. Woods, C. Fulton, C. Masson, C. Haeggstroem, C. Fitzgerald, D.A. Nicholson, D.R. Hagen, D.V. Pasechnik, E. Olivetti, E. Martin, E. Wieser, F. Silva, F. Lenders, F. Wilhelm, G. Young, G.A. Price, G.-L. Ingold, G.E. Allen, G.R. Lee, H. Audren, I. Probst, J.P. Dietrich, J. Silterra, J.T. Webber, J. Slavic, J. Nothman, J. Buchner, J. Kulick, J.L. Schoenberger, J.V. de Miranda Cardoso, J. Reimer, J. Harrington, J.L.C. Rodriguez, J. Nunez-Iglesias, J. Kuczyński, K. Tritz, M. Thoma, M. Newville, M. Kuemmerer, M. Bolingbroke, M. Tartre, M. Pak, N.J. Smith, N. Nowaczyk, N. Shebanov, O. Pavlyk, P.A. Brodkorb, P. Lee, R.T. McGibbon, F. Feldbauer, S. Lewis, S. Tygier, S. Sievert, S. Vigna, S. Peterson, S. More, T. Pudlik, T. Oshima, T.J. Pingel, T.P. Robitaille, T. Spura, T.R. Jones, T. Cera, T. Leslie, T. Zito, T. Krauss, U. Upadhyay, Y.O. Halchenko, Y. Vazquez-Baeza, SciPy 1.0: fundamental algorithms for scientific computing in Python, Nature Methods 17 (3) (2020) 261-272, http://dx.doi.org/10.1038/s41592-019-0686-2.

[86] C.M. Bishop, Neural Networks for Pattern Recognition, Oxford University Press, 1995.

[87] I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, 2016.

[88] M. Kuhn, K. Johnson, Applied Predictive Modeling, Springer New York, New York, NY, 2013, http://dx.doi.org/10.1007/978-1-4614-6849-3.

[89] A. de Myttenaere, B. Golden, B. Le Grand, F. Rossi, Mean absolute percentage error for regression models, Neurocomputing 192 (2016) 38-48, http://dx.doi.org/10.1016/j.NEUCOM.2015.12.114.

[90] F. Murtagh, P. Legendre, Ward's hierarchical agglomerative clustering method: Which algorithms implement ward's criterion? J. Classification 31 (3) (2014) 274-295, http://dx.doi.org/10.1007/S00357-014-9161-9/METRICS.

[91] L. Breiman, Random forests, Mach. Learn. 45 (1) (2001) 5-32, http://dx.doi.org/10.1023/A:1010933404324/METRICS.

[92] P.T. Boggs, J.E. Rogers, Orthogonal distance regression, Contemp. Math. 112 (1990) 183-194.

[93] Y.-c.K. Chen, Y.S. Chu, J. Yi, I. McNulty, Q. Shen, P.W. Voorhees, D.C. Dunand, Morphological and topological analysis of coarsened nanoporous gold by x-ray nanotomography, Appl. Phys. Lett. 96 (4) (2010) 043122, http://dx.doi.org/10.1063/1.3285175.
[94] A. van Herwijnen, J. Gaume, E.H. Bair, B. Reuter, K.W. Birkeland, J. Schweizer, Estimating the effective elastic modulus and specific fracture energy of snow-pack layers from field experiments, J. Glaciol. 62 (236) (2016) 997–1007, http://dx.doi.org/10.1017/jog.2016.90.

[95] I. Reiweger, J. Schweizer, Weak layer fracture: facets and depth hoar, Cryosphere 7 (5) (2013) 1447–1453, http://dx.doi.org/10.5194/tc-7-1447-2013.

[96] C. Chandel, P.K. Srivastava, P. Mahajan, Determination of failure envelope for faceted snow through numerical simulations, Cold Reg. Sci. & Technol. 116 (2015) 56–64, http://dx.doi.org/10.1016/j.coldregions.2015.04.009.

[97] T. Mede, G. Chambon, P. Hagenmuller, F. Nicot, Snow failure modes under mixed loading, Geophys. Res. Lett. 45 (24) (2018) 13351–13358, http://dx.doi.org/10.1029/2018GL080637.

[98] E. Wang, X. Fu, H. Han, X. Liu, Y. Xiao, Y. Leng, Study on the mechanical properties of compacted snow under uniaxial compression and analysis of influencing factors, Cold Reg. Sci. & Technol. 182 (2021) 103215, http://dx.doi.org/10.1016/j.coldregions.2020.103215.

[99] C. Chandel, P.K. Srivastava, V. Kumar, P. Datt, R. Sheoran, P. Satayawali, Laboratory set-up for surface hoar layer growth over rounded grain snow, Cold Reg. Sci. & Technol. 205 (2023) 103705, http://dx.doi.org/10.1016/j.coldregions.2022.103705.

[100] I. Reiweger, J. Gaume, J. Schweizer, A new mixed-mode failure criterion for weak snowpack layers, Geophys. Res. Lett. 42 (5) (2015) 1427–1432, http://dx.doi.org/10.1002/2014GL062780.

[101] P. Hagenmuller, G. Chambon, M. Naaim, Microstructure-based modeling of snow mechanics: a discrete element approach, Cryosphere 9 (5) (2015) 1969–1982, http://dx.doi.org/10.5194/tc-9-1969-2015.

[102] W.L. Roque, K. Arcaro, A. Alberich-Bayarri, Mechanical competence of bone: A new parameter to grade trabecular bone fragility from tortuosity and elasticity, IEEE Trans. Biomed. Eng. 60 (5) (2013) 1363–1370, http://dx.doi.org/10.1109/TBME.2012.2234457.

[103] S. Torquato, Random heterogeneous media: Microstructure and improved bounds on effective properties, Appl. Mech. Rev. 44 (2) (1991) 37–76, http://dx.doi.org/10.1115/1.3119494.

[104] A.P. Roberts, E.J. Garboczi, Elastic properties of model porous ceramics, J. Am. Ceram. Soc. 83 (12) (2000) 3041–3048, http://dx.doi.org/10.1111/j.1151-2916.2000.TB01680.X.

[105] B. Lawn, Fracture of Brittle Solids, Cambridge University Press, 1993, http://dx.doi.org/10.1017/CBO9780511623127.

[106] P. Guo, Y. Wu, R. Chen, B. Lu, Y. Jia, W. Zhang, B. Cao, S. Yu, J. Wei, Development of strong, stiff and lightweight compression-resistant mechanical metamaterials by refilling tetrahedral wireframes, Virtual Phys. Prototyp. 19 (1) (2024) http://dx.doi.org/10.1080/17452759.2024.2365852.

[107] L.R. Meza, A.J. Zelhofer, N. Clarke, A.J. Mateos, D.M. Kochmann, J.R. Greer, Resilient 3D hierarchical architected metamaterials, Proc. Natl. Acad. Sci. USA 112 (37) (2015) 11502–11507, http://dx.doi.org/10.1073/pnas.1509120112.

[108] J. Schöttner, M. Walet, E. Marques, R. Carbas, V. Adam, M. Hohl, P. Weißgraeber, P. Rosendahl, S. Kraus, F. Rheinschmidt, M. Hedvard, L.F.M. da Silva, J. Schweizer, A. van Herwijnen, Developing two multiaxial testing machines to link strength and microstructure of weak snow layers, in: International Snow Science Workshop Proceedings 2024, Tromsø, Norway, Copernicus GmbH, 2024.

[109] C. Soyarslan, S. Bargmann, M. Pradas, J. Weissmüller, 3D stochastic bicon tinuous microstructures: Generation, topology and elasticity, Acta Mater. 149 (2018) 326–340, http://dx.doi.org/10.1016/J.ACTAMAT.2018.01.005.