# Two decades of glacier mass loss along the Andes

I. Dussaillant, E. Berthier, F. Brun, M. Masiokas, R. Hugonnet, V. Favier, A.

Rabatel, P. Pitte, L. Ruiz

► To cite this version:

I. Dussaillant, E. Berthier, F. Brun, M. Masiokas, R. Hugonnet, et al. Two decades of glacier mass loss along the Andes. Nature Geoscience, 2019, 12 (10), pp.802-808. 10.1038/s41561-019-0432-5. hal-02365259

HAL Id: hal-02365259

https://hal.science/hal-02365259v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Two decades of glacier mass loss along the Andes

I. Dussaillant<fnref n="1" />,²*, E. Berthier<fnref n="1" />, F. Brun<fnref n="1" />,², M. Masiokas<fnref n="3" />, R. Hugonnet<fnref n="1" />, V. Favier<fnref n="1" />,², A. Rabatel<fnref n="1" />,²,
P. Pitte<fnref n="3" /> and L. Ruiz<fnref n="3" />

Andean glaciers are among the fastest shrinking and largest contributors to sea level rise on Earth. They also represent crucial water resources in many tropical and semi-arid mountain catchments. Yet the magnitude of the recent ice loss is still debated. Here we present Andean glacier mass changes (from 10°N to 56°S) between 2000 and 2018 using time series of digital elevation models derived from ASTER stereo images. The total mass change over this period was −22.9 ± 5.9 Gt yr⁻<fnref n="1" /> (−0.72 ± 0.22 m w.e. yr⁻<fnref n="1" /> (m w.e., metres of water equivalent)), with the most negative mass balances in the Patagonian Andes (−0.78 ± 0.25 m w.e. yr⁻<fnref n="1" />) and the Tropical Andes (−0.42 ± 0.24 m w.e. yr⁻<fnref n="1" />), compared to relatively moderate losses (−0.28 ± 0.18 m w.e. yr⁻<fnref n="1" />) in the Dry Andes. Subperiod analysis (2000–2009 versus 2009–2018) revealed a steady mass loss in the tropics and south of 45°S. Conversely, a shift from a slightly positive to a strongly negative mass balance was measured between 26 and 45°S. In the latter region, the drastic glacier loss in recent years coincides with the extremely dry conditions since 2010 and partially helped to mitigate the negative hydrological impacts of this severe and sustained drought. These results provide a comprehensive, high-resolution and multidecadal data set of recent Andes-wide glacier mass changes that constitutes a relevant basis for the calibration and validation of hydrological and glaciological models intended to project future glacier changes and their hydrological impacts.

The Andes contain a wide variety of topographic and climatic conditions from ~10°N to 56°S that results in a large diversity of ice masses<fnref n="1" />,². It hosts the largest glacierized area in the Southern Hemisphere outside Antarctica. Glaciers cover a wide range of altitudes from the highest peaks and volcanoes higher than 6,000 m above sea level in the Tropical and Equatorial Andes to sea level in Patagonia and Tierra del Fuego. All the recent global estimates of glacier contribution to sea level rise (SLR) place the Andean glaciers among the highest contributors<fnref n="3" />⁻⁸. This region experienced the largest cumulative specific mass losses since 1961.<fnref n="1" /> Many Andean glaciers represent important water resources along an extensive portion of the Tropical and Dry Andes where glacier melt acts as a buffer during periods of drought¹⁰⁻¹². Together with Central Asia, the Andes is a region in which the snow and ice melt contribution can reach 50% or more of the total runoff in some basins¹³, and in which the largest reductions in glacier runoff are projected for the end of the twenty-first century<fnref n="1" />. These glaciers also constitute important tourist attractions that generate significant revenues for local and regional economies¹⁴.

Although the recent shrinkage of Andean glaciers is relatively well documented, the magnitude of the total mass loss remains debated. Glacier mass change estimates have been obtained by four main strategies: (1) extrapolation of the scarce glaciological (field) measurements¹⁵⁻¹⁸; (2) using geodetic estimates derived from digital elevation models (DEMs) for specific glacierized areas, such as the North Patagonian Icefield (NPI, for example, Rignot et al.¹⁹), South Patagonian Icefield (SPI, for example, Malz et al.²⁰) and Cordillera Darwin Icefield (for example, Melkonian et al.²¹) and, only recently, for the entire Andes²²; (3) using low-resolution remote sensing methods⁴,⁵,⁷ that provide regional estimates but do not resolve individual glaciers nor regions with a small concentration of glaciers (ICESat-based estimates are restricted to 2003–2009, and Gravity Recovery and Climate Experiment (GRACE)-based estimates need to account for leakage of the hydrological signal and the strong signal from glacio-isostatic adjustment²³,²⁴) and (4) using models

of different complexity⁶,²⁵,²⁶. Among others, uncertainties in these models arise from the lack of measurements to calibrate both the surface mass balance (SMB) (accumulation is particularly uncertain over the Patagonian icefields²⁷,²⁸) and frontal ablation (calving) when included in the glacier model.

Recently, a study²² provided a comprehensive, high-resolution assessment of recent changes in the South American glaciers by differencing radar-based DEMs. However, it remains unclear whether the penetration of the radar signal into snow and firn can be neglected for all the latitudes and altitudes across the Andes. Further, Braun et al.²² examined a single time interval (from 2000 to 2011/15) and thus did not resolve the decadal mass balance variability.

In this study, we provide an Andes-wide mass balance estimate that extends from 10°N to 56°S based on the 'Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) monitoring of ice towards extinction' geodetic method (hereafter called ASTERIX). This method (Methods and Brun et al.²⁹ give more details) exploits the extensive and freely available archive of optical stereo images acquired by the ASTER sensor, previously used in Alaska, the high mountains of Asia and the NPI²⁹⁻³¹. One of the strengths of the ASTERIX method is that it relies on optical imagery, which avoids the uncertainties associated with the signal penetration into snow and ice of radar-based DEMs³⁰. Compared to most of the earlier assessments, our mass loss estimates have a significantly higher resolution (down to individual glaciers), cover nearly two decades (2000–2018) and capture the shift in mass balance between two subperiods, 2000–2009 and 2009–2018. The results are based on the processing of 90 1° by 1° (latitude by longitude) tiles, that altogether contain 99.6% (31,505 km²) of the Andes glacierized area, for which the overall ASTERIX coverage with valid elevation change data is close to three-quarters (Table 1).

## Spatial variability of glacier changes

Results are reported as average mass balance rates over grid cells of 1° by 1° to help with their visualization and interpretation

<footnote n="1">¹LEGOS, Université de Toulouse, CNES, CNRS, IRD, UPS, Toulouse, France. ²Université Grenoble Alpes, CNRS, IRD, Grenoble INP, IGE, Grenoble, France.</footnote>

<footnote n="3">³IANIGLA, CCT-CONICET, Mendoza, Argentina. *e-mail: ines.dussaillant@gmail.com</footnote>
<table><caption>Table 1 | Region-wide mass balances for different subregions of the Andes from 2000 to 2018</caption><thead><tr><th rowspan="2">Region</th><th rowspan="2">Total glacier area (km²)</th><th rowspan="2">Glacier area measured (%)</th><th rowspan="2">Mean number of DEMs</th><th rowspan="2">Survey period (yyyy/mm)<sup>a</sup></th><th colspan="2">Mass balance</th></tr><tr><th>(m w.e. yr⁻¹)</th><th>(Gt yr⁻¹)</th></tr></thead><tbody><tr><td>Inner Tropics</td><td>174</td><td>80</td><td>11.5</td><td>2001/08-2016/11</td><td>-0.37 ± 0.25</td><td>-0.1 ± 0.0</td></tr><tr><td>Outer Tropics</td><td>1,967</td><td>94</td><td>17.3</td><td>2001/06-2017/03</td><td>-0.42 ± 0.23</td><td>-0.9 ± 0.5</td></tr><tr><td>Desert Andes</td><td>387</td><td>97</td><td>31.5</td><td>2000/11-2017/09</td><td>-0.12 ± 0.17</td><td>-0.1 ± 0.1</td></tr><tr><td>Central Andes</td><td>1,795</td><td>97</td><td>28.3</td><td>2001/08-2017/11</td><td>-0.31 ± 0.19</td><td>-0.6 ± 0.3</td></tr><tr><td>North Patagonia</td><td>1,753</td><td>96</td><td>18.2</td><td>2001/12-2017/07</td><td>-0.57 ± 0.22</td><td>-1.0 ± 0.4</td></tr><tr><td>South Patagonia<sup>b</sup></td><td>21,362</td><td>70</td><td>11.3</td><td>2000/10-2016/09</td><td>-0.86 ± 0.27</td><td>-18.3 ± 5.7</td></tr><tr><td>Fuegian Andes<sup>b</sup></td><td>4,053</td><td>62</td><td>6.9</td><td>2000/06-2016/03</td><td>-0.48 ± 0.27</td><td>-1.9 ± 1.1</td></tr><tr><td>Total Andes</td><td>31,631</td><td>74</td><td>12.8</td><td>2000/11-2016/10</td><td>-0.72 ± 0.22</td><td>-22.9 ± 5.9</td></tr><tr><td>Tropical Andes</td><td>2,271</td><td>93</td><td>16.9</td><td>2001/06-2017/03</td><td>-0.42 ± 0.24</td><td>-1.0 ± 0.5</td></tr><tr><td>Southern Andes</td><td>29,360</td><td>73</td><td>12.4</td><td>2000/11-2016/10</td><td>-0.75 ± 0.23</td><td>-21.9 ± 5.8</td></tr></tbody></table>

Supplementary Tables 1 and 2 give the subperiods. <sup>a</sup>For a specific subregion, the survey period varies because, for each pixel, the date of the first/last DEM can be different. The survey period is thus the interval from the mean first date to the mean last date of all the DEMs used to extract rate of elevation change (Supplementary Figs. 6 and 7 show the distribution of these start/end dates for each subregion). <sup>b</sup>SRTM DEM is included when computing the final trend.

(Fig. 1). Selected maps of the elevation change rates are shown in Supplementary Figs. 2–4 to illustrate the variety of patterns and the level of noise off-glacier. Region-wide thinning extends across the entire Andes for the period 2000–2018. A contrasted latitudinal pattern is also revealed with rapidly thinning small glaciers in the northernmost latitudes (10°N to 20°S) and over the largest glacierized area of the Patagonian Andes. An intermediate region from 26°S to 37°S exhibits comparatively smaller thinning rates (Fig. 1).

Mass balance estimates were computed for 17,832 glaciers out of 18,799 in the Randolph Glacier Inventory (RGI) v6.0. Region-wide mass balance rates were calculated for seven subregions, defined by latitude and climatology<sup>32,33</sup> (Fig. 2, Table 1), as the area-weighted sum of all glaciers belonging to the sub-region. The Andes-wide mass loss rate (area weighted sum of the subregional rates) totalled −22.9 ± 5.9 Gt yr⁻¹ (−0.72 ± 0.22 m w.e. yr⁻¹) from March 2000 to April 2018 (not considering the subaqueous losses). These dates correspond to the first and last available DEMs in our time series, but the exact temporal sampling varies from one subregion to another (Table 1).

Over the entire study period, the most negative mass balance rate was found in South Patagonia (−0.86 ± 0.27 m w.e. yr⁻¹), the subregion that contains the two largest icefields (NPI and SPI) and is characterized by large calving outlet glaciers. The glaciers of North Patagonia and the Fuegian Andes also presented negative mass balance rates of −0.57 ± 0.22 m w.e. yr⁻¹ and −0.48 ± 0.27 m w.e. yr⁻¹, respectively. The tropics, which correspond to rather small and spatially isolated glaciers in Ecuador and Colombia (Inner Tropics) and larger glaciers in Peru and Bolivia (Outer Tropics), also presented negative mass balance rates of −0.37 ± 0.25 and −0.42 ± 0.23 m w.e. yr⁻¹, respectively. Slightly less negative mass balance rates of −0.31 ± 0.19 m w.e. yr⁻¹ were detected in the Central Andes. The least negative value was measured in the Desert Andes, with −0.12 ± 0.17 m w.e. yr⁻¹. The ASTERIX coverage with valid glacier elevation change data exceeded 94% of the glacierized area in most subregions, with the exception of the Inner Tropics (80%), South Patagonia (70%) and the Fuegian Andes (62%) (Table 1). For South Patagonia this is explained by the homogeneity of the snow surface on the flat NPI and SPI plateaus, which prevents the correlation of the stereo pairs and leads to missing values at high elevations in many of the DEMs, and consequently to data gaps in the final maps of the elevation changes (Supplementary Fig. 4). Despite these data gaps, our total mass losses for the SPI and NPI are in good agreement with independent estimates based on an almost complete coverage<sup>20,34</sup>. In the Fuegian Andes, the relatively poorer ASTERIX coverage is

explained mostly by the high frequency of cloud cover. Nonetheless, the good agreement with radar-based estimates, together with the validation over Alaska icefields using airborne laser altimetry<sup>30</sup>, and a further sensitivity analysis to data gaps performed here (Supplementary Information) proves the relevance of our method, even in the featureless accumulation area of large icefields.

Consistently among the different subregions, the highest thinning rates were observed at the lowest elevations (Fig. 2). The strongest thinning rates were observed at the lowest elevations of the Patagonian subregions (up to 4 m yr⁻¹ in South Patagonia). A record-breaking thinning rate of 44 m yr⁻¹ was measured at the tongue of glacier HPS12 (SPI). In the glacier upper reaches, the rates of elevation change tend toward zero in most subregions. Exceptions to this are the glaciers in North Patagonia and the Inner Tropics, and to a lesser extent in the Outer Tropics, in which negative rates of elevation change were observed even at their highest elevation bands. This suggests that glaciers in these subregions are in disequilibrium with the current climate and in an advanced state of decline<sup>14,35–37</sup>.

## Decadal variability of glacier changes

To detect temporal variations in the glacier mass change, we computed the elevation change rates from March 2000 to March 2009 and from January 2009 to April 2018 for the different subregions (Fig. 1b,c and Supplementary Tables 1 and 2). For most subregions, the glacier mass changes do not differ significantly from one period to the other. Conversely, a noticeable difference was observed for the Desert Andes, Central Andes and North Patagonia. In these subregions, the rates of elevation changes during the first subperiod were close to zero or slightly positive. In contrast, strongly negative elevation change rates were detected during the second subperiod. The weighted-average mass balance rate for these three subregions decreased from 0.04 ± 0.18 m w.e. yr⁻¹ (2000–2009) to −0.63 ± 0.18 m w.e. yr⁻¹ (2009–2018). These estimates are robust, as these subregions present the largest percentage of valid data in the elevation change maps, with only 3% of data gaps for the 2000–2018 period (11% for 2000–2009 and 6% for 2009–2018). Moreover, the few available glaciological and local geodetic measurements support our findings of a shift from slightly positive to strongly negative mass balances after 2009<sup>38–40</sup>.

## Comparison with previous mass balance estimates

Previous studies provided mass balance estimates for the Southern Andes and low latitude glaciers (Fig. 3 and Supplementary Table 2),
<figure><img src="image_2.png" /></figure>

**Fig. 1** | Andes-wide glacier mass balance rates averaged by 1° latitude and 1° longitude tiles. a, March 2000 to April 2018. b, January 2000 to March 2009. c, January 2009 to April 2018. Histograms represent the glacierized area (grey bars) and the mean rate of elevation change (red bars) as a function of latitude. Empty values (Nan cells) correspond to tiles for which the mean coverage with rate-of-elevation-change data is less than 20%, that is, tiles situated in the Southern and Fuegian Andes, near to the Pacific coast, where cloud cover is persistent (Supplementary Fig. 1). Note that the colour scale intervals are irregular to help visualization.

but in most cases the comparison is not straightforward because of
differences in the time periods and/or different geographical cover-
age. The annual mass change estimates reported here for the entire
Andes (−22.9±5.9 Gt yr⁻¹) is slightly more negative than the recent
Andes-wide radar-based geodetic estimate from Braun et al.²²
(−18.4±0.6 Gt yr⁻¹). Both studies agree over the large Patagonian
icefields, where the coverage with radar DEMs is superior to ours.
However, our ASTERIX estimates are consistently more negative
than those from Braun et al.²² over all other Andean subregions
(Supplementary Table 5 and Supplementary Fig. 9). The differences
stem mostly from methodological aspects (that is, the choice of the
filtering technique applied to the elevation-change maps leads to
biases in the subregion mass losses) and the differential penetration
of the C-band (Shuttle Radar Topography Mission (SRTM)) and
X-band (TanDEM-X) radar signals used in Braun et al.²². A compre-
hensive comparison between our results and those of Braun et al.²²
is given in the Supplementary Information.

For low latitude Andean glaciers, the annual mass changes
reported here (−1.0±0.5 Gt yr⁻¹ for 2000–2018 and −1.1±0.5 Gt
yr⁻¹ for 2000–2009) are significantly less negative than those esti-
mated using the interpolation of glaciological and local geodetic
measurements for 2003–2009 (−4±1 Gt yr⁻¹) (ref. ⁴), 2002–2014
(−2.2±1.3 Gt yr⁻¹) (ref. ⁷) and 2000–2016 (−2.0±0.6 Gt yr⁻¹) (ref. ⁴
update of ref. ¹⁵). A better agreement was found with Wouters et al.⁴¹,
based on GRACE (−1.2±3.5 Gt yr⁻¹) for the period 2002–2016.

For the Southern Andes, our mass changes of -21.9±5.8 Gt-yr⁻¹ for 2000-2018 and -21.0±6.6 Gt yr⁻¹ for 2000-2009 are in good agreement with the Jacob et al.⁵ GRACE estimates for 2003-2010 (-23±9 Gt yr⁻¹) but are less negative than other GRACE estimates from Wouters et al.⁴¹ (-30±11 Gt yr⁻¹), Gardner et al.⁴ (-29±11 Gt yr⁻¹) and Reager et al.⁷ (-29±11 Gt yr⁻¹) and with the glaciological and geodetic interpolated estimates from Zemp et al.⁹ (-34±11 Gt yr⁻¹). Overall, these differences suggest a strong sensitivity of the region-wide values to very small number of glaciers with direct measurements in the Andes and their limited representativity. Mass changes estimated by GRACE are generally more negative in both RGI regions. This overestimation of the mass loss is possibly due to the strong influence of the glacio-isostatic adjustment in Patagonia²⁴,⁴²,⁴³. The magnitude of the hydrological corrections also varies greatly depending on the model⁴¹. A detailed comparison of our ASTERIX results with published local glaciological and geodetic observations is presented in the Supplementary Information.

The modelled SMB rate from Marzeion et al.⁶ is in good agree-
ment with our mass balance rate for low latitudes but is significantly
less negative for the Southern Andes. This is probably because
their model does not account for ice mass loss due to frontal
ablation, an important component of mass loss in the Southern
Andes. We separated calving and non-calving glaciers and com-
pared the modelled SMB rates from Marzeion et al.⁶ (updated to
<figure><img src="image_3.png" /><figcaption>**Fig. 2** | Delimitation of seven subregions and altitudinal distribution of glacier elevation change. **a**, Seven subregions defined according to latitude and climate. The Patagonian and Dry Andes correspond to the RGI Southern Andes region. **b**, Rate of elevation change of all the Andes subregions from 2000 to 2018 as a function of elevation. **c**, Same as **b**, but expressed as the normalized elevation defined as $(h-h_1) = (h_{99}-h_1)$, where $h$ is the elevation and $h_1$ and $h_{99}$ are the elevations of the 1st and 99th percentiles of the glacierized area, respectively.</figcaption></figure>

2000–2015) with the ASTERIX mass loss rates for the same period. The ASTERIX (−5.5 ± 3.3 Gt yr⁻¹) and modelled (−4.8 Gt yr⁻¹) estimates agree for non-calving glaciers, which gives confidence in the modelling approach. For calving glaciers, the ASTERIX mass loss (−15.9 ± 4.8 Gt yr⁻¹) is considerably larger than the modelled one (−4.3 Gt yr⁻¹). Assuming that the Marzeion et al.⁶ SMB rate estimates are also valid for calving glaciers, we can provide a first-order estimate of the frontal ablation of Andes calving glaciers by subtracting the SMB from our geodetic mass loss summed with the subaqueous mass loss (approximated to be −3.1 ± 1.2 Gt yr⁻¹ in Patagonia for 2000–2016)²². Our first-order estimate of frontal ablation is therefore 14.7 ± 6.0 Gt yr⁻¹ for 2000–2015, assuming that the basal and internal processes are negligible.

## Influence of glacier mass loss on river runoff

To estimate the contribution of glacier mass loss to surface runoff¹⁰,¹², we selected seven major Andean river basins and compared the total discharge with the observed glacier changes for the subperiods 2000–2009 and 2009–2018 (Fig. 4 and Supplementary Table 4). The glacier contribution to streamflow can be split into two terms: the seasonally delayed component and the imbalance (or excess discharge) contribution⁴⁴,⁴⁵,⁴⁶. Our decadal mass loss estimates only provide an indication of the imbalance

contribution, which is a lower bound of the glacier contribution to the streamflow, in particular for winter accumulation regimes in which the seasonally delayed fraction due to snow and ice melt is the dominant contribution.

Assuming that all the measured glacier mass loss contributed to the discharge measured in the outlet gauging stations, the glacier imbalance in the Santa River basin contributed to 5 m³ s⁻¹ of the measured streamflow in each decade (that is, around 4% of the total annual discharge). For the four catchments analysed in the Dry Andes of Chile and Argentina (the San Juan, Mendoza, Aconcagua and Maipo rivers), the river discharge decreased markedly after 2009, with reductions that ranged between 28 and 46%. This worrisome river runoff reduction in these populated, semi-arid basins coincides with a substantial increase in glacier mass loss after 2009. These changes reflect the very dry conditions experienced in this region since 2010, which have collectively been termed the 'megadrought' and can also be observed in terms of rainfall, snow accumulation, lake levels, water table depths and vegetation productivity, among other indicators⁴⁷,⁴⁸. Our results suggest that the decrease in perennial ice storage in these basins since 2009 effectively helped to mitigate the magnitude of this widespread dry period by contributing 3–8% of the total discharge recorded in these rivers. As indicated in previous analyses⁴⁹,⁵⁰, the glacier contribution
<figure><img src="image_4.png" /><figcaption>**Fig. 3** | Glacier mass change rate estimates from multiple methodologies. **a**, Low latitude (10°N–25°S, 99% Andean glaciers) regions. **b**, Southern Andes (25°S–56°S) regions. The results from this study do not account for subaqueous mass loss. In **a**, the estimates from Braun et al.<sup>22</sup> and this study correspond to the Andean low latitude glaciers only.</figcaption></figure>

was probably substantially larger during the warm months, and in particular during the several extremely dry late summers that have occurred in this region since 2010<sup>47,48</sup>. The two Patagonian basins show a slight increase on the decadal glacier imbalance contribution (from 2 to 3% in the Baker river that drains the NPI and from 14 to 16% in that Santa Cruz river that drains the SPI). These results reflect the more negative glacier mass balance rates observed during the most recent decade.

## Contribution of Andean glaciers to sea level rise

The total Andes glacier mass loss of −22.9 ± 5.9 Gt yr⁻¹ (−26.0 ± 6.0 Gt yr⁻¹ including the subaqueous loss estimated in Braun et al.<sup>22</sup>) during 2000–2018 corresponds to 0.06 ± 0.02 mm yr⁻¹ of SLR (subaqueous losses are not taken into account), slightly more than 10% of the global glacier contribution (0.55 mm yr⁻¹) during 2002–2016<sup>1</sup>. This value is 50% larger than the SLR contribution of all the High Mountain Asia glaciers<sup>29</sup>, although the latter cover an area three times larger. Our results also show that the total contribution of the Andean glaciers to the SLR was relatively steady during the first two decades of the twenty-first century, largely due to the unabated glacier mass loss in Southern Patagonia, where the largest glaciers are located. Nonetheless, within this large-scale pattern of unabated ice mass loss, we found an extensive region (between ~26° and 45°S) in which glaciers remained relatively stable or slightly gained mass from 2000 to 2009 (Fig. 1). This behaviour is similar to that observed in the West Kunlun and Karakoram

<figure><img src="image_5.png" /><figcaption>**Fig. 4** | The seven main river basins used in the hydrological analyses. Santa basin (Peru), Aconcagua, Maipo and Baker basins (Chile) and San Juan, Mendoza and Santa Cruz basins (Argentina). The blue colour scale represents the relative decadal glacier imbalance contribution (%) to the river runoff for each subperiod. Subregions as in Fig. 2.</figcaption></figure>

regions of the high mountains of Asia<sup>29,51</sup>, and puts these regions among the few glacierized areas on Earth in which the glaciers were stable or slightly gained in mass at the beginning of the twenty-first century.

Our new multiperiod and high-resolution data set provides original, up-to-date information on Andes-wide glacier changes for the first two decades of the twenty-first century. These estimates have numerous potential applications for local, regional and larger-scale assessments, and will contribute to better understand the response of Andean glaciers to climate change. They also have important implications for the management and conservation of Andean water resources, especially in the semi-arid regions. The highly resolved maps of glacier elevation changes for the entire Andes will also provide useful information to validate glaciological and hydrological models throughout South America.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, statements of code and data availability and associated accession codes are available at https://doi.org/10.1038/s41561-019-0432-5.

Received: 11 December 2018; Accepted: 22 July 2019; Published online: 16 September 2019
## References

1. Clapperton, C. M. The glaciation of the Andes. *Quat. Sci. Rev.* **2**, 83–155 (1983).

2. Garreaud, R. D., Vuille, M., Compagnucci, R. & Marengo, J. Present-day South American climate. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **281**, 180–195 (2009).

3. Bamber, J., Westaway, R., Marzeion, B. & Wouters, B. The land ice contribution to sea level during the satellite era. *Environ. Res. Lett.* **13**, 063008 (2018).

4. Gardner, A. S. et al. A reconciled estimate of glacier contributions to sea level rise: 2003 to 2009. *Science* **340**, 852–857 (2013).

5. Jacob, T., Wahr, J., Pfeffer, W. T. & Swenson, S. Recent contributions of glaciers and ice caps to sea level rise. *Nature* **482**, 514–518 (2012).

6. Marzeion, B., Leclercq, P. W., Cogley, J. G. & Jarosch, A. H. Brief communication: global reconstructions of glacier mass change during the 20th century are consistent. *Cryosphere* **9**, 2399–2404 (2015).

7. Reager, J. T. et al. A decade of sea level rise slowed by climate-driven hydrology. *Science* **351**, 699–703 (2016).

8. Glasser, N. F., Harrison, S., Jansson, K. N., Anderson, K. & Cowley, A. Global sea-level contribution from the Patagonian icefields since the Little Ice Age maximum. *Nat. Geosci.* **4**, 303–307 (2011).

9. Zemp, M. et al. Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016. *Nature* **568**, 382–386 (2019).

10. Buytaert, W. et al. Glacial melt content of water use in the tropical Andes. *Environ. Res. Lett.* **12**, 114014 (2017).

11. Huss, M. & Hock, R. Global-scale hydrological response to future glacier mass loss. *Nat. Clim. Change* **8**, 135–140 (2018).

12. Kaser, G., Großhauser, M. & Marzeion, B. Contribution potential of glaciers to water availability in different climate regimes. *Proc. Natl Acad. Sci. USA* **107**, 20223–20227 (2010).

13. Huss, M. et al. Toward mountains without permanent snow and ice. *Earth's Future* **5**, 418–435 (2017).

14. Vuille, M. et al. Rapid decline of snow and ice in the tropical Andes—impacts, uncertainties and challenges ahead. *Earth-Sci. Rev.* **176**, 1–95–213 (2018).

15. Cogley, J. G. Geodetic and direct mass-balance measurements: comparison and joint analysis. *Ann. Glaciol.* **50**, 96–100 (2009).

16. Mernild, S. H. et al. Mass loss and imbalance of glaciers along the Andes Cordillera to the sub-Antarctic islands. *Glob. Planet. Change* **133**, 109–119 (2015).

17. Rabatel, A. et al. Current state of glaciers in the tropical Andes: a multi-century perspective on glacier evolution and climate change. *Cryosphere* **7**, 81–102 (2013).

18. Zemp, M. et al. Historically unprecedented global glacier decline in the early 21st century. *J. Glaciol.* **61**, 745–762 (2015).

19. Rignot, E., Rivera, A. & Casassa, G. Contribution of the Patagonia Icefields of South America to sea level rise. *Science* **302**, 434–437 (2003).

20. Malz, P. et al. Elevation and mass changes of the Southern Patagonia Icefield derived from TanDEM-X and SRTM data. *Remote Sens.* **10**, 188 (2018).

21. Melkonian, A. K. et al. Satellite-derived volume loss rates and glacier speeds for the Cordillera Darwin Icefield, Chile. *Cryosphere* **7**, 823–839 (2013).

22. Braun, M. H. et al. Constraining glacier elevation and mass changes in South America. *Nat. Clim. Change* **9**, 130–136 (2019).

23. Ivins, E. R. et al. On-land ice loss and glacial isostatic adjustment at the Drake Passage: 2003–2009. *J. Geophys. Res.* **116**, B02403 (2011).

24. Richter, A. et al. Crustal deformation across the Southern Patagonian Icefield observed by GNSS. *Earth Planet. Sci. Lett.* **452**, 206–215 (2016).

25. Huss, M. & Hock, R. A. A new model for global glacier change and sea-level rise. *Front. Earth Sci.* **3**, 54 (2015).

26. Mernild, S. & Wilson, R. The Andes Cordillera. Part III: glacier surface mass balance and contribution to sea level rise (1979-2014). *Int. J. Climatol.* **37**, 3154–3174 (2016).

27. Collao-Barrios, G. et al. Ice flow modelling to constrain the surface mass balance and ice discharge of San Rafael Glacier, Northern Patagonia Icefield. *J. Geol.* **64**, 568–582 (2018).

28. Schaefer, M., Machguth, H., Falvey, M., Casassa, G. & Rignot, E. Quantifying mass balance processes on the Southern Patagonia Icefield. *Cryosphere* **9**, 25–35 (2015).

29. Brun, F., Berthier, E., Wagnon, P., Kääb, A. & Treichler, D. A spatially resolved estimate of High Mountain Asia glacier mass balances from 2000 to 2016. *Nat. Geosci.* **10**, 668–673 (2017).

30. Berthier, E., Larsen, C., Durkin, W. J., Willis, M. J. & Pritchard, M. E. Brief communication: unabated wastage of the Juneau and Stikine icefields (southeast Alaska) in the early 21st century. *Cryosphere* **12**, 1523–1530 (2018).

31. Dussaillant, I., Berthier, E. & Brun, F. Geodetic mass balance of the Northern Patagonian Icefield from 2000 to 2012 using two independent methods. *Front. Earth Sci.* **6**, 8 (2018).

32. Masiokas, M. H. et al. Glacier fluctuations in extratropical South America during the past 1000 years. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **281**, 242–268 (2009).

33. Sagredo, E. A. & Lowell, T. V. Climatology of Andean glaciers: a framework to understand glacier response to climate change. *Glob. Planet. Change* **86**, 101–109 (2012).

34. Jaber, W. A., Floricioliu, D. & Rott, H. Geodetic mass balance of the Patagonian icefields derived from SRTM and TanDEM-X data. In 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS) 342–345 (IEEE, 2016); https://doi.org/10.1109/IGARSS.2016.7729082

35. Rabatel, A. et al. Toward an imminent extinction of Colombian glaciers? *Geogr. Annal. A* **100**, 75–95 (2017).

36. Cook, S. J., Kougkoulos, I., Edwards, L. A., Dortch, J. & Hoffmann, D. Glacier change and glacial lake outburst flood risk in the Bolivian Andes. *Cryosphere* **10**, 2399–2413 (2016).

37. Paul, F. & Mölg, N. Hasty retreat of glaciers in northern Patagonia from 1985 to 2011. *J. Glaciol.* **60**, 1033–1043 (2014).

38. Burger, F. et al. Interannual and decadal variability in glacier contribution to runoff from high-elevation Andean catchments: understanding the role of debris cover in glacier hydrology. *Hydrol. Process.* **33**, 214–229 (2019).

39. Farías-Barahona, D. et al. Geodetic mass balances and area changes of Echaurren Norte glacier (Central Andes, Chile) between 1955 and 2015. *Remote Sens.* **11**, 260 (2019).

40. Masiokas, M. H. et al. Reconstructing the annual mass balance of the Echaurren Norte glacier (Central Andes, 33.5°S) using local and regional hydroclimatic data. *The Cryosphere* **10**, 927–940 (2016).

41. Wouters, B., Gardner, A. S. & MoHoldt, G. Global glacier mass loss during the GRACE satellite mission (2002–2016). *Front. Earth Sci.* **7**, 96 (2019).

42. Dietrich, R. et al. Rapid crustal uplift in Patagonia due to enhanced ice loss. *Earth Planet. Sci. Lett.* **289**, 22–29 (2010).

43. Ivins, E. R. et al. On-land ice loss and glacial isostatic adjustment at the Drake Passage: 2003–2009. *J. Geophys. Res.* **116**, B02403 (2011).

44. Pritchard, H. D. Asia's shrinking glaciers protect large populations from drought stress. *Nature* **569**, 649–654 (2019).

45. Radić, V. & Hock, R. Glaciers in the Earth's hydrological cycle. Assessments of glacier mass and runoff changes on global and regional scales. *Surv. Geophys.* **35**, 813–837 (2014).

46. Lambrecht, A. & Mayer, C. Temporal variability of the non-steady contribution from glaciers to water discharge in western Austria. *J. Hydrol.* **376**, 353–361 (2009).

47. Garreaud, R. D. et al. The 2010–2015 megadrought in central Chile: impacts on regional hydroclimate and vegetation. *Hydrol. Earth Syst. Sci.* **21**, 6307–6327 (2017).

48. Rivera, J. A., Penalba, O. C., Villalba, R. & Araneo, D. C. Spatio-temporal patterns of the 2010–2015 extreme hydrological drought across the Central Andes, Argentina. *Water* **9**, 652 (2017).

49. Gascoin, S. et al. Glacier contribution to streamflow in two headwaters of the Huasco River, Dry Andes of Chile. *Cryosphere* **5**, 1099–1113 (2011).

50. Ragetti, S. & Pelliccioli, F. Calibration of a physically based, spatially distributed hydrological model in a glacierized basin: on the use of knowledge from glaciometerological processes to constrain model parameters. *Water Resour. Res.* **48**, W03509 (2012).

51. Kääb, A., Treichler, D., Nuth, C. & Berthier, E. Brief communication: contending estimates of 2003–2008 glacier mass balance over the Pamir–Karakoram–Himalaya. *Cryosphere* **9**, 557–564 (2015).

## Acknowledgements

We acknowledge the French Space Agency (CNES) and the Région Occitanie for PhD fellowship to I.D. E.B. also acknowledges the support from the CNES and the French Programme National de Télédétection Spatiale grant PNTS-2016-01. L.R., P.P. and M.M. acknowledge the support from Agencia de Promoción Científica (projects PICT 2010-1438 and PICT 2014-1794) and CONICET. ASTER images are courtesy of NASA/ METI/AIST/Japan Space systems and the US/Japan ASTER Science Team. This work was only possible thanks to the GLIMS project (www.glims.org/) which allowed the population of a vast archive of freely available ASTER stereo images over the glaciers. We thank the Ames Stereo Pipeline support and developer teams for their help. Most of the computation was performed using the Froggy platform of the CIMENT infrastructure (http://ciment.ujf-grenoble.fr), supported by the Rhone-Alpes region (grant CPER07_13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the Programme Investissements d'Avenir supervised by the French Agence National pour la Recherche. We are grateful to K. Moxham for reviewing the English.

## Author contributions

I.D. and E.B. designed the study. I.D. performed the ASTER DEM analysis with inputs from F.B., R.H. and E.B. M.M. made the streamflow analysis. All the authors contributed to the interpretation of the results. I.D. led the writing of the paper and all the other co-authors contributed to it.

## Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41561-019-0432-5.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Correspondence and requests for materials** should be addressed to I.D.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2019
## Methods

The data and methodology (ASTERIX) applied here to the entire Andes was inspired by previous studies<sup>22–54</sup>. Below we briefly summarize the processing steps and highlight the specificity for the application to the Andes. Brun et al.<sup>29</sup> gives further details.

The ASTER DEMs were processed based on 1° by 1° (latitude by longitude) tiles. The entire processing was done in Universal Transverse Mercator at a resolution of 30 m (Supplementary Information). Tiles that contained less than 10 km² of glacier were not processed and were substituted by their regional average (Supplementary Information). We derived about 30,000 DEMs from raw ASTER level L1A (freely available on Earthdata)<sup>55</sup> acquired from March 2000 to April 2018 using the open source Ames Stereo Pipeline routine<sup>56</sup>. After co-registration of individual ASTER DEMs to the SRTM DEM<sup>57</sup>, we fitted a two-step linear regression through the time series of the available ASTER DEMs to estimate the rate of elevation change for each pixel. These trends were extracted after excluding unrealistic outliers from the ASTER DEMs, that is, pixels for which the absolute elevation difference was larger than 25 m plus 60 m per year as the last valid elevation (initialized by the elevation value of the SRTM DEM). The latter value was chosen based on the highest observed thinning rate in the Andes of approximately 50 m yr<sup>-1</sup>. A first linear trend was fitted to the elevation time series from ASTER DEMs and the SRTM DEM. Then, a second linear regression was fitted excluding the SRTM DEM and all the data outside the 99% confidence interval of the first fit. Exclusion of the SRTM DEM is generally recommended to avoid the biases due to radar penetration into dry snow and ice<sup>30</sup>. Exceptionally for the 1° by 1° tiles south of 48° S (SPI to the South), we kept the SRTM DEM when computing the final trend. This choice is justified by the limited number of ASTER DEMs in these subregions due to a higher occurrence of clouds and reduced image acquisition and the lack of measurable penetration of the SRTM C-band radar signal on the wet surface of the Patagonian icefields in February 2000, the summer period in the Southern Hemisphere<sup>31,58</sup>. The linear regression was computed for the entire period (March 2000–April 2018) and two subperiods of about nine years (March 2000–March 2009 and January 2009–April 2018).

We then integrated the rates of elevation change using the local hypsometric method<sup>29</sup> with 100 m elevation bands. Data gaps were thus filled on an individual glacier basis using the mean rate of elevation change of the corresponding elevation bin, proved to be one of the most robust methods of void interpolation<sup>30</sup>. The mean elevation change rate was computed after excluding all the values that lay further than three normalized median absolute deviations from the median of the elevation bands. We also excluded pixels on slopes larger than 45°, pixels with uncertainties in the linear fit larger than 3.0 m yr<sup>-1</sup> (at the 95% confidence level) and unrealistic absolute dh/dt values larger than 30 m yr<sup>-1</sup>. Excluding glacier HPS12, the most negative elevation change rates (−28 m yr<sup>-1</sup>) were observed on the Upsala and Jorge Mont glaciers (SPI). Glacier HPS12 was the only glacier that presented stronger elevation change rates (up to −44 m yr<sup>-1</sup>) and was processed separately. The total individual glacier volume change rate was then calculated as the sum of the mean elevation change rates multiplied by the area of the elevation band to which it corresponded. In the rare cases in which no elevation change rates were available for an altitude band, we interpolated hypometrically using a linear function forced to 0 for the highest band, which corresponded to the glacier accumulation zone. Finally, the mass balance rate was computed assuming a volume to mass conversion factor of 850 ± 60 kg m<sup>-3</sup> (ref. 60). Mass balance rates for large regions (that is, the entire Andes, Andes subregions, river basins, 1° by 1° tiles) were computed as the weighted sum of all individual glaciers. In the specific case of 1° by 1° tiles, all the glaciers whose geometric centroids were located inside the tile limits were considered. Glaciers that presented mass balance errors higher than 1 m w.e. yr<sup>-1</sup> and a measured area lower than 20% of their total surface area were excluded and their mass balance rate replaced by the regional mean.

For the DEM adjustment over stable terrain and mass balance calculations, we used a glacier mask derived from the Southern Andes and low latitudes regions of RGI v6.0 (excluding glaciers in Africa and Indonesia, less than 1% of the total glacierized area in the low latitudes region). Water-body outlines from the Global Lakes and Wetlands Database, levels 1 and 2<sup>61</sup>, and ocean water polygons<sup>62</sup> were also used to define stable terrain. We selected the RGI inventory for three reasons: RGI is available for the entire Andes<sup>63</sup>, it represents the glacier extent around 2000, the year of the earliest ASTER DEMs, and allows the comparison to model results also based on RGI. Areal changes were not considered in our mass balance change estimates due to the lack of a more recent inventory for the entire Andes. The sensitivity of our mass balance estimates to areal changes and to the imperfections inherent to the RGI are discussed in the Supplementary Information.

**Uncertainty assessment.** The total uncertainty of the mass balance rates for our seven subregions, the selected river basins and the large icefields was calculated as the quadratic sum of the random ($\sigma_{M_{\text{rdm}}}$) and the triangulation error ($\sigma_{M_{\text{trg}}}$) (refs. 29,57):

$$ \sigma_{M_{\text{tot}}} = \sqrt{(\sigma_{M_{\text{rdm}}})^2 + (\sigma_{M_{\text{trg}}})^2} \quad (1) $$

The random and triangulation errors were calculated by considering subregions as a single ice body. The random error has three main sources, assumed to be independent of one another: the uncertainty in the rate of elevation change, the uncertainty in the volume-to-mass conversion factor and the uncertainty in the glacierized area. This calculation of the uncertainty in the rate of elevation change is derived from Fischer et al.<sup>64</sup> and Rolstad et al.<sup>65</sup>, with details given in Brun et al.<sup>29</sup>.

The triangulation error was assessed as the absolute value of the triangulation residual, r, between two equally long subperiods:

$$ r = \left| \frac{(\dot{M}_{2000-2009} + \dot{M}_{2009-2018})}{2} - \dot{M}_{2000-2018} \right| \quad (2) $$

where $\dot{M}_{y_1-y_2}$ is the mass balance rate (m w.e. yr<sup>-1</sup>) for the period between year y<sub>1</sub> and year y<sub>2</sub>. Similar to Menounos et al.<sup>66</sup>, we calculated the residuals for every individual glacier in the Andes and found that r is dependent on the number of DEMs used to calculate $\dot{M}_{2000-2018}$, increasing when mass balances were calculated using fewer DEMs. We derived the 68th percentile of the residual distribution in a moving window of two DEMs. The triangulation error was then approximated as the linear relation of the observed trend of the 68th percentile of residuals for each window (Supplementary Information and Supplementary Fig. 5). Our triangulation error varies from 0.27 m w.e. yr<sup>-1</sup> when the trend is fitted on three DEMs to 0.08 m w.e. yr<sup>-1</sup> for 60 DEMs. The subregional triangulation error was then calculated using this relation and the mean number of DEMs for every subregion.

The total uncertainty of our Andes-wide glacier mass change ($G_{\text{tyr}}$<sup>-1</sup>) was calculated as the quadratic sum of the total uncertainty on the glacier mass change for each Andes subregion. The total uncertainty of the specific mass balance (m w.e. yr<sup>-1</sup>) was then calculated by conserving the relative uncertainty.

## Data availability

Rate of elevation change maps are distributed through PANGAEA platform (https://doi.pangaea.de/10.1594/PANGAEA.903618). Individual glacier mass balances will be provided through the World Glacier Monitoring Service database (https://wgms.ch/data_exploration/). The code used to generate ASTER DEMs from the freely available 1A images can be downloaded at https://github.com/FannyBrun/ASTER_DEM_from_L1A. Individual ASTER DEMs are also available upon request to the corresponding author (I.D.).

## References

52. Naimura, T., Fujita, K., Yamaguchi, S. & Sharma, R. R. Elevation changes of glaciers revealed by multitemporal digital elevation models calibrated by GPS survey in the Khumbu region, Nepal Himalaya, 1992-2008. *J. Glaciol.* **58**, 648–656 (2012).

53. Wang, D. & Kääb, A. Modeling glacier elevation change from DEM time series. *Remote Sens.* **7**, 10117–10142 (2015).

54. Willis, M. J., Melkonian, A. K., Pritchard, M. E. & Ramage, J. M. Ice loss rates at the Northern Patagonian Icefield derived using a decade of satellite remote sensing. *Remote Sens. Environ.* **117**, 184–198 (2012).

55. Earthdata (NASA, accessed 3 August 2018); https://earthdata.nasa.gov/

56. Shean, D. E. et al. An automated, open-source pipeline for mass production of digital elevation models (DEMs) from very-high-resolution commercial stereo satellite imagery. *ISPRS J. Photogram. Remote Sens.* **116**, 101–117 (2016).

57. Nuth, C. & Kääb, A. Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change. *Cryosphere* **5**, 271–290 (2011).

58. Jaber, W. A., Floricioiu, D., Rott, H. & Eineder, M. Surface elevation changes of glaciers derived from SRTM and TanDEM-X DEM differences. In *2013 IEEE International Geoscience and Remote Sensing Symposium (IGARSS)* 1893–1896 (IEEE, 2013); https://doi.org/10.1109/IGARSS.2013.6723173

59. McNabb, R., Nuth, C., Kääb, A. & Girod, L. Sensitivity of glacier volume change estimation to DEM void interpolation. *Cryosphere* **13**, 895–910 (2019).

60. Huss, M. Density assumptions for converting geodetic glacier volume change to mass change. *Cryosphere* **7**, 877–887 (2013).

61. Lehner, B. & Döll, P. Development and validation of a global database of lakes, reservoirs and wetlands. *J. Hydrol.* **296**, 1–22 (2004).

62. OpenStreetMap (OpenStreetMap, accessed 8 December 2018); https://www.openstreetmap.org/copyright

63. Pfeffer, W. T. et al. The Randolph Glacier Inventory: a globally complete inventory of glaciers. *J. Glaciol.* **60**, 537–552 (2014).

64. Fischer, M., Huss, M. & Hoelzle, M. Surface elevation and mass changes of all Swiss glaciers 1980–2010. *Cryosphere* **9**, 525–540 (2015).

65. Rolstad, C., Haug, T. & Denby, B. Spatially integrated geodetic glacier mass balance and its uncertainty based on geostatistical analysis: application to the western Svartisen ice cap, Norway. *J. Glaciol.* **55**, 666–680 (2009).

66. Menounos, B. et al. Heterogeneous changes in western North American glaciers linked to decadal variability in zonal wind strength. *Geophys. Res. Lett.* **46**, 200–209 (2019).