# Increased mass loss of glaciers in Volcán Domuyo (Argentinian Andes) between 1962 and 2020, revealed by aerial photos and satellite stereo imagery

Daniel Falaschi, Etienne Berthier, Joaquín Belart, Claudio Bravo, Mariano

Castro, Marcelo Durand, Ricardo Villalba

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

► To cite this version:

Daniel Falaschi, Etienne Berthier, Joaquín Belart, Claudio Bravo, Mariano Castro, et al. In-creased mass loss of glaciers in Volcán Domuyo (Argentinian Andes) between 1962 and 2020, re-vealed by aerial photos and satellite stereo imagery. Journal of Glaciology, 2023, 69 (273), pp.40-56. 10.1017/jog.2022.43 . hal-03719615

HAL Id: hal-03719615

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

Submitted on 11 Jul 2022

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Article

**Cite this article:** Falaschi D, Berthier E, Belart JMC, Bravo C, Castro M, Durand M, Villalba R (2022). Increased mass loss of glaciers in Volcán Domuyo (Argentinian Andes) between 1962 and 2020, revealed by aerial photos and satellite stereo imagery. *Journal of Glaciology* 1–17. https://doi.org/10.1017/jog.2022.43

Received: 16 December 2021

Revised: 22 May 2022

Accepted: 23 May 2022

**Key words:**

Glacier mass balance; glacier modelling;
glacier volume; mountain glaciers

**Author for correspondence:**

Daniel Falaschi,

E-mail: dfalaschi@mendoza-conicet.gob.ar

Increased mass loss of glaciers in Volcán
Domuyo (Argentinian Andes) between 1962 and
2020, revealed by aerial photos and satellite
stereo imagery

Daniel Falaschi¹,²,³ ID, Etienne Berthier⁴ ID, Joaquín M. C. Belart⁵,⁶ ID,
Claudio Bravo⁷ ID, Mariano Castro¹, Marcelo Durand⁸ and Ricardo Villalba¹ ID

¹Instituto Argentino de Nivología, Glaciología y Ciencias Ambientales, Mendoza, Argentina; ²Departamento de
Geografía, Facultad de Filosofía y Letras, Universidad Nacional de Cuyo, Mendoza, Argentina; ³School of
Geography and Sustainable Development, University of St. Andrews, St. Andrews, Scotland, UK; ⁴Laboratoire
d'Etudes en Géophysique et Océanographie Spatiales, Université de Toulouse, CNES, CNRS, IRD, UPS, Toulouse,
France; ⁵National Land Survey of Iceland, Akranes, Iceland; ⁶Institute of Earth Sciences, University of Iceland,
Reykjavik, Iceland; ⁷Centro de Estudios Científicos, Valdivia, Chile and ⁸Departamento de Geodesia y Topografía,
Facultad de Ciencias Exactas y Tecnología, Universidad Nacional de Tucumán, San Miquel de Tucuman, Argentina

Abstract

We present geodetic mass-balance estimates for ten glaciers (22.6 ± 1.1 km²) around Volcán Domuyo between 1962 and 2020 (and 46 glaciers covering 29 ± 1.5 km² between 1984 and 2020), derived from airborne, ASTER and Pléiades imagery. Overall, we find a slightly negative mass balance (−0.15 ± 0.09 m w.e. a⁻¹) for the entire 1962–2020 time span. A closer inspection of sub-periods reveals, however, an increasingly negative mass balance with time. The Domuyo glaciers shifted from a moderately positive mass balance of +0.28 ± 0.13 m w.e. a⁻¹ between 1962 and 1984, to a strongly negative mass balance of −0.99 ± 0.19 m w.e. a⁻¹ between 2012 and 2020. An increase in summer temperatures and a decrease in winter precipitation during the last four decades are likely drivers of the observed glacier changes. We support this finding by implementing a minimal glacier model, relying solely on monthly precipitation and air temperatures. The mass-balance evolution detected in Volcán Domuyo is consistent with other sites spread across the Central Andes, suggesting rapidly increasing glacier wastage impacts are occurring at a geographically wider scale.

**1. Introduction**

The Andes of South America host a great variety of glaciers in relation to their morphology, size, altitudinal distribution, climatic regime and response time to climatic perturbations (Sagredo and Lowell, 2012). They are relevant for their role as water resources in semi-arid regions, hydro-electric power generation, human consumption and various industrial activities (Masiokas and others, 2020 and references therein). Andean glaciers have also been acknowledged as one of the major contributors to sea level rise (Gardner and others, 2013; Zemp and others, 2019; Hugonnet and others, 2021). Traditionally, glacier mass balance is determined by means of the glaciological method, which uses a combination of stakes and snow pits to measure the accumulation and ablation components (Cogley, 2009). Owing to the inherent difficulties in sustaining continuous, long-term surveys, glaciological mass-balance records in the Southern Andes are infrequent, especially those that surpass a full decade of measurements (WGMS, 2021).

In recent years, the geodetic method, which retrieves glacier surface topography changes over different years from multi-sourced digital elevation models (DEMs), has allowed for significant advances in regard to understanding glacier mass-balance trends at different latitudes of the Southern Andes in Chile and Argentina, and their relationship to climate variability, at local to regional scales (e.g. Melkonian and others, 2013; Falaschi and others, 2017, 2019; Dussaillant and others, 2018, 2019; Braun and others, 2019; Kinnard and others, 2020). These studies have compensated, to a certain degree, for the overall paucity of direct glaciological measurements. The geodetic method can serve as a validation or calibration tool for the glaciological method (Zemp and others, 2013; Andreassen and others, 2016; Klug and others, 2018). Unfortunately, despite a few exceptions (e.g. Rignot, 2003; Falaschi and others, 2019; Farías Barahona and others, 2019, 2020a, 2020b), most geodetic surveys in the Southern Andes usually span <20 years, and do not extend well back into the 20th century, hampering a long-term, comprehensive insight into glacier trends and the glacier-climate relationships.

One of the least sampled regions of the Andes in glaciological research is between the Central, Lake District and Patagonian Andes, between ~36° S and 42° S (Lliboutry, 1998). This region is characterized by a significant drop in average elevation of the mountain peaks relative to the theoretical equilibrium line altitude (ELA) (Condom and others, 2007) and consequently, the scarce glacier ice extent is restricted to the highest volcanic edifices (Reinthaler and others, 2019) of the Southern Volcanic Zone (33–46° S; Cembrano and

© The Author(s), 2022. Published by
Cambridge University Press. This is an Open
Access article, distributed under the terms of
the Creative Commons Attribution-
NonCommercial-NoDerivatives licence
(https://creativecommons.org/licenses/by-nc-
nd/4.0/), which permits non-commercial re-
use, distribution, and reproduction in any
medium, provided the original work is
unaltered and is properly cited. The written
permission of Cambridge University Press
must be obtained for commercial re-use or in
order to create a derivative work.

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Overview of the Volcán Domuyo area and major glaciers shown in the insets (a1–a3). Cv, Covunco; T, Turbio; MC, Manchana Covunco; Ch, Chadileo. Background: Pléiades © CNES 2020, Distribution Airbus D&S (a). Location of the study site in Northern Patagonia (b). Location of geodetic mass-balance sites with >20 year records along the Chilean–Argentinean Andes (c).</figcaption></figure>

Lara, 2009). It is therefore not surprising that glaciological studies on this portion of the Andes are sparse, focusing mainly on glaciier area changes (Rivera and others, 2012; Rivera and Bown, 2013; Reinthaler and others, 2019). Only recent large-scale remote-sensing studies have looked at glacier mass changes in this region (Braun and others, 2019; Dussaillant and others, 2019; Ferri and others, 2020) but after 2000 only.

In this study we aim to generate a geodetic mass-balance record between 1962 and 2020 for the glaciers of Volcán Domuyo (~36° S), based on the differencing of DEMs derived from archives of aerial photographs, ASTER and Pléiades stereo-imagery. Our analysis covers four time intervals during the last six decades, enabling for the investigation of possible climatic drivers of glacier change on a sub-decadal to bidecadal basis. We also use the above-mentioned imagery to produce a high-resolution, updated glacier inventory and to derive glacier area changes during each sub-period. We analyse our mass and area change results for Volcán Domuyo in a broader geographical context in the Central and Northern Patagonian Andes of Chile and Argentina.

A second, main goal of this study was to apply a surface mass balance (SMB) model for four of the Volcán Domuyo glaciers during a period analogous to the geodetic results and to assess the observed differences in glacier mass balance. Lastly, we estimate climate trends for the Volcán Domuyo area using available nearby instrumental records and gridded climate products to evaluate the possible influence of climatic trends in the region on the observed glacier changes.

## 2. Study area

Volcán Domuyo (36°38′ S, 70°25′ W, 4706 m a.s.l.) lies in a relatively remote location ~50 km east of the main Andean axis

(Fig. 1). Although from a strict geographic definition Volcán Domuyo is located at the northern tip of (but still within) Patagonia, the high aridity of the area and the extensive debris-covered and rock glaciers (Falaschi and others, 2016) are more typical of the Central Andes of Argentina and Chile (Dussaillant and others, 2019). Unsurprisingly, the Andean range where the volcano is located is often termed *Transition Cordillera*. Volcán Domuyo, nicknamed ‘The Roof of Patagonia’, is the highest peak in Patagonia, and forms part of a protected natural conservation area located in the Argentinean Neuquén Province at the head of the Rio Neuquén catchment. Volcán Domuyo sits within a vast Middle Miocene to Early Pliocene Volcanic complex. A striking feature of the area is the widespread occurrence of hydrothermal activity, including boiling fluids that enter creeks, maintaining >30°C temperatures (Chiodini and others, 2014). Snow and ice from the Domuyo district nourish the Neuquén river, where construction of hydropower plants is scheduled for the near future. The Rio Neuquén catchment has a pluvio-nival regime, and is characterized by a seasonally bimodal discharge.

Flow is concentrated from May to December, with a major rainfall peak in July–August and a minor snowfall peak in October–November. The lowest flows are recorded between February and April, representing only 8.5% of the total annual discharge (AIC, 2010). The glacier aerial views are courtesy of Hector Valdez.

Beyond the pan-Andes studies of Braun and others (2019) and Dussaillant and others (2019), a small amount of glaciological research has been conducted on ice-capped volcanoes between 36° S and 42° S. In a survey aimed at awareness-raising about glaciier-volcanic interactions and associated hazards, Rivera and others (2015) estimated the overall thinning since 1961 in Volcán Villarrica (39° S), whereas Schaefer at al. (2017) carried
out in situ measurements and also modelled the mass balance of the Mocho Glacier atop Volcán Mocho Choshuenco (40° S) between 2006 and 2015.

In spite of its dome-shaped appearance, Volcán Domuyo does not host a typical ice cap with approximately radial ice flow as in many glacierized volcanoes in Patagonia (Reinthaler and others, 2019), rather it hosts mountain glaciers and ice aprons. Falaschi and others (2016) produced an inventory of Domuyo glaciers accounting for ~25.4 km² of ice for 2009 and an estimated 26% (~1.4% a⁻¹) reduction since 1990. Another unusual characteristic of these glaciers is that debris-covered ice is not restricted to the flatter glacier snouts, but has been proportionally increasing in their upper (and steeper) reaches as well in recent years. These authors also hinted at significant glacier thinning in recent decades as well as increasingly larger internal rock outcrops in accumulation areas, both of which are indicative of rising ELAs for out-of-balance glaciers. However, a more comprehensive and in-depth glaciological investigation of glacier changes in the Volcán Domuyo is still needed.

From a climatic point of view, the Andean mountains south of ~30° S are permanently under the influence of the mid-latitude westerly flow. According to Garreaud (2009), the main characteristic of the entire region is the passage of frontal systems that are driven eastwards from the Pacific Ocean, producing precipitation maxima during the austral winter months (July–September). In particular, the Domuyo area exhibits a higher amount of moisture compared to the Central Andes to the north, but has drier conditions in comparison with the Patagonian Andes immediately to the south (hence the *Transition Cordillera* term). There is also a strong northsouth gradient of the ELA, decreasing from ~5000 m a.s.l. at 31° S to <2000 m a.s.l. at 49° S (Sagredo and Lowell, 2012). Around 36° S, the latitude of Volcán Domuyo, a rough estimate of 3200–3300 m a.s.l. for the ELA can be inferred from Condom and others (2007).

# 3. Data and methods

## 3.1 DEMs derived from optical satellite imagery

DEMs and orthoimages were generated using Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) (2003) and Pléiades (2012 and 2020) stereo imagery. We generated a 30 m resolution DEM and 15 m resolution orthoimage from the 9 April 2003 ASTER stereopair using the NASA Ames Stereo Pipeline (ASP, Shean and others, 2016) using the Semi-Global-Matching algorithm following Shean and others (2020). Altimetric precision of individual ASTER DEMs ranges from 5 to 15 m (Girod and others, 2017; Shean and others, 2020). The two most recent DEMs were derived from two Pléiades stereopairs with 0.7 m spatial resolution, acquired on 7 March 2012 and 18 February 2020, respectively. Each Pléiades stereopair was processed using ASP to derive DEMs with a ground-sampling distance (GSD) of 4 m (DEM₂₀₁₂ and DEM₂₀₂₀ hereafter) and 0.5 m resolution orthoimages, using the

Semi-Global-Matching algorithm and the set of processing parameters derived in Deschamps-Berger and others (2020) without ground control points (GCPs). Even without GCPs, Pléiades-derived DEMs are located within ~10 m (90% circular error) and have an altimetric precision of ~1 m (Berthier and others, 2014). They have been already used to derive glacier elevation changes in diverse-glacierized regions worldwide (Kutuzov and others, 2019; Abermann and others, 2020; Andreassen and others, 2020; Kapitsa and others, 2020; Denzinger and others, 2021).

## 3.2 Generation of DEMs and orthoimages from aerial photographs

Vertical aerial photographs covering the Volcán Domuyo area are available from the Instituto Geográfico Nacional of Argentina (IGN) for 1962 and 1984 (Table 1). Both photogrammetric flights have a 60–70% ‘along track’ overlap approximately, show no cloud cover and good contrast on glaciers. The original negative films and camera calibration certificate are available for the 1984 flight. They were scanned at 10 µm resolution in a Vexcel UltraScan 5000 unit. Unfortunately, the original 1962 negative films were unavailable from the IGN, and instead contact prints were scanned at 2400 dpi resolution using a Hewlett Packard Designjet 4500 scanner.

The aerial photographs were processed with MicMac (IGN, France, Pierrot Deseilligny and Clery, 2011; Rupnik and others, 2017), following Belart and others (2019). The MicMac routines ‘tapioca’, ‘tapas’, ‘GCPBascule’, ‘Campari’ and ‘Malt’ perform the sequence of steps needed in the photogrammetric processing, consisting of the computation of tie point generation, solving of relative and absolute image orientation, and performing the dense matching and creation of the DEM and orthoimages. The processing is semi-automated: a series of 23 GCPs (1984) and 28 GCPs (1962) were manually digitized in the 2020 Pléiades orthoimage and in the corresponding pairs or triplet of aerial photographs (see Fig. S1). The elevation of each GCP was extracted from the 2020 Pléiades DEM. This scheme had been previously used by Papasodoro and others (2015), and resulted in photogrammetric DEMs and orthomosaics co-registered to the Pléiades 2020 dataset.

## 3.3 Glacier mapping

The basis for glacier mapping at the Volcán Domuyo was the outlines produced by Falaschi and others (2016), which were derived from 2.5 m spatial resolution ALOS PRISM scenes acquired in April 2009. These polygons were manually corrected and adjusted by visual interpretation for 1962, 1984, 2003, 2012 and 2020 using all processed orthoimages in a UTM Zone 19S projection. We then calculated glacier area change rates for each time interval and for the full study period. We report the glacier inventory and area change (and mass balance, see Section 3.4) for 46 glaciers for the 1984–2020 period. The 1962 images did not cover some of these

<table><caption>**Table 1.** Summary of aerial and satellite data and characteristics of the derived DEMs/orthoimages</caption><thead><tr><th>Dataset</th><th>Acquisition date (dd/mm/yyyy)</th><th>Number of images</th><th>DEM name and spatial resolution m</th><th>Orthoimage spatial resolution m</th></tr></thead><tbody><tr><td>WILD RC5/AG</td><td>08/04/1962</td><td>5</td><td>DEM<sub>1962</sub> (5)</td><td>0.5</td></tr><tr><td>WILD RC10/UAG II</td><td>11/03/1984</td><td>8</td><td>DEM<sub>1984</sub> (5)</td><td>0.5</td></tr><tr><td>ASTER</td><td>09/04/2003</td><td>2</td><td>DEM<sub>2003</sub> (30)</td><td>15</td></tr><tr><td>Pléiades</td><td>07/03/2012</td><td>2</td><td>DEM<sub>2012</sub> (4)</td><td>0.5</td></tr><tr><td>Pléiades</td><td>18/02/2020</td><td>2</td><td>DEM<sub>2020</sub> (4)</td><td>0.5</td></tr></tbody></table>

The footprint of the aerial (and satellite) imagery used in this study is displayed in Fig. S1.
glaciers, Manchana Covunco being the largest one (1.7 ± 0.1 km²). Still, from 1962 to 2020, we report area and mass changes for the ten largest glaciers (>0.5 km²) which total 22.6 km² in 1984, i.e. 78% of the total glacier area.

Surface conditions (cloud and snow cover) allowed for the clear delineation of glacier margins in the 1984 aerial photos, the 2003 ASTER images and Pléiades 2020 scenes. The Pléiades 2012 scene and 1962 aerial photos do nevertheless display a slight amount of seasonal snow. In the first of these two cases, glacier mapping was performed with the aid of a WorldView-2 4th April 2012 image available from GoogleEarth with very dry conditions. Finding an alternative dataset that can be used to validate the 1962 glacier outlines is far more problematic. We are nevertheless confident that our acute knowledge of the terrain involved, based on at least four field campaigns and the 1984 aerial photos, has allowed us to avoid much of the seasonal snow in the glacier mapping procedure. The glacier area uncertainty was calculated following Wagnon and others (2021) as the product of the glacier outline initial perimeter and two times the GSD of the satellite imagery used (e.g. 15 m for ASTER VIS bands and 0.5 m for the Pléiades panchromatic scene).

As for glacier denomination, there are no official glacier names in the Domuyo district. We therefore termed the largest glaciers following the names of the mountain creeks (and consecutive numbers in case of various glaciers) that originate from them (Fig. 1a).

### 3.4 DEM differencing and mass-balance calculation

The full 1962–2020 study period was divided into four time intervals: 1962–84 (21.94 years), 1984–2003 (19.09 years), 2003–12 (8.92 years) and 2012–20 (7.95 years), spanning 57.9 years in total. Since all DEMs and orthoimages are coregistered to the most recent Pléiades DEM$_{2020}$, any remaining shifts (which we term triangulation following Nuth and Kaab, 2011) among the various DEMs is minimized. For each individual time period, we subtracted the later DEM from the initial DEM, and produced elevation change (dh/dt) maps. The internal consistency of our dh/dt was examined by comparing the weighted sum of the individual dh/dt values for each time interval against the accumulated 1962–2020 dh/dt value.

dDEM maps can be affected by residual errors in the interior and exterior sensor orientations, which can cause localized biases in elevation. This effect of sensor misorientation can be particularly problematic for older frame cameras, where calibration data are often not well defined or directly missing (e.g. Belart and others, 2019). Previous geodetic mass-balance studies using both ASTER and Pléiades data have shown low-frequency biases in dh/dt residual maps, which are commonly oriented in the westerst direction (e.g. Girod and others, 2017; Dussaillant and others, 2019), and can be attributed to satellite attitude oscillations along-track (jitter) (Deschamps-Berger and others, 2020). Similar undulations have been previously found in other satellite products as well, including WorldView DEMs (Shean and others, 2016; Bessette-Kirton and others, 2018). Based on the above, for the elevation change dh/dt grids involving ASTER and Pléiades data, we performed an along-track correction to remove the abovementioned low-frequency undulation pattern. This was done by fitting a spline to the residual off glaciers in the along-track direction and then removing this average residual in the entire dh/dt grids (e.g. Gardelle and others, 2013).

Examination of dh/dt grids including the DEM$_{1962}$ reveals spatially dependent errors, spatially more restricted and more randomly distributed in comparison with the ASTER and Pléiades undulations mentioned above. This is likely because the contact prints were exposed to deformations from storage, and they were

scanned with a non-photogrammetric scanner. For refinement of these results, the photogrammetric DEM$_{1962}$ was cropped around each glacier unit with a 1 km buffer. The clipped DEM was co-registered to the Pléiades 2020 DEM using the Nuth and Kääb (2011) co-registration algorithm implemented by Shean and others (2016). We then accounted for non-linear deformations, applying a fifth-degree polynomial fit from the stable areas in the buffer.

Overall, the individual dh/dt elevation change grids show a very small amount of data gaps on glacier area (~6% maximum of data voids in the 2003–12 grid). Only DEM$_{1962}$ missed 14% of the uppermost part of Turbio I Glacier. For void filling, we followed the local hypsometric approach of McNabb and others (2019), calculating the mean elevation change on elevation bins of 50 m.

For a given time period and glacier unit, the total volumetric change Δv (in m³) is calculated as the product of the average elevation change Δh (m) and the initial glacier area S$_i$ (m²):

$$ \Delta v = \Delta h \times S_i \quad (1) $$

Then, the glacier specific geodetic mass balance Δm (m w.e. a⁻¹) can be calculated as:

$$ \Delta m = \frac{\Delta v \times \rho}{(S_i + S_f)/2} \times t^{-1} \quad (2) $$

ρ being the density conversion factor, S$_f$ the final glacier area and t the time interval in years. For all the calculations, we considered a mean density of 850 kg m⁻³ (Huss, 2013).

We determined the uncertainty associated with the volumetric mass balance as the quadratic sum of the volumetric uncertainties on glacier area and mean elevation change as:

$$ \sigma \Delta v = \sqrt{(S_i \times \sigma_{dh})^2 + (dh \times \sigma_{S_i})^2} \quad (3) $$

where σ$_{dh}$ is the uncertainty on the rate of elevation change and σ$_{S_i}$ is the area uncertainty. The uncertainty on the mean elevation change σ$_{dh}$ is calculated from the patch method (Berthier and others, 2016). This approach aims to empirically determine the uncertainty associated with the mean elevation change dh by sampling patches (or tiles) of the stable terrain of various sizes, in order to constrain the decay of the error with the averaging area (see Supplementary material of Miles and others, 2018; Wagnon and others, 2021 for details).

Finally, the uncertainty on the area-averaged mass balance σ$_m$ for a given glacier is again calculated as a quadratic sum, considering also the uncertainty σ$_{fΔv}$ of ±60 kg m⁻³ in the density assumption for volume to mass conversion (Huss, 2013):

$$ \frac{\sqrt{(f_{\Delta v} \times \sigma \Delta v)^2 + (\Delta v \times \sigma f_{\Delta v})^2}}{(S_i + S_f)/2} \times t^{-1} \quad (4) $$

Further uncertainty sources (apart from the glacier mapping and density assumption), such as seasonality corrections, may be taken into account in geodetic mass-balance determinations. Seasonality corrections can have a great effect on maritime glaciers with a large annual mass turnover (e.g. Andreassen and others, 2005; Belart and others, 2019). However, the Volcán Domuyo area is located in a fairly continental environment where mass turnover is expected to be smaller than in maritime regions. We further note that our DEMs are all acquired at almost the same time of the year. We therefore have excluded seasonality corrections to our mass-balance estimations.
<table><caption>**Table 2.** Validation of the ASTER, Pléiades and photogrammetric (1984) DEM elevations compared to 14 field-surveyed GCPs</caption><thead><tr><th>DEM</th><th>Mean<br>m</th><th>STDV<br>m</th><th>RMSE<br>m</th><th>SE<br>m</th></tr></thead><tbody><tr><td>Pléiades (DEM<sub>2020</sub>)</td><td>1.6</td><td>1.4</td><td>7.8</td><td>0.6</td></tr><tr><td>Pléiades (DEM<sub>2012</sub>)</td><td>4.3</td><td>1.8</td><td>13.4</td><td>0.7</td></tr><tr><td>ASTER (DEM<sub>2003</sub>)</td><td>-8.7</td><td>4.8</td><td>13.8</td><td>1.8</td></tr><tr><td>Photogrammetric (DEM<sub>1984</sub>)</td><td>5.7</td><td>2.0</td><td>17.3</td><td>0.8</td></tr></tbody></table>

## 3.5 Comparison of DEMs with in situ global navigation satellite system (GNSS) measurements

Validation of all DEMs, except for DEM<sub>1962</sub>, was carried out by comparing field-surveyed check points (CPs) off glacier against each DEM (Fig. S1). This was not possible for DEM<sub>1962</sub> because the buffered coverage (for DEM co-registration purposes) around individual glaciers did not coincide with CPs (see Section 3.4). A total of 14 CPs were surveyed during field campaigns in 2012 and 2016 using Trimble 5700 receivers and spanning 2400–3800 m a.s.l. in elevation. CP coordinates were differentially post-processed using data from the Chos Malal GNSS continuous station (37°22′ S–70°16′ W, 875 m, ellips.). The mean horizontal and vertical precision of these CPs were 0.02 and 0.06 m, respectively.

First, we determined the offset of the most recent DEM<sub>2020</sub> relative to the CPs. This yielded a mean XY shift of 6.8 m ($\sigma$ = 0.9 m) and a mean vertical (Z) shift of 1.6 m ($\sigma$ = 1.4 m). Then, after co-registration of all derived orthoimages and DEMs with respect to the DEM<sub>2020</sub>, we applied the systematic shift above and compared the Z coordinate of each DEM against the CPs (Table 2). Particularly, in the cases of DEM<sub>2020</sub> and DEM<sub>2012</sub>, the elevation differences are along the lines of other Pléiades-derived DEMs with no use of GCPs, that were tested on stable (off-glacier) terrain over several evaluation sites worldwide (Berthier and others, 2014; Błaszczyk and others, 2019).

## 3.6 Analysis of climate records

We analysed summer (October–March) air temperatures and winter (April–September) precipitation records available from the Pampa de Chacaico meteorological station (36°28′56.4″ S, 70°36′9.3″ W, 2580 m a.s.l.), located 17 km north of Domuyo. Monthly data from this station were available between 1997 and 2010. These records closely reflect the climate conditions at

Domuyo, as shown when monthly temperature and precipitation records from Pampa de Chacaico are compared with ERA5 reanalysis data using the KNMI climate explorer platform (www.climexp.knmi.nl). Both the temperature and precipitation records from Pampa de Chacaico depict a far-reaching spatial correlation field with $r > 0.6$ over the Central and Northern Patagonian Andes between ~31° S and 43° S (Fig. 2). We extended the monthly time series of air temperature and precipitation at Pampa de Chacaico using ERA5-Land reanalysis data for the period between 1981 and 2020. ERA5-Land is an enhanced dataset derived from ERA5 for the land component (Muñoz-Sabater and others, 2021). We corrected the ERA5-Land data at the location of Pampa del Chacaico through a linear regression, using the meteorological observations available for the 1997–2010 period and ERA5-Land outputs for the same period. Regression models are a simple, appropriate and straightforward method to apply for monthly temperature considering the linearity of the relationship and the normal distribution (Wilby and others, 2004). An RMSE of 0.9°C was obtained after the correction.

The regression parameters were used to bias-adjust the ERA5-Land time series for the entirety of the analysed period (1981–2020). We then calculated the anomalies of the mean summer temperature and total winter precipitation for each year using the 1981–2020 period as a baseline. Lastly, to determine the relation between the glacier balance and regional climate variability, we compared the glacier-wide mass balance with the seasonal variations of temperature and precipitation, averaging the climate records for analogous sub-periods in which the geodetic mass balance was considered. Since ERA5-Land data are available from 1981, the initial 1962–84 sub-period could not be included in this particular analysis.

Using data from Paso de los Indios gauge station (38°31′ S, 69°24′ W; 498 m a.s.l.; not visible in Fig. 1b), we assessed the inter-annual variations of Río Neuquén streamflow. Annual runoff was calculated for the 1903–2019 period (116 years), considering the hydrological year from April to March of the following calendar year. We then evaluated the streamflow variability in relation to the observed geodetic mass balances for similar sub-periods: 1962–83 (21 years), 1984–2002 (18 years), 2003–10 (7 years) and 2011–19 (8 years).

## 3.7 SMB modelling for Volcán Domuyo glaciers

To determine the extent to which the observed geodetic mass balance is linked to climate forcing, we reconstructed the SMB of the

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Location of the Pampa de Chacaico (PC) meteorological station used in this study and 1981–2020 spatial correlation field between monthly anomalies of (a) temperature and (b) precipitation, and the corresponding anomalies in ERA5.</figcaption></figure>
Volcán Domuyo glaciers between the hydrological years 1981/82 and 2019/20, using the corrected ERA5-Land air temperature and precipitation data (previous section) and a minimal glacier mass balance model (MGM). This model was originally presented by Marzeion and others (2012) for glaciers in the alps and has been subsequently applied to reconstruct the annual SMB of Echaurren Norte Glacier (Masiokas and others, 2016) and El Morado Glacier in central Chile (Farias-Barahona and others, 2020b). The model is defined as:

$$MB = \sum_{i=1}^{12} (\alpha P_i - \mu(\max(0, T_i - T_{\text{melt}}))) \quad (5)$$

where MB is the annual specific mass balance. $P_i$ is the monthly precipitation and $\alpha$ is a scaling parameter for precipitation, $T_i$ is the mean monthly air temperature at the front of the glacier, $T_{\text{melt}}$ is the monthly mean temperature above which melt occurs and $\mu$ is the factor to calculate the melt (in mm °C). The maximum operator ensures that melting occurs during months with mean air temperature above $T_{\text{melt}}$. Since no mass-balance observations exist for the Volcán Domuyo glaciers, here we used the same parameters estimated by Masiokas and others (2016) in Echaurren Norte Glacier (~340 km away from Domuyo, the closest glacier where the MGM has been previously used). Hence, we considered a mean $\alpha = 3.9$, $\mu = 90.1 \text{ mm °C}^{-1}$, and we prescribed $T_{\text{melt}} = 0^{\circ}\text{C}$. Although the climate setting between both locations attains a correlation coefficient of over 0.6 (Fig. 2) suggesting similar climate controls, model parameters could vary considerably between glaciers depending on their characteristics and local-scale climate conditions. The uncertainty of the model is relatively high considering the lack of available glaciological observations to calibrate the model. To account for some of the uncertainty of the MGM, we assumed a $\pm 1$ value for the $\alpha$ parameter. Marzeion and others (2012) found a great range of values for the $\alpha$ parameter in the alps. Their mean value of 2 is lower than the $\alpha = 3.9$ value estimated by Masiokas and others (2016), which we also use here.

We applied the model to Turbio I, Chadileo I, Chadileo II and Chadileo IV glaciers (Fig. 1), the four largest glaciers around Volcán Domuyo. Turbio I and Chadileo II have a continuous debris-covered layer at the tongue, whereas Chadileo I and Chadileo IV also present areas with debris-cover, though spatially discontinuous. In view of the continuous debris-covered layer on the tongues that can reduce the mass-balance sensitivity (e.g. Anderson and Mackintosh, 2012), we also applied the model using the elevation at which bare-ice prevails over debris-covered ice in the cases of Turbio I and Chadileo II. Based on the above, we discarded from the analysis Covunco Glacier (also one of the largest glaciers), as the debris-cover is predominant across its surface.

To extrapolate the air temperature data from Pampa de Chacaico to the elevation of the front and bare-ice area of each selected glacier (input to the model), we derived monthly lapse rates using the available records from stations in the vicinity (<37 km) of Volcán Domuyo during the time interval between 1997 and 2010. These stations were Varvarco (1180 m a.s.l.), Nehuén (1225 m a.s.l.), Cajón de los Chenques (1533 m a.s.l.) and Pampa de Chacaico (2580 m a.s.l.) (Fig. 1b). This way, a catchment-scale temperature lapse rate was generated, and a mean of $4.5 \pm 1.6^{\circ}\text{C km}^{-1}$ was obtained. Steeper lapse rates were computed during summer months, reaching values close to the environmental lapse rate (6.5°C km$^{-1}$), whereas shallower lapse rates were computed for the winter months (Fig. S2). The extended air temperature time series (April 1981–March 2020) was extrapolated using the monthly lapse rates to the front of

the selected glaciers, which was estimated using the glacier outlines of 1984, 2003, 2012 and 2020. The elevation of the glacier fronts (i.e. minimum elevation of each glacier) was extracted using the 2020 Pléiades DEM, and the glacier fronts were linearly interpolated between the initial and final years of survey periods analogous to the geodetic mass balance with the inventoried glacier polygons.

## 4. Results

### 4.1 Glacier inventory and area changes between 1962 and 2020

When computing the glacier changes for the full sampled time period (1962–2020), we considered the ten largest glaciers (Table 3). Accordingly, the glacier area loss between 1962 and 2020 was $4.7 \text{ km}^2$ (20.8%) at a rate of $-0.08 \text{ km}^2 \text{ a}^{-1}$ (0.36% a$^{-1}$). The glacier with the highest relative area loss was Chadileo III, losing half of its 1962 area. Turbio I was the glacier that lost the greatest area ($-1.3 \text{ km}^2$; ~27%) between 1962 and 2020, owed to the fast retreat of its narrow debris-covered tongue. Interestingly, the rate of glacier shrinkage was limited for almost four decades (the glacier area changes between 1962 and 2003 was negligible), but picked up pace after 2003, when it attained $0.8 \text{ a}^{-1}$.

Based on the 1984 aerial photographs, we identified and mapped 46 glaciers larger than $0.01 \text{ km}^2$ that represented a total area of $29.0 \pm 0.2 \text{ km}^2$. The largest inventoried glacier is Chadileo II ($5.1 \pm 0.02 \text{ km}^2$). By 2020, five glacier units mapped in the 1984 orthophotos had disappeared, and the total glacier area had decreased to $20.5 \pm 0.2 \text{ km}^2$. The total glacier area loss of $8.5 \text{ km}^2$ corresponds to 29.3% of the 1984 glacier area, and the area loss rate between 1984 and 2020 was $\sim 0.25 \text{ km}^2 \text{ a}^{-1}$ ($0.8 \text{ a}^{-1}$).

### 4.2 Geodetic mass balance

We determined the area-averaged geodetic mass balance of 46 glaciers for the 1984–2020 period, while we also look at 11 individual glaciers in greater temporal detail. From these 11 glaciers, we additionally determined the mass balance of ten (owing to Manchana Covunco Glacier missing from the 1962 aerial photographs) between 1962 and 2020 (Fig. 3 and Table 4).

Overall, the cumulative mass balance was $-15.12 \pm 3.60 \text{ m w.e.} (-0.42 \pm 0.10 \text{ w.e. a}^{-1})$ during 1984–2020 (46 glaciers), and the total volume loss was $12.2 \pm 2.7 \times 10^6 \text{ m}^3 \text{ a}^{-1}$. During this period, Turbio I Glacier ($-0.55 \pm 0.11 \text{ m w.e. a}^{-1}$) had the most negative mass balance, whereas Covunco ($-0.23 \pm 0.09 \text{ m w.e. a}^{-1}$) had the least negative mass loss rate. As mentioned in Section 3.4, we performed a mass-balance triangulation test to evaluate the residuals between the sum of the individual time interval $dh/dt$ values and the accumulated 1984–2020 $dh/dt$ balance. This yielded a mean difference of $+0.01 \text{ m w.e. a}^{-1}$ for the total study

<table><caption>Table 3. Glacier inventory and area changes</caption><thead><tr><th rowspan="2">Year</th><th rowspan="2">Number of glaciers</th><th rowspan="2">Glacier area km<sup>2</sup></th><th colspan="2">Area change %</th></tr><tr><th>Initial time = 1984</th><th>Initial time = 1962</th></tr></thead><tbody><tr><td>1962</td><td>10</td><td>(22.6 ± 0.2)</td><td>-</td><td>-</td></tr><tr><td>1984</td><td>46</td><td>29 ± 0.2 (23.1 ± 0.2)</td><td>-</td><td>+2.2</td></tr><tr><td>2003</td><td>46</td><td>27.9 ± 7.0 (22.6 ± 4.5)</td><td>-3.8</td><td>0</td></tr><tr><td>2012</td><td>46</td><td>23.4 ± 0.2 (20.1 ± 0.2)</td><td>-19.3</td><td>-11.1</td></tr><tr><td>2020</td><td>41</td><td>20.5 ± 0.2 (17.9 ± 0.2)</td><td>-29.3</td><td>-20.8</td></tr></tbody></table>

Values in brackets indicate glacier changes for the ten glaciers mapped for 1962. Note also the disappearance of five small glaciers in the 2020 inventory
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Annual glacier elevation changes between 1962 and 2020 (a) and 1984 and 2020 (b), and for four sub-periods between 1962 and 2020 (c-f). For comparability the same colour scale is used for all panels. Note the much noisier signal on stable terrain in the 1984–2003 and 2003–12 panels (where the ASTER DEM comes to play), compared to the 1962–84 and 2012–20 panels (where DEMs are derived from repeated Pleiades or aerial imagery).</figcaption></figure>

area, while for individual glaciers residuals ranged between ±0.03 m w.e. a⁻¹.

In spite of the pronounced general mass loss throughout the 1984–2020 period and a fairly homogeneous mass-balance response among individual glaciers within individual sub-periods

<table><caption>**Table 4.** Accumulated 1984–2020 and 1962–2020 volume change (Δv) and mass balance (Δm) for 11 specific glaciers and total glacier area</caption><thead><tr><th rowspan="2">Glacier</th><th colspan="2">1984–2020</th><th colspan="2">1962–2020</th></tr><tr><th>Δv (10<sup>6</sup> m<sup>3</sup> a<sup>-1</sup>)</th><th>Δm (10<sup>6</sup> m<sup>3</sup> a<sup>-1</sup>)</th><th>Δv (10<sup>6</sup> m<sup>3</sup> a<sup>-1</sup>)</th><th>Δm (m w.e. a<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Domuyo</td><td>-0.9 ± 0.2</td><td>-0.5 ± 0.11</td><td>-0.3 ± 0.2</td><td>-0.15 ± 0.12</td></tr><tr><td>Norte I</td><td></td><td></td><td></td><td></td></tr><tr><td>Domuyo</td><td>-0.3 ± 0.1</td><td>-0.52 ± 0.12</td><td>-0.1 ± 0.1</td><td>-0.20 ± 0.15</td></tr><tr><td>Norte II</td><td></td><td></td><td></td><td></td></tr><tr><td>Chadileo I</td><td>-2.1 ± 0.4</td><td>-0.46 ± 0.10</td><td>-0.7 ± 0.5</td><td>-0.15 ± 0.11</td></tr><tr><td>Chadileo II</td><td>-2.1 ± 0.5</td><td>-0.38 ± 0.10</td><td>-0.8 ± 0.5</td><td>-0.14 ± 0.10</td></tr><tr><td>Chadileo III</td><td>-0.3 ± 0.1</td><td>-0.46 ± 0.13</td><td>-0.1 ± 0.1</td><td>-0.11 ± 0.16</td></tr><tr><td>Chadileo IV</td><td>-0.6 ± 0.2</td><td>-0.32 ± 0.10</td><td>-0.3 ± 0.2</td><td>-0.16 ± 0.12</td></tr><tr><td>Covunco</td><td>-0.6 ± 0.2</td><td>-0.23 ± 0.09</td><td>-0.1 ± 0.3</td><td>~0 ± 0.10</td></tr><tr><td>Manchana</td><td>-0.8 ± 0.2</td><td>-0.51 ± 0.12</td><td>-</td><td>-</td></tr><tr><td>Covunco</td><td></td><td></td><td></td><td></td></tr><tr><td>Turbio I</td><td>-2.7 ± 0.5</td><td>-0.55 ± 0.11</td><td>-1.1 ± 0.5</td><td>-0.23 ± 0.11</td></tr><tr><td>Turbio II</td><td>-0.3 ± 0.1</td><td>-0.49 ± 0.12</td><td>-0.1 ± 0.1</td><td>-0.19 ± 0.14</td></tr><tr><td>Turbio III</td><td>-0.3 ± 0.1</td><td>-0.41 ± 0.13</td><td>-0.1 ± 0.1</td><td>-0.13 ± 0.15</td></tr><tr><td>Total glacier area</td><td>-12.2 ± 2.7</td><td>-0.42 ± 0.10</td><td>-</td><td>-</td></tr><tr><td>Total glacier area<sup>a</sup></td><td>-10.4 ± 2.2</td><td>-0.43 ± 0.10</td><td>-3.6 ± 2.1</td><td>-0.15 ± 0.09</td></tr></tbody></table>

(Figs 3, 4), mass-balance conditions varied significantly among the three sub-periods.

From 1984 to 2003, the studied glaciers were close to an equilibrium state (−0.06 ± 0.09 m w.e. a⁻¹). Turbio II (−0.27 ± 0.20 m w.e. a⁻¹) and Turbio III (−0.23 ± 0.19 m w.e. a⁻¹) were the glaciers with the most negative mass balance over this study period. From 2003 to 2012, Domuyo glaciers shifted to strongly negative mass balances (−0.54 ± 0.34 m w.e. a⁻¹), with the greatest mass loss measured for Domuyo Norte I (−0.67 ± 0.38 m w.e. a⁻¹) and Turbio I (−0.66 ± 0.31 m w.e. a⁻¹). The most negative mass-balance conditions throughout the study period were attained between 2012 and 2020, with a strongly negative glacier-wide mass balance of −0.99 ± 0.19 m w.e. a⁻¹. During this period, Domuyo Norte II (−1.46 ± 0.33 m w.e. a⁻¹) had the largest mass loss rate.

Looking at the ten glaciers sampled from 1962 to 1984, we found an overall moderately positive mass balance of +0.28 ± 0.13 m w.e. a⁻¹. Domuyo Norte I (+0.37 ± 0.20 m w.e. a⁻¹) had the most positive mass balance for this time interval. The positive mass balances from 1962 to 1984 and the close to equilibrium conditions between 1984 and 2003 partly counterbalanced the losses after 2003 (Table 5). Consequently, the overall mass balance for these ten glaciers for the full 1962–2020 period was only slightly negative (−0.15 ± 0.09 m w.e. a⁻¹), and the mass balance among glaciers varied to a limited extent between 0 and −0.23 m w.e. a⁻¹ (Table 4).

The mean elevation change averaged over 50 m elevation bands (Fig. 5) reveals that for the 2003–12 and 2012–20 intervals nearly all of the glacier area thinned, including the uppermost

<footnote><sup>a</sup>Refers to the listed glaciers in the table (minus Manchana Covunco) only.</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Glacier specific area reduction (a) and mass balance (b) per time interval. Also plotted in panel (b) are the geodetic mass balances of Echaurren Norte and El Morado glaciers by Farias-Barahona and others (2019, 2020b).</figcaption></figure>

reaches. The 1962–84 curve shows a strongly positive gradient between 2800 and 3000 m (shifting from −1.5 to 1.3 m a⁻¹), and remains under moderately positive dh/dt conditions (<0.5 m a⁻¹) higher up-glacier, in line with the overall positive mass balance. Elevation gains are visible below 3100 m for the 1984–2003 and 2003–12 periods respectively, which are a consequence of the slight advances of the Covunco and Chadileo II tongues around those time intervals (Figs 3d, e). The rates of surface elevation change were similar between 1962–84 and 1984–2003 above 3700 m, and generally show a much gentler gradient in comparison with the 2012–20 curve, where the gradient is much steeper.

## 4.3 Climate and hydrological trends

During the 1981–2019 period, we observed an increase in summer air temperatures (0.13°C decade⁻¹) and a decrease in winter precipitation (−50 mm decade⁻¹). A correlation analysis over the different sub-periods showed high *r* values between the glacier-wide mass balance, summer temperatures (*r* = 0.76; *p* = 0.44), and winter precipitation (*r* = 0.91; *p* = 0.27) for the 1984–2020 time interval (Fig. 6). The *p*-values of these correlations indicate that they are statistically not significant, which is due to the low number of mass-balance observations (only three periods). Consequently, these findings must be approached with care, and motivate the implementation of the MGM to better understand

the climate forcing on glaciers. The 1984–2003 period, when glaciers were close to an equilibrium state (−0.06 ± 0.09 m w.e. a⁻¹), coincides with a period of very slightly positive precipitation anomalies (+2%) and near-zero temperature anomalies. Between 2003 and 2012, the glacier-wide mass balance became more negative (−0.54 ± 0.34 m w.e. a⁻¹), concomitant with near-zero temperature and precipitation anomalies. The final 2012–20 interval represents the peak in negative mass balance for the last 58 years (−0.99 ± 0.19 m w.e. a⁻¹), and was concurrent with a mean summer temperature and winter precipitation anomalies of +0.23°C and −0.1, respectively.

The mean streamflow of the Rio Neuquén at Paso de los Indios gauge station was 299.24 m³ s⁻¹ over the long-term 1903–2019 period. We found a slightly negative trend of −0.5 m³ s⁻¹ a⁻¹ over this interval (Fig. 8), representing an annual streamflow decrease of 58.0 m³ s⁻¹ during the last 116 years. In the context of this long-term trend of streamflow, shorter periods with considerably higher or lower flows occurred. Considering the sub-periods analysed in the geodetic mass-balance assessment, the two earlier intervals (1962–84 and 1984–2003) showed a mean streamflow similar to the full 1903–2019 period (302.4 and 296.8 m³ s⁻¹, respectively). In addition, throughout these two periods, we observed streamflow increases of ~140.0 m³ s⁻¹ (i.e. 6.4 m³ s⁻¹ a⁻¹) and 47.0 m³ s⁻¹ (i.e. 2.5 m³ s⁻¹ a⁻¹), respectively. The 1984–2003 interval depicts a slight streamflow increase
<table><caption>**Table 5.** Time-interval volume change (Δv) and mass balance (Δm) for 11 specific glaciers and total glacier area</caption><thead><tr><th rowspan="2">Glacier</th><th colspan="2">1962-84<sup>a</sup></th><th colspan="2">1984-2003</th><th colspan="2">2003-12</th><th colspan="2">2012-20</th></tr><tr><th>Δv (10<sup>6</sup> m<sup>3</sup> a<sup>-1</sup>)</th><th>Δm (m w.e. a<sup>-1</sup>)</th><th>Δv (10<sup>6</sup> m<sup>3</sup> a<sup>-1</sup>)</th><th>Δm (m w.e. a<sup>-1</sup>)</th><th>Δv (10<sup>6</sup> m<sup>3</sup> a<sup>-1</sup>)</th><th>Δm (m w.e. a<sup>-1</sup>)</th><th>Δv (10<sup>6</sup> m<sup>3</sup> a<sup>-1</sup>)</th><th>Δm (m w.e. a<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Domuyo Norte I</td><td>+0.7 ± 0.4</td><td>+0.37 ± 0.20</td><td>-0.2 ± 0.3</td><td>-0.08 ± 0.16</td><td>-1.2 ± 0.7</td><td>-0.67 ± 0.38</td><td>-1.6 ± 0.4</td><td>-1.03 ± 0.27</td></tr><tr><td>Domuyo Norte II</td><td>+0.2 ± 0.2</td><td>+0.22 ± 0.24</td><td>0.1 ± 0.1</td><td>0.19 ± 0.20</td><td>-0.4 ± 0.3</td><td>-0.61 ± 0.45</td><td>-0.8 ± 0.2</td><td>-1.46 ± 0.33</td></tr><tr><td>Chadileo I</td><td>+1.7 ± 0.8</td><td>+0.33 ± 0.17</td><td>0.1 ± 0.7</td><td>0.02 ± 0.13</td><td>-2.9 ± 1.5</td><td>-0.59 ± 0.30</td><td>-5 ± 1</td><td>-1.15 ± 0.24</td></tr><tr><td>Chadileo II</td><td>+1.6 ± 1</td><td>+0.27 ± 0.16</td><td>0.4 ± 0.8</td><td>0.07 ± 0.13</td><td>-3 ± 1.6</td><td>-0.51 ± 0.28</td><td>-5.3 ± 1.1</td><td>-0.98 ± 0.22</td></tr><tr><td>Chadileo III</td><td>+0.3 ± 0.2</td><td>+0.35 ± 0.23</td><td>-0.1 ± 0.2</td><td>-0.09 ± 0.19</td><td>-0.3 ± 0.4</td><td>-0.35 ± 0.45</td><td>-0.7 ± 0.2</td><td>-1.29 ± 0.34</td></tr><tr><td>Chadileo IV</td><td>+0.2 ± 0.4</td><td>+0.09 ± 0.20</td><td>0.3 ± 0.4</td><td>0.13 ± 0.16</td><td>-0.7 ± 0.7</td><td>-0.37 ± 0.35</td><td>-2.2 ± 0.4</td><td>-1.23 ± 0.26</td></tr><tr><td>Covunco</td><td>+0.9 ± 0.5</td><td>+0.34 ± 0.19</td><td>&lt;-0.1 ± 0.4</td><td>-0.03 ± 0.15</td><td>-0.9 ± 0.8</td><td>-0.35 ± 0.32</td><td>-1.2 ± 0.6</td><td>-0.47 ± 0.23</td></tr><tr><td>Manchana Covunco</td><td>-</td><td>-</td><td>-0.4 ± 0.3</td><td>-0.2 ± 0.17</td><td>-0.9 ± 0.6</td><td>-0.57 ± 0.39</td><td>-1.5 ± 0.4</td><td>-1.08 ± 0.27</td></tr><tr><td>Turbio I</td><td>+1.3 ± 0.9</td><td>+0.23 ± 0.17</td><td>-0.9 ± 0.7</td><td>-0.16 ± 0.13</td><td>-3.8 ± 1.6</td><td>-0.66 ± 0.31</td><td>-5 ± 1</td><td>-1.13 ± 0.23</td></tr><tr><td>Turbio II</td><td>+0.2 ± 0.2</td><td>+0.28 ± 0.24</td><td>-0.3 ± 0.1</td><td>-0.27 ± 0.20</td><td>-0.4 ± 0.3</td><td>-0.53 ± 0.48</td><td>-0.6 ± 0.2</td><td>-1.34 ± 0.30</td></tr><tr><td>Turbio III</td><td>+0.3 ± 0.2</td><td>+0.30 ± 0.22</td><td>-0.1 ± 0.2</td><td>-0.23 ± 0.19</td><td>-0.3 ± 0.4</td><td>-0.54 ± 0.52</td><td>-0.6 ± 0.2</td><td>-1.13 ± 0.29</td></tr><tr><td>Total glacier area</td><td>-</td><td>-</td><td>-2 ± 3.0</td><td>-0.06 ± 0.09</td><td>-16.2 ± 7.2</td><td>-0.54 ± 0.34</td><td>-25.5 ± 3.9</td><td>-0.99 ± 0.19</td></tr></tbody></table>

compared to the former interval, and is still heavily influenced by the larger streamflow during the last 3 years (2000–02, see Fig. 7b). In contrast to the 1962–2003 period, the two following intervals (2003–12 and 2012–19) show lower mean streamflow values compared to the overall mean (280.0 and 179.4 m³ s⁻¹), and strongly decreasing trends of −13.4 and −4.6 m³ s⁻¹ a⁻¹, respectively. Again, for all these four time intervals, the correlation coefficient between the glacier-wide mass balance and mean streamflow are high (r ~ 0.9), although not significant (p = 0.12) at the 95% confidence level due to the low number of periods compared (n = 4).

## 4.4 Modelled SMB and comparison against the geodetic mass results

We applied the minimum glacier mass-balance model to Chadileo I, Chadileo IV, Turbio I and Chadileo II glaciers (Table 6, Fig. 8). The mass balances of Chadileo II and Turbio I were modelled using two different elevation inputs: the elevation of the glacier front (including the debris-covered area), and the estimated elevation at which bare-ice prevails over debris-covered ice. The modelled glaciers show a predominance of negative SMB. The mean SMB shows that just 6 hydrological years are positive (i.e. 1982, 1984, 1997, 2001,

2005 and 2006). Notably, the mass balances for the last 14 years are all negative. As expected, modelled SMB values are more negative when the glacier front (including debris-covered areas) is assumed as input. In this scenario, Chadileo I and IV show some positive SMB years (up to 1 m w.e. a⁻¹), while the SMB of Turbio I and Chadileo II remains negative throughout all of the time series. When the elevation of the bare-ice area is considered, both Turbio I and Chadileo II show a larger number of positive SMB years.

We compared the results of the modelled SMB and the geodetic assessment for the two glaciers where bare-ice prevails (Chadileo I and IV) and for the two other glaciers where debris-covered ice is more abundant (Chadileo II and Turbio I). For the two former glaciers, we found a good agreement between the geodetic and modelled values for some individual periods, although the model can either under- or overestimate the geodetic mass balance by up to 0.6 m w.e. a⁻¹ for other periods (Fig. 9). In the case of the two latter glaciers, there was a significant underestimation (up to 1.2 m w.e. a⁻¹) of the SMB (more negative) when the model uses the front of the glacier as input. Correspondence between the modelled and geodetic results improves when the elevation of the bare-ice area is used as input to the SMB model, yet differences of ~1 m w.e. a⁻¹ persisted for certain periods.

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Hypsometric curve of surface elevation changes over the Domuyo glaciers. The grey bars represent the glacier area distribution on 50 m elevation bins for the full 1984 glacier area (right axis).</figcaption></figure>

<footnote><sup>a</sup>Refers to the listed glaciers in the table (minus Manchana Covunco) only.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Time-interval evolution of geodetic glacier mass balance compared to (a) summer temperature and (b) winter precipitation. In the upper (a) panel, the mass-balance values on the right axis have been reversed for a better interpretation.</figcaption></figure>

# 5. Discussion

## 5.1 Glacier area changes at Volcán Domuyo and nearby glacierized areas

Previously, glacier area loss rates at Volcán Domuyo were estimated to be 1.4% a⁻¹ between 1990 and 2009 by Falaschi and others (2016), who used medium-to-high spatial resolution Landsat™ (30 m) and ALOS PRISM (2.5 m) images. With an extended analysis period and improved glacier outlines thanks to higher spatial resolution imagery, we can separate an initial period (1962–84) of slight area increase (+2.2%, i.e. 0.1% a⁻¹ on average), and a later period (1984–2020) when glacier area loss rates increased to ~0.6% a⁻¹. Overall, the glacier area change rate was -0.36% a⁻¹ between 1962 and 2020. Interestingly, the scatter in area change percentage is relatively low between 1962 and 2003, but becomes much larger afterwards (Fig. 4a). Reinthaler and others (2019) investigated area changes of ice-caped volcanoes along the South American Andes and report a 1.5% a⁻¹ reduction for the 1986–2015 period for the Domuyo district. For a similar period (1984–2012), we report a much lower reduction rate of 0.6% a⁻¹, and conclude that differences are owed to the fact that Reinthaler and others (2019) refer mostly to the debris-free ice (Fig. S3). Our 0.36% a⁻¹ recession rate during 1962–2020 is similar to the ~0.5% a⁻¹ reported by Malmros and others (2016) for the Central Andes of Chile between 1955 and 2013/14, though highly variable rates among individual glaciers

have also been observed by Barcaza and others (2017) in this region.

## 5.2 The 1962–2020 Volcán Domuyo Glacier mass-balance record in a broader geographic perspective

### 5.2.1 Mass balance pre-2000

In the Central Andes of Chile, glacier volume and mass changes of the Maipo catchment (33–34° S, Fig. 1c) were investigated by Farias-Barahona and others (2020a) using topographic maps and the SRTM DEM over the 1955–2000 period. The geodetic mass balance for the entire catchment (361.4 km² of glacier area) was -0.12 ± 0.06 m w.e. a⁻¹. At the sub-catchment scale, however, these authors found an in-balance to slightly positive mass-balance trend for the southernmost Upper Maipo and Volcán sub-basins, in coincidence with our results for the 1962–2003 period. Incidentally, these geographical differences in glacier mass balance, which show more negative (positive) mass balance north (south) of 34° S, have been previously linked to lower (higher) precipitation over two hydroclimatic regimes at this region of the Andes (González-Reyes and others, 2017; Ayala and others, 2020).

Farias-Barahona and others (2019, 2020b) also assessed individual geodetic mass balances of the World Glacier Monitoring Service (WGMS) reference Echaurren Norte (33°34′ S–70°7′ W; 0.17 km²) and El Morado (33°44′ S–70°3′ W; 1.05 km²) glaciers
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Streamflow time series from the Paso de los Indios in the Rio Neuquén catchment for the 1903–2019 period (a) and detailed time-interval streamflow variation in relation to the mean geodetic mass balance in Volcán Domuyo (b).</figcaption></figure>

within the Maipo catchment. In the case of Echaurren Norte, an average mass balance of –0.85 ± 0.08 m w.e. a⁻¹ was obtained for the 1955–2000 period, again based on topographic maps and

SRTM. On the other hand, El Morado lost mass at a rate of –0.14 ± 0.19 m w.e. a⁻¹ between 1955 and 1996, as revealed from photogrammetric surveys. No mass-balance data before 2000 for

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Time series of the reconstructed annual SMB computed using the MGM. Lines within the grey area show the values for the individual glaciers. The blue bars represent the mean SMB of all four glaciers and the six experimental runs, and the black bars correspond to the mass balance estimated with the glaciological method in Echaurren Norte Glacier (33° S) (WGMS, 2021). Values in parentheses are the correlation coefficients between each time series and the glaciological mass balance of Echaurren Norte.</figcaption></figure>
<table><caption>**Table 6.** Modelled glacier mass balance (1984–2020)</caption><thead><tr><th rowspan="3">Glacier</th><th colspan="4">SMB (m w.e. a<sup>-1</sup>)</th></tr><tr><th colspan="2">1984–2003</th><th colspan="2">2003–12</th><th colspan="2">2012–20</th></tr><tr><th>MGM</th><th>Geodetic</th><th>MGM</th><th>Geodetic</th><th>MGM</th><th>Geodetic</th></tr></thead><tbody><tr><td>Chadileo I</td><td>-0.03</td><td>0.02 ± 0.13</td><td>-0.05</td><td>-0.59 ± 0.30</td><td>-0.72</td><td>-1.15 ± 0.24</td></tr><tr><td>Chadileo II</td><td>-2.16 (–0.84)</td><td>0.07 ± 0.13</td><td>-2.12 (–0.66)</td><td>-0.51 ± 0.28</td><td>-2.87 (–1.32)</td><td>-0.98 ± 0.22</td></tr><tr><td>Chadileo IV</td><td>-0.50</td><td>0.13 ± 0.16</td><td>-0.45</td><td>-0.37 ± 0.35</td><td>-1.17</td><td>-1.23 ± 0.26</td></tr><tr><td>Turbio I</td><td>-1.39 (0.02)</td><td>-0.16 ± 0.13</td><td>-0.90 (0.28)</td><td>-0.66 ± 0.31</td><td>-1.50 (–0.16)</td><td>-1.13 ± 0.23</td></tr></tbody></table>

Values in brackets indicate mass-balance values referred to the glacier bare-ice area only.

the Northern Patagonian Andes are available, so that no parallel can be drawn with Volcán Domuyo.

Beyond the geodetic mass-balance determinations available for the Central Andes, the Chilean Departamento General de Aguas has carried out in situ glaciological investigations on Echaurren Norte since 1975. These surveys yielded a mean mass balance of –0.03 m w.e. a⁻¹ between 1975 and 1984, and –0.35 m w.e. a⁻¹ between 1984 and 2003 (WGMS, 2021). While some survey years are missing from our 1962–84 period, these data nevertheless confirm a shift towards more negative mass balance in the Central Andes from this initial period to the follow-up 1984–2003 interval. On the contrary, the reconstructed Echaurren Norte 1962–84 (–0.53 m w.e. a⁻¹) and 1984–2003 (–0.40 m w.e. a⁻¹) average mass

balance (Masiokas and others, 2012) do not replicate a similar shift as indicated by the glaciological record. Thus, the positive balance of 1962–84 observed over Volcán Domuyo cannot be consistently and accurately tracked elsewhere in the Central or Northern Patagonian Andes based on geodetic or modelled mass-balance records available so far.

### 5.2.2 Mass balance post-2000

Our results for Volcán Domuyo, showing negative mass-balance conditions between 2003 and 2012, and very strongly negative conditions between 2012 and 2020, are in line with the findings of Farias-Barahona and others, (2019) for Echaurren Norte Glacier (~33.56° S), who found a geodetic mass balance of

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Comparison of the geodetic mass balance and the modelled SMB for the three periods analysed and the four selected glaciers. Error bars of the geodetic mass balance are explained in Section 3.4 and the error bars in the modelled SMB represent the std dev. during the corresponding period.</figcaption></figure>
precipitation at the Pampa de Chacaico site is highly correlated with the regional snow accumulation series based on eight snow-pack records spread along the Central Andes (Masiokas and others, 2012). The overall picture does not suggest a clear winter precipitation trend since 1951, though shorter periods of drier and more humid conditions can be observed. If our initial 1962–84 geodetic mass-balance survey period is analysed within this time series, the snowpack-derived record shows below average winter precipitation between 1962 and 1973, and above average precipitation between 1974 and 1984. Dry conditions were also present from the late 1980s through the early 1990s.

As described above, the mass balance during 2012–20 was strongly negative (−0.99 ± 0.19 m w.e. a⁻¹), and coincides with highly negative precipitation anomalies from 2010 onwards (Fig. 6b). This can be interpreted as the onset of the spatially and temporally severe mega-drought that affects both western and eastern flanks of the Central and Northern Patagonian Andes of Chile and Argentina (Garreaud and others, 2017, 2019; Rivera and others, 2017; Dussaillant and others, 2019).

Further insight is provided by the Río Neuquén streamflow recorded at Paso de los Indios. The catchment’s streamflow has a predominantly pluvio-nival regime and the gauge station is located in the lower part of the catchment (~220 km from Domuyo), meaning that contributions from other sources, beyond the Domuyo glaciers, dominate the streamflow at Paso de los Indios. We correlated the Paso de los Indios streamflow with the much shorter 1974–2009 record of Puente Andacollo gauge station, located in the upper Río Neuquén catchment (~60 km from Domuyo; Fig 1b). As we found a high correlation (r = 0.86; p < 0.001) between these two time series, we put forward that the Paso de los Indios record is representative of the upper-most Rio Neuquén catchment. With these caveats in mind, we observed a close relationship between the geodetic glacier mass balance and the mean streamflow at Paso de los Indios. In particular, the 1962–84 geodetic mass balance (+0.28 m w.e. a⁻¹) is consistent with a marked streamflow increase of ~140.0 m³ s⁻¹ over this interval, and a mean streamflow value of 302.4 m³ s⁻¹, higher than the mean streamflow for the entire 1903–2019 period. On the contrary, the 2003–11 and 2012–19 periods showed highly negative mass balances, concurrent with lower mean values of 280.0 and 179.0 m³ s⁻¹, and pronounced runoff reductions of ~140.0 and ~41.0 m³ s⁻¹ over these intervals, respectively (Fig. 8b). Based on this evidence, it can be argued that glacier melt is not able to compensate for the mega-drought-driven deficit in the precipitation-fed Río Neuquén catchment, because of the limited fraction (<1%, Zalazar and others, 2020) of the catchment covered by glaciers.

## 5.4 Insight from the SMB model in Volcán Domuyo and the benchmark Echaurren Norte Glacier in the Central Andes

The choice of the minimal glacier model of Marzeion and others (2012) lies in the minimal input requirements, using solely monthly precipitation and air temperatures as climate forcing. Although glaciological surveys can be used to calibrate the model, there are no in situ measurements available for Volcán Domuyo glaciers. Calibration to geodetic results has been performed for large-scale glacier models recently, even when few time periods are available (e.g. Zekollari and others, 2019). This, however, is likely to result in equifinality. Due to the lack of calibration of the model parameters, we suggest that absolute values of estimated SMB must be taken with caution.

We observed significant discrepancies in absolute values between the geodetic and modelled SMB, especially in the (partly) debris-covered glaciers (Chadileo II and Turbio I). Also, the modelled SMB shows periods of either more positive or less negative

mass budget compared to the geodetic results, meaning that the model is not fully able to reproduce the evolution of the mass balance of individual glaciers. We interpret that the differences found between the modelled and geodetic determinations may be owed to a number of factors. On the one hand, there are factors associated with the uncertainty of the model parameters (as mentioned in Section 3.7) and also with the use of extrapolated lapse rates derived from low/medium elevations to the high and complex topography of Domuyo glaciers. On the other hand, other intervening factors that are not accounted for in our modelling include (1) the thick debris-cover on some portions of the sampled glaciers in Volcán Domuyo (which, as explained in Section 3.7, decreases glacier sensitivity to climate forcing) and (2) volcanic unrest and surface warming acting at a decadal to pluri-decadal scale (Astort and others, 2019; Lundgren and others, 2020) may come into play in the geodetic determination but not in the modelled SMB. We have nevertheless neglected the possible effect of the high geothermal flux (present in the Domuyo area) on glacier mass balance, as it would not be trivial to assess from an instrumental point of view (see e.g. Trombotto Liaudat and others, 2014).

The most interesting aspect of the modelled annual SMB lies in that the overall interannual variability of the four selected glaciers was similar to the mass balance of Echaurren Norte Glacier, determined by the Departamento General de Aguas of Chile using the glaciological method (WGMS, 2021) (Fig. 8), and located 340 km north of Volcán Domuyo. The correlation coefficients between both the annual time series are ~0.6, and all are statistically significant (p-value <0.05), indicating a strong large-scale climate forcing. Although this is an interesting development (especially considering the distance between Echaurren Norte and Domuyo), we stress here that the SMB results may be biased, and that this is not a calibrated model.

As a final validation of the modelled SMB using the MGM, we compared our results against those provided by ten forcing datasets of the ensemble-based model of Malles and Marzeion (2021) for the four modelled glaciers (Fig. S5). Again, the correlation coefficient between the ensemble mean and our MGM values is high (r = 0.8), indicating that the interannual variability is well represented, though large differences (~1 m w.e. a⁻¹) can still persist in absolute terms.

## 6. Conclusion

In this study we have calculated the 1984–2020 (1962–2020) geodetic mass balance of 46 (10) glaciers in Volcán Domuyo in the northern Patagonian Andes by differencing DEMs from aerial photographs, ASTER and Pléiades satellite stereo images. This assessment, which encompasses almost six decades of glacier mass change, represents the longest mass-balance record for a region where in-depth glaciological studies are few and far between. Our results for the entire 1962–2020 study period show an increase in ice loss rates with time, shifting from an initial positive budget (+0.28 ± 0.13 m w.e. a⁻¹) between 1962 and 1984, to a strongly negative mass balance of −0.99 ± 0.19 m w.e. a⁻¹ from 2012 to 2020. Comparisons with other geodetic determinations for the Domuyo glaciers found in the literature show that the main differences can be linked to the usage of different DEMs, the diverse duration of sampled time periods and the contrasting quality of glacier inventories. The latter also highlights the importance of using highly accurate glacier outlines for geodetic assessments at a more local scale.

The observed mass-balance evolution of the Domuyo glaciers is related to an increase in summer air temperatures (0.13°C decade⁻¹) and a decrease in winter precipitation (−50 mm decade⁻¹) during the 1981–2019 period. In particular, we hypothesize that
the so-called mega-drought of the Central Andes is the main dri-
ver of the unusually high mass loss rates found for the 2012–20
interval.

The variability in mass balance throughout the study period
that we observed using the DEM differencing method can be,
in general terms, effectively tracked among previously published
research for other catchment scale or individual glacier sites
within the Central Andes, especially for the post-2000 period.
More data are nevertheless necessary to give further support to
the positive mass-balance conditions found for the initial 1962–
84 period.

Despite the limitations of the modelled SMB in terms of abso-
lute values, the interannual variability showed similarity to the
glaciological records of Echaurren Norte Glacier, confirming
that the mass-balance evolution of the Domuyo glaciers do bear
resemblance to those located in the Central Andes further
north. With no signs of a decline in the mega-drought in sight,
similar or even higher glacier mass loss rates can be expected in
the near future for Volcán Domuyo in particular and the
Central Andes in general.

**Supplementary material.** The supplementary material for this article can be found at https://doi.org/10.1017/jog.2022.43

**Data.** The newly produced glacier outlines and glacier-specific mass-balance data will be made available in the WGMS and/or GLIMS databases. In the meantime, these datasets are available upon request to the corresponding author. The DEM co-registration algorithm from Nuth and Kääb (2011) is available on https://github.com/dshean/demcoreg (Shean and others, 2016). The NGI glacier outlines are available from https://www.glacairesargentinos.gob.ar; the RGI v6.0 outlines are available from the GLIMS initiative at https://www.glims.org/RGI/rgi60_dl.html. The elevation change grids from Dussaillant and others (2019) are distributed through PANGEA from https://doi.pangaea.de/10.1594/PANGAEA.903618. The Pléiades imagery used in this study was provided by the Pléiades Glacier Observatory (PGO) initiative of the French Space Agency (CNES). Pléiades data © CNES 2012 and 2020, Distribution Airbus D&S. The climate data were obtained for the Rio Neuquén catchment from the Autoridad Interjurisdiccional de las Cuencas de los ríos Limay, Neuquén y Negro (www.aic.gov.ar), and the streamflow data from Paso de los Indios gauge station were downloaded from the Subsecretaría de Recursos Hídricos de Argentina (https://www.argentina.gob.ar/obras-publicas/infraestructura-y-politica-hidrica).

**Acknowledgements.** We acknowledge funding by the Agencia Nacional de Promoción Científica y Técnica (Grant PICT 2007-0379 and PICT 2016-1282) and CONICET (Grant PIP 11220110100200). EB acknowledges support from the French Space Agency (CNES). We thank Sergio Cimbaro (IGN), Martín Matula (AEROMAPA), Mario Rosas, Paola Banegas, Ricardo Manzur (SEGEMAR Mendoza), Luis Lenzano, Pierre Pitte (IANIGLA) and David from Taller 4 (Mendoza City) for handling and scanning the photographic material. We also thank Tobias Bolch and Owen King (St. Andrews University) for their assistance on the mapping of debris-covered ice and language editing, respectively. The Check Points were surveyed and postprocessed using equipment and software lent by the Geomatics division of IANIGLA. We are grateful to Jan-Hendrik Malles and Ben Marzeion (University of Bremen, Germany) for kindly providing the ensemble-based model data. Special thanks go to the late Dr Eduardo Llambías (CIG-La Plata), Mariana Aubone, Héctor Valdéz, Daniel and Heraldo Castillo (Área Natural Protegida del Sistema Domuyo), for providing field support and photographs. Finally, we thank for the comments and suggestions provided by two anonymous reviewers and the Editors Evan Miles and Nicolas Cullen.

**References**

Abermann J, Steiner JF, Prinz R, Wecht M and Lisage P (2020) The Red Rock ice cliff revisited – six decades of frontal, mass and area changes in the Nunatarssuaq area, Northwest Greenland. *Journal of Glaciology* **66**(258), 567–576. doi: **10.1017/jog.2020.28**

AIC (2010) Estadísticas Climáticas e Hidrológicas. Periodo 2001-2010. Cuencas de los Ríos Limay, Neuquén y Negro. Secretaría Operativa y

Fiscalización. Cipolletti, Rio Negro, Argentina. Available at www.aic.gob.ar/sitio/archivos/201801/parte%201.pdf.

Anderson B and Mackintosh A (2012) Controls on mass balance sensitivity of maritime glaciers in the Southern Alps, New Zealand: the role of debris cover. *Journal of Geophysical Research* 117(F1), F01003. doi: 10.1029/2011JF002064

Andreassen LM, Elvehøy H, Kjøllmoen B and Belart JMC (2020) Glacier change in Norway since the 1960s – an overview of mass balance, area, length and surface elevation changes. *Journal of Glaciology* 66(256), 313–328. doi: 10.1017/jog.2020.10

Andreassen LM, Elvehøy H, Kjøllmoen B and Engeset RV (2016) Reanalysis of long-term series of glaciological and geodetic mass balance for 10 Norwegian glaciers. *The Cryosphere* 10(2), 535–552. doi: 10.5194/tc-10-535-2016

Andreassen LM, Elvehøy H, Kjøllmoen B, Engeset RV and Haakensen N (2005) Glacier mass-balance and length variation in Norway. *Annals of Glaciology* 42, 317–325. doi: 10.3189/172756405781812826

Astort A and 6 others (2019) Unrest at Domuyo Volcánó, Argentina, detected by geophysical and geodetic data and morphometric analysis. *Remote Sensing* 11(18), 2175. doi: 10.3390/rs11182175

Ayala Á and 5 others (2020) Glacier runoff variations since 1955 in the Maipo river basin, in the semiarid Andes of Central Chile. *The Cryosphere* 14(6), 2005–2027. doi: 10.5194/tc-14-2005-2020

Barcaza G and 7 others (2017) Glacier inventory and recent glacier variations in the Andes of Chile, South America. *Annals of Glaciology* 58, 166–180. doi: 10.1017/aog.2017.28

Belart MCJ and 5 others (2019) The geodetic mass balance of Eyjafjallajökull Ice Cap for 1945–2014: processing guidelines and relation to climate. *Journal of Glaciology* 65(251), 395–409. doi: 10.1017/jog.2019.16

Berthier E and 10 others (2014) Glacier topography and elevation changes derived from Pléiades sub-meter stereo images. *The Cryosphere* 8(6), 2275–2291. doi: 10.5194/tc-8-2275-2014

Berthier E, Cabot V, Vincent C and Six D (2016) Decadal region-wide and glacier-wide mass balances derived from multi-temporal ASTER satellite digital elevation models. Validation over the Mont-Blanc area. *Frontiers in Earth Science* 4, 63. doi: 10.3389/feart.2016.00063

Bessette-Kirton EK, Coe and JA and Zhou (2018) Using stereo satellite imagery to account for ablation, entrainment, and compaction in volume calculations for rock avalanches on glaciers: application to the 2016 Lamplugh rock avalanche in Glacier Bay National Park, Alaska. *Journal of Geophysical Research: Earth Surface* 123, 622–641. doi: 10.1002/2017JF004512

Błaszczyk M and 10 others (2019) Quality assessment and glaciological applications of digital elevation models derived from space-borne and aerial images over two tidewater glaciers of Southern Spitsbergen. *Remote Sensing* 11(9), 1121. doi: doi.org/10.3390/rs11091121

Braun MH and 8 others (2019) Constraining glacier elevation and mass changes in South America. *Nature Climate Change* 9(2), 130–136. doi: 10.1038/s41558-018-0375-7

Burger F and 7 others (2019) Interannual variability in glacier contribution to runoff from a high-elevation Andean catchment: understanding the role of debris cover in glacier hydrology. *Hydrological Processes* 33(2), 214–229. doi: 10.1002/hyp.13354

Cembrano J and Lara L (2009) The link between volcánism and tectonics in the southern volcánic zone of Chilean Andes: a review. *Tectonophysics* 471, 96–113. doi: 10.1016/j.tecto.2009.02.038

Chiodini G and 8 others (2014) The Domuyo Volcánic system: an enormous geothermal resource in Argentine Patagonia. *Journal of Volcanology and Geothermal Research* 274, 71–77. doi: 10.1016/j.jvolgeores.2014.02.006

Cogley GJ (2009) Geodetic and direct mass-balance measurements: comparison and joint analysis. *Annals of Glaciology* 50(50), 96–100. doi: 10.3189/17275640978769744

Condom T, Coudrain A, Sicart JE and Théry S (2007) Computation of the space and time evolution of equilibrium-line altitudes on Andean glaciers (10°N–55°S). *Global and Planetary Change* 59(1–4), 189–202. doi: 10.1016/j.gloplacha.2006.11.021

Denzinger F and 7 others (2021) Geodetic mass balance of Abramov Glacier from 1975 to 2015. *Journal of Glaciology* 67(262), 331–342. doi: 10.1017/jog.2020.108

Deshamps-Berger C and 7 others (2020) Snow depth mapping from stereo satellite imagery in mountainous terrain: evaluation using airborne lidar data. *The Cryosphere* 14, 2925–2940. doi: 10.5194/tc-14-2925-2020

Dussaillant I and 8 others (2019) Two decades of glacier mass loss along the Andes. *Nature Geoscience* 12(10), 802–808. doi: 10.1038/s41561-019-0432-5
Dussaillant I, Berthier E and Brun F (2018) Geodetic mass balance of the northern Patagonian icefield from 2000 to 2012 using two independent methods. *Frontiers in Earth Science* **6**, 8. doi: 10.3389/feart.2018.00008

Falaschi D and 6 others (2017) Mass changes of alpine glaciers at the eastern margin of the northern and southern Patagonian icefields between 2000 and 2012. *Journal of Glaciology* **63**(238), 258–272. doi: 10.1017/jog.2016.136

Falaschi D and 5 others (2019) Six decades (1958–2018) of geodetic glacier mass balance in Monte San Lorenzo, Patagonian Andes. *Frontiers in Earth Science* 7, 326. doi: 10.3389/feart.2019.00326

Falaschi D, Masiokas M, Tadono T and Couvreux F (2016) ALOS-derived glacier and rock glacier inventory of the Volcán Domuyo region (~36° S), southernmost Central Andes, Argentina. *Zeitschrift für Geomorphologie* **60**(3), 195–208. doi: 10.1127/zfg/2016/0319

Farías-Barahona D and 8 others (2019) Geodetic mass balances and area changes of Echaurren Norte Glacier (Central Andes, Chile) between 1955 and 2015. *Remote Sensing* **11**(3), 260. doi: 10.3390/rs11030260

Farías-Barahona D and 9 others (2020a) 60 Years of glacier elevation and mass changes in the Maipo river basin, Central Andes of Chile. *Remote Sensing* **12**(10), 1658. doi: 10.3390/rs12101658

Farías-Barahona D and 13 others (2020b) A near 90-year record of the evolution of El Morado Glacier and its proglacial lake, Central Chilean Andes. *Journal of Glaciology* **66**(259), 845–860. doi: 10.3390/rs11030260

Ferri L and 9 others (2020) Ice mass loss in the Central Andes of Argentina between 2000 and 2018 derived from a new glacier inventory and satellite stereo-imagery. *Frontiers in Earth Science* **8**, 530997. doi: 10.3389/feart.2020.530997

Gardelle J, Berthier E, Arnaud Y and Kääb A (2013) Region-wide glacier mass balances over the Pamir–Karakoram–Himalaya during 1999–2011. *The Cryosphere* **7**(4), 1263–1286. doi: 10.5194/tc-7-1263-2013

Gardner AS and 15 others (2013) A reconciled estimate of glacier contributions to sea level rise: 2003 to 2009. *Science* **340**(6134), 852–857. doi: 10.1126/science.1234532

Garread RD (2009) The Andes climate and weather. *Advances in Geosciences* **22**, 3–11. doi: 10.5194/adgeo-22-3-2009

Garread RD and 8 others (2017) The 2010–2015 megadrought in Central Chile: impacts on regional hydroclimate and vegetation. *Hydrology and Earth System Sciences* **21**, 6307–6327. doi: 10.5194/hess-21-6307-2017

Garread RD and 5 others (2019) The Central Chile mega drought (2010–2018): a climate dynamics perspective. *International Journal of Climatology* **40**(1), 421–439. doi: 10.1002/joc.6219

Girod L, Nuth C, Kääb A, McNabb R and Galland O (2017) MMASTER: improved ASTER DEMs for elevation change monitoring. *Remote Sensing* **9**(7), 704. doi: 10.3390/rs9070704

González-Reyes Á and 8 others (2017) Spatiotemporal variations in hydroclimatic across the Mediterranean Andes (30°–37°S) since the early twentieth century. *Journal of Hydrometeorology* **18**(7), 1929–1942. doi: 10.1175/JHM-D-16-0004.1

Hugonnet R and 10 others (2021) Accelerated global glacier mass loss in the early twenty-first century. *Nature* **592**(7856), 726–731. doi: 10.1038/s41586-021-03436-z

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

Kapitsa V and 5 others (2020) Assessment of changes in mass balance of the Tuyuksu group of glaciers, Northern Tien Shan, between 1958 and 2016 using ground-based observations and Pléiades satellite imagery. *Frontiers in Earth Science* **8**, 259. doi: 10.3389/feart.2020.00259

Kinnard C and 8 others (2020) Mass balance and climate history of Haig-altitude glacier, Desert Andes of Chile. *Frontiers in Earth Science* **8**, 40. doi: 10.3389/feart.2020.00040

Klug C and 8 others (2018) Geodetic reanalysis of annual glaciological mass balances (2001–2011) of Hintereisferner, Austria. *The Cryosphere* **12**(3), 833–849. doi: 10.5194/tc-12-833-2018

Kutuzov S, Lavrentiev I, Smirnov A, Nosenko G and Petrakov D (2019) Volume changes of Elbrus glaciers from 1997 to 2017. *Frontiers in Earth Science* **7**, 153. doi: 10.3389/feart.2019.00153

Lliboutry L (1998) Glaciers of Chile and Argentina. In Williams RS and Ferrigno JG (eds), *Satellite Image Atlas of Glaciers of the World. South America*, Vol. 1386–I.U.S. Denver, CO: Geological Survey, p. 1103.

Lundgren P and 8 others (2020) The dynamics of large silicic systems from satellite remote sensing observations: the intriguing case of Domuyo Volcáno, Argentina. *Scientific Reports* **10**(1), 11642. doi: 10.1038/s41598-020-67982-8

Malles J-H and Marzeion B (2021) Twentieth century global glacier mass change: an ensemble-based model reconstruction. *The Cryosphere* **15**(7), 3135–3157. doi: 10.5194/tc-15-3135-2021

Malmros JK, Mernild SH, Wilson R, Yde JC and Fensholt R (2016) Glacier area changes in the central Chilean and Argentinean Andes 1955–2013/14. *Journal of Glaciology* **62**(232), 391–401. doi: 10.1017/jog.2016.43

Malz P and 5 others (2018) Elevation and mass changes of the southern Patagonia icefield derived from TanDEM-X and SRTM data. *Remote Sensing* **10**(2), 188. doi: 10.3390/rs10020188

Marzeion B, Hofer M, Jarosch AH, Kaser G and Mölg T (2012) A minimal model for reconstructing interannual mass balance variability of glaciers in the European Alps. *The Cryosphere* **6**, 71–84. doi: 10.5194/tc-6-71-2012

Masiokas MH and 7 others (2012) Snowpack variations since AD 1150 in the Andes of Chile and Argentina (30°–37°S) inferred from rainfall, tree-ring and documentary records: Andean snowpack variations since AD 1150. *Journal of Geophysical Research: Atmospheres* **117**, D05112. doi: 10.1029/2011JD016748

Masiokas MH and 11 others (2016) Reconstructing the annual mass balance of the Echaurren Norte Glacier (Central Andes, 33.5°S) using local and regional hydroclimatic data. *The Cryosphere* **10**(2), 927–940. doi: 10.5194/tc-10-927-2016

Masiokas MH and 11 others (2020) A review of the current state and recent changes of the Andean cryosphere. *Frontiers in Earth Science* **8**, 99. doi: 10.3389/feart.2020.00099

McNabb R, Nuth C, Kääb A and Girod L (2019) Sensitivity of glacier volume change estimation to DEM void interpolation. *The Cryosphere* **13**(3), 895–910. doi: 10.5194/tc-13-895-2019

Melkonian AK and 5 others (2013) Satellite-derived volume loss rates and glacier speeds for the Cordillera Darwin icefield, Chile. *The Cryosphere* **7**(3), 823–839. doi: 10.5194/tc-7-823-2013

Miles ES and 8 others (2018) Glacial and geomorphic effects of a supraglacial lake drainage and outburst event, Everest region, Nepal Himalaya. *The Cryosphere* **12**(12), 3891–3905. doi: 10.5194/tc-12-3891-2018

Muñoz-Sabater J and 16 others (2021) ERA5-Land: a state-of-the-art global reanalysis dataset for land applications. *Earth Syst. Sci. Data* **13**(9), 4349–4383.

Nuth C and Kääb A (2011) Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change. *The Cryosphere* **5**(1), 271–290. doi: 10.5194/tc-5-271-2011

Papasodoro C, Berthier E, Royer A, Zdanowicz C and Langlois A (2015) Area, elevation and mass changes of the two southernmost ice caps of the Canadian Arctic Archipelago between 1952 and 2014. *The Cryosphere* **9**(4), 1535–1510. doi: 10.5194/tc-9-1535-2015

Pfeffer T and 19 others (2014) The Randolph glacier inventory: a globally complete inventory of glaciers. *Journal of Glaciology* **60**(221), 537–552. doi: 10.3189/2014JoG13J176

Pierrot Deseilligny M and Clery I (2011) Apero, an open source bundle adjustment software for automatic calibration and orientation of set of images. *The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences* **3816**, 269–276. doi: 10.5194/ispsarchives-XXXVIII-5-W16-269-2011

Reinthaler J, Paul F, Granados HD, Rivera A and Huggel C (2019) Area changes of glaciers on active Volcánoes in Latin America between 1986 and 2015 observed from multi-temporal satellite imagery. *Journal of Glaciology* **65**(252), 542–556. doi: doi.org/10.1017/jog.2019.30

Rignot E (2003) Contribution of the Patagonia icefields of South America to sea level rise. *Science* **302**(5644), 434–437. doi: 10.1126/science.1087393

Rivera A and 7 others (2015) Recent changes in total ice volume on Volcán Villarrica, Southern Chile. *Nat. Hazards* **75**(1), 33–55.

Rivera A and Bown F (2013) Recent glacier variations on active ice capped Volcánoes in the southern volcánic zone (37°–46°S), Chilean Andes. *Journal of South American Earth Sciences* **45**, 345–356. doi: 10.1016/j.jsames.2013.02.004

Rivera A, Bown F, Carrión D and Zenteno P (2012) Glacier responses to recent volcánic activity in southern Chile. *Environmental Research Letters* **7**(1), 014036. doi: 10.1088/1748-9326/7/1/014036

Rivera J, Penalba O, Villalba R and Araneo D (2017) Spatio-temporal patterns of the 2010–2015 extreme hydrological drought across the Central Andes, Argentina. *Water* **9**(9), 652. doi: doi.org/10.3390/w9090652

Robson BA and 5 others (2022) Glacier and rock glacier changes since the 1950s in the La Laguna catchment, Chile. *The Cryosphere* **16**(2), 647–665. doi: 10.5194/tc-16-647-2022
**Rupnik E, Daakir M and Pierrot Deseilligny M** (2017) MicMac – a free, open-source solution for photogrammetry. *Open Geospatial Data, Software and Standards* **2**, 14. doi: 10.1186/s40965-017-0027-2

**Sagredo EA and Lowell TV** (2012) Climatology of Andean glaciers: a framework to understand glacier response to climate change. *Global and Planetary Change* **86–87**, 101–109. doi: 10.1016/j.gloplacha.2012.02.010

**Schaefer M, Rodriguez JL, Scheiter M and Casassa G** (2017) Climate and surface mass balance of Mocho Glacier, Chilean Lake District, 40°S. *Journal of Glaciology* **63**(238), 218–228. doi: 10.1017/jog.2016.129

**Shean DE and 6 others** (2016) An automated, open-source pipeline for mass production of digital elevation models (DEMs) from very-high-resolution commercial stereo satellite imagery. *ISPRS Journal of Photogrammetry and Remote Sensing* **116**, 101–117. doi: 10.1016/j.isprsjprs.2016.03.012

**Shean DE and 5 others** (2020) A systematic, regional assessment of high mountain Asia glacier mass balance. *Frontiers in Earth Science* **7**, 363. doi: 10.3389/feart.2019.00363

**Wagnon P and 10 others** (2021) Reanalysing the 2007–19 glaciological mass-balance series of Mera Glacier, Nepal, Central Himalaya, using geodetic mass balance. *Journal of Glaciology* **67**(261), 117–125. doi: 10.1017/jog.2020.88

**WGMS** (2021) Global Glacier Change Bulletin No. 4 (2018–2019). In Zemp M, Nussbaumer SU, Gärtner-Roer I, Bannwart J, Paul F and Hoelzle M, eds. ISC(WDS)/IUGG(IACS)/UNEP/UNESCO/WMO, World Glacier Monitoring Service, Zurich, Switzerland. doi: 10.5904/wgms-fog-2021-05

**Wilby RL and 5 others** (2004) Guidelines for use of climate scenarios developed from statistical downscaling methods. *IPCC Task Group on Data and Scenario Support for Impact and Climate Analysis (TGICA)*, Geneva, Switzerland.

**Zalazar L and 9 others** (2020) Spatial distribution and characteristics of Andean ice masses in Argentina: results from the first national glacier inventory. *Journal of Glaciology* **66**(260), 938–949. doi: 10.1017/jog.2020.55

**Zekollari H, Huss M and Farinotti D** (2019) Modelling the future evolution of glaciers in the European Alps under the EURO-CORDEX RCM ensemble. *The Cryosphere* **13**(4), 1125–1146. doi: 10.5194/tc-13-1125-2019

**Zemp M and 16 others** (2013) Reanalysing glacier mass balance measurement series. The Cryosphere **7**(4), 1227–1245. doi: 10.5194/tc-7-1227-2013

**Zemp M and 14 others** (2019) Global glacier mass changes and their contributions to Sea-level rise from 1961 to 2016. *Nature* **568**(7752), 382–386. doi: 10.1038/s41586-019-1071-0