Spatio-temporal variability of snow water equivalent in the extra-tropical Andes Cordillera from distributed energy balance modeling and remotely sensed snow cover

Seasonal snow cover is the primary water source for human use and ecosystems along the extratropical Andes Cordillera. Despite its importance, relatively little research has been devoted to understanding the properties, distribution and variability of this natural resource. This research provides high-resolution (500 m), daily distributed estimates of end-of-winter and spring snow water equivalent over a 152 000 km domain that includes the mountainous reaches of central Chile and Argentina. Remotely sensed fractional snow-covered area and other relevant forcings are combined with extrapolated data from meteorological stations and a simplified physically based energy balance model in order to obtain melt-season melt fluxes that are then aggregated to estimate the end-of-winter (or peak) snow water equivalent (SWE). Peak SWE estimates show an overall coefficient of determination R of 0.68 and RMSE of 274 mm compared to observations at 12 automatic snow water equivalent sensors distributed across the model domain, withR values between 0.32 and 0.88. Regional estimates of peak SWE accumulation show differential patterns strongly modulated by elevation, latitude and position relative to the continental divide. The spatial distribution of peak SWE shows that the 4000– 5000 m a.s.l. elevation band is significant for snow accumulation, despite having a smaller surface area than the 3000– 4000 m a.s.l. band. On average, maximum snow accumulation is observed in early September in the western Andes, and in early October on the eastern side of the continental divide. The results presented here have the potential of informing applications such as seasonal forecast model assessment and improvement, regional climate model validation, as well as evaluation of observational networks and water resource infrastructure development.


Introduction
Accurately predicting the spatial and temporal distribution of snow water equivalent (SWE) in mountain environments remains a significant challenge for the scientific community and water resource practitioners around the world.The Andes Cordillera, a formidable mountain range that constitutes the backbone of the South American continent, remains one of the relatively least studied mountain environments due to its generally low accessibility and complex topography.The extratropical stretch of the Andes, extending south from approximately latitude 27 • S, is a snow-dominated hydrological environment that provides key water resources for a majority of the population in Chile and Argentina.Until now, a very sparse network of snow courses and automated snow measuring stations (snow pillows) has been the only source of information about this key resource.In a context of sustained climate change characterized by warming trends and likely future precipitation reductions (Vera et al., 2006;Vicuña et al., 2011), it becomes ever more relevant to understand the past dynamics of the seasonal snowpack in order to validate predictive models of future snow water re-Published by Copernicus Publications on behalf of the European Geosciences Union.
sources.This research presents the first spatially and temporally explicit high-resolution SWE reconstruction over the snow-dominated extratropical Andes of central Chile and Argentina based on a physical representation of the snowpack energy balance (Kustas et al., 1994) and remotely sensed snow extent (Dietz et al., 2012) between years 2001 and 2014.A key advantage of the presented product is its independence from notoriously scarce and unreliable precipitation measurements at high elevations.Estimates of maximum SWE accumulation and depletion curves are obtained at 500 m resolution, coincident with the MODIS MOD10A1 Fractional Snow Cover product (Hall et al., 2002).
Patterns of hydroclimatic spatio-temporal variability in the extratropical Andes have been studied with increased intensity over the last couple of decades, as pressure for water resources has mounted while at the same time rapid changes in land use and climate have highlighted the societal need for increased understanding of water resource variability and trends under present and future climates.The vast majority of studies have relied on statistical analyses of instrumental records and regional climate models to present synopticscale summaries of precipitation (e.g., Aravena and Luckman, 2009;Falvey and Garreaud, 2007;Garreaud, 2009), temperature (Falvey and Garreaud, 2009), snow accumulation (Masiokas et al., 2006) and streamflow variability (Cortés et al., 2011;Núñez et al., 2013).Currently, no highresolution, large-scale distributed assessments of snow water equivalent are available for the Andes region.
The SWE reconstruction method seeks to estimate endof-winter accumulation by back accumulating melt energy fluxes during the depletion season.The methods and assumptions required for SWE reconstruction have been tested and refined since initial development (Cline et al., 1998).Applications across a variety of scales have been presented in recent years.In the Sierra Nevada, Jepsen et al. (2012) compared SWE reconstructions to distributed snow surveys in a 19.1 km 2 basin (R 2 = 0.79), while Guan et al. (2013) obtained good correlation with SWE observations from an operational snow sensor network across the entire Sierra Nevada (R 2 = 0.74).In the Rocky Mountains, Jepsen et al. (2012) obtained an R 2 value of 0.61 when comparing reconstructed SWE to spatial regression from snow surveys, and Molotch (2009) estimated SWE with a mean absolute error (MAE) of 23 % compared to intensive study areas.A useful discussion on the uncertainties of the SWE reconstruction method -albeit one based on temperature-index melt equations -was presented by Slater et al. (2013), who demonstrated that errors in forcing data are at least, if not more, important than snow-covered area data availability.The vast majority, if not all, of SWE reconstruction exercises have been developed in the northern hemisphere, under environmental conditions quite different from those predominant in the extratropical Andes Cordillera.Here, snow distribution and properties have been analyzed in a few local studies (e.g., Ayala et al., 2014;Cortés et al., 2014b;Gas-coin et al., 2013), but no large-scale estimations at a relevant temporal and spatial resolution for hydrologic applications have been presented.In fact, the Andes of Chile and Argentina display near-ideal conditions for the SWE reconstruction approach due to (1) the near absence of forest cover over a large fraction of the domain where snow accumulation is hydrologically significant; (2) the sharp climatological distinction between wet (winter: June through August) and dry (spring/summer: September through March) seasons, with most of annual precipitation falling during the former; and (3) the low prevalence of cloudy conditions during the spring and summer months over the mountains, which afford a high availability of remotely sensed snow cover information.Conversely, the SWE reconstruction presented here is certainly subject to a series of uncertainty sources, such as the sparseness of the hydrometeorological observational network, which limits both the availability of forcing and validation data.
However, this is the first estimation of peak SWE and snow depletion distribution at this scale and spatial resolution for the extratropical Andes, and the information shown here can be useful for several applications such as understanding yearto-year differential accumulation patterns that may impact the performance of seasonal streamflow forecast models that rely on point-scale data only.Also, the SWE reconstruction can be used to validate output from global or regional climate models and reanalysis, which are being increasingly employed to estimate hydrological states and fluxes in ungauged regions.By analyzing the spatial correlation of snow accumulation and hydrometeorological variables, distributed SWE estimates can inform the design of improved climate observation networks.Likewise, from analyzing the obtained SWE estimates in light of the necessary modeling assumptions and data availability we are able to highlight future research directions aimed at quantifying and reducing these uncertainties.
The objectives of this research include the following: (1) to assess the dominant patterns of spatio-temporal variability in snow water equivalent of the snow-dominated extratropical Andes Cordillera; and (2) to explicitly evaluate the strengths and weaknesses of the SWE reconstruction approach in different sub-regions of the extratropical Andes using snow sensors and distributed snow surveys.

Study area
Figure 1 shows the study area, which includes headwater basins in the Andes Mountains of central Chile and Argentina, between 27 and 38 • S.These basins supply freshwater to low valleys located on both sides of the cordillera, a topographic barrier more than 5 km high that strongly controls the spatial variability in atmospheric processes (Garreaud, 2009;Montgomery et al., 2001).In Chile, runoff from the Andes Mountains benefits 75 % of the popula- tion (http://www.ine.cl) as well as most of the country's agricultural output, hydropower and industrial activities.In the case of Argentina, 7 % of the population is located in the provinces of La Rioja, San Juan, Mendoza and northern Neuquén (http://www.indec.gov.ar/), with primary water uses in agriculture and hydropower.The selected watersheds have unimpeded streamflow observations and a snowdominated hydrologic regime (Fig. 2).River basins included in this study have been grouped in eight clusters, or hydrologic response units, based on the seasonality of river flow, numbered C1 to C8 in Fig. 1b.Due to differences in topography and locations of stream gages, the number of headwater basins contained within clusters differs markedly on both sides of the cordillera, with larger watersheds on the Argentinean side.
The hydro-climate is mostly controlled by orographic effects on precipitation (Falvey and Garreaud, 2007) and interannual variability associated with the Pacific Ocean through the El Niño-Southern Oscillation and Pacific Decadal Oscillation (Masiokas et al., 2006;Newman et al., 2003;Rubio-Álvarez and McPhee, 2010).Precipitation is concentrated in winter months on the western slope (Aceituno, 1988) and sporadic spring and summer storms occur on the mountain front plains of the eastern slope.The vegetation cover presents a steppe-type condition on the western slope up to 33 • S, transitioning to the south into tall bushes and sparse mountain forest.On the eastern slope the steppe vegetation prevails until 37 • S, with an intermittent presence of mountain forests in the Patagonian plains (Eva et al., 2004).
Figure 2 summarizes the dominant climatology and associated hydrological regime of rivers in the study region.The temperature seasonality (upper left panel) is typical of a temperate, Mediterranean climate, and precipitation is strongly concentrated in the fall-winter months of May through August (upper right panel).The hydrological regime is markedly snow-dominated in the northern part of the domain, which can be seen from the sharp increase in river flow from October and into the summer months of Dec, Jan and Feb (lower right panel) that follows the seasonal melt of snow (lower left panel).Only rivers in the southern subregion display a significant rainfall-dominated seasonal hydrograph.The importance of SWE for the region is demonstrated by the fact that for the studied basins, ablation-season (September-March) river flow accounts for two-thirds of av-E.Cornwell et al.: Spatio-temporal variability of snow water equivalent in the extra-tropical Andes Cordillera erage annual streamflow.Maximum SWE accumulation is reached between the months of August and September on the western side and between late September and early October on the eastern side (Fig. S4 in the Supplement).Scattered snow showers in mid-spring (September through November) affect the study area, but they do not affect significantly the decreasing trend of snow-covered area during the melt season (see timing of peak SWE and fractional snow-covered area (fSCA) analysis in online supplementary material).This feature is essential for choosing the SWE reconstruction methodology used in this work, which is most applicable to snow regimes with distinct snow accumulation and snow ablation seasons.
By and large, the existing network of high-elevation meteorological stations does not include appropriately shielded solid precipitation sensors.Some climate reanalysis products exist, but their representation of Andean topography is crude, and their spatial resolution is not readily amenable to hydrological applications without significant bias correction (Krogh et al., 2015;Scheel et al., 2011).Previous attempts at estimating precipitation amounts at high-elevation reaches in the Andes suggest uncertainties on the order of 50 % (Castro et al., 2014;Falvey and Garreaud, 2007;Favier et al., 2009).In some basins, runoff is partially dictated by glacier contributions, which occur in summer.According to the Randolph Glacier Inventory (http://www.glims.org/RGI/), the central Andes Cordillera has a glacier area of 2245 km 2 between 27 and 38 • S, which is equivalent to 1.5 % of the modeling domain surface area (∼ 152 000 km 2 ).

SWE reconstruction model
A retrospective SWE reconstruction model based on the convolution of the fSCA depletion curve and time-variant energy inputs for each domain pixel is implemented.For each year, the model is run at a daily time step between 15 August (end of winter) and 15 January (mid-summer).This time window ensures that the most likely time at which peak SWE occurs is captured -which itself is variable from year to year -and the almost complete depletion of the seasonal snowpack.Isolated pixels with non-negative fSCA values may remain after 15 January at glacier and perennial snowpack sites.However, the relative area that these pixels represent with respect to the entire model domain is very low (< 1.5 %), and can be neglected in the context of this work.
The energy balance model adopted here derives from the formulation proposed by Brubaker et al. (1996), which considers explicit net shortwave and longwave radiation terms and a conceptual, pseudo-physically based formulation for turbulent fluxes that depends only on the degree-day air temperature: where M p is potential melt; Q nsw is the net shortwave energy flux; Q nlw is the net longwave energy flux; T d is the degree-day temperature, a r (mm • C −1 day −1 ) is the restricted degree-day factor, and f B is the energy-to-mass conversion factor with a value of 0.26 (mm W −1 m 2 day −1 ).Actual melt is obtained by multiplying potential melt by fractional snow cover area: where fSCA fc is the fSCA MOD10A1 estimate adjusted to forest cover correction by a vegetation fractional f veg (0 to 1) from the MOD44B product (Hansen et al., 2003): The SWE for each pixel is computed for each year by accumulating the melt fluxes back in time during the melt season, starting from the day on which fSCA reaches a minimum value, and up to a date such that winter fSCA has plateaued, according to the relations where SWE 0 is end-of-winter or initial maximum SWE accumulation, SWE n is a minimum or threshold value.The model was run retrospectively until 15 August, an adequate date before which little melt can be expected for most of the winter seasons within the modeling period in this region (please see Fig. S5).

Fractional snow-covered area and land use data
Spatio-temporal evolution of snow-covered area was estimated using the fSCA product from the Moderate Resolution Imaging Spectroradiometer (MODIS) on-board the Terra satellite (MOD10A1 C5 Level 3).The MOD10A1 product provides daily fSCA estimates at 500 m resolution.Percentages of snow extent (i.e., 0 to 100 %) are derived from an empirical linearization of the Normalized Difference Snow Index (NDSI), considering the total MODIS reflectance in the visible range (0.545-0.565 µm; band 4) and shortwave infrared (1.628-1.652µm; band 6) (Hall et al., 2002;Hall and Riggs, 2007).Binary and fractional MODIS fSCA estimates are limited by the use of an empirical NDSI-based method.These errors are notoriously sensitive to surface features such as fractional vegetation and surface temperature (Rittger et al., 2013).Arsenault et al. (2014) reviewed MODIS fSCA accuracy estimates from several studies under different climatic conditions, and report a range between 1.5 and 33 % in terms of absolute error with respect to ground observations and operational snow cover data sets.Errors stem mainly from cloud masking and detection of very thin snow (<10 mm depth), forest cover and terrain complexity.In general, commission and omission errors are greatest in the early and late portions of the snow cover season (Hall and Riggs, 2007) and decrease with increasing elevation (Arsenault et al., 2014).Molotch and Margulis (2008) compared MODIS and Landsat Enhanced Thematic Mapper performance in the context of SWE reconstruction, showing that significant differences in SWE estimates were a result of SCA estimation accuracy and less so of model spatial resolution.The latter conclusion supports the feasibility of using the snow-covered area products at a 500 m spatial resolution for regional-scale studies.In order to minimize the effect of cloud cover on the temporal continuity and extent of the fSCA estimates, the MOD10A1 fSCA product was post-processed by a modified algorithm for non-binary products, based on the algorithm proposed by Gafurov and Bárdossy (2009).Their method is adapted here to the fractional snow cover product, applying a three-step correction consisting of: (1) a pixel-specific linear temporal interpolation over 1, 2 or 3 days prior and posterior to a cloudy pixel; (2) a spatial interpolation over the eight-pixel kernel surrounding the cloudy pixel, retaining information from lower-elevation pixels only; and (3) assigning the 2001-2014 fSCA pixel specific average when steps (1) and ( 2) where not feasible.This step minimized the effect of cloud cover on data availability over the spatial domain, yielding cloud cover percentages ranging from 21 % in September to 8 % in December.
The Normalized Difference Vegetation Index (NDVI) (Huete et al., 2002) derived from the MOD13Q1 v5 MODIS Level 3 product (16 days -250 m) is used to classify forest presence for each model pixel.For pixels classified as forested, both fSCA and energy fluxes where corrected: fractional SCA was modified on the basis of percentage forest cover (Molotch, 2009;Rittger et al., 2013), using the average of the forest percentage product from MOD44B V51.Forest attenuation (below canopy) of energy fluxes at the snow surface was estimated from forest cover following the method from Ahl et al. (2006) assuming invariant NDVI over each melt season.The selected NDVI pattern is obtained by averaging the four NDVI scenes available in the December-January time window through 14 study years.This time window displays the average state of evergreen forest with the maximum amount of data.

Model forcings
Spatially distributed forcings are required at each grid element in order to run the SWE reconstruction model.In order to ensure the tractability of the extrapolation process, we divided the model domain into sub-regions or clusters, composed of one or more river basins.The river basins were grouped using a clusterization algorithm (please see Sect.S2 in the Supplement) based on melt-season river flow volume as described in Rubio-Álvarez and McPhee (2010).Then, spatially distributed variables (surface temperature, fSCA, global irradiance) are combined with homogeneous variables for each cluster (e.g., cloud cover index) and point data from meteorological stations in order to obtain a distributed product as described below.A further benefit of the clustering process is that it allows us to analyze distinct regional features of the SWE reconstruction parameters, input variables and output estimates.
Net shortwave radiation, Q nsw is estimated as a function of incoming solar radiation based on the equation where α s is snow surface broadband albedo; G ↓ is incoming solar radiation (global irradiance); and τ a is the shortwave transmissivity as a function of LAI for mixed forest cover (Pontailler et al., 2003;Sicart et al., 2004), which in turn is estimated as with κ = 0.52 for mixed forest species (DeWalle and Rango, 2008).Equation ( 7) is valid for NDVI values between 0.16 and 0.87.Global irradiance under cloudy sky conditions is estimated considering a daily distributed spatial pattern of clear sky irradiance G c↓ derived by the r.sun GRASS GIS module (Hofierka and Suri, 2002;Neteler et al., 2012) and the clear sky index K c derived from the insolation incident on a horizontal surface from the "Climatology Resource for Agroclimatology" project in the NASA Prediction Worldwide Energy Resource "POWER" (http://power.larc.nasa.gov/) 1 • × 1 • gridded product.
In Eq. ( 8), G r↓ and G c↓ are spatial averages over each hydrologic response unit (cluster) of the POWER and r.sun-derived products, respectively.A snow-age decay function based on snowfall detection is implemented to estimate daily snow surface albedo (Molotch and Bales, 2006) constrained between values of 0.85 and 0.40 (Army Corps of Engineers, 1960).Snowfall events were diagnosed using a unique minimum threshold for fSCA increments of 2.5 % for each hydrologic unit area.
Net longwave radiation estimates are derived using where T a is air temperature, T s is the snow surface temperature, ε s is the snow emissivity (i.e., 0.97), ε sf is the canopy emissivity (i.e., 0.97), f sv is the sky-view factor (i.e., assumed equal to shortwave transmissivity; Pomeroy et al., 2009;Sicart et al., 2004), σ is the Stefan-Boltzmann constant, and L ↓ is the incoming longwave radiation.Air vapor pressure (e a ) required for longwave radiation estimates was derived from air temperature and relative humidity, which in turn was assumed constant throughout the melt period and equal to 40 % based on observations at selected high-elevation meteorological stations.The multiplying factor (1 + a c C 2 ) represents an increase in energy input relative to clear sky conditions due to cloud cover, where a c equals 0.17 and C = 1 − K c is an estimate of the cloud cover fraction (DeWalle and Rango, 2008).Spatially distributed air temperature is generated by combining daily air temperature recorded at index meteorological stations and a weekly spatial pattern of skin temperature derived from the MODIS Land Surface Temperature product (MOD11A1.V5) (Wan et al., 2002(Wan et al., , 2004)).The product MOD11A1 V5 Level 3 estimates surface temperature from thermal infrared brightness temperatures under clear sky conditions using daytime and nighttime scenes and has been shown to adequately represent measurements at meteorological stations (R 2 ≥ 0.7), displaying moderate overestimation in spring and underestimation in fall (Neteler, 2010).Other studies have reported similar accuracies, with RSME values around 4.5 • K in cold mountain environments (Williamson et al., 2014).Taking into account the high correlation between air temperature and LST (Benali et al., 2012;Colombi et al., 2007;Williamson et al., 2014), we define where T a base is daily air temperature at an index station for each cluster and T a is the difference in air temperature between any pixel and the pixel where the index station is located.To determine T a we use a linear regression between MODIS LST data and T a considering pairs of stations located at high-altitude and valley (base) sites, taking into account the melt season average values over the 2001-2014 period.In Eq. ( 11), LST − LST base denotes the difference between skin temperatures from any pixel and the index station pixel.The linear regression between skin temperature and air temperature differences has a slope µ of 0.65, an intersect ν of −0.5 and R 2 of 0.93 (Fig. S3).Estimation of LST during cloudy conditions is done as follows: (1) a pixel-specific linear temporal interpolation is performed over 1 and 2 days prior and posterior to the cloudy pixel; and (2) estimation of remaining null values by an LST-elevation linear regression (Rhee and Im, 2014).
This spatial extrapolation method was preferred over more traditional methods -for example, based on vertical lapse rates (Minder et al., 2010;Molotch and Margulis, 2008)after initial tests showed that the combined effect of the relatively low elevation of index stations and the large vertical range of the study domain resulted in unreasonably low air temperatures at pixels with the highest elevations.Likewise, the scarcity of high-elevation meteorological stations and the large spatial extent of the model domain precluded us from adopting more sophisticated temperature estimation methods (e.g., Ragettli et al., 2014).
Snow surface temperature and degree-day temperature are estimated (Brubaker et al., 1996) as where T is the difference between air and snow surface temperature.To the best of our knowledge, no direct, systematic values of snow surface temperature exist in this region, so for the purposes of this paper we adopt an average value T = 2.5 [ • C], following the suggestion in Brubaker et al. (1996).Slightly higher values ranging from 3 to 6 • C are shown for continental and alpine snow types (Raleigh et al., 2013), indicating an additional source of uncertainty over net longwave radiation computations.More sophisticated parametrizations for T s , for example based on heat flow through the snowpack, have been proposed (e.g., Rankinen et al., 2004;Tarboton and Luce, 1996) but those require explicit knowledge about the snowpack temperature profile and/or more complex model formulations to estimate the internal snowpack heat and mass budgets simultaneously.The a r coefficient in the restricted degree-day energy balance equation was computed using a combination of station and reanalysis data, and assumed spatially homogeneous within each of the clusters that subdivide the model domain.Brubaker et al. (1996) propose a scheme in which this parameter can be explicitly computed from air and snow surface temperature, air relative humidity, and atmospheric pressure and wind speed.Wind speed was obtained from the NASA POWER reanalysis described previously.A correction for atmospheric stability is applied on the bulk transfer coefficient C h according to the formulation presented by Kustas et al. (1994), assuming a surface roughness of 0.0005 m: where Ri is the Richardson number, g is the gravity acceleration (9.8 [m s −2 ]), z is the standard air temperature measurement height (2 m) and u is wind speed.The calculation of Ri and a r is based on the standard assumptions of T s at the freezing point and a water vapor saturated snow surface over all high-elevation meteorological stations with available air temperature and relative humidity records (Molotch and Margulis, 2008).Further in the text, we discuss some implications of these assumptions and of the input data used on the ability of the model of simulating relevant components of the snowpack energy exchange.
Table 1 shows the main cluster characteristics and regionalized model parameters.It can be seen that for those clusters located in the southern and middle reaches of the model domain, the a r parameter values range from 0.10 to 0.23 (cm • C −1 day −1 ), which is similar to values reported in previous studies performed in other mountain ranges in the Northern Hemisphere (0.20-0.25 in Martinec, 1989; 0.17 in Kustas et al., 1994;0.20 in Brubaker at al., 1996;0.15 in Molotch and Margulis, 2008).However, values associated with the northernmost clusters of our study area are quite low, reaching under 0.02 for the C1 cluster in northern Chile.Clear sky index (K c ) values range between 0.78 and 0.89, which is similar to values reported by Salazar and Raichijk (2014), who estimate K c values on the order of 0.90 for a single location at 1200 m a.s.l. in northern Argentina.A 5 to 6 • C difference can be observed in mean air temperature at index stations between the northern and southern edge of the domain.Temperatures for the C4 cluster are subject to greater uncertainty, because no high-elevation climate station data were available for this study (Fig. S4).Forest cover values are lower than 6 % throughout the model domain, with the exception of cluster C3, with a value of 13.8 %.The difference in forest cover between clusters C3 and C8 can be attributed to the precipitation shadow effect induced by the Andes ridge.Forest corrections applied to MODIS fSCA resulted in a 17 % increase with respect to the original values over the southern sub-domain (C3).

Evaluation data: SWE, snow depth and river flow observations
Operational daily snow-pillow data from stations maintained by government agencies in Chile and Argentina were available for this study (Table 2).Only stations with 10 or more years of record were included, and manual snow course data were neglected because of their discontinuous nature.Approximately 10 % of observed maximum SWE accumulation values were discarded due to obvious measurement er- rors and data gaps.An analysis of the seasonal variability of snow-pillow records on the western and eastern slopes of the Andes suggests that the peak-SWE date is somewhat delayed on the latter, by approximately 1 month.Therefore, peak-SWE estimates for Chilean and Argentinean stations are evaluated on 1 September and 1 October, respectively, although in the results section we show values for 15 September in order to use a unique date for the entire domain.Manual snow depth observations were taken in the vicinity of selected snow-pillow locations in order to evaluate the representativeness of these measurements at the MODIS grid scale during the peak-SWE time window.These depth observations were obtained in regular grid patterns within an area the approximate size of a MODIS pixel (500 m), centered about the snow-pillow location.On average, 120 depth observations spaced at approximately 50 m increments were obtained at each snow-pillow site.Snow density was estimated by a depth-weighted average of snow densities measured in snow pits with a 1000-cc snow cutter.Samples where obtained either at regular 10 cm depth intervals along the snow pit face, or at the approximate mid depth of identifiable snow strata for very shallow snow pack conditions.Weights were computed as the fraction of total depth represented by each snow sample.
Distributed snow depth observations were available from snow surveys carried out during late winter between 2010 and 2014 at seven study catchments on the western side of the Andes, between latitudes 30 and 37 • S (Fig. 1, Table 3).Snow depths were recorded with 3 m graduated avalanche probes inserted vertically into the snow pack.Depending on the terrain conditions, between three and five individual point snow depth measurements were obtained at each location, from which a mean snow depth and standard error are calculated; i.e., three-point observations are made forming a line with a spacing of 1 m and five-point observations are made forming a cross with an angle of 90 • and a spacing of 1 m.Pixel-scale SWE estimates are obtained by averaging all depth observations within the limits of MODIS pixels and multiplying them by density observations from snow pits excavated at the time of each snow survey, i.e., two or three snow pits per field campaign.After this, individual depth observations are converted into SWE for model validation.Modeled SWE values are averaged at all MODIS pixels where manual depth observations are available, and their summary statistics are compared to those of SWE estimated from manual depth observations at the same pixels, multiplied by average density from snow pits.
Spring and summer season (September to March) total river flow volume (SSRV) for the 2001-2014 period is obtained from unimpaired (no human extractions) streamflow records at river gauges located in the mountain front along the model domain.Data were pre-selected leaving out series that showed too many missing values, and verified through the double mass curve method (Searcy and Hardison, 1960) in order to discard anomalous values and to ensure homogeneity throughout the period of study.Regional consistency was verified through regression analysis, only including streamflow records with R 2 values greater than 0.5 among neighboring catchments.Missing values constituted about 3.7 % of the entire period and were filled through linear regression.

Model validation
Figure 3 compares reconstructed peak SWE (gray circles) to observed values at three snow-pillow locations (black diamonds) where additional validation sampling at the MODIS pixel scale was conducted (box plots).At the Cerro Vega Negra site (CVN), located in cluster C1, the model overestimates peak SWE (1 September) with respect to the snowpillow value by 97 % in 2013 and by 198 % in 2014.At the Portillo site (POR, cluster C2), reconstructed SWE underestimates recorded values by 51 % in 2013 and 72 % in 2014.At the Laguna Negra site (LAG, also C2), reconstructed peak SWE slightly overestimates recorded values (8 %) (Table 4).However, reconstructed SWE compares favorably to distributed manual SWE observations obtained in the vicinity of the snow pillows at the POR and LAG sites.At POR, model estimates approach upper (2012) and lower (2013 and 2104) quartiles, while at LAG the model estimates are closer to the minimum value observed in 2013 and very similar to the observed mean in 2014.
Figure 4 depicts the comparison between reconstructed SWE and snow surveys carried out at pilot basins throughout the model domain.From left to right, it can be seen that the model slightly overestimates SWE with respect to observations at CVN (i.e., 18 % overestimation).Further south, there is a very good agreement at ODA-MAR (i.e., 4 % underestimation), with less favorable results at MOR-LVD (i.e., 39 % underestimation) and OB-RBL (i.e., 36 % underestimation).At CHI the model significantly underestimates SWE (i.e., by 67 %); note that this site is heavily forested.For the 2013a and 2014a boxes (Fig. 4) -which correspond to clearing sites -there is still underestimation, but of lesser magnitude (20 %).Summarizing, we detect model overestimation with respect to snow survey medians in four cases and underestimation in fifteen cases.In 11 out of 19 cases, reconstructed SWE lies within the snow survey data uncertainty bounds (standard deviation).Figure 5 shows a comparison between model estimates of peak (15 September) SWE and corresponding observations at snow-pillow sites.In general, directly contrasting pixel-based estimates with sensor observations should be attempted with caution.In areas with complex topography, slight variations in the position of the sensor with respect to the model grid, combined with high spatial variability in snow accumulation could lead to large differences between model estimates and observations.Also, small-scale variations in snow accumulation near the sensor, for example induced by protective fences, could introduce bias to the results (e.g., Meromy et al., 2013;Molotch and Bales, 2006;Rice and Bales, 2010).Taking the above into consideration, Fig. 5 suggests that the model tends to overestimate observed peak SWE at the two northernmost sites on the Chilean side (QUE and CVN); the equivalent cluster on the Argentinean side (C4) lacks SWE observations.The R 2 values indicated below refer to the best linear fit; regression line slope and intercept coefficients are provided in Table 4. Overall, we find a better agreement at the eastern slope sites (i.e., R 2 = 0.74) than at their western counterparts (i.e., R 2 = 0.43), with a combined R 2 value of 0.61.Individually, the worst and best linear agreements are obtained at POR (R 2 = 0.32) and LOA (R 2 = 0.88), respectively.Time series of observed SWE and model estimates for these two extreme cases are shown in the supplementary online material, and indicate a significant degree of inter-annual variability in model discrepancies in terms of peak SWE, but less in terms of, for instance, snow cover duration.Average standard error, SE x is 284 mm (SE x = 242 mm at the western slope; SE x = 302 mm at the eastern slope), with a range between 72 mm (TOS) and 378 mm (ATU) (Table 4).Relative errors display some variability, with overestimation higher than 30 % at the two northernmost (QUE and CVN) and at the southernmost (PEH) snow pillows.For all other snow  pillows, the model estimates are lower than the sensor observation; the range of relative errors for those sites with underestimation goes from −52 to −5 %.

Correlation with melt-season river flows
Under the assumption of unimpaired flows, peak SWE and seasonal flow volume should show some degree of correlation, even though no assumptions can be made here about other relevant hydrologic processes, such as flow contributions from glaciated areas, subsurface storage carryover at the basin scale and influence of spring and summer precipitation.Differences can be expected due to losses to evapotranspiration and sublimation affecting the snowpack and soil water throughout the melt season.Hence, basin-averaged peak SWE should always be higher than melt season river volume.A clear regional pattern emerges when inspecting the results of this comparison in Fig. 6.Correlation between peak SWE and melt season river flow is higher in clusters C1 and C4 with R 2 values of 0.84 and 0.86, respectively.The re- sult for Cluster C4 indicates that liquid precipitation during the melt season (Fig. 2) does not result in decreased correlation between peak SWE and river flow.Clusters C2, C5, C6 and C7 display a somewhat lower correlation, with some individual years departing more significantly from the overall linear trend.R 2 values range between 0.46 and 0.78 in these cases.Finally, not only are correlation coefficients lower for the southern clusters C3 (R 2 = 0.56) and C8 (R 2 = 0.48), but also estimated peak SWE is always lower than river flow, which indicates the importance of spring and summer precipitation in determining streamflow variability.In fact, Castro et al. (2014) analyze patterns of daily precipitation in this area and document average spring and summer rainfall amounts of approximately 520 mm in C3 and 85 mm in C8.A promising avenue for further research in this region emerges when comparing the correlation between meltseason river flow and the spatially distributed reconstructed product versus that of river flows and snow-pillow data.Table 5 shows values of R 2 for the linear regression between these variables.It can be seen that for two of the three clusters on the western side of the continental divide, the end-ofwinter distributed reconstruction has more predictive power than observed SWE.Only for central Chile the Laguna Negra (LAG, with a value of 0.82) site has a better correlation with river flows, but the reconstructed product has a value of 0.78, which lies in between those found for LAG and for Portillo (POR, with a value of 0.68).For the eastern side of the continental divide, the distributed product shows similar skill than that of snow pillows except for Atuel, which has a very high correlation (R 2 of 0.87) with cluster C6 river flows, and for cluster C7, in which the reconstruction shows higher predictive power (R 2 of 0.89) than the available SWE observations (VAL and PEH).

Regional SWE estimates
Figure 7 shows the 15 September SWE average over the 2001-2014 period obtained from the reconstruction model, and the percent annual deviations (anomalies) from that average.Steep elevation gradients can be inferred from the climatology, as well as the latitudinal variation expected from precipitation spatial patterns.For the northern clusters (C1 and C4), the peak SWE averaged over snow-covered areas is on the order of 300 mm, while in the middle of the domain (C2, C5, C6), it averages approximately 750 mm.The southern clusters (C3, C7, C8) do show high accumulation averages (≈ 650 mm), despite the sharp decrease in the Andes elevation south of latitude 34 • S. The anomaly maps convey the important degree of inter-annual variability, as well as distinct spatial patterns associated with it.Between 2001Between and 2014Between , years 2002Between and 2005 stand out for displaying large positive anomalies throughout the entire mountainous region of the model domain, with values 2000 mm and more above the simulation period average.Other years prior to 2010 show differential accumulation patterns, where either the northern or southern parts of the domain are more strongly affected by positive or negative anomalies.Overall, the northern clusters (C1 and C4) show above-average accumulation in only 3 (2002, 2005 and 2007) of the 14 simulated years, whereas the other clusters show above-average accumulation for 6 years (2001, 2002, 2005, 2006, 2008 and 2009).In particular, years 2007 and 2009 show a bimodal spatial structure, with excess accumulation (deficit) in the northern (southern) clusters during the former, and the inverse pattern in the latter year.
A longitudinal pattern in the distribution of negative anomalies can be discerned from Fig. 7, whereby drought conditions tend to be more acute on one side of the divide versus the other.Conversely, during positive anomaly years, both sides of the Andes seem to show similar behavior.Further research on the mechanisms of moisture transport during below-average precipitation years may shed light on this result.
Figure 8 provides a different perspective on the region's peak SWE climatology by presenting our results aggregated into elevation bands for each hydrologic unit.Elevation bands are defined at 1000 m increments starting from 1000 m a.s.l.Crosses indicate average peak SWE for each band (mm), and circle areas are proportional to the surface area covered by each elevation band.From north to south, hydrologic unit C4 shows slightly higher SWE than C1 between 3000 and 5000 m a.s.l., but much larger surface areas (∼ 32 000 vs. ∼ 17 000 km 2 ), indicating a larger water resource potential.C2 stands out as having the greatest areaweighted cluster SWE and the greatest SWE for each elevation band.Compared to its counterpart on the eastern side of the Andes range (C5), C2 shows higher accumulations (up to ∼ 1800 mm) at all elevations.The area included between 2000 and 4000 m a.s.l.(∼ 13 000 km 2 ), which shows an estimated peak SWE accumulation on the order of 600 mm, represents the most predominant snow volume accumulation zone.Although the 4000-5000 m a.s.l.elevation band contributes approximately half the 2000-4000 band surface area in C2, its average peak SWE is roughly twice that of the 3000-4000 band (∼ 6000 km 2 ).This makes this subregion interesting for future research, because most snow observations in the area are obtained below 4000 m a.s.l.; the same is true for unit C5.Further to the south, the barrier effect of the Andes is also suggested by the displacement of the SWE- elevation distribution in C6 and C7 when compared to C3.
On the eastern side of the model domain, it is interesting to see a steepening of the average peak SWE elevation profile between C6 and C8, suggesting that C8 is less affected by Andes blockage than its northern counterparts.Estimated net energy inputs (Fig. 9) shows a decrease from the northern (C1 and C4) into the mid-range clusters (C2, C5 and C6), with increases again in the southern reaches of the domain (C3, C7 and C8).This is a result of a combination of an increasing trend in net shortwave radiation in the southnorth direction and a reverse spatial trend in net longwave radiation exchange, which increases (approaches less negative values) in the north-south direction.Modeled turbulent energy fluxes (Eq. 1) are negligible in the northern clusters, but their contribution to the net energy exchange increases with latitude as a result of the spatial variation in the a r parameter.
Figure 10 shows the temporal (seasonal) variation in average fSCA and SWE for each cluster, and Table 6 shows peak SWE at the watershed scale, averaged both over the entire basin and over the snow-covered area.Maximum fSCA increases in the north-south direction, consistent with the climatological increase in winter precipitation and decrease in temperature.A dramatic increase in snow coverage is observed between the northern (i.e., C1 and C4) and adjacent southern clusters (i.e., C2 and C5), with average peak fSCA increasing from 20 to 50 %.The highest average snow coverage is observed for cluster C8, with more than 60 %.Snow water equivalent displays a similar regional variability with lower seasonal variability than snow cover for all clusters except for C2, where fSCA and SWE variability throughout the melt season are identical.Mean peak SWE in northern Chile is the lowest among the eight clusters, with approximately 100 mm SWE over the 2001-2014 period.The largest estimate is for cluster C2, central Chile, where mean peak SWE exceeds 500 mm.The rain shadow effect of the Andes range is apparent in the comparison of SWE and fSCA in C2 and C5-C6-C7.Fractional snow-covered area is lower on the eastern side because of the larger basin sizes, which increases the proportional area of lower elevation terrain.In addition, peak SWE is approximately 25 % lower on the eastern side, with less than 400 mm SWE for the eastern clusters.Cluster C4 is not affected by this phenomenon, showing higher snow coverage and water equivalent accumulation than its counterpart, C1.Cluster C8 represents an interesting exception in that its average fSCA is the largest within the model domain, but peak SWE is not significantly higher than the estimates in the other clusters on the Argentinean side of the Andes.

Sensitivity analysis
The Andes Cordillera, on the one hand, displays ideal conditions for SWE reconstruction, including low cloud cover, infrequent snowfall during spring and summer, and very low forest cover.On the other hand, the scarcity of basic climate data poses challenges that would affect any modeling exercise.A local sensitivity analysis is implemented in order to gain insights regarding the influence of some of the assumptions required for SWE modeling (Fig. 11).The influence of the clear sky factor (K c ), snow surface albedo (∝ s ), the slope of the LST vs. T a relationship (µ), the a r parameter, and the difference between air and snow surface temperature are explored.Results are shown for the model pixels corresponding to two of the snow-pillow sites, each located at the northern and southern sub-regions of the model domain respectively.The clear sky factor, snow albedo and LST vs. T a slope are the most sensitive parameters at the northern (CVN) site.Increasing the slope in the LST vs. T a relationship results in decreasing temperature at pixels with higher elevations than the index station, thus lowering longwave cooling and resulting in higher SWE estimates.The impact of increasing slope values decreases progressively, because an increasing slope results in increased pixel air temperature, but snow surface temperature cannot exceed 0 • C. The influence of snow albedo is analyzed by perturbing the entire albedo time series for each season from the values predicated by the USACE model.Increasing albedo values restricts the energy available for melt therefore decreasing peak SWE estimates.Again, a nonlinear effect is observed, constrained by a minimum albedo value of 0.4.The sensitivity of the clear sky factor, on the other hand, is monotonic, with increasing values generating more available solar energy, resulting in higher SWE estimates.At the southern site (ALT), the shape of the sensitivity functions is the same as at CVN, but the magnitude of SWE variations as a function of parameter perturbations is smaller.This is likely related to the fact that turbulent fluxes constitute a larger fraction of the simulated overall energy balance at the southern sites; a r parameter values are greater in the southern portions of the domain.There-   fore, perturbations of the other terms account for a smaller fraction of the energy exchange at the southern sites.

Model performance and conceptual energy balance representation
Among the many factors that influence model performance, the sub-region delineation involves the selection of index meteorological stations for extrapolating input data at the domain level.Thus, for example, two adjacent pixels that are part of different sub-regions may be assigned input data derived from two different meteorological stations that are many kilometers apart.It would be preferable to use distributed inputs only, but these were not available for this do-main.Future research is needed to explore alternative strategies for domain clustering.
Overall, the model performance, evaluated against SWE observations, is comparable to that achieved in other mountain regions of the world.Our average coefficient of determination R 2 of 0.68 is lower than that obtained by Guan et al. (2013) in the Sierra Nevada (0.74) when comparing operational snow-pillow observations, although this value is affected by three stations with much lower agreement (POR, LAG, ATU); the median R 2 in our study, on the other hand, is 0.73, which we consider satisfactory in light of the scarcity of forcing data and direct snow properties observations available in this region.The overall relative error is −2 % for observations from snow pillows within our study region, but this value is strongly affected by two stations where we observed significant overestimation (QUE and CVN).When including the remaining ten snow pillows only, relative error increases to −16 %.Given that forest cover is minimal in our modeling domain, we can attribute this bias to either weaknesses in the simplified energy balance model formulation or to errors in the MOD10A1 fSCA product.Previous work in the northern hemisphere (Rittger et al., 2013) has shown that MODIS can underestimate fractional snow cover during the snowmelt season.On the one hand, land cover heterogeneity at spatial resolutions lower than the MODIS scale (i.e., 500 m) results in mixed-pixel detection problems.On the other hand, spectral unmixing based on the NDSI approach tends to underestimate fSCA under patchy snow distributions.In addition, surface temperatures greater than 10 • C -more likely to exist during late spring -induce MODIS fSCA underestimation.Molotch and Margulis (2008) tested the SWE reconstruction model using Landsat ETM and MOD10A1 and found that maximum basin-wide mean SWE estimates were significantly lower when using MOD10A1.More recently, Cortés et al. (2014a) showed that a similar pattern can be seen for the extratropical Andes, whereby MODIS fSCA consistently underestimated LANDSAT TM fSCA retrievals.MODIS fSCA underestimation during spring combined with increased net energy fluxes over the snowpack can result in a marked underestimation (∼ 20 %) for available energy flux for snowpack melting and consequently (∼ 45 %) for maximum SWE (Molotch and Margulis, 2008).Comparisons against spatial interpolations from intensivestudy areas in the Sierra Nevada or Rocky Mountains (e.g., Erxleben et al., 2002;Jepsen et al., 2012) are not directly applicable, because in this study we do not employ interpolation methods to derive our manual snow survey SWE estimates.However, the average overestimation found with respect to snow survey data could be explained by the fact that manual surveys are limited by site accessibility and sampling procedures.For example, snow probes utilized are only 3.0 m long, which precludes observation of deeper snowpack; likewise, deep snow is expected in sites exposed to avalanching, which were generally avoided in snow survey design due to safety considerations.On the other hand, manual snow surveys do not visit steep snow-free areas where snow depth is expected to be lower than the 500 m pixel reconstruction.The combined effect of these two contrasting effects is the subject of further research in this region.
Another possible explanation for model errors is the simplified formulation of the energy balance equation, which may be problematic when applied over a large, climatically variable model domain.To explore the implications of the simplified energy balance with respect to model errors, we focus on the representation of turbulent energy fluxes, represented here through a linear temperature-dependent term.Figure 12 describes the spatial distribution of the a r parameter, and its dependence on air temperature and relative humidity observed at index meteorological stations.The implication for energy balance modeling is that turbulent fluxes would account for a very small portion of the snowpack energy and mass balance in the northern area (C1 and C4), which is characterized by low air temperatures and relative humidity, which yield very low a r values.The reader must recall that a r values were computed based on index station data and assumed spatially homogeneous over each cluster.The simplified model formulation used in this research, however, although pseudo-physically based -compared to degree-day or fully calibrated models -allows only for positive net turbulent fluxes, because both the a r and the degree-day temperature index are positive values.However, previous studies in this region (Corripio and Purves, 2005;Favier et al., 2009) have suggested that latent heat fluxes have a relevant role because of high sublimation rates favored by high winds and low relative humidity conditions predominant in the area.
In order to diagnose differential performance of the model across the hydrologic units defined in this study, we compute the Bowen ratio (β) at the point scale from data available only at the few high-elevation meteorological stations in the region with recorded relative humidity.The calculations show that at stations located within cluster C1, latent heat fluxes are opposite in sign and larger in magnitude than sensible heat fluxes (Fig. S6).While this results in net turbulent cooling of the snowpack, this energy loss is not considered in our simplified energy balance approach.Note that for the clusters C5, C6, C7 and C8, all located on the eastern (Argentinean) slope of the Andes, sensible and latent heat fluxes are positive, compared to negative latent heat fluxes for all the index stations within clusters C2 and C3 on the Chilean side.This result is consistent with Insel et al. (2010), who applied a regional circulation model (RegCM3) in the area and showed a significant difference in relative humidity (∼ 70 % on the eastern side vs. ∼ 40 % on the western side).The fact that we extrapolate the a r parameter value based on relatively low elevation meteorological observations throughout the southern Argentinean hydrologic units may result in a yet not quantified overestimation of seasonal energy inputs and peak SWE for those clusters.

Conclusions
Snow water equivalent is the foremost water source for the extratropical Andes region in South America.This paper presents the first high-resolution distributed assessment of this critical resource, combining instrumental records with remotely sensed snow-covered area and a physically based snow energy balance model.Overall errors in estimated peak SWE, when compared with operational station data, amount to −2.2 %, and correlation with observed melt-season river flows is high, with an R 2 value of 0.80.MODIS fractional SCA data proved adequate for the goals of this study, affording high temporal resolution observations and an appropriate spatial resolution given the extent of the study region.These results have implications for evaluating seasonal water supply forecasts, analyzing synoptic-scale drivers of snow accumulation, and validating precipitation estimates from regional climate models.In addition, the strong correlation between peak SWE and seasonal river flow indicates that our results could be useful for the evaluation of alternative water resource projects as part of development and climate change adaptation initiatives.Finally, the regional SWE and anomaly estimates illustrate the dramatic spatial and temporal variability of water resources in the extratropical Andes, and provide a striking visual assessment of the progression of the drought that has affected the region since 2009.These results should motivate further research looking into the climatic drivers of this spatially distributed phenomenon.
The Supplement related to this article is available online at doi:10.5194/hess-20-411-2016-supplement.

Figure 1 .
Figure 1.Study area and model domain: (a) river basins, stream gages (red circles) and sites where snow survey data are available (green circles); (b) hydrologic units (C1 to C8) and snow-pillow stations (white circles).

Figure 2 .
Figure 2. Summarized hydro-climatology of the model domain.Data from meteorological stations located within zones C1, C4, C3 and C8 summarized the hydro-climatological regime of the northwestern, northeastern, southwestern and southeastern zones, respectively.Total SWE is SWE measured at selected snow-pillow stations.

Figure 3 .
Figure 3. Reconstructed SWE validation at selected snow-pillow sites.Black diamonds are instrumental records, gray circles are model estimates, and box plots summarize the manual verification data set around the pillow site.Upper and lower box limits are the 75 and 25 % quartiles, the horizontal line is the median, the white box is the mean, upper and lower dashes represent plus and minus 2.5 SD from the mean, and crosses are outlying values.

Figure 4 .Figure 5 .
Figure 4. Reconstructed SWE validation at pixels with snow survey data.Box plots summarize all individual measurements at pixels colocated with SWE reconstruction.Symbology analogous to Fig. 3.

Figure 6 .
Figure 6.Area-specific spring-summer runoff volume (SSRV) versus peak SWE.Clusters 1 through 3 include rivers on the Chilean (western) slope of the Andes range; clusters 4 through 8 correspond to Argentinean (eastern) rivers.Solid line represents the 1 : 1 line.C4 and C8 SSRV were estimated by the area-transpose method.

Figure 8 .
Figure 8. Maximum SWE through 1000 m elevation bands (EB).Crosses are mean values within EB, lines are the estimated SWE-elevation profile.Circle radius indicates EB area (km 2 ) scaled by 0.05 and takes values from the SWE axis.

Figure 9 .
Figure 9.Time series of energy fluxes over snow surface (average over 14 years) and global average per cluster.Unique axes scale for all plots.

Figure 10 .
Figure 10.Average seasonal evolution of fSCA and SWE in the study region.Lower right panel shows the spatial correlation between time-averaged fSCA, SWE and specific melt-season river discharge.

Figure 11 .
Figure 11.Sensitivity of peak SWE estimates to model forcings and parameters.Average over the 2001-2014 period at selected snowpillow sites.x represents the percentage change over each parameter studied with respect to the base case.

Figure 12 .
Figure 12.Restricted degree-day factor as a function of space (basin cluster) and climatological properties.Bowen (β) coefficient shown between parentheses in the legend.

Table 1 .
Study area subdivision, relevant characteristics and model parameters.

Table 2 .
Snow-pillow measurements available within the study domain.

Table 3 .
Summary of snow depth and density intensive study campaigns.
a Without forest cover (upper part of basin).b With forest cover (lower part of basin).

Table 4 .
Model validation statistics against intensive study area observations around snow pillows and at catchment scale.

Table 5 .
Coefficient of determination R 2 between river melt season flows (SSRV) and estimated and observed SWE (end-of-winter).2014 flows in Argentina unavailable to us at the time of writing.* * 2009 is considered an outlier year for the reconstruction at Argentinean sites. *

Table 6 .
Peak SWE 2001-2014climatology for river basins within the study region.Basin-wide averages, SCA-wide averages and basinwide water volumes shown.