Analysis of spatial and temporal distribution

Seasonal snow cover is analyzed based on DEMs with a grid size of 1 m. The surface elevation at ALS(t 1) is subtracted from the surface elevation at ALS(t 2) to calculate surface elevation change, Δz _ALS:

(1)

Values of Δz _ALS are interpreted as snow depth, HS. HS values in the range −10 m ≤ HS ≤ 15 m are used for the analysis. Negative HS, which are an artifact of the measurement principle (Reference Hopkinson, Hardy and FrankensteinHopkinson and others, 2001), are not removed, to maintain the original distribution of the ALS-derived information. The temporal variability of the mean annual snow depth, , and its spatial variability in terms of the standard deviation, σ, and the corresponding CV are analyzed. A threshold HS of 0.15 m is used to distinguish between snow-free and snow-covered areas. According to the assumed relative vertical accuracy (Reference Bollmann, Sailer, Briese and FritzmannBollmann and others, 2011), this threshold accounts for the distribution of zero snow depths at ±0.15 m, i.e. HS < 0.15 m can result from zero snow depth and DEM error. Areas of HS < 0.15 m are also assumed to feature a scattered snow cover. The interannual persistence of HS and snow-covered area, SCA, with altitude is analyzed in 25 m elevation zones.

Pearson correlation coefficients, r, of HS are calculated between each combination of the five HS distributions. HS are sampled to coarser grid sizes (10, 30, 50 and 100 m) using the area-weighted mean of the 1 m HS values. Variability is known to decrease with increase of spatial aggregation (e.g. Reference Gutknecht and KirnbauerBlöschl, 1999). The spatial scale of 10 m is chosen to reduce the noise of the original 1 m data without losing too much information about snow-cover variability. A spatial resolution of 30 m was used by Reference Sturm and WagnerSturm and Wagner (2010) and Reference Winstral and MarksWinstral and Marks (2014), who analyzed manually sampled snow depth data. A typical scale of 50 m is used for the simulation of high-mountain snow cover and mountain hydrology (e.g. Reference StrasserStrasser, 2008; Schöber and others, in press). Spatial scales of >100 m show a reduction in variability of high-resolution snow depth data (e.g. Grünewald and others, 2013) and cause a distinct loss in information on spatial snow depth patterns (e.g. Reference Melvold and SkaugenMelvold and Skaugen, 2013).

Spatially distributed standardized snow depths, HS _s, are calculated according to Reference Sturm and WagnerSturm and Wagner (2010) for the interannual comparison of seasonal snow depths, using

(2)

where σ_HS is the annual areal standard deviation of HS in one winter. Squared errors of HS _s (SE_s) are calculated for all possible combinations of two available snow depth distributions (1 and 2) using

(3)

The resulting maps of SE_s (ten in total) are used to identify regions of high interannual snow depth variability and areas of interannually persistent HS _s. Areal means of the squared errors of standardized snow depths, MSE_s, are used as a metric of overall similarity of the snow depth distributions.

Residual analysis

Spatial snow-cover variability is expressed by residuals of HS (r_HS ) to their vicinity. We use a moving-window approach to calculate the residual of each HS to the mean HS of its squared vicinity of 5 × 5 raster cells:

(4)

where n is the number of cells of the window and i, j are the distances in the x – and y –directions from the center. While a homogeneous snow cover is characterized by almost zero r_HS , a heterogeneous snow cover shows distinctly negative and positive values of r_HS . The magnitude of r_HS is controlled by spatial variations in the underlying topography and variability of the snow on its surface. Values of r_HS are calculated at a raster resolution of 10 m within a 50 m square. Thus, the residual analysis covers the scale for which snow typically shows a high spatial autocorrelation (e.g. Reference Shook and GrayShook and Gray, 1996; Reference Deems, Fassnacht and ElderDeems and others, 2008; Reference Schirmer and LehningSchirmer and Lehning, 2011). We use a k –means clustering algorithm (Reference ForgyForgy, 1965; Reference MacQueen, Le Cam and NeymanMacQueen, 1967; Reference Hartigan and WongHartigan and Wong, 1979) integrated in SAGA (System for Automated Geoscientific Analyses) to automatically group observations with similar characteristics in magnitude and temporal variability of r_HS . A set of ~36 000 spatially distributed observations of r_HS is used. Each observation consists of a five-dimensional vector representing r_HS values over the five years. This clustering method does not create spatially connected clusters (e.g. as image segmentation algorithms). Rather, the high number of observations distributed over the entire basin are grouped into classes. We clustered the observations to a pre-set of ten classes. Frequency distributions are calculated for the annual r_HS of all observations grouped in each of the ten classes. The advantage of clustering the r_HS is the spatially distributed information obtained on the temporal persistence of snow depth residuals to mean snow depth in their vicinity.

Geostatistical analysis

Geostatistical analyses are used in cryospheric sciences for scaling issues (e.g. Reference Shook and GrayShook and Gray, 1996; Blöschl, 1999), to investigate topographic and wind control on spatial snow depth distribution and to analyze regional and temporal variability of fractal dimensions and typical scale breaks of mountain snow cover (e.g. Reference Deems, Fassnacht and ElderDeems and others, 2006, Reference Deems, Fassnacht and Elder2008; Reference TrujilloTrujillo, 2007; Reference Lehning and SchirmerLehning and others, 2011; Reference Mott, Schirmer and LehningMott and others, 2011a; Reference Schirmer and LehningSchirmer and Lehning, 2011).

Experimental variograms are used to analyze the spatial autocorrelation of data pairs as a function of their separation distance h (Reference OleaOlea, 1999; Reference Sun, Xu, Gong and LiangSun and others, 2006). The semivariance, γ (h_k ), of HS is given for the mean lag distance h of bin k as

(5)

where HS are snow depths at locations i and i + h, and n_k is the number of data pairs in bin k. The spatial characteristics of HS can be explored by the properties of the variogram, i.e. nugget, sill and range (Reference OleaOlea, 1999). The nugget indicates small-scale variability that has not been captured by the measurement. The sill is the overall variability, γ (h_k ), in the dataset, i.e. the higher the sill, the higher the spatial variability in the data. The distance between sampling points (range) where the increase in γ (h_k ) with increasing h changes considerably is called the scale break distance, L. In the case of flattening of the increase in γ (h_k ) on a log-log plot, L marks the separation of the sampling points at which there is no further autocorrelation between data pairs.

The fractal (Hausdorff–Besicovich) dimension D was found to be useful to analyze the irregularity of objects and surfaces (Reference BurroughBurrough, 1993; Reference Gao and XiaGao and Xia, 1996). Nearly planar surface features have D = 2, while a rugged surface has a limiting D = 3.

Fractal dimensions are estimated from the slope of a power law fitted to the log-spaced data

(6)

where α is the ordinal intercept and the power β is the slope of the linear fit in the log-transformed variogram. The fractal dimension D is calculated from β (Reference Gao and XiaGao and Xia, 1996):

(7)

The application and interpretation of geostatistics is strongly limited by the measurement scale and the sampling design. For example, if the spacing between the sampling points is large compared with the extent of the sampled area, the true variability in the data might not be captured. The high spatial resolution of the ALS data supports the applicability of the semivariance analysis. However, the results of L and D are highly influenced by the methods used. Sources of different estimates for scale breaks and fractal dimensions are discussed by Reference Xu, Moore and GallantXu and others (1993) and Reference Sun, Xu, Gong and LiangSun and others (2006).

For this study, semivariance is calculated for HS over the entire ice-free terrain with respect to the largest glacier extent in 2001 (17.75 km²), separately for HS on Hintereisferner (HEF; 6.99 km²) and on Kesselwandferner (KWF; 3.65 km²) and for HS on the remaining small glaciers (2.00 km²). With respect to the results of other studies analyzing geostatistical properties of Alpine HS (e.g. Reference Deems, Fassnacht and ElderDeems and others, 2008; Reference MottMott and others, 2011b; Reference Schirmer and LehningSchirmer and Lehning, 2011), scale breaks in the order of meters to tens of meters are expected. Considering this length characteristic we calculate omni-directional variograms with 2 m grid spacing. The lag distance h is subdivided into 20 log-spaced bins for each power of ten. The maximum h is set to 200 m for ice-free terrain. With respect to the more homogeneous snow cover on glaciers, semivariance is calculated with a maximum h = 1000 m for small glaciers and of h = 3000 m for the two large glaciers, HEF and KWF.

For ice-free terrain, a preliminary first break distance L _p is calculated by automatically fitting an exponential function to h and using a threshold of 95% of γ (h_k ) (Reference Webster and OliverWebster and Oliver, 2007). This automatically fitted L is only used to separate the data of the used bins into a short-range set and a long-range set. Power laws are fitted to the two datasets in the log-transformed variogram, solving least squares to calculate the short-range fractal dimension, D _s, for h ≤ L _p and the long-range fractal dimension, D _l, for h > L _p. The intersection of the two fitted lines presents the final L, which depends on the slope of the power laws and differs from L _p.

In contrast to Reference Shook and GrayShook and Gray (1996), who used a fixed D _l = 3 (completely random spatial distribution), we also calculated D _l for the long-range section.

Spatial and temporal snow distribution

Mean snow depth at the end of each accumulation season, , and its spatial and temporal variability are shown in Table 3. ranged from 1.31 m in 2011 to 2.58 m in 2009. This range corresponds to 68% of the temporal mean, , of 1.88 m. Whereas the temporal CV of and the long-term temporal CV of winter precipitation measured at cumulative gauges in the investigated catchment (Table 2) are similar, the observed range between minimum and maximum winter precipitation is higher than the range of in the five investigated years.

Table 3. Mean snow depth, HS, the standard deviation of spatially distributed snow depth, σ_HS , the corresponding spatial CV and the snow-covered area, SCA, for each accumulation season calculated from 1 m gridded snow depth data. The interannual mean and interannual variability in terms of range, σ and temporal CV are also shown

Snow-covered area, SCA, varied interannually by 16% of its mean value. The temporal mean snow-covered area fraction was 87% of the total area. The standard deviation of spatially distributed HS (σ_HS ) was lowest in the year of minimum , 2011. σ_HS was similar in 2002, 2003 and 2012. In these years, the spatial CV varies due to interannual variations in . Maximum σ_HS occurred in 2009, whereas the spatial CV was lowest in this year of highest . This shows the great influence of on the spatial CV of one year compared with the influence of σ_HS . The mean spatial CV of was 0.72. The temporal CV shows that was interannually more variable than σ_HS and the spatial CV.

Pearson correlation coefficients (r; Table 4) of spatially distributed HS show a high correlation of HS between all accumulation seasons. At the high spatial resolution of 1 m, r is highest for the two seasons with a high fraction of snow-free area in 2011 and 2012 and second highest for the two seasons of the lowest in 2011 and the highest in 2009. The correlations are lower for combinations that include the HS distributions of 2002 and 2003. Thus, there is no clear dependence of correlation on . In general, correlations of HS distributions decrease with increasing time lag between the two seasons. At a spatial resolution of 100 m, correlations between the HS distribution of 2003 and later years are similar, but correlations including HS distributions of 2002 still show lower values. In general, an increase in correlation with decrease in raster resolution is obvious, due to a reduction in small-scale variability by spatial aggregation, as shown by, for example, Reference Melvold and SkaugenMelvold and Skaugen (2013).

Table 4. Correlation coefficients, r, of spatially distributed HS calculated in raster resolutions from 1 to 100 m for each possible combination of the individual accumulation seasons

A high Pearson correlation (r = 0: 94) was found between and SCA with lower , indicating less SCA. This is clear from Figure 3, where mean HS and SCA fraction of the investigated years are shown for 25 m elevation zones. In this figure, large triangles represent elevation zones with a large area and, thus, a high contribution to total snow-cover volume compared to elevation zones marked with small triangles. In these elevations of largest contribution to total catchment area, HS varied by a factor of two, whereas SCA fraction varied by only 11% (between 87% in 2011 and 98% in 2009).

Fig. 3. Mean snow depth, HS (solid lines), and SCA fraction (triangles) of 25 m elevation zones. The size of the triangles represents the area of the corresponding elevation zone scaled by the total investigation area, i.e. small triangles represent only small areas with less contribution to total snow-cover volume compared with the contribution from areas marked with large triangles.

In general, the annual elevation distribution of HS and SCA fraction seem to run parallel (Fig. 3). HS increased to an elevation of 3100 m a.s.l. and decreased above 3200 m a.s.l. Interannual differences of HS and SCA fraction can be found at the lowest elevations and at mountain ridges >3400 m a.s.l. Differences in mean HS of up to 0.5 m occurred for interannually similar SCA fractions <2800 m a.s.l. (i.e. 2002 and 2009).

The remarkably high fraction of snow-free areas in 2011 and 2012 is clear at HS = 0 ± 0: 15 m in Figure 4. The year with lowest (2011) shows a peak of most frequent HS at 2 m. In contrast, the year with highest (2009) has a peak of most frequent HS at 3.2 m. Frequency distributions of standardized snow depths, HS _s, were found to be inter-annually more consistent than frequency distributions of HS. Hence, HS _s was chosen to distinguish between areas of persistent snow cover and areas of high interannual snow depth variability.

Fig. 4. Frequency distribution of snow depth, HS, and standardized snow depth, HS _s (Eqns (1) and (2)), for the five accumulation seasons. The uncertainty of ALS data (±0.15 m) is shaded in blue at HS = 0.

The temporal mean of HS _s (Fig. 5a) represents the standard snow pattern in the investigation area using the available ALS data. South-facing slopes at lower elevations showed almost no snow cover at the time of ALS(t 2). At higher elevations, low snow coverage can be found at sheer rock faces and along mountain crests. In general, HS _s on the glacier surfaces was higher than HS _s in ice-free areas. Increased HS _s can be found at the base of steep slopes and along the glacier margins. These are the areas of high temporal standard deviation of HS _s (σ HS _s) with respect to the available ALS data.

Fig. 5. Spatial distribution of (a) temporal mean standardized snow depth, HS _s, and (b) its interannual standard deviation, σHS _s, for all five accumulation seasons. A, B, C and D mark regions of high interannual variation of HS _s.

Areas A and B in Figure 5b are examples of regions showing high σ HS _s as a result of avalanche release and accumulation. Corresponding profiles of HS are shown in Figure 6. In both areas, HS for 2002 is low at steep slopes and high on glacier surfaces within ~100 m of the slope base. In profile A, high HS can be found directly at the slope base in the other years. Peaks of HS at a distances between 0 and 100 m in profile A and at 200 m in profile B are caused by crevasses at the ALS(t 1) reference surface. A substantial snow volume was kept on the slope in area B (Fig. 6b) and not released to the glacier surface prior to ALS(t 2) in 2009. Box C highlights a typical dead-ice area next to the tongue of Hintereisferner (investigated by, e.g., Reference Bollmann, Sailer, Briese and FritzmannBollmann and others, 2011; Reference Sailer, Bollmann, Hoinkes, Rieg and SprossSailer and others, 2012). In dead-ice zones, the high temporal variability of HS _s results from spatially differing surface lowering within the 10 year period. The retreat of the glacier tongue of Kesselwandferner in the 10 year period caused a high interannual variability of HS _s in box D (see also next subsection). In contrast to zones A–D, large parts of the investigation area have low σ HS _s, which is indicated by the dominant blue color graduation in Figure 5b. An area fraction of 75% corresponds to a threshold σ HS _s of 0.4 (Fig. 7). This large section is characterized by an inter-annually persistent HS _s. Conversely, about a quarter of the total area shows a high interannual variability of HS _s. These areas were found to be not necessarily located in steep terrain, and are more frequent at the base of steep slopes and along glacier margins (Figs 1 and 5).

Fig. 6. Snow depth HS along profiles A and B (Fig. 5b) for the five accumulation seasons. The dotted lines show the corresponding mean elevation of the ALS(t 1) surface (right axis).

Fig. 7. Cumulative probability density function (PDF) of the temporal standard deviation of standardized snow depths, σHS _s.

The multi-annual cross-combination of squared errors of HS _s (SE_s; Eqn (3)) helps in analysis of the temporal frequency of the processes causing high interannual variability of snow depths, and to explain the magnitude of spatial means of SE_s (MSE_s; Fig. 8). Most parts of the areas featuring SE_s <0.5 (masked white in Fig. 8) correspond to areas of σ HS _s < 0. 4 (Fig. 5b).

Fig. 8. Maps of squared errors of standardized snow depths, SE_s (Eqn (3)), between all possible combinations of the five accumulation seasons. In the lower-left corner, mean squared errors, MSE, of the cross-combinations are shown for different spatial resolutions.

High SE_s are clear within areas A and B (Fig. 8, upper row), where snow variability is a result of the release and accumulation of avalanches (Fig. 6). These interannual differences of HS _s are especially pronounced in combinations that include 2002. High SE_s between the years 2002, 2003 and the more recent surveys 2009, 2011 and 2012 are clear for areas C and D along the glacier margins. In these areas, high values of SE_s are caused by continuous changes in topography, due to melt of dead ice, collapse of moraines and glacier retreat within at least 6 years between the surveys (e.g. Reference FischerFischer, 2011; Reference Sailer, Bollmann, Hoinkes, Rieg and SprossSailer and others, 2012). Thus, MSE_s (Fig. 8) and correlation coefficients (r; Table 4) are highest for cross-combinations including 2002 and lower for shorter time lags between the years. Even in the 2 year period 2009–11 and the 3 year period 2009–12, high SE_s can be found along the glacier margins (Fig. 5b, box C). SE_s values along the glacier margins are small for the 1 year periods 2002/03 and 2011/12, which reduces MSE_s for these cross-combinations.

Residual analysis

While high interannual variability was detected for areas along the glacier margins, it has not yet been shown what causes the interannual variability of SE_s in these zones. The temporal clustering of the residuals of each HS value from its surrounding mean (r_HS ; Eqn (4)) is applied to show more detail of the differences in the areas highlighted in the previous subsection. The clustering of r_HS resulted in seven classes. Three of the ten initial classes are not discussed further because they originate from extreme values of individual seasons and cover only a small relative fraction (0.1%) of the total area. A map of the clusters with a zoom to the area of the tongue of KWF is shown in Figure 9.

Fig. 9. Clusters derived by k –means clustering of snow depth residuals, r_HS . The letters pair classes with similar temporal characteristics in magnitude of r_HS . Index 1 indicates negative r_HS and index 2 indicates positive r_HS . The detailed view highlights areas of changes in snow accumulation as a result of glacier retreat (black glacier outlines) in the 10 year period 2001–11.

Figure 10 shows the frequency distribution of r_HS within the classes E–H and for the total area. Sixty-three percent (22.5 km²) of the investigated area is characterized by a homogeneous snow cover, predominantly located on glaciers and in areas with no snow in any of the five years (class H). In these areas the mean r_HS values are close to zero (maximum σ = 0.18 m) in all five accumulation seasons. The transition zone from homogeneous (r_HS ≠ 0) to heterogeneous (r_HS 0) snow distribution (class F1 and F2) covers 25% of the total area. Frequency distributions of r_HS in F1and F2 were similar for the individual accumulation seasons (Fig. 10). More interannual variability of r_HS is clear in classes E1 and E2 of distinct positive and negative r_HS (Fig. 10). Maximum and minimum r_HS were both found in 2009. This shows that the spatial variability of HS increases in years of high , i.e. HS increases in sinks and channels, whereas crests and bumps remain snow-free. In 2002, the frequency distribution of r_HS shifted towards lower values in areas of positive r_HS (Fig. 10). This resulted from a snow cover accumulated at the mountain crests prior to ALS(t 1). Classes G1 and G2 are characterized by a change in r_HS within the five accumulation seasons. These two classes are located at the glacier margins (Fig. 9) and together cover 6% of the total area. They show homogeneous snow cover in 2002 and 2003, but distinct positive or negative r_HS in 2009, 2011 and 2012 (Fig. 10). This change is a result of glacier retreat, which exposes the more heterogeneous surface of the former glacier bed compared with the homogeneous glacier surface. However, differences in these two classes are not visible in the frequency distribution of r_HS for the total area, due to the large fraction of homogeneous snow cover (class H) within the investigation area.

Fig. 10. Histograms of annual r_HS for the clusters derived by k –means clustering of snow depth residuals, r_HS . Letters E–H are the classes shown in Figure 9. Note that the scales on the area axes vary between the classes.

Geostatistical analysis

Variograms were calculated for HS in ice-free terrain, for HS on Hintereisferner (HEF) and Kesselwandferner (KWF) and for HS on the remaining small glaciers (Fig. 11). The scale break, L, the corresponding sill, γ (h_L ), the short-range fractal dimension, D _s, and the long-range fractal dimension, D _l, were calculated for ice-free terrain (Table 5).

Fig. 11. Variograms for HS of the five accumulation seasons. (a) In ice-free terrain, (b) on small glaciers, (c) on Hintereisferner and (d) on Kesselwandferner. Vertical lines in (a) show the location of the calculated scale breaks, L.

Table 5. Annual scale break, L, the corresponding sill, γ (L), short-range fractal dimension, D _s, and long-range fractal dimension, D _l, calculated for HS in ice-free terrain

An interannually consistent L was found in ice-free terrain at ~20 m, whereas γ (h_L ) ranged considerably between the two years of lowest and highest (Table 5). The short-range variability of HS is highest in 2002 (Fig. 11a). This results in a higher D _s than the interannually persistent D _s of the other four years (Table 5). In all five years, the increase of γ (h) with h flattens out to a consistent D _l of 2.9 for HS in ice-free terrain. γ (h) for h > L depends on the spatial mean snow depth, , of the individual years with increasing γ (h) for increasing .

In contrast, no distinct flattening of γ (h) with increasing h could be found for h > 10 m on the glacier surface. The overall spatial variability of HS in ice-free terrain was found to be higher than γ (h) on glaciers. A break in scaling behavior of HS on glaciers is indicated at L = 10 m (Fig. 11b–d). However, γ (h) is of the magnitude of a few centimeters to one decimeter at h = 10 m and, thus, in the range of the assumed accuracy of Δz _ALS. At h = 10 m, γ (h) is slightly higher on HEF than on small glaciers and on KWF.

At scales of h < 10 m, γ (h) is highest in 2003 and 2012 on all glaciers. Although γ (h) of HS is lowest in 2011, differences can be found between γ (h) in different accumulation seasons on the glaciers. In particular, HS in 2002 shows low variability at shorter scales, but γ (h) increases towards scales of ~100 m on small glaciers and on HEF. However, variability of HS in 2002 stays at comparatively low values on KWF. The variability of HS on KWF in 2003 does not increase over a large range of h.

Spatial and temporal snow distribution

The interannual analysis of the snow cover in this partly glacierized catchment shows a higher range of (68% of its temporal mean) compared with the interannual range of SCA (16% of its temporal mean; Table 3). Nevertheless, annual SCA and of the total area have a high correlation (r = 0: 94). With respect to the total snow-cover volume, and hence , elevations of the largest area present a roughly parallel shift in elevation distribution of HS and SCA (Fig. 3). But differences in SCA fraction are small at these elevations. This is a result of the large fraction of glacierized area that is totally snow-covered in all years. At lower elevations, snowmelt caused different HS and variable SCA, due to the seasonal course of weather conditions. However, variable HS at these elevations has less influence on , due to their small area fraction. The small variation of SCA compared with the variation of might be challenging for the determination of relationships between snow-covered area fraction and snow depths, such as those used to calculate subgrid variability of snow-cover products from satellite data (e.g. Reference Swenson and LawrenceSwenson and Lawrence, 2012) or to model spatial variability of snow depths (e.g. Reference Luce, Tarboton and CooleyLuce and others, 1999; Reference ListonListon, 2004; Reference ClarkClark and others, 2011). At the time of peak snow accumulation, very detailed surveys of SCA over the largest-area elevation zones would be necessary to obtain mean snow depth of the partly glacierized catchment from small changes in SCA.

As a result of ice flow and densification of firn, HS can be underestimated in typical firn areas of the glaciers. This is more relevant for elevations >3250 m a.s.l. Reference Helfricht, Kuhn, Keuschnig and HeiligHelfricht and others (2014) compared ALS data and ground-penetrating radar measurements on the investigated glaciers and found a mean underestimation of HS of 0.4 m (σ = 0: 34 m) by ALS data in the firn areas of the glaciers. The decrease of mean HS in glacierized areas >3250 m a.s.l. (Fig. 3) can partly result from firn densification and ice dynamics. Areas of emergent ice flow were found to be small, having no considerable influence on the interpretation of HS. However, we assumed these dynamical processes to be constant in time, so no corrections were applied. Further, in many studies rough and steep surfaces were found to store less snow on average (Reference KirnbauerBlöschl and Kirnbauer, 1992; Reference LehningGrünewald and Lehning, 2011; Reference Helfricht, Seiser, Fischer and KuhnHelfricht and others, 2012; Grünewald and others, 2013; Schöber and others, in press).

An increase in annual σ_HS with increasing was found. This was also noted by Reference Sturm and WagnerSturm and Wagner (2010), who analyzed derived from snow depth samples taken on a regular 30 m grid within a 1 km × 1 km Arctic test area. Reference Melvold and SkaugenMelvold and Skaugen (2013) also found a decreasing spatial CV with increasing HS with values between 0.86 (= 1.42 m) and 0.62 ( = 2.89 m). Reference Winstral and MarksWinstral and Marks (2014) discussed the variation of spatial CV when changes at a faster rate than σ. In their study, the mean spatial CV over 11 years was 0.71, with a slightly higher value for snow surveys in spring (0.81). This finding is in agreement with the variability in mean spatial CV found in our study. Reference Winstral and MarksWinstral and Marks (2014) also found a high spatial CV, of 0.89, for 1 year ALS data. Reference Winstral and MarksWinstral and Marks (2014) related higher spatial CV of HS to higher wind speeds and higher melt rates. To some degree this might be true for the spatial CV of HS reported here, in 2011 and 2012, when potential melt conditions appeared prior to ALS(t 2).

Marchand and Killingtveit (2005) showed that the spatial CV is in the range 0.72–0.96 for open sites in Norway. Reference Jepsen, Molotch, Williams, Rittger and SickmanJepsen and others (2012) analyzed a total of 5300 snow depth measurements taken at the end of 12 accumulation seasons. They found the spatial CV of the snow cover ranged from 0.73 to 1.09 in a continental Alpine watershed. In their study, the temporal variability of the spatial CV was 0.15, which is higher than the temporal variability of CV of 0.1 found here. Reference Schirmer, Mott and LehningGrünewald and others (2010) derived a spatial CV of HS between 0.6 and 2.0 from terrestrial lidar data, with the lowest values at the beginning of the ablation season.

In general, the temporal mean spatial CV of ALS-derived HS in this study (0.72) is at the low end of values derived for open and high-elevation alpine areas in previous studies. Schöber and others (in press) present a map of spatially distributed CVs of HS in 2011, which includes the investigation area of the present study and shows low CV values on glaciers. Thus, with respect to the homogeneous snow cover on the glacier surface, annual and multi-annual mean spatial CVs of HS in partly glacierized catchments tend to be lower than the CVs in ice-free terrain.

Here we find the maximum correlation across two years is 0.84 (2011/2012) for the original spatial resolution of 1 m. This is lower than interannual correlations found for HS by Reference Deems, Fassnacht and ElderDeems and others (2008) (r = 0: 88–0.92), Reference Melvold and SkaugenMelvold and Skaugen (2013) (r = 0: 95) or Reference Schirmer and LehningSchirmer and others (2011) (up to r = 0: 97). This results from the large area and complex terrain investigated in our study. The wide range of snow depths (snow-free to HS > 5 m in avalanche accumulation areas), the high number of different topographic features (glaciers, rock walls, slopes, mountain crests, etc.) and the corresponding dynamical processes (glacier retreat, rockfall, collapse of moraines, etc.) reduce the interannual correlation compared with smaller, temporally more consistent test sites.

Since HS is scaled for statistical models to simulate snow-cover distribution (e.g. Reference Sturm and WagnerSturm and Wagner, 2010; Reference Lehning and SchirmerLehning and others, 2011; Grünewald and others, 2013), different measures of snow depth have been used. We found standardized snow depth, HS _s, best suited to describe the interannual frequency distribution of HS of this area. Thus, we used HS _s to distinguish between areas with interannually persistent snow cover and areas with high interannual variability. In the same way as Reference Sturm and WagnerSturm and Wagner (2010), we showed a climatological snow distribution pattern in terms of spatially distributed, temporal mean HS _s of the investigated snow surveys. Although Reference Sturm and WagnerSturm and Wagner (2010) found that interannually persistent areas were visually striking, they did not present a measure to differentiate between interannually variable HS _s and areas of persistent snow cover. We used the temporal standard deviation of HS _s (HS _s) (Fig. 5b) to demonstrate that 75% of the investigated area showed interannually consistent HS _s for a threshold σHS _s of 0.4.

In general, σHS _s is low on glaciated areas and on gentle slopes (Fig. 5b). Especially in these areas, lidar measurements of HS have an appropriate accuracy and most of the snow-cover volume is accumulated there. Thus, ALS is suitable for recording characteristic snow depth distributions for the calibration and validation of hydro-meteorological models (e.g. Schöber and others, in press). High HS values were found consistently on glaciers characterized by flat, sheltered surfaces at high elevations. The preferential accumulation in these areas is also the reason why glaciers exist where they do. For instance, Reference Dadic, Mott, Lehning and BurlandoDadic and others (2010a) related the location of glaciers to differences in the wind field. Redistribution of snow to glacier surfaces has already been included in hydro-meteorological models adapted to the investigated catchment (Reference KuhnKuhn, 2003).

The reduction of MSE_s with increasing aggregation of raster resolution (Fig. 8) is distinctly higher for MSE_s values that include 2002 or 2003 compared with MSE_s calculated between surveys since 2009. Small-scale features of high SE_s are smoothed, especially in narrow areas along sheer rock walls and glacier margins. Two studies have calculated mean-squared errors of HS _s (MSE_s) for multi-annual data. Sturm and Wagner (2010) found HS _s for snow surveys in March and April compared with one survey in March in the range 0.31–1.34. Lower MSE_s values were found by Reference Winstral and MarksWinstral and Marks (2014). Based on the snow depth data over 11 years at a spatial resolution of 30 m, their calculated MSE_s was in the range 0.06–0.6, with a mean MSE_s of 0.31. This range is comparable with the mean MSE_s found in our study (0.22–0.5) with a mean of 0.34 for a spatial resolution of 30 m (Fig. 8).

Areas of high interannual snow depth variability (i.e. high σHS _s) are not restricted to mountain crests and steep slopes, where ALS data are affected by large errors. In general, 82% of all σHS _s ≥ 0.4 were located at slopes of <40°. High σHS _s values were found at gently inclined slope toes and on homogeneous glacier surfaces (Fig. 7). These are typical accumulation areas of gravitational snow transport by avalanches and snow sloughs.

Geomorphological processes continuously change the mountains’ surface topography (Reference Bollmann, Sailer, Briese and FritzmannBollmann and others, 2011; Reference Sailer, Bollmann, Hoinkes, Rieg and SprossSailer and others, 2012). Thus, the variation in the reference surface at ALS(t 1) affects snow depth variability. In addition, glacier retreat alters the surface relief from a more homogeneous surface to heterogeneous landforms featuring ridges and sinks, which lead to a heterogeneous snow depth distribution (Figs 9 and 10). This coincides with the findings of Reference Lehning and SchirmerLehning and others (2011) and Reference Schirmer and LehningSchirmer and others (2011), who demonstrated that topographic roughness controls snow-cover variability. The shift of the frequency distribution of r_HS towards more positive and negative values in the year of highest mean snow depth (2009) shows that differences of HS between preferential accumulation areas and areas with less or no snow increase with an increasing .

Scaling behavior of snow depths

We assessed the multi-scaling behavior of HS separated by a scale break distance, L, into a short-range fractal dimension, D _s, and a long-range fractal dimension, D _l. In ice-free terrain, L was found to be interannually consistent at ~20 m, comparable with the L = 17.3 and 20.4 m found by Schirmer and Lehning (2011) for a cross-loaded slope (i.e. variable wind direction) in complex mountain terrain. However, L found by Reference Deems, Fassnacht and ElderDeems and others (2008) for the HS distribution of two seasons in a forested, moderate-elevation terrain with low rolling topography (L = 15.5 and 15.9 m) and at an Alpine site above the tree line (L = 31.5 and 26.4 m) are slightly lower and slightly higher, respectively. Reference Mott, Schirmer and LehningMott and others (2011a) presented very similar values to those found in this study for peak HS calculated from terrestrial laser scanning measurements in a small Alpine area (L = 20: 4, D _s = 2: 31), and Reference Schirmer and LehningSchirmer and Lehning (2011) showed a reduction in D _s towards values between 2.3 and 2.5 at the end of the accumulation season. D _l was found to be consistently close to 2.9 in all studies.

Here we found that L determined for HS in ice-free terrain marked the beginning of a curve-like reduction in the increase of γ (h), rather than a distinct bend in γ (h), unlike, for instance, those found by Reference Mott, Schirmer and LehningMott and others (2011a) and Reference Schirmer and LehningSchirmer and Lehning (2011). This is a result of the large area (~14 km²) and the high number of different surface types (rocks, boulders, talus, etc.) in a complex topography (sheer walls, gentle slopes, flat valley floors) with different exposure to meteorological conditions (radiation, wind). Since values of L depend strongly on the bin spacing, the grid resolution, the area considered for the variogram and the maximum range of h used for calculation, they give an order of magnitude (here a few tenths of a meter) rather than exact values.

We found γ (h) of HS at h = 2 m is highest in 2002 (0.23 m) and stays high for h lower than the scale break, L (Fig. 11a). Whereas γ (h) at h < 4 m might be influenced by the accuracy of the ALS data (Section 2), similarities between snow-cover distribution in 2002 and 2009 can also be found in SCA and HS at lower elevations (Fig. 3). In addition, the avalanche activity in 2002 (Fig. 7) increased the small-scale variability of HS. No distinct, interannually persistent scale break could be detected on glaciers. Nevertheless, the increase in γ (h) appeared to be higher at scales below 10^-1m compared with the increase in γ (h) at larger h. Even though this variability is clear in some years, γ (h) below 10^–1 m is close to the accuracy of Δz _ALS, of ±0.15 m (Section 2). Thus, an interpretation of this scale break would be speculative. However, years without considerable snowfall prior to ALS(t 1) (2003, 2012) show higher γ (h) at short scales compared with years where snowfall occurred prior to the survey in autumn (see Section 3). This is a result of crevasses and surface structures, such as meltwater channels and moraines on the glacier surface which, if they are not covered by snow at ALS(t 1), increase small-scale variability.

Especially in 2011, γ (h) of HS was reduced over the whole range of h due to snow accumulation prior to ALS (t 1), the comparatively low precipitation over this season and the fact that melt had begun at the snow surface on the glaciers prior to ALS(t 2). The lower values of γ (h) on glaciers in comparison with γ (h) in ice-free terrain highlight the more homogeneous snow depth distribution on the ice surfaces.