the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Snow mechanical property variability at the slope scale – implication for snow mechanical modelling

### Francis Gauthier

### Alexandre Langlois

Snow avalanches represent a natural hazard to infrastructure and backcountry recreationists. Risk assessment of avalanche hazard is difficult due to the sparse nature of available observations informing on snowpack mechanical and geophysical properties and overall stability. The spatial variability of these properties also adds complexity to decision-making and route finding in avalanche terrain for mountain users. Snow cover models can simulate snow mechanical properties with good accuracy at fairly good spatial resolution (around 100 m). However, monitoring small-scale variability at the slope scale (5–50 m) remains critical, since slope stability and the possible size of an avalanche are governed by that scale. To better understand and estimate the spatial variability at the slope scale, this work explores links between snow mechanical properties and microtopographic indicators. Six spatial snow surveys were conducted in two study areas across Canada. Snow mechanical properties, such as snow density, elastic modulus and shear strength, were estimated from high-resolution snow penetrometer (SMP) profiles at multiple locations over several studied slopes, in Rogers Pass, British Columbia, and Mt. Albert, Québec. Point snow stability metrics, such as the skier crack length, critical propagation crack length and a skier stability index, were derived using the snow mechanical properties from SMP measurements. Microtopographic indicators, such as the topographic position index (TPI), vegetation height and proximity, wind-exposed slope index, and potential radiation index, were derived from unoccupied aerial vehicle (UAV) surveys with sub-metre resolution. We computed the variogram and the fractal dimension of the snow mechanical properties and stability metrics and compared them. The comparison showed some similarities in the correlation distances and fractal dimensions between the slab thickness and the slab snow density and also between the weak layer strength and the stability metrics. We then spatially modelled snow mechanical properties, including point snow stability, using spatial generalized additive models (GAMs) with microtopographic indicators as covariates. The use of covariates in GAMs suggested that microtopographic indicators can be used to adequately estimate the variation in the snow mechanical properties but not the stability metrics. We observed a difference in the spatial pattern between the slab and the weak layer that should be considered in snow mechanical modelling.

- Article
(8308 KB) - Full-text XML
- BibTeX
- EndNote

Snow avalanches represent a natural hazard to infrastructure and backcountry recreationists in mountainous areas all over the world (Stethem et al., 2003; Techel et al., 2016). Snow avalanches can be divided into different types: wet, dry, non-cohesive or slab avalanches. However, dry-snow slab avalanches are the most difficult to predict and are responsible for the most fatalities (Techel et al., 2016). They require a shear crack, usually initiated by a person or new snowfall in a weak porous layer underneath a cohesive snow slab. Then, the crack must reach a critical size in order to self-propagate across the slope for a slab avalanche to occur. Practitioners and forecasters estimate the probability and size of an avalanche from point-scale information on weak layers and slab properties across different scales. However, the sparse nature of available observations on snowpack properties makes the forecasting of slab avalanches difficult (Hägeli and McClung, 2004). The snow spatial variability at different scales also adds complexity to this challenging task by adding uncertainty to whether the properties measured in the field are representative of the slab and weak layer system (Schweizer et al., 2008a).

The spatial variability of snow properties is well documented in climate studies (e.g. Harper and Bradford, 2003), glacier dynamics (e.g. Pulwicki et al., 2018), snow hydrology (e.g. Deems et al., 2006), mountain meteorology (e.g. Mott et al., 2011), permafrost (e.g. Wirz et al., 2011) and snow (e.g. Schweizer et al., 2008a). Numerous studies have investigated the spatial distribution of snow depth and its water equivalent to support hydrological models (e.g. Deems et al., 2006; Grünewald et al., 2010; Schirmer et al., 2011; Winstral et al., 2002). Some researchers went further to estimate and analyze the spatial pattern of snow depth (Deems et al., 2006; Mott et al., 2011; Schirmer and Lehning, 2011; Trujillo et al., 2007). They analyzed the scaling properties and the fractal dimension of the snow depth, which can be estimated with the slope of a log–log variogram or with the periodogram of the spatial signal. The idea behind the scaling properties and fractal dimension is that many scales can define a spatial pattern instead of one scale like the correlation length in a variogram. Fractal dimension also characterizes the roughness or smoothness of a spatial pattern across multiple scales. These researchers compared the fractal dimension of snow depth with the fractal dimension of topographic indicators and vegetation. However, no studies have explored the fractal dimension of snow mechanical properties. Most studies have relied on lidar or manual snow probe surveys to estimate snow depth. However, snow depth is not a sufficient indicator of the conditions required for snow avalanches to occur.

There are more effective indicators, such as snow stability tests, to estimate the conditions for snow avalanches. These tests are widely used in the avalanche industry to assess snow stability and, ultimately, snow avalanche hazard. These tests provide a qualitative evaluation of the mechanical interaction between the cohesive slab and the weak layer. Some studies have investigated the variability of several snow stability tests on an avalanche-prone slope (Kronholm and Schweizer, 2003; Birkeland, 2001; Campbell and Jamieson, 2007). They demonstrated a variation in the test results and spatial patterns with variograms and correlation distances of around 5–20 m. However, these snow stability tests do not provide information on the snow mechanical properties of the slab and the weak layer. Additionally, these tests are time-consuming, leading to limited spatial sampling density and extent for statistical analysis, around 30 m measurements covering 20 m. To address this limitation, the high-resolution snow micro-penetrometer (SMP) was used to characterize the mechanical and structural properties of the snow, including slab and weak layer thickness, density, elastic modulus, and microstructural strength of the weak layer (Proksch et al., 2015; Löwe and van Herwijnen, 2012; Johnson and Schneebeli, 1999). Several studies have characterized stability based on snow mechanical properties of the slab and the weak layer (Föhn, 1987; Gaume and Reuter, 2017; Reuter et al., 2015b; Monti et al., 2016; Schweizer and Reuter, 2015; Reuter and Schweizer, 2018; Rosendahl and Weißgraeber, 2020). Gaume and Reuter (2017) proposed a stability index that accounts for both failure initiation and propagation propensity, using an analytical method applicable to SMP profiles.

The SMP was used in snow spatial studies because it can rapidly and accurately measure the mechanical properties of the snow relevant to snow stability on a slope prone to avalanche (Bellaire and Schweizer, 2011; Feick et al., 2007; Kronholm and Schweizer, 2003; Landry et al., 2004; Lutz et al., 2007; Lutz and Birkeland, 2011). These studies reported spatial patterns of weak layer properties with correlation distances ranging from 0.5 to 20 m. However, the sampling density of the survey was between 20 to 50 SMP measurements depending on the studies, and the spatial extent covered 20 to 50 m. Reuter et al. (2016) used stability metrics based on SMP-derived snow mechanical properties to show spatial patterns of snow stability with a higher sampling density and extent compared to the other studies. The correlation distance obtained from this study was still in the same range as the others with some exceptions between 40 and 60 m. The differences in spatial patterns of snow instability among surveys were attributed to various meteorological processes interacting with the terrain and snow cover (e.g. Schweizer et al., 2008a; Reuter et al., 2016).

Based on these findings, several studies have simulated artificial spatial patterns of the weak layer in mechanical models to understand the effect of the spatial variability of the weak layer on the slope stability, given the likelihood of an avalanche (Gaume et al., 2014, 2013; Kronholm and Birkeland, 2005; Fyffe and Zaiser, 2004). Gaume et al. (2015) used the same method to estimate the propensity for tensile failure in the slab and the relationship with the size of avalanche release. These studies typically assumed that the spatial patterns of the weak layer ranged from 0.5 to 10 m, with the other parameters being constant for simplicity. Kronholm (2004) and Bellaire and Schweizer (2011) demonstrated that the spatial patterns of the weak layer and the slab could have different correlation distances for the same survey, in some cases resulting in a smoother slab variation than the weak layer or the opposite. However, the spatial extent of the snow sampling was relatively small, only twice that of the measured correlation length, and could affect the estimation of the correlation length (e.g. Dale and Fortin, 2014; Skøien and Blöschl, 2006). This matter should be further explored with a spatial sampling extent greater than 20 m in order to improve the implementation of snow variability in mechanical models.

Spatial patterns of snow properties can be explained and estimated by statistical models with exploratory spatial variables. In the past, environmental variables were mapped using a linear regression model and kriging with external drift. Several studies have used kriging to map point snow stability, such as snow stability test results, SMP-derived mechanical properties, and stability metrics (Birkeland, 2001; Mullen and Birkeland, 2008; Reuter et al., 2015a; Schweizer and Kronholm, 2007). These studies have demonstrated that point snow stability can be spatially estimated using topographic indicators, such as aspect, elevation and slope angle, on the regional scale. These indicators capture the complex interactions between meteorological processes and terrain features, such as snow deposition by wind and the influence of solar radiation on the snow surface between different slopes (Reuter et al., 2016). However, despite the use of statistical models incorporating topographic indicators, spatially autocorrelated residuals persisted. This residual spatial variability could be attributed to other spatial phenomena at a smaller scale.

In studies focused on the slope scale, researchers successfully explained and estimated the spatial variability of snow depth, even in cases where slope angle, aspect and elevation remained relatively constant (e.g. Deems et al., 2006; Grünewald et al., 2010; Pulwicki et al., 2018; Revuelto et al., 2020; Meloche et al., 2022; Trujillo et al., 2007; Winstral et al., 2002). They used in their studies microtopographic indicators, such as the shape of the slope (topographic position index, TPI), vegetation index, and microclimate indexes such as wind exposure (Winstral index) or the potential of solar radiation. Guy and Birkeland (2013) related terrain parameters to potential trigger zones, but the relationships were not unique, and their study was limited to the presence of depth hoar layers. However, the presence of depth hoar crystals is insufficient to characterize snow stability, which requires more information on snow mechanical properties for the slab and the weak layer. These mechanical properties can be accurately measured with the SMP (Reuter et al., 2019). Reuter et al. (2016) have linked snow stability from SMP-derived snow mechanical properties with microtopography indicators at the basin scale. While previous spatial studies have explored linear relations between point snow stability and topographic indicators, Reuter et al. (2016) suggested that the relation between point snow stability and topographic indicators could be non-linear.

The snow mechanical variability can also affect the overall slope stability with the so-called knockdown effect (Fyffe and Zaiser, 2004; Gaume et al., 2014; Kronholm and Schweizer, 2003; Schweizer et al., 2008a). This effect denotes that variations in weak layer strength can cause the slope to fail before the load reaches the corresponding average strength, and this effect is more prominent with a longer correlation length. Additionally, spatial variation in snow can influence the size of the avalanche release (Gaume et al., 2015). Small-scale variation can promote slab tensile failure and smaller avalanches.

It is necessary to spatially explain and estimate the mechanical properties of snow and snow stability with microtopography indicators at the slope scale. This study is based on the limitations and suggestions of Reuter et al. (2016), who modelled the spatial patterns of two stability metrics at the basin scale with terrain-based indicators such as slope angle, aspect and elevation. This work aims to estimate spatial variation at a smaller scale using microtopographic indicators through non-linear regression. As such, the first objective of this paper is to compare the scaling effect of the snow mechanical properties and the stability metrics for slopes prone to avalanches with different characteristics. The second objective is to spatially estimate snow spatial variability using microtopography indicators. An additional objective is to compare our dataset with two empirical power law fits from the literature (Bažant et al., 2003; McClung, 2009), which estimate the shear strength of the weak layer and slab density from the slab thickness.

## 2.1 Study sites

In order to spatially estimate the spatial variability of snow using microtopography indicators, we selected four study sites based on their specific microtopography and microclimate context. The first study site is located on Mount Albert in Gaspésie National Park, Québec, Canada (Fig. 1a). The winter climate of the region is characterized by extreme changes caused by (1) low-pressure continental systems that bring heavy snowfall of up to 100 cm in 48 h followed by Arctic cold air masses with strong northwesterly winds and (2) warm and wet air masses coming from the south, creating rain-on-snow events (Meloche et al., 2018). The study site is named Arête de Roc (AR) and is located in a subalpine/tundra area heavily affected by wind and snow transport compared to the other sites. This site has a high ground roughness with large boulders and small trees (1 m high). The slope angle is uniform (33°) with a convex roll at the top and a concavity at the bottom (Fig. 1). Two other surveys in Mt. Albert at Épaule du Mur (EP), where the snow slabs are thicker and denser, were added for our additional objective, namely to verify the parametrization of snow mechanical properties based on slab thickness (Bažant et al., 2003; McClung, 2009). However, these two surveys were not used in the variogram analysis and spatial modelling due to their insufficient spatial density and extent compared to the other surveys. They were added to the study to increase the data range for our additional objective.

Two study sites are in Glacier National Park, located in Rogers Pass, British Columbia, Canada (Fig. 1). Our study sites are on Mount Fidelity, which receives heavy snow precipitation (Hägeli and McClung, 2003) and has a snow cover of around 2–3 m and sometimes up to 4 m. The Mount Fidelity area is classified as a transitional snow and avalanche climate influenced by warm and wet air masses from the Pacific that bring heavy snowfall and cold air masses from the north, leading to the development of persistent weak layers (Hägeli and McClung, 2003). This study area annually experiences several persistent weak layers consisting of buried surface hoars or facets, relevant for stability assessment. The first study site at Mount Fidelity is located just above the tree line at 2020 m a.s.l. on a shoulder named Round Hill (RH). This site is an alpine area with low soil roughness (Fig. 1). The slope angle is relatively low (near 25°), with long and smooth convex rolls around 20–30 m. The last study site, Jim Bay Corner (JBC), is located below the tree line at 1830 m a.s.l. It is an open forested area with relatively low ground roughness with small shrubs. The site has 10 m tall trees which create some shaded areas, and the slope angle is relatively constant (near 20°) with small convex rolls around 5–10 m (Fig. 1).

## 2.2 Data collection and sampling strategies

For the spatial analysis, this study presents four snow spatial surveys collected during the winter of 2021–2022 (Fig. 1): 25 February 2022 at the Arête de Roc site (AR22-PP), 27 January 2022 at the Round Hill site (RH22-PP), 19 January 2022 at Jim Bay Corner (JBC22-SH) and 24 January 2022 at Jim Bay Corner (JBC22-PP). Two more surveys were added for the comparison of different parametrizations of snow mechanical properties: 24 January 2019 at Épaule du Mur (EP19-FC) and 29 February 2020 at Épaule du Mur (EP20-DF). Snow mechanical properties were measured using the high-resolution SMP. To compare the spatial patterns of snow mechanical properties and snow stability, each SMP measurement was made following a sampling scheme, according to the concept of the scale triplet, which is the support, spacing, and extent described by Blöschl and Sivapalan (1995). The support is the diameter of the SMP tip, which is 5 mm, guaranteeing a proper estimation of the microstructural properties of the snow. A minimum spacing of 2 m and a study site extent covering around 60 to 100 m were chosen to allow the spacing to be at least half of the expected correlation length and the extent to be 2 to 5 times the expected correlation length. The expected correlation length has been reported to be around 5–20 m from several studies (Bellaire and Schweizer, 2011; Lutz et al., 2007; Reuter et al., 2016; Schweizer and Reuter, 2015). This method ensures a reliable estimate of the spatial pattern, defined by both spatial variance and autocorrelation distance (Skøien and Blöschl, 2006; Dale and Fortin, 2014). Our sampling scheme also needs to be adequate for the second objective, which is the spatial estimation of snow mechanical properties and stability metrics using microtopographic indicators. Therefore, the sampling scheme was adjusted for each specific study site to obtain a representative distribution of microtopographic indicator values while respecting the scale triplets mentioned above. The sampling was conducted by randomly traversing the study site while adhering to the minimum spacing and also by characterizing the down and cross-slope for an isotropic sampling. The sampling was stopped when the study site was almost covered with 60 to 80 SMP measurements. The resulting sampling is shown in Fig. 1. Random sampling contributes to the obtainment of a robust estimation of the correlation length with limited SMP measurements (Kronholm and Birkeland, 2007; Skøien and Blöschl, 2006).

To ensure an accurate interpretation of the SMP signal, the weak layer needed to be identified and characterized from a snow profile. Full characterization of the snow stratigraphy was not needed for our analysis, so a shorter version of snow profile was used to optimize the time in the field. Two or three snow profiles were conducted per snow spatial survey, spaced at least 20 m apart and positioned next to SMP measurements (Fig. 1). In each test snow profile, we first performed two compression tests to identify the weak layer (Canadian Avalanche Association, 2016). The weak layer was attributed to the uppermost compression test results consistent in both compression tests. Then, we visually characterized the types and sizes of the snow grains of the weak layer. Finally, a propagation saw test was performed to measure the critical crack length of the weak layer (Gauthier and Jamieson, 2008). Layers situated above the weak layer were considered part of the slab. This assessment allowed us to accurately identify the weak layer in the nearest SMP profile and subsequently in the remaining SMP profiles. Each snow measurement, SMP or snow profile, was georeferenced using a GNSS receiver with centimetre accuracy. Furthermore, for each study site, aerial imagery was captured by a quad-rotor UAV with an RGB sensor to characterize the topography in both the summer and the winter on the same day as the spatial snow survey. Ground/surface models were generated using a *structure from motion (sfm)* photogrammetry algorithm (Westoby et al., 2012) with ground and snow control points, ensuring georeferenced models with centimetre accuracy (<2 cm in *x* and *y* and <5 cm in *z*).

## 2.3 Snow mechanical properties and stability metrics

This section describes the workflow used to process every SMP profile, extracting several snow mechanical properties needed for stability assessment. Three stability metrics were derived from these snow mechanical properties. Figure 2 presents the summary of this workflow.

### 2.3.1 SMP signal processing and snow properties

Each SMP signal was visually interpreted to identify distinct layers. First, the weak layer was identified on the SMP signal next to the snow profile, based on the failure depth in the corresponding compression test. Then homogeneous layers above the weak layer were classified into slab layers (S_{1}, S_{2}, …, S_{i}). This procedure was repeated for the remaining SMP signal. To obtain the macroscopic mechanical properties for each snow layer, the SMP signal was analyzed using a Poisson shot noise model with a moving window of 2.5 mm (Löwe and van Herwijnen, 2012). This analysis was used to recover microstructural parameters, including the peak force *F* the deflection at rupture *δ*, and the element length *L* (Löwe and van Herwijnen, 2012). Then, each structural and macroscopic snow mechanical property necessary for estimating the stability metrics could be retrieved: the slab thickness *D*, the weak layer thickness *D*_{wl}, the slab density *ρ*, the weak layer density *ρ*_{wl}, the elastic modulus of the slab *E* and the shear strength of the weak layer *τ*_{p}. Specifically, the slab thickness *D* and the weak layer thickness *D*_{wl} are directly extracted from the SMP profile. Slab density *ρ* and weak layer density *ρ*_{wl} are derived using the *F* and *L* parameters based on the method proposed by Proksch et al. (2015):

where *ρ* is in kilograms per cubic metre, *L* in millimetres and *F* in N. The coefficients were obtained by Calonne et al. (2020). The slab density *ρ* was determined as the mean value of all sub-slab layers above the weak layer, while *ρ*_{wl} is the mean value of the signal inside the weak layer. The effective macroscale elastic modulus of the slab (*E*) was derived with the new formulation recently adapted by Reuter et al. (2019), originally developed by Johnson and Schneebeli (1999):

The SMP cannot specifically measure the shear strength of the weak layer due to the mixed-mode loading on the weak layer caused by the slope angle in the field. Reuter et al. (2015a) previously assumed that the shear strength of the weak layer *τ*_{p} is approximately equal to the microstructural strength of the element defined by ${\mathit{\sigma}}_{\mathrm{micro}}^{\mathrm{th}}=F/{L}^{\mathrm{2}}$. We retained the same assumptions, but we used the macroscale strength ${\mathit{\sigma}}_{\mathrm{macro}}^{\mathrm{th}}$ to estimate *τ*_{p}. The formulation is similar but scaled with the number of active contacts $\frac{\mathit{\delta}}{L}$ over the 2.5 mm processing moving window of the SMP, following the formulation of Johnson and Schneebeli (1999):

### 2.3.2 Stability metrics

The skier propagation index (SPI) proposed by Gaume and Reuter (2017) was used to describe the skier stability. The SPI is defined as the ratio of two lengths: the skier crack length *l*_{sk} and the critical crack length *a*_{c}. The skier crack length represents the length of the crack in the weak layer induced by the weight of a skier on top of a slab. The critical crack length is the length of the crack required for the onset of crack propagation. The skier crack length is computed by solving the following equation: $\mathit{\tau}+\mathrm{\Delta}\mathit{\tau}={\mathit{\tau}}_{p}$, where *τ*=*ρ**g**D*sin *ψ* is the shear stress due to the slab weight with *g* as the gravitational acceleration. The stress due to the skier Δ*τ* was originally defined by Föhn (1987) and refined by Monti et al. (2016):

where *R* is the skier load set to 780 N, and *ψ* is the snow surface slope angle derived from UAV imagery. The angle *α* is defined as the angle between the point at the snow surface under the skier load to the point of maximum induced shear stress at the weak layer. Additionally, *D*_{e} is the new multilayered slab thickness proposed by Monti et al. (2016), considering that slabs are often made up of multiple layers with different properties, influencing stress redistribution (Habermann et al., 2008). The computation of *D*_{e} follows Eqs. (2), (3) and (4) in Monti et al. (2016), based on each layer elastic modulus *E* that composed the slab. In order to determine *l*_{sk}, the roots of the equation are found where $\mathit{\tau}+\mathrm{\Delta}\mathit{\tau}={\mathit{\tau}}_{p}$. The roots define two angles, *α*_{1} and *α*_{2}, describing the area of stress from the surface beneath the skier to the weak layer. From these two angles, the skier crack length is calculated (*l*_{sk}) with the following equation (Gaume and Reuter, 2017):

It is important to note that *D*_{e} was used exclusively in Eqs. (4)–(5), and the slab thickness *D* was used in the *a*_{c} formulation (explained below) and in both spatial analysis and estimation.

The critical crack length is computed using the formulation from Gaume et al. (2017):

where *σ*=*ρ**g**D*cos *ψ*, and Λ is a characteristic length of the system defined by

where ${E}^{\prime}=E/(\mathrm{1}-{v}^{\mathrm{2}})$, *v* is Poisson's ratio set to 0.3, *D*_{wl} is the weak layer thickness and *G*_{wl} is the shear modulus of the weak layer. However, Richter et al. (2019) proposed to change the formulation of Λ by excluding *D*_{wl} due to its sensitivity in snow cover modelling (SNOWPACK), which is also challenging to visually interpret in an SMP profile. Instead, they proposed a parametrization based on the weak layer density and optical grain size, replacing the ratio $\frac{{D}_{\mathrm{wl}}}{{G}_{\mathrm{wl}}}$ by *F*_{wl} (Eq. 8) into the characteristic length $\mathrm{\Lambda}=\sqrt{{E}^{\prime}D{F}_{\mathrm{wl}}}$. They normalized the optical grain size with a critical grain size (1.25 mm) from Schweizer et al. (2008b). The critical grain size of 1.25 mm was determined with a statistical analysis comparing weak layer properties from profiles classified as stable or unstable. We adapted this approach by replacing the optical grain size with the SMP parameter *L*. Following Pielmeier and Marshall (2009), we used *L*_{0}=1.09 mm. Consequently, we obtained the following formulation for *F*_{wl}:

where *ρ*_{wl} is the weak layer density, and *L*_{wl} is the element length *L* of the SMP signal analysis averaged over the thickness of the weak layer. The values are slightly different from those reported by Richter et al. (2019). Additionally, critical crack lengths were obtained in the field with the propagation saw test (PST) conducted next to the snow profile for each snow sampling survey. We compare the critical crack length *a*_{c} from the SMP with the critical crack length from the PST tests to assess the precision of our approach. However, we do not intend to accurately predict the stability metrics but to model their spatial variation. Finally, the skier propagation index (SPI) is defined as the ratio of the critical crack length (*a*_{c}) and the skier crack length (*l*_{sk}) (Gaume and Reuter, 2017):

A snowpack loaded by a skier is considered stable for SPI >1 and unstable for SPI <1.

## 2.4 Analysis of spatial pattern

The first objective of this paper is to compare the scaling effect on snow mechanical properties and stability metrics for slopes prone to avalanches with different characteristics. We choose three mechanical properties: the slab thickness *D*, the slab density *ρ* and the shear strength of the weak layer *τ*_{p}, as well as the three stability metrics described above, which are the skier crack length *l*_{sk}, the critical crack length *a*_{c} and the skier propagation index SPI. The spatial patterns of each snow mechanical property and stability metric were compared between the snow spatial surveys as an exploratory analysis. The omnidirectional experimental variogram *γ* was computed following the equation for a variable *y* (Chilès and Delfiner, 1999):

where *N* is the number of observations, and *h* is the distance between observations. The experimental variogram is defined by three parameters: the nugget or the non-spatial variance, the sill (which is the spatial variance), and the range or correlation length (which is the distance where the variance levels out). While the sill is difficult to compare across properties due to differing units, the correlation length is comparable as it shares the same unit. The correlation length provides insight into the scale of spatial variation. Four types of covariance models (Gaussian, exponential, spherical, Matérn) were fitted to the experimental variogram using iterative reweighted least squares estimation with the function fit.variogram from the *gstat* package (Pebesma, 2004) in RStudio (R Core Team, 2013). Furthermore, the fractal dimension, which expresses the roughness or complexity of a surface (2D–3D) in a non-integer dimension (Gao and Xia, 1996), was estimated from the variogram. We estimated the slope *ϕ* of the transformed log–log variogram and then obtained the fractal dimension (Gao and Xia, 1996):

## 2.5 Spatial modelling

### 2.5.1 Covariate processing

The second objective of this study is to explore the link between microtopographic indicators and snow mechanical properties and stability metrics in order to estimate snow spatial variability. The scale of these microtopographic indicators is defined by the size of the moving window used to derive them. Different sizes of moving windows were used to allow for a multiscale approach describing the spatial process (e.g. Revuelto et al., 2020; Meloche et al., 2022; Veitinger et al., 2014). The choice of different window sizes used in this study is based on the literature and is developed further below. Microtopography indicators are used as exploratory spatial variables and are referred to as covariates in the spatial model. These covariates were derived from a digital terrain and surface model (DTM/DSM) generated through photogrammetry using the UAV imagery. The classification between the ground and the vegetation was performed manually through visual inspection, given the small extent of the study site. Additionally, canopy models were generated for each snow study site by differentiating the DSM from the DTM. Snow depth maps were generated using a snow surface model (DSM_{snow}) and compared to the DTM model to retrieve the snow depth for each spatial snow survey.

All covariates were raster data with an original spatial resolution below 0.1 m and were upscaled to a spatial resolution of 0.5 m. The final resolution of the spatial model was the same as the covariates. The choice of covariates was based on multiple studies that focus on spatial variation in snow depth and is described below. Three groups of covariates, terrain shape, vegetation and microclimate, are presented in Table 1. Two indicators were chosen to describe the terrain shape, the topographic position index (TPI) and the vector ruggedness measure (VRM). The TPI is a slope descriptor indicating ridges, valleys or slopes at a given scale, referencing the position in elevation relative to neighbouring cells (Weiss, 2001). The TPI was measured between a minimum radius and a maximum radius with weighted distance from the maximum radius (Table 1). The VRM indicates the ruggedness of the terrain independently of slope angle and aspect. The ruggedness was derived as the sum of elevation differences with neighbouring cells but then decoupled with slope angle and aspect, which means that a flat and a steep slope could be homogeneous with low ruggedness (Sappington et al., 2007). These two indicators were used to explain and estimate snow depth (e.g. Revuelto et al., 2020; Meloche et al., 2022; Veitinger et al., 2014). The sizes of the different moving windows were chosen based on the values used in these studies to have a multiscale approach (Table 1). We used the slope angle and convexity of the terrain as exploratory variables. Vegetation also has an impact on the spatial variation of snow depth (Deems et al., 2006); we choose to use the canopy height for the influence of shrubs (around 0.3 and 0.5 m) and small trees (around 1 or 2 m) because snow cover can be up to 3 or 4 m in some areas of JBC and RH. Only trees above 5 m were masked from the study sites. We used the radial proximity to vegetation greater than 2 m to represent proximity to trees. Some authors also found that solar radiation (e.g. Lutz and Birkeland, 2011) and wind exposure (e.g. Winstral et al., 2002) were important in spatially estimating snow properties. We selected the potentially incoming solar radiation as covariates: the algorithm simulated over a DSM (including trees), the trajectory of the sun in the sky based on the time of the year and the latitude of the study site. The covariate represents direct insolation (shade and sunshine areas), calculated over a month prior to the survey. The Winstral index or upwind maximum slope parameter *S*_{x} represents the shelter or exposure areas provided by the terrain upwind of each pixel (Winstral et al., 2002). The upwind terrain was defined with the maximum search distance and the prevalent wind direction based on the mean wind direction from the nearest weather station of the study sites over the winter. The snow depth values from the DSM_{snow} were taken as covariates. The last covariates used were the spatial coordinates (easting and northing). The fitting of a smooth function to spatial coordinates, explained in the following section, takes the residual spatial autocorrelation into account (Nussbaum et al., 2017). The processing of the covariates involved the use of the geoprocessing library SAGA (Conrad et al., 2015), QGIS 3.14 and a Python implementation of the Winstral index *S*_{x} according to Winstral et al. (2002).

### 2.5.2 General additive model

General additive models (GAMs) can represent non-linear relationships between the covariates and the response variable. GAMs have been used in the past for spatial estimation of environmental variables (Nussbaum et al., 2017). They produce good results while remaining easy to interpret compared to more complex tree classification methods and machine learning algorithms (Nussbaum et al., 2017). A GAM can be described as a generalized linear model with a linear prediction involving a sum of a smooth function *s* of covariate *x* (Wood, 2017):

where *f* is a link function to a family distribution, *Y*_{i} is a response variable from some exponential family distribution *K* and *X*_{i} is a row of the model matrix for any strictly parametric component with vector parameter *θ*. Each smooth function or spline *s*_{j} can be expressed through a basis expansion *b* with a weight parameter *β* and *k* defining the order of the basis expansion.

Each smooth function represents a combination of linear terms fitted to a covariate *x*_{j}. The order of the smooth function determines the non-linear degree or the *wigliness* of the fitted GAM. We kept a low order (*k*=3) to avoid overfitting and non-realistic variation. While stepwise procedures are commonly used, they lack stability compared to newer methods such as shrinkage and boosting procedures (Hesterberg et al., 2008). We used the double-penalty approach as a shrinkage method proposed by Marra and Wood (2011), which adds a smoothing parameter for each covariate spline function. This method is implemented in the *mgcv* package in R. We applied this method for six response variables *Y*: the three snow mechanical properties (slab thickness *D*, slab density *ρ*_{slab} and the shear strength of the weak layer *τ*_{p}) and the three stability metrics (skier crack length *l*_{sk}, critical crack length *a*_{c} and skier propagation index). The estimation of these response variables used GAMs with the 13 covariates listed in Table 1.

The performance of our models was evaluated with the root mean square error (RMSE) and the mean absolute error (MAE) using a 10-fold cross-validation approach. This involves randomly splitting the sample into 10 subsets, fitting the model to the 9 subsets, comparing it to the remaining subset and repeating this procedure 10 times. The percentage of deviance explained (sum of squared errors) was computed to demonstrate the amount of total variance accounted for by the model. This metric is more suited to non-linear models compared to *R*^{2}, which is still shown in the results for comparison. Once our model was fitted (and cross-validated) and the covariates were selected, the response variable was estimated for every location at each study site on a 0.5 m resolution grid. A smaller resolution will not be in line with the assumption of homogeneous snowpack for the computation of the skier crack *l*_{sk} and the critical crack length *a*_{c}. All statistical computations were performed in R (R Core Team, 2013).

## 3.1 Summary of spatial snow surveys

The first spatial snow survey was conducted at the AR site. A weak layer of precipitation particles with an observed grain size of 0.5–1 mm was investigated on 25 February 2022 (AR22-PP), with 45 SMP measurements and a spatial extent of 71 m. The average slab thickness was 0.28 m and the mean slab density was relatively high: 252 kg m^{−3} (Table 2).

At the RH site (RH22-PP), a weak layer of precipitation particles with an observed grain size of 0.5 to 1 mm was found beneath a relatively soft snow slab. The mean slab thickness was 0.19 m, and the mean density was 171 kg m^{−3}. This survey, conducted on 27 January 2022, included 64 SMP measurements and covered a spatial extent of 116 m. The slab consisted of one homogeneous layer of storm snow, and both the slab and the weak layer originated from the same meteorological event.

We conducted two spatial snow surveys at the JBC site in two different areas of the site. The first survey at this site took place on 19 January 2022 (JBC22-SH) when there was a persistent weak layer of buried surface hoar of 1–2 mm. The slab was composed of multiple layers with a mean slab thickness of 0.39 m and a mean density of 188 kg m^{−3} above the surface hoar crystals. This survey consisted of 53 SMP measurements, covering a spatial extent of 102 m. The second survey (JBC22-PP) was characterized by a weak layer of precipitation particles buried under a fresh snow slab of 0.21 m thickness and an average slab density of 166 kg m^{−3}, deposited by the same meteorological event as RH22-PP. This survey included 55 SMP measurements, and the spatial extent was 74 m (Table 2).

The last two surveys presented in Table 2 were added to the study to increase the data range in Fig. 3. The snow spatial survey EP20-DF had a mean slab thickness of 0.32 m and slab density of 241 kg m^{−3}, similar to AR22-PP. The snow spatial survey EP19-FC recorded the highest mean slab thickness of 0.85 m and the highest mean slab density of 333 kg m^{−3}. Although the number of SMP measurements and spatial extent were not sufficient for spatial analysis, these surveys provided valuable data points characterized by larger slab thickness *D*, contributing to a fair assessment of the two empirical power law fits (Bažant et al., 2003; McClung, 2009).

Figure 3 shows slab density *ρ* and weak layer shear strength *τ*_{p} in relation to slab thickness *D*. These relations are often established, as snow density and snow strength should increase as the slab load increases. We fitted two power laws to our SMP-derived dataset and compared them with two other empirical power laws commonly used in the literature (Bažant et al., 2003; McClung, 2009). Figure 3 indicates a poor fit for both parameters (*ρ* and *τ*_{p}). The power law from McClung (2009) was better suited for the two surveys characterized with relatively low density (*ρ*<250 kg m^{−3}), which were conducted at Mount Fidelity (Fig. 3a). The surveys with higher density (*ρ*>250 kg m^{−3}) were on Mount Albert, which is a heavily wind-exposed area that could explain these highly dense slabs. Figure 3b shows some surveys aligned with the two power laws, especially the surveys from Mount Fidelity (circles). However, the Mount Albert surveys exhibited more variability compared to the Mount Fidelity surveys. In general, our dataset fitted poorly with the power laws from the literature, and a lot of variability remained in each survey. The intra-survey variability and implications for snow mechanical modelling are discussed in Sect. 4.1.

## 3.2 Comparison of spatial patterns

For all spatial snow surveys, the empirical variogram showed smaller correlation lengths for the slab thickness compared to other properties, ranging from 5 to 10 m (Fig. 4). The variograms for the slab density exhibited a correlation length in the same range as for the slab thickness, particularly for JBC22-PP and RH22-PP, with 5 and 8 m, respectively. These two spatial snow surveys had the same weak layer and slab meteorological deposition event characterized by a new snow instability. The correlation length for the slab thickness and slab density at AR22-PP was 10 m, with the same type of new snow instability. The variogram of the slab density at JBC22-SH was the only survey that had a longer correlation length of 34 m. Variograms of the slab density from JBC22-SH, JBC22-PP and AR22-PP appeared to exhibit fractal characteristics. These variograms showed a distinct plateau of variance around 10–20 m, followed by an increase in variance around 30–40 m, indicating a multiscale pattern around these two distances (10 and 40 m). Variograms of the weak layer shear strength indicated a longer correlation length around 20 m compared to the ones of slab properties, which were around 10 m. In the JBC22-PP and RH22-PP surveys, which shared the same meteorological deposition event, the variance stabilized at 20 m without any further increase in variance. The other surveys (JBC22-SH and AR22-PP) had longer correlation lengths and showed fractal characteristics with no stabilization in variance as the sampling distance increased. The primarily used variogram models were spherical and exponential, characterized by a rapid increase in variance for short distances. These models tend to be less smooth than Gaussian models, which have a smaller variance for short distances. Gaussian models were fitted for slab thickness at JBC22-SH and slab density at JBC22-PP. In general, the correlation lengths tended to be shorter for the thickness and density of the slab compared to the shear strength of the weak layer in each snow spatial survey.

At first glance, all the correlation lengths of the stability metrics were longer than those of the slab properties. Surveys at the Jim Bay Corner (JBC22-SH and JBC22-PP) showed correlation lengths around 20 m (Fig. 5). The other two surveys (AR22-PP and RH22-PP) exhibited an empirical variogram that did not show a clear plateau of variance to determine a correlation length. These surveys either had a longer correlation length than the spatial extent of the sampling or showed a fractal behaviour over multiple scales. The correlation lengths of the stability metrics ranged from 10 to 20 m, which is longer compared to the slab properties (Fig. 5). The most frequently used variogram model was spherical, but Gaussian models were also applied for the skier crack length (JBC22-PP, RH22-PP, AR22-PP) and skier index (JBC22-SH, JBC22-PP). Gaussian models were more frequently fitted to stability metrics than to snow properties, suggesting smoother spatial patterns for the stability metrics. The variogram of the stability metrics exhibited more similarities with the variogram of the weak layer shear strength than with the slab properties.

The fractal dimensions for the snow properties indicated a difference in surface complexity between the slab properties, the weak layer properties and the stability metrics (Fig. 6). The slab properties had higher fractal dimensions, around 2.85, indicating a higher surface complexity compared to the weak layer and the stability metrics, which had a similar fractal dimension around 2.7. Despite the stability metrics being computed from both slab mechanical properties and weak layer properties, their fractal dimension values were closer to those of the weak layer than to those of the slab. These results suggest that the spatial patterns of the stability metrics were more similar to those of the weak layer than to those of the slab properties.

## 3.3 Spatial modelling

The spatial models created by the GAMs explained the variance of the response variable but not entirely. The *R*^{2} and the percentage of deviance explained before ranged from 0.17 to 0.84 and from 22 % to 84 % (Tables 3–4). On average for all models, the *R*^{2} was approximately 0.5, and the percentage of deviance was 55 %. The average *R*^{2} was 0.47 for snow properties and 0.55 for stability metrics, and the average percentage of deviance explained previously was the same at 55 %. The performance of the models was assessed with a 10-fold cross-validated RMSE and MAE. The cross-validated RMSE and MAE for the slab thickness *D* were mostly 1–2 cm, except for 12 cm at AR22-PP, and around 4 to 27 kg m^{−3} for the slab density. The RMSE and MAE for the shear strength ranged from 30 to 128 Pa, except for 752 Pa for AR22-PP, but this snow spatial survey was also the one which had the highest variance (500 to 3500 Pa).

The spatial surfaces estimated by the GAMs in JBC22-SH for the snow mechanical properties are presented in Fig. 7. The estimated surfaces for the slab thickness and density exhibited a similar variation with comparable maximum and minimum areas. However, the estimated surface for the shear strength of the weak layer differed slightly from the slab properties. This finding reinforces the above results, indicating that the spatial pattern of the weak layer differed from the slab properties in our dataset. Estimation errors for the critical crack length ranged from 0.11 to 0.60 m, except for 1.2 m for AR22-PP. The RMSE and MAE for the skier propagation index ranged from 0.27 to 4, showing significant variability and relatively high values for an index. The estimation errors for the stability metrics were notably high, demonstrating that the model estimations were not reliable compared to the snow mechanical properties. However, Fig. 8 suggests that some outliers might have contributed to overestimating the RMSE, particularly with low values of *l*_{sk} and high SPI values (SPI ≈10). The spatial patterns of the stability metrics revealed two major weak spots represented by two clusters of low SPI values near 0, located on the north side (right) and northwest (upper middle). These weak spots corresponded to areas with lower shear strength values and slightly thicker and higher-density slabs.

There are no clear covariates selected by the model for every site, snow properties or stability metrics. However, some covariates were selected more frequently by the spatial models than others. The covariates most frequently used by the models, for both snow properties and stability metrics, were multiscale TPI and VRM, but their usage varied depending on the scale (Fig. 9). Spatial models for the shear strength of the weak layer appeared to select mainly TPI2550 and VRM5, whereas for slab density, VRM15 and convexity were chosen predominantly. Canopy height was selected in the snow property models but rarely in the stability metric models. The easting and northing coordinates (*x*,*y*) were widely used in the models, indicating the presence of spatially autocorrelated residuals. Surprisingly, snow depth was not used as frequently as other covariates. Convexity was selected numerous times, especially for the slab density, but almost never for the slab thickness. Overall, these results indicate that there are no universal covariates or specific covariates for snow properties or stability metrics that could be extrapolated to other sites. The selection of covariates by the spatial models was site-specific and also specific to different snow properties. The spatial models presented using microtopography indicators were fairly reliable for estimating absolute values of snow properties and not reliable for the stability metrics but rather for capturing relative spatial variability.

## 4.1 Snow mechanical parametrization

Our study aligns with the well-known relationship between slab thickness and slab density, attributed to snow settlement. The comparison of spatial patterns between surveys indicated that these two properties exhibited similar trends in their variogram, the fractal dimension and their covariates used for spatial modelling. For further research, the empirical power law fit $\mathit{\rho}\sim \mathrm{100}+\mathrm{135}{D}^{\mathrm{0.4}}$, suggested by McClung (2009), provides a simple approach to obtain average values that represent the interaction between these two properties for mechanical simulation (e.g. Gaume and Reuter, 2017). The power law fitted to our SMP-derived dataset appears to yield better average values for denser slabs in wind-exposed areas. However, it is important to note that these power laws fitted poorly with our dataset, indicating that significant variability remains. Nevertheless, these power laws could be used in a snow mechanical model to generate a slab density variation according to the spatial pattern of the slab thickness. Until now in snow mechanical modelling research, the spatial variation in snow properties was limited to the weak layer. Our study shows a distinction between the spatial variation in the slab properties and the weak layer, already observed by Kronholm (2004) and Bellaire and Schweizer (2011). We propose accounting for both slab property variation and weak layer variation since spatial patterns can differ between them.

Weak layer variations exhibited longer correlation lengths (smoother spatial pattern) than slab variations, and the increase in shear strength did not necessarily match the increase in the slab thickness. In general, shear strength should increase with slab thickness due to the slab load, but some variation was still present in our dataset (Fig. 3). The interaction between slab thickness and shear strength can be described with a power law ${\mathit{\tau}}_{p}\sim c+\mathrm{1370}{D}^{\mathrm{1.3}}$ (Bažant et al., 2003), reported according to the Mohr–Coulomb relation with initial cohesion *c* (300 Pa in Fig. 3) (Gaume et al., 2014). This power law represented the average values of the survey from Mount Fidelity well, but our fitted power law could also be used for thicker (denser) slabs in wind-exposed areas. However, the four power laws tested did not adequately capture the variability in values for a specific spatial survey. The constant parameter must be adjusted for each spatial survey to fit the values. Overall, these power laws should be used with caution when estimating the average snow values (strength and density) if only the slab thickness is available.

Gaume et al. (2013) proposed a method to generate a weak layer with spatial heterogeneity. The method generates a random field with a specified mean, variance and correlation length for the cohesion of the weak layer, where the shear strength of the weak layer is defined by a Mohr–Coulomb relation. The friction term of the Mohr–Coulomb relation, which incorporates the slab load, was added to the cohesion to obtain the shear strength. Although their friction term was constant due to a constant slab thickness, the method can easily be adapted to accommodate a variable friction term, reflecting a variation in slab thickness. This adaptation would enable the specification of two distinct random fields for the properties of the slab and the weak layer while ensuring consistency with the load of the slab. This method still requires input values for mean, variance and correlation length. The empirical power law can estimate mean values, but according to our dataset, the variance is not well represented (Fig. 3). Future work should explore the possibility of estimating the variance and correlation length of snow properties using the covariance of microtopography combined with distributed snow cover modelling. Such approaches could contribute to more realistic simulations in avalanche modelling, enhancing forecasting capabilities for both the probability of skier triggering and the size of avalanche releases.

## 4.2 Spatial modelling

This study gathers a unique dataset characterizing the spatial variation in snow mechanical properties and stability metrics at four different study sites. The comparison of variograms and fractal dimensions highlights differences in scale between slab properties and, on the other hand, weak layer properties and stability metrics (smoother patterns). Spatial GAMs were used to estimate with fair accuracy the snow mechanical properties using microtopography. However, the spatial modelling of the stability metrics was poor and not reliable. Additionally, a portion of spatial variances remained unexplained by the models, potentially due to non-spatial variances, such as instrument error or our processing data strategy. This strategy included a visual interpretation of the layer in the SMP resistance profile, as misclassification or misidentification of the weak layer boundaries can impact the result. Nevertheless, the modification of using the parametrization *F*_{wl} proposed by Richter et al. (2019) instead of the weak layer thickness for the computation of the critical crack length makes the method less dependent on weak layer thickness, enhancing its robustness. While the cross-validated RMSE of snow mechanical properties suggests sufficient precision, the high RMSE of stability metrics indicates that the spatial modelling of these metrics is not reliable (Table 3). Future work could use spatial estimations of the snow mechanical properties to compute the stability metrics from the spatial field of snow properties.

The cross-validation procedure was performed by randomly selecting 10 subsets. Future work should consider the correlation length during the random selection of subsets in cross-validation procedures to ensure complete independence between subsets. This could improve the reliability of RMSE and MAE estimations. However, our 10-fold cross-validation (repeated 10 times) still provides a reliable estimation of the performance of our models.

## 4.3 Microtopographic covariates

This study aimed to use microtopographic covariates for spatial estimations of snow spatial variability and stability. Our spatial generalized additive modelling did not reveal a universal covariate that could spatially estimate both snow mechanical properties and stability metrics. The study of Reuter et al. (2016), based on larger-scale terrain-based covariates, did not find a consistent covariate in all surveys to estimate instability at the basin scale. They reported that the slope aspect was selected as a estimator by the model in all of their surveys, but each survey used a different combination of covariates. In the present study, the selection of covariates was specific to each survey with no clear trend or takeaway regarding the choice of covariates. Notably, snow depth was not a reliable spatial estimator of snow mechanical properties and stability metrics, a finding consistent with the study by Reuter et al. (2016). The limited selection of snow depth as an estimator in our study might be attributed to the homogeneity in the dataset regarding snow depth or the weak layer's spatial variation being unrelated to the snow accumulation process. It is also noteworthy that, despite AR22-PP being a wind-exposed site, the GAM did not select the Winstral index *S*_{x} as a predictor. This could be related to the research distance in *S*_{x} being too large (100 m), and adjusting the scale of this indicator, similar to TPI and VRM, could reveal *S*_{x} as a significant covariate, especially at the wind-exposed site (AR22-PP).

Unfortunately, no link could be made between our only persistent weak layer survey consisting of surface hoar crystals (JBC22-SH) and the remaining non-persistent weak layer surveys. A bigger dataset is needed to demonstrate clear differences between persistent vs. non-persistent weak layers, as well as between alpine vs. forested areas. The covariates TPI and VRM emerged as the most significant covariates for estimating snow properties; this was also observed by previous studies using spatial models (random forest) for snow depth estimation (Meloche et al., 2022; Revuelto et al., 2020). The optimal scale or window size for TPI and VRM varied depending on the study site, snow properties and stability metrics. Future work with a more extensive dataset should investigate whether the optimal scale is linked to the specific scale of terrain features at each site, the scale of the meteorological process affecting the slab and the weak layer, or a combination of both factors.

The transferability of our results to different sites is not feasible. The selection of covariates by the model was specific to each site, snow properties and stability metrics. As demonstrated by Reuter et al. (2016), the interaction between meteorological processes and terrain leads to distinct spatial variations in snow properties across different surveys. These micrometeorological processes vary between sites, and differences emerge not only between slab deposition patterns, but, crucially, between different types of weak layers. More spatial snow surveys are needed to gather a robust dataset to highlight trends in spatial patterns between different types of weak layers, slab deposition, microtopographic contexts and microclimatic contexts. To obtain a more robust dataset, future research should aim for an equivalent or higher sampling density and extent presented in this study (60 and more SMP profiles covering 80 m extent). Lowering the sampling density and extent could compromise the estimation of the experimental variogram and the spatial modelling. An alternative approach to sampling with fewer SMP measurements could be to incorporate distributed 3D snow cover modelling tools like ALPINE3D. This avenue was explored by Reuter et al. (2016), but they acknowledged the need for improving performance in distributed snow cover modelling. Implementing 3D snow cover modelling has the potential to capture a portion of these site-specific micrometeorological processes without requiring an extensive spatial survey of SMP measurements.

The study provides insights into the spatial variability of snow mechanical properties and stability metrics. First, we show that in our dataset, the slab properties exhibit spatial patterns that were different from the weak layer spatial patterns. In fact, the slab properties, both the slab thickness and the density, had smaller correlation lengths in their variograms than the weak layer strength. The slab properties had higher fractal dimensions than the weak layer strength, which demonstrates a more “rough” spatial surface. Secondly, spatial modelling using microtopography variables allows for the estimation of snow mechanical properties with reasonable accuracy, although the reliability of stability metric estimations was poor and not reliable. We also show the usefulness of using microtopography to estimate snow spatial variability, but the selection of the indicators was specific to each study site and the snow properties. The spatial models did not predominantly select a microtopographic indicator, indicating that there is no possible extrapolation to other sites or advice to backcountry recreationists. Future research could explore the capability of multiscale microtopographic indicators, like the topographic position index (TPI) and vector ruggedness measure (VRM), to estimate spatial patterns of snow mechanical properties with 3D snow cover modelling. This may contribute to the development of predictive methods for operational avalanche forecasting services, potentially estimating avalanche release sizes through snow cover modelling and mechanical models. Additional work is needed to gather a robust dataset regarding the spatial pattern of snow mechanical properties in order to elucidate trends between different types of weak layers and terrain features.

The code and the data are available upon request.

FM conceptualized and led the research, wrote the code for the processing and analysis of the data, and drafted the paper. FG and AL conceptualized the research and reviewed the paper. AL provided the major part of the funds for the project.

At least one of the (co-)authors is a member of the editorial board of *The Cryosphere*. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

The authors would like to thank Jeff Goodrich and the Mount Revelstoke Park and Glacier National Park staff for their support. This research was also possible with the help of Claude Isabel and the Gaspésie National Park (SEPAQ) as well as Dominic Boucher and Avalanche Québec staff. The authors would also like to thank Jean-Benoît Madore, Julien Meloche, Antoine Rolland, Alex Blanchette, Jacob Laliberté and William Durand for their help in the field. We want to thank the two anonymous reviewers for their helpful and constructive comments, which significantly improved the quality of our paper. Lastly, we want to thank Jürg Schweizer for his useful comments that improved the presentation of this paper.

This project was funded by the Search and Rescue New Initiatives Fund (SAR NIF) from Public Safety Canada, the Natural Sciences and Engineering Research Council of Canada (NSERC), and the Quebec Research Funds – Nature and Technologies (FRQNT), and the Canada Foundation for Innovation (CFI) provided funding for the Station d’études montagnardes des Chic-Chocs (SEM).

This paper was edited by Jürg Schweizer and reviewed by two anonymous referees.

Bažant, Z. P., Zi, G., and McClung, D.: Size effect law and fracture mechanics of the triggering of dry snow slab avalanches, J. Geophys. Res.-Sol. Ea., 108, 2119, https://doi.org/10.1029/2002jb001884, 2003. a, b, c, d, e, f

Bellaire, S. and Schweizer, J.: Measuring spatial variations of weak layer and slab properties with regard to snow slope stability, Cold Reg. Sci. Technol., 65, 234–241, https://doi.org/10.1016/J.COLDREGIONS.2010.08.013, 2011. a, b, c, d

Birkeland, K. W.: Spatial patterns of snow stability throughout a small mountain range, J. Glaciol., 47, 176–186, https://doi.org/10.3189/172756501781832250, 2001. a, b

Blöschl, G. and Sivapalan, M.: Scale issues in hydrological modelling: A review, Hydrol. Process., 9, 251–290, https://doi.org/10.1002/hyp.3360090305, 1995. a

Calonne, N., Richter, B., Löwe, H., Cetti, C., ter Schure, J., Van Herwijnen, A., Fierz, C., Jaggi, M., and Schneebeli, M.: The RHOSSA campaign: multi-resolution monitoring of the seasonal evolution of the structure and mechanical stability of an alpine snowpack , The Cryosphere, 14, 1829–1848, https://doi.org/10.5194/tc-14-1829-2020, 2020. a

Campbell, C. and Jamieson, B.: Spatial variability of slab stability and fracture characteristics within avalanche start zones, Cold Reg. Sci. Technol., 47, 134–147, https://doi.org/10.1016/j.coldregions.2006.08.015, 2007. a

Canadian Avalanche Association: Observations guidelines and recording standards for weather, snowpack and avalanches, Tech. rep., Revelstoke, 2016. a, b

Chilès, J.-P. and Delfiner, P.: Geostatistics: Modelling Spatial Uncertainty, John Wiley & Sons, Ltd, New-York, https://doi.org/10.1002/9781118136188, 1999. a

Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L., Wehberg, J., Wichmann, V., and Böhner, J.: System for Automated Geoscientific Analyses (SAGA) v. 2.1.4, Geosci. Model Dev., 8, 1991–2007, https://doi.org/10.5194/gmd-8-1991-2015, 2015. a

Dale, M. R. T. and Fortin, M.-J.: Spatial Analysis: A guide for Ecologists, Cambridge University Press, 2nd edn., https://doi.org/10.1017/CBO9780511978913, 2014. a, b

Deems, J. S., Fassnacht, S. R., and Elder, K. J.: Fractal distribution of snow depth from lidar data, J. Hydrometeorol., 7, 285–297, https://doi.org/10.1175/JHM487.1, 2006. a, b, c, d, e

Feick, S., Kronholm, K., and Schweizer, J.: Field observations on spatial variability of surface hoar at the basin scale, J. Geophys. Res.-Earth Surf., 112, 1–16, https://doi.org/10.1029/2006JF000587, 2007. a

Föhn, P.: The stability index and various triggering mechanisms, IAHS, 162, 195–214, 1987. a, b

Fyffe, B. and Zaiser, M.: The effects of snow variability on slab avalanche release, Cold Reg. Sci. Technol., 40, 229–242, https://doi.org/10.1016/j.coldregions.2004.08.004, 2004. a, b

Gao, J. and Xia, Z. G.: Fractals in physical geography, Prog. Phys. Geogr., 20, 178–191, https://doi.org/10.1177/030913339602000204, 1996. a, b

Gaume, J. and Reuter, B.: Assessing snow instability in skier-triggered snow slab avalanches by combining failure initiation and crack propagation, Cold Reg. Sci. Technol., 144, 6–15, https://doi.org/10.1016/j.coldregions.2017.05.011, 2017. a, b, c, d, e, f

Gaume, J., Chambon, G., Eckert, N., and Naaim, M.: Influence of weak-layer heterogeneity on snow slab avalanche release: Application to the evaluation of avalanche release depths, J. Glaciol., 59, 423–437, https://doi.org/10.3189/2013JoG12J161, 2013. a, b

Gaume, J., Schweizer, J., Herwijnen, A., Chambon, G., Reuter, B., Eckert, N., and Naaim, M.: Evaluation of slope stability with respect to snowpack spatial variability, J. Geophys. Res.-Earth Surf., 119, 1783–1799, https://doi.org/10.1002/2014jf003193, 2014. a, b, c, d

Gaume, J., Chambon, G., Eckert, N., Naaim, M., and Schweizer, J.: Influence of weak layer heterogeneity and slab properties on slab tensile failure propensity and avalanche release area, The Cryosphere, 9, 795–804, https://doi.org/10.5194/tc-9-795-2015, 2015. a, b

Gaume, J., van Herwijnen, A., Chambon, G., Wever, N., and Schweizer, J.: Snow fracture in relation to slab avalanche release: critical state for the onset of crack propagation, The Cryosphere, 11, 217–228, https://doi.org/10.5194/tc-11-217-2017, 2017. a

Gauthier, D. and Jamieson, B.: Evaluation of a prototype field test for fracture and failure propagation propensity in weak snowpack layers, Cold Reg. Sci. Technol., 51, 87–97, https://doi.org/10.1016/J.COLDREGIONS.2007.04.005, 2008. a

Grünewald, T., Schirmer, M., Mott, R., and Lehning, M.: Spatial and temporal variability of snow depth and ablation rates in a small mountain catchment, The Cryosphere, 4, 215–225, https://doi.org/10.5194/tc-4-215-2010, 2010. a, b

Guy, Z. M. and Birkeland, K. W.: Relating complex terrain to potential avalanche trigger locations, Cold Reg. Sci. Technol., 86, 1–13, https://doi.org/10.1016/j.coldregions.2012.10.008, 2013. a

Habermann, M., Schweizer, J., and Jamieson, J. B.: Influence of snowpack layering on human-triggered snow slab avalanche release, Cold Reg. Sci. Technol., 54, 176–182, https://doi.org/10.1016/j.coldregions.2008.05.003, 2008. a

Hägeli, P. and McClung, D. M.: Avalanche characteristics of a transitional snow climate-Columbia Mountains, British Columbia, Canada, Cold Reg. Sci. Technol., 37, 255–276, https://doi.org/10.1016/S0165-232X(03)00069-7, 2003. a, b

Hägeli, P. and McClung, D. M.: Hierarchy theory as a conceptual framework for scale issues in avalanche forecast modeling, Ann. Glaciol., 38, 209–214, https://doi.org/10.3189/172756404781815266, 2004. a

Harper, J. T. and Bradford, J. H.: Snow stratigraphy over a uniform depositional surface: Spatial variability and measurement tools, Cold Reg. Sci. Technol., 37, 289–298, https://doi.org/10.1016/S0165-232X(03)00071-5, 2003. a

Hesterberg, T., Choi, N. H., Meier, L., and Fraley, C.: Least angle and l1 penalized regression: A review, Statistical Surveys, 2, 61–93, https://doi.org/10.1214/08-SS035, 2008. a

Johnson, J. B. and Schneebeli, M.: Characterizing the microstructural and micromechanical properties of snow, Cold Reg. Sci. Technol., 30, 91–100, https://doi.org/10.1016/S0165-232X(99)00013-0, 1999. a, b, c

Kronholm, K.: Spatial Variability of Snow Mechanical Properties with regard to avalanche formation, Ph.D. thesis, University of Zurich, Zurich, 2004. a, b

Kronholm, K. and Birkeland, K. W.: Integrating spatial patterns into a snow avalanche cellular automata model, Geophys. Res. Lett., 32, L19504, https://doi.org/10.1029/2005GL024373, 2005. a

Kronholm, K. and Birkeland, K. W.: Reliability of sampling designs for spatial snow surveys, Comput. Geosci., 33, 1097–1110, https://doi.org/10.1016/j.cageo.2006.10.004, 2007. a

Kronholm, K. and Schweizer, J.: Snow stability variation on small slopes, Cold Reg. Sci. Technol., 37, 453–465, https://doi.org/10.1016/S0165-232X(03)00084-3, 2003. a, b, c

Landry, C., Birkeland, K., Hansen, K., Borkowski, J., Brown, R., and Aspinall, R.: Variations in snow strength and stability on uniform slopes, Cold Reg. Sci. Technol., 39, 205–218, https://doi.org/10.1016/j.coldregions.2003.12.003, 2004. a

Löwe, H. and van Herwijnen, A.: A Poisson shot noise model for micro-penetration of snow, Cold Reg. Sci. Technol., 70, 62–70, https://doi.org/10.1016/j.coldregions.2011.09.001, 2012. a, b, c

Lutz, E. and Birkeland, K. W.: Spatial patterns of surface hoar properties and incoming radiation on an inclined forest opening, J. Glaciol., 57, 355–366, https://doi.org/10.3189/002214311796405843, 2011. a, b

Lutz, E., Birkeland, K. W., Kronholm, K., Hansen, K., and Aspinall, R.: Surface hoar characteristics derived from a snow micropenetrometer using moving window statistical operations, Cold Reg. Sci. Technol., 47, 118–133, https://doi.org/10.1016/j.coldregions.2006.08.021, 2007. a, b

Marra, G. and Wood, S. N.: Practical variable selection for generalized additive models, Computational Statistics and Data Analysis, 55, 2372–2387, https://doi.org/10.1016/j.csda.2011.02.004, 2011. a

McClung, D. M.: Dimensions of dry snow slab avalanches from field measurements, J. Geophys. Res.-Earth Surf., 114, F01006, https://doi.org/10.1029/2007JF000941, 2009. a, b, c, d, e, f, g

Meloche, F., Gauthier, F., Langlois, A., and Boucher, D.: The Northeastern Rainy Continental snow climate: A snow climate classification for the Gaspé Peninsula, Québec, Canada, in: International Snow Science Workshop, Innsbruck, Austria, 1025–1029, 2018. a

Meloche, J., Langlois, A., Rutter, N., McLennan, D., Royer, A., Billecocq, P., and Ponomarenko, S.: High-resolution snow depth prediction using Random Forest algorithm with topographic parameters: A case study in the Greiner watershed, Nunavut, Hydrol. Process., 36, e14546, https://doi.org/10.1002/HYP.14546, 2022. a, b, c, d

Monti, F., Gaume, J., van Herwijnen, A., and Schweizer, J.: Snow instability evaluation: calculating the skier-induced stress in a multi-layered snowpack, Nat. Hazards Earth Syst. Sci., 16, 775–788, https://doi.org/10.5194/nhess-16-775-2016, 2016. a, b, c, d

Mott, R., Schirmer, M., and Lehning, M.: Scaling properties of wind and snow depth distribution in an Alpine catchment, J. Geophys. Res.-Atmos., 116, 1–8, https://doi.org/10.1029/2010JD014886, 2011. a, b

Mullen, R. S. and Birkeland, K. W.: Mixed Effect and Spatial Correlation Models for Analyzing a Regional Spatial dataset, International snow science workshop, 21–27 September 2008, Whistler, Canada, 8, https://arc.lib.montana.edu/snow-science/item/65 (last access: 10 July 2023), 2008. a

Nussbaum, M., Walthert, L., Fraefel, M., Greiner, L., and Papritz, A.: Mapping of soil properties at high resolution in Switzerland using boosted geoadditive models, SOIL, 3, 191–210, https://doi.org/10.5194/soil-3-191-2017, 2017. a, b, c

Pebesma, E. J.: Multivariable geostatistics in S: the gstat package, Comput. Geosci., 30, 683–691, https://doi.org/10.1016/j.cageo.2004.03.012, 2004. a

Pielmeier, C. and Marshall, H. P.: Rutschblock-scale snowpack stability derived from multiple quality-controlled SnowMicroPen measurements, Cold Reg. Sci. Technol., 59, 178–184, https://doi.org/10.1016/j.coldregions.2009.06.005, 2009. a

Proksch, M., Löwe, H., and Schneebeli, M.: Density, specific surface area, and correlation length of snow measured by high-resolution penetrometry, J. Geophys. Res.-Earth Surf., 120, 346–362, https://doi.org/10.1002/2014JF003266, 2015. a, b

Pulwicki, A., Flowers, G. E., Radic, V., and Bingham, D.: Estimating winter balance and its uncertainty from direct measurements of snow depth and density on alpine glaciers, J. Glaciol., 64, 781–795, https://doi.org/10.1017/JOG.2018.68, 2018. a, b

R Core Team: R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, http://www.R-project.org/ (last access: 10 July 2023), 2013. a, b

Reuter, B. and Schweizer, J.: Describing Snow Instability by Failure Initiation, Crack Propagation, and Slab Tensile Support, Geophys. Res. Lett., 45, 7019–7027, https://doi.org/10.1029/2018GL078069, 2018. a

Reuter, B., Herwijnen, A. V., and Schweizer, J.: Simple drivers of snow instability, Cold Reg. Sci. Technol., 120, 168–178, https://doi.org/10.1016/j.coldregions.2015.06.016, 2015a. a, b

Reuter, B., Schweizer, J., and van Herwijnen, A.: A process-based approach to estimate point snow instability, The Cryosphere, 9, 837–847, https://doi.org/10.5194/tc-9-837-2015, 2015b. a

Reuter, B., Richter, B., and Schweizer, J.: Snow instability patterns at the scale of a small basin, J. Geophys. Res.-Earth Surf., 121, 257–282, https://doi.org/10.1002/2015JF003700, 2016. a, b, c, d, e, f, g, h, i, j, k

Reuter, B., Proksch, M., Löwe, H., Van Herwijnen, A., and Schweizer, J.: Comparing measurements of snow mechanical properties relevant for slab avalanche release, J. Glaciol., 65, 55–67, https://doi.org/10.1017/jog.2018.93, 2019. a, b

Revuelto, J., Billecocq, P., Tuzet, F., Cluzet, B., Lamare, M., Larue, F., and Dumont, M.: Random forests as a tool to understand the snow depth distribution and its evolution in mountain areas, Hydrol. Process., 34, 5384–5401, https://doi.org/10.1002/hyp.13951, 2020. a, b, c, d

Richter, B., Schweizer, J., Rotach, M. W., and van Herwijnen, A.: Validating modeled critical crack length for crack propagation in the snow cover model SNOWPACK, The Cryosphere, 13, 3353–3366, https://doi.org/10.5194/tc-13-3353-2019, 2019. a, b, c

Rosendahl, P. L. and Weißgraeber, P.: Modeling snow slab avalanches caused by weak-layer failure – Part 2: Coupled mixed-mode criterion for skier-triggered anticracks, The Cryosphere, 14, 131–145, https://doi.org/10.5194/tc-14-131-2020, 2020. a

Sappington, J. M., Longshore, K. M., and Thompson, D. B.: Quantifying Landscape Ruggedness for Animal Habitat Analysis: A Case Study Using Bighorn Sheep in the Mojave Desert, J. Wildlife Manage., 71, 1419–1426, https://doi.org/10.2193/2005-723, 2007. a

Schirmer, M. and Lehning, M.: Persistence in intra-annual snow depth distribution: 2.Fractal analysis of snow depth development, Water Resour. Res., 47, 1–14, https://doi.org/10.1029/2010WR009429, 2011. a

Schirmer, M., Wirz, V., Clifton, A., and Lehning, M.: Persistence in intra-annual snow depth distribution: 1.Measurements and topographic control, Water Resour. Res., 47, 1–16, https://doi.org/10.1029/2010WR009426, 2011. a

Schweizer, J. and Kronholm, K.: Snow cover spatial variability at multiple scales: Characteristics of a layer of buried surface hoar, Cold Reg. Sci. Technol., 47, 207–223, https://doi.org/10.1016/j.coldregions.2006.09.002, 2007. a

Schweizer, J. and Reuter, B.: A new index combining weak layer and slab properties for snow instability prediction, Nat. Hazards Earth Syst. Sci., 15, 109–118, https://doi.org/10.5194/nhess-15-109-2015, 2015. a, b

Schweizer, J., Kronholm, K., Jamieson, J. B., and Birkeland, K. W.: Review of spatial variability of snowpack properties and its importance for avalanche formation, Cold Reg. Sci. Technol., 51, 253–272, https://doi.org/10.1016/j.coldregions.2007.04.009, 2008a. a, b, c, d

Schweizer, J., McCammon, I., and Jamieson, J. B.: Snowpack observations and fracture concepts for skier-triggering of dry-snow slab avalanches, Cold Reg. Sci. Technol., 51, 112–121, https://doi.org/10.1016/J.COLDREGIONS.2007.04.019, 2008b. a

Skøien, J. O. and Blöschl, G.: Sampling scale effects in random fields and implications for environmental monitoring, Environ. Monitor. A., 114, 521–552, https://doi.org/10.1007/s10661-006-4939-z, 2006. a, b, c

Stethem, C., Jamieson, B., Schaerer, P., Liverman, D., Germain, D., and Walker, S.: Snow avalanche hazard in Canada – A review, Nat. Hazards, 28, 487–515, https://doi.org/10.1023/A:1022998512227, 2003. a

Techel, F., Jarry, F., Kronthaler, G., Mitterer, S., Nairz, P., Pavšek, M., Valt, M., and Darms, G.: Avalanche fatalities in the European Alps: long-term trends and statistics, Geogr. Helv., 71, 147–159, https://doi.org/10.5194/gh-71-147-2016, 2016. a, b

Trujillo, E., Ramírez, J. A., and Elder, K. J.: Topographic, meteorologic, and canopy controls on the scaling characteristics of the spatial distribution of snow depth fields, Water Resour. Res., 43, W07409, https://doi.org/10.1029/2006WR005317, 2007. a, b

Veitinger, J., Sovilla, B., and Purves, R. S.: Influence of snow depth distribution on surface roughness in alpine terrain: a multi-scale approach, The Cryosphere, 8, 547–569, https://doi.org/10.5194/tc-8-547-2014, 2014. a, b

Weiss, A.: Topographic position and landforms analysis, Poster presentation, ESRI User Conference, San Diego, CA, 64, 227–245, 2001. a

Westoby, M. J., Brasington, J., Glasser, N. F., Hambrey, M. J., and Reynolds, J. M.: “Structure-from-Motion” photogrammetry: A low-cost, effective tool for geoscience applications, Geomorphology, 179, 300–314, https://doi.org/10.1016/J.GEOMORPH.2012.08.021, 2012. a

Winstral, A., Elder, K., and Davis, R. E.: Spatial snow modeling of wind-redistributed snow using terrain-based parameters, J. Hydrometeorol., 3, 524–538, https://doi.org/10.1175/1525-7541(2002)003<0524:SSMOWR>2.0.CO;2, 2002. a, b, c, d, e, f

Wirz, V., Schirmer, M., Gruber, S., and Lehning, M.: Spatio-temporal measurements and analysis of snow depth in a rock face, The Cryosphere, 5, 893–905, https://doi.org/10.5194/tc-5-893-2011, 2011. a

Wood, S. N.: Generalized additive models: an introduction with R, 2nd edn., Chapman and Hall/CRC, https://doi.org/10.1201/9781315370279, 2017. a