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, Reference Sagredo and Lowell2012). 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, Reference Masiokas2020 and references therein). Andean glaciers have also been acknowledged as one of the major contributors to sea level rise (Gardner and others, Reference Gardner2013; Zemp and others, Reference Zemp2019; Hugonnet and others, Reference Hugonnet2021). 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, Reference Cogley2009). 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, Reference Melkonian2013; Falaschi and others, Reference Falaschi2017, Reference Falaschi2019; Dussaillant and others, Reference Dussaillant, Berthier and Brun2018, Reference Dussaillant2019; Braun and others, Reference Braun2019; Kinnard and others, Reference Kinnard2020). 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, Reference Zemp2013; Andreassen and others, Reference Andreassen, Elvehøy, Kjøllmoen and Engeset2016; Klug and others, Reference Klug2018). Unfortunately, despite a few exceptions (e.g. Rignot, Reference Rignot2003; Falaschi and others, Reference Falaschi2019; Farías Barahona and others, Reference Farías-Barahona2019, Reference Farías-Barahona2020a, Reference Farías-Barahona2020b), 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, Reference Lliboutry, Williams and Ferrigno1998). 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, Reference Condom, Coudrain, Sicart and Théry2007) and consequently, the scarce glacier ice extent is restricted to the highest volcanic edifices (Reinthaler and others, Reference Reinthaler, Paul, Granados, Rivera and Huggel2019) of the Southern Volcanic Zone (33–46° S; Cembrano and Lara, Reference Cembrano and Lara2009). 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, Reference Rivera, Bown, Carrión and Zenteno2012; Rivera and Bown, Reference Rivera and Bown2013; Reinthaler and others, Reference Reinthaler, Paul, Granados, Rivera and Huggel2019). Only recent large-scale remote-sensing studies have looked at glacier mass changes in this region (Braun and others, Reference Braun2019; Dussaillant and others, Reference Dussaillant2019; Ferri and others, Reference Ferri2020) 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, Reference Falaschi, Masiokas, Tadono and Couvreux2016) are more typical of the Central Andes of Argentina and Chile (Dussaillant and others, Reference Dussaillant2019). 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, Reference Chiodini2014). 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.
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 (Reference Braun2019) and Dussaillant and others (Reference Dussaillant2019), 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 ( Reference Rivera2015) estimated the overall thinning since 1961 in Volcán Villarrica (39° S), whereas Schaefer at al. (Reference Schaefer, Rodriguez, Scheiter and Casassa2017) 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, Reference Reinthaler, Paul, Granados, Rivera and Huggel2019), rather it hosts mountain glaciers and ice aprons. Falaschi and others (Reference Falaschi, Masiokas, Tadono and Couvreux2016) produced an inventory of Domuyo glaciers accounting for ~25.4 km2 of ice for 2009 and an estimated 26% (~1.4% a−1) 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 (Reference Garreaud2009), 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, Reference Sagredo and Lowell2012). 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 (Reference Condom, Coudrain, Sicart and Théry2007).
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, Reference Shean2016) using the Semi-Global-Matching algorithm following Shean and others (Reference Shean2020). Altimetric precision of individual ASTER DEMs ranges from 5 to 15 m (Girod and others, Reference Girod, Nuth, Kääb, McNabb and Galland2017; Shean and others, Reference Shean2020). 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 (DEM2012 and DEM2020 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 (Reference Deshamps-Berger2020) 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, Reference Berthier2014). They have been already used to derive glacier elevation changes in diverse-glacierized regions worldwide (Kutuzov and others, Reference Kutuzov, Lavrentiev, Smirnov, Nosenko and Petrakov2019; Abermann and others, Reference Abermann, Steiner, Prinz, Wecht and Lisage2020; Andreassen and others, Reference Andreassen, Elvehøy, Kjøllmoen and Belart2020; Kapitsa and others, Reference Kapitsa2020; Denzinger and others, Reference Denzinger2021).
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.
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, Reference Pierrot Deseilligny and Clery2011; Rupnik and others, Reference Rupnik, Daakir and Pierrot Deseilligny2017), following Belart and others (Reference Belart2019). 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 (Reference Papasodoro, Berthier, Royer, Zdanowicz and Langlois2015), 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 (Reference Falaschi, Masiokas, Tadono and Couvreux2016), 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 km2). Still, from 1962 to 2020, we report area and mass changes for the ten largest glaciers (>0.5 km2) which total 22.6 km2 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 (Reference Wagnon2021) 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 DEM2020, any remaining shifts (which we term triangulation following Nuth and Kaab, Reference Nuth and Kääb2011) 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, Reference Belart2019). 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, Reference Girod, Nuth, Kääb, McNabb and Galland2017; Dussaillant and others, Reference Dussaillant2019), and can be attributed to satellite attitude oscillations along-track (jitter) (Deschamps-Berger and others, Reference Deshamps-Berger2020). Similar undulations have been previously found in other satellite products as well, including WorldView DEMs (Shean and others, Reference Shean2016; Bessette-Kirton and others, Reference Bessette-Kirton, Coe and Zhou2018). 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, Reference Gardelle, Berthier, Arnaud and Kääb2013).
Examination of dh/dt grids including the DEM1962 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 DEM1962 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 (Reference Nuth and Kääb2011) co-registration algorithm implemented by Shean and others (Reference Shean2016). 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 DEM1962 missed 14% of the uppermost part of Turbio I Glacier. For void filling, we followed the local hypsometric approach of McNabb and others (Reference McNabb, Nuth, Kääb and Girod2019), 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 m3) is calculated as the product of the average elevation change Δh (m) and the initial glacier area S i (m2):
Then, the glacier specific geodetic mass balance Δm (m w.e. a–1) can be calculated as:
ρ 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, Reference Huss2013).
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:
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, Reference Berthier, Cabot, Vincent and Six2016). 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, Reference Miles2018; Wagnon and others, Reference Wagnon2021 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, Reference Huss2013):
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, Reference Andreassen, Elvehøy, Kjøllmoen, Engeset and Haakensen2005; Belart and others, Reference Belart2019). 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 DEM1962, was carried out by comparing field-surveyed check points (CPs) off glacier against each DEM (Fig. S1). This was not possible for DEM1962 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 DEM2020 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 DEM2020, we applied the systematic shift above and compared the Z coordinate of each DEM against the CPs (Table 2). Particularly, in the cases of DEM2020 and DEM2012, 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, Reference Berthier2014; Błaszczyk and others, Reference Błaszczyk2019).
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, Reference Muñoz-Sabater2021). 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, Reference Wilby2004). An RMSE of 0.9°C was obtained after the correction.
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 (Reference Marzeion, Hofer, Jarosch, Kaser and Mölg2012) for glaciers in the alps and has been subsequently applied to reconstruct the annual SMB of Echaurren Norte Glacier (Masiokas and others, Reference Masiokas2016) and El Morado Glacier in central Chile (Farias-Barahona and others, Reference Farías-Barahona2020b). The model is defined as:
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 (Reference Masiokas2016) 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−1, 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 (Reference Marzeion, Hofer, Jarosch, Kaser and Mölg2012) 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 (Reference Masiokas2016), 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, Reference Anderson and Mackintosh2012), 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 km2 (20.8%) at a rate of −0.08 km2 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 km2; ~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.
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 km2 that represented a total area of 29.0 ± 0.2 km2. The largest inventoried glacier is Chadileo II (5.1 ± 0.02 km2). By 2020, five glacier units mapped in the 1984 orthophotos had disappeared, and the total glacier area had decreased to 20.5 ± 0.2 km2. The total glacier area loss of 8.5 km2 corresponds to 29.3% of the 1984 glacier area, and the area loss rate between 1984 and 2020 was ~0.25 km2 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).
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).
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 × 106 m3 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.
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 (Reference Farías-Barahona2019, Reference Farías-Barahona2020b).
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).
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.
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−1) and a decrease in winter precipitation (−50 mm decade−1). 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.
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 m3 s–1 over the long-term 1903–2019 period. We found a slightly negative trend of −0.5 m3 s–1 a–1 over this interval (Fig. 8), representing an annual streamflow decrease of 58.0 m3 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 m3 s–1, respectively). In addition, throughout these two periods, we observed streamflow increases of ~140.0 m3 s–1 (i.e. 6.4 m3 s–1 a–1) and 47.0 m3 s–1 (i.e. 2.5 m3 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 m3 s–1), and strongly decreasing trends of −13.4 and −4.6 m3 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).
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.
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.
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.
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 (Reference Falaschi, Masiokas, Tadono and Couvreux2016), 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−1 on average), and a later period (1984–2020) when glacier area loss rates increased to ~0.6% a−1. 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 (Reference Reinthaler, Paul, Granados, Rivera and Huggel2019) 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 (Reference Reinthaler, Paul, Granados, Rivera and Huggel2019) 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 (Reference Malmros, Mernild, Wilson, Yde and Fensholt2016) 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 (Reference Barcaza2017) 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 (Reference Farías-Barahona2020a) using topographic maps and the SRTM DEM over the 1955–2000 period. The geodetic mass balance for the entire catchment (361.4 km2 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, Reference González-Reyes2017; Ayala and others, Reference Ayala2020).
Farias-Barahona and others (Reference Farías-Barahona2019, Reference Farías-Barahona2020b) also assessed individual geodetic mass balances of the World Glacier Monitoring Service (WGMS) reference Echaurren Norte (33°34′ S–70°7′ W; 0.17 km2) and El Morado (33°44′ S–70°3′ W; 1.05 km2) 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, Reference Masiokas2012) 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, (Reference Farías-Barahona2019) 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, Reference Farías-Barahona2020b). 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, Reference Farías-Barahona2020a). In line with the above, Burger and others (Reference Burger2019) 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−1 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, Reference Pfeffer2014) in the Northern Patagonian and Central Andes area are available from Dussaillant and others (Reference Dussaillant2019). 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 (Reference Dussaillant2019). 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.
2000–09 (a) and 2009–18 (b) elevation change maps derived from ASTER DEMs from Dussaillant and others (Reference Dussaillant2019) and our own 2003 and 2012 glacier outlines for the Volcán Domuyo area, later used in Ferri and others (Reference Ferri2020). (c) Elevation change as retrieved from Braun and others (Reference Braun2019) 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 (Reference Ferri2020) to derive glacier mass balance in the Central Andes, using the same ASTER DEM series of Dussaillant and others (Reference Dussaillant2019). 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 (Reference Dussaillant2019), 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 (Reference Dussaillant2019) and Ferri and others (Reference Ferri2020) 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 (Reference Ferri2020) and the present study.
Braun and others (Reference Braun2019) 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, Reference Dussaillant, Berthier and Brun2018; Malz and others, Reference Malz2018; Robson and others, Reference Robson2022).
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 (Reference Falaschi, Masiokas, Tadono and Couvreux2016) 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, Reference Masiokas2012). 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−1), 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, Reference Garreaud2017, Reference Garreaud2019; Rivera and others, Reference Rivera, Penalba, Villalba and Araneo2017; Dussaillant and others, Reference Dussaillant2019).
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 m3 s–1 over this interval, and a mean streamflow value of 302.4 m3 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 m3 s–1, and pronounced runoff reductions of ~140.0 and ~41.0 m3 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, Reference Zalazar2020) 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 (Reference Marzeion, Hofer, Jarosch, Kaser and Mölg2012) 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, Reference Zekollari, Huss and Farinotti2019). 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, Reference Astort2019; Lundgren and others, Reference Lundgren2020) 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 (Reference Malles and Marzeion2021) 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−1) and a decrease in winter precipitation (−50 mm decade−1) 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.
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 (Reference Nuth and Kääb2011) is available on https://github.com/dshean/demcoreg (Shean and others, Reference Shean2016). The NGI glacier outlines are available from https://www.glaciaresargentinos.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 (Reference Dussaillant2019) 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 of 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.
Open access
