Fractional snow-covered area parameterization over complex topography

Fractional snow-covered area (SCA) is a key parameter in large-scale hydrological, meteorological and regional climate models. Since SCA affects albedos and surface energy balance fluxes, it is especially of interest over mountainous terrain where generally a reduced SCA is observed in large grid cells. Temporal and spatial snow distributions are, however, difficult to measure over complex topography. We therefore present a parameterization of SCA based on a new subgrid parameterization for the standard deviation of snow depth over complex topography. Highly resolved snow depth data at the peak of winter were used from two distinct climatic regions, in eastern Switzerland and in the Spanish Pyrenees. Topographic scaling parameters are derived assuming Gaussian slope characteristics. We use computationally cheap terrain parameters, namely, the correlation length of subgrid topographic features and the mean squared slope. A scale dependent analysis was performed by randomly aggregating the alpine catchments in domain sizes ranging from 50 m to 3 km. For the larger domain sizes, snow depth was predominantly normally distributed. Trends between terrain parameters and standard deviation of snow depth were similar for both climatic regions, allowing one to parameterize the standard deviation of snow depth based on terrain parameters. To make the parameterization widely applicable, we introduced the mean snow depth as a climate indicator. Assuming a normal snow distribution and spatially homogeneous melt, snow-cover depletion (SCD) curves were derived for a broad range of coefficients of variations. The most accurate closed form fit resembled an existing fractional SCA parameterization. By including the subgrid parameterization for the standard deviation of snow depth, we extended the fractional SCA parameterization for topographic influences. For all domain sizes we obtained errors lower than 10 % between measured and parameterized SCA.


Introduction
At the peak of winter, a snow cover resembles a sparkling, smooth blanket.However, it is well known that the spatial distribution of snow depths underneath is heterogeneous.Complex topography adds extra spatial variability due to spatial patterns of wind (sheltering/exposure), precipitation (e.g., mountain luv/lee), shortwave radiation (shading, sky view, terrain reflections) and longwave radiation (sky view, terrain emission).Furthermore, in complex topography, snow relocation can occur due to snow avalanches.To complicate matters, these processes operate at different spatial scales (cf.Liston, 2004).The result is a patchy snow cover consisting of snow-free and snow-covered areas (SCAs).In various scientific and operational applications, knowledge about spatial snow depths plays a key role.Hydrologists are interested in predicting the timing of snowmelt runoff as well as the overall amount of snow in a catchment to estimate the water stored, allowing one to forecast available water resources.This is a relevant issue, e.g., in controlling the drinking water supply, in hydropower production planning or in warning of spring flooding.Climatologists, studying present and future climates, are interested in the snow coverage in a largescale model grid cell which forms a key parameter in general circulation models (e.g., Roesch et al., 2001).For instance, from fractional SCAs, coarse-scale surface albedos can be derived by weighting snow-free and snow-covered albedos (Liston and Hiemstra, 2011).Since snow has a high surface albedo, it alters the energy and moisture fluxes on the earth and thus the surface energy budget (Dingman, 1994).Knowing the actual spatial snow depth distribution, especially in mountainous terrain, is therefore a relevant topic in large-scale hydrological, meteorological and regional climate models.Due to computational constraints, large-scale models often have to simplify physical processes over snow surfaces and within snow.Frequently, they lack a subgrid snow distribution representation which is a shortcoming that deteriorates atmospheric interaction simulations (cf.Liston et al., 1999).In general, the purpose of subgrid parameterizations is to account for subgrid scale processes, i.e., unresolved processes, with analytical approximations in largescale model systems.The IPCC (2007) considered subgrid snow distributions as important for simulating observations of seasonal snow cover.
A few studies previously tackled subgrid snow distributions.Liston (2004) improved a regional climate model by performing separate surface energy balance calculations over snow-covered and snow-free fractions of each model grid cell.Similar, Ménard et al. (2014) calculated vertical and horizontal energy fluxes between the atmosphere and snow, snow-free and vegetation grid cell portions and found a warming feedback through decreases in surface albedo and increases in sensible heat fluxes to the atmosphere.Liston (2004) computed SCAs by assuming lognormally distributed snow depth and by introducing a dichotomous key for coefficient of variations for snow depth (CV is standard deviation divided by mean) depending on topographic variability, air temperature and wind speed.Liston and Hiemstra (2011) introduced a snow cover protruding vegetation fraction for grid cell portions covered by shrubs or grass.Essery and Pomeroy (2004) validated previously published ad hoc closed forms of SCA over non-forested terrain with those derived from a peak of winter lognormal distribution that undergoes homogeneous melt.They found the closest snow-cover depletion (SCD) curves using a functional form proportional to tanh, similar to what was proposed by Yang et al. (1997) and Roesch et al. (2001).Instead of a roughness length of the surface (Yang et al., 1997) or the standard deviation of the summer digital surface model (DSM) (Roesch et al., 2001), Essery and Pomeroy (2004) included the peak of winter standard deviation of snow depth in the SCA parameterization.However, peak of winter standard deviations of snow depth are rarely available.
Numerous studies analyzed catchment snow depth distributions by relating measured snow depth data to small-scale terrain parameters (for a recent literature overview see Clark et al., 2011).Until now, multiple linear regressions were frequently applied to relate mean snow depth, standard deviation of snow depth or deviations of the mean to smallscale terrain parameters such as elevation, slope or aspect.Others found linear (Pomeroy et al., 2004) or power-law (Egli and Jonas, 2009;Egli et al., 2011) relationships for the accumulation period, solely between standard deviation of snow depth and mean snow depth using constant fit param-eters.While the CVs presented by Liston (2004) depend on topographic variability, the relationships of Pomeroy et al. (2004), Egli and Jonas (2009) and Egli et al. (2011) result in CVs which neglect varying complexities of terrain.Even though previous parameterizations for the snow distribution parameters provide good descriptions for the investigated regions, they might easily fail in a different geographic region with other terrain characteristics.Recently, Grünewald et al. (2013) analyzed snow depth data from seven mountainous catchments around the world.For each catchment, their developed multiple regression equations for the relative snow depth (HS -catchment mean) using subgrid topographic parameters showed good performance.However, a similar performance for a global model, based on all data sets, could not be achieved, and Grünewald et al. (2013) argue that the snow depth and topography are less universally related than hypothesized by Lehning et al. (2011).
A poorer performance of a subgrid parameterization for the snow distribution can also arise from the different scales on which the spatial variability of snow depths is created in complex topography.Recently, Grünewald et al. (2013) and Melvold and Skaugen (2013) therefore investigated the influence of scale on aggregated snow depth data.By analyzing snow depth data in differently sized grid cells up to 800 m for several catchments, Grünewald et al. (2013) found a lower limit of 400 m for the grid cell size to explain most of the remaining larger-scale spatial variability.By analyzing snow depth data from a large mountainous area in Norway in grid cell sizes up to 1 km, Melvold and Skaugen (2013), however, determined a larger lower limit of 1 km to eliminate most of the spatial variability such that the mean adequately represents the average grid cell snow depth.A reason why a global parameterization might not be derivable at one certain horizontal resolution is that too many different snow-cover shaping processes are still active, at that scale, making it a challenge to parameterize the subgrid snow distribution.
How can we acquire snow depth data spatially in order to better investigate subgrid snow depth distributions?Measuring snow distribution, both temporally and spatially, is a challenging task in mountainous terrain.To overcome the limitations of point measurements of automated stations or hand probing, terrestrial laser scanning (TLS) was introduced to continuously measure snow depths in very high resolutions (Prokop et al., 2008;Grünewald et al., 2010).Airborne laser scanning (ALS) can cover larger regions in a shorter time without the limitations of TLS (Hopkinson et al., 2004;Deems et al., 2006;Grünewald and Lehning, 2011).ALS measurements are, however, quite expensive and for larger regions they require large investments to gather snow depths in adequate temporal and spatial resolutions (e.g., NASA's Airborne Snow Observatory, 2013).Visible satellite remote sensing provides information on snow coverage in various horizontal and temporal resolutions.However, the interpretation of satellite signals is difficult and requires complex algorithms extracting clouds (e.g., Hüsler et al., 2012) and the influence of topography on the signal (e.g., Stöckli, 2013).Small-scale distributed snow surface modeling (e.g., Lehning et al., 2006) over complex topography could fill the gap of missing temporal and spatial snow depth data.However, for large regions this is rarely feasible due to computational constraints and/or the lack of small-scale input data.Erroneous input data could easily blur modeled distributed snow depths.For now, we therefore prefer spatially and temporally measured snow depth data to investigate subgrid snow depth distributions.
To our knowledge, a systematic analysis of snow depth data from a large region, aggregated in grid sizes comparable to those of large-scale models, is still missing.Here, we are aiming for grid cell sizes where the subgrid variability is deducible from the underlying characteristic terrain lengths.We assume that the smoothing out of smallscale snow depth heterogeneities originating from processes such as snowdrift or avalanches reveals the large-scale topographic influences on precipitation and the shortwave radiation balance.Our hypothesis is motivated by the observation of Liston et al. (1999), in that, in contrast to summer convective-precipitation systems, the spatial distribution of winter precipitation is more influenced by topographic distributions.Furthermore, it is motivated by the results of Grünewald et al. (2013) and Melvold and Skaugen (2013), which confirmed that the snow depth distribution is dominated by topography at scales of several hundred meters.
In this study our principal goal is thus to develop a subgrid parameterization of SCA for large-scale model grid cell sizes of a few kilometers that account for varying levels of complex, treeless topography.For this, we relate snow depth data to terrain parameters in view of a subgrid parameterization of the standard deviation of snow depth.We use easily accessible, computationally cheap terrain parameters calculated from the summer DSM.We employ highly resolved spatial snow depth data from alpine terrain of two large areas in the eastern Swiss Alps as well as from one in the eastern part of the Spanish Pyrenees, i.e., from two distinct climates.The snow depth data resolves for all small-scale variability of the snow cover.We analyze the probability density functions (pdf) of snow depth and the two defining parameters, mean and standard deviation, and examine the data both within and between domain sizes of various dimensions.Finally, we point out the limitations of our subgrid parameterizations originating from using measured snow depth data sets.

Site descriptions
To account for the influence of different climates on the spatial snow distribution, we used snow depth data from three large alpine areas in two distant geographical regions.Two alpine areas, called Wannengrat and Dischma, are located in eastern Switzerland around Davos (Fig. 1a).Wannengrat covers about 30 km 2 and Dischma about 120 km 2 .In the Wannengrat area, elevations range from 1517 to 2781 m and in the Dischma area elevations range from 1516 to 3227 m.The mean slope angle, which was computed from 2 m elevation differences, is 26 • for Wannengrat and 28 • for Dischma.
The third alpine catchment, called Val de Núria, is located in the eastern part of the Spanish Pyrenees (Fig. 1b) showing a dryer snow climate.Val de Núria covers about 28 km 2 of treeless mountainous terrain (Moreno Banos et al., 2009).Elevations range from 1910 to 2910 m for Val de Núria.The mean slope angle, which was computed from 1 m elevation differences, is 24 • .

Digital photogrammetry
For the Wannengrat and Dischma sites, spatial snow depth data were obtained using an opto-electronic line scanner (Sensor ADS80, Leica Geosystems) mounted on a plane.Photogrammetric image correlation techniques were applied for summer and winter aerial imagery to calculate DSMs in 2 m horizontal resolution (Bühler et al., 2015).Spatial snow depths were obtained by subtracting the summer from the winter DSM.The winter DSM of the Wannengrat area (cf.Fig. 1a) shows a root mean square error (RMSE) of approximately 33 cm with snow depths from simultaneously conducted TLS measurements and a RMSE of approximately 19 cm with snow depths from snow probing in plots with 5 by 5 probes per plot (Bühler et al., 2015).The winter DSM of the Dischma area (cf.Fig. 1a) shows a larger RMSE of approximately 43 cm with snow depths obtained from ground penetrating radar measurements at the valley bottom.The snow depth data sets were acquired at approximate peak of winter on 20 March 2012.The mean snow depth at Wannengrat was 1.72 and 2.07 m at the Dischma area.

Airborne laser scanning
For the Val de Núria site, point clouds of snow depth values were obtained by ALS measurements (Moreno Banos et al., 2009).Based on this data, Grünewald et al. (2013) calculated summer and winter DSM in 1 m horizontal resolution, which they then subtracted to obtain the spatial snow depth data.The mean accuracy in vertical direction is 30 cm which is similar to the ADS80 data.The ALS campaign took place at the approximate peak of winter on 9 March 2009.The mean snow depth at Val de Núria was 1.07 m.

Preprocessing
For Wannengrat and Dischma, we neglected all measurements that coincided with trees, buildings, rivers and glaciers.Negative snow depth values were set to zero.In total we obtained about 6 × 10 6 usable snow depth measurements for Wannengrat and about 22×10 6 for Dischma.The data set of Val de Núria was preprocessed as described in Grünewald et al. (2013) resulting in about 28 × 10 6 usable snow depth measurements.Figure 2 shows the pdf of all measured snow depths for the three areas.

Aggregating snow depth data
Analyzing a sufficiently large number of differently sized domains from a large mountainous region allows one to study snow distributions at different scales.By randomly selecting different grid origins, we aggregated the snow depth data sets in different squared domain sizes L. Note that L can be seen as a coarse grid cell size x in a large-scale model (cf.Fig. 1a).We chose domain sizes of L = 50, 100,200,500,750,1000,1250,1500,1750,2000,2500 and 3000 m covering the range of typical grid cell sizes from hydrologic measurement campaigns to the smallest grid cell sizes in meteorological models.For each domain size we used 50 realizations allowing for overlap between domain sizes L (cf.Fig. 1a).In total we generated ensembles of 600 snow depth grids for each Swiss site.In Val de Núria we could not aggregate snow depth data in domain sizes L larger than 1500 m, resulting in 400 snow depth grids at this site.
For building domain averages, all data points were spatially averaged in a domain size L.However, we only used domain sizes L with at least 75 % valid snow depth measurements (including zero values).For larger domain sizes L ≥ 1 km in Val de Núria we had to allow for a maximum of 40 % of missing values due to the irregular perimeter of that catchment (cf.Fig 1b).In the following, we omit the normally used overbars for domain-averaged variables.

Terrain characteristics
To relate the snow depth distribution parameters to topographic features, we computed several terrain parameters from the summer DSMs.For selecting terrain parameters, we exploited the fact that real topographic slope characteristics are reasonably well described by Gaussian statistics (Helbig and Löwe, 2012).Gaussian random fields with a Gaussian covariance such that topography is reduced to only two underlying large length scales in a model domain of size L, were previously used to systematically investigate radiative transfer in complex terrain via the radiosity approach (Helbig et al., 2009;Helbig and Löwe, 2012;Löwe and Helbig, 2012) as well as to develop a parameterization for domain-averaged sky view factors in complex terrain (Helbig and Löwe, 2014).
Assuming a Gaussian covariance for the summer topography, the two underlying characteristic length scales are a valley-to-peak elevation difference σ (typical height of topographic features), which is the standard deviation of the elevation model, and a lateral extension ξ (typical width of topographic features), which is the correlation length of the elevation model.We use a terrain parameter µ = √ 2σ/ξ , which is related to the mean squared slope and which can be derived from first partial derivatives ∂ x z and ∂ y z (slope components) in orthogonal directions: ) 2 as outlined by Löwe and Helbig (2012).We also use the L/ξ ratio where a large ratio indicates that more topographic features are included in a domain size L. Note that the typical width of topographic features ξ in a domain size L can be obtained via ξ = √ 2σ z /µ, with the standard deviation of the summer DSM σ z .Helbig et al. (2009) showed that to minimize influences of (subgrid) grid size x and domain size L on domain-averaged shortwave terrain reflected radiation, the condition x ξ L must be fulfilled.The relevance of including enough terrain in a domain, here L × L, was confirmed by Helbig and Löwe (2014), where errors of a subgrid parameterization for the sky view factor over complex topography decreased with increasing L/ξ ratio.We believe that in complex terrain for domain-averaged snow depths, the above condition should always be met in order to accurately capture the predominant subgrid processes shaping the snow distribution at the corresponding scale.Consequently, we need to detrend the summer DSMs in order to obtain the correct characteristic length scales for the corresponding domain size L. Linearly detrending reveals the dominant processes that shape the scale dependent characteristic snow depth distribution by shifting the scaling parameters.For small domain sizes L this leads to smaller correlation lengths ξ and thus to larger L/ξ ratios.

Parameterizing spatial variability of snow depth
In order to specify the spatial variability of snow depth over mountainous, treeless topography for large-scale grid cells, we first need to define the pdf of snow depths in a domain size L. Commonly applied snow depth distributions at the peak of winter range from lognormal for complete snow cover (Donald et al., 1995;Pomeroy et al., 1998;DeBeer and Pomeroy, 2009) to gamma (Skaugen, 2007;Egli et al., 2012) to normal in forests (Marchand and Killingtveit, 2005).Second, we need to scale the defining parameters mean and standard deviation of the snow depth distribution, HS and σ HS , respectively, with the underlying subgrid terrain characteristics.Previously published linear (Pomeroy et al., 2004) or power-law (Egli and Jonas, 2009) relationships, solely be- tween σ HS and HS, lead to snow depth coefficients of variation CV which do not depend on varying topography.Yet, we computed a mean CV for L ≥ 1 km of 0.63 for the Wannengrat and 0.48 for the Dischma region.The CV for the catchment in the eastern Spanish Pyrenees for L ≥ 1 km is 1.04, i.e., considerably larger than for the two large areas in the eastern Swiss Alps.Deriving the CVs from the powerlaw relationship (via σ HS = HS 0.84 ) results in overall larger but similar CV values among the three regions: 0.91 for Wannengrat, 0.89 for Dischma and 1.01 for Val de Núria.The CV of the eastern Swiss Alps compares well to the CV categories of the dichotomous key in that geographic region of 0.5 to 0.7, which was based on topographic variability, air temperature and wind speed (Liston, 2004).However, for the area in the eastern Spanish Pyrenees the CV of the dichotomous key of Liston ( 2004) is about 0.06, i.e., completely different to our 1.04.Given that we use snow depth data sets from two distinct climate regions, we can focus on the development of a subgrid parameterization of the standard deviation of snow depth σ HS which is not constrained to one specific geographic area but is more widely applicable.For this, we employ the mean snow depth HS as a climate indicator variable for each domain size L.However, mean HS is generally not easily measured.We therefore investigate if mean snow depth HS can be approximated by averaged flat field measurements HS flat .A flat field was defined as a 22 m × 22 m (for Wannengrat and Dischma) or a 11 m × 11 m (for Val de Núria) area where each slope angle was lower than or equal to 10 • .We computed the average flat field snow depth from all snow depth values within a flat field.To obtain an average flat field snow depth HS flat for each domain size L, we averaged all mean snow depths of flat fields within each L. Note that in the following we will use the superscript m for measured, mean quantities when opposed to parameterized quantities.

Snow depth distribution
We found mostly unimodal distributions of snow depths in all domain sizes L ranging from 50 m to 3 km in all three areas (Fig. 3).We tested three, previously published theoretical pdfs on our ensembles of gridded snow depth data: normal, lognormal and gamma density functions.While for small domain sizes a gamma distribution best described the measured snow depth distributions, for larger domain sizes (L ≥ 500 m) a normal distribution worked as well or better (Fig. 4).The mean RMSE between theoretical pdfs and measured snow depths decreased with increasing domain size L for all three areas.A comparison of computed quantiles for the theoretical and measured snow depth distribution also resulted in decreasing mean RMSE with increasing L. Note that our domain sizes do include subgrid snow-free values.

Scaling of snow depth data grids
We analyzed our ensemble of snow depth data grids to relate mean and standard deviation of each snow depth distribution, HS and σ HS , to terrain parameters.An interesting result is that the mean of σ HS increased with increasing L. For domain sizes of L ≥ 1 km the overall changes in the mean of σ HS became small (Fig. 5).Similar to Grünewald et al. (2013) and Melvold and Skaugen (2013) we found that overall, with larger domain size L, the scatter in standard deviation of snow depth σ HS decreased (Fig. 5).However, in comparison to Grünewald et al. (2013) and to Melvold and Skaugen (2013), we also included L > 1 km and found that for L ≥ 1 km the scatter in σ HS still somewhat decreased.Note that we obtained similar trends and magnitudes of σ HS as a function of domain size L for both climates, which allowed us to pool the data of all three areas.Furthermore, similar trends in σ HS were found with terrain parameters in all three areas, suggesting that a parameterization can be developed which can be applied under a broad range of topographic characteristics.For example, Fig. 6 shows the standard deviation of snow depth σ HS of the three areas as function of the standard deviation of the summer DSM, σ z .In all areas σ HS increased similarly with increasing σ z and with increasing domain size L. Furthermore, the scatter or the standard deviation of σ HS among the same domain sizes L decreased with increasing L and σ HS .A correlation analysis between terrain characteristics and standard deviation of snow depth σ HS revealed significant Pearson correlation values ranging from 0.22 to 0.65 for pooled snow depth data from all catchments (Table 1).The overall larger scatter in snow depths for all L in the Dischma catchment (cf.Figs. 5 and 6) resulted in  lower correlation values r when looking at the correlations coefficients of each area separately (cf.Table 1).We found weaker correlations between mean snow depth HS and terrain parameters, than between σ HS and terrain parameters (Table 1).For the correlation between terrain parameters and pooled snow depth data from all catchments, the significance was marginally lower than for σ HS .However, the correlation analyses between HS and terrain parameters conducted for each catchment separately often showed statistically insignificant correlations, i.e., p values ≥ 0.05 (Table 1).Yet, we observed an approximately linear relationship between HS and mean flat field snow depths HS flat when we pooled snow depth data of all areas, especially for domain sizes larger than 1500 m (Fig. 7).The overall deviations between HS and HS flat decreased with increasing domain size L.For the overall relationship of HS and HS flat , we obtained a Pearson correlation coefficient r of 0.86, a squared correlation coefficient R 2 of 0.65, a RMSE of 36.7 cm, a normalized root mean square error NRMSE of 5.4 % and a mean squared error (MSE) of 13.4 cm.

Parameterization for the standard deviation of snow depth
In order to develop a parameterization for σ HS , we pooled the snow depth data of all three areas.We derived the following subgrid parameterization for the standard deviation of snow depth σ HS over mountainous terrain from snow depth data aggregated in domain sizes ranging from L = 50 m to 3 km: with a = 0.549 and b = 0.309 and HS, ξ and L in meters.
When fitting for each area separately, the parameters changed slightly.The standard deviation of snow depth σ HS in Eq. ( 2) has three scaling parameters: a terrain parameter µ (Eq.1), related to the mean squared slope in each domain size L, the mean snow depth HS and the L/ξ ratio, roughly describing how many subgrid topographic features are in a domain size L. The functional form of our subgrid parameterization was motivated by the result that we consistently obtained the largest correlation coefficients for σ HS with the terrain parameter µ (cf.sizes L with varying correlation lengths of topographic features ξ the condition L/ξ 1 is not always fulfilled and corrections are required.Naturally, the correction factor decreases with increasing L/ξ ratio.We chose a Gaussian factor e −(ξ/L) 2 based on our result that in large-scale grid sizes the snow depth distribution can be described by a Gaussian distribution.Assuming that topography is the major driver for the snow distribution, the Gaussian factor is also a consequence of previously found Gaussian slope statistics for real topographies (Helbig and Löwe, 2012).Mean snow depth HS has to be included in a parameterization of σ HS to account for varying surface climates.We performed the nonlinear regression analysis to optimize the parameters in Eq. ( 2) by robust M-estimators using iterated reweighted least squares; see R v2.15.2 statistical programming language (R Core Team, 2012) and its robustbase v0.9-7 package (Rousseeuw et al., 2012).Our subgrid parameterization, as in Eq. ( 2), predicts the observed σ HS well (cf.Fig. 8a).The performance of the parameterization improves with increasing domain sizes L. Our subgrid parameterization for the standard deviation of snow depth σ HS is statistically significant (Pearson r = 0.70, p value < 0.001, R 2 of 0.45, RMSE of 22.9 cm, NRMSE of 7.6 % and MSE of 5.2 cm).The performance of parameterized σ HS (Eq.2) also improved compared to previously published parameterizations of σ HS , which did not explicitly account for subgrid topography (Fig. 8b and c).Note, that the subgrid parameterization for σ HS was developed for peak of winter snow depth data.

Parameterization of fractional snow-covered area
Snow-covered area is an important parameter in the energy balance of large-scale models, e.g., to weight energy flux components and surface albedos for snow-covered and snowfree fractions.Fractional SCA f in a large-scale grid cell is, however, reduced due to subgrid topographic effects on the snow depth distribution.Here, we showed that the standard deviation of snow depth σ HS at the peak of winter over complex topography scales with the underlying terrain characteristics combining previously published observations.We therefore suggest including σ HS , as in Eq. ( 2), in a closed form parameterization of the fractional SCA f .When deriving a functional form for f , Essery and Pomeroy (2004) concentrated on homogeneous surface units where the peak of winter snow depth distribution could be described by a lognormal distribution.We are focussing on large-scale grid cell sizes over complex topography where we employ our result that the simpler normal distribution describes the snow depth distribution equally well or better (cf.Fig. 4).We start the derivation from a normal distribution at the peak of winter over alpine terrain (including snow-free sub-pixels):  2), (b) parameterized via Egli and Jonas (2009) and (c) parameterized via Pomeroy et al. (2004).Colors indicate corresponding domain size L. NRMSEs are given for each parameterization.
with σ HS 0 as the standard deviation of snow depth and HS 0 as the mean snow depth at the peak of winter, both indicated here with the subscript 0. The SCA f is obtained by assuming a homogeneous melt amount M and by integrating over the peak of winter snow depth distribution from M to ∞: The mean snow depth HS is obtained from leading to We followed the procedure of Essery and Pomeroy (2004) to derive a more practical closed form of f than Eq. ( 4).For this we also assumed homogeneous melt for our peak of winter normal snow depth distribution (Eq.3).Note that Egli and Jonas (2009) showed that the concept of spatially uniform melt can even be applied over mountainous terrain when starting from a measured snow distribution.In contrast to Essery and Pomeroy (2004) we included a larger range for coefficients of variations CV to derive a closed form of f .We chose the CV values of Liston (2004) defining snow distribution categories around the world but added a maximum CV value of 1: 0.06, 0.09, 0.12, 0.17, 0.4, 0.5, 0.6, 0.7, 0.85, 1. Dashed lines in Fig. 9 show the fitted f (HS) by means of f (HS) = tanh 1.30 HS σ HS 0 (7) using σ HS 0 , the standard deviation of snow depth at the peak of winter.We obtain the closest fit with the same functional form as Essery and Pomeroy (2004) who started with a lognormal snow depth distribution (Fig. 9).Our pre-factor in Eq. ( 7) varies slightly from the one presented by Essery and Pomeroy (2004).For our data and the fit parameter we computed a 95 % confidence interval ranging from 1.27 to 1.35.For the fit in Eq. ( 7) we obtain a mean RMSE of 0.02, and a mean NRMSE of 2.5 % for all CV.Similar to the fit of Essery and Pomeroy (2004) our RMSEs increase with increasing CV with the largest RMSE of 0.04 for a CV = 1.We extend the fractional SCA f (HS) of Eq. ( 7) to complex topography by employing standard deviation of snow depth at the peak of winter parameterized for complex subgrid topography (cf.Eq. 2). Figure 10a shows that the mean errors between parameterized and observed SCA f for all our areas decrease with increasing domain size L. Also, the scatter per L decreases with increasing L (cf. error bars in Fig. 10a).Note that the largest mean errors are still below 10 %.When using previously derived parameterizations for σ HS in parameterized f (HS) (Eq.7) both mean errors and scatter also decrease with increasing L, however, the overall errors are larger and mean errors do not approach zero for the largest domain sizes L ≥ 1750 m (Fig. 10b and c).

Discussion and Conclusion
Scaling snow depth distribution parameters is a relevant issue for various applications in large-scale hydrological, meteorological and regional climate models.In this study, we derived a parameterization for the fractional SCA over complex, treeless topography for large-scale models with grid cell sizes of a few kilometers.This required developing a subgrid parameterization for the standard deviation of snow depth over mountainous terrain.For the parameterization we chose easy to derive subgrid terrain parameters and the mean snow depth as a climate indicator variable.We derived the subgrid parameterization from highly resolved snow depth data sets in large areas gathered at the peak of winter.
Investigating a spatial distribution entails studying the distribution parameters, mean and standard deviation.Furthermore, measured mean and standard deviation of snow depths require to be analyzed as a function of scale in order to reveal the scale at which the dominant shaping processes can be reliably parameterized, i.e., when small-scale snow depth variations no longer resolved.We performed a scale dependent analysis by creating data sets from randomly selecting differently sized squared domain sizes L (equivalent to a coarse grid cell size of a large-scale model) ranging from  4) as function of mean HS normalized with the peak of winter mean snow depth HS 0 (indicated by the subscript 0) (solid lines).Dashed lines represent parameterized fractional SCA f via Eq.( 7).Coefficient of variations CV vary between 0.06 (upper left) and 1 (lowest one).
Figure 10.Error in fractional SCA f between measured f m and parameterized f (Eq.7) as function of the domain size L for all three areas.(a) Parameterized using σ HS from Eq. ( 2), (b) parameterized using σ HS from Egli and Jonas (2009) and (c) parameterized using σ HS from Pomeroy et al. (2004).Mean values are indicated by squares.Error bars show the standard deviation of the error per L. NRMSEs are given for each parameterization.
better performance (Fig. 4).We therefore conclude that over alpine terrain, in large-scale grid cells, the snow depth distribution can be approximated by a simple normal distribution.We also found a strong dependency of the distribution parameter, the standard deviation of snow depth σ HS , as a function of coarse grid cell size L for domain sizes L ≤ 1 km in each of the three data sets separately (Fig. 5).This indicated that there should be a lower limit for large-scale grid cell sizes to minimize scatter, which we suggest being ≥ 1 km, similar to Melvold and Skaugen (2013).A scale analysis of domainaveraged snow depth values with domain-averaged terrain parameters revealed similar trends and magnitudes with terrain characteristics for each of the three catchments (Fig. 6 and Table 1).Scattering within a domain size L consequently decreased with increasing L. We therefore concluded that a parameterization using terrain parameters can possibly predict the standard deviation of snow depth at the peak of winter.Furthermore, this allowed us to create a pooled data set in order to derive a subgrid parameterization independently of one geographic region.Note that, despite similar trends between the three catchments, the scatter in the standard deviation of snow depth σ HS varied, for which we assume two reasons.First, there was increased overlap of the randomly picked domains in the smaller catchments of Wannengrat area and Val de Núria (cf.Fig. 1).Second, an overall larger scatter in the Dischma area data set might stem from a larger flight height resulting in higher measurement uncertainties which was, however, necessary due to local topographic features.
We developed a subgrid parameterization of snow depth distributions based on spatial snow depth data sets acquired by aerial imagery and photogrammetric image correlation techniques.Even though measurement errors can reach up to 33 cm (cf.Bühler et al., 2015) compared to small-scale modelings of spatial snow depths, which require detailed input data and which sometimes even rely on parameterizations, errors are clearly defined.Three snow depth data sets from large, alpine areas were analyzed to develop the subgrid parameterization of snow depth distributions.Two areas were located in eastern Switzerland and one catchment in the eastern part of the Spanish Pyrenees showing a somewhat dryer snow climate than the other areas (Fig. 1).We focussed on developing a subgrid parameterization for the standard deviation of snow depth σ HS independent of one specific geographic area or winter season.For this, we introduced the mean snow depth as a climate indicator variable.By analyzing flat field snow depth measurements, gathered at the peak of winter, we found that the mean of all average flat field snow depth measurements in a domain size L was approximately linearly correlated with the mean snow depth in the same L (Fig. 7).This was especially true for domain sizes L larger than about 1.5 km.It also has the interesting practical advantage that deriving the mean snow depth for a large domain at the peak of winter can be conducted by measuring snow depths on several flat field sites which are representa-tive for a specific geographic region.Since measuring snow depths on a few flat fields within each domain size for largescale models covering a wide area is generally not feasible, we suggest that those can be replaced by an automated flat field measurement, showing good climate representativeness for the corresponding large-scale domain size L. Though the linear relationship might have to be further verified in other geographic catchment areas and during other seasons, using measured flat field snow depths as an easily accessible climate descriptor allows one to develop a parameterization for the standard deviation of snow depth independently of its geographical region.The three snow depth data sets were gathered in two different winters, each time at approximately the peak of winter (Fig. 2).Until now we do not have measurements during other seasons and a re-evaluation of the subgrid parameterization for the standard deviation of snow depth σ HS (Eq.2) during other seasons might be necessary.However, in principle, using the mean snow depth as a climate indicator variable, Eq. ( 2) should also capture seasonal differences.
To relate snow depth distributions, measured at the peak of winter, to terrain characteristics we chose Gaussian statistics to approximate slope characteristics of real summer topographies.Assuming that real topographies can be described by a Gaussian covariance (cf.Helbig and Löwe, 2012) topography is reduced to two underlying characteristic length scales, namely, a typical height of topographic features σ (standard deviation of the summer DSM σ z ) and a typical width of topographic features ξ .From these we computed the L/ξ ratio indicating how many topographic features are included in a domain size L as well as a terrain parameter µ, which is related to the mean squared slope (Eq.1).Before deriving the terrain parameters we linearly detrended the summer DSM to reveal the correct characteristic terrain length associated with the shaping process of the snow depth distribution at the corresponding scale.Detrending all summer DSMs then resulted in reasonably large L/ξ ratios ranging from 2.7 to 15 for all domain sizes L. Without detrending overall smaller L/ξ ratios prevailed, with the smallest L/ξ ratio of 1.7 for L = 50 m.In grid cells with small ratios the relevant shaping processes might not be accurately resolved, and a subgrid parameterization could be flawed.Domain averages were built by spatially averaging the data in a domain size.
Overall, our subgrid parameterization for the standard deviation of snow depth σ HS (Eq.2) describes measured snow distributions in the three different alpine areas very well (Fig. 8a).As expected, the accuracy of parameterized σ HS increased with increasing domain size L (Fig. 8a).This is partly because at small scales the shaping processes are more diverse (or random) which are, however, smoothed out at larger scales, here for L ≥ 1 km.The parameterization in Eq. ( 2) describes the processes dominating at larger scales.On the other hand, the accuracy of the subgrid parameterization of σ HS also increases with increasing L/ξ ratios, i.e., the subgrid topographic features and their impact on snow depth distributions are represented more accurately.In the following we discuss the three scaling parameters in the subgrid parameterization of σ HS .First, it includes a terrain parameter, related to the mean squared slope µ (Eq.1), describing the influence of topography due to varying incident shortwave radiation and precipitation.This terrain parameter was motivated by the result that we consistently obtained the largest correlation coefficients for σ HS with µ (cf.Table 1).For now, we assume that parameterized σ HS approaches zero for mean squared slopes µ of zero.Even though, mean slope angles of all domains range from 2 to 58 • , the lowest domainaveraged slope angles only coincide with the smallest domain sizes.Equation ( 2) can be extended for large-scale grid cells showing slopes of zero, once the necessary snow depth data become available.As for the second scaling parameter, the parameterization includes the L/ξ ratio, a correction term for finite grid sizes which can show a range of correlation lengths of subgrid topographic features ξ (cf.Helbig and Löwe, 2014) that might or might not be captured by the domain size L. As a consequence of the overall good agreement of the pdf of snow depths with a normal distribution at larger scales we used a Gaussian factor e −(ξ/L) 2 .The Gaussian factor also follows from the assumption that topography has a large impact on the snow distribution in large-scale grid cell sizes and from the previously found Gaussian slope characteristics for slope characteristics of real topographies (cf.Helbig and Löwe, 2012).The third parameter in the σ HS parameterization includes the mean snow depth which accounts for climate or seasonal differences.
Since the snow depth data sets were only acquired at approximately the peak of winter slight hysteresis phenomena of the alpine, seasonal snow depth distribution (Egli and Jonas, 2009) were introduced (cf.Fig. 8a).With snow depth data gathered exactly at the peak of winter, constant parameters a and b in Eq. ( 2) might change but overall errors are expected to decrease.Note that we optimized a and b in Eq. ( 2) with a nonlinear regression analysis.Our parameterization performed better than previously published parameterizations for σ HS , which did not account for subgrid topography (Fig. 8b and c).Since the averaged coefficient of variation for snow depth CV of all domain sizes in our catchments of 0.63 resembles the one for alpine tundra of 0.43 which Pomeroy et al. (2004) used in a linear relationship, this parameterization shows an overall better performance among the tested parameterizations (Fig. 8).
By employing the new subgrid parameterization for the standard deviation of snow depth σ HS we developed a parameterization for the fractional SCA over complex mountainous terrain (Eq.7).For this large-scale model application we reevaluated a previously presented functional closed form for homogeneous landscapes (Essery and Pomeroy, 2004).To obtain a parametric fractional SCA f (Eq.4) we similarly integrated the snow distribution assuming uniform melt but started from a normal snow depth distribution (Fig. 4).Fitting the resultant parametric f (Eq.4) we obtained the same functional closed form fit as Essery and Pomeroy (2004) which is proportional to tanh (Eq.7).We assume that the slightly differing pre-factor stems from our broader range for CV stretching from 0.06 to 1 compared to the one used by Essery and Pomeroy (2004) with CV values from 0.1 to 0.5.We stress that the parameterization for σ HS (Eq.2) as a function of terrain characteristics coincides well with previously presented dependencies of f on terrain parameters such as the roughness length of the surface (Yang et al., 1997) or the standard deviation of the summer DSM (Roesch et al., 2001).Overall, we found decreasing errors between parameterized and measured f , for our three areas at the peak of winter, with increasing domain size L with the largest errors being below 10 % (Fig. 10a).When applying previously derived parameterizations for σ HS we also found decreasing errors between parameterized and measured f with increasing L. However, we obtained overall larger errors and errors did not approach zero for the largest domain sizes L ≥ 1750 m (Fig. 10b and c).We emphasize that applying Eq. ( 7) with parameterized σ HS leads to a NRMSE of only 4 % more than when applying measured σ m HS in Eq. ( 7).Note that in line with replacing exhausting snow depth measurements in large domain sizes L by parameterized σ HS via Eq.( 2), we investigated the increase of error in the parameterization for f when applying averaged flat field snow depths instead of mean snow depth HS.Applying measured HS flat instead of mean snow depth HS, but using measured snow depth distribution σ m HS , in the SCA parameterization increased the NRMSE only by about 7 %.
We believe that the parameterization for fractional SCA is also applicable during the accumulation or melt season, during other winters and in a different geographic region.However, our assumption requires verification once highly resolved spatial snow depth data become available, preferably in different snow climates, at times other than at the peak of winter, and from less topographical influenced regions.By performing snow depths measurements over several winter seasons, persistent snow depth distributions at the peak of winter were already found (Luce et al., 1999;Schirmer et al., 2011).These findings suggest that our parameterization for σ HS should be applicable during other winters, and motivates to investigate the evolution of spatial distributions of snow depth throughout the (melting) season.
Regarding grid cell size, horizontal resolutions of largescale meteorological and regional climate models can be much larger than our largest tested grid cell size of 3 km.However, at these larger scales, the presented parameterizations should also be applicable.To mimic the dominant snow-cover shaping processes in a domain size L, the domain size has to be substantially larger than the subgrid topographic correlation length ξ , i.e., L ξ .Note that every grid cell shows a different terrain correlation length ξ due to different subgrid topographies.A reliable subgrid parameterization for σ HS (Eq.2) was therefore derived by including a scaling parameter that corrects for finite L/ξ ratios (as in Helbig and Löwe, 2014).
We summarize that the subgrid parameterization for σ HS depends on both terrain length scales and on mean snow depth, which allowed developing a parameterization for the SCA over complex topography independent of a specific geographic region.A parameterization for the SCA over mountainous terrain has several practical applications.For instance, from SCA surface albedos can be derived to improve radiation balance estimates in large-scale meteorological models.Moreover, a SCA parameterization can be used to improve simulations of averaged snowmelt fluxes in large grid cells, which is a relevant issue for flood warnings.

Figure 1 .
Figure 1.Maps of (a) measured snow depths at Wannengrat and in the Dischma area in the eastern Swiss Alps and (b) hillshade at Val de Núria in the eastern part of the Spanish Pyrenees.The black squares illustrate examples of our randomly selected domain sizes of varying size.The underlying pixelmap (1 : 2 000 000) in (a) stems from swisstopo ©2008.The picture in (b) is taken from Grünewald et al. (2013).

Figure 2 .
Figure 2. Probability density functions (pdf) of measured snow depths are shown for the three areas.

Figure 3 .
Figure 3.One example probability density function (pdf) of measured snow depths HS for each domain size L (in color) in each area.

Figure 4 .
Figure 4. Mean root mean square errors (RMSE) between theoretical probability density functions (pdf) and measured pdfs as function of domain size L. Error bars indicate standard deviation of RM-SEs.

Figure 5 .
Figure 5.Standard deviation of snow depth σ HS as function of domain size L for all three areas.The squares represent mean σ HS .

Figure 6 .
Figure 6.Standard deviation of snow depth σ HS as a function of detrended valley-to-peak elevation difference σ (indicated by σ z ) of the underlying topographic features.Colors indicate corresponding domain size L.

Figure 7 .
Figure 7. Measured mean snow depth HS as function of mean measured flat field HS flat for all three areas.Colors indicate corresponding domain size L.

Figure 8 .
Figure 8. Measured standard deviation of snow depth σ m HS as function of parameterized standard deviation of snow depth σ HS for all three areas.(a) Parameterized via Eq.(2), (b) parameterized viaEgli and Jonas (2009) and (c) parameterized viaPomeroy et al. (2004).Colors indicate corresponding domain size L. NRMSEs are given for each parameterization.

Figure 9 .
Figure9.Snow-cover depletion curves derived assuming normally distributed snow depth and homogeneous melt via Eq.(4) as function of mean HS normalized with the peak of winter mean snow depth HS 0 (indicated by the subscript 0) (solid lines).Dashed lines represent parameterized fractional SCA f via Eq.(7).Coefficient of variations CV vary between 0.06 (upper left) and 1 (lowest one).

Table 1 .
Pearson correlation coefficients r for mean snow depth HS and standard deviation of snow depth σ HS with terrain parameters for pooled data of all three catchments as well as for each catchment separately.Gaussian covariance parameters σ (σ z ) and ξ are obtained as described in Sect.3.2.For mean slope µ, see Eq. (1).Values in bold indicate statistically significant correlations (p values < 0.05).