Interannual variability in glacier contribution to runoff from a high‐elevation Andean catchment: understanding the role of debris cover in glacier hydrology

We present a field‐data rich modelling analysis to reconstruct the climatic forcing, glacier response, and runoff generation from a high‐elevation catchment in central Chile over the period 2000–2015 to provide insights into the differing contributions of debris‐covered and debris‐free glaciers under current and future changing climatic conditions. Model simulations with the physically based glacio‐hydrological model TOPKAPI‐ETH reveal a period of neutral or slightly positive mass balance between 2000 and 2010, followed by a transition to increasingly large annual mass losses, associated with a recent mega drought. Mass losses commence earlier, and are more severe, for a heavily debris‐covered glacier, most likely due to its strong dependence on snow avalanche accumulation, which has declined in recent years. Catchment runoff shows a marked decreasing trend over the study period, but with high interannual variability directly linked to winter snow accumulation, and high contribution from ice melt in dry periods and drought conditions. The study demonstrates the importance of incorporating local‐scale processes such as snow avalanche accumulation and spatially variable debris thickness, in understanding the responses of different glacier types to climate change. We highlight the increased dependency of runoff from high Andean catchments on the diminishing resource of glacier ice during dry years.


Introduction
Seasonal snow and glacier melt in the semiarid Chilean Andes provide water to more than two thirds of Chile's population as well as maintaining key economic activities, ecosystems and ecosystem services (Favier et al., 2009). Central Chile is characterised by warm and dry summers, and humid cold winters, and ice melt provides a key contribution to runoff in dry periods and during late summer and autumn, in a water balance otherwise dominated by snowmelt (Ragettli and Pellicciotti, 2012;Ohlanders et al., 2013;Ragettli et al., 2013;Ayala et al., 2016;Rodriguez et al., 2016). While glaciers in the region have been receding and losing mass over the past few decades (Rivera et al., 2002;Bown et al., 2008;Mernild et al., 2015;Malmros et al., 2016;Ragettli et al., 2016;Barcaza et al., 2017), the runoff response to climate and glacier changes is still poorly understood. Recent trends of decreasing runoff from high elevation catchments (Casassa et al., 2009;Mernild et al., 2016;Ragettli et al., 2016) suggest that the peak runoff, corresponding to the maximum contribution from a catchment (Pellicciotti et al., 2010;Huss and Hock, 2018), was reached at some time in the past. Results obtained from advanced glacio-hydrological modelling at relatively high resolutions Ragettli et al., 2016) and trend analysis (Casassa et al., 2009) are in agreement with global-scale models that also show declining trends in runoff from central Andean catchments (Bliss et al., 2014;Huss and Hock, 2018). However, these analyses, have been based on either a few intensively studied glaciers (Ragettli et al., 2016) or obtained from large scale studies with grid resolutions too coarse to capture differences caused by important local-scale processes (Mernild et al., 2016, Huss andHock, 2018).
Multidecadal studies focusing on the processes generating glacier streamflow response to a changing climate are needed to bridge this scale gap.
Recent studies in the region have advanced our understanding of the spatial patterns of ablation and glacier mass balance (Ayala et al., 2016(Ayala et al., , 2017Bravo et al., 2017), but they are generally limited to a maximum of a few seasons of data and none has investigated decadal changes of mass balance and runoff. Time series analysis of observational records (Casassa et al., 2009;Burger et al., 2018) is a useful tool to establish general data trends, but is of limited use when observations are scarce in space and time, and cannot provide insights into which processes drive observed changes. Additionally, whilet satellite-based glacier inventories (Rivera et al., 2002;Nicholson et al., 2010;Rabatel et al., 2011;Malmros et al., 2016;Barcaza et al., 2017) aided the establishment of baseline areal changes, they do not generally assess mass balance or volumes change, and cannot be used to explain the causes of observed changes. Therefore, there is a need for an integrated approach to understand the mid-and long-term changes in glaciers and glacier runoff in the high elevation catchments of the central Andes that combines both high resolution glacio-hydrological modelling and mass balance estimates from remote sensing (Pellicciotti et al., 2014), which are increasingly used to evaluate model simulations.
Determining glacier mass balance regimes and glacier hydrological contribution in the Andes is further complicated by the presence of debris-covered and rock glaciers, which account for approximately 3,200 km 2 of the 23,700 km 2 of inventoried ice (Barcaza et al., 2017). While the contribution of these glaciers to high elevation streamflow is poorly understood (Ayala et al., 2016), increasing evidence from other mountain regions shows a very distinct response of debris-covered glaciers to climate change compared to clean ice glaciers (Benn et al., 2012;Rowan et al., 2015). In one of the first glacio-hydrological modelling studies to explicitly include debris-covered glaciers, Ayala et al., (2016) showed that the contribution of a debriscovered glacier to total runoff over two years was of a similar magnitude to that of two debris-free glaciers in the same catchment.
Our main aims are to (1) reconstruct the glacier changes for the period 2000-2015 and (2) estimate the corresponding glacier contribution to runoff for the period. These aims are addressed through application of a physically-oriented and fully distributed glaciohydrological model, in-situ data, and the first geodetic mass balance estimates for the region.
We use this combination of modelling, field data and satellite observations to compare the hydrological contributions of a debris-covered glacier and two debris-free glaciers in the study catchment.
2 Study site and data

Study site
The Rio del Yeso catchment (33.55°S, 69.91°W, 3000-5230 m a.s.l, 62 km 2 ) is located ~70 km east of Santiago in the semiarid Andes of central Chile (Figure 1), and contains three principal glaciers: Bello, Yeso and Piramide. The former are small, debris-free valley-type glaciers (4.6 and 2.2 km 2 , respectively), while the debris-covered Piramide Glacier covers a larger elevation gradient although has a total area similar to that of Bello Glacier (4.7 km 2 ).
Piramide has a typical reverse ablation gradient and an estimated debris thickness ranging from 0.01 to 0.6 m (Ayala et al., 2016, colour-scale in Figure 1a). Mixed snow-debris avalanches typically feed the highest elevations of the glacier from local headwalls.

Meteorological data
The model was forced with meteorological data (temperature and precipitation) from Yeso Embalse Automatic Weather station (AWS) from the Chilean Water Directorate (Dirección General de Aguas, DGA) meteorological network (Figure 1b). This AWS recorded daily maximum and minimum temperatures and daily precipitation for the entire simulation period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). Additional temperature data from AWSs in the catchment were also used (Table 1). Hourly and daily maximum and minimum temperature from the Laguna Negra AWS (2780 m a.s.l.) ( Figure 1b) were used to identify the best disaggregation approach and evaluate the performance of the selected approach to derive hourly data from the temperature time series at Yeso Embalse AWS. Lapse rates used to extrapolate air temperature from Yeso Embalse AWS to the entire catchment were calculated using Yeso Embalse (2475 m a.s.l.), Yeso off-glacier (4300 m a.s.l.), Piramide off-glacier (3020 m a.s.l.), and Piramide on-glacier (3494 and 3655 m a.s.l.) AWSs (Figure 1c) between 2013 and 2015, using common data periods for each month (Table 2). Air temperature lapse rates were validated against hourly temperatures measured at AWSs installed on Bello and Yeso glaciers. Finally, daily cloud transmissivity coefficients were derived from reanalysis ERA-Interim data by considering constant values during the day.

Digital Elevation Models (DEMs)
Digital Elevation Models (DEMs) were used as both a basis for model runs (2000 SRTM at 30 m resolution), as well as to quantify glacier thinning through the 16-year simulation period. For the DEM differencing used to validate the model simulations we used the bistatic TanDEM-X and SRTM-C for the period 2000-2013, and two repeated airborne LiDAR surveys for 2012-2015. TerraSAR/X and TanDEM-X (TDX) correspond to an ongoing satellite constellation launched by the German Aerospace Center (DLR) and Airbus Defense and Space. TDX has a swath width of 30 km with a ground resolution of 0.4 arcsec. The Shuttle Radar Topography Mission (SRTM) was an interferometric synthetic aperture radar (InSAR) mission carried out simultaneously in the C-and X-band frequencies between 11 and 22 February 2000 between 56° S and 60° N (Farr et al., 2007). We used the void-filled LPDAAC NASA version of the SRTM DEM at 1-arcsec resolution. The two airborne LiDAR surveys carried out by the DGA have an estimated precision of ±0.30 m with an average of four points per square meter (ppm2) (DGA, 2012(DGA, , 2015. to April 2015. The system included a Canon EOS Rebel T3 camera with a resolution of 12.2 MP and a focal length of 18 mm. The camera was programmed to take photos every day at 13 h (local time). Photos were georeferenced following Corripio (2004) and albedo was derived from the terrestrial photos as explained in Ayala et al. (2016).

TOPKAPI ETH model
We used the TOPKAPI-ETH model to simulate glacier mass balance, glacier changes and runoff generation in the Río del Yeso catchment. TOPKAPI-ETH is an extended version of the original rainfall-runoff model TOPKAPI (TOPographic Kinematic APproximation and Integration) (Liu andCiarapica and and it has been applied in glacierized catchments from a few tens to more than 30,000 km 2 in the semiarid Andes (Ragettli and Pellicciotti, 2012;Ragettli et al., 2014a;Ayala et al., 2016), the Swiss Alps (Finger et al., 2011, 2012, Fatichi et al., 2014, 2015 and the Himalaya (Ragettli et al., 2013(Ragettli et al., , 2014b(Ragettli et al., , 2016. TOPKAPI-ETH offers a compromise between a detailed representation of high-mountain hydrological processes and computational efficiency. The model incorporates physically-based parameterizations of most relevant hydrological processes in high-mountain catchments, such as snow and ice melt (Pellicciotti et al., 2005), ice melt under debris (Carenzo et al., 2016), glacier area and elevation changes (Huss et al., 2010), snow albedo decay (Brock et al., 2000), gravitational redistribution of snow (Bernhardt and Schulz, 2010), and englacial storage and release of snow and ice meltwater (Hock and Noetzli, 1997).
In this study, we used the same TOPKAPI-ETH model setup as described in Ayala et al. (2016). Ayala et al. (2016) extensively calibrated and validated the model for the Rio del Yeso catchment using two years (2013)(2014)(2015) of field data that included manual snow depth measurements, ablation stakes, meteorological data from four AWSs, albedo time series from radiation measurements at Bello and Yeso AWSs, distributed fields of daily albedo derived from optical photos, and streamflow measurements at the outlet of Bello and Yeso glaciers.
To avoid error compensation and parameter ambiguity, the model was calibrated in a stepwise approach, in which each main parameter set was calibrated individually against specific field observations (see Figure 2 in Ayala et al., 2016). Here we perform an additional calibration step to account for the uncertainty in precipitation over the longer period of record of this study (see Section 3.5 below). Using the calibrated model setup, we then simulate glacier mass balance, elevation and areal changes for the period 2000 to 2015. Glacier volume and geometry changes were simulated using the Δh-parametrization developed by Huss et al. (2010). The Δh-parametrization is an empirical method that quantifies ice thickness changes as a function of previously observed elevation changes, and was developed based on a large datasets of glaciers in the Swiss Alps (Huss et al., 2010). Given the lack of repeated DEMs or ice volume observations for our region, we have used the set of parameters originally calibrated by Huss et al. (2010) for glaciers with an area smaller than 5 km 2 . The model was run for 10 years in a spin up mode to produce initial conditions of albedo and snow height and simulations started in 2000 with these initial conditions.

Extrapolation of meteorological variables
To drive the glacio-hydrological model, both air temperature and precipitation measurements are required at hourly resolution. While there was relatively good spatial coverage of AWSs in the study catchment in 2013-15, only a few stations were available over the complete study period (Table 1), and so data were extrapolated in both space and time. Hourly temperature time series were calculated from daily minimum and maximum values recorded at the Yeso Embalse AWS using the method suggested by Wilkerson et al. (1983), and adapted by Reicosky et al. (1989) (subroutine WCALC). The approach uses a sinusoidal function to interpolate between extreme values by dividing the day into three time periods: midnight to sunrise plus two hours; daylight hours; and sunset to midnight. It assumes that the minimum value occurs two hours after sunrise and the maximum halfway between sunrise and sunset time, obtained from the matlab subroutine Sunset, based on Montenbruck and Pfleger, (2000).
To evaluate the results of this approach in the study catchment, the method was tested at Laguna Negra AWS (LN, 2780 m a.s.l.; Figure 1b) using the data available for 2013-2015.
This site was selected because daily minimum and maximum data as well as independently recorded hourly data were available, and the station was at a similar elevation to the Yeso Embalse AWS (Table 1). The evaluation resulted in a Nash-Sutcliffe model efficiency (NSE) criterion of 0.93.
The interpolated hourly temperature data were distributed from the Yeso Embalse (YE) AWS ( Figure 1b) to the rest of the catchment using monthly means of hourly lapse rates, calculated from the available meteorological data in the catchment for each month and year with concurrent data (Table 2). In order to represent the effect of the glacier boundary layer (Greuell and Böhm, 1998;Brock et al., 2010), we use a parameter to decrease air temperature over glacier surfaces (Tmod) by subtracting 1°C for debris-free areas and a Tmod debris of 0.3°C to increase temperature for debris cover grid cells, calibrated for glacier-covered areas by Ayala et al. (2016). Figures 2 and 3 show the disaggregated and extrapolated temperatures at the Bello ( Figure 2) and Yeso (Figure 3) AWSs. In general, calculated temperatures correspond well to measured values, especially during the daytime in the summer months; however the disaggregation and extrapolation method performs less well on glacier surfaces before sunrise when there snow was present, and during winter time when air temperatures were below 0° C.
Precipitation is not evenly distributed across the study catchment, and according to previous studies in the region, a logarithmic model can be used to represent precipitation spatial variability at high elevations (Vicuña et al., 2010;Ragettli et al., 2014a). Thus, we extrapolated the hourly precipitation measurements from Yeso Embalse AWS for the period 2000 to 2015 using a logarithmic model as follows (Ragettli et al., 2014a): Where P(Yeso Embalse) is precipitation measured at Yeso Embalse, and P(z) the precipitation at the elevation "z". The values of the coefficients in Equation 1  reproduced. To represent local effects on precipitation and its accumulation over glaciers, we used a local scaling factor that modifies precipitation over each glacier (Huss et al. (2008); Magnusson et al. 2011;Farinotti et al., 2012). This factor was calibrated to match the simulated long-term glacier elevation changes with those derived from the geodetic mass balance to account for local processes governing snow accumulation on glaciers. Use of a local scaling factor for each glacier was supported by evidence from field observations of preferential deposition, scouring, and snow removal by wind, that cannot be captured by a regional precipitation gradient. The local factor was calibrated against the elevation change of the period 2000-2013 obtained from the DEM differencing. The estimated factors were 0.74 for Bello and Yeso, and 1.88 for Piramide. A similar approach has been used by Magnusson et al. (2011) andFarinotti et al., (2012), and in light of the lack of local observations, it is a plausible way of preventing precipitation uncertainty from dominating the modelling exercise. Finally, precipitation was disaggregated to hourly values by distributing daily values homogeneously during the day.

Debris Thickness estimation
We used the debris thickness map derived by Ayala et al. (2016) for the debris-covered areas on the Piramide, Bello and Yeso glaciers ( Figure 1). The map was derived by solving the distributed energy balance of the debris-covered areas at the moment of acquisition of a Landsat 8 thermal image of the area at 90-m spatial resolution. This method was originally presented by Foster et al., (2012) as a physically-oriented alternative to empirical relationships between surface temperature and debris thickness (e.g. Mihalcea et al., 2008).
Different versions of the method have been subsequently presented (e.g. Rounce and McKinney, 2014) but several uncertainties are still associated with this approach, such as the debris temperature profiles, heat storage rate and turbulent heat fluxes (Schauwecker et al., 2015;Ayala et al., 2016). Technical details regarding the development of the debris thickness map can be found in Appendix 1 in Ayala et al. (2016).

Geodetic elevation change
In recent years, the geodetic method has been widely used to obtain glacier changes over short-or long-term periods (e.g. Rankl and Braun, 2016;Bolch et al., 2017). The method, based on the differencing of digital elevation models (DEMs) can provide glacier changes over several years for large remote areas.
We followed the TDX processing scheme in Malz et al. (2018) using the GAMMA software.
We acquired co-registered single look complex images (CoSSC) in HH polarization provided by DLR. All the TanDEM-X scenes along one track were concatenated and the strips processed by differential interferometry (DInSAR) (e.g. Malz et al., 2018). We used the SRTM-C DEM as a reference and subtracted it from the bistatic TanDEM-X interferogram.
For this, the topographic phase was simulated from the SRTM-C DEM using the TanDEM-X orbit parameters. We unwrapped the interferograms using a Minimum Cost Flow algorithm.
In order to remove the phase noise from the differential interferogram we applied a Goldstein filter, and areas with low coherence (coherence < 0.2) were masked out (Goldstein and Werner, 1998;Vijay et al.;Rankl and Braun, 2016). The unwrapped differential phase was converted into absolute differential heights. These differential heights were added back to the topographic heights from SRTM-C DEM to generate a TanDEM-X DEM. The resulting TanDEM-X DEMs were geocoded with the SRTM-C DEM to maintain planimetric consistency (e.g. Vijay et al., 2016;Malz et al., 2018).
The post-processing comprises the mosaicking of all raw DEMs resulting from the InSAR strip processing. We used a stable ground mask derived from optical data (Landsat OLI 2013) and corrected vertical biases between the strips applying a polynomial fitting (Malz et al., 2018). TanDEM-X DEMs were iteratively co-registered (vertically and laterally) to the SRTM-C DEM using the approach of Nuth and Kääb (2011).
Uncertainties in the geodetic elevation changes for the period 2000-2013 were estimated by calculating the Median Absolute Deviation (MAD) for the elevation differences on stable areas. Since the deviation is known to be slope-dependent , the area of interest was divided in 5° slope-bins and the total MAD was calculated by weighting the area of each bin (e.g. Malz et al., 2018). We discarded any significant bias associated to the radar signal penetration in snow and ice of SRTM-C and TanDEM-X, as previous studies have shown that summer DEM acquisitions in the Southern Hemisphere (during melting conditions) reveal negligible penetration (Jaber et al. 2013;Jaber et al. 2016;Falaschi et al., 2017;Dussaillant et al., 2018). MODIS. We compared the SCA derived from a daily MODIS MOD10A1 product (Hall et al., 2002) to the catchment-wide snow cover area derived from the TOPKAPI-ETH on a daily timestep. The daily MODIS data were discarded if more than 10% of the total area was covered by cloud, and any remaining cloud covered cells were filtered using a linear interpolation of SCA quantity based upon a temporal search window of two days either side of the cloud cover observation at the given cell. The MODIS grids were resampled to a 30 m grid and clipped to the same area as the model domain (see blue line in Figure 1c), and for each day, the catchment average SCA was extracted and compared to that of the TOPKAPI-ETH model simulations. The validation procedure of albedo and SCA using the terrestrial camera is detailed in Ayala et al. (2016) and not repeated here.

Model validation results
The ground-based and satellite validations of the TOPKAPI-ETH model are given in Figures   5 and 6, respectively. The model captures the variability in albedo measured at Bello and Yeso AWSs and both the albedo and SCA derived from the terrestrial camera ( Figure 5).
Specifically, modelled albedo at Bello Glacier AWS follows the measured decay rates and albedo increase after spring storm events (Figure 5b)  (Figures 9a-b). However, geodetic values were more negative than the 2000-2013 simulations between 4200 and 4800 m a.s.l. (Figures 9a-b).
The differences between the debris-free and debris-covered glaciers are also evident in the spatial distribution of the geodetic mass balance, with a profile that is quite smooth for Bello (and slightly less so for Yeso) while the spatial variability at Piramide is very high, a feature also well captured by the model. According to the geodetic mass balance between 2000 and 2013, simulations in Bello and Yeso overestimate melt below ~4300 m a.s.l., and overestimate accumulation above that elevation, while at Piramide the geodetic mass balance pattern for that period is similar to the simulations, with a small melt overestimation between 3500 and 3800 m a.s.l.

Runoff and runoff components
Total runoff and the annual contribution to runoff from snowmelt, icemelt and precipitation are shown in Figure 10. A clear overall decline in runoff can be observed, which is marked for the period after 2009. The main runoff contribution comes from snowmelt, accounting for is highest in dry years with low total runoff (such as 2014-2015, which have the highest proportional contribution from icemelt). The liquid precipitation contribution is consistently small and does not exceed 6% in any given year. While annually snowmelt represents the main water input to the system, in summer (January and February) icemelt becomes increasingly dominant, contributing equally with snowmelt by March (Figure 11). The contribution of icemelt to total runoff from each of the three glaciers is distinct (Figure 12). snowmelt contribution from all three glaciers has reduced compared to the previous years (orange and red lines on the right hand panels in Figure 12). The seasonality is also distinct from 2013 onwards, with snowmelt and icemelt both occurring earlier at Piramide (Figure   12e and f) than on the other two glaciers (Figure 12a-d).

Mass balance and runoff contribution
Geodetic mass balance estimates show a very distinct behaviour in the two analysed periods (2000-2013 and 2013-2015). While an almost neutral elevation change was obtained for the first period (-0.01 to -0.14 m yr -1 ), a high ice thinning rate is evident for the second period (-  (Garreaud et al., 2017), termed the "Mega-drought" (Boisier et al., 2016) which has been characterized by historically low precipitation levels, shallow seasonal snowpacks (Cornwell et al., 2016;Cortés and Margulis, 2017) and high temperatures (Garreaud et al., 2017), especially in spring and autumn (Burger et al., 2018).
There is a good agreement between the simulated and geodetic mass balance at Piramide Glacier ( Figure 8). The model does not fully reproduce the values from the geodetic mass balance at Bello and Yeso glaciers in 2015, though they are within the range of uncertainties of the geodetic measurements. Simulations do show a decreasing trend starting in 2010, but the geodetic mass loss in the period 2010-2015 is larger than that from the model. Despite these differences, a more detailed analysis of the patterns of glacier surface change shows that the model is able to reproduce some of the elevation-dependent differences evident in the geodetic elevation changes (Figure 9). Importantly, the model is able to resolve much of the spatial variability of glacier surface change in relation to the differences in debris thicknesses, as well as those related to the avalanches that feed the debris-covered Piramide glacier.
The glaciers analysed in this study provide an excellent example of the different processes that affect the long-term evolution of debris-free and debris-covered glaciers within one catchment. While the mass balance of Bello and Yeso is more strongly controlled by temperature gradients affecting the precipitation phase and ablation components of the model (Ayala et al., 2016) and snow removal on the mid glacier (evidence from unpublished data on Bello), the mass balance of Piramide Glacier is controlled principally by debris thickness and snow accumulation from avalanches ( Figure 9). A sensitivity analysis revealed that ignoring the contribution of avalanching in the model simulations produces a total cumulative elevation change more than three times more negative for Piramide compared to our reference model (Figure 13). Ignoring this process has a greater impact on the total absolute glacier elevation change than artificially providing 10 centimetres more surface debris for the glacier (red triangles in Figure 13). The strong control that debris thickness exerts over the glacier mass balance has been extensively demonstrated in the literature (e.g. Ragettli et al., 2016;Vincent et al., 2016) and is shown by the results from our study to be important for long term modelling of debris-covered glaciers mass balance.
The removal of snow from the central sections of Bello Glacier would contradict a hypothesis that glacier areas act as net sinks of snow during the accumulation season (Dadic. et al., 2010;Gascoin et al., 2013). However, the mass losses that seem to occur on Bello and Yeso glaciers might be compensated by the contribution of avalanches from the surrounding upper slopes, which is still shown to be an important process for these glaciers ( Figure 13).
Evidence of negative surface changes in the central sections of Bello ( Figure 7a) reveals a potential local process of wind effects on snow which agrees with field observations, though modelling snow redistribution on high elevation glaciers remains an important area for future studies, even if outside of the scope of this paper.
Runoff generation from debris-free and debris-covered glaciers also exhibit a distinct behaviour. Particularly interesting is the different interannual spread in the runoff generated by icemelt from Bello and Piramide glaciers ( Figure 12). Icemelt from Bello Glacier shows a large interannual variability, which most likely indicates a sensitivity to climatic variability of precipitation and temperature. In turn, icemelt from Piramide Glacier is largely insulated from climatic changes due to its supraglacial debris. Interestingly, the opposite pattern is evident for snowmelt, i.e. large interannual variability on Piramide Glacier in comparison with Bello Glacier (Figure 12f). In this case, the predominant low elevation of Piramide Glacier makes the snowpack sensitive to air temperature and the proportional amounts of solid winter precipitation. Further still, the mass loss characteristics of a debris-covered glacier such Piramide Glacier is strongly governed by the presence of ice cliffs (Steiner et al., 2015;Buri et al., 2016) which are currently not considered in TOPKAPI-ETH. Future modelling studies of debris-covered glaciers may therefore benefit from physically-based or parameterised representations of these processes.
As demonstrated by previous studies for individual years, snowmelt is the largest contributor to runoff in Andean catchments of central Chile with an outlet point around 3000 m a.s.l. (Ragettli and Pellicciotti, 2012;Ohlanders et al., 2013). In this study, we showed that this result holds for a long time period, although with important interannual variability ( Figure   10). Particularly important is the evidence that icemelt contribution becomes more relevant during drought periods. We observed an increase in the relative contribution to runoff of ice melt from 11.6% in the period 2000-2009 to 20.5% in the period 2010-2015. Quantifying the shift in relative streamflow contributions for a drier climate is highly relevant to the socioeconomy of central Chile, and our study provides a new insight into the longer term evolution of these contributions, and most significantly a shift to increasing dependency on the declining resource of glacier ice.

Sources of uncertainty
As suggested by Ayala et al. (2016), the spatial distribution of forcing variables and the debris thickness are relevant controls on the glacier mass balance in this catchment.
Particularly challenging is the modelling of snow accumulation on debris-free glaciers and wind-exposed locations. As our simulations show, the cumulative elevation change of Bello and Yeso glaciers was consistently overestimated when using regional precipitation gradients that provide a good representation of annual water balances in the region. More precise regional estimates of precipitation and the simulation or parameterization of snow transport would likely improve the simulation of glacier mass balance. On the other hand, more accurate estimates of the current debris thickness distribution and its time evolution will become necessary for future studies on long-term glacier mass balances (Rowan et al., 2015).
However, changes in debris thickness over a period of 16 years are likely minor or restricted to specific areas.

Conclusion
In this paper, we have used a combination of distributed glacio-hydrological modelling, new estimates of geodetic mass balance for the region and a large amount of field data over two seasons to investigate the mass balance and runoff contribution of the glaciers of the Rio del Yeso catchment. Our main conclusions are: • • The spatial distribution of the mass balance over the debris-covered glacier is distinct to that of the debris-free glaciers, and its elevation profile is related most strongly to the debris thickness variability and avalanches contribution. There is also a contrast in behaviour in terms of runoff contribution between the debris-covered and debris-free glaciers. The interannual variability of snowmelt contributions from Piramide Glacier is much larger and more sensitive to the 0°C isotherm than the debris-free Bello and Yeso glaciers. By contrast, there is a less variable interannual sub-debris icemelt on Piramide Glacier that is decoupled from high frequency climate variability.
• We witness a clear decrease in runoff over the period of study, with very low runoff during the years of the mega-drought, but a decline which is evident from the start of the study period in 2000. Our period of record is too short to confirm a general longterm trend, such as those that have been modelled or suggested by other studies in the region. Dry years show an increased dependency of runoff on the declining resource of glacier ice.
Given this result and those of the global studies, there is a clear need to extend model simulations and reconstruction of geodetic elevation changes for debris-free and debriscovered glaciers to longer time periods, in order to establish whether the peak water has been reached already and what the contribution of distinct types of glaciers is.   (2015)