Satellite image archives are increasingly vast, as well as open-access, and computational advances have led to methodological frameworks that allow for robust analyses of glacier change that we employ here in the Cordilleras Vilcanota, Vilcabamba and Urubamba. The Landsat archive provides a record of ice area change almost continuously from 1972 to the present day, interrupted only by periods of heavy cloud and snow and short discontinuities in satellite operation. The small size of some glaciers in this region presents a challenge to delineating ice areas, but for decadal rates of change the measured rates of recession frequently exceed any uncertainties in the data.

Deriving volumetric changes from these stereoscopic sensors relies on there being sufficient contrast over snow and ice-covered surfaces for feature matching to be successful. Reconnaissance images (e.g. Corona KH-4; Bolch et al., 2008) acquired over the high-elevation accumulation areas of the Cordilleras Vilcanota, Vilcabamba, and Urubamba are largely featureless however, meaning surface elevations cannot be extracted. For this reason, we focussed our mass balance estimations on the ASTER archive, which provides a continuous record from 2000 to the present day of stereoscopic scenes across the world. ASTER acquires images every 16 days through a 14-band nadir sensor and a single band backward looking sensor to acquire stereo images. These data are collected at 15 and 30 m resolution to produce 30 m-resolution digital elevation models (DEMs). Cloud and snow also interfere with the production of digital elevation models, but the huge data redundancy afforded by the archive means this issue is relatively easy to overcome.

Orthorectified Landsat-2 scenes (Level L1TP, 60 m spatial resolution; Supplementary Table S1) were used to manually delineate glacier outlines for the 1970s, conducted in ArcMap. All Landsat-2 scenes were acquired between the 25th June—31st July 1975 to minimise the influence of snow cover. All other glacier outlines were produced from 30 m imagery from the full-archive Landsat-5 Tier 1 Top of Atmosphere (1984–2012) and Landsat-8 Tier 1 Top of Atmosphere (2013–2020) datasets using a Normalised Difference Snow/Water Index (NDSI, NDWI) in Google Earth Engine. A cloud-free composite was produced for each year by extracting the median pixel value from the driest months (July 1st to September 30th) to minimise snow and cloud cover. Due to heavy cloud in these regions, there are some years within the time series with no usable imagery. The resulting images were split into their respective boundaries from the Randolph Glacier Inventory 6.0 (RGI) dataset (RGI Consortium, 2017) and manually extended outwards using the ASTER GDEM. On inspection, a small number of RGI outlines in these regions seemed highly improbable, and so we manually delineated flowlines where it was clear the RGI database was merging multiple glaciers together. All shapefiles were clipped to remove polygons smaller than 0.01 km² in area and manually inspected to ensure no ice-free area was being incorrectly classified.

When determining glacier area, we calculated a conservative estimate of uncertainty for Landsat-2 scenes by assuming all glacier periphery pixels were misclassified, and calculated the total area of these periphery pixels, representing an error of 11%. Studies which have calculated uncertainty in a similar way produce comparable estimates of uncertainty (e.g. 11% for Kochtitzky et al., 2018). Rocky outcrops between glaciers and nunataks were particularly problematic to determine whether they were glacierised or snow-covered. Glacier boundaries from 1984 onwards were calculated via an automated method. Where this method is applied over predominantly clean-ice, associated error is 2–5% (Paul et al., 2013), and so we use a conservative estimate of 5% as the uncertainty value for glacier area delineated in this way, as the Southern Peruvian Andes are mostly clean-ice. We acknowledge that this method excludes any small debris cover patches that exist, but this is consistent with other studies who also exclude these areas due to their scarcity in this region (Salzmann et al., 2013).

To calculate glacier volume change, we extracted all L1A ASTER tiles collected over these regions in the 20-year timespan from 1st January 2000 to 31st December 2020. We derived 631 DEMs using the AMES Stereo Pipeline with a 7 pixel correlation kernel, using a 12 m TanDEM-X tile seed. To ensure accurate co-registration, they were corrected for planimetric and altimetric shifts against the TanDEM-X tile with a third-order polynomial transformation in stable, off-ice areas using the method of Nuth and Kääb (2011). Using a higher resolution DEM (12 m) compared to Hugonnet et al. (2021) (90 m) for co-registration should result in more accurate elevation models. DEMs that passed this correction phase were uploaded to Google Earth Engine using geeup (Roy, 2021), and clipped to elevations between 2000 and 6500 m a.s.l. (the approximate range of elevation within these Cordilleras) to exclude grossly erroneous pixels. A single value per pixel per year was extracted from all DEM tiles by collecting the median value (excluding values above and below the median ±100 m). Using these single values per year, a linear fit through the surface elevation was then calculated over three time periods: 2000–2010; 2010–2020; and 2000–2020 to represent elevation change per year during these intervals. Hugonnet et al. (2021) produce elevation change at 5-year intervals, but the associated error is often far larger than in the decadal intervals, and so we consider only the 10-year products to be usable. More detailed methodology including the parameters used at each stage is given in Supplementary Datasheet S2.

For the conversion of elevation change to mass change, there are three key sources of uncertainty:

1) Error in glacier area ( σ A )

We followed the method presented by Brun et al. (2017) to calculate overall uncertainty for mass balance calculations. We calculate uncertainty for the region-wide mass loss and use this fraction to attribute uncertainty in mass balance. Uncertainty in volume change is calculated using the standard deviation of elevation change on stable ground ( σ Δ h ) . We cast a 2 km buffer outwards from the glacier margin into stable ground and used the standard deviation of all pixels in this buffer area as σ Δ h . We follow Brun et al. (2017) in using a decorrelation length of 500 m, but use a glacier area uncertainty ( σ A ) of 5% as we have greater confidence that this region is predominantly clean-ice. Finally, we use an ice density uncertainty ( σ f Δ V ) of 60 kg m⁻³ following Huss (2013).

Geodetic mass balance was calculated by taking into account the elevation change per pixel, pixel size (30 m), and ice density of 850 ± 60 kg m³ (Huss, 2013). Calculation of surface mass balance, which would include estimates of glacier flow from modelled ice thickness and velocity (Miles et al., 2021), was not possible in this region owing to lack of available flow data. This precluded a robust estimate of the glacier ELA using our surface elevation data alone. Therefore, in the absence of any in-situ data, we investigated previously published methods to quantify ELA change through time. These included quantifying changes to the freezing-line altitude (e.g. Schauwecker et al., 2017; Yarleque et al., 2018), the end of summer snowline altitude (Condom et al., 2007; Rabatel et al., 2013b), and the use of the median glacier elevation as an ELA proxy (e.g. Braithwaite and Raper, 2009; King et al., 2017). We adopted this final approach in the absence of higher spatial resolution climate data or reliable cloud-free satellite data to quantify snowlines, though we acknowledge this is likely to underestimate actual ELAs as the glaciers are out of balance with the climate. We calculated median glacier elevation for every glacier, and herein summarise an area-weighted average region-wide median glacier elevation for the Cordilleras Vilcabamba and Urubamba, and separate the Cordillera Vilcanota into five sub-regions reflecting the differing altitudes and topographies of each (Supplementary Figure S3).

To identify the most significant controls on ice loss in each region, we correlated the variables of change (volume and area loss) against topographic and climate data. Topographic data (aspect, elevation, slope) were derived from the ASTER GDEM and averaged over the shapefiles of each glacier. Climate data were created using the Weather Research and Forecasting (WRF) model and bias-corrected against daily precipitation and maximum and minimum daily temperature data in the region, from 1980 to 2018. Frost days were calculated as the number of days in the year where the minimum daily temperature fell below 0°C. For full details of the WRF model setup and bias-correction, see (Fyffe et al., 2021).

Climate data were correlated to glacier area change from 1984 (the first year of overlap between the Landsat record and climate dataset) to 2018 (the end of the climate dataset) for every glacier individually. These coefficients were then transformed to z-values using Fisher’s Z to ensure a normal distribution of data. A mean was taken of these Fisher’s Z values, which was subsequently back-transformed to achieve an overall region-wide correlation coefficient between ice area and climatic variables. Significance for these values was calculated using the Harmonic mean p-value, weighted to the size of the glacier, as a robust measure of significance in large datasets (Wilson, 2019). These steps ensure each glacier is statistically analysed against only the climate data that are of most relevance to that particular glacier.

Overall, glacierised area of the Southern Peruvian Andes shrunk from 833.6 ± 91.7 km² in 1975 to 312.6 ± 15.6 km² in 2020 (Figure 2). This loss has predominantly occurred at low elevations (Figure 3). The Cordillera Vilcanota shrunk from a regional total of 540.6 ± 59.4 km² in 1975 to 246.8 ± 12.3 km² in 2020. This represents a loss of 54.3% of the glacier ice that was present in 1975 over a 45-year time period. 82 of the 257 glaciers present in 1975 had been lost by 2020 (or shrunk smaller than 0.01 km² in size). There has been consistent and steady recession of ice across the region in this time. The Quelccaya ice cap, which is included within the area for the Cordillera Vilcanota, shrunk from 57.4 ± 5.4 km² in 1975 to 41.6 ± 2.1 km² in 2020, representing a loss (27.5%) that is below the regional average (Figure 2). There has, in contrast to the region-wide Cordillera Vilcanota, been an observable slowdown of ice area loss in recent years.

The Cordillera Vilcabamba shrunk from 226.1 ± 24.9 km² in 1975 to 98.9 ± 4.9 km² in 2020, a loss of 56% of glacierised area. 25 of the 169 glaciers present in 1975 had been lost by 2020. The Cordillera Urubamba shrunk from a glacierised area of 66.9 ± 7.4 km² in 1975 to 24.2 ± 1.2 km² in 2020; a loss of 64%. This represents the largest relative decline in glacierised area of the regions as this is the lowest-lying of the three Cordilleras. Of the 39 glaciers in the Cordillera Urubamba in 1975, seven had been lost by 2020.

These area changes, which predominantly affect lower-lying zones, have driven the glaciers to higher elevations over the 45-year study period (Figure 3). It is notable that in all three Cordilleras, glacier ice has all but disappeared from elevations below 4,700 m a.s.l. but is largely preserved above elevations of 5,500 m a.s.l. In each case, the hypsometry is skewed towards lower elevations. In contrast, the hypsometry of Quelccaya is skewed towards higher elevations, and consequently the bins with greatest area have been largely preserved. For Quelccaya, there is no observed change in glacier area above 5,550 m a.s.l. from 1975 to 2020.

Across the Cordillera Vilcanota, there is a moderate correlation between mean annual air temperature (which has risen by ∼0.2°C per decade from 1980 to 2020) and rate of ice area change (R = -0.514, p < 0.01), and maximum annual temperature and ice area loss (R = -0.576, p < 0.01) (Figure 3; Table 1). There is no statistically significant link (at p < 0.01) between ice area change and precipitation, but there is a corresponding correlation between the number of frost days in the year and ice area loss. The key topographic control on rate of ice area change is maximum elevation (R = −0.394, p < 0.01), with a second-order control of slope (R = −0.269, p < 0.01). There is no significant (p < 0.01) co-linearity between maximum elevation and slope. The prevalence of cloud or snow-covered scenes in our record limits our ability to analyse the role that ENSO plays in affecting glacier area change (shown in Figure 2).

Over the Cordilleras Vilcabamba and Urubamba, there is also a moderate correlation between maximum annual temperature and rate of ice area change (R = −0.338, p < 0.01 and R = −0.576, p < 0.01 respectively). There is no statistically significant link between ice area change and precipitation. Similar to the Cordillera Vilcanota, the primary topographic control on ice area change in the Cordilleras Vilcabamba and Urubamba is maximum elevation (R = −0.241, p < 0.01; R = −0.677, p < 0.01 respectively) though this driver is much weaker in the Cordillera Vilcabamba. Slope is not a statistically significant control in the Cordillera Vilcabamba, but is in the Cordillera Urubamba (R = −0.301, p < 0.01), though there is co-linearity between elevation and slope in the Urubamba (R = 0.515, p < 0.01).

Between 2000 and 2020, the Cordillera Vilcanota shrunk at a rate of −0.159 ± 0.022 Gt yr⁻¹ (−0.48 ± 0.07 m w.e. yr⁻¹). There was no significant change in the mass balance between the two sub-decadal periods; from −0.49 ± 0.17 m w.e. yr⁻¹ between 2000–2010 to −0.44 ± 0.16 m w.e. yr⁻¹ between 2010–2020 (Figure 4). Owing to the cloudiness of ASTER data, smaller glaciers, and steeper slopes over the Cordilleras Vilcabamba and Urubamba, the uncertainty in our mass balance calculation was much larger than it was for the Vilcanota. Between 2000 and 2020, the Cordillera Vilcabamba shrunk at a rate of −0.405 ± 0.451 m w.e. yr⁻¹, while the Cordillera Urubamba shrunk at a rate of −0.559 ± 1.575 m w.e. yr⁻¹. We present these figures here for completeness, but do not analyse them further given the extent to which the uncertainty exceeds any signal in the data. Thinning of debris-covered ice, which is not included in the glacier area polygons but only makes up a minimal percentage of ice in the region, can be observed in the surface elevation change (e.g. Figures 4B,C).

There is a moderate correlation between volume change over the full 20-year period and elevation (R = 0.38, p < 0.001), suggesting that low-lying areas are thinning, as well as receding, the most. There is a slight, but significant, correlation between volume change and slope (R = 0.16, p < 0.001) and aspect (R = −0.04, p < 0.001). This is particularly noticeable over the Nevado Ausangate cluster of glaciers, where South-facing glaciers are largely shadowed from solar radiation, while North-facing glaciers are experiencing large thinning (Figure 4B).

The median elevation (which we use as a proxy for ELA) of the Cordilleras Vilcabamba and Urubamba rose at the fastest rates across the full 1975–2020 time series (3.28 m yr⁻¹ and 3.27 m yr⁻¹ respectively), from 4,966 m a.s.l. to 5,110 m a.s.l. for the Cordillera Vilcabamba and 4,930 m a.s.l. to 5,077 m a.s.l. for the Cordillera Urubamba (Figure 5). However, the rate of change between the regions through time varies. The median elevation of the Cordillera Urubamba is accelerating upwards (from 2.88 m yr⁻¹ for 2000–2010 to 4.59 m yr⁻¹ for 2010–2020), while the rate of change in the Cordillera Vilcabamba is decelerating (from 3.69 m yr⁻¹ for 2000–2010 to 3.06 m yr⁻¹ for 2010–2020) (Supplementary Figure S5).

The rate of median elevation change over the Cordillera Vilcanota is lower in all sub-regions than the Cordilleras Vilcabamba and Urubamba. The northernmost region of the Cordillera Vilcanota (Nevado Qolquepunco) has the lowest median elevation of the region (at 5,163 m a.s.l. in 2020; compared to 5,451 m a.s.l. at Quelccaya), but the rate of change has remained steady in recent decades. Quelccaya exhibits the lowest rate of median elevation change of all the regions at 1.59 m yr⁻¹ from 1975 to 2020, likely due to its high elevation and flat topography. This rate of change has decelerated to 1.06 m yr⁻¹ between 2010 and 2020. The key driver for the rate of median elevation change across the Cordillera Vilcanota is temperature; there is a significant correlation between maximum annual temperature and median elevation (R = 0.603, p < 0.001). The role of precipitation and frost in driving median elevation was not significant at p < 0.01.

Published estimates of glacier area in these regions are scarce, and many studies limit their efforts to studying specific watersheds within the wider region (e.g. Veettil and de Souza, 2017; Drenkhan et al., 2018). Furthermore, definitions of the boundaries of the Cordilleras are not fixed, with differing assessments of how many smaller periphery glaciers are included in calculations of the glacierised area of the region (Hanshaw and Bookhagen, 2014). Nevertheless, being cognisant of the slight variability in study areas, our estimations of ice recession fit well with the trend presented by other publications (Figure 6). To our knowledge, only the National Inventory of Glaciers presents information on glacier area change in the Cordilleras Vilcabamba and Urubamba. Our results for the area of the Cordillera Vilcabamba in 2017 (100.07 ± 5.00 km²) are comparable to the INAIGEM data (95.54 km²), but are slightly higher for the Cordillera Urubamba (25.27 ± 1.26 km² to the INAIGEM’s 24.92 km²). Note there is a discrepancy between the INAIGEM published estimate for the Cordillera Urubamba of 17.89 km² and their dataset of 24.92 km², likely due to the way boundaries are drawn between Cordilleras (INAIGEM, 2017; ANA, 2021).

There have been no comparable estimates of mass balance over the Cordillera Vilcanota specifically; however, the global ASTER dataset of glacier change (Hugonnet et al., 2021) provides reasonable context for our findings. The primary difference between the two datasets is their spatial resolution: the current study reports data at 30 m spatial resolution, while Hugonnet et al. (2021) is downsampled to 100 m resolution. Our mass balance estimate over the 20-year time period (−0.477 ± 0.067 m w.e. yr⁻¹) is nearly identical to using their dataset (−0.479 ± 0.048 m w.e. yr⁻¹). We attribute our slightly higher level of uncertainty to our finer spatial resolution which introduces some small artefacts into our dataset. However, the trade-off is that we are able to observe much more closely the changes of individual glaciers (an example comparison is shown over Chumpe glacier, in the Cordillera Vilcanota in Figure 7). The problem of very large uncertainty of mass balance change over the Cordilleras Vilcabamba and Urubamba was not unique to our dataset as Hugonnet et al. (2021) had large data gaps over this region, rendering their data unusable. Our high errors were due to cloud rather than steep terrain. Over the Cordilleras Vilcabamba and Urubamba, we produced far fewer DEMs per pixel (∼20 DEMs over the 20-year time period) in comparison to the Cordillera Vilcanota (∼75 DEMs per pixel).

Our 2020 median glacier elevation estimate of 5,357 m a.s.l. for the central sub-region of the Cordillera Vilcanota is comparable to the 5,399 m a.s.l. (2016 figure) calculated by Drenkhan et al. (2018), with the difference attributed to different spatial areas covered (their study includes only glaciers in the Vilcanota-Urubamba-Vilcabamba basin). Our 1975–2020 rate of median elevation increase over the Quelccaya ice cap (15.9 m per decade) is comparable to the 1980–2017 estimate of 16.3 m per decade over the Quelccaya ice cap estimated by Yarleque et al. (2018). Our 1975–2020 rate of median elevation increase over the Central region of the Cordillera Vilcanota (24.8 m per decade) is slightly higher than the 22.8 m per decade between 1988 and 2010 estimated by Hanshaw and Bookhagen (2014). Our changes in the rate of median glacier elevation change show high variability across regions and decadal intervals (Supplementary Figure S6). This is not inconsistent with the interannual variability of measured ELA change over glaciers in the Cordillera Blanca (Kaser and Georges, 1997; Condom et al., 2007) and Cordillera Ampato (Veettil et al., 2016) which have been attributed to fluctuations in air humidity and the role of ENSO respectively. The interannual variation we exhibit in the Cordillera Vilcanota correlated to maximum air temperature, though this warrants further investigation given the high variability between sub-regions.

Ice area retreat is prolific across all areas of the Southern Peruvian Cordilleras. In the Cordillera Vilcanota, 82 of 257 glaciers (which represented 11.2% of the 1975 glacierised area) disappeared completely over the study period–most notably a cluster of glaciers to the south-east of the Quelccaya ice cap (Figure 8). However, deglaciation is occurring at differing rates between regions. Between 1988 and 2020, the Cordillera Vilcabamba lost 34% of its glacial area, compared to 45% for the Cordillera Vilcanota between 1984 and 2020. This is likely due to the differing morphology of glaciers between the regions–glaciers in the Vilcabamba are generally smaller and at lower elevations (Drenkhan et al., 2018) in comparison to the Vilcanota, which means less relative area is exposed to rapid melt. A differing climate setting between the regions is also likely to affect area change: the Cordilleras Vilcabamba and Urubamba are generally wetter than the Vilcanota (Kaser, 2001; Fyffe et al., 2021) though, unlike Fyffe et al. (2021), we did not observe precipitation to be a statistically significant driver of ice area change. This is likely due to the temporal resolution of our dataset; we examined annual precipitation and often with gaps between years in the dataset. Fyffe et al. (2021) studied the influence of precipitation at daily to weekly timescales, and found high variability in melt due to these precipitation dynamics. Further work is needed at a higher temporal resolution to better understand the role precipitation plays in regulating ice area change over long time scales.

Glacier thinning has been consistent throughout the two decadal intervals studied. It is difficult to make an authentic comparison of our mass balance rate (−0.477 ± 0.067 m w.e. yr⁻¹) to other estimates as differing proportions of the regions are studied. Using thickness-volume scaling, Salzmann et al. (2013) calculated an approximate thinning of −0.39 m a⁻¹ from 1985–2006 (Seehaus et al., 2019), which would suggest that thinning has accelerated in recent decades when compared to our estimate of −0.563 m a⁻¹ from 2000–2020. Limited field-based point measurements of mass balance over individual glaciers of the Cordillera Vilcanota have largely focussed on understanding the role micro-topography plays on the spatial differentiation of mass balance; ranging from −4 to 1.4 m w.e. yr⁻¹ across the surface of Suyuparina glacier in the Cordillera Vilcanota in 2013 (Molina et al., 2015). While we do not have exact location information on their stakes, our results for the 2010–2020 time interval show ranges of −4.79 to 2.28 m w.e. yr⁻¹ over the same glacier. Our estimates of mass balance are slightly larger than the −0.42 ± 0.24 m w.e. yr⁻¹ calculated by Dussaillant et al. (2019), which suggests this region is also thinning faster than their Andes-wide average. A Peru-wide calculation of glacier thinning from 2000–2016 was estimated to be −0.359 ± 0.068 m a⁻¹ (Seehaus et al., 2019), which suggests that thinning in the Cordillera Vilcanota is generally faster than the Peru-wide average.

In this study, we did not explicitly consider the change in glacial lakes across these regions and how this affects glacier loss. Previous studies (e.g. Drenkhan et al., 2018; Wood et al., 2021), have shown glacial lake expansion across each of these Cordilleras in line with climatic changes. The role of glacial lakes may be important in evaluating future change for tropical mountain glaciers, as lakes have been shown to accelerate ice recession when in direct contact (King et al., 2019). As more lakes develop in overdeepenings left by glacial recession (Colonia et al., 2017), this opens up the opportunity for accelerated ice loss. Similarly, in areas where glacial lakes do currently exist and ice is currently in contact with these lakes (such as the Western face of the Quelccaya ice cap), there may be a slow-down in areal recession in future decades as ice retreats up-slope and out of these lakes. This is particularly important in the Cordillera Vilcanota where 117 new large (>10,000 m²) lakes are projected to develop by 2050; the greatest of any Peruvian Cordillera (Guardamino et al., 2019), which is likely to accelerate glacier recession.

Future projections of glacier change in the Andes are stark. Even under a low emissions pathway, Schauwecker et al. (2017) suggest a rise in the freezing line of the Cordillera Vilcanota of 230 ± 190 m by 2100, which would lead to greater committed ice loss. The close relationship between temperature and ELA, projected 2.5–4.5°C warming over these glaciers by 2100 (Yarleque et al., 2018; Fyffe et al., 2021; Veettil and de Souza, 2017), and presence of low-lying glacier tongues across the regions suggest there will be significant melt of the Southern Peruvian Andes this century. Studies projecting future increase in ELA over the Vilcanota and Quelccaya (Schauwecker et al., 2017; Yarleque et al., 2018) also suggest this rate is likely to rise sharply in the coming decades. Low-lying regions within the Cordillera Vilcanota (South-East and North; Supplementary Figure S3) will likely experience the most glacier loss in the near future. Predictions of remaining glacial ice in 2100 for the Cordillera Vilcanota range from a pessimistic 13 km² to an optimistic 155 km² (5.4–63.0% of 2020 ice extent) (Schauwecker et al., 2017).

For people living in close proximity to this glacial ice, this represents a major loss of a freshwater resource. Following the approach of Drenkhan et al. (2018), we can contextualise this by noting that in our 2000–2020 study of volume change, the loss of 3.44 km³ of water from the Cordillera Vilcanota represents ∼84 years of water supply for the city of Cusco. Up to 33% of the population of Cusco rely on glacial meltwater during dry seasons (3% annually) (Martínez et al., 2020). Future ice melt is also unlikely to be even. Key glaciers supplying the Laguna Sibinacocha would lose over 90% of their glacierised area by 2100 (Drenkhan et al., 2018), while the Quelccaya ice cap would likely be below the ELA by mid-century (Yarleque et al., 2018). The low-lying Cordillera Urubamba is particularly at threat of severe ablation over the coming decades owing to the high rate of ELA change, posing an imminent risk of loss of water for irrigation, hydropower, and a drinking source in dry seasons (Chevallier et al., 2011). Ice melt is also unlikely to be even in time. While few studies have examined these specific Cordilleras for the timing of “peak water” (a maximum river discharge from the glacial melt contribution), indications from other Cordilleras of Peru show that peak water is likely to pass in the near-future (Drenkhan et al., 2014; Mark et al., 2017; Huss and Hock, 2018; Drenkhan et al., 2019).

An immediate concern is also the increased risk of outburst flood events and hazards to communities downstream of these destabilising glaciers (Vilca et al., 2021). Urgent adaptation measures, such as alternative agricultural practices (e.g. sustainable irrigation systems) and glacial lake outburst flood monitoring systems, are needed across the Southern Peruvian Andes to minimise the impact of this loss on communities (Veettil and Kamp, 2019). This must be integrated with a social perspective to recognise the cultural and spiritual loss marked by the recession of these glaciers, and connection of indigenous communities to the ice (Motschmann et al., 2020).

The drivers of ice loss across all three regions are largely due to warming maximum air temperatures (which correlated with ice area loss and ELA changes in all regions) and topographic controls, most notably, elevation. This is as expected; freezing levels are rising in line with temperature rises and low-lying glacier tongues, many of which are low-gradient and thus susceptible to even minor temperature shifts, are now exposed to melt processes (Thompson et al., 2011). However, we did not have the data available to explore all the predicted key climatic drivers of ice loss, such as atmospheric humidity and radiation (Kaser, 1999; Vuille et al., 2008). Furthermore, the role of ENSO as a driver of tropical glacier mass loss is becoming realised over individual glaciers across the Andes (Maussion et al., 2015; Veettil et al., 2016), but we could not explore this phenomenon within the available datasets. This signal is thought to be lesser over the Cordilleras Vilcanota, Vilcabamba, and Urubamba as they receive moisture from the Amazon during La Niña; thereby dampening the ENSO signal (Perry et al., 2014). This highlights the need to complement remote sensing observations with in-situ data to bridge such gaps.

The optical imagery from Landsat and ASTER used in this study provides useful insights into the dynamics of changes to mountain glaciers, and the long archives of both have been particularly beneficial to generate these results. However, their limitation is highlighted by the lower numbers of usable images we were able to acquire over the Cordilleras Vilcabamba and Urubamba and steep terrain. Mass balance estimates over the Peruvian Andes are few and can be improved using satellite data (Vuille et al., 2018), but the level of uncertainty driven by lack of usable DEM data from ASTER images remains too high. We have demonstrated in this study that the 100 m resolution Hugonnet et al. (2021) dataset provides comparable estimates of mass balance over small mountain glaciers as a higher resolution 30 m product. Radar satellites are able to penetrate cloud, which overcomes the primary limitation, but are challenged by steep topography, which often limits their usefulness to glaciers far larger than those in the Southern Peruvian Andes (Lee et al., 2013). Furthermore, with the steady decay of ASTER as the Terra satellite drifts further out of orbit, and with no current plans for an alternative non-commercially operated stereo satellite, there is an urgent need to find alternative sources to derive regional-level and multi-temporal digital elevation data.