Fractional snow-covered area: scale-independent peak of winter parameterization

The spatial distribution of snow in the mountains is significantly influenced through interactions of topography with wind, precipitation, shortwave and longwave radiation, and avalanches that may relocate the accumulated snow. One of the most crucial model parameters for various applications such as weather forecasts, climate predictions and hydrological modeling is the fraction of the ground surface that is covered by snow, also called fractional snow-covered area (fSCA). While previous subgrid parameterizations for the spatial snow depth distribution and fSCA work well, performances were scale-dependent. Here, we were able to confirm a previously established empirical relationship of peak of winter parameterization for the standard deviation of snow depth σ_HS by evaluating it with 11 spatial snow depth data sets from 7 different geographic regions and snow climates with resolutions ranging from 0.1 to 3 m. An enhanced performance (mean percentage errors, MPE, decreased by 25 %) across all spatial scales ≥ 200 m was achieved by recalibrating and introducing a scale-dependency in the dominant scaling variables. Scale-dependent MPEs vary between −7 % and 3 % for σ_HS and between 0 % and 1 % for fSCA. We performed a scale- and region-dependent evaluation of the parameterizations to assess the potential performances with independent data sets. This evaluation revealed that for the majority of the regions, the MPEs mostly lie between ±10 % for σ_HS and between −1 % and 1.5 % for fSCA. This suggests that the new parameterizations perform similarly well in most geographical regions.

1 Introduction

Whenever there is snow on the ground, there will be large spatial variability in snow depth. In mountainous terrain, this spatial distribution of snow is significantly influenced by topography due to corresponding spatial variations in wind, precipitation, and shortwave and longwave radiation and in steep terrain due to avalanches that may relocate the accumulated snow. As a result, the snow-covered landscape can consist of a complex mix of snow-free and snow-covered areas, especially in steep terrain or during snowmelt. A parameter which describes how much of the ground is covered by snow is the fractional snow-covered area (fSCA). Most of the time, fSCA is tightly linked to snow depth (HS) and in particular to its spatial distribution. A fSCA is able to bridge the spatial mean HS and the actual observed snow coverage. Sound fSCA models are therefore crucial since for the same spatial mean HS in early winter and in late spring, the associated fSCA can be completely different (e.g., Luce et al., 1999; Niu and Yang, 2007; Magand et al., 2014).

A fSCA plays a key role in modeling physical processes for various applications such as weather forecasts (e.g., Douville et al., 1995; Doms et al., 2011), climate simulations (e.g., Roesch et al., 2001; Mudryk et al., 2020) and avalanche forecasting (Bellaire and Jamieson, 2013; Horton and Jamieson, 2016; Vionnet et al., 2016). As climate warms, fSCA is an highly relevant indicator for spatial snow-cover changes in climate projections (e.g., Mudryk et al., 2020). A decrease in spatial snow extent prominently changes surface characteristics, such as albedo in mountain landscapes, leading to changes in surface radiation, which is a primary component of the surface energy balance. A fSCA is also a parameter in hydrological models to scale water discharges appropriately to help manage basin water supply (e.g., Luce et al., 1999; Thirel et al., 2013; Magnusson et al., 2014; Griessinger et al., 2016). Errors in fSCA estimates directly translate into errors of snowmelt rates and melt water discharge (Magand et al., 2014). Thus, accurately describing fSCAs is of key importance for multiple model applications in mountainous terrain where highly variable spatial snow distributions occur.

A fSCA can be obtained from satellite remote sensing products using optical imagery with varying spatiotemporal resolutions. For instance, Sentinel-2 gathers data at a spatial resolution of 10 to 20 m at frequent global revisit intervals (< 5 d, cloud permitting) (Drusch et al., 2012; Gascoin et al., 2019). The availability of satellite-derived fSCA remains, however, inconsistent due to time gaps between satellite revisits, data delivery and the frequent presence of clouds which obscure the ground, especially in winter in mountainous terrain, thus reducing the availability of images drastically (e.g., Parajka and Blöschl, 2006; Gascoin et al., 2015). Satellite-derived fSCAs can also not be used directly for forecasting. Alternatively, fSCAs can be obtained from spatially averaging by using snow models at subgrid scales. While such snow-cover models are available (e.g., Tarboton and Luce, 1996; Marks et al., 1999; Lehning et al., 2006; Essery et al., 2013; Vionnet et al., 2016), up until now they could not be used at very high spatial resolutions over very large regions in part due to a lack of detailed input data such as fine-scale surface wind speed and precipitation, as well as due to high computational cost. Often they are limited by model parameters and structure requiring calibration. Integrating data assimilation algorithms in snow models is able to mitigate some of these limitations, which has led, for instance, to improvements in runoff simulations (e.g., Andreadis and Lettenmaier, 2006; Nagler et al., 2008; Thirel et al., 2013; Griessinger et al., 2016; Huang et al., 2017; Griessinger et al., 2019). However, uncertainties inherently present in the input or assimilation data still remain, which are generally accentuated over snow-covered catchments (Raleigh et al., 2015). Today, fSCA parameterizations describing the subgrid snow depth variability therefore remain unavoidable for complex model systems and for complementing the assimilation of satellite-retrieved fSCA products especially over mountainous terrain.

A parameterization of fSCAs describes the relationship between fSCA and grid-cell-averaged HS or snow water equivalent (SWE) by a so-called snow-cover depletion (SCD) curve. SCD curves were originally introduced in models without taking into account subgrid topography or vegetation. In principle, there are two commonly applied forms: so-called closed functional forms and parametric probabilistic SCD curve formulations (Essery and Pomeroy, 2004). Parametric SCD curves have disadvantages for practical applications such as numerical stability, computational efficiency and assuming an unimodal distribution which might be less appropriate for large grid cells covering heterogeneous surface such as mountainous terrain (e.g., Essery and Pomeroy, 2004; Swenson and Lawrence, 2012). Various closed functional forms for fSCAs are therefore applied in land surface and climate models (e.g., Douville et al., 1995; Roesch et al., 2001; Yang et al., 1997; Niu and Yang, 2007; Su et al., 2008; Swenson and Lawrence, 2012). Most of these parameterizations use simple relationships between fSCA and spatial mean HS or SWE. Since topography strongly determines the spatial snow depth or snow water equivalent distribution (Clark et al., 2011), in the past, terrain characteristics were mostly heuristically introduced in closed form curves to account for subgrid terrain influences on fSCA (e.g., Douville et al., 1995; Roesch et al., 2001; Swenson and Lawrence, 2012). To verify the commonly applied closed forms of fSCA, Essery and Pomeroy (2004) integrated log-normal SWE distributions and fitted the parametric SCD curves. The best fit obtained resulted in a function proportional to tanh, which is a previously derived closed form from Yang et al. (1997). By using a normal probability density function (pdf), Helbig et al. (2015) obtained the same form fit for fSCA as Essery and Pomeroy (2004). The functional form for fSCA from Yang et al. (1997) could thus be inferred from integrating normal, as well as log-normal, snow depth distributions with subsequent fitting of the parametric SCD curves. The main difference between the form of Yang et al. (1997) and Essery and Pomeroy (2004) is the variable in the denominator. Yang et al. (1997) used the aerodynamic roughness length, whereas Essery and Pomeroy (2004) obtained the standard deviation of snow depth (σ_HS) at the peak of winter in the denominator. The advantage of introducing σ_HS in the closed form for fSCAs is that subgrid terrain characteristics contributing to shape the dominant spatial snow depth distribution can be used to parameterize σ_HS and thus to extend the fSCA parameterization of Essery and Pomeroy (2004) for mountainous terrain (Helbig et al., 2015).

Until recently, it was not possible to derive an empirical parameterization for σ_HS based on high-resolution snow depth data due to the lack of such high-resolution spatial data. New measurement methods such as terrestrial laser scanning (TLS), airborne laser scanning (ALS) and airborne digital photogrammetry (ADP) nowadays provide a wealth of spatial snow data at fine-scale horizontal resolutions. Since recently, digital photogrammetry can also be applied to high-resolution optical satellite imagery (Marti et al., 2016; Deschamps-Berger et al., 2020; Eberhard et al., 2021; Shaw et al., 2020). Snow depth data at these high resolutions now enable statistical analyses of spatial snow depth patterns for various purposes (e.g., Melvold and Skaugen, 2013; Grünewald et al., 2013; Kirchner et al., 2014; Grünewald et al., 2014; Revuelto et al., 2014; Helbig et al., 2015; Voegeli et al., 2016; López-Moreno et al., 2017; Helbig and van Herwijnen, 2017; Skaugen and Melvold, 2019). Based on spatial snow depth data sets, σ_HS could be related to terrain parameters. For instance, Helbig et al. (2015) parameterized σ_HS at the peak of winter using spatial mean HS and subgrid terrain parameters, namely a squared-slope-related parameter and terrain correlation length, and Skaugen and Melvold (2019) parameterized σ_HS for the accumulation season using current spatial mean HS and stratifications according to landscape classes and standard deviations in squared slope. Though both approaches are promising and also somehow similar, e.g., both use the squared slope as a significant scale variable, they also differ, e.g., in the considered horizontal scale lengths at the development of the parameterization. While the parameterization of Helbig et al. (2015) was developed for squared grid cell sizes from 50 m to 3 km, Skaugen and Melvold (2019) presented parameterizations for 0.5 km × 1 km grid cells. Helbig et al. (2015) observed improved performances for larger scales (> 1000 m), Skaugen and Melvold (2019) observed the same performances when validating it for 0.5 km × 10.25 km grid cells. This can be explained by the physical processes shaping the complex mountain snow cover which is predominantly interacting with topography at different length scales, e.g., precipitation, wind and radiation (Liston, 2004). A multiscale behavior has been found in various studies using different spatial coverages and measurement platforms (e.g., Deems et al., 2006; Trujillo et al., 2007; Schirmer et al., 2011; Mendoza et al., 2020), but a thorough analysis of spatial autocorrelations using many spatial snow depth data sets up to several kilometers at horizontal resolutions far below the first estimated scale break of about 10 to 20 m has not been presented so far. Such an analysis could reveal a scale range from which the spatial snow distribution in mountainous terrain can be parameterized with consistent accuracy. Using the newly available wealth of spatial snow data, we now have the opportunity to better understand the differences in previous empirically developed closed form fSCA parameterizations by adding variability in the evaluation of data sets, i.e., by using data from different geographic regions, as well as by taking into account the spatial scale in scaling parameters.

This article presents a new fSCA parameterization for mountainous terrain for various snow model applications. Since snow model applications operate at different spatial scales, a fSCA parameterization should work across spatial scales, as well as for various snow climates. Two important points were therefore tackled compared to a previous fSCA parameterization. (1) We derived the empirical parameterization for σ_HS from a large pool of spatial snow depth data sets at the peak of winter from various geographic sites and validated it scale- and region-dependently. (2) Based on a spatial scale analysis, we introduced scale-dependent parameters in peak of winter parameterization of Helbig et al. (2015) for σ_HS such that the new fSCA parameterization can be reliably applied for grid cell sizes starting at 200 m and increasing to 5 km. While a seasonal fSCA model algorithm can be built using parameterized σ_HS at the peak of winter, we need additional information on the history of previous HS and SWE values to mimic the real seasonal fSCA evolution. We will present a seasonal fSCA algorithm separately.

2 Data

We compiled 11 spatial snow depth data sets from 7 different geographic sites in mountainous regions of Switzerland, France and the US, i.e., from two continents (Fig. 1). These data sets have horizontal grid cell resolutions Δx between 0.1 and 3 m and cover areas from 0.14 to 280 km². In addition to that, the snow depth data sets were acquired by five different remote sensing methods, i.e., using different platforms. The diversity of the data sets can be seen in Fig. 2, which shows the pdfs for snow depth, elevation and the squared-slope-related parameter μ (Helbig et al., 2015) which is described in Sect. 3.1. All snow depth data were gathered at the local approximate point in time when snow accumulations had reached their annual maximum. Except for the two snow depth data sets shown in Fig. 3, the data sets have been published before, or the geographic location is described elsewhere. In the following, all snow depth data sets are listed and grouped according to their mountain range.

Figure 1 The map shows the approximate location of the 11 spatial snow depth data sets. The colors of the trays indicate the region, measurement platforms or acquisition date as presented in Fig. 2.

Figure 2 Probability density functions for fine-scale (a) snow depth, (b) elevation and (c) squared-slope-related parameter per observation data set in its original horizontal resolution, i.e., between 0.1 and 3 m. Densities were normalized with the corresponding maximum density of all data sets. Note that for elevation (b) the y axis was cut for better visibility. Colors represent the different geographic regions, measurement platforms or acquisition dates (number) of the compiled data set as indicated in Sect. 2.1 to 2.4.

Figure 3 Snow depth maps of the eastern Swiss Alps: (a) in the Dischma region (ALS data) and (b) at Gaudergrat (UAS data) at the peak of winter. The red dot in the inset map for Switzerland shows the location of the two sites. Pixmap © 2020 Swisstopo (5704000000), reproduced with the permission of Swisstopo (JA100118).

2.1 Eastern Swiss Alps

We used snow depth data sets acquired by three different platforms at four different alpine sites in the eastern Swiss Alps.

The first platform was airborne digital scanning (ADS) using an optoelectronic line scanner on an airplane. Data were acquired from the Wannengrat and Dischma areas near Davos in the eastern Swiss Alps (Bühler et al., 2015). ADS-derived snow depth data sets were used from 20 March 2012 (“ads-CH²”) to 9 March 2016 (“ads-CH¹”), together with summer digital elevation models (DEMs) (Marty et al., 2019). The data set covers about 150 km² at 2 m resolution. Bühler et al. (2015) validated the 2 m ADS-derived snow depth data among others with TLS data. They obtained a root mean square error (RMSE) of 33 cm and a normalized median absolute deviation (NMAD) of the residuals (Höhle and Höhle, 2009) of 26 cm.

The second platform was an unmanned aerial system (UAS) recording optical imagery with real-time kinematic (RTK) positioning of the image acquisition points of the snow cover with a standard camera over two different smaller regions near Davos in the eastern Swiss Alps (Bühler et al., 2016; Eberhard et al., 2021). These images were photogrammetrically processed into a digital surface model (DSM). By subtracting the snow-free DSM from the winter flight, the HS values were obtained (Bühler et al., 2017). An UAS-derived snow depth data set was used from 7 April 2018 (“uav-CH⁹”) from Schürlialp, together with a UAS-acquired summer DEM (Eberhard et al., 2021). The Schürlialp data set covers about 3.2 km² which we used at 30 cm resolution. A second UAS-derived snow depth data set was used from 29 March 2019 (“uav-CH⁸”) from Gaudergrat, together with a UAS-acquired summer DEM. The Gaudergrat data set covers about 0.8 km² at 10 cm resolution (Fig. 3b). Compared to snow depth data from snow probing, Eberhard et al. (2021) obtained an RMSE of 16 cm and a NMAD of 11 cm for UAS-derived snow depth data at 9 cm horizontal resolution from Schürlialp.

The third platform was airborne laser scanning (ALS) above the Dischma region near Davos in the eastern Swiss Alps (Fig. 3a). This acquisition was a Swiss partner mission of the Airborne Snow Observatory (ASO) (Painter et al., 2016). For consistency reasons, the same lidar setup was used, and similar processing standards to the ASO campaigns in California were applied (Sect. 2.2). ALS-derived snow depth data were used from 20 March 2017 (“als-CH³”), together with a summer DEM from 2017. The ALS data set from Switzerland used here covers about 260 km² at 3 m resolution. Details on the derivation of the ALS data can be found in Mazzotti et al. (2019), though this study focused on three 0.5 km² forested sub-data sets. Validation of 1 m ALS-derived snow depth grids from 20 March 2017 against data from snow probing within the forest but outside the canopy (i.e., not below a tree) resulted in an RMSE of 13 cm and a bias of −5 cm.

2.2 The Sierra Nevada, CA, US

We used data sets acquired by two different platforms above Tuolumne basin in the Sierra Nevada (California) in the US.

The first platform was ALS performed by ASO (Painter et al., 2016). ALS-derived snow depth data were used from 26 March 2016 (“als-US⁷”) and 2 May 2017 (“als-US⁶”), together with a summer DEM (Painter, 2018). The second platform was a Pléiades product from 1 May 2017 (“plei-US⁶”). A detailed data description of the Pléiades data set derivation is given in Deschamps-Berger et al. (2020).

We used the ASO summer DEM for the Pléiades, as well as the ALS snow depth data sets. Given that the extent of the Pléiades snow depth data set was much smaller than the ALS domain, we cropped the ALS data sets to the Pléiades data set extension resulting in a coverage of about 280 km². The horizontal resolution used here was 3 m for both data sets. Compared to snow probe measurements in relatively flat areas, ALS snow depth data at 3 m horizontal resolution was found unbiased with an RMSE of 8 cm (Painter et al., 2016). Pléiades-derived snow depth data were recently validated with ASO data over 137 km² at 3 m resolution above Tuolumne basin (Deschamps-Berger et al., 2020). An RMSE of 80 cm, a NMAD of 69 cm and a mean bias of 8 cm were obtained for the Pléiades data set.

2.3 Eastern French Pyrenees

A Pléiades product was acquired over the Bassiès basin in the northeastern French Pyrenees. Pléiades-derived snow depth data were used from 15 March 2017 (“plei-FR⁴”), together with a summer DEM (Marti et al., 2016). The data set we used covers about 113 km² at 3 m resolution. Marti et al. (2016) derived a median of the bias between 2 m Pléiades data and snow probe measurements of −16 cm and with UAS measurements of −14 cm. They further obtained a NMAD of 45 cm with snow probe measurements and a NMAD of 78 cm with UAS measurements.

2.4 Southeastern French Alps

TLS-derived snow depth data were acquired at two alpine mountain passes in the southeastern French Alps. One snow depth data set was acquired over Col du Lac Blanc on 9 March 2015 (“tls-FR¹⁰”) (Revuelto et al., 2020). A site and data description can be found in Naaim-Bouvet et al. (2010), Vionnet et al. (2014), and Schön et al. (2015, 2018). We used a UAS-acquired summer DEM (Guyomarc'h et al., 2019). The data set covers about 0.6 km² at 1 m resolution. The second TLS-derived snow depth data set was acquired over Col du Lautaret at 27 March 2018 (“tls-FR⁵”) (Revuelto et al., 2020). We used a TLS-acquired summer DEM. The data set covers about 0.14 km² at 1 m resolution. Previously, mean biases between 4 and 10 cm for TLS laser target distances up to 500 m were obtained between TLS-derived and reference tachymetry measurements (Prokop, 2008; Prokop et al., 2008; Grünewald et al., 2010).

2.5 Preprocessing

In all data sets, grid cells Δx with forest, rivers, glaciers or buildings were masked out. In order to avoid introducing any biases, we consistently neglected fine-scale snow depth values in all data sets that were lower than 0 m or larger than 15 m. We used a snow depth threshold of 0 m to decide whether or not a fine-scale grid cell was snow-covered.

3 Methods

We parameterize the standard deviation of snow depth σ_HS to reassess the validity of the fSCA parameterization for complex topography from Helbig et al. (2015) for a range of spatial scales, in particular for sub-kilometer spatial scales.

3.1 Fractional snow-covered area parameterization

Helbig et al. (2015) derived an fSCA parameterization by integrating a normal pdf assuming spatially homogeneous melt. Subsequent fitting over a range of coefficients of variation CV (standard deviation divided by its mean) between 0.06 and 1.00 resulted in a similar closed form fit for fSCA as Essery and Pomeroy (2004) obtained by integrating a lognormal pdf,

$$\begin{array}{}\text{(1)}& \mathrm{fSCA}=\tanh \left(\mathrm{1.3}{\displaystyle \frac{\mathrm{HS}}{{\mathit{\sigma }}_{\mathrm{HS}}}}\right),\end{array}$$

using current HS and the standard deviation of previous maximum snow depth or peak of winter. The standard deviation of snow depth at the peak of winter was derived by relating the peak of winter high-resolution spatial snow depth data from Switzerland and Spain to underlying summer terrain parameters (Helbig et al., 2015):

$$\begin{array}{}\text{(2)}& {\mathit{\sigma }}_{\mathrm{HS}}={\mathrm{HS}}^a{\mathit{\mu }}^b\exp \left[-(\mathit{\xi }/L{)}^{\mathrm{2}}\right],\end{array}$$

with a=0.549, b=0.309, HS and terrain correlation length ξ are in meters, and ξ and μ are summer terrain parameters, where μ is related to the mean squared slope via $\mathit{\mu }={\left\{\stackrel{\mathrm{‾}}{\left[\right({\partial }_{x}z{)}^{\mathrm{2}}+\left({\partial }_{y}z{)}^{\mathrm{2}}\right]}/\mathrm{2}\right\}}^{\mathrm{1}/\mathrm{2}}$ using partial derivatives of subgrid terrain elevations z, i.e., from a DEM. The correlation length ξ or typical width of topographic features in a domain size L was derived via $\mathit{\xi }=\sqrt{\mathrm{2}}{\mathit{\sigma }}_z/\mathit{\mu }$ with the standard deviation of elevations σ_z. The $L/\mathit{\xi }$ ratio indicates how many characteristic topographic features of length scale ξ are included in each L. Figure 2 in Helbig et al. (2009) shows a transect of a topography indicating the described characteristic length scales. Similar to Helbig et al. (2015), we linearly detrended the summer DEM before deriving the terrain parameters to unveil the correct terrain characteristics associated with the shaping process of the snow depth distribution at the corresponding scale. Using Eq. (1), fSCA can thus be derived with grid cell mean snow depth from a snow model and grid cell mean subgrid terrain parameters derived from a fine-scale summer DEM. With the σ_HS formulation shown in Eq. (2), Helbig et al. (2015) extended the fSCA parameterization (Eq. 1) for mountainous terrain.

3.2 Aggregating and pooling of data sets

Pooling all snow depth data sets yields a data pool with a vast variety in snow climates, topographic characteristics and thus snow depth distributions. We first aggregated all snow data in squared domain sizes L in regular grids between 3 m and 5 km covering each geographic site. Our choice of the smallest applicable L in a data set was defined by a large enough $L/\mathrm{\Delta }x$ ratio (here ≥ 20) to minimize the influence of grid cell resolutions when spatially averaging (Helbig et al., 2009). When aggregating, we required at least 70 % valid data in a domain size which was the maximum threshold to obtain a sufficient number of domains for the largest domain sizes L of 3 m to 5 km. In addition to that, we excluded L with spatial mean slope angles larger than 60^∘ and spatial mean snow depth HS lower than 5 cm. By applying these limitations and since horizontal resolutions Δx, as well as the overall extent of the data sets, vary, the full range of L values, consisting of 41 different L values, was not represented by each data set. Overall, this resulted in a pool of 367 643 domains with L values between 3 m and 5 km. We obtain a decreasing number of domains for increasing L with a range between 59 376 for L=90 m and 17 for L=5000 m (Fig. 4). Spatial averages and standard deviations were built for each L. The resulting pooled data set shows a large variability in summer terrain characteristics. Spatial average slope angles range from 4 to 60^∘ (μ from 0.05 to 1.22; Fig. 5c), terrain correlation lengths ξ from 6 to 775 m and $L/\mathit{\xi }$-ratios from 3 to 40. Thus, typical summer terrain characteristics captured by coarse climate model grid cells are well represented. The diversity of the remaining domains with regards to snow depth is shown by means of the pdfs for spatial mean HS and σ_HS as a function of domain size L in Fig. 5a and b. Since the data pool also covers a broad range in spatial mean HS (from 5 cm to 12.4 m) and spatial variability in snow depth σ_HS (from 1 cm to 4.6 m) (Fig. 5a and b), we assume interannual snow depth variability is well described. In the following, overbars are neglected for spatial averages, i.e., for instance HS represents spatial mean snow depth exclusively.

Figure 4 Total number of valid domain sizes L per domain size L in log–log scale.

Figure 5 Probability density functions for (a) snow depth, (b) standard deviation of snow depth and (c) squared-slope-related parameter per domain size L after preprocessing, pooling all data sets and building spatial averages.

3.3 Autocovariances for scale breaks

The spatial autocovariance allows us to find spatial scale breaks up to which snow depth values are highly correlated, i.e., up to which length scale the snow depth distribution is strongly dictated by local topographic interactions of the snow cover with wind, precipitation and radiation. Below this scale, break process models should ideally explicitly resolve these interactions to reliably describe the spatial snow depth distribution. Above this scale break, we assume that dominant wind or precipitation patterns due to larger-scale topography impacts dictate spatial snow depth distributions. At this scale range, the normalized standard deviations of snow depth σ_HS start leveling out (Fig. 6a), as well as the normalized variability in σ_HS among similarly sized L (Fig. 6b).

Figure 6 (a) Normalized standard deviation of HS as a function of domain size L for each data set separately. (b) Normalized percentage of standard deviation of panel (a) among each L. L ranges from 3 m to 5 km. Colors represent the different geographic regions, measurement platforms or acquisition dates (number) of the compiled data set as indicated in Sect. 2.1 to 2.4.

We calculated spatial autocovariances for snow depth data sets with the fast Fourier transform (FFT), which allows for the computing of spatial autocovariances up to large distances by keeping the fine grid cell resolutions. We used the R function fft() of the “stats” package (see R Core Team, 2020). For each autocovariance, we then determine scale breaks using the R function uik() of the “inflection” package (R Core Team, 2020).

3.4 Deriving a new scale-independent fractional snow-covered area parameterization

Helbig et al. (2015) showed that fSCA performances increased with spatial scale and yielded their best performance for spatial scales larger than 1000 m. Since the fSCA parameterization was empirically developed on snow depth data from two geographic regions, here we reevaluated the scaling variables for the spatial variability in snow depth σ_HS, as well as the functional form of the parameterization using the much larger compiled HS data set of this study. Various scaling variables were previously employed to capture σ_HS in mountainous terrain. Helbig et al. (2015) selected HS, the squared-slope-related parameter μ and the $L/\mathit{\xi }$ ratio (Eq. 2). Skaugen and Melvold (2019) used HS and the standard deviation of the squared slope, with sqS being derived using sqS=tan ²ζ, where ζ is the slope angle in radians. Since tan ²ζ is the same as 2μ² (e.g., Löwe and Helbig, 2012), we here derive sqS from 2μ². Several other studies used σ_z as terrain parameter (e.g., Roesch et al., 2001). Here, we were interested in finding dominant scaling variables that correlate consistently across scales with σ_HS. We therefore analyzed the Pearson correlation coefficient r between various candidate parameters and σ_HS as a function of spatial scale, i.e., domain size L. Based on the results of previous studies, we selected the following candidate parameters: HS, μ, sqS, σ_sqS, $L/\mathit{\xi }$ and σ_z.

3.5 Performance measures

The performance in this article is evaluated by the following measures: the root mean square error (RMSE), normalized root mean square error (NRMSE; normalized by the range of measured data – max-min – or the mean of the measurements for fSCA), mean absolute error (MAE), the mean absolute percentage error (MAPE; absolute bias with measured minus parameterized and normalized with measurements), the mean percentage error (MPE; bias with measured minus parameterized and normalized with measurements) and the Pearson correlation coefficient r as a measure for correlation. We also evaluate the performances by deriving the two-sample Kolmogorov–Smirnov test (K-S test) statistic values D (Yakir, 2013) for the pdfs and by computing the NRMSE for quantile–quantile plots (NRMSE_quant; normalized by the range of measured quantiles, max-min) for probabilities with values in the range of [0.1,0.9].

4 Results

4.1 Spatial correlation range from snow depth data

We derived a total of 40 autocovariances for the available domain sizes L of 3 km with grid cell sizes Δx of 2 or 3 m. We obtained scale breaks between 183 and 681 m with a mean of 284 m (±σ 86 m) (Fig. 7). The zero crossings for each autocovariance were between 402 and 1815 m with a mean of 1011 m (±σ 402 m). For the mean autocovariance, we obtained a scale break at about 279 m and a zero crossing at about 1212 m. Based on the observed scale breaks, we selected a minimum length scale of 200 m for deriving a new scale-dependent fSCA parameterization for all larger scales. In the following, all results are therefore restricted to L≥ 200 m, leaving a pool of 41 249 domain sizes L with L between 200 m and 5 km for the development of the parameterization.

Figure 7 FFT-derived autocovariances for spatial snow depth. Individual ranges, mean range and mean autocovariance zero crossing are shown.

4.2 Scaling variables for σ_HS

Correlation coefficients varied differently across spatial scales (Fig. 8a). For all scales, we obtained the largest correlation coefficients for HS ranging from 0.48 to 0.98 with a mean of 0.79. From correlations with the various subgrid terrain parameters, the largest correlations across all scales were reached for the squared-slope-related parameter μ ranging from 0.22 to 0.61 with a mean of 0.36. Similar consistent correlation coefficients across scales but slightly smaller for L≤ 1800 m resulted in the squared slope sqS having an overall mean of 0.33. The correlation coefficients for the standard deviation of sqS (σ_sqS) and σ_z were much less consistent across scales than for μ and sqS and were overall lower. The mean correlation for σ_sqS is 0.15, for $L/\mathit{\xi }$ 0.21 and for σ_z 0.01. Though the mean correlation between σ_HS and $L/\mathit{\xi }$ is rather low, the correlation remains more consistent across scales up to about 2500 m and increases for larger scales considerably up to 0.67 (cf. Fig. 8a).

Figure 8 (a) Correlation coefficients between σ_HS and various parameters as a function of domain size L. (b) Standard deviation of snow depth σ_HS as a function of HS and μ.

We selected HS, μ and $L/\mathit{\xi }$ as main scaling parameters for σ_HS across spatial scales from 200 m to 5 km (Fig. 8b).

4.3 Scale-independent fSCA parameterization

The correlation analysis across scales revealed the same dominant correlation parameters as in Helbig et al. (2015). We therefore kept the functional form for σ_HS at the peak of winter suggested by Helbig et al. (2015) using the three scaling variables HS, μ and $L/\mathit{\xi }$. The new σ_HS parameterization at the peak of winter thus has the same functional form than the one suggested by Helbig et al. (2015) which was presented in Eq. (2). However, the fit parameters a and b therein are replaced by the new parameters c and d which we specify below. To derive the new parameters c,d, we fitted nonlinear regression models by robust M-estimators using iterated reweighed least squares; see R (R Core Team, 2020) and its robustbase package version 0.93-6 (Maechler et al., 2020). We started at the scale length of 200 m, defined by the scale break which we derived before from spatial snow depth autocovariances.

Fit parameters were first derived for the entire data pool and L≥200 m yielding c=0.6589 (±0.0037) and d=0.5638 (±0.0043) with the 90 % confidence intervals of the fit parameters given in parentheses. These “new” constant parameters c,d are larger than the previously derived constants a,b in Eq. (2) (cf. Fig. 9).

Figure 9 Fit parameters for Eq. (2) as a function of domain sizes L to scale variables (a) HS and (b) μ. Colored lines show the fit parameters derived for each of the 500 random 80 % samples of each of the 25 sub-data pools. The dark blue dots depict the ensemble median per L. Previously obtained constant parameters of Helbig et al. (2015) (light blue dots) and newly fitted constant (red dots), as well as newly fitted scale-dependent (pink circles) parameters are shown.

In addition to fitting over the entire data pool and L≥ 200 m, we split the entire data pool into 25 sub-data pools for any available domain size between 200 m and 5 km (cf. Fig. 4). Thereby, each sub-data pool comprised all domains larger than or equal to the corresponding domain size, i.e., L≥ 200 m, L≥ 240 m, etc. Fitting over such a sub-data pool should allow us to describe the combined larger-scale topography–wind–precipitation impacts on the spatial snow depth distribution in mountainous terrain acting at scales larger than the observed scale break of about 200 m. From each of the 25 created sub-data pools, we randomly took 500 subsamples in which each subsample was restricted to 80 % data of the sub-data pool. Each of the 500 subsamples per sub-data pool was unique. Scale-dependent parameter values were derived for each of the 500 subsamples drawn from each of the 25 sub-data pools (cf. individual colored lines in Fig. 9). Given that the values of c,d clearly increase with spatial scale L (Fig. 9), we introduced L in c,d to improve the application of Eq. (2) across scales. By fitting the ensemble medians of all scale-dependent fit parameters (dark blue dots in Fig. 9) across all domain sizes between 200 m and 5 km, we obtained scale-dependent parameters c(L) and d(L). Thus, Eq. (2) using the following scale-dependent parameters c(L) and d(L) assembles our new σ_HS parameterization for L≥ 200 m:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& c\left(L\right)=\mathrm{0.5330}{L}^{\mathrm{0.0389}},\\ & d\left(L\right)=\mathrm{0.3193}{L}^{\mathrm{0.1034}},\end{array}\end{array}$$

with the 90 % confidence intervals of ±0.0097, ±0.0026 and ±0.0183, and ±0.0079 in the order of introduced constants in Eq. (3).

The new σ_HS parameterization using c(L) and d(L) (Eq. 2 with Eq. 3) is applied in the previously derived fSCA parameterization (Eq. 1). To demonstrate the resulting differences when using scale-dependent versus scale-independent fit parameters in parameterized σ_HS (Eq. 2), we will also validate the performance using constant c,d in the previously derived fSCA parameterization, as well as in the σ_HS parameterization.

4.4 Evaluation

4.4.1 Evaluation for σ_HS and fSCA for all L values

Parameterized σ_HS and fSCA perform well for all domain sizes, i.e., for L≥ 200 m, of the entire data pool. Very similar performance measures are obtained for the parameterizations using the newly derived constant fit parameters c,d and the parameterizations using the scale-dependent parameters c(L),d(L) (cf. Table 1 and Ia and IIa). We obtain a slightly better MPE for σ_HS when using scale-dependent fit parameters (−4 % versus -5 %); however, for fSCA, MPEs are the same (0.2 %). The same rather low NRMSE results for σ_HS (8 %) and for fSCA (2 %) are obtained when using constant or scale-dependent fit parameters.

Helbig et al. (2015)Skaugen and Melvold (2019)Helbig et al. (2015)Skaugen and Melvold (2019)

Table 1 Performance measures for all L values between measurement and parameterization of (I) standard deviation of snow depth σ_HS with (a) Eq. (2) and constant or L-dependent fit parameters c,d (Eq. 3) and (b) σ_HS as in Helbig et al. (2015) and Skaugen and Melvold (2019) and of (II) fSCA with (a) Eq. (1) and (1a) and (b) fSCA as in Helbig et al. (2015) and Skaugen and Melvold (2019) using Eq. (1).

4.4.2 Scale-dependent evaluation for σ_HS and fSCA

While mean performance measures of the σ_HS and fSCA parameterization are almost uninfluenced by using constant or scale-dependent fit parameters (cf. Table 1 and Ia and IIa), we found diverging performances when analyzing performance measures as a function of scale (Fig. 10). Across scales, improved or similar performances were achieved when using scale-dependent fit parameters in parameterized σ_HS, especially for larger scales. Maximum performance improvements for σ_HS of 4 % and for fSCA of 0.7 % occurred for L of 2500 m when using scale-dependent fit parameters. Thus, introducing scale-dependent fit parameters enhanced the σ_HS parameterization for application across scales.

Figure 10 Mean percentage error (MPE) as a function of L for (a) σ_HS and (b) fSCA. MPEs are shown for the σ_HS and fSCA parameterizations using Eq. (1) to (3) with scale-dependent c(L),d(L), as well as for constant c,d. The second y axis shows the number of valid domains per L on a logarithmic scale.

4.4.3 Scale- and region-dependent evaluation for σ_HS and fSCA

A large data set from various geographic regions allows us to develop a more reliable empirical parameterization than being limited to the characteristics of a few data sets. Here, we not only compiled data sets from various geographic regions, but the data sets were also acquired by different measurement platforms coming with a range of inaccuracies between below 10 and 80 cm. As a consequence, larger scatter in performances appears when performance measures are depicted not only as a function of spatial scale but also regionwise, including platformwise. While most of the MPEs are still between −20 % and 10 %, some regions stand out because they have much larger MPEs when binned by scale, as well as by region (Fig. 11). For instance, a MPE of up to 60 % for σ_HS was obtained for TLS data from the southeastern French Alps, and overall larger MPEs, though consistent across scales, for the Pléiades data from the northeastern French Pyrenees were obtained. MPEs for fSCA on the other hand do not show a similar large spread among the regions and are low between −1 % to 2 % (Fig. 11b).

Figure 11 Mean percentage error (MPE) as a function of L for the compiled data set for (a) σ_HS and (b) fSCA using Eqs. (1) to (3) with scale-dependent c(L),d(L). Colors represent the different geographic regions, measurement platforms or acquisition dates (number) of the compiled data set as indicated in Sect. 2.1 to 2.4.

4.4.4 Evaluation of previous closed form parameterizations

To increase our understanding of the performances achieved with the new parameterizations, we also tested two previously derived empirical parameterizations. Specifically, we investigated how parameterized σ_HS using Eq. (2) (Helbig et al., 2015) and using the recently published formulation of Skaugen and Melvold (2019) compare to observed σ_HS of our compiled data set (Fig. 12a). We further tested both σ_HS parameterizations in the fSCA parameterization (Eq. 1; Fig. 12b). The parameterization of Helbig et al. (2015) works well. The performance measures for all L values are only slightly worse compared to the new parameterizations using both constant and scale-dependent fit parameters (Table 1). However, compared to the performance measures for the parameterization of Skaugen and Melvold (2019), the performances of Helbig et al. (2015) are clearly improved. Though MPEs of both previous σ_HS parameterizations are scale-dependent, the MPEs of Skaugen and Melvold (2019) reveal a larger scale-dependency of the performances compared to Helbig et al. (2015) (Fig. 12a). In particular, individual MPEs vary a lot from MPEs for all L values given in Table 1.

Figure 12 Mean percentage error (MPE) as a function of L for the compiled data set for (a) σ_HS and (b) fSCA. MPEs are shown for the σ_HS and fSCA parameterizations of Helbig et al. (2015), as well as for the σ_HS parameterization of Skaugen and Melvold (2019) and the σ_HS parameterization of Skaugen and Melvold (2019) applied in the fSCA parameterization of Helbig et al. (2015) (Eq. 1).

5 Discussion

5.1 Spatial correlation range

While multiscale behavior for spatial snow depth data has been found in various studies, observed scale breaks depend on the extent and horizontal resolution of the investigated snow depth data sets. A first scale break of spatial snow depth data in treeless, alpine terrain has been observed between 10 to 20 m (e.g., Deems et al., 2006; Trujillo et al., 2007; Schweizer et al., 2008; Schirmer and Lehning, 2011; Helfricht et al., 2014; Mendoza et al., 2020), and a second scale break has been observed at around 60 m (Trujillo et al., 2009). By computing spatial autocovariances starting with domain sizes L of 200 m at 0.1 to 1 m resolution and increasing up to 3 km at 2 to 3 m resolution, we also detected the two previously found scale breaks (not shown). However, by additionally covering larger spatial extents than have been previously investigated, we also detected a third scale break with a mean at about 280 m (Fig. 7). Similar long-range scale breaks between 185 and 300 m were very recently reported from analyzing 24 TLS-derived snow depth data sets acquired during six snow seasons in a subalpine catchment in the Spanish Pyrenees (Mendoza et al., 2020). Furthermore, a similar scale break at around 200 m was recently found by analyzing performance decreases in distributed snow modeling in various grid cell sizes, together with a semivariogram analysis of subgrid summer terrain slope angles in the same catchment in the High Atlas (Baba et al., 2019). While for other application studies, such as in avalanche forecasting, the smaller-scale breaks are decisive for explicitly describing the relevant snow-cover processes, here we are more interested in the largest detected scale break. Above scale lengths of 200 m, the longer-range processes of precipitation, wind and radiation interactions with topography most dominantly influence the spatial snow distribution in mountainous terrain, while we believe there are different physical processes which establish the smaller-scale breaks at around 10 to 20 and 60 m. The results presented here indicate that the model described by Eqs. (1) and (3) is reliably parameterizing the spatial snow distribution shaped by the longer-range precipitation, wind and radiation interactions with topography for spatial scales between 200 m and 5 km. Above the detected scale range of around 200 m, not only the spatial autocorrelations approach zero (Fig. 7), but normalized σ_HS clearly start leveling out, as well as the normalized variability in σ_HS among similarly sized L (Fig. 6). Thus, even though we could not verify the fSCA parameterization for length scales larger than 5 km, we believe that as long as grid cell mean slope angles are larger than zero, Eqs. (1) and (3) might also hold for larger grid cell sizes than 5 km.

5.2 Scaling parameter

We not only investigated dominant correlations between the spatial snow depth distribution and terrain parameters, but we also analyzed these correlations as a function of spatial scale. For some commonly applied scaling parameters, this revealed large variations in correlations across scales such as for σ_z (Fig. 8a). Similar to our results, Skaugen and Melvold (2019) also obtained large correlations between σ_HS and mean squared slope sqS for spatial snow depth data sets acquired at the peak of winter in Norway, though this was only analyzed for grid cells of 0.5 km² (0.5 km × 1 km). Nevertheless, this confirms our findings since mean squared slope is related to the slope-related parameter μ used here by sqS=2μ². However, Skaugen and Melvold (2019) obtained a slightly improved correlation for the standard deviation of squared slope and therefore selected this parameter to stratify the topography for parameterizing σ_HS. Across spatial scales, as well as for all L values, we obtained lower correlations between the standard deviation of squared slope and σ_HS, though we observed cross-correlations between the mean and the standard deviation of squared slope of 0.71, indicating that both parameters correlate well with σ_HS.

5.3 Scale-independent fSCA parameterization

The closed form fractional snow-covered area parameterization fSCA given in Eq. (1) got enhanced by recalibration and introducing scale-dependent fit parameters (Eq. 3) to make the performance consistent across spatial scales.

We developed the parameterization on a large snow depth data set. Large variability in the snow depth data set was gained by compiling 11 individual data sets from varying geographic regions, as well as various measurement platforms. While the latter might explain the remaining performance differences discussed below, the first led to large variability in summer terrain characteristics and snow climates and consequently spatial snow depth distributions (cf. Fig. 2). Though our presented parameterization for σ_HS was empirically derived, it is reassuring that for a new empirical derivation on a much larger and more diverse snow depth data set, the same underlying functional form could be used. Furthermore, larger (about 17 % and 45 %) but overall consistent constant fit parameters were obtained compared to those from Helbig et al. (2015) based on a more limited number of data sets and just two geographic regions (cf. a and b in Eq. 2 and c,d presented in Sect. 4.3 or Fig. 9).

In addition to deriving constant fit parameters across spatial scales, we took 500 random 80 % subsamples from each of the 25 sub-data pools (Sect. 4.3). Scale-dependent constants considerably increased with increasing scale from L = 200 m to L = 5 km by at most 12 % and 38 %, respectively (Fig. 9). This demonstrates that accounting for scale-dependent constants in the fSCA parameterization (Eq. 1 with Eqs. 2 and 3) had to be done. While we did not split our data set in development and validation subsets, fitting over the ensemble median of all scale-dependent parameters to derive c(L),d(L) ensures confidence in the resulting fit parameters.

An increase in scatter among all c(L) and d(L) with increasing domain scale L (Fig. 9) can be most likely explained by a concurrent decrease in available valid data in larger L values. Though we required at least 70 % valid data per L when aggregating fine-scale snow depth data in domain sizes L, the maximum threshold of 70 % was more often required for the larger L values than for smaller L values.

5.4 Evaluation

5.4.1 Evaluation for σ_HS and fSCA

Upon deriving performance measures on parameterized and observed σ_HS and fSCA for all L values (i.e., the pooled performance), we obtained very similar performances when using newly derived constant or scale-dependent fit parameters, i.e., c,d or c(L),d(L) (Table 1). Despite considerable differences up to 12 % for c and up to 38 % for d between constant and scale-dependent fit parameters (Fig. 9), pooled performances for all L values for σ_HS and fSCA were similar (Table 1). An explanation for this is that the number of available domains is strongly decreasing with increasing L. For L≥ 3000 m, we have only about 0.33 % (137 in total) valid domains available compared to the total of 41 249 for L< 3000 m (Fig. 4). This emphasizes the need for a scale-dependent evaluation.

5.4.2 Scale-dependent evaluation for σ_HS and fSCA

The largest improvement in MPE for all L values seems to originate from the recalibration using the new compiled data set with a reduction in MPE from −30 % to −5 % compared to a reduction from −5 % to −4 % when introducing scale-dependent fit parameters (Table 1). However, MPEs as a function of scale clearly demonstrated the improved behavior when using scale-dependent c(L),d(L) instead of constant fit parameters c,d in the σ_HS and fSCA parameterization (Fig. 10). Given that constant c,d were fitted over the entire data set as have been c(L),d(L), any performance improvement using c(L),d(L) instead of constant c,d for parameterized σ_HS and fSCA originates in introducing scale-dependent parameters. For the parameterizations using the constant fit parameters c,d, errors varied slightly more across scales than when using the scale-dependent c(L),d(L) version. Individual scale-dependent errors were in part larger than the MPEs for all L values given in Table 1. Unequal numbers of valid domains per L most likely also contributed to this.

5.4.3 Scale- and region-dependent evaluation for σ_HS and fSCA

While we did not perform an evaluation using independent snow depth data sets, studying regionwise performances reveals the spread in errors we can expect when the new parameterizations are applied on an individual independent data set (Fig. 11). We obtain much larger positive MPEs for σ_HS at lower spatial scales of L= 200 m and L= 300 m for the two TLS data sets in the southeastern French Alps and overall larger MPEs between 20 % and 30 %, though consistent across scales, for the Pléiades data from the Bassiès basin in the northeastern French Pyrenees. It is unclear if these larger MPEs originate in uncertainties of the data acquisition, i.e., are platform specific, or if they are linked to spatial snow depth distributions which could not be captured by the proposed new parameterizations. RMSEs for the various remote sensing platforms and data sets used here (Sect. 2) decrease from 80 cm for Pléiades data from the Sierra Nevada to 33 cm for the ADS, to 16 cm for UAS, to 13 cm for ALS data from Switzerland, to 8 cm for ALS data from the Sierra Nevada and to 4 to 10 cm for TLS data in general. Given the rather low errors typically obtained for TLS data compared to the other remote sensing platforms, the reason for the large deviations in the TLS data sets might not originate in inaccuracies of the data acquisition. On the contrary, the observed bias in the Pléiades data from the northeastern French Pyrenees might indeed be attributed to the rather large inaccuracies of the platform with NMADs of 45 to 78 cm (Marti et al., 2016). However, Pléiades data from the Sierra Nevada come with a similarly large NMAD of 69 cm, but σ_HS can be parameterized very well with MPEs lower than ±3 % across spatial scales (Fig. 11a). Observed σ_HS from the TLS, as well as from the Pléiades data in France, was considerably larger than parameterized σ_HS, but mean slope angles alone can also not explain this behavior (between 6 and 23^∘ for the TLS data and between 13 and 50^∘ for the Pléiades data).

While we were not able to clearly relate some of the poorer regionwise performances to uncertainties related to the platform, other studies entirely focused on performing extensive intercomparisons between platforms for large-scale snow depth mapping in alpine terrain (e.g., Bühler et al., 2015, 2016; Eberhard et al., 2021; Deschamps-Berger et al., 2020).

5.4.4 Evaluation of previous closed form parameterizations

Though we developed a new σ_HS parameterization (Eq. 3), empirically derived parameterizations can only describe the variability inherent in the data set used to derive the parameterization. In addition to the regionwise evaluation, analyzing performances of previous empirically derived parameterizations may therefore allow us to estimate expected performance sensitivity to independent data sets. While both tested parameterizations of σ_HS (Helbig et al., 2015; Skaugen and Melvold, 2019) showed a worse performance than the new parameterizations and less consistency as a function of scale, the model performances of Helbig et al. (2015) were only slightly worse than the new parameterizations (Table 1). The parameterization for σ_HS of Skaugen and Melvold (2019) was developed on mean domain sizes L of 750 m, whereas Helbig et al. (2015) used L values between 50 m to 3 km. This difference might be one reason for the overall poorer performances of Skaugen and Melvold (2019) compared to Helbig et al. (2015) across spatial scales. Since only 1 out of the 11 data sets used in this study was previously used to develop the parameterization of Helbig et al. (2015), an overall similar performance of Helbig et al. (2015) (Fig. 12) with the large compiled data set of this study clearly confirms the underlying functional form of σ_HS suggested by Helbig et al. (2015) which was reapplied here.

6 Conclusions

We presented an empirical peak of winter parameterization for the standard deviation of snow depth σ_HS for treeless, mountainous terrain, describing the spatial snow depth distribution in a grid cell for various model applications. The scaling variables in the new parameterization of σ_HS and fSCA are the same as in Helbig et al. (2015) which are spatial mean snow depth, a squared-slope-related parameter and a terrain correlation length. All subgrid terrain parameters can be easily derived from fine-scale summer DEMs for each coarse grid cell.

By introducing spatial scale dependencies in the variables in the formulation for σ_HS of Helbig et al. (2015), σ_HS can be consistently parameterized across spatial scales starting at scales ≥200 m. The spatial snow depth variability or σ_HS is the important variable to parameterize the fractional snow-covered area or fSCA (Helbig et al., 2015). Performance improvements across spatial scales of the σ_HS parameterization therefore directly enhanced the fSCA parameterization. Between length scales of 200 m and 5 km, mean percentage errors (MPE) were between −7 % and 3 % for σ_HS and between 0 % and 1 % for fSCA.

The subgrid parameterization of σ_HS was developed from 11 spatial snow depth data sets from seven different geographic regions at high spatial resolutions between 0.1 to 3 m and with spatial coverage between 0.14 to 280 km². An evaluation of two previously presented empirical σ_HS parameterizations confirmed the functional form of the parameterization of Helbig et al. (2015), as well as the need to enhance its performance across scales. By analyzing data from the large pool of spatial snow depth data sets, we were able to recalibrate the subgrid parameterization of σ_HS and achieve improved performances using new constant fit parameters. Additionally introducing a scale dependency in the dominant scaling variables further improved the performance across spatial scales. Mean MPEs of σ_HS over all scales (i.e., pooled performance) reduced from −30 % using Helbig et al. (2015) to −5 % after recalibration to −4 % after introducing scale-dependent fit parameters (Table 1). Individual scale-dependent improvements in MPEs reached up to 4 % when using newly derived scale-dependent fit parameters compared to newly derived constant fit parameters for σ_HS from the large data pool. This shows the improvement thanks to introducing scale-dependent parameters (Fig. 10). Towards estimating the possible spread in performances when applying empirically derived σ_HS and fSCA for independent geographic regions, we validated the parameterizations separately for each region and scale. While this clearly increased MPEs for three data sets, the majority of the region- and scale-dependent MPEs were between ±10 % for σ_HS and between −1 % and 1.5 % for fSCA, indicating that the parameterizations perform similarly well in most geographical regions.

By parameterizing peak of winter σ_HS in mountainous terrain, Helbig et al. (2015) extended the fSCA parameterization of Essery and Pomeroy (2004) for topography. Here, we extended peak of winter σ_HS parameterization, and thus the fSCA parameterization, to be applicable over a large range of snow climates and topographic characteristics, as well as across spatial scales. Since fSCA is an essential model parameter, a seasonal fSCA algorithm describing the variability throughout a winter season is of relevance for various model applications. A peak of winter parameterization of σ_HS that describes the maximum spatial snow depth variability can not, however, be used alone to describe that seasonal fSCA evolution. A reliable model application of any fSCA parameterization requires an implementation accounting for alternating snow accumulation and melt events during the season, i.e., to describe the SCD. Especially at lower elevations, the separation of the SCD in only one accumulation period followed by a melting period is no longer valid (Egli and Jonas, 2009). A description of an algorithm for a seasonal fSCA model implementation which uses the new scale-independent peak of winter parameterization of σ_HS in the fSCA parameterization presented here is currently in preparation. Extending the empirical peak of winter fSCA parameterization to a broader range of spatial scales and snow climates was thus a meaningful step towards accounting for spatiotemporal variability in snow depth for multiple snow model applications.