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

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^–1) 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^–1 between 1962 and 1984, to a strongly negative mass balance of −0.99 ± 0.19 m w.e. a^–1 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 Lara, 2009). It is therefore not surprising that glaciological studies on this portion of the Andes are sparse, focusing mainly on glacier 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.

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).

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 glacier–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.

Table 1. Summary of aerial and satellite data and characteristics of the derived DEMs/orthoimages

The footprint of the aerial (and satellite) imagery used in this study is displayed in Fig. S1.

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 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₂₀₂₀, 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 westeast 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₁₉₆₂ 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₁₉₆₂ 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₁₉₆₂ 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²):

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

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

(2)$$\Delta m = \displaystyle{{\Delta v \times \rho } \over {( {S_{\rm i} + S_{\rm f}} ) /2}} \times t^{{-}1}$$

ρ 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^–3 (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:

(3)$$\sigma \Delta v = \sqrt {{( {S_{\rm i} \times \sigma_{{\rm d}h}} ) }^2 + {( {{\rm d}h \times \sigma_{S_{\rm i}}} ) }^2} $$

where σ _dh is the uncertainty on the rate of elevation change and $\sigma _{S_{\rm 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 σƒ_Δv of ±60 kg m^–3 in the density assumption for volume to mass conversion (Huss, 2013):

(4)$$\displaystyle{{\sqrt {{( {\,f_{\Delta v} \times \sigma \Delta v} ) }^2 + {( {\Delta v \times \sigma {\rm f}_{\Delta v}} ) }^2\;} \;} \over {( {S_{\rm i} + S_{\rm f}} ) /2\;}} \times \;t^{{-}1}$$

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.

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

Validation of all DEMs, except for DEM₁₉₆₂, was carried out by comparing field-surveyed check points (CPs) off glacier against each DEM (Fig. S1). This was not possible for DEM₁₉₆₂ 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₂₀₂₀ relative to the CPs. This yielded a mean XY shift of 6.8 m (σ = 0.9 m) and a mean vertical (Z) shift of 1.6 m (σ = 1.4 m). Then, after co-registration of all derived orthoimages and DEMs with respect to the DEM₂₀₂₀, 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₂₀₂₀ and DEM₂₀₁₂, 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).

Table 2. Validation of the ASTER, Pléiades and photogrammetric (1984) DEM elevations compared to 14 field-surveyed GCPs

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.

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.

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 interannual 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 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:

(5)$${\rm MB} = \mathop \sum \limits_{i = 1}^{12} \;( {\alpha \;P_i-\;\mu ( {{\rm max}( {0, \;\;T_i-T_{{\rm melt}}} ) \;} ) } ) $$

where MB is the annual specific mass balance. P _i is the monthly precipitation and α is a scaling parameter for precipitation, T _i is the mean monthly air temperature at the front of the glacier, T _melt is the monthly mean temperature above which melt occurs and μ 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 _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 α = 3.9, μ = 90.1 mm °C⁻¹, and we prescribed T _melt = 0°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 ±1 value for the α parameter. Marzeion and others (2012) found a great range of values for the α parameter in the alps. Their mean value of 2 is lower than the α = 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 ± 1.6°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 km² (20.8%) at a rate of −0.08 km² 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 km²; ~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% a^–1.

Table 3. Glacier inventory and area changes

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

Based on the 1984 aerial photographs, we identified and mapped 46 glaciers larger than 0.01 km² that represented a total area of 29.0 ± 0.2 km². The largest inventoried glacier is Chadileo II (5.1 ± 0.02 km²). By 2020, five glacier units mapped in the 1984 orthophotos had disappeared, and the total glacier area had decreased to 20.5 ± 0.2 km². The total glacier area loss of 8.5 km² corresponds to 29.3% of the 1984 glacier area, and the area loss rate between 1984 and 2020 was ~0.25 km² a^–1 (0.8% 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).

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).

Table 4. Accumulated 1984–2020 and 1962–2020 volume change (Δv) and mass balance (Δm) for 11 specific glaciers and total glacier area

^a Refers to the listed glaciers in the table (minus Manchana Covunco) only.

Overall, the cumulative mass balance was −15.12 ± 3.60 m w.e. (–0.42 ± 0.10 w.e. a^–1) during 1984–2020 (46 glaciers), and the total volume loss was 12.2 ± 2.7 × 10⁶ m³ a^–1. During this period, Turbio I Glacier (−0.55 ± 0.11 m w.e. a^–1) had the most negative mass balance, whereas Covunco (−0.23 ± 0.09 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 m w.e. a^–1 for the total study area, while for individual glaciers residuals ranged between ±0.03 m w.e. a^–1.

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 (Figs 3, 4), mass-balance conditions varied significantly among the three sub-periods.

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).

From 1984 to 2003, the studied glaciers were close to an equilibrium state (−0.06 ± 0.09 m w.e. a^–1). Turbio II (−0.27 ± 0.20 m w.e. a^–1) and Turbio III (−0.23 ± 0.19 m w.e. a^–1) 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^–1), with the greatest mass loss measured for Domuyo Norte I (−0.67 ± 0.38 m w.e. a^–1) and Turbio I (−0.66 ± 0.31 m w.e. a^–1). 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^–1. During this period, Domuyo Norte II (−1.46 ± 0.33 m w.e. a^–1) 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^–1. Domuyo Norte I (+0.37 ± 0.20 m w.e. a^–1) 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^–1), and the mass balance among glaciers varied to a limited extent between 0 and −0.23 m w.e. a^–1 (Table 4).

Table 5. Time-interval volume change (Δv) and mass balance (Δm) for 11 specific glaciers and total glacier area

^a Refers to the listed glaciers in the table (minus Manchana Covunco) only.

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 reaches. The 1962–84 curve shows a strongly positive gradient between 2800 and 3000 m (shifting from −1.5 to 1.3 m a^–1), and remains under moderately positive dh/dt conditions (<0.5 m a^–1) 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.

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).

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^–1), 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^–1), 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^–1), and was concurrent with a mean summer temperature and winter precipitation anomalies of +0.23°C and −0.1, respectively.

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.

The mean streamflow of the Rio Neuquén at Paso de los Indios gauge station was 299.24 m³ s^–1 over the long-term 1903–2019 period. We found a slightly negative trend of −0.5 m³ s^–1 a^–1 over this interval (Fig. 8), representing an annual streamflow decrease of 58.0 m³ s^–1 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^–1, respectively). In addition, throughout these two periods, we observed streamflow increases of ~140.0 m³ s^–1 (i.e. 6.4 m³ s^–1 a^–1) and 47.0 m³ s^–1 (i.e. 2.5 m³ s^–1 a^–1), respectively. The 1984–2003 interval depicts a slight streamflow increase 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^–1), and strongly decreasing trends of −13.4 and −4.6 m³ s^–1 a^–1, 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).

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).

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^–1), 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.

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.

Table 6. Modelled glacier mass balance (1984–2020)

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

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^–1 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^–1) 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^–1 persisted for certain periods.

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.

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^–1 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^–1 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^–1 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^–1, 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^–1 recession rate during 1962–2020 is similar to the ~0.5% a^–1 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^–1. 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 within the Maipo catchment. In the case of Echaurren Norte, an average mass balance of −0.85 ± 0.08 m w.e. a^–1 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^–1 between 1955 and 1996, as revealed from photogrammetric surveys. No mass-balance data before 2000 for 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^–1 between 1975 and 1984, and −0.35 m w.e. a^–1 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^–1) and 1984–2003 (−0.40 m w.e. a^–1) 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 −0.36 ± 0.09 m w.e. a^–1 for the 2000–13 period and −1.20 ± 0.09 m w.e. a^–1 between 2009 and 2015. This time-variable pattern is confirmed at El Morado Glacier (~33.75° S) by a thinning rate of 1.36 ± 0.07 m a^–1 between 2000 and 2013, and subsequently by even higher rates (up to 4.1 m a^–1) between 2013 and 2015 (Farias-Barahona and others, 2020b). In general, glaciers in the Maipo catchment, which includes both the Echaurren Norte and El Morado glaciers, lost mass (from −0.15 ± 0.13 to −0.32 ± 0.13 m w.e. a^–1) at the sub-basin scale between 2000 and 2013 (Farias-Barahona and others, 2020a). In line with the above, Burger and others (2019) investigated three glaciers in the Río del Yeso catchment (~33.55° S), and found almost no elevation changes between 2000 and 2013, yet an increase in glacier thinning of <0.75 m a⁻¹ between 2013 and 2015.

Geodetic mass-balance records from 2000 to 2018 for all glaciers of the Randolph Glacier Inventory (RGI 6.0; Pfeffer and others, 2014) in the Northern Patagonian and Central Andes area are available from Dussaillant and others (2019). From their glacier-specific mass-balance database, we extracted an area-averaged mass balance of −0.04 ± 0.13 m w.e. a^–1 between 2000 and 2009 and −0.40 ± 0.07 m w.e. a^–1 between 2009 and 2018 for Volcán Domuyo (Figs 10a, b). We interpret that the differences between these values with our much more negative mass-balance signals (−0.54 ± 0.34 m w.e. a^–1 [2003–12] and −0.99 ± 19 m w.e. a^–1 [2012–20], respectively), likely lie in the usage of the RGI glacier outlines in Dussaillant and others (2019). For geodetic mass-balance determinations at a local scale, the RGI outlines are limited, as they miss great amounts of debris-covered ice, often have poorly defined drainage divides and overestimate glacier area by including seasonal snow and highly reflective, gypsum-rich rock outcrops.

Fig. 10. 2000–09 (a) and 2009–18 (b) elevation change maps derived from ASTER DEMs from Dussaillant and others (2019) and our own 2003 and 2012 glacier outlines for the Volcán Domuyo area, later used in Ferri and others (2020). (c) Elevation change as retrieved from Braun and others (2019) using the RGI (6.0) glacier outlines.

Improved glacier outlines from the National Glacier Inventory (NGI) of Argentina (Zalazar and other, 2020) (Fig. S4) were used by Ferri and others (2020) to derive glacier mass balance in the Central Andes, using the same ASTER DEM series of Dussaillant and others (2019). The mass balances for Volcán Domuyo were +0.11 ± 0.23 m w.e. a^–1 between 2000 and 2009, and −1.2 ± 0.23 m w.e. a^–1 between 2009 and 2018, which shows good agreement with our 2012–20 (−0.99 ± 19 m w.e. a^–1) estimate, but differs significantly from our 2003–12 results (−0.54 ± 0.34 m w.e. a^–1). If our 2003 and 2012 glacier outlines are used in combination with the ASTER-derived rate of elevation changes of Dussaillant and others (2019), the glacier-wide mass balance is −0.09 ± 0.23 m w.e. a^–1 between 2000 and 2009, and −1.11 ± 0.23 m w.e. a^–1 between 2009 and 2018.

The ‘ASTER monitoring of ice towards extinction “ASTERIX”’ dh/dt grids of Dussaillant and others (2019) and Ferri and others (2020) depict spatially dependent errors due to the ASTER satellite along-track oscillations described in Section 3.4, especially in the 2000–09 grids. Such errors can be observed as crests and valleys of predominantly negative or positive residuals over on- and off-glacier terrains (compare upper and lower part of Fig. 10a). We hence identify here a second source of mass-balance difference (the first one being the different glacier inventories) when comparing the ASTERIX-derived elevation change maps against our Domuyo estimates for a similar period. A third possible cause for the observed differences stem from the differences in the periods studied in Ferri and others (2020) and the present study.

Braun and others (2019) also produced an assessment of the geodetic mass balance of most glaciers in South America. From their available elevation change maps, we extracted an average elevation change of −0.04 ± 0.08 m a^–1 between 2000 and 2011/2015 for the Domuyo glaciers, which represents a much smaller mass loss compared to our results (Fig. 10c). Again, we attribute the difference to the usage of the RGI glacier outlines, and the possible penetration of the radar pulse in the snowpack, a common bias in the retrieved elevation of the synthetic aperture radar-derived DEMs (Dussaillant and others, 2018; Malz and others, 2018; Robson and others, 2022).

5.3 Climate trends and the geodetic record in the Volcán Domuyo area

While the temperature and precipitation evolution in the Domuyo region prior to 1981 cannot be tracked from the Pampa de Chacaico records, Falaschi and others (2016) showed that the precipitation at the Pampa de Chacaico site is highly correlated with the regional snow accumulation series based on eight snowpack 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 Rio 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 uppermost 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^–1) is consistent with a marked streamflow increase of ~140.0 m³ s^–1 over this interval, and a mean streamflow value of 302.4 m³ s^–1, 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^–1, and pronounced runoff reductions of ~140.0 and ~41.0 m³ s^–1 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 Rio 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^–1) 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^–1) between 1962 and 1984, to a strongly negative mass balance of −0.99 ± 0.19 m w.e. a^–1 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 driver 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 absolute 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.