Topographic control of snowpack distribution in a small catchment in the central Spanish Pyrenees : intra-and interannual persistence

Introduction Conclusions References


Introduction
Assessing the snow distribution in mountain areas is important because of the number of processes in which snow plays a major role, including erosion rates (Pomeroy and Gray, 1995), plant survival (Keller et al., 2000;Wipf et al., 2009), soil temperature and moisture (Groffman et al., 2001), and the hydrological response of mountain rivers (Bales and Harrington, 1995;Barnett et al., 2005;Liston, 1999;Pomeroy et al., 2004).As mountain areas are highly sensitivity to global change (Beniston, 2003), snow accumulation and melting processes are likely to be subject to marked changes in coming decades, affecting all processes influenced by the presence of snow (Caballero et al., 2007;López-Moreno et al., 2011, 2012b;Steger et al., 2012).For these reasons, much effort has been devoted to understanding the main factors that control the spatial and temporal dynamics of snow (Egli et al., 2012;López-Moreno et al., 2010;Mott et al., 2010;Schirmer et al., 2011).
One of the main difficulties in snow studies is obtaining reliable information of the variables that describe snow distribution, including snow depth (SD), snow water equivalent (SWE) and snow covered area (SCA).Manual measurements Published by Copernicus Publications on behalf of the European Geosciences Union.
have traditionally been used to provide information on the distribution of snowpack, with different sampling strategies having been applied at various spatial scales (Jost et al., 2007;López-Moreno et al., 2012a;Watson et al., 2006).However, manual sampling is not feasible for large areas because of the time involved, especially when SWE measurements are also acquired.In the last decade the use of airborne laser scanners (ALSs) (Deems et al., 2006) and terrestrial laser scanners (TLSs) (Prokop, 2008), both of which are based on lidar (light detection and ranging) technology, have provided for major advances in obtaining data on the SD distribution at unprecedented spatial resolutions.These developments have enabled studies of several factors that in the past have been only marginally considered, including scaling issues (Fassnacht and Deems, 2006;Mott et al., 2011;Schirmer and Lehning, 2011;Trujillo et al., 2007), the detailed dynamics of snow accumulation and ablation (Grünewald et al., 2010;Schirmer et al., 2011;Scipión et al., 2013), and snow transport processes (Mott et al., 2010).In addition, the high density measurements provided by lidar technologies are a valuable resource for detailed investigation of the linkage between snow distribution and topography.In the past, this linkage has mostly been studied using manual measurements, and hence with generally limited spatial and temporal resolution (López-Moreno et al., 2010).
Previous studies have highlighted the marked control of topography on snow distribution in mountain areas (Anderton et al., 2004;Erickson et al., 2005;Lehning et al., 2011;Mott et al., 2013), and the importance of vegetation and wind exposure (Erxleben et al., 2002;Trujillo et al., 2007).The most commonly used approach has been to develop digital elevation models (DEMs) that describe the spatial distribution of elevation, from which other terrain variables are derived such as slope, terrain aspect, curvature, wind exposure or sheltering, and potential solar radiation.This enables one to analyze the linear or nonlinear relation of these variables to punctual SD or SWE values to be established (Grünewald et al., 2010;Schirmer et al., 2011).Various statistical methods have been applied for this purpose, including linear regression models (Fassnacht et al., 2003;Hosang and Dettwiler, 1991), generalized additive models (GAMs) (López-Moreno and Nogués-Bravo, 2005), and binary regression trees (BRTs) (Breiman, 1984) which have been widely applied in a diversity of regions (Elder et al., 1991;Erxleben et al., 2002;McCreight et al., 2014).
The extent to which topographic variables explain snow distribution can change during the snow season; the variability of terrain characteristics can drive processes related to the spatial variability of snow accumulation (snow blowing, terrain curvature) (Lehning et al., 2008), or affect the energetic exchange between terrain and the snowpack (temperature, incoming solar radiation), so the importance of topographic variables is modified during the season (Molotch et al., 2005).In addition, during a snow season the terrain changes markedly (is smoothed) by snow accumulation (Schirmer et al., 2011).However, few studies have systematically analyzed the intra-and inter-annual persistence of the relation between snow distribution and topography.Recent studies have assessed whether the influence of topography is constant among different years, e.g., the similarities observed at the end of the accumulation season (Schirmer and Lehning, 2011;Schirmer et al., 2011), or the consistent fractal dimensions in two analyzed years (Deems et al., 2008); in both cases there was a relation with the dominant wind direction, which highlights the predictive ability of topographic variables.
The main focus of this study was to assess the influence of topography on the spatial distribution of snowpack and its evolution over time.The high temporal and spatial density of the data set collected during the study enabled analysis of the main topographic factors controlling snow distribution, and assessment of whether topographic control of the snowpack varied during the snow season and between years having very contrasting climatic conditions.For this purpose, we conducted 12 surveys over 2012 (6) and 2013 (6) in a small mountain catchment representing a typical subalpine environment in the central Spanish Pyrenees, and obtained highresolution SD measurements using lidar technology with a TLS.

Study area and snow and climatic conditions
The Izas experimental catchment (42 • 44 N, 0 • 25 W) is located in the central Spanish Pyrenees (Fig. 1).The catchment is on the southern side of the Pyrenees, close to the main divide (Spain-France border), in the headwaters of the Gállego river valley, and ranges in elevation from 2000 to 2300 m above sea level.The catchment is predominantly east-facing, with some areas facing north or south, and has a mean slope of 16 • .There are no trees in the study area, and the basin is mostly covered by subalpine grasslands dominated by Festuca eskia and Nardus stricta, with rocky outcrops in the steeper areas.Flat, concave and convex areas occur in the basin.
The climatic conditions are influenced by the proximity of the Atlantic Ocean, with the winters being humid compared with zones of the Pyrenees more influenced by Mediterranean conditions.The mean annual precipitation is 2000 mm, of which snow accounts for approximately 50 % (Anderton et al., 2004).The mean annual air temperature is 3 • C, and the mean daily temperature is < 0 • C for an average of 130 days each year (del Barrio et al., 1997).Snow covers a high percentage of the catchment from November to the end of May.
The two years analyzed in the study represent climatic extremes during recent decades.Severe drought occurred during 2012, leading to snow accumulation well below the longterm average.The thickness of the snowpack, measured at the automatic weather station (AWS, Fig. 1), during winter 37 10.Figures  in this year was less than the 25th percentile of the available historical data series of this AWS (1996-2011) (Fig. 2).Only at the end of spring did late snowfall events increase the amount of snow, but this rapidly melted.The opposite occurred in 2013, a year in which the deepest snowpack and the longest snow season of recent decades were recorded.Winter and spring in 2013 were extremely humid, with temperatures mostly between the 25th and 75th percentiles of the AWS historical series.SD accumulation was very high between February and June (exceeding the 75th percentile).
In some areas of the basin snow lasted until late July, which is one month longer than in most of the preceding years for which records are available.Regarding net solar radiation data (shortwave), no measurements were available before December 2011, However, the annual evolution has been tracked on Fig. 2   show the TLS survey days.Note that during some surveys no snow was present at the AWS, but some areas of the Izas experimental catchment were covered by snow.(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011).The vertical dashed lines show the TLS survey days.Note that during some surveys no snow was present at the AWS, but some areas of the Izas experimental catchment were covered by snow.

Snow depth measurements
During the study period high-resolution SD maps (Fig. 3) were generated using a long range TLS (Riegl LPM-321), which enables safe acquisition of SD information with short acquisition times from remote areas, compared with measurements obtained manually.This technique has been extensively tested (Prokop et al., 2008;Revuelto et al., 2014;Schaffhauser et al., 2008), and systematically applied to the study of snow distribution in mountain terrain (Egli et al., 2012;Grünewald et al., 2010;Mott et al., 2013;Schirmer et al., 2011).In a previous study, the mean absolute error in the most distant areas of the catchment was less than 10 cm (Revuelto et al., 2014), which is consistent with errors reported in previous studies (Grünewald et al., 2010;Prokop, 2008;Prokop et al., 2008;Schaffhauser et al., 2008).
TLS provides high-resolution three-dimensional information on the terrain.Nevertheless, error sources need to be considered because they can have large effects on the measurements.To reduce the influences of instability on the TLS (originated through small displacements of the tripod be-cause the TLS vibrates while it is operating), which leads to misalignment with reference points and atmospheric change, a well-defined protocol must be applied.The protocol applied in this study for generating high-resolution SD maps with a 1 m cell size was described by Revuelto et al. (2014).This protocol is based on the following main points: data collection; which includes experimental setup design and information acquisition by the scanning procedure; and data processing, when data is filtered, quality checked and the SD maps generated.Mainly, the methodology was based on differences between DEMs obtained with snow coverage in the study area and a DEM taken at 18 July 2012, when the catchment had no snow cover.Twelve SD maps at a spatial resolution of 5 m were generated for the 2012 and 2013 snow seasons (Fig. 3).In each year three surveys were undertaken from February to April (2012: 22 February, 2 April, 17 April;2013: 17 February, 3 April, 25 April), and three were undertaken from May to June when intense melting conditions dominated (2012: 2, 14 and 24 May; 2013: 6, 12 and 20 June).The average SD and SCA, and the maximum SD are shown in Table 1.It shows that much lower SD and SCA were observed in 2012 compared to 2013.

Digital elevation model and topographic variables
From the two scan stations located in the study area (Fig. 1), 86 % of the total area of the catchment was surveyed using the TLS.DEMs of 1 m grid size were initially obtained from point clouds of varying density in different areas, but always with a minimum of 1 point m −2 (Revuelto et al., 2014).Some of the predictor variables cannot be calculated where data gaps occur in the DEM (e.g., the topographic position index), and others require a DEM with a greater surface than the area scanned during the study (e.g., to calculate the potential solar radiation, including the shadow effect from surrounding topography, or to calculate the maximum upwind slope parameter, topographic information for distances up to 1200 m is included from the exterior limit of the DEM obtained with the TLS).Thus, a DEM having a 5 m grid size, available from the Geographical National Institute of Spain (Instituto Geográfico Nacional, http://www.ign.es), was combined with the snow-free DEM obtained using the TLS resampled from 1 to 5 m resolution (the empty raster of the Geographical National Institute was used for the resampling, averaging all values within each cell).The 1 m grid-size SD maps were also resampled to 5 m to enable matching of the two different data sources.
Elevation was obtained directly from the DEM, while the other variables were calculated using ArcGIS10.1 software.This software calculates Slope as the maximum rate of change in value from a specific cell to that of its neighbors (15 m × 15 m window size), and Curvature is determined from the second derivative of the fitted surface to the DEM in the direction of maximum slope of the terrain for the neighbors cells (15 m × 15 m window size too).Radiation was obtained using the algorithm of Fu and Rich (2002) and reported in Wh m −2 based on the average conditions for the 15day period prior to each snow survey.This algorithm calcu-lates the potential incoming solar radiation (shortwave) under clear sky conditions, which may strongly differ from the real radiation as a consequence of cloud cover.This measure provided the relative difference in the extraterrestrial incoming shortwave solar radiation among areas of the catchment for a given period under given topographical conditions (Fassnacht et al., 2013).In this way, Radiation can be considered as a good proxy of the spatial distribution of incoming solar energy within the study area.Easting and Northing exposure were calculated directly as the sine and cosine, respectively, of the angle between direction north and terrain orientation or aspect.It provided information on the east (positive)/west (negative) exposure and the north (positive)/south (negative) exposure.
The TPI provides information on the relative position of a cell in relation to the surrounding terrain at a specific spatial scale.Thus, this index compares the elevation of each cell with the average cell elevation at specific radial distances as follows (De Reu et al., 2013;Weiss, 2001): where z o is the elevation of the cell in which TPI is calculated and z is the average elevation of surrounding cells obtained from Eq. ( 2) for a radial distance R. For each pixel the TPI was calculated for 5, 10, 15, 25, 50, 75, 100, 125, 150 and 200 m radial distances (scale factors).
For specific wind directions, the maximum upwind slope parameter, averaged for 45 • upwind windows ( Sx; Winstral et al., 2002) provided information on the exposure or sheltering of individual cells at various distances, resulting from the topography.Rather than considering the contribution for the dominant wind directions (Molotch et al., 2005), Sx(Sx further on) values for eight directions were selected and directly related to the SD.The directions were 0 • for north (N), 45 • for northeast (NE), 90 • for east (E), 135 • for southeast (SE), 180 • for south (S), 225 • for southwest (SW), 270 • for west (W) and 315 • for northwest (NW).For Sx, the searching distances (Winstral et al., 2002) considered were 100, 200, 300 and 500 m.These distances were selected to enable assessment of the range at which Sx exhibited greatest control on SD dynamics, as has occurred in previous studies (Schirmer et al., 2011;Winstral et al., 2002).

Statistical analysis
The 12 SD maps at 5 m spatial resolution were related to each of the topographic variables considered (including the 40 Sx combinations, and the 9 distances for TPI).The large number of cells for which SD data were available enabled robust correlations between topography and snow distribution to be obtained, and provided a very large data set for training and validating the SD distribution models.
Pearson's r coefficients were obtained between SD and each topographic variable.Using the whole data set, each variable was correlated, for all available points, with the SD value for the specific survey day.Given the large amount of data for surveys, the degrees of freedom for correlation analyses were very high and hence they can account for statistically significant correlations even with very low correlation coefficients.Moreover, the use of a very dense data set of observations may have associated problems derived from spatial autocorrelation (Koenig, 1999).For reducing effects derived from spatial autocorrelation, we followed a Monte Carlo procedure, in which 1000 random samples of 100 SD cases were extracted from the entire data set (an average of 20 000 SD measurements for each day) and correlated with topographic variables for assessing significance.A threshold 95 % confidence interval (α < 0.05) was used to assess the significance of correlations (r = ±0.197,based on 100 cases).The spatial scales of Sx and TPI for which SD showed a higher correlation, 200 and 25 m respectively, were selected for further analysis (not presented in the manuscript).
To assess the explanatory capacity when all topographic variables were considered simultaneously, two statistical models were used: (1) multiple linear regressions (MLRs) and (2) binary regression trees (BRTs).A wide variety of regression analyses for interpretation of much more complex spatial data are available with greater capacity than MLRs and BRTs to deal with spatial autocorrelation issues and the nonlinear nature of the relationship between predictors and the response variable (Beale et al., 2010).However, in this study we used MLRs and BRTs because these methods have been and are still widely used in snow studies, and because both enable to isolate accurately the weight of each independent variable within the model, which was the main objective of this research, rather than deriving models with maximum predictive capacity.Prior to running the models, a principal component analysis (PCA) was applied to the entire data set for detecting correlations between independent variables that could originate multicollinearity in MLR and BRT.This analysis (not shown) grouped the topographic variables in three components, showing that TPI and Curvature are highly correlated with PCA component one, and also Northing and Radiation (but in this case with opposite signs) presented high correlation with component two of the PCA.TPI and Northing showed both higher correlations with their respective components and in general higher Pearson's r coefficients with SD than Curvature and Radiation (see result section).Therefore, Curvature and Radiation were discarded as predictors in MLR and BRT analyses.
1. Multiple linear regression estimates the linear influence of topographic variables on SD.Despite its simplicity and the rather limited capability under nonlinear conditions (López-Moreno et al., 2010), MLR was used to quantify the relative contribution of each variable to the entire SD distribution model.SD was calculated from the topographic variables at a specific location and day.
The threshold for a variable to enter in the model was set at α < 0.05.Beta coefficients (obtained dividing the standardized units of the coefficients by the mean value of each variable) were used to compare the weight of each variable within the regression models.As a first step, a reduced data set (1000 cases) was randomly extracted to avoid an excessive number of observations that may lead to spurious identification of statistically significant predictor variables.A stepwise procedure was used to obtain these variables from the random extraction.The variables determined for each survey were used to obtain the final model, but using the entire data set (except 5000 cases for model validation), forcing variables entrance in models.
2. Binary regression trees have been widely used to model snowpack distribution from topographic data (Erxleben et al., 2002;Molotch et al., 2005).These are nonparametric models that recursively split the data sample, based on the predictor variable that minimizes the square of the residuals obtained (Breiman, 1984).One BRT was created for each sampling date.The BRTs were run until a new split was not able to account for 1 % of the explained variance, or when a node had less than 500 cases; a maximum of 15 terminal nodes was set to reduce tree complexity.As there were no overfitting problems associated with sample size, 15 000 cases were used to grow the trees and 5000 cases were used for validation.By scaling the explained variance of each variable introduced into each BRT (based on the % of the total explained variance by the BRT), we were able to compare the relative importance of each topographic variable between the different models.
Coefficients of determination (r 2 ) and Willmott's D statistic were used to assess the ability of each model to predict SD over an independent random sample of 5000 cases.Willmott's D was determined using Eq.(3) (Willmott, 1981) where N is the number of observations, O i is the observed value, P i is the predicted value, and Ō is the mean of the observed values.The index ranges from 0 (minimum) to 1 (maximum predictive ability).

Single correlations
Table 2 shows the correlation between SD and Sx for the eight wind directions at a distance of 200 m (identified as the best correlated searching distance in previous analysis).
Despite differences in magnitude, the correlations for surveys carried out at the beginning of the season (22 February 2012 and 17 February 2013) in each year showed that SD was clearly affected by N and NW wind directions.The contribution of N and NW wind directions is clearly evident for the surveys on 17 February 2013 (Fig. 4, where wind roses with average wind speeds and direction, for the 15-day period before each survey, are presented), when greater SD was recorded in the leeward slopes from a northerly direction (Fig. 3, northerly areas of the maps).In the two years of the study, a correlation with W and SW wind directions was observed to increase progressively during the snow season (Fig. 4 and Table 2 correlations).In 2013 this phenomenon was less marked because of the greater SD accumulation at the beginning of the snow season accompanied with NW direction winds, which resulted in only moderate changes in the Sx for the most strongly correlated wind directions.It was also observed that in both study years once the snow had started to melt (the last three surveys in each season) the snow distribution did not change in relation to Sx directions.The best correlated Sx directions for each survey are in good agreement with wind roses main directions (Fig. 4).These directions are for the following survey days: 315 Correlations between the best correlated Sx direction for each day and SD were compared with correlations between SD and the other topographic variables (Table 3).This showed that Sx had one of the greatest coefficient of correlation with SD (range 0.22-0.56).The correlations were higher during the accumulation periods, especially in the 2013 snow season, with a reduction in correlations values occurring during the melt period at the end of each snow season.
The TPI at 25 m showed the highest correlation with SD for the 12 sampled days.During 2012 the mean correlation values ranged from −0.32 to −0.58 for those surveys during which snow accumulation dominated in the days preceding the surveys.The r values were closer to the significance level for the surveys where the preceding days were dominated by melting conditions (14 and 24 May).In 2013, the TPI was more highly correlated with SD than in 2012, with Pearson's r coefficients < −0.6 for all survey days.Curvature also had a high correlation with SD, and similar to TPI with a 25 m searching distance was significantly correlated on all the survey dates, but unlike the TPI, the correlation of Curvature with SD did not decrease during the snowmelt periods.The significant correlations of TPI and Curvature with SD highlight the importance of terrain curvature on the SD distribution.The importance of terrain curvature at different scales for SD distribution is clearly evident in Fig. 3, which shows that higher SD values were usually found for concave areas, which showed snow presence until the end of each snow season.
The correlation between Elevation and SD varied among survey days (Table 3).The correlations were usually positive, but only statistically significant (or approaching significance) for days when melting dominated (the last two surveys in 2012 and 2013).Slope was relatively weakly correlated with SD during the 2012 snow season.In 2013 the correlation was greater, and was statistically significant for all days.Similarly to Elevation, the correlation between Slope and SD was variable between the two study years, and showed a similar temporal pattern to Easting, probably because of the presence of steeper areas on the east-facing slopes.
The correlation between Northing and SD was rarely statistically significant, highly variable and contributed to explaining SD in very different ways in 2012 and 2013.In 2012 no correlation between SD and Northing was found during the accumulation period, but during the melting period a slight positive correlation was observed, as snow remained longer on north-facing slopes.The 2013 snow season started with a large precipitation event dominated by strong winds from a northerly direction, leading to high levels of snow accumulation on the south-facing slopes.This explains the strong and statistically significant negative correlation of SD with Northing for 17 February 2013.This event influenced the rest of the season (as evident in Table 2 for 2013), but a progressive decrease in its influence was evident for the following survey days.Radiation had an almost opposite influence on SD to that observed for Northing.During the melting period, for each year, the Pearson's r correlation between SD and Radiation was negative, indicating a thinner snowpack on the most irradiated slopes.This relation was statistically significant at the end of the 2013 snow season.However, during the accumulation period in 2013, statistically significant and positive correlations were observed with Northing and Radiation, which are connected to the strong snow redistribution by winds from N-NW directions.

Multiple linear regression and binary regression tree models
Figure 5 shows the Willmott's D values and the coefficients of determination (r 2 ) obtained in the comparison of observed and predicted values using MLRs and BRTs for a data set reserved for validation (5000 cases).The MLRs produced r 2 values ranging from 0.25 to 0.65 and Willmott's D values ranging from 0.60 to 0.88, while the BRTs produced r 2 values ranging from 0.39 to 0.58 and Willmott's D values ranging from 0.72 and 0.85.For both methods the relations between the observed and predicted values were stronger for 2013.Accuracy decreased at the end of the snow season, when the snowpack was mostly patchy across the basin; this was particularly the case for the end of the 2012 season.
Overall, the performance of the MLRs was more variable than that of the BRTs, which were more constant amongst the various snow surveys.For those days on which the models were most accurate in predicting SD variability, the MLRs showed slightly better scores than the BRTs.However, for days on which the accuracy between predictions and observations was lower, the BRTs provided better estimates than the MLRs.For 2012, slightly better results were obtained using MLRs, while the opposite occurred in 2013.Neverthe-less, only large differences in the accuracy of each model were evident by the end of 2012 snow season, in the two last surveys, which were characterized by thin and patchy snowpack.
As shown for single correlations, the TPI variable explained most of the variance in MLR models developed for all analyzed days (Table 4).The contribution of the other variables varied markedly among surveys, particularly when the two years were compared.In most cases, Elevation was the second most important variable explaining the SD distribution in 2012, followed by Sx and Slope.The other variables made a much smaller contribution, or were not included in the models.The contribution of Elevation was much less in 2013, and it was not included in three of the six surveys, whereas in 2012 it was included in all surveys.For the entire 2013, Sx was the second most important variable, followed by Easting, which had an almost negligible influence in 2012.Northing was only included in the models for the surveys carried out during periods dominated by snow accumulation, and was not included in the models during the periods dominated by melting.
Figure 6   The relative importance (scaled from 0 to 100) of each topographic variable in each BRT is shown in Table 5.This shows that TPI was the first most important variable explain-ing SD for all survey days.For the 2012 snow season, TPI explained more than 67 % of the total explained variance in all BRTs, and 75 % during the accumulation period (the first three surveys).Thus, for most of the survey days the variance explained by the other variables was < 30 %. Furthermore, TPI was in all cases the first split variable (which accounted from a 23 to a 30 % of the explained variance), with a critical value that ranged from −0.67 m to −0.The final nodes (with ellipses) show the predicted SD in the zone having the specified terrain 859 characteristics.At each branch point, one topographic variable is considered; if the value is 860 less than the specified value, the left branch is selected, but if it is equal to or greater than the 861 specified value, the right branch is selected.The final nodes (with ellipses) show the predicted SD in the zone having the specified terrain characteristics.At each branch point, one topographic variable is considered; if the value is less than the specified value, the left branch is selected, but if it is equal to or greater than the specified value, the right branch is selected.
contributor to the total explained variance, exceeding 50 % for almost all survey days, and approaching or > 70 % during the snowmelt period.For this year, also TPI was the first split variable in nearly all BRT, with critical values ranging from −0.47 to −0.15 m and an average value of −0.28 m; except for the 13 February 2013, in which Sx was the first split variable.The influence of Sx was more important in 2013 than in the previous year.At the beginning of 2013 the contribution of Sx to the total explained variance was almost 46 %, and remained > 20 % for the rest of the snow season.An exception was the last survey, when melting dominated and its effect declined to 12 %.When snow was not mobilized for long periods by wind (no changes on the best correlated wind direction of Sx are observed), the SD distribution was more dependent on variables related to terrain curvature (TPI and Curvature).During 2013, Elevation contributed approximately 5 % to the total explained variance during the entire snow season.Northing made a significant contribution to the model (14.7 %) only one day (3 April 2013), and a smaller contribution on the following survey day (25 April 2013).
When included in the BRTs, the other variables (Easting, Radiation) made minor contributions to the total explained variance.
Figure 7 shows the mean contribution of each topographic variable vs. the coefficient of variation from the twelve surveys for the different statistical approaches considered in this study (Pearson's r coefficients, beta coefficients of the MLRs and the contribution to the explained variance for each BRT).Clearly, TPI is the most important variable to explain the snow distribution in the catchment, but it is also the variable that exhibits a lower variability between the different surveys (CV < 0.2).Besides it has been introduced as predictor for MLRs and BRTs in all studied days.Sx is the next variable in importance to explain snow distribution, being introduced as predictor in the majority of the modelled days (11 and 10 out of 12 days for MLRs and BRTs, respectively).It shows a low temporal variability when correlation's coefficients are calculated (CV = 0.24), but the variability in its contribution to MLRs and BRTs increases noticeably, with CV values of 0.35 and 0.59, respectively.The rest of the variables show a much lower mean contribution for explaining snow distribution and a high temporal variability in their explanatory role.Lower CV values are observed for MLRs, ranging the majority between 0.3 and 0.4, than for BRTs models, ranging the majority between 0.4 and 0.8.

Discussion
The distribution of snow in mountain areas is highly variable in space and time, as shown for Izas experimental catchment during two consecutive years.Many meteorological and topographic parameters affect the snow distribution and its evolution through time with different weights subjected to several factors.In this context, we demonstrated that terrain characteristics significantly affect SD distribution in a subalpine catchment.Also we have shown that its effect evolved during the snow accumulation and melting periods over two years having highly contrasting climatic conditions and snow accumulation amounts Many studies have analyzed the spatial distribution of SD in mountain areas considering both, intra-and inter-annual variability of the topographic control on the snowpack distribution (Anderton et al., 2004;Erickson et al., 2005;López-Moreno et al., 2010;McCreight et al., 2014).Other researches have also focused their attention in long-term interannual snow distribution analyses (Jepsen et al., 2012;Sturm and Wagner, 2010;Winstral and Marks, 2014).The results of these previous works have highlighted the difficulties in fully explaining the distribution of snow in complex mountainous terrain.In addition, results differ among studies, and suggest that different variables govern the distribution of snowpack among areas as consequence of their different characteristics and geographical settings.These differences include surface extension, the altitudinal gradients, the importance of wind redistribution, the presence or absence of vegetation and the topographic complexity as concluded by Grünewald et al. (2013) in a study where seven study sites across the world were considered.
Most of the topographic variables investigated in this study have been included in previous studies, including Elevation, Slope, Radiation, Curvature and Sx.Other variables, in particular TPI, have received little attention in previous research (López-Moreno et al., 2010).We showed that TPI at a scale of 25 m had the greatest capacity to explain the SD distribution in the study catchment.Curvature (which refers to a smaller spatial scale of terrain curvature when compared with TPI) is also highly correlated with the SD distribution, but not as much as TPI.This reinforces the importance of considering terrain curvature at various scales for explaining the SD distribution in mountain environments.The correlation between snowpack and the TPI decreased during melting periods, whereas the correlation with Curvature remained constant.This suggests that snow accumulates more in small deep concavities, but is shallower at the end of the season in wider concave areas that were identified by the 25 m TPI scale.This effect was evident at the end of the snow season, when snow was present only in deep concavities, as shown in Fig. 3. To explain the snow distribution, Anderton et al. (2004) compared the relative elevation of a cell with the terrain over a 40 m radius, and observed that this had a major role on SD distribution, which sustain curvature importance at different scales.
The maximum upwind slope (Sx; Winstral et al., 2002) has also been identified as a key variable explaining snow distribution, improving the results obtained when it is introduced into models.Our results, 200 m searching distance for Sx, is comparable with those of other studies that have shown that the optimum searching distance for correlating Sx with the SD distribution is 300 m (Schirmer et al., 2011), which is not a large difference for the considered distances in this work reaching 500 m.As it is observed from the reported wind information, Izas experimental catchment has W-NW dominant wind direction what is consistent with the best correlated Sx directions.For this reason, the Sx preferred direction for each date was selected, and showed that there were intra-annual shifts in the most highly correlated direction.The change in the most important Sx direction was similar between the 2012 and 2013 snow seasons; it started with a northerly component and evolved to a dominant westerly direction.We also found a decrease in the correlation between Sx and the snow distribution at the end of each snow season, Each graph point is accompanied with its variable, and in the case of MLR and BRT, in brackets the number of days in which each variable was included in models.when melting conditions dominated.This is consistent with the findings of previous studies (Winstral and Marks, 2002).
Sx parameter takes into account sheltering effects with topographic origin in relation to wind directions.SD distribution maps show higher SD amounts in leeward slopes, located in E-SE slopes.TPI is not able to explain snow drifts, because this index considers the topographic characteristics in all directions.Nevertheless, terrain characteristics at the study site in relation to SD distribution have shown a higher importance of TPI when compared to Sx.The most plausible explanation accounting for this result is that the basin has a rather reduced size, shows the same general aspect (SE facing) and topography is relatively gentle.Under such conditions, during wind blowing events snow is accumulated in all wide concavities of the basin (represented by TPI) independently of its specific location.Nonetheless, wind redistribution will be affected by a combination of local topography and main wind directions; which makes it necessary to consider the Sx parameter.As it has been observed, this effect lasts in time until the melting season has advanced.
Only for two days (22 February 2012 and 2 April 2012) was there no contribution (or it was minor) of Sx to explain SD distribution, according to BRTs and MLRs.On these days Northing was introduced into the models, and was found to explain some of the variance of Sx from the northerly direction, the best correlated direction for these days (Table 2).
Although Elevation has been found to largely explain the snow distribution in areas having marked altitudinal differences (Elder et al., 1998;Erxleben et al., 2002;Molotch and Bales, 2005), in our study no strong association was found between SD and Elevation, with significant correlations occurring only during the snowmelt period.This is because of the low elevation range of the study area (300 m).During the accumulation period the entire catchment is generally above the freezing height.However, during spring the 0 • C isotherm shifts to higher elevations, which may lead to different melting rates within the basin.Despite the relatively weak correlation between Elevation and SD, this variable was introduced as a predictor in the MLRs and BRTs for most of the days analyzed.Similarly, López-Moreno et al. ( 2010) reported that elevation was of increasing importance as the grid size increased.Anderton et al. (2004) also informed about the importance of elevation to explain snowpack distribution in the same study area.The results of the present study suggest the increase in importance of Elevation at the end of the snow season, and particularly when it is considered in combination with other topographic variables in MLR and BRT models.
Slope has a weak explanatory capacity for snow distribution, probably because the slope in most of the catchment is not steep enough to trigger gravitational movements including avalanches and slush during the snowmelt period, which could thin the snowpack on the steepest slopes (Elder et al., 1998).Most likely, some of Slope explanatory capacity is included on Radiation explanatory capacity, because it affects solar light incident angle, besides the steeper areas of the catchment are in south facing zones.Nevertheless, quantifying such kind of effects is highly difficult due to the high complexity of SD dynamic in mountain terrain.
Radiation, Northing and Easting showed no close correlation with the snowpack distribution; their relationships with SD were variable over time, with statistically significant correlations occurring on some days and only weak correlations on other days.The results suggested that Radiation and Northing (which showed almost opposite patterns) may be related to SD for two different reasons.During the accumulation period in 2013 heavy snowfalls associated with northerly winds led to the accumulation of deep snow on south-facing areas (more irradiated), whereas during the snowmelt period the greater exposure of the southern slopes to solar energy led to a positive (negative) correlation with Northing (Radiation).This phenomenon was also observed by López-Moreno et al. (2013), using a physically based snow energy balance model in the same study area.Moreover, the high and opposite correlation between Northing and Radiation obtained in PCA results (not shown in the manuscript), was showing a potential problem of multicollinearity.Thus, only Northing was considered for MLRs and BRTs (the same occurred with TPI and Curvature, being only considered in statistical models the TPI).Although Northing did not show a significant correlation with SD during accumulation periods; when the surveys were closer to the snowmelt period, the negative correlation of this variable with SD was more evident, possibly due to the increase of the difference in the energetic exchange between sun exposed and shaded areas.The importance of Northing in MLR models, combined with the contribution of Easting during the accumulation period may be related to the high snow redistribution originated by wind directions from N-NW directions.Therefore, terrain aspect relation with SD distribution (considered with Northing and Easting) in winter is tightly related to the accumulation patterns resulting from wind redistribution, whereas in spring they were associated with the unequal distribution of solar radiation, which leads to higher melting rates on the most irradiated slopes, which has shown better explanatory capacity than Radiation at Izas Experimental catchment.
The MLRs and BRTs provided reasonably high accuracy scores when observed and predicted SD data were compared.The scores were comparable, and in some cases better, to values reported in previous researches using similar methods.As an example, Molotch et al. (2005) reported r 2 values between 0.31 and 0.39 using BRT; and Winstral et al. (2002), who considered different number of terminal nodes of BRT, obtained an optimal tree size of 16 nodes data set with an r 2 value close to 0.4.Moreover results presented here were obtained from a separate data set, and data used to create the models were not considered for testing, thanks to the large available data set.One reason for the improvement may be the use of the TPI as a SD predictor, as this variable has not been considered in previous studies.Nevertheless, it should be noted that the study sites considered in other studies, could differ in terms on complexity of terrain, and also in SD accumulation amounts.For the 12 survey days the TPI had the greatest explanatory capacity in both approaches.However, based on comparison of the different dates and surveys, the other variables made more varying contributions, as a result of the different roles they play during the snow accumulation and melting periods, and the wind conditions during the main snowfall events.The models had less capacity to explain spatial variability of the snowpack when the snow was thinner and patchy.The BRT and MLR approaches were consistent with respect to error estimates.The results obtained using each approach were comparable, so the trends in the variable ranking with both models for each survey day were similar.Only during conditions of snow scarcity did the BRT approach demonstrate better capability to relate SD to topography.This is probably a consequence of the greater capacity of BRTs to take account of the nonlinear response of the snowpack to topography, and the occurrence of sharp thresholds typical of days when the snowpack is patchy (López-Moreno et al., 2010;Molotch et al., 2005).
Despite model results differ between survey days and years, the most important variable, TPI, is always present in the models and their contribution to the total explained variances show very low CV values.Other variables with an also important role to explain SD distribution (i.e., Sx) are included in most of the models as predictors showing their influence on snowpack distribution, although their contribution to the final models changes noticeably amongst different surveys.Moreover for 2012 and 2013 a consistent inter-annual distribution of the snowpack in the catchment is observed; the areas of maximum SD and the location of snow free zones were consistent between both years of the study, and more importantly there is a strong consistency of the effect of topography on SD.This spatial consistency of snowpack has implications for soil dynamics and plant cycles, because some parts of the basin will tend to remain free of snow cover during longer periods favoring the presence of temporary frozen soils, and reducing the isolation effect of snowpack to the plants (Keller et al., 2000;Pomeroy and Gray, 1995).Besides, it suggests that the information acquired from TLS during several years could be useful to design long-term monitoring strategies of SD in the basin based on few manual measurements in representative points according their terrain characteristics.

Conclusions
The TPI at a 25 m searching distance was the best topographic variable, and the most persistent in time, for explaining SD distribution in the Izas experimental catchment.This suggests the importance of including this index in future snow studies, and the need to establish the best searching distance for relating this variable to SD distribution at other study sites.The maximum upwind slope (Sx) at a searching distance of 200 m was also an important variable explaining the SD distribution, but its influence varied markedly between years and surveys, depending of the specific wind conditions during and after main snowfall events.Nevertheless, Sx has shown a similar evolution pattern for the best correlated direction in the two analyzed snow seasons.The influence of the other topographical variables on the spatial distribution of SD was lower, and showed higher intraand inter-annual variability.The total variance explained by BRTs and MLRs clearly decreased for periods on which the snowpack was thinner and more patchily.The results from BRTs and MLRs models were consistent in terms of variables importance ranking, and the explanatory capacities of the main variables were similar for all surveys.Except TPI, that showed very low coefficient of variations for the two approaches, the variability of the contribution of each to-pographic variable for the different surveys was noticeably lower for MLRs than for BRTs.

Figure 1 :
Figure 1: Location of the Izas experimental catchment, and the digital elevation model

Figure 1 .
Figure 1.Location of the Izas experimental catchment, and the digital elevation model showing the positions of the scan stations and the automatic meteorological station.The two images in the bottom part of the figure, from Scan Station 1, show the terrain characteristics with (1) and without snow cover (2).
(bottom)  showing a clear increase of incoming solar radiation while snow season advance, with high variability due to meteorological factors.

Figure 2 :
Figure 2: Daily average temperature, snow depth and net solar radiation (short wave) at the
shows two examples of BRTs, obtained for 2 May 2012 (upper panel) and 3 April 2013 (bottom panel), which accounted for the largest amount of snow accumulation in www.the-cryosphere.net/8

Figure 4 .
Figure 4. Wind roses from the automatic weather station placed at the catchment obtained for a 15-day period.

Figure 5 :Figure 5 .
Figure 5: Willmott's D and r 2 values between the observed and predicted SD, based on the 854

Figure 6 .
Figure 6.Binary regression tree obtained for 2 May 2012 (top) and 3 April 2013 (bottom).The final nodes (with ellipses) show the predicted SD in the zone having the specified terrain characteristics.At each branch point, one topographic variable is considered; if the value is less than the specified value, the left branch is selected, but if it is equal to or greater than the specified value, the right branch is selected.

Figure 7 :
Figure 7: Mean contribution of topographic variables to models and Pearsons'r coefficients

Figure 7 .
Figure 7. Mean contribution of topographic variables to models and Pearsons' r coefficients vs. the coefficients of variation for all considered surveys.Upper panel shows Pearson's r coefficients; middle panel shows beta coefficients of the multiple linear regression; and the bottom panel shows the contribution to the explained variance to the binary regression trees.Each graph point is accompanied with its variable, and in the case of MLR and BRT, in brackets the number of days in which each variable was included in models.

Table 1 .
Summary statistics of the snowpack distribution and the snow covered area of the basin.Note that the snow covered area is expressed as a % of the total area surveyed by the TLS, and the mean SD is the average of all SDs not including zero values.
• for 22 February 2012, 270 • for 2 and 17 April 2012 and 225 • for the three surveys in May 2012; in 2013, 315 • was the best correlated direction for 17 February, and 270 • for the other five surveys of the snow season.

Table 2 .
Pearson's r coefficients between SD and Sx, calculated for the eight studied wind directions over the survey days.Correlations that were statistically significant (α < 0.05) in at least half of the samples (500 out of 1000 samples) from the Monte Carlo approach, and bold r coefficients represent the best correlated Sx direction for a specific survey day. *

Table 3 .
Pearson's r coefficients between SD and the topographic variables.
* Correlations that were statistically significant (α < 0.05) in at least half of the samples (500 out of 1000 samples) from the Monte Carlo approach, and bold r coefficients represent the best correlated topographic variable for a specific survey day.

Table 4 .
Multiple linear regression beta coefficients for each independent variable and sampled day.
amongst survey days.Thus, Sx was the second most influential variable during May (except for 24 May 2012), following the largest snowfall in the season (which occurred the 1 May 2012), and Elevation was the most important variable in the other surveys during 2012.Northing also had an evident influence during the two first surveys of the year, but subsequently had minimal explanatory capacity, as was the case for all the other variables.In 2013 TPI was also the main

Table 5 .
Contribution of the various topographic variables to the explained variance of SD distribution in the binary regression tree models for 2012 and 2013.Values have been rescaled from 0 to 100.