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

. The spatial distribution of snow in the mountains is signiﬁcantly inﬂuenced 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 con-ﬁrm 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 recalibrat-ing 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 sug-gests that the new parameterizations perform similarly well in most geographical regions.


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 simula-N. Helbig et al.: Fractional snow-covered area tions (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 fS-CAs 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 finescale 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 wa-ter 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: socalled 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 lognormal 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;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 vali- dated 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.

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 2 . 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 squaredslope-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.

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). ADSderived snow depth data sets were used from 20 March 2012 ("ads-CH 2 ") to 9 March 2016 ("ads-CH 1 "), together with summer digital elevation models (DEMs) (Marty et al., 2019). The data set covers about 150 km 2 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 me-dian 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 UASderived snow depth data set was used from 7 April 2018 ("uav-CH 9 ") from Schürlialp, together with a UAS-acquired summer DEM (Eberhard et al., 2021). The Schürlialp data set covers about 3.2 km 2 which we used at 30 cm resolution. A second UAS-derived snow depth data set was used from 29 March 2019 ("uav-CH 8 ") from Gaudergrat, together with a UAS-acquired summer DEM. The Gaudergrat data set covers about 0.8 km 2 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 UASderived 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 3 "), together with a summer DEM from 2017. The ALS data set from Switzerland used here covers about 260 km 2 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 2 forested sub-data sets. Validation of 1 m ALSderived 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.

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 7 ") and 2 May 2017 ("als-US 6 "), together with a summer DEM (Painter, 2018). The second platform was a Pléiades product from 1 May 2017 ("plei-US 6 "). 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 2 . 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 2 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.

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 4 "), together with a summer DEM (Marti et al., 2016). The data set we used covers about 113 km 2 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.

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 10 ") . A site and data description can be found in Naaim-Bouvet et al. (2010), Vionnet et al. (2014), and Schön et al. (2015Schön et al. ( , 2018. We used a UAS-acquired summer DEM (Guyomarc'h et al., 2019). The data set covers about 0.6 km 2 at 1 m resolution. The second TLS-derived snow depth data set was acquired over Col du Lautaret at 27 March 2018 ("tls-FR 5 ") . We used a TLS-acquired summer DEM. The data set covers about 0.14 km 2 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).

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.

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

Fractional snow-covered area parameterization
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): 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 µ = 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 ξ = √ 2σ z /µ with the standard deviation of elevations σ z . The L/ξ 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.

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/ 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/ξ -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.

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 dom- inant 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).
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).

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/ξ ratio (Eq. 2). Skaugen and Melvold (2019) used HS and the standard deviation of the squared slope, with sqS being derived using sqS = tan 2 ζ , where ζ is the slope angle in radians. Since tan 2 ζ is the same as 2µ 2 (e.g., Löwe and Helbig, 2012), we here derive sqS from 2µ 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/ξ and σ z .

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 twosample 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].

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

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/ξ 0.21 and for σ z 0.01. Though the mean correlation between σ HS and L/ξ 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). We selected HS, µ and L/ξ as main scaling parameters for σ HS across spatial scales from 200 m to 5 km (Fig. 8b).

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/ξ . 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.  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). 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: 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 pa-624 N. Helbig et al.: Fractional snow-covered area rameters 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.

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.

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.

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

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).  (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.

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;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 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)   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.

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 spa-tial 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 2 (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µ 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 spa-tial 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 .

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

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 scaledependent evaluation.

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

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 N. Helbig et al.: Fractional snow-covered area 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 largescale snow depth mapping in alpine terrain (e.g., Bühler et al., 2015Bühler et al., , 2016Eberhard et al., 2021;Deschamps-Berger et al., 2020).

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.

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 snowcovered 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 2 . 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 scaledependent 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.
Data availability. All data used in this study are described in the data section. The data can be downloaded from the referenced repositories, or data availability is described in the referenced publications.
Author contributions. YB, LE, CDB, SG, MD, JR, JSD and TJ were responsible for data acquisition and processing, experiment and platform design, and setup. NH was responsible for the development of the parameterization. All authors contributed to the paper with discussions and ideas. NH wrote the paper with contributions from all coauthors.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. We thank Vincent Vionnet, Yannick Deliot, Gilbert Guyomarc'h and Aymeric Richard for their effort in acquiring and processing TLS observations at Col du Lac Blanc. CNRM/CEN is part of Labex OSUG@2020 (investissement d'avenir -ANR10 LABX56). Pléiades imagery was acquired through DINAMIS (Dispositif Institutionnel National d'Approvisionnement Mutualisé en Imagerie Satellitaire). We further thank the two anonymous reviewers for the helpful comments on the manuscript. Financial support. Nora Helbig was funded by a grant of the Swiss National Science Foundation (SNF) (grant no. IZSEZ_186887), as well as partly funded by the Federal Office of the Environment (FOEN). Simon Gascoin was partly supported from CNES Tosca and Programme National de Télédétection Spatiale (PNTS) (grant no. PNTS-2018-4). Jesus Revuelto was funded by an AXA grant (grant no. CNRM 3.2.01/17) and by ANR JCJC EBON (grant no. ANR-16-CE01-0006).
Review statement. This paper was edited by Carrie Vuyovich and reviewed by two anonymous referees.