Glacier runoff variations since 1955 in the Maipo River basin, in the semiarid Andes of central Chile

As glaciers adjust their size in response to climate variations, long-term changes in meltwater production can be expected, affecting the local availability of water resources. We investigate glacier runoff in the period 1955–2016 in the Maipo River basin (4843 km2, 33.0–34.3 S, 69.8–70.5W), in the semiarid Andes of Chile. The basin contains more than 800 glaciers, which cover 378 km2 in total (inventoried in 2000). We model the mass balance and runoff contribution of 26 glaciers with the physically oriented and fully distributed TOPKAPI (Topographic Kinematic Approximation and Integration)-ETH glacio-hydrological model and extrapolate the results to the entire basin. TOPKAPI-ETH is run at a daily time step using several glaciological and meteorological datasets, and its results are evaluated against streamflow records, remotely sensed snow cover, and geodetic mass balances for the periods 1955–2000 and 2000–2013. Results show that in 1955–2016 glacier mass balance had a general decreasing trend as a basin average but also had differences between the main sub-catchments. Glacier volume decreased by one-fifth (from 18.6±4.5 to 14.9±2.9 km3). Runoff from the initially glacierized areas was 177±25 mm yr−1 (16±7 % of the total contributions to the basin), but it shows a decreasing sequence of maxima, which can be linked to the interplay between a decrease in precipitation since the 1980s and the reduction of ice melt. Glaciers in the Maipo River basin will continue retreating because they are not in equilibrium with the current climate. In a hypothetical constant climate scenario, glacier volume would reduce to 81± 38 % of the year 2000 volume, and glacier runoff would be 78± 30 % of the 1955–2016 average. This would considerably decrease the drought mitigation capacity of the basin.

69.8-70.5 • W), in the semiarid Andes of Chile. The basin contains more than 800 glaciers, which cover 378 km 2 in total CE1 (inventoried in 2000). We model the mass balance and runoff contribution of 26 glaciers with the physically oriented and fully distributed TOPKAPI-ETH CE2 glacio-10 hydrological model and extrapolate the results to the entire basin. TOPKAPI-ETH is run at a daily time step using several glaciological and meteorological datasets, and its results are evaluated against streamflow records, remotely sensed snow cover, and geodetic mass balances for the periods

Introduction
Most glaciers on Earth have retreated due to global atmospheric warming during the 20th century (Zemp et al., 2019). Glaciers that are still out of balance with the present cli-35 mate are committed to lose part of their mass in the coming decades, even without further warming (Zemp et al., 2015;Marzeion et al., 2018), and major changes in their meltwater production can be anticipated Huss and Hock, 2018;IPCC, 2019). In the absence of precipitation 40 changes, a temporary increase in meltwater generation from a retreating glacier occurs as a consequence of higher air temperatures and enhanced ablation, but after this transient phase melt amounts decrease due to the reduction of the available snow, firn and ice volumes (Jansson et al., 2003). The pe- 45 riod in which the annual melt volume reaches its long-term maximum has been termed "peak water" (Gleick and Pala-2 Á. Ayala et al.: Glacier runoff variations since 1955niappan, 2010Baraer et al., 2012). Global-scale studies indicate a large heterogeneity in the geographical distribution of peak water, while several catchments in the Himalayas and Alaska are expected to increase their glacier runoff due to the enhanced ablation in the next decades and reach a maximum at some point of the 21st century, other regions in the world, such as the semiarid Andes, central Europe, and western Canada, have already reached a regional maximum, and thus glacier runoff will only decrease in the future Huss and Hock, 2018). While these studies provide 10 global trends that are key for macro-regional assessments, studies focusing on the catchment scale can provide more specific information about local hydro-glaciological changes to communities and stakeholders for the generation of mitigation and adaptation strategies. Additionally, catchment- 15 scale studies place glacier runoff in the context of other components of the water cycle and evaluate the impacts of glacier changes on downstream areas.
In this study, we focus on glacier changes and their impacts on long-term glacier runoff contribution in the semiarid 20 Andes. Meltwater originating in the Andes is key for Chile and the western areas of Argentina, as it represents the main source of water for drinking water, agriculture, industry, mining, and natural ecosystems. The climate of this region is characterized by its strong interannual variability of precip- 25 itation linked to periodic atmosphere-ocean variations over the Pacific Ocean (Montecinos and Aceituno, 2003;Falvey and Garreaud, 2007) and a sustained air temperature increase during the last decades (Carrasco et al., 2005;Burger et al., 2018). A few studies have estimated the present (Ragettli and 30 Pellicciotti, 2012; Ayala et al., 2016;Burger et al., 2019) and future (Ragettli et al., 2016;Huss and Hock, 2018) glacier runoff contribution in the semiarid Andes, but its past variations have not been analyzed in detail, mostly due to the lack of long-term glaciological data. As future climate sce-35 narios anticipate a decrease in glacier runoff (e.g., Ragettli et al., 2016), the question of whether peak water has already occurred still remains open. The assessment of long-term changes in glacier runoff is particularly useful for water planners because it provides reference information for the role of 40 glacier meltwater in river flows, and the impacts that can be anticipated in the absence of its contribution.
Glaciers in the semiarid Andes underwent a major retreat in the 20th century (Le Quesne et al., 2009;Malmros et al., 2016) and the last 2 decades (Braun et al., 2019;Dussaillant 45 et al., 2019). Historical documents, aerial photographs, and dendrochronological studies suggest that the general retreat trend started around the mid-19th century, but it has been interrupted by occasional periods of positive mass balance accompanied by glacier advances (Le Quesne et al., 2009;50 Masiokas et al., 2009). Masiokas et al. (2016) performed a reconstruction of the annual mass balance of Echaurren Norte Glacier (3650-3900 m a.s.l.) since 1909 using a simple glacier mass balance model forced with monthly precipitation and air temperature. The model was verified against 55 streamflow records and direct mass balance measurements on the glacier, where the first mass balance monitoring program in the southern Andes started in 1975. Masiokas et al. (2016) found a general retreat interrupted by three periods of sustained positive mass balance in the 1920-1930s, 1980s, 60 and 2000s. The latter positive or balanced mass budget in the semiarid Andes in the 2000-2009 period has been recently verified by geodetic mass balances (Braun et al., 2019;Dussaillant et al., 2019;Farías-Barahona et al., 2019) and has been supported by independent modeling results (Burger et 65 al., 2019). As the findings by Masiokas et al. (2016) are based on a relatively simple model applied to only one glacier at low elevation (< 4000 m a.s.l.), they cannot be extrapolated to other glaciers. This is especially true due to the large spatial variability of response times and retreat rates reported 70 for this region (Malmros et al., 2016). Thus, a more detailed analysis based on the specific characteristics of each glacier is needed to complement these results and estimate regional changes of ice volume and glacier runoff. From a climatic perspective, glacier retreat in the semiarid Andes 75 has been driven by a general temperature increase and modulated by a strong temporal variability of precipitation. Air temperature showed an increasing trend of about 0.25 • C per decade in the period 1979-2006(Falvey and Garreaud, 2009, mostly explained by a spring and autumn warming 80 (Burger et al., 2018), which can be used to explain an increase in the 0 • C isotherm and the regional equilibrium line altitude (ELA) (Carrasco et al., 2005(Carrasco et al., , 2008. Precipitation, on the other hand, exhibited an average decrease of −65 mm (−7.1 %) per decade in the period 1979-2016 (Boisier et al., 85 2016), although it has a large interannual and inter-decadal variability (Montecinos and Aceituno, 2003).
Our main objectives are to reconstruct glacier changes (area and volume) during the last 6 decades in one of the main catchments of the semiarid Andes, the Maipo River 90 basin; analyze the role of glaciers in the regional hydrology; and identify the main trends in glacier runoff. Glacier runoff is defined as the water originating from ice melt, snowmelt and rain over a given glacier. Additionally, we estimate glacier changes under synthetic scenarios of committed 95 ice loss, in which air temperature, precipitation and cloudiness are assumed to stay at their current levels until the end of the century. We use these scenarios for (i) understanding how far the glaciers are from an equilibrium after the climatic changes that took place in the period 1955-2016 and 100 (ii) providing a baseline for the future changes in hydrology that the basin will experience in any case, i.e., even in the hypothetical case that climate change was to stall. They are thus highly conservative and do not correspond to a realistic projection for the future. The calculation of glacier 105 changes and runoff contribution is carried out for a subset of the largest glaciers using the physically oriented and fully distributed TOPKAPI-ETH glacio-hydrological model (Ayala et al., 2016;Ragettli et al., 2016), and the resulting mass balances are extrapolated to the entire basin (Huss, 2012). We 110 set up the glacier model using glacier inventories, digital elevation models (DEMs), and estimates of ice thickness, and we force it with a combination of local meteorological stations and reanalysis data for precipitation, air temperature, and solar radiation. The model is calibrated and validated 5 using remotely sensed snow cover, streamflow records, and geodetic mass balances covering the periods 1955-2000and 2000(Braun et al., 2019Farías-Barahona et al., 2020).

Study area
The study focuses on the headwaters of the Maipo River 10 basin (hereafter we refer to these areas as the Maipo River basin for simplicity). The basin is located in central Chile (∼ 33 • S, ∼ 70 • W), to the east of the Chilean capital city, Santiago (Fig. 1a), to which it provides about 70 % of its drinking water (DGA, 2004). The basin outlet is the Maipo en El Manzano gauging station, which roughly marks the boundary between rural mountain areas and Santiago urban districts. The selected basin has an area of 4843 km 2 , its elevation ranges from 850 to 6570 m a.s.l., and more than 800 glaciers covering about 378 km 2 (7.8 % glacierized) 20 were inventoried in 2000 (Barcaza et al., 2017). The Maipo River and its tributaries are the primary source for drinking water, agriculture, hydropower, and industry in the region, which concentrates about 40 % of the country's population. The region has a Mediterranean-type climate, with 25 a strong seasonality characterized by cold and wet winters and hot and dry summers. Average precipitation in Santiago was 308 mm yr −1 in the period 1950-2018, but values as low as 69 mm yr −1 and as high as 712 mm yr −1 have been registered, with a coefficient of variation of 0.45. Recurrent 30 droughts have been reported since the beginning of hydrometeorological records. Precipitation amounts are, in general, larger towards the south and towards higher elevations. An early study estimated that glacier runoff in the Maipo River basin represents about 34 % of the total discharge in Febru-35 ary and up to 67 % during summer months of dry years, such as 1968-1969(Peña and Nazarala, 1987. There are five major sub-catchments in the study area, from north to south: Olivares, Colorado, Yeso, Volcán, and Upper Maipo (Fig. 1b). According to the national Chilean 40 inventory (Barcaza et al., 2017) (described in next section), the highest glacierized sites are in the Olivares and Colorado sub-catchments, with mean elevations between 4200 and 4500 m a.s.l., and some glaciers reaching elevations higher than 5500 m a.s.l. (Fig. 1c). The Upper Maipo sub- 45 catchment, on the other hand, has the lowest-lying glaciers, with mean elevation varying between 3500 and 4000 m a.s.l., and several glaciers reaching elevations below 3000 m a.s.l. (Fig. 1c). Glacierized areas vary from 40 km 2 in the Volcán sub-catchment to 99 km 2 in Colorado. Upper Maipo 50 has the largest number of individual glaciers (348), and most of them correspond to low-elevation, rock, and debris-covered glaciers (Fig. 1d). In general, glacier size tends to decrease towards the south, with the largest glaciers being located in Olivares (Juncal Sur, Olivares Gama, and Olivares 55 Beta glaciers) and on the slopes of Tupungatito and Marmolejo volcanoes (Volcán Tupungatito, Azufre, and Marmolejo glaciers) in the Colorado sub-catchment. Another series of relatively large glaciers corresponds to debris-covered ones, such as Pirámide, Loma Larga, and Cerro Castillo 60 glaciers.

Geographic and topographic information
Glacier outlines are extracted from the national Chilean inventory (Barcaza et al., 2017) and the Marangunic inventory 65 (Marangunic, 1979). While the information for the Maipo River basin in the national inventory was produced using two satellite images from the Landsat Enhanced Thematic Mapper ETM+ of 2003, the Marangunic inventory was mostly based on aerial photographs taken in 1955 during a national 70 geodetic program and maps presented by Lliboutry (1956) for the few missing areas. For consistency with the DEM obtained from the Shuttle Radar Topography Mission (SRTM), we assume that the outlines in the national inventory from 2003 are also valid for 2000. Additionally, the glacierized ar-75 eas in the national inventory that are not identified as such in 1955 (mostly rock glaciers and debris-covered areas) are added to the 1955 inventory. In this study, we assign an error of 5 % to the year 2000 inventory, which is a common choice for glacier inventories (Paul et al., 2013), and has been used 80 for this inventory in particular (Barcaza et al., 2017). As the inventory of 1955 suffers from additional errors (such as the presence of snow patches that likely made the interpretation of glacierets difficult, and the use of Lliboutry maps to fill missing areas), we assume an error of 10 % for that year. 85 Based on the resulting inventories, we estimate that the total glacier area changed from 532 km 2 in 1955 to 378 km 2 in 2000 (−28.9 %) and that the number of individual glaciers decreased from 861 to 854. Although some small glaciers might have effectively disappeared, the decreasing trend in 90 the total number of glaciers is also balanced by the fragmentation of large glaciers, such as the Olivares Alfa glacier complex, into several smaller units (Malmros et al., 2016).
In addition to the glacier inventory, we generate a mask of debris-covered glacier areas from the same Landsat im-95 ages that were used to produce the Chilean glacier inventory. For this, we use the semiautomatic method based on band ratio segmentation of TM4 and TM5 Landsat bands (Paul et al., 2004), and we manually correct the results using Google Earth imagery. For the year 1955, we maintain the same de-100 bris cover maps as in the 2000-2010 period, i.e., assuming that no major changes have occurred in the extension of de- Here, we provide a brief description of the derivation of these datasets, but more details are included in the Supplement. The 1955 DEM was calculated 15 from digitized 50 m contour lines of the 1 : 50 000 official Chilean cartography product, which was also obtained from the 1955 geodetic program. While the DEMs for the year 2000 were extracted from the SRTM product, the DEM of Maipo River basin for 2013 was derived from TanDEM-X 20 post-processed products (which for this region correspond to the year 2013). The DEMs were co-registered following Nuth and Kääb (2011). Errors from the geodetic mass balances were assessed over stable ground, and calculated us-ing a standard error propagation procedure, including typical 25 error sources such as radar penetration signal. Two glaciers (San Francisco and Mirador del Morado) were discarded from the geodetic mass balance because the original SRTM product was not available for those areas (only the void-filled product). As rock glaciers exhibit changes that are smaller 30 than the estimated uncertainties, they were also discarded from the geodetic mass balance.

Ice thickness
Distributed glacier ice thickness in 2000 is estimated for all individual glaciers using the method of Huss and 35 Farinotti (2012) with the glacier outlines and the SRTM DEM. Standard model parameters are used except for glaciers classified as debris-covered or rock glaciers. For these two types of ice bodies, the parameter prescribing ice flux is substantially reduced to obtain thicknesses compara-40 ble to the direct thickness observations on the debris-covered Pirámide Glacier and its neighboring rock glaciers (DGA, 2012). The obtained ice thickness estimates compare well with ground-penetrating radar (GPR) measurements (DGA, 2014) on Volcán Tupungatito (1685 data points) and Mar-45 molejo (1544 data points) glaciers extracted from the Glacier Thickness Database (GlaThiDa) (Gärtner-Roer et al., 2014), for which we find a root-mean-square error (RMSE) of 9.8 and 8.5 m, respectively.
Once the distributed ice thickness is calculated for every glacier for the year 2000, we use the geodetic mass bal-5 ance in the period 1955-2000 to estimate the ice thickness distribution in 1955. In this procedure, we find the problem that for some grid cells showing a positive elevation change from the geodetic mass balance for the 1955-2000 period, the ice thickness in year 2000 is too small, resulting in an 10 inferred negative thickness. To avoid this and obtain meaningful 1955 ice thicknesses that are consistent with both the geodetic mass balance and the glacier inventory, we assign the year 2000 ice thickness to 1955 in these grid cells and add the estimated positive elevation change. In this way, we 15 obtain a corrected ice thickness value in 2000 for 4.8 % of the glacierized area. As no geodetic mass balance was calculated for rock glaciers, they are assumed to have the same thickness in 1955 as in year 2000. A similar result was found in the study of Bodin et al. (2010) for rock glaciers near San-20 tiago. Finally, to calculate the 1955 ice thickness of small glaciers that are not included in the year 2000 inventory, we use the 1955 glacier areas from the glacier inventory and a scaling relation to calculate mean ice thickness (h) as a function of the glacier area (S) and assume average thickness to 25 be valid for every grid cell in these glaciers: where γ = 1.357 and c = 28.5 are standard parameters in the volume-area scaling theory (Chen and Ohmura, 1990). At the basin scale, we find a total ice volume of 18.6 ± and 2013 is larger than in 2000 since it also includes the uncertainty from the geodetic mass balances. In the calculation of glacier volumes, we implicitly assume that no basal melting takes place. The error introduced by neglecting this process is much less than the uncertainty associated with the ice 45 thickness estimates and the geodetic mass balance.

Hydrometeorological data
Precipitation and temperature data for the period 1979-2016 are derived from daily gridded products developed by the Centre for Climate and Resilience Research in Santi- 50 ago, Chile (CR2, http://www.cr2.cl TS3 ). These products were generated for a national water balance study led by the Chilean directorate of water resources (DGA) (DGA, 2017;Álvarez-Garretón et al., 2018). The CR2 daily precipitation product was generated by means of a statistical downscaling 55 of precipitation and moisture fluxes from the ERA-Interim reanalysis. The downscaling procedure is based on multiple linear regressions with topographic parameters, which were calibrated with quality-controlled precipitation records. The CR2 temperature product was obtained using near-surface 60 temperature from ERA-Interim and land surface temperature (LST) from the Moderate Resolution Imaging Spectroradiometer (MODIS), by means of multiple regression models using LST as the explanatory variable and validated with local observations. For our study, while the CR2 precipitation 65 product is linearly interpolated from its original resolution (0.05 • ) to the spatial resolutions of our glacio-hydrological models (1 km and 100 m; see Sect. 4.1.2 and 4.1.3) to generate monthly average maps, the CR2 temperature product is used to generate basin-scale daily temperature lapse rates. 70 Daily cloud transmissivity of solar radiation is calculated from the Chilean solar radiation database (http://www. minenergia.cl/exploradorsolar/TS4) for the period 2004-2016 at the location of Embalse El Yeso meteorological station, which is placed close to the centroid of the Maipo River 75 basin, and assumed to be uniform over the catchment. The solar radiation database was derived using reanalysis data to force a radiative transfer model for clear-sky solar irradiance and an empirical model based on satellite data for cloudy conditions (Molina et al., 2017).

80
In addition to the information from the CR2 products, we use local records of air temperature and precipitation from Embalse El Yeso and Quinta Normal (located in Santiago) meteorological stations, respectively, as a base for extrapolating these variables during the period 1955-1978 (see 85 Sect. 4.1.2 and 4.1.3). Values for air temperature gradients and cloud transmissivity in the study periods without information from CR2 and the Chilean solar radiation database (1955 to 1978 and 1955 to 2003, respectively) are randomly selected from a pool of values recorded in the same day of 90 the year in the periods with available information. Finally, streamflow data for the Maipo River basin are available as monthly mean records at the gauging station of Maipo en El Manzano. These time series were already corrected for extractions and reservoirs to approximate the natural flow in 95 the study of CONIC-BF (2008).

Additional datasets
To calibrate and validate the snow processes in the study area, we use two products: (i) post-processed MODIS snow cover area (SCA), downloaded from an online platform (http:// 100 www.dgf.uchile.cl/rene/MODIS/TS5) that automatically calculates SCA from MODIS Terra and Aqua satellite products at a spatial resolution of 500 m in several Chilean river basins and (ii) daily basin-scale snow water equivalent (SWE) estimates for the period 1984-2014, extracted from 105 6 Á. Ayala et al.: Glacier runoff variations since 1955 the Chilean version of the Catchment Attributes and Meteorology for Large Sample Studies (CAMELS-CL) database (http://camels.cr2.cl/TS6). These basin-scale SWE estimates were aggregated by Alvarez-Garretón et al. (2018) from a daily gridded SWE reconstruction for the Andes Cordillera generated by Cortés and Margulis (2017) at a 180 m resolution. The SWE reconstruction was obtained from a data assimilation framework that integrates a land surface and depletion model, the assimilation of Landsat imagery, and the Modern Era Retrospective Analysis for Research and Appli-10 cations (MERRA) reanalysis as a forcing dataset (Cortés et al., 2016;Cortés and Margulis, 2017). Although not all physical processes are included in the assimilation process (for example, blowing snow sublimation), the dataset has been validated at several sites across the southern Andes (Cortés 15 et al., 2016;Cortés and Margulis, 2017), and it should provide a good estimate of snow on the ground that can be used for hydrological modeling.
For modeling evapotranspiration and subsurface water fluxes, we generate land use and soil types maps, respec-20 tively. The land use maps are extracted from the National Forest Corporation (CONAF) database (CONAF, 2013), and the same maps are used to estimate the spatial distribution of soil types in the basin. For simplicity, and due to the absence of more detailed data, we define only two soil types based 25 on the presence or absence of vegetation. The vegetated soil type dominates areas at low elevations and close to streams, whereas the non-vegetated soil type dominates on mountain slopes. To our knowledge, there are too few detailed datasets to evaluate changes in land use throughout the study period, 30 and we keep land use and soil types constant in our simulations.

35
TOPKAPI-ETH is a physically oriented, fully distributed, glacio-hydrological model that was adapted from a rainfall runoff model (Ciarapica and Todini, 2002) to simulate snow cover evolution and glacier mass balance in high mountain areas. The model has been used successfully in the semiarid 40 Andes (Ragettli et al., 2014;Ayala et al., 2016), the Alps (Fatichi et al., 2014(Fatichi et al., , 2015, and the Himalayas (Ragettli et al., 2013(Ragettli et al., , 2015; can be run at different spatial and time steps (typically hourly or daily); and is well suited for long-term simulations (Ragettli et al., 2016). 45 TOPKAPI-ETH is forced with time series of precipitation, air temperature, and cloud transmissivity of solar radiation. The model simulates snowfall at a given grid cell when precipitation occurs and air temperature is below a threshold parameter. If air temperature is above that threshold, precipita-50 tion is considered rain. When snow accumulation exceeds a slope-dependent threshold of a given grid cell (snow holding depth, S hd ), excess snow is moved to a lower grid cell based on the SnowSlide gravitational transport model (Bernhardt and Schulz, 2010): where SGR C (m) and SGR a are empirical parameters and SLP is the slope of the grid cell. Snow and ice melt is calculated with the Enhanced Temperature Index (ETI) model (Pellicciotti et al., 2005), depending on the net solar radia-60 tion and near-surface air temperature: where M is melt (mm h −1 ), SRF is the shortwave radiation factor (mm m 2 h −1 W −1 ), S in is the incoming shortwave radiation (W m −2 ), α is surface albedo, TF is the temperature fac-65 tor (mm h −1 • C), T a is air temperature ( • C), and T T is the air temperature threshold parameter for the onset of melt ( • C). TOPKAPI-ETH internally converts the units of the ETI variables and parameters to a daily time step. TOPKAPI-ETH does not compute sublimation. To calculate ice melt under 70 supra-glacial debris we also use the ETI model but with reduced melt factors (see Sect. 4.1.3). Although TOPKAPI-ETH includes a melt module that accounts for debris thickness in the computation of sub-debris ice melt, we did not use it due to the lack of debris thickness information in the region 75 and the large uncertainties that are present in large-scale estimates of debris thickness (Rounce and McKinney, 2014;Schauwecker et al., 2015). As a result of our assumptions, we expect that some of the spatial patterns of glacier ablation induced by the spatial variability of supraglacial debris 80 thickness are not accurately represented in our simulations. Once snow accumulation and melt are integrated to calculate the annual glacier surface mass balance, TOPKAPI-ETH translates it to elevation changes at the end of each hydrological year (from April to March) by means of the h approach 85 (Huss et al., 2010). This is done by using the originally proposed, glacier-size-dependent parameters (cf. CE4 Fig. 3b in Huss et al., 2010). Negative annual mass balances can result in glacier area reductions, but no area increases due to positive mass balances are prescribed. Area changes are ap-90 plied at the end of March. While snow melt over a nonglacierized grid cell is added to the respective soil layers, snow and ice melt over glaciers are added to a conceptual water reservoir for each glacier, which releases its water by means of a linear reservoir equation (Jansson et al., 2003). 95 In non-glacierized grid cells, the model simulates subsurface water flow, evapotranspiration, and water routing (Ciarapica and Todini, 2002).

Model setup for the Maipo River basin
We set up an instance of the TOPKAPI-ETH model for the entire Maipo River basin at a spatial resolution of 1 km that does not include glaciers. Glaciers and their runoff contribution are accounted for separately in the next section, but their 5 ice melt contribution is included in this section for the calibration of the subsurface parameters. The objective of the 1 km resolution setup for the entire Maipo River basin is to simulate snowmelt and rain, which account for the largest runoff volumes in the basin, at a resolution that allows for 10 multiple model runs and the automatic calibration of the subsurface flow parameters. The model was run continuously from 1955 to 2016 at a daily time step.
We spatially distribute daily precipitation over the basin using monthly mean maps derived from the CR2 pre- Embalse El Yeso using basin-scale daily temperature lapse rates (see Sect. 3.3). Periods with no direct information of daily mean air temperature at Embalse El Yeso are filled using correlation with records of daily extreme temperatures at the same station (mainly the period [1962][1963][1964][1965][1966][1967][1968][1969][1970][1971][1972][1973][1974][1975][1976][1977] or at Quinta 25 Normal station . The calibration of the Maipo River basin model was performed for the period April 2003 to March 2016, and consists of two steps: (i) the snow parameters are varied in order to fit SCA and SWE aggregated at the scale of the entire basin 30 from the MODIS and CAMELS-CL datasets (Sect. 3.4), and (ii) the parameters controlling subsurface fluxes are varied in order to fit monthly mean streamflow records at Maipo en El Manzano. While parameters in step (i) are manually calibrated and largely correspond to default values from previ-35 ous studies using TOPKAPI-ETH, parameters in step (ii) are automatically calibrated by minimizing three different evaluation metrics (Nash-Sutcliffe, NS; root-mean-square error, RMSE; and mean bias, BIAS). As the SWE reconstruction data are available starting in 1984, we use the period During the calibration procedure, we find that the use of the precipitation amounts derived from the CR2 product leads to an underestimation of SCA and SWE over the basin area, and streamflow at the basin outlet. This under-45 estimation of precipitation by the CR2 product was already identified by Alvarez-Garretón et al. (2018) when analyzing runoff ratios across Chile and attributed to a limitation of satellite-derived precipitation estimates over high-elevation areas. Similar results have been found in this region using a 50 regional climate model driven by ERA-Interim (Bozkurt et al., 2019) and the MERRA reanalysis (Cortés et al., 2016). Although the CR2 precipitation product corrects the ERA-Interim values by comparing them with ground data, these data are available only below 3000 m a.s.l. in this region, and 55 have not been corrected for gauge undercatch (DGA, 2017), which can also contribute to the underestimation of precipitation at the highest elevations (Rasmussen et al., 2012). We obtain a precipitation correction factor by manually fitting the observed and modeled curves of SCA and SWE and at 60 the same time closing the water balance of the basin. We obtain a value of +50 %. This correction generates precipitation amounts on the order of 3 to 4 times larger than that registered on low-lying areas. This value is larger than those estimated by previous studies on the west side of the semi-65 arid Andes (Falvey and Garreaud, 2007;Viale et al., 2011;Cortés et al., 2016), which estimated that the orographic effect results in a precipitation enhancement on the order of 2 to 3. The spread of precipitation amounts estimates over the semiarid Andes (and in general over mountain areas) is in 70 fact large, and previous hydrological studies have performed different types of corrections to close the water balance at the basin scale (Vicuña et al., 2011;Ragettli and Pellicciotti, 2012;Burger et al., 2019).
An additional aspect of model simplifications identified 75 during the model calibration is that air temperature over areas above 5000 m a.s.l. (about 5 % of the basin) is most of the time lower than the air temperature threshold parameter for melt onset, generating large snow accumulation that is not seen in the SWE reconstruction product. As snow on these 80 high-elevation areas is in reality removed by wind transport and sublimation, we reset the SWE in the model to zero at the beginning of each hydrological year. Although this implies that the model is not strictly mass-conserving, we verify that the discarded snow is, on average, 34 mm yr −1 over 85 the entire basin (or 688 mm yr −1 = 1.9 mm d −1 over the areas above 5000 m a.s.l.), which is similar to the estimates of sublimation amounts for this region (Corripio, 2003;Ayala et al., 2017a, b) and is on the order of the model uncertainties (49.9 mm w.e. in Fig. 2a). As elevation decreases south, 90 the discarded snow varies from about 121 mm w.e. yr −1 over the Colorado sub-catchment to about 10 mm w.e. yr −1 over Upper Maipo. Figure 2 shows the results of the model calibration for daily time series of SWE (Fig. 2a); monthly time series 95 of streamflow (Fig. 2b); and seasonal variations in SCA (Fig. 2c), SWE (Fig. 2c), and streamflow (Fig. 2d). The final calibrated snow parameters for this setup are shown in Table 1, whereas values for the subsurface flux parameters are shown in the Supplement (Table S1). The quality metrics 100 for snow and streamflow variables show very good results in both the calibration and validation periods. Extreme values are well captured, except for the humid winter of 1988, in which the model underestimates snow accumulation and streamflow.

Model setup for individual glaciers
In addition to the basin-scale model, we set up an instance of TOPKAPI-ETH for each one of the glaciers larger than 1 km 2 in the catchment (about 59 glaciers). These instances have a spatial resolution of 100 m, which is more adequate to 5 simulate the processes governing glacier mass balance. The domain of these models runs correspond approximatively to the smallest catchment that contains the 1955 glacier extent of each glacier. The models are run at a daily time step starting in the year 1955 and are then restarted in 2000 using the 10 topographic and geographic information from that year. The models are forced using daily precipitation at the location of the centroid of each glacier, linearly interpolated from basinaveraged precipitation (including the 50 % precipitation correction) (Alvarez-Garretón et al., 2018), and assumed to be 15 uniform over each corresponding domain. Air temperature is extrapolated from the Embalse El Yeso meteorological sta-tion using a constant air temperature gradient equal to the environmental lapse rate (−6.5 • C km −1 ). For the study period in which no CR2 precipitation products are available, 20 the Quinta Normal and Embalse El Yeso stations are used.
We choose a set of model parameters typically used in the literature for this region (Ragettli and Pellicciotti, 2012;Ayala et al., 2016;Burger et al., 2019) for all individual glacier models and keep parameter calibration at a minimum 25 level. For each glacier, we vary only the ETI model parameters within ranges suggested in the literature (Finger et al., 2011;Ragettli and Pellicciotti, 2012;Ayala et al., 2017b) to fit the glacier-wide mass balance as derived from the geodetic mass balances. Glacier-wide mass balance is considered 30 fitted when the difference between the simulated and observed balance is smaller than a certain threshold. We find that choosing a threshold equal to half of the uncertainty in the geodetic mass balance allows for reliable simulations while keeping an acceptable computation time. The uncer-35 Although we set up a TOPKAPI-ETH model for all glaciers with an area above 1 km 2 in 2000 (equivalent to 59 glaciers), we find that staying within the selected ranges for the ETI parameters only allows us to fit the geodetic mass balances in 26 cases. Among the discarded glaciers, about 5 half of them are smaller than 3 km 2 , and the rest correspond to those lying on the slopes of the Tupugatito Volcano and San José volcanic complex (Volcán Tupungatito, Azufre, and Marmolejo glaciers). We suspect that this is an expression of the fact that some of the processes not included in TOPKAPI-10 ETH (namely permafrost, sublimation, snow dynamics, or geothermal fluxes) may play a role governing the mass balance of these glaciers. However, it might also be related to local deficiencies in the spatial distribution of air temperature and precipitation. No rock glaciers are included in this 15 subset of glaciers. The location and main properties of the 26 modeled glaciers in comparison with those of the total sample are shown in Fig. 3 and Table 2, respectively. The simulated glaciers are spread over the entire basin, and their mean elevations are in the middle range of the total sample. 20 Glaciers smaller than 1 km 2 , from which 85 % correspond to rock glaciers or glacierets, are less well represented by the sample of 26 glaciers. The sample of 26 glaciers is mostly oriented towards the south (aspect > 90 • ) and does not include the steepest glaciers. In Fig. 3, we also highlight the 25 areas with the discarded large glaciers on Tupungatito Volcano and the San José volcanic complex.
The results of the calibration of the TOPKAPI-ETH models for the 26 modeled glaciers are shown in Fig. 4a (  and Fig. 4b (2000-2013). The calibration results are 30 very good for both periods with area-weighted RMSEs of 1 and 0.2 m w.e. for the 1955-2000 and 2000-2013 periods, respectively. These errors are well within the uncertainty bounds of the geodetic mass balance. Figure 4c shows the resulting cumulative glacier mass balance for all simulated 35 glaciers, their area-weighted average, and the comparison with the glaciological mass balance measured on the Echaurren Norte Glacier since 1975. The fastest declining line of the sample corresponds to the Olivares Alfa Glacier, which has been previously identified as one of the glaciers with the 40 largest retreating rates in the basin (Malmros et al., 2016). Interestingly, several of the glaciers show a positive or nearneutral mass balance over the entire period, which might be an indication that these glaciers have already retreated close to a new equilibrium. However, this is not the general trend 45 in the basin (as shown by the average values in Fig. 4), and it is limited to some specific cases where glaciers have retreated to elevations above the basin-average ELA, or have been covered by thick debris. 50 We extrapolate the mass balance of the 26 modeled glaciers to the entire basin based on the methodology described by Huss (2012). In that work, a set of in situ glacier mass bal- Figure 3. Location of the 26 glaciers modeled with TOPKAPI-ETH. We highlight the volcanic areas in which some large glaciers were discarded from the modeled sample. We include the name of the main glaciers in this sample.

Extrapolation
ance measurements for Switzerland were used to calculate the mass balance of all glaciers in the European Alps. Here, 55 we calculate the annual surface mass balance B (m w.e.) of glacier g in year y with the following equation: whereB(gp) is the average annual mass balance in the study period p and B(sy) is the glacier annual mass bal-60 ance anomaly in the sub-catchment s, where glacier g is located. WhileB(gp) is extracted from the geodetic mass balance, B(sy) is derived from the TOPKAPI-ETH simulations. Equation (4) where g * ,s is the subset of modeled glaciers ( * ) in sub-70 catchment s. The time series of annual mass balance B (g, y) are then used to estimate the volume changes of each glacier throughout the study periods. Glacier areas (S) are updated due to negative changes in glacier volume (V ) (we do not prescribe 75 increases in glacier area due to positive annual mass balance)  by means of the volume-area scaling formula: where γ and c are the scaling parameters. In line with recommendations of the volume-area scaling theory (Bahr et al., 2015), the parameter γ is kept constant in all periods at 5 a value of 1.357, and we let c vary in order to fit the total glacier volume in the basin in years 1955 and 2000 (calculated in Sect. 3.2). Parameter c is calculated as 28.1 for 2000 (this value is also used afterwards), but a value of 21.1 is the one that fits best to our estimates of ice thickness in 1955. In 10 between these 2 years we use a linear interpolation of c.
In the calculation of area and volume evolution, we account for the uncertainties in the annual mass balance, inventoried glacier areas in 1955 and 2000, and the parameter c, by disturbing each variable with a random variation. These ran-15 dom variations are 1000 realizations of three normal probability distributions of mean 0 and standard deviations equivalent to the typical errors of each variable. From the uncertainties in the geodetic mass balance, we estimate a typical error in the annual mass balance of 0.08 m w.e. yr −1 for the pe- estimates, and results in a value of 4.1. The uncertainty in parameter c should indirectly account for the different boundary conditions (such as basal sliding or surface geometry) that are found at each glacier (Bahr et al., 2015).
Glacier runoff, including all its components (i.e., ice melt, snow melt and rain), is extrapolated directly from the TOPKAPI-ETH results for the 26 modeled glaciers to the rest of the glacierized areas and does not depend on the extrapolated time series of glacier mass balance. At a particular year, the uncertainty in glacier runoff is estimated as a frac-10 tion of the same variable. That fraction is the same as that between glacier volume and its uncertainty in that year. As in Huss and Hock (2018), we compute glacier runoff as the water originating from the initially glacierized area (1955 in our case), i.e., independent of the glacier area in a particular 15 year. This allows for the evaluation of changes in total headwater runoff due to glacier retreat. However, in our study we also evaluate specific variations in the ice melt component. Throughout this paper, glacier runoff and its components are presented as normalized by the area of the entire Maipo River 20 basin.

Committed ice loss estimates
We estimate the committed glacier ice loss caused by the temperature increase in the last decades by conducting a set of 10 additional TOPKAPI-ETH simulations, and by extrap- 25 olating them using the same analysis as described in the previous section. The additional simulations are run under different synthetic climate scenarios in which the climate of the last 2 decades is stochastically repeated for a 100-year period. The meteorological inputs are built by repeating 1-year-30 long blocks of the input variables (precipitation, temperature, and cloud transmissivity) corresponding to a randomly selected year between 1993 and 2016 (23 years). We select this period because air temperature was relatively stable in the basin and precipitation showed the characteristic interannual 35 variability of this region. While the anomaly term B(sy) is calculated in the same way as for the period 1955-2016 (i.e., from the TOPKAPI-ETH simulations), as no geodetic mass balances are available for the synthetic scenarios, we calculateB (g, p) using 40 two different approximations depending on glacier size. For glaciers that are larger than the size of the smallest modeled glacier (1.1 km 2 ), we use a multiple linear regression of the mass balance of the modeled glaciers in each scenario with their topographic parameters in year 2000: 45B (g, p) = a 1 · x 1 + . . . + a n · x n , where a i are calibrated coefficients and x n are topographic parameters. In average for the 10 synthetic scenarios, the best results are given by glacier area, median glacier elevation, percentage of debris cover, mean sky view factor, and mean 50 aspect. Together, these five variables explain 52 % of the total variance, which is in the range of the original application of this methodology (obtained 35 % using three variables and 51 % using six Huss, 2012). Results of this procedure are summarized in Table 3. For glaciers smaller than 1.1 km 2 , 55 we use the average mass balance of modeled glaciers in the corresponding sub-catchment. As in the 1955-2016 period, rock glaciers are assumed to have a balanced mass budget.
Once the time series of mass balance for the 10 synthetic scenarios are calculated, we compute area and volume evolution 60 of each glacier and their associated uncertainties, using the same methodology as for the 1955-2016 period.

75
In Fig. 6, we present the temporal variability of precipitation (Fig. 6a), air temperature (Fig. 6a), the equilibrium line altitude (ELA) (b), and cumulative mass balance ( Fig. 6c and d) in the Maipo River basin since 1955. While the large interannual variability of the basin's mean precipitation (Fig. 6a, 80 blue bars) directly relates to the El Niño-Southern Oscillation (ENSO) phenomenon, a 3-year moving average of this variable exposes a sequence of dry (e.g., 1967-1969, 2010-2016) and wet periods (e.g., 1978-1987 and 2000-2008). This sequence has been related to other climatic in-85 dices, such as the Pacific Decadal Oscillation (PDO) or the Interdecadal Pacific Oscillation (IPO) (Boisier et al., 2016;González-Reyes et al., 2017). From 2010 on, precipitation has decreased due to a severe drought across Chile (Garreaud et al., 2017). Air temperature over the basin shows a 90 sustained increase in the long term but with relatively stable values since the mid-1990s. Since the 1960s, air temperature has increased in about 2 • C. Figure 6b shows the annual and decadal variability of the ELA of the 26 modeled glaciers. The ELA is calculated as the average elevation of 95 all grid cells with an annual mass balance of ±10 cm, and the estimated range (in light red) corresponds to the standard deviation. Since the 1960s, the elevation of the ELA has increased by 239 m (or 39 m per decade). These estimates of the ELA change are larger than those calculated by Carrasco 100 et al. (2005), who estimated an increase in the elevation of   Table S1.
the 0 • C isotherm of about 160 m for central Chile in the period 1975-2001. Figure 6c and d integrate the results of TOPKAPI-ETH and the extrapolation procedure. Figure 6c presents the cumulative surface mass balance of glaciers in the Maipo River 5 basin since 1955, including the 26 glaciers modeled with TOPKAPI-ETH, and the remaining glaciers in the basin for which extrapolation was used. The cumulative mass balance shows a decreasing trend interrupted by short periods of positive or near-neutral mass balance, with a more negative fi-10 nal value for the 26 modeled glaciers than for all glaciers in the basin. The more negative value for modeled glaciers might be caused by their larger area in comparison to the rest of the glaciers, as large glaciers have shrunk more extensively (Malmros et al., 2016). For comparison with the longterm glacier mass balance reference in the region, we include the direct measurements on Echaurren Norte Glacier, which presents a more negative trend, most likely due to its low elevation (3650 to 3900 m a.s.l.). In Fig. 6d, we present the surface mass balance of glaciers in each sub-catchment, where 20 relatively large differences can be seen. In general, glaciers in southern catchments show more positive mass balance than those in northern catchments. This can be explained by larger precipitation amounts and a higher proportion of both debriscovered and rock glaciers. Most notably, glaciers in Olivares 25 show the most negative mass balance throughout the study period, whereas those in Volcán present a positive mass balance until the mid-2000s. However, after the start of the current drought in 2010, negative glacier mass balances dominate across the entire Maipo River basin. The information in-30 cluded in Fig. 6d is summarized in Table 4, which shows the simulated glacier mass balance in each sub-catchment for the 1955-2016 period. For comparison, we include the mean elevation, mean latitude, and the 1955 glacierized area of each sub-catchment.

35
In Fig. 7 we show the variations in glacier runoff and its components (ice melt, snow melt and rain) in the initially glacierized areas over the period 1955-2016. While the annual and summer interannual variability is presented in Fig. 7a and b, Fig. 7c presents the average seasonal curve 40 and the percentage of each contribution. The summer period is chosen as January to March. Glacier runoff was 177 ± 25 mm yr −1 over the entire period and shows a sequence of three decreasing maxima (1968-1969, mid-1980s,    (c) cumulative glacier mass balance for the modeled glaciers (simulated with TOPKAPI-ETH), the entire basin (extrapolation) and its associated uncertainty, and the measurements on Echaurren Norte Glacier; and (d) cumulative glacier mass balance for each sub-catchment. In (b) the difference between the ELA in the last 10 years (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) and the first 10 years (1955)(1956)(1957)(1958)(1959)(1960)(1961)(1962)(1963)(1964)(1965) of the study periods is indicated, as is the equivalent ELA increase rate. The shadowed area in (b) shows the standard deviation of the elevation of grid cells with a mass balance between −0.1 and 0.1 m w.e.
ical year in record), and it averages 158±27 mm yr −1 during the current drought (2010)(2011)(2012)(2013)(2014)(2015)(2016). Figure 7a shows that the interannual variability of ice melt is very large (with a coefficient of variation of 0.57) and that its share in total glacier runoff can vary from less than 10 % (as in 1982-1983, 1997-5 1998, and 2002-2003) to more than 90 % (as in 1968-1969). Except for 1968-1969, snow melt in these areas is consistently the largest runoff contributor at the annual scale, but the contribution during summer is very variable. In Fig. 7c, we show the summary of runoff contributions at the annual 10 scale. Runoff contribution is dominated by snowmelt (60 %), with ice melt representing 37 % of the annual total. Rain represents about 3 %, but these amounts have increased since 1955 ( Fig. 7a and b).
In Fig. 8, we quantify the role that glacier runoff has 15 played in the entire Maipo River basin over the study period. At the annual scale, glaciers provide 16 ± 7 % of the total runoff, but this contribution can increase up to 59 ± 23 % drought is close to the average value over the entire study period.  1968-1969, 1990-1991, and 2011-2012). We estimate that if glaciers reached equilibrium with the current climate, the peaks would be 50 considerably lower than those in the 1955-2016 period. In this equilibrium situation, the peaks are close to 40 % of the largest ice melt runoff contribution in the past (1968)(1969). The information presented in Figs. 9c and 10 is also summarized in Table 4. In that table, we present the average glacier 55 runoff contribution per sub-catchment for the 1955-2016 period and the last 20 years of the committed ice loss scenarios.

Glacier changes
Our results indicate that the total glacier volume in the Maipo 60 River basin decreased by about one-fifth in the period 1955-2016. The cumulative glacier mass balance in the Maipo River basin shows variations that are similar to those registered on Echaurren Norte Glacier. These variations consist of a general decreasing trend, concurrent with an increase 65 in the ELA (Carrasco et al., 2005(Carrasco et al., , 2008, which has been interrupted by periods of slightly positive or neutral mass balance. Since the mid-1980s, there has been a strong mass loss, interrupted only by a positive period in the beginning of the 2000s. In fact, from 2000 on, we observe a 10-year 70    and the committed ice loss scenarios assuming a constant climate. In (a) and (b) we use results from the glacier inventories and the combination of ice thickness estimates and geodetic mass balances, respectively, as observations of glacier area and volume. Glacier runoff in (c) is computed for the Maipo River basin (i.e., runoff units are normalized by the basin area). While the uncertainty bars of the observations are shown in green, those of the simulations are shown in blue for the 1955-2016 period and in red for the committed ice loss scenarios. For visual purposes, we present the committed ice loss scenarios using the period 2000-2100 in the x axis. period with positive or nearly neutral mass balance. This was also described by glaciological observations (Masiokas et al., 2016) and geodetic mass balances (Braun et al., 2019;Dussaillant et al., 2019). Following this period, strongly negative mass balances have been observed (Masiokas et al., 2016;Figure 10. Variations in ice melt in the Maipo River basin for the past period (1955-2016, in blue) and each one of the 10 committed ice loss scenarios (light red). The peaks of ice melt over the past period and those at the final decade of the committed ice loss scenarios are highlighted in red. On the right axis, we set the ice melt estimated for the severe drought of 1968-1969 as 100 %. Ice melt is computed for the Maipo River basin (i.e., runoff units are normalized by the basin area). For visual purposes, we present the committed ice loss scenarios using the period 2000-2100 in the x axis.
central Chile, unprecedented in extension and duration (Garreaud et al., 2017). For the period before 1975, when the mass balance measurements on Echaurren Norte Glacier started, we compare our results to the reconstruction obtained by Ma-10 siokas et al. (2016). The latter estimated a strongly negative mass balance in the period 1955-1975, while we obtain a nearly neutral mass balance between 1955 and 1968 and a more negative balance from 1968 to 1975. These differences might either correspond to differences between the basin-15 averaged mass balance and that on Echaurren Norte Glacier or be a consequence of the different methodologies between our study and that of Masiokas et al. (2016). While Masiokas et al. (2016) relied on the correlation between hydrometeorological records and the measured surface mass balance on 20 Echaurren Norte, our model is calibrated to the geodetic mass balances.
The different trends of glacier mass balance in the subcatchments of the Maipo River basin are an expression of the diverse climatic and morphological characteristics that dom-25 inate across the basin. For example, the positive and near-neutral glacier mass balances in Volcán and Upper Maipo might be related to higher precipitation towards the south or that several glaciers have retreated close to a new equilibrium state. Within the Olivares sub-catchment, the geodetic mass balance is in line with the large areal changes found by Malmros et al. (2016). This might be explained by a strong imbalance of the large glaciers in that catchment and/or by the impacts of nearby mining activities (Los Bronces and Andina mines), especially dust deposition on Olivares Alfa Glacier and its neighbor glaciers. However, more specific 10 studies addressing albedo changes are necessary to obtain more conclusive results.
Our estimates of committed ice loss show that glaciers will continue to shrink if the climate remains stable, with an estimated committed ice loss of 20 % relative to the vol-15 ume in the year 2000 (30 % relative to 1955). We stress that these estimates are an indication of the glacier changes that past climate will produce in any case. Future projections under emission scenarios will certainly show more dramatic reductions of glacier area and volume. In the context 20 of future projections, we highlight that most projections for glacier changes in the Central Andes (Marzeion et al., 2012;Radić et al., 2014;Huss and Hock, 2015) are included in the macro-region of the "southern Andes", which also contains the Patagonian Ice Fields. Future projections of glacier 25 changes in the region are thus strongly influenced by these large ice masses with their peculiar climate and physical processes (e.g., calving) that are not representative of the small mountain glaciers along the semiarid Andes (Mernild et al., 2015). 30

Glacier runoff
Despite the reduction of glacier volume in the period 1955-2016, our estimates of glacier runoff do not show (Figs. 8 and 9c) the typical increasing or decreasing phases of peak water observed or projected for other catchments across the 35 world (Baraer et al., 2012;Farinotti et al., 2012). This result is similar to that obtained by Casassa et al. (2009), who did not find significant trends in an analysis of Maipo streamflow records. Peaks in glacier runoff have reduced their magnitude over the last decades, due to a combination of a decrease in 40 precipitation (Boisier et al., 2016) and the reduction of ice volume, but it is difficult to identify if there was an increasing phase of glacier runoff in the period 1955-2016.
We suggest that the strong interannual and inter-decadal climatic variability observed in the semiarid Andes (Monte-45 cinos and Aceituno, 2003; Masiokas et al., 2006;Falvey and Garreaud, 2007) is also transferred to the glacier runoff time series, modifying or masking the typical trends associated with glacier retreat. Once an extended time period is considered (in this case a committed ice loss scenario, Fig. 9c), 50 peak water emerges more clearly. Huss and Hock (2018) estimated that glacier runoff in the Rapel River basin (south of the Maipo River basin) experiences peak water in the current decade (2010-2020), but their analyses also show a strong interannual glacier runoff variability from 1980 to 2010. This 55 makes peak water evident only when compared to the future projections under different emission scenarios, in which the glacier runoff is considerably lower than present levels.
6.3 Uncertainties in the modeling of glacier changes in data-scarce regions 60 We identify four main sources of errors and uncertainties in our study: (i) the glaciological datasets, i.e., the geodetic mass balance, glacier outlines, ice thickness, and debris cover areas; (ii) the spatial distribution of meteorological inputs; (iii) modeling limitations in TOPKAPI-ETH; and (iv) limita-65 tions of the extrapolation methodology.
In general, the uncertainties of the elevation changes and glacier properties are well quantified and explicitly stated in the confidence bounds of the geodetic mass balance (Fig. 9). However, a few properties were not explicitly quantified, 70 such as the ice content in rock glaciers, which is in fact a key problem in the semiarid Andes (Schaffer et al., 2019). Due to this, there might be an overestimation in the ice content due to the presence of rock glaciers. In the future, more geophysical measurements to acquire information on ice content in 75 rock glaciers of the semiarid Andes (e.g., Croce and Milana, 2002) could improve the estimates of runoff generation from these landforms.
The accuracy in the spatial distribution of meteorological inputs is particularly difficult to evaluate, and it likely cor-80 responds to a major source of uncertainty (especially precipitation). This is because of the relatively sparse network of meteorological stations installed in the basin, the difficulties of atmospheric models to represent precipitation processes over the Andes (Bozkurt et al., 2019), and the un-85 derestimation of satellite-based precipitation products over high-elevation areas (Alvarez-Garretón et al., 2018). However, the indirect evaluation of precipitation amounts through snow cover products and the basin's water balance increase the confidence in the results of this study. An additional sim-90 plification in the meteorological distribution is the extrapolation of air temperature from one single station. Nevertheless, we are confident that air temperature variability is well constrained over the catchment because it usually correlates well over long distances, daily lapse rates are derived from 95 the basin-wide CR2 temperature dataset, and the timing of snow disappearance is well simulated by TOPKAPI-ETH.
Although our study has benefited from a series of new meteorological and glaciological datasets presented for the southern Andes in recent years (Cortés and Margulis, 2017;100 Alvarez-Garretón et al., 2018;Farías-Barahona et al., 2020), the lack of field data in the Maipo River basin is something that needs to be taken into account in glacio-hydrological modeling studies in the region, particularly at remote highelevation sites. In this study, we alleviate the difficulties 105 posed by the lack of basin-wide field data, and its impact on the TOPKAPI-ETH results, by deriving most of the model parameters from data collected in previous field campaigns in this region starting in 2008 (Pellicciotti et al., 2008;Ragettli and Pellicciotti, 2012;Ayala et al., 2016). These previous studies have also shown that many of the parameters required by the model are fairly stable, in the sense that they can be extrapolated from one glacierized area to another with a reasonable degree of confidence (Ragettli et al., 2014;Ayala et al., 2017b;Burger et al., 2019). In addition, the representation of processes driving the mass balance at some specific sites re-10 quires more fundamental work, and additional parameterizations or more physically based representations are required. Such sites correspond mainly to sublimation-dominated sites above 5500 m a.s.l., where we had to correct our simulations of snow depth, and the temperature index modeling is inaccu-15 rate (Ayala et al., 2017a, b), debris-covered areas with complex distributions of debris thickness , or steep glacierized slopes such as the volcanoes in the Maipo River basin.
As the average rates of annual mass balance in the 1955-20 2016 period are calculated from the geodetic mass balance, the long-term glacier changes derived from the extrapolation methodology should be well simulated, but the year-to-year variations in mass balance depend on the representativity of the modeled glaciers. As discussed by Huss (2012), low rep-25 resentativity of the reference glaciers could lead to large errors in the mass balance of individual glaciers and years, but these errors should be lower at the mountain range scale and over long time periods. In the committed ice loss scenarios, the uncertainty of glacier changes is higher because it also 30 relies on the multiple linear regression analysis. In any case, we explicitly accounted for the uncertainties in the extrapolated mass balance (at least partly) by means of the random perturbation described in Sect. 4.2.

35
We have reconstructed the changes that glaciers in the Maipo River basin experienced over the last 6 decades, with a focus on glacier runoff and the impacts of its long-term variations on the basin's hydrology. These results add a missing piece to the current hydro-climatological knowledge of the 40 semiarid Andes and can be useful for water managers and stakeholders to develop adaptation or mitigation strategies. Although some uncertainties still remain, our results successfully take into account a number of independent datasets, including snow cover area variations, snow water equivalent 45 reconstructions, streamflow records, glacier inventories, and geodetic mass balances.
Our main conclusions are as follows.
a. Over the period 1955-2016, the total glacier volume in the Maipo River basin has decreased by 20±14 % (from 50 18.6 ± 4.5 to 14.9 ± 2.9 km 3 ). In agreement with other studies, our results show that the cumulative glacier mass balance over the study period had a general decreasing trend, interrupted by short periods of positive or near-neutral mass balance. This might be an in-55 dication that some glaciers temporarily retreated to a new equilibrium state. Strongly negative mass balances have dominated since the start of the current drought in 2010. Despite the general trend, there are important differences between the glacier mass balances of 60 the sub-catchments, with the southern sub-catchments (Volcán and Upper Maipo) showing positive or nearneutral mass balances until 2000 and the Olivares subcatchment showing a strongly negative mass balance over the entire period. 65 b. The average glacier contribution to runoff in the Maipo River basin (i.e., the runoff contribution from liquid precipitation, snowmelt ,and ice melt from the areas that were glacierized in 1955) was 177 ± 25 mm yr −1 in the period 1955-2016. Instead of a clear peak wa-70 ter, we identify a decreasing sequence of runoff maxima that can be linked to both a decrease in precipitation since the 1980s and a reduction of ice melt. The exact occurrence of peak water will also depend on future changes (e.g., more precipitation or more ice 75 melt), which are not addressed in our article. Glacier runoff has decreased since the severe drought of 1968-1969, when glacier runoff peaked at 245 ± 62 mm yr −1 (49 % of the basin's total runoff). During the current drought, which started in 2010, the contribution has 80 been 158 ± 27 mm yr −1 (17 % of the total runoff).
c. If climate was to stabilize at the level of the past 2 decades, we estimate a committed glacier ice mass loss of 19 ± 38 %. This would cause glacier runoff to reduce by 22±30 % when compared to the 1955-2016 average 85 or by 39 ± 21 % when compared to [1968][1969]. Based on these numbers, we anticipate that the future capacity of the basin to mitigate severe droughts will be reduced.
Our results shed light on the glacier runoff evolution in the semiarid Andes and complement recent studies that as-90 sessed regional-scale glacier changes (Braun et al., 2019;Dussaillant et al., 2019). Some topics deserving further attention that should be addressed are the drivers behind the positive mass balance in the southern catchments (Volcán and Upper Maipo), the processes governing mass balance 95 on glaciers on active volcanoes, the possible anthropogenic impacts on glaciers in the Olivares sub-catchment, and the quantification of the hydrological role of rock glaciers. While our simulations of committed ice loss provide estimates of the minimum changes that glaciers will experience due to 100 past changes of the climate, future studies driven by climate model simulations and emission scenarios should provide more realistic projections for the future of the region's glaciers.
Data availability. Data used in this study are available upon request to the authors.
Author contributions. ÁA designed the study with contributions 5 from DF and DFB and performed the calculations with TOPKAPI-ETH and the extrapolation method. DFB calculated the geodetic mass balances. MH calculated the ice thicknesses for 2000. All coauthors provided scientific advice during the study. ÁA prepared the manuscript with contributions from all co-authors.
Eduardo Muñoz for his help in the calibration of TOPKAPI-ETH. DFB acknowledges CONICYT Becas Chile. The work benefited from the financial support of a seed grant within ETH Zurich's Research for Development (R4D) funding scheme and was conducted within the research collaboration "Glacier runoff contributions 20 in the Maipo River catchment" between the Fundación para la Transferencia Tecnológica (UNTEC), Santiago, Chile, and ETH Zurich, Switzerland.
Financial support. This research has been supported by the NAME OF FUNDER (grant no. GRANT AGREEMENT NO). TS10

25
Review statement. This paper was edited by Elizabeth Bagshaw and reviewed by Francisca Bown and one anonymous referee.
Remarks from the language copy-editor CE1 Please confirm the change.

CE2
Does this need to be defined?

CE3
Please note that, following the dominant usage, this paper has been edited according to the standards of American English.

CE4
Please check the use of "cf." (i.e. "compare") throughout. Do you mean "see" or "compare"? If you mean "see", then it should be changed to "see" throughout.
Remarks from the typesetter TS1 Please provide city and country.

TS3
Please provide date of last access.

TS4
Please provide date of last access.

TS5
Please provide date of last access.

TS6
Please provide date of last access.

TS7
Please provide year.

TS8
Please confirm change.

TS9
Please confirm change.

TS10
Please note that there is funding information given in the acknowledgements but you have not indicated any funding upon manuscript registration. Therefore, we were not able to complete the financial support statement. Please provide the missing information and double-check your acknowledgements to see whether repeated information can be removed from the acknowledgements. Thanks.

TS11
Please provide volume with article number or page range.

TS12
Please provide date of last access.

TS13
Please provide article number or page range.

TS14
Please provide journal.

TS15
Please provide article number or page range.

TS16
Please provide update if available.

TS17
Please provide an update if available.

TS18
Please provide date of last access.

TS19
Please provide date of last access.

TS20
Please provide article number or page range.