Processes governing snow ablation in alpine terrain – detailed measurements from the Canadian Rockies

The spatial distribution of snow water equivalent (SWE) and melt are important for estimating areal melt rates and snow-cover depletion (SCD) dynamics but are rarely measured in detail during the late ablation period. This study contributes results from high-resolution observations made using large numbers of sequential aerial photographs taken from an unmanned aerial vehicle on an alpine ridge in the Fortress Mountain Snow Laboratory in the Canadian Rocky Mountains from May to July in 2015. Using structurefrom-motion and thresholding techniques, spatial maps of snow depth, snow cover and differences in snow depth (dHS) during ablation were generated in very high resolution as proxies for spatial SWE, spatial ablation rates and SCD. The results indicate that the initial distribution of snow depth was highly variable due to overwinter snow redistribution; thus, the subsequent distribution of dHS was also variable due to albedo, slope/aspect and other unaccountable differences. However, the initial distribution of snow depth was 5 times more variable than that of the subsequent dHS values, which varied by a factor of 2 between the north and south aspects. dHS patterns were somewhat spatially persistent over time but had an insubstantial impact on SCD curves, which were overwhelmingly governed by the initial distribution of snow depth. The reason for this is that only a weak spatial correlation developed between the initial snow depth and dHS. Previous research has shown that spatial correlations between SWE and ablation rates can strongly influence SCD curves. Reasons for the lack of a correlation in this study area were analysed and a generalisation to other regions was discussed. The following questions were posed: what is needed for a large spatial correlation between initial snow depth and dHS? When should snow depth and dHS be taken into account to correctly model SCD? The findings of this study suggest that hydrological and atmospheric models need to incorporate realistic distributions of SWE, melt energy and cold content; therefore, they must account for realistic correlations (i.e. not too large or too small) between SWE and melt in order to accurately model SCD.

suggests that in alpine terrain the question of relative contribution of spatially variable melt 28 rates or snow redistribution on SCD can be reduced to the question of whether such a 29 correlation between melt and SWE exists and how large it is. 30 Besides theoretical considerations, there are a number of existing field and modelling 31 studies on the relative importance of spatially variable melt or snow redistribution on SCD. 32 There are studies that have found the temporal progression of snow-cover depletion (SCD) to 1 be governed primarily by the variability caused by snow redistribution rather than the 2 variability caused by melt rate differences (Shook and Gray, 1994;Luce et al., 1998;3 Anderton, 2004;Egli et al., 2012). These studies show with spatial observations that snow 4 cover depletion (SCD) can be modelled with statistics derived during peak accumulation. measurements of spatial melt rate and SWE via snow depth (HS) and the change in HS (dHS) 10 by terrestrial LiDAR and applying a few measured bulk densities to estimate SWE and 11 ablation rates. They found that SWE and melt rate were spatially uncorrelated over most of 12 their ablation season, except for a correlation coefficient of r = -0.4 for one sub period. They 13 noted that the variability of SWE was much larger than the variability of melt rates. In the 14 same study area over an additional winter season Egli et al. (2012) calculated SCD curves that 15 assumed correlations between HS and the change in HS (dHS), however these curves deviated 16 substantially from observations, suggesting that such correlations did not exist. Neither study 17 examined why such correlations were absent. 18 On the other side, spatially varying melt ratescaused by differences in insolation due 19 to aspect (Marks and Dozier, 1992), net solar irradiance due to albedo differences (Skiles et 20 al., 2015), internal energy storage due to deep, cold snow (DeBeer and Pomeroy, 2010) and 21 advected energy due to bare ground (Mott et al., 2011; can alter this pre-melt SWE 22 distribution and, when correlated to SWE, result in substantial changes to SCD curves 23 (Pomeroy et  to consider also the spatial correlation between melt and SWE to correctly model SCD. One 25 attempt in this direction is shown in DeBeer and Pomeroy (2010)  This study wants to assess the the influence of peak SWE variability, melt rate 3 variability and especially their spatial correlations on SCD in alpine terrain using high-4 resolution distributed measurements rather than sparse manual sampling or relying on model 5 results. The use of high resolution measurements is potentially important because peak alpine 6 snow depth, and thus also peak SWE, is known to vary most substantially below the scale of 7 tens of metres in alpine environments. small scale aspect variations can be found in the slopes on both sides of the ridge. Extreme 26 south and north aspects are underrepresented in the terrain snowcovered at the beginning of 27 the study period (Fig. 2b). The snowcovered area is reasonably steep with two peaks in the 28 slope distribution at 10 o and 25 o (Fig. 2c). On both sides of the ridge the slope steepens to up 29 to 40 o . Vegetation played a role in snow deposition patterns, mainly in the lee of shrubs and 30 clusters of small trees in krummholz with heights up to 2 m. Areas within these vegetation 31 clusters were excluded from the study as vegetation degraded the digital surface models 1 (DSMs) derived from UAV SfM photogrammetry (see section 2.4). The included area only 2 covered bare or sparsely vegetated ground so that vegetation effects can be excluded. 3 Two weather stations are located at the ridge, one on top of the ridge (Fortress Ridge -4 FRG) and one in a south facing slope (Fortress Ridge South -FRS, Fig. 2a). The local 5 Fortress Mountain Snow Laboratory within the regional Canadian Rockies Hydrological 6 Observatory provided five more weather stations within less than 2 km distance of the ridge, 7 which were used to interpret weather situations and for quality control. 8

Unmanned aerial vehicle (UAV) Data acquisition 9
A "Sensefly Ebee RTK" fixed wing UAV was used with a modified consumer-grade 10 Canon Elph compact RGB camera. As a base station a Leica GS15 differential GPS system 11 communicated with the UAV to tag captured images with corrected geolocations. 12 Additionally, ground control points were measured with this differential GPS system, which 13 we had to restrict analysis to two melt periods. 28

Accuracy evaluation and manual measurements 1
The accuracy assessment of this rather new method to determine snow depth was 2 given a high priority and is described in detail by Harder et al. (2016) for this environment 3 and others. In short, 100 differential GPS measurements on bare ground were taken. 4 Approximately 60% of the area was bare at the beginning of the study period which allowed 5 distribution of GPS ground measurements over a large part of the study area (Fig. 3a) and 6 thus widespread detection of any general misalignment of DSMs or local tilts. These points 7 could be used for all available flights. Differential GPS measurements were also taken at the 8 snow surface on the day of the specific flights, but technical problems often allowed only 9 limited additional time for these surveys. For most of the days up to 20 differential GPS 10 measurements on snow could be taken. At these GPS measurement points, snow depth was 11 also manually measured and snow density was measured at approximately each third of these 12 points. Density measurements were not sufficient to confidently estimate SWE from snow 13 depth into SWE and ablation rates from differences in snow depth. As such HS and dHS are 14 detail. In short, from 100 measurements on bare ground, the root mean square errors of bare 19 ground surface elevation ranged between 4 and 15 cm with a mean of less than 9 cm. Over 20 snow with fewer measurements an increase in these error measures could not be detected. A 21 signal-to-noise ratio (SNR) was used to ensure that the signal of the UAV was sufficiently 22 larger than the error, defined as mean of the signal divided by the standard deviation of the 23 error. The potential impact of this error on the results presented is discussed in section 3.2. 24

Data analysis 22
To diagnose reasons for spatial differences in snow depth change (dHS), Pearson 23 correlation coefficients were calculated with several potential explaining variables as Slope, 24 Deviation from North, Solar Irradiance, Brightness, current snow depth (HS) and snow depth 25 at the beginning of the study period (HS0). Scatter plots also were visually inspected to detect 26 reasons for strong or weak correlations or non-linear dependencies. The scatter plots were too 27 dense to interpret visually because of the high resolution and so instead of plotting point pairs, 28 the density of point pairs in a limited area of the plot was visualized (e.g. in Fig. 5a). 29 Spatial dependencies of the spatial structure of dHS and its correlation with explaining 30 variables were analysed with variograms and correlograms. Variograms were calculated with 31 Correlations between two variables x and y in a certain lag distance h were calculated with the 3 cross-variogram as an estimator of the covariance (Webster and Oliver, 2007). 4 (2) 5 This covariance was scaled with estimators of the variance x ˆ(Eq. 1) using 6  Due to high wind speeds, large parts of the ridge were snow-free during most of the 24 exceptionally warm and dry winter season. After a large late November 2014 snow storm, the 25 FRG station rarely documented snow on the ground and shallow snowpacks that did form 26 were regularly eroded by wind within a few days. The snow covered area (SCA) reached the 1 seasonal maximum in late November after this substantial snowfall (80 mm) with light winds 2 and dropped dramatically due to subsequent wind redistribution from blowing snow storms. 3 When aerial measurements began on 19 May 2015, SCA was slightly below its typical winter 4 value as spring ablation was under way. Without excluding any areas (see section 2.4) SCA 5 was approximately 0.45 in both subareas (Fig. 3a). 6 Dust-on-snow was an obvious feature in late winter and the beginning of the melt 7 season (Fig. 3b). It had not been observed to such extent in over a decade of observations in 8 the region. This dust was locally eroded from the fine frost-shattered and saltation-pulverized 9 shale particles at the ridge-top and was transported by wind to adjacent lee slopes and into 10 gullies, similarly to wind-transported snow. Hence dust was deposited preferentially to snow 11 drifts. Subsequent snow accumulation and melt processes led to a dust-on-snow pattern of 12 high small-scale variability. The lower albedo from dust deposition may have influenced 13 snowmelt energetics, but its high variability is different from the large scale, areally uniform An example of reductions in snow depth (dHS) due to ablation over a period of 13 26 days is shown in Fig. 3d. At the first glance differences between aspects are obvious, as well 27 as smaller scale impact of albedo variations (cf. Fig. 3b). The driving forces to differences in 28 ablation inferred from the observed differences in depth change will be examined in 29 section 3.3. 30 The study covered the late melt period, when the highest ablation rates occurred. Peak 31 SWE of 500 mm was measured with a weighing snow lysimeter (Sommer "Snow Scale") in a 32 nearby forest clearing on 20 April 2015. By the start of the study period on 19 May, SWE had 1 gradually decreased to 300 mm, often interrupted by snowfall. During the study period after 2 19 May no significant (>3 cm) snowfall was observed. The much higher ablation rates 3 compared to the previous weeks caused the snow to disappear at this station on 30 May. A 4 very similar development could be observed at two other stations using snow depth sensors 5 within the Fortress Mountain Snow Observatory, including the FRS station (c.f. Fig. 2a). On 6 30 May a SCA of 0.2 was measured from the UAV over the whole flight domain. Considering 7 a typical pre-melt SCA of approximately 0.45, the presence of a significant SCA illustrates 8 the value of spatially distributed measurements of snow ablation and cover, when all seven 9 meteorological stations in the ~3 km 2 region were snow-free. 10 A meteorological overview during the study period is given in Fig. 4 Table 1 shows the Pearson's correlation coefficient for above mentioned melt periods 2 and different subareas. This univariate analysis shows clearly two driving factors for the 3 earlier melt period, P1, albedo and solar radiation differences, expressed respectively with 4 Brightness and either Deviation from North or with Solar Irradiance. The sign of the 5 correlations is mainly as expected: More southerly and darker pixels showed larger dHS 6 values. Exceptions (e.g. during P2 in the southern subarea) may be explained with observable 7 differences between a few remaining snow patches with different albedo values, slope, snow The correlations of dHS with Brightness, Deviation from North and Solar Irradiance 23 were often strong. dHS increased from 5 to 7 cm/d (nearly 60% increase) as aspect shifted 24 about 115 deg from north to south or snow from clean to dusty (c.f. Fig. 5b). This shows the 25 importance of spatial variation in net solar irradiance to melt energeticsas exemplified by 26 the modelled energy budget shown in Fig. 4b. The impact of dust on albedo and slope on 27 solar irradiance is well established in the snow literature and so this is expected. 28 What is a more interesting finding here is that dHS was not correlated with initial HS0, 29 13 SCD curves (Pomeroy et al., 2001), which will be highlighted in section 3.5. Figure 5c shows 1 the areal mean values for HS0 and dHS for flat areas (slope < 5 o ) and areas on both sides of 2 the ridge (threshold aspect is 235 o , slope ≥ 5 o ). The hypothesis for this study period was that 3 large drifts on south-facing parts of the ridge cause a correlation between melt energy and 4 SWE. Indeed, the southeast part showed larger HS0 and dHS compared to the flat and 5 northwest part of the study area. This suggests a correlation between HS0 and dHS, which 6 was not apparent when analysing all pixels. In each subarea the range of snow depth was 7 large, which diminished the observed correlation. More importantly, on the south-eastern face 8 a mild negative correlation of r = -0.35 developed (Fig. 5d), which may be explained by a 9 remaining cold content in deep drifts. This negative correlation is not apparent for smaller 10 dHS values, in the northwest part of the ridge (Fig. 5e). The lack of correlation in the Fig. 5c  Given these scenarios some guidelines for modelling areal SCD can be provided. 30 Models must be able to represent realistic correlations between SWE and melt in order to 31 model the effect of this correlation on SCD (Essery and Pomeroy, 2004). Potential pitfalls are 32 14 incomplete modelling representations that might neglect a governing process. To capture the 1 spatial correlations, models need to include snow redistribution, internal snowpack energetics 2 and melt rate variability on slopes at fairly fine scales (<100 m) in complex terrain. Semi-3 distributed models with homogenous snow distribution over large areas or distributed models 4 that neglect blowing snow redistribution may misrepresent spatial correlations of SWE and 5 melt. 6 Another reason for models misrepresenting spatial correlations between HS0 and dHS 7 is discussed in section 3.6, in which the mismatch of scales of dHS and HS0 patterns is 8 discussed. 9 Table 2 shows mean, standard deviation and CV values of HS and dHS in different 11 periods and subareas. Throughout the melt season CV values of dHS were about five times 12 smaller than those of HS. At the start of the study period, the variability of dHS was smaller 13 than that of HS by a factor 3.7 to 6.7. 14 For the whole area only a weak temporal correlation (r = 0.36) was found in a pixel-15 by-pixel analysis between ablation patterns over the two long periods P1 and P2. Larger HS at the start of the study period (HS0) and a spatially uniform dHS value was 30 applied for each pixel. This value was determined with observed mean ablation 1 values shown in Table 3. Each pixel was reduced by this mean value and any 2 negative values in HS were set to 0. SCA was defined as the ratio of the number of 3 grid points with HS > 0 to all pixels. 4 2. Uniform HS0/variable dHS: In this scenario, the mean initial snow depth as shown 5

Variability of dHS in relation to HS0 and temporal persistence 10
in Table 3 was uniformly distributed in the whole snow-covered area. Spatially 6 variable dHS values as measured with the UAV were applied to each pixel. To 7 obtain the exact melt-out time this scenario was calculated daily using a temporally 8 constant dHS value between flights. No exact dHS amounts were available for 9 pixels which melted out between flights. For those pixels the mean areal dHS 10 value was applied. The general shape of SCD curves obtained when this scenario 11 was also calculated on the time resolution of the UAV flights. 12 3. Uniform HS0/uniform dHS: Similar to scenario 2, but a spatially uniform dHS 13 value was applied to each pixel, each of which had a uniform HS0. This scenario 14 was also calculated on a daily resolution . 15 In all scenarios, SCA was set to 1 for the area which was snow-covered at the start of 16 the study period. Figure 6 shows mean HS ablation and SCD curves for the whole area and 17 the northern subarea (top), for which more flights are available. Differences between 18 measured development and the first scenario of uniform dHS and variable HS0 were not 19 large. However, a large difference between measurements and the second and third scenarios 20 with uniform HS0 and either variable or uniform dHS is obvious. Areal dHS in those 21 scenarios was overestimated before modelled melt-out because of the overestimation of SCA. 22 Later, areal dHS was underestimated (or zero) since most or all snow disappeared too early. The main reason why the observed dHS differences, which were substantial and partly 29 persistent, did not influence SCD curves can be found in the small to negligible spatial 30 correlation between dHS and HS0 (cf. section 3.3 and Table 1 Where correlation is insignificant, spatial melt differences can be quite large without 3 affecting SCD curves. In this case, spatially variable melt can be viewed as a nearly random 4 processit introduces noise into the log-normal frequency distribution of HS, but does not 5 affect the emergent behaviour of the SCD curve. Here, with a much larger variability of HS0 6 compared to dHS (see section 3.4) and only small spatial correlations between them (see 7 Table 1), HS0 controls the SCD. variance. Indeed, the correlogram (Fig. 8b) confirms that the largest correlation with dHS to 30 ρ xy = 0.4 was achieved at those larger distances. 31 The same analysis for initial snow depth (HS0) can be seen in Fig. 8c and d. Most of 1 the variance for snow depth is at length scales less than 100 m. The periodic behaviour shown 2 beyond that scale may be due to the patchy snow cover which has long snow-free patches. No 3 substantive correlation with dHS is observable on all scales (Fig. 8d). 4 This analysis offers further explanation why dHS and HS0 were not spatially 5 correlated in these observations. dHS variance was related to large scale aspect changes on 6 both slopes and medium scale albedo change, whilst snow depth was variable mainly at much 7 smaller scales. This scale mismatch leads to a larger scatter between dHS and HS0 values and 8 thus prevented a substantive spatial correlation. 9 Two processes were previously discussed and described in Fig. 5c which could drive 10 compensating correlations between HS0 and dHS; cold content and melt energy. Cold content 11 likely acts on a similar scale as HS0, since it depends mainly on snow depth. As shown in Fig.  12 5d and 5e a negative correlation driven by cold content is not uniformly present. Melt energy 13 differences, i.e. differences in net shortwave radiation, turbulent fluxes, and net longwave 14 radiation, are not directly dependent on snow depth, but need to spatially coincide by chance 15 (e.g. by direction of redistribution). Acknowledging that Solar Irradiance is a simple proxy of 16 melt energy, spatial coincidences between accumulation and melt energy are only present 17 over larger distances (Fig 8b). The large scatter between HS0 and dHS results from the 18 observation that most of the variance of HS0 occurs at much smaller scales (Fig. 8c). Figure  19 8d illustrates variability in the compensating correlations. At small scales below 50 m, the 20 differences in Solar Irradiance are small and the cold content is responsible for slight 21 negative correlation between HS0 and dHS. This is counteracted by Solar Irradiance until the 22 distance of 250 m (cp. Fig 8a). 23 There needs to be a match in scaling behaviour between SWE and melt rate for these 24 variables to develop spatial correlations. Assuming melt is primarily driven by aspect and 25 slope differences as in the proxy Solar Irradiance, SWE must vary on similar scales for a 26 correlation to develop. This may be achieved if SWE varies primarily over larger scales, e.g. 27 in a simple topography of a ridge without gullies and with one predominant wind direction 28 during blowing snow, in which one slope face has much larger SWE values than the other. 29 This may also be achieved if Solar Irradiance acts on a smaller scale similar to HS0. This 30 might be possible in highly complex terrain in which most slope/aspects differences can be 31 found on scales below 100 m but this does not correspond to the "ridge" in our study site. 32

Conclusions and outlook 1
The aim of this study was to determine factors which influence areal snow ablation 2 patterns in alpine terrain using spatially intensive observation. The dependency of snow 3 accumulation and topographic variables with spatial melt rates were analysed for an alpine 4 ridge in the Fortress Mountain Snow Laboratory located in the Canadian Rocky Mountains. 5 Detailed maps of snow depth, snow depth change and snow-covered area were generated 6 during late season ablation with UAV-based orthophotos, photogrammetry and Structure-7 from-Motion techniques. Snow depth and its change served as proxies for snow accumulation 8 and melt rates. Snow depth change values were found to be spatially variable and mainly 9 dependent on variation in solar irradiance and albedo, and likely on the cold content of the 10 snowpack which is a function of snow depth. Local and small-scale dust variations, which had 11 not previously been observed in the area, increased the variability of ablation. 12 Snowcover depletion curves were mostly dominated by the variability of initial snow 13 depth at the start of this study rather than the variability in snow depth change. Initial snow 14 depth variability was approximately five times larger than the variability in snow depth 15 change in this windswept environment. The scales of variability of snow depth and snow 16 depth change were mismatched, with snow depth variability occurring at small scales (<10 m) 17 and snow depth change associated with the medium scale (50 m) of albedo variation or the 18 slope scale (100s of m) of solar irradiance variation. As a result, the initial snow depth and 19 changes in snow depth were not strongly correlated over space, and so only initial snow depth 20 influenced snowcover depletion. 21 The observations collected here show the prerequisites for strong correlations that can 22 impact snowcover depletion curves. Correlation between melt and snow accumulation may 23 be driven by cold content and melt energy distributions. Whilst cold content can create a 24 negative correlation between melt and snow accumulation, melt energy variations can create 25 either positive or negative correlations. In order to not compensate for each other, one process 26 needs to be dominant, or the both processes need to create a similar negative correlations. It 27 is also important that these variations occur at the same spatial scales. 28 To further investigate these arguments, longer time series of spatially detailed 29 snowpack and snowcover observations need to be made in order to further test and examine 30 the temporal evolution of the spatial covariance and variance of ablation and accumulation in 31 various global alpine environments. The results of such a study could suggest how to 32