4.1 Glacier inventory

Since the manual delineation of glacier extents is laborious, time-consuming, and subjective, several methods have been developed to automatically map glacier outlines based on multispectral images (Veettil and Kamp, 2017). The most widely used and robust approaches are the computation of the normalized difference snow index (NDSI) or the band ratio (BR) and the application of a threshold value to differentiate between on- and off-glacier areas (GLIMS algorithm working group, 2004; Paul et al., 2013). In this study, we first used the NDSI to classify glacier areas and combined it with BR information to improve the mapping in areas affected by shadows. NDSI maps generated from top-of-atmosphere reflectance values show a better performance than NDSI maps based on digital number values. For the BR computation, digital number values are taken. The threshold value selection is supported by high-resolution satellite imagery (Google Earth) from the respective dates. A NDSI threshold value of 0.8 is selected, which is higher than the thresholds of 0.5–0.6 applied by other studies in this region (e.g. Silverio and Jaquet, 2005; Veettil et al., 2017). This offset might be induced by the application of top-of-atmosphere reflectance values instead of digital number values. The threshold for the BR data is set to 1.7 for Landsat 5 TM and 1.5 for Landsat 8 OLI data. Finally, polygons of the glacier outlines are generated from the computed glacier masks.

The detection of the debris-covered glacier termini extents in the Cordillera Blanca is difficult using multi-spectral imagery. Therefore, we generated SAR coherence maps from repeat-pass SAR acquisitions to distinguish the debris-covered ice from the surrounding ice-free areas (Atwood et al., 2010; Lippl et al., 2018). The surface structure of the debris-covered glacier areas changes over time due to the dynamics and melting of the underlying ice. This leads to a temporal decorrelation of the backscattered SAR signal of repeat-pass SAR imagery and thus to lower coherence compared to the surrounding ice-free areas. This difference in coherence facilitates the delineation of the debris-covered ice areas. Data from the Sentinel-1 mission are used to map the debris-covered areas in 2016. No suitable repeat-pass SAR acquisitions are available for the Cordillera Blanca in 2013. Thus, we had to rely on TerraSAR-X and Sentinel-1 data from 2014 to map the outlines of the debris-covered glacier tongues. In 2000 (and ±1 year), only repeat-pass SAR data are available from the European Remote Sensing (ERS) satellite with a repetition cycle of 1 d. Due to this short temporal baseline, a separation between debris-covered ice and surrounding ice-free areas is unfeasible. Consequently, we combined our outlines from 2000 with the manually delineated debris-cover masks available from the Global Land Ice Measurements from Space (GLIMS) database from 2003 based on SPOT imagery (Racoviteanu, 2005).

The catchment discriminations of the RGI are applied to split the resulting polygons into individual glacier basins. In the next step, the glacier inventories are visually inspected and misclassified areas are manually corrected. These manual corrections are supported by high-resolution imagery from the respective years (Google Earth). According to the RGI 6.0 Technical Report (RGI Consortium, 2017), topographic parameters (minimum, maximum and median elevation, mean slope, and aspect) of the individual glacier basins are computed using the void-filled SRTM DEM as an elevation reference. Finally, the areas (S) of the complete inventories and of each glacier are measured in UTM projection (UTM Zone 18S for subregion R1 and UTM Zone 19S for subregion R2 and R3). The uncertainties of the area gauging (δ_S) are calculated following the approach of Malz et al. (2018) based on an error evaluation of 3 % for alpine glacier outlines derived from Landsat images (Paul et al., 2013). This estimate is scaled by the area-to-perimeter ratio of the studied subregion compared to the area-to-perimeter ratio of Paul et al. (2013) in order to account for differences in the shape of the glacierized areas.

4.2 Elevation change

Surface elevation change information is computed by differencing DEMs from SRTM and TDX data. DEMs are derived from the bistatic TDX imagery following the differential interferometric approach (e.g. Malz et al., 2018; Seehaus et al., 2015; Vijay and Braun, 2016), which is briefly summarized in the following.

First, acquisitions from the same relative orbit and date are concatenated in the along-track direction. A differential interferogram is computed using the void-filled SRTM DEM as elevation reference. In the next steps, the interferogram is filtered and unwrapped by applying the branch cut and minimum cost flow algorithm, and the unwrapped differential phase is transferred into differential elevations. Subsequently, the topographic information of the SRTM DEM is added to obtain absolute height information and finally the product is geocoded and orthorectified. The DEMs are visually checked for phase jumps and the best results of both phase-unwrapping methods are selected for further processing. Areas affected by remaining phase jumps are masked out.

The TDX DEMs need to be precisely horizontally and vertically coregistered to the respective reference DEM (SRTM for 2000, TDX for 2013) in order to accurately map elevation changes on the glacierized areas. Figure S2 illustrates the applied processing chain used to perform this coregistration. First, smooth stable reference areas are defined by masking out vegetation, water, and glacier areas. The vegetation and water masks are derived from region-wide cloud-free Landsat 8 mosaics and using a normalized difference vegetation index (NDVI) threshold of 0.3 and a normalized difference water index (NDWI) threshold of 0.1. Additionally, a slope threshold of 15^∘ (of the respective reference DEM) is applied. Thereafter, the TDX DEMs are bilinearly vertically corrected for offsets to the reference DEM, which are measured for the defined stable regions. Subsequently, a horizontal coregistration between the reference DEM and the TDX DEMs is carried out following the widely used approach of Nuth and Kääb (2011). Afterwards, a second bilinear vertical coregistration of the TDX DEMs to the reference DEM is run to reduce any biases that remain. Finally, the coregistered TDX DEMs are merged to a regional DEM mosaic, which include a date stamp for each grid cell.

To obtain elevation change rates Δh∕Δt of the respective study periods, the SRTM DEM and the TDX DEM mosaics are differentiated. The mean date of the 11 d SRTM mission (16 February 2000) is assigned to the SRTM DEM. Since data voids in the SRTM DEM are filled with data from other sources (no date information available), the non-SRTM data values are masked out using the coverage information provided by LP DAAC NASA. Numerous studies have revealed that the glaciers in Peru are in general retreating (Sect. 1). Thus, the glacier inventory from the beginning of the respective observation period is employed to create surface elevation change maps for on- and off-glacier areas. The average regional and glacier-wise elevation change rates are obtained by integration of Δh∕Δt over the respective areas. Slopes steeper than 50^∘ are rejected (5.7 % of the glacier area in 2000) since major ice aggregation is quite unlikely there (avalanche slopes, backed up by field observations) and DEMs are less accurate on these steep slopes (Toutin, 2002). To account for data voids in the elevation change fields on glacierized areas, the measured Δh∕Δt values are area weighted based on the hypsometric area distribution using 100 m elevation bins in order to calculate regional mean values. This is one of the recommended methods to obtain regional estimates and provides reliable results for datasets with up to 60 % voids (McNabb et al., 2019).

The average Δh∕Δt of individual glaciers is calculated according to McNabb et al. (2019) by using elevation bins of 10 % of the glacier elevation range, if it is < 500 m, and bins of 50 m for glaciers with elevation ranges > 500 m. A coverage of more than two-thirds of the elevation bins and < 60 % data gaps is used as a criterion for exclusion. Voids in the hypsometric Δh∕Δt distribution are filled by applying a 3rd-order polynomial fit. Outliers in the respective elevation bins (regional and glacier-wise analysis) are sorted out using 3 times the normalized median absolute deviation (NMAD) (Brun et al., 2017). For all hypsometric analyses of elevation changes (SRTM to TDX, TDX to TDX), the void-filled SRTM DEM is utilized.

Figure 2Off-(red) and on-glacier (light blue) area and off-glacier elevation change (blue dots) distributions as a function of slope in subregion R1 for the period 2000–2016. Error bars represent NMAD of Δh∕Δt values in the individual slope interval. The dotted line indicates the applied slope threshold (see Sect. 4.2). Glacier area measurements are based on the glacier outlines from 2000. Note: for better representation, on-glacier areas are scaled by a factor of 10. Plots for other subregions are provided in the Supplement.

Download

The uncertainties of the generated elevation change rates are assessed by evaluating the elevation change rates on non-vegetated stable off-glacier areas (water and vegetation masks; see above). The lowest and highest 2 % quantiles of the change rates are rejected to suppress the impact of processing artefacts and outliers. To account for the dependency of the offsets on the slope (Figs. 2 and S3 and S4), the deviations are binned in slope intervals of 5^∘. Remaining outliers are removed by employing a 3*NMAD filter for each slope bin. Finally, the area-weighted standard deviations σ_AW based on the offsets in off-glacier areas and the slope distribution in glacier areas are calculated.

Since we integrate elevation change information over the glacierized area, spatial autocorrelation of the elevation change fields must be considered in the accuracy assessment. We estimated the uncertainty of the computed average elevation change rates (σ_Δh∕Δt) according to the approach of Rolstad et al. (2009):

where $A_{\mathrm{cor}}=\mathit{\pi }\cdot d_{\mathrm{cor}}$ is the correlation area, A_gl is the analysed glacier area, and σ_AW is the assessed accuracy of the elevation change rates (explained in the next paragraph). The correlation length (d_cor) is obtained by generating semivariograms with 100 000 random samples of Δh∕Δt values on the off-glacier areas. A binning in 30 m distance intervals and a maximum distance of 20 km are applied. Spherical semivariogram functions are fitted to the data and an average correlation length of 387 m results from the analysed elevation change fields. Equation (1) is applied for each continuous glacierized area (icecap or connected glaciers) and the area-weighted average of the individual ice-covered areas is taken as the region-wide σ_Δh∕Δt.

The hypsometric extrapolation of elevation change information leads to an additional uncertainty that is hard to quantify. We employed the approach of Berthier et al. (2014). A scaling factor (we selected a factor of 2) is applied to σ_Δh∕Δt for the area fraction with hypsometric extrapolation of Δh∕Δt, in order to obtain the uncertainty of our region-wide average elevation change rates δ_Δh∕Δt.

4.3 Mass balances

The geodetic mass balances ΔM∕Δt of the analysed regions are computed according to Fountain et al. (1997) by multiplying the integrated elevation change rates by the average ice density (volume-to-mass conversion factor). We applied two density scenarios. The first scenario follows the suggestion of Huss (2013) for alpine glaciers and we applied an average density (ρ) of 850 kg m⁻³. For the second scenario, two different conversion factors of 600 and 900 kg m⁻³ are applied for ablation and accumulation areas, respectively (Kääb et al., 2012; Gardelle et al., 2012). The average equilibrium line altitude (ELA) (see Table S3 and below) of each subregion is used to distinguish between both glacier sections. As revealed by Huss (2013), the conversion factor can vary strongly. Following his findings, we applied an uncertainty of ±300 kg m⁻³ for the conversion factor for observation intervals shorter than 10 years and mass budgets lower than ±200 kg m⁻² a⁻¹ and an uncertainty of ±60 kg m⁻³ for all other cases (see Braun et al., 2019). Mass budgets of individual glaciers are only calculated for the constant density scenario.

In order to estimate the accuracy of the geodetic mass balances, the following error contributions are considered:

• accuracy of the elevation change rates δ_Δh∕Δt

• accuracy of the glacier outlines δ_S

• uncertainty of the applied average ice density δ_ρ

• potential bias due to different SAR signal penetration V_pen∕Δt.

This leads to the following formula to calculate the error of ΔM∕Δt:

The assessments of δ_S and δ_Δh∕Δt are depicted in Sect. 4.1 and 4.2. The uncertainty contribution by potentially different SAR signal penetration in the glacier surfaces is evaluated following the approach of Malz et al. (2018). No difference in the SAR signal penetration in the ablation areas below the equilibrium line altitude (ELA) is assumed. These areas experience melt throughout the year in the tropical Andes (Kaser, 2001) and differences in the radar signal penetration in wet glacier surfaces are small (Casey et al., 2016; Rossi et al., 2016; Ulaby et al., 1984). A linear increase in the penetration depth bias towards 5 m between the ELA and the highest peaks is employed in the accumulation areas to calculate V_pen. Depending on the applied scenario, a volume-to-mass conversion factor of 850 or 600 kg m⁻³ is taken. The late summer snow line altitude (SLA) derived from optical satellite images is a good proxy for the ELA in the tropical Andes (Rabatel et al., 2012). Thus, the average ELA in the study period for the subregions is estimated based on published ELA and SLA values. Table S3 provides an overview of the considered ELA and SLA values. Uncertainties in the estimated ELA positions affect the penetration depth bias estimation and thus the error budget. A sensitivity analysis of an ELA mismatch on the error budget is provided in Braun et al. (2019). Specific mass balances are calculated according to the UNESCO Glossary of Glacier Mass Balance (Cogley et al., 2011) by applying the average glacier area of each observation period (mean glacier extent with slopes < 50^∘).

Figure 3Polar plot of relative area changes (dot colour) in subregion R1 in the period 2000–2016 of individual glaciers. Dot size: glacier size in 2000; radius: median elevation; orientation: mean aspect. Red circle: equilibrium line altitude (ELA); see also Table S3. Plots for other subregions and time intervals are provided in the Supplement.

Download

5.1 Glacier inventory

The glacier outlines from 2000, 2013, and 2016 of some example mountain ranges from all three subregions are shown in Fig. S5. A clear recession of the glaciers in all subregions is obvious. The measured extents of the glaciers for all subregions and time steps are presented in Table 1. Since cloud and snow cover do not allow for a complete coverage of subregions R1 and R2 in 2013 and no suitable SAR data cover the debris-covered ice areas in the Cordillera Blanca in 2013, a small fraction of the glacier area (subregion R1: 2.1 %; subregion R2: 4.4 %) is mapped using imagery from 2014. For subregion R1, a nearly complete mapping (99.5 %) of the glacier areas in 2014 is possible. The obtained average area change rate between 2013 and 2014 at glaciers delimited in both years is employed to correct the area change measurement in subregion R1 for 2013 at glaciers only delimited in 2014 and vice versa. In subregion R2, the change rate between 2013 and 2016 is applied to perform a similar correction of the area change measurements in 2013. In total, ice covered areas of 1916.6±48.3 km² in 2000, 1571.9±43.1 km² in 2013, and 1368.1±45.5 km² in 2016 are estimated for Peru. Based on the catchment divides of the RGI6.0, a reduction in the number of glaciers from 1973 in 2000 to 1803 in 2016 is observed. The variations in the glacier count and area in the subregions and in the whole country are summarized in Tables 1 and 2. It should be noted that the uncertainty of the area changes (Table 2) results from the uncertainty of the individual inventories (quadratic sum), assuming independence of the individual area measurements. The revealed area changes of individual glaciers in the different subregions and study periods are correlated with topographic parameters (glacier area, median elevation, mean aspect) and plotted in Figs. 3 and S6–S13.

Table 1Measured glacier extents for different years and regions.

S: measured glacier area. δ_S: uncertainty of measured glacier area. n: number of glacier catchments (delineations based on RGI 6.0).

Download Print Version | Download XLSX

5.2 Elevation changes

The obtained unfiltered elevation change rates on ice-covered areas at some example mountain ranges are illustrated exemplarily in Figs. 4 and 5 for the periods 2000–2013 and 2013–2016, respectively. In the interval 2000–2013, a clear thinning of the lower glacier parts in the Cordillera Vilcanota and Apolobamba is observed. The glaciers in the Cordillera Blanca and Huayhuash show a more balanced pattern and Coropuna experienced thinning throughout most of its ice cap. An increase in the surface lowering rates at most mountain ranges is obvious in the interval 2013–2016. Only the ice caps in subregion R3 show reduced lowering rates. Since the TDX data in 2013 do not cover the whole glacier area in subregion R3, the voids are filled with TDX data from 2014. The glacier area covered by Δh∕Δt measurements using TDX data from 2014 amounts to only 1.4 % (glacier outlines from 2000). The impact of this void filling is considered negligible since the analysis uses change rates. The average measured and extrapolated surface elevation change rates of all subregions for different observation periods are listed in Table 2. The fraction of glacier area covered by Δh∕Δt measurements is lowest in subregion R1, when using the SRTM DEM as a reference (Table 2). Large data voids in the SRTM data lead to a partial coverage of only 46 % of the ice-covered area in subregion R1 in 2000 (61 % in subregion R2, 89 % in subregion R3). The number of measurements on the glacierized areas clearly increases when deriving the elevation change solely from TDX datasets (Table 2). In the period 2013–2016, a coverage of up to 80 %, 69 %, and 89 % is obtained in subregions R1, R2, and R3, respectively. Since the differential InSAR approach is used to generate the TDX DEMs, the proportion of areas affected by phase jumps in the unwrapped interferogram is small (∼1 % of the total glacier area). The major area affected by phase jumps (at the Cordillera Vilcanota in 2016) is highlighted in Fig. 5. Figures 4 and 5 indicate that the elevation change rates are altitude dependent. The hypsometric distributions of the measured elevation change rates and the respective ice-covered areas are plotted in Figs. 6, S14, and S15 for each subregion in the interval 2000–2016. Surface lowering is found on areas below 5700 m in subregion R1, 5800 m in subregion R2, and 5900 m in subregion R3. Considerable positive elevation changes are observed in the higher reaches of the glaciers in subregion R1. The obtained elevation change rates of the glaciated areas are available via the PANGAEA database (https://doi.org/10.1594/PANGAEA.906211, Seehaus et al., 2019).

Figure 4Surface elevation changes in the period 2000–2013 at mountain ranges in all three subregions (unfiltered). The frame colour of the panels indicates the subregions (see Fig. 1). Background: SRTM DEM hillshade © NASA.

Figure 5Surface elevation changes in the period 2013–2016 at mountain ranges in all three subregions (unfiltered). The frame colour of the panels indicates the subregions (see Fig. 1). White dashed line: masked out areas affected by phase jumps in the unwrapped interferogram. Background: SRTM DEM hillshade © NASA.

Table 2Measured area, surface, and mass changes for different periods and regions. ^a Constant volume-to-mass conversion factor of 850 kg m⁻³. ^b Volume-to-mass conversion factor of 600 kg m⁻³ and 900 kg m⁻³ for region below and above the ELA (Table S3), respectively. ^c Mean glacier area of observation interval is used to calculate specific mass balances (see Sect. 4.3)

dt: observation period (mean time lag between DEM dates). dS: glacier area change. dS∕dt: glacier area change rate (mean time lag between glacier inventories is used). dh_M∕dt: average measured surface lowering rate. dh_E∕dt: average extrapolated surface lowering rate. S_M: fraction of glacier area covered by dh_M∕dt measurements (below slope threshold). ΔM∕Δt: mass balance (average and specific). ΔM: total mass change in the observation period.

Download Print Version | Download XLSX

5.3 Mass balances

The geodetic mass balances derived from the elevation change information of the individual subregions are listed in Table 2. In total, a mass loss of $-\mathrm{7.62}\pm \mathrm{1.05}$ or $-\mathrm{9.18}\pm \mathrm{1.10}$ Gt is found for Peru in the period 2000–2016, depending on the applied density scenario (Sect. 4.3). The mass loss rates show an increase between the observation periods 2000–2013 and 2013–2016. Only in subregion R3 is a slight reduction of the mass loss observed. The most prominent increase in glacier wastage is revealed in subregion R1. The mass budget changed from nearly stable conditions of 45±97 kg m⁻² a⁻¹ in 2000–2013 to a specific mass loss rate of $-\mathrm{752}\pm \mathrm{452}$ kg m⁻² a⁻¹ in 2013–2016 (second density scenario). The computed mass balances of individual glaciers in the different subregions and study periods are correlated with topographic parameters (glacier area, median elevation, mean aspect) and plotted in Figs. 7 and S16–S20.

Figure 6Hypsometric distribution of measured (red bars) and total (light blue bars) glacier area with elevation in subregion R1 in the interval 2000–2016. Blue dots represent the mean Δh∕Δt value in each elevation interval. Error bars indicate NMAD of Δh∕Δt for each hypsometric bin. Grey areas mark the lower and upper 1 % quantile of the total glacier area distribution. Black dashed line: mean glacier elevation; red dashed line: equilibrium line altitude (ELA); see also Table S3. Area measurements are based on the glacier outlines from 2000, considering only regions with slopes below applied slope threshold (50^∘; see Sect. 4.2). Plots for other subregions are provided in the Supplement.

Download

Figure 7Polar plot of specific mass balance (dot colour) of individual glaciers in subregion R1 in the period 2000–2016 of individual glaciers. Dot size: glacier size in 2000; radius: median elevation; orientation: mean aspect. Red circle: equilibrium line altitude (ELA); see also Table S3. Note: only glaciers with elevation change information > 50 % are included. Plots for other subregions and time intervals are provided in the Supplement.

Download

6.1 Glacier retreat

We observed a dramatic recession of the glaciers throughout Peru of −29 % ($-\mathrm{548.5}\pm \mathrm{65.7}$ km²; −1.8 % a⁻¹) between 2000 and 2016. Our total mapped glacier extent of 1368.1±44.5 km² in 2016 is comparable to the reported coverage of 1298.6 km² by the recent Peruvian glacier inventory (UGRH, 2014), considering that UGRH did not include the glacier areas of the Bolivian Cordillera Apolobamba (∼70 km² ). The glacier area mapped in the RGI6.0 amounts to 1602.96 km², which is in the range of our measurements for 2000 and 2013. A direct comparison is complex since the RGI6.0 is a blended product with outlines discriminated between 2000 and 2009 in Peru. However, visual inspection of our outlines and the RGI6.0 revealed that numerous small glaciers (especially in the southern section of subregion R1 and throughout the whole subregion R3) mapped in the RGI6.0 are artefacts, which are caused most probably by temporary snow cover. This issue has already been mentioned in the RGI Technical Report (RGI Consortium, 2017). Thus, our glacier inventory of Peru in 2000 has 350 fewer features than the RGI6.0. By contrast, the Peruvian glacier inventory (UGRH, 2014) lists 2679 glaciers. This difference can be explained by the different basin delineations applied.

The comparison of our area measurement of 1916.6±48.3 km² in 2000 and the first Peruvian glacier inventory (Hidrandina SA, 1989) in 1970 (2041.85 km²) results in a retreat of −7 % (−139.9 km²; 0.2 % a⁻¹). However, not all cordilleras were mapped completely by the first Peruvian glacier inventory. When considering only the glacierized cordilleras which were completely mapped in the first Peruvian glacier inventory (UGRH, 2014), the area changes amount to −23 %. Figure 8 illustrates the temporal evolution of the glacier area changes and highlights that Peru's glaciers have experienced long-term shrinkage since 1970 with increasing rates in recent years. This observed trend is in accordance with the findings of previous studies at individual mountain ranges. Silverio and Jaquet (2017) summarized area measurements of several studies in the Cordillera Blanca (subregion R1) and revealed an increase in the loss rate from about −5 km² a⁻¹ between 1971 and 1996 towards −23 km² a⁻¹ between 2014 and 2016. Burns and Nolin (2014) reported a ∼3.5 times higher area loss rate in 2004–2010 compared to 1970–2003 in the Cordillera Blanca. Area loss rates of ∼1.8 % a⁻¹ (> 50 %) since 1975 are reported by the Instituto Geofisico del Peru (2010) and López-Moreno et al. (2014) for the Cordillera Huaytapallana, which is higher than our finding of −1.1 % a⁻¹ in subregion R1 in the period 2000–2013. However, their study period includes strong El Niño events in the 1980s and 1990s (Fig. 8), which typically lead to increased glacier melt in the tropical Andes (Wagnon et al., 2001). Retreat measurements in subregion R2 are carried out at the Cordillera Vilcanota (Hanshaw and Bookhagen, 2014; Salzmann et al., 2013; Veettil and Souza, 2017), the Cordillera Apolobamba (Cook et al., 2016; Veettil et al., 2017a), the Cordilleras Carabaya, Urubamba, and Vilcabamba (Veettil et al., 2017d), and at glaciers draining into the Vilcanota–Urubamba basin (Drenkhan et al., 2018). All studies revealed significant glacier shrinkage since the 1980s with retreat rates in the range of −0.9 to −1.7 % a⁻¹. These findings are comparable with our observed rate of −1.3 % a⁻¹ in subregion R2 in the period 2000–2013. Drenkhan et al. (2018) revealed reduced shrinkage rates for 2010–2016 compared to 2004–2014. This is contradictory to our observed strong increase in glacier recession after 2013. They analysed only glaciers in the Vilcanota–Urubamba basin, of which most face south. We measured the highest retreat for glaciers facing north in subregion R2 in the period 2013–2016 (Fig. S22) and only low retreat for south-facing glaciers. Thus, different spatial extents and representativeness of topographic attributes in the analysed regions leads to the mismatch of the observed retreat trends. Moreover, this suggests that the representativeness of topographic settings needs to be considered when doing region-wide upscaling of sampled glacier change measurements. At the ice cap of the Coropuna volcano (subregion R3), various ice coverage estimates are available going back to 1955 (Peduzzi et al., 2010; Racoviteanu et al., 2007; Silverio and Jaquet, 2012; Ubeda, 2011; Veettil et al., 2016). The average loss rate of ∼1.5 % a⁻¹ (1955–2015) is lower than our estimated rate of −2.3 % a⁻¹ for subregion R3 in the period 2000–2013. In addition to differences in the study periods, we attribute this deviation to the fact that our estimate for subregion R3 includes numerous small glaciers located at lower altitudes (Fig. S10). These glaciers are in general more sensitive to climate change (Francou Bernard et al., 2003; Vuille et al., 2008) in comparison to the more elevated glaciers of the Coropuna volcano. Moreover, Veettil et al. (2016) discovered an increased retreat and uplift of the SLA after ∼2000 at Coropuna's ice cap, supporting our findings.

Figure 8Temporal evolution of (a) relative glacier area changes, (b) observed specific mass balance (spMB; colour shaded areas represent the uncertainty of the corresponding measurements), and (c) Oceanic Nino Index (ONI) in the period 1970–2016. ^* Area measurement in 1970 taken from the first Peruvian glacier inventory (Hidrandina SA, 1989) and only mountain ranges with complete coverage by the inventory in 1970 are considered in the area change computation (see also Sect. 6.1).

Download

For the period 2013–2016, we discovered a 4 times higher countrywide retreat rate compared to the period 2000–2013. The strongest increase is found in subregion R2, whereas in subregion R3 (only ∼5 % of total glacier area) the glacier area remained quite stable after 2013. The increased retreat rates in subregions R1 and R2 after 2013 can be attributed to the strong ENSO activities in the years 2015–2016. An average ONI of 0.42 is reported for 2013–2016 and the maximum ONI of 2.6 in December 2015 indicates distinct El Niño conditions. On the other hand, an average ONI of −0.17 is revealed for 2000–2013, indicating that La Niña conditions dominated this period. Since El Niño periods typically lead to increased glacier wastage in the tropical Andes (e.g. Vuille et al., 2008; Wagnon et al., 2001) (Sect. 2), our observed increased shrinkage after 2013 can be attributed to the ENSO activities in this period. This pattern also fits the findings by Morizawa et al. (2013) (Condoriri Glacier, Bolivia) and Silvero and Jaquet (2017) (Cordillera Blanca, subregion R1). Both studies reported enhanced recession during El Niño events and even area gains during La Niña epochs. The latter study also discovered a linear relation between glacier retreat and ONI (R²=0.8) and reports a change rate of −5 % a⁻¹ for 2014–2016 that is equal to our change rate at subregion R1 for 2013–2016. Moreover, the glaciers in the tropical Andes are on average the thinnest worldwide (Farinotti et al., 2019). Therefore, the increased melt will lead to more pronounced changes in the glacier area compared to other glacier regions. The stagnation of the glacier retreat in subregion R3 in the interval 2013–2016 is difficult to explain. During the strong El Niño in 1997/98, increased precipitation was observed at the Coropuna volcano in subregion R3 (Herreros et al., 2009; Silverio and Jaquet, 2012), whereas El Niño usually leads to reduced precipitation rates. Veettil et al. (2016) reported positive precipitation anomalies after 2011 and clearly negative anomalies in the period 2009–2011 for the Coropuna volcano. However, the total precipitation data from ERA-Interim indicate lower precipitation rates after 2013 (Fig. S24) and do not clearly indicate an increase in precipitation during El Niño 1997/98. Since the small glacier areas in subregion R3 cover mainly volcano peaks, the revealed precipitation values from the spatially coarse ERA-Interim data do not necessarily reflect the local precipitation pattern at the prominent, high-altitude volcano peaks. Additionally, the mean glacier altitude in 2013 shifted ∼100 m above the ELA (Fig. S25), whereas the mean glacier altitude in 2000 was nearly similar to the ELA (Fig. S15). A significant temporary snow cover at the non-glacierized Cordillera del Barroso (subregion R3, close to the boarder to Bolivia) was observed during the dry seasons in 2015 and 2016 (Lèon et al., 2019), which fits to the suggested increased precipitation at high elevations after 2013. Moreover, snowfall events during the dry season strongly affect the glacier albedo and thus lead to reduced ablation. However, the temporary snow cover during the dry season in 2016 in subregion R3 might have also affected the mapping of the glacier outlines. Although only imagery with no or only minimal snow coverage is selected, it is quite likely that some remaining snow cover was located at the glacierized peaks, leading to a slightly larger glacier area in 2016 compared to 2013. This bias is not quantifiable and certainly within the range of applied uncertainty. Thus, we conclude that the glacier area did not significantly change in subregion R3 after 2013 and attribute this to the allocation of the remaining ice at higher altitudes and increased precipitation, especially during the dry season, even though a strong El Niño event occurred in this period.

The analysis of area fluctuations of individual glaciers revealed in all three subregions indicates higher recession for glaciers with lower median elevation and for small glaciers (Figs. 3 and S6–S13). This is in accordance with findings reported in previous studies (e.g. Kaser and Osmaston, 2002; Mark and Seltzer, 2005; Ramirez et al., 2001). Small glaciers have in general a more narrow altitude range compared to larger glaciers, which can maintain the ELA below the maximum glacier elevation. A rise in SLA (a proxy of the ELA in the tropical Andes, Sect. 4.3) is observed throughout the Peruvian Andes by various studies (Hanshaw and Bookhagen, 2014; López-Moreno et al., 2014; McFadden et al., 2011; Veettil et al., 2016, 2017c, d). This corresponds to our observed retreat pattern. Figures 3 and S6–S13 suggest that glaciers with slopes facing on average the south-south-west direction experienced in general higher relative retreats. The total amount of lost glacier area repeats this general pattern (Figs. S21–S23). This can be attributed to the fact that more low-lying, small glaciers with mean aspects facing southwards still exist (Figs. 3 and S6–S13). Higher retreat rates in general were observed for north-orientated glaciers before 2000 (Veettil et al., 2017a; Veettil, 2018), leading already to the disappearance of small north-facing glaciers before the start of our observation periods.

In total, we discovered that 177 glaciers had disappeared in our observation period, of which most disappeared after 2013 and were south facing (Table 1 and Figs. S11–S13). At the Artesonraju Glacier in the Cordillera Blanca (subregion R1), Vuille et al. (2018) projected an uplift of the ELA by 300–700 m until 2100 based on the CMIP5 scenarios RCP4.5 and RCP8.5. Thus, the proceeding climate change will lead to the further disappearance of numerous small low-lying glaciers in the tropical Andes within the next decades, as predicted by Ramirez et al. (2001) and Huh et al. (2017).

The gain and formation of proglacial lakes are consequences of glacier recession (Cook and Quincey, 2015) and increase the GLOF risk of downstream areas. Veettil et al. (2017a) discovered an increase in the number of glacial lakes from 697 to 903 in the Cordilleras Apolobamba and Carabaya between 1985 and 2015. In the Cordillera Vilcanota, Hanshaw and Bookhagen (2014) observed stable or increasing extents in 77 % of the lakes connected to glacial watersheds. Colonia et al. (2017) compiled an inventory of 201 potential future glacier lakes based on modelled glacier bed overdeepenings. Considering the revealed findings by this study as well as previously reported values and projections, a further region-wide monitoring of glacier retreat and lake development is highly advisable to identify potential GLOF risks in the tropical Andes, as also suggested by Cook et al. (2016).

6.2 Surface elevation changes and mass balances

The average countrywide glacier surface elevation change between 2000 and 2016 amounts to $-\mathrm{0.359}\pm \mathrm{0.068}$ m a⁻¹, which corresponds to a specific mass budget of $-\mathrm{296}\pm \mathrm{41}$ or $-\mathrm{357}\pm \mathrm{43}$ kg m⁻² a⁻¹ depending on the applied density scenario. The second density scenario leads to 8 %–21 % lower countrywide mass balances. As pointed out by Kääb et al. (2012), this scenario is suitable for glaciers with no dh∕dt due to ice dynamics and when dh∕dt is clearly driven by melt or increased accumulation. Since these conditions are typical for the Peruvian glaciers (see Sect. 1 and further down), we used the results of the second density scenario for further discussion and analyses. On the other hand, the results from the first scenario provide suitable information for comparison with other studies. Moreover, using the second scenario reduces the potential bias by the SAR signal penetration depth differences (Sect. 3) on the mass budget by ∼40 % compared to the first scenario. A computation of individual glacier mass balances using the second scenario is not performed due to the lack of ELA information on glacier scales.

The highest average surface lowering is revealed for subregion R1 in the period 2013–2016 (Table 2). However, the stable surface elevations before 2013 in subregion R1 suppress the long-term average value. The nearly balanced budget contradicts the observed glacier retreat of −15 % in subregion R1 in the interval 2000–2013 at first glance. The mean surface elevation change in the retreat areas amounts to −0.28 m a⁻¹, clearly indicating mass loss in the deglacierized areas. However, elevation gain is found at high altitudes. The more La Niña-like conditions of ENSO in the period 2000–2013 (Sect. 6.1) and an increase in precipitation in this region due to stronger upper-tropospheric easterlies (Schauwecker et al., 2017) have most probably led to higher accumulation rates. ERA-Interim reanalysis data also show an increase in total precipitation in this period, especially around 2007 (La Niña event). Thus, the accumulation gain in the upper reaches balanced the ice losses at the termini, even though temperatures increased (Schauwecker et al., 2017). In all subregions an increase in temperature is found in the reanalysis data for 2000–2013, but the strongest positive precipitation anomaly is found in subregion R1 (Figs. S24 and S26). This explains the mass losses in subregions R2 and R3 in this period.

Skin temperature was still above the long-term average in subregions R1 and R2 after 2013 (Fig. S26). The downward surface thermal radiation shows an increase in subregion R1 (Fig. S27), whereas total precipitation decreased in subregion R2 and remained nearly stable in subregion R1 (Fig. S24). Those climatic variations enhance the ablation and facilitate a positive feedback that further increases the glacier melt. The higher temperature as well as reduced and delayed precipitation, which are typical during El Niño, lead to liquid precipitation in the ablation regions and a reduced glacier albedo, enhancing the short-wave radiation absorption (e.g. Vuille et al., 2018; Maussion et al., 2015). Thus, the climatic variations explain the more negative mass balances in both subregions in the period 2013–2016, which also correlate with the strong El Niño activities in this interval (Fig. 8).

Only at subregion R3 did the thinning rate reduce after 2013, although El Niño conditions dominated. We attribute this, like the stable glacier area, to higher precipitation rates, especially during the dry season, at the ice-capped volcanos and the allocation of the remaining ice masses at high elevations (see Sect. 6.1 for more details). Moreover, the revealed spatial pattern of mass balances after 2013, with the highest mass loss rates at the northernmost subregion R1 (Table 2), matches the observed trend of higher river runoff towards northern Peru during strong El Niño events (Casimiro et al., 2012, 2013).

The analysis of the mass balance of individual glaciers and topographic parameters in subregions R2 and R3 reveals a trend towards higher mass losses for glaciers facing north-north-east in the period 2000–2016 (Figs. 7, S16, and S17). This trend agrees with observations by Soruco et al. (2009a) in the Bolivian southern wet outer tropics. In the interval 2000–2016, no clear dependency of the specific glacier mass balances on the median elevation or aspect is obvious in subregion R1 (Fig. 7). However, Fig. 5 suggests a trend towards higher surface lowering rates on the western slopes of the Cordillera Blanca in the period 2013–2016. The subregion-wide analysis of glacier elevation changes as a function of aspect does not reveal any clear trend (Fig. S28); however, when analysing only the Cordillera Blanca an east-to-west gradient is obvious (Fig. S29). We attribute this to the changed precipitation pattern during El Niño. Higher precipitation rates are typically found on east-facing slopes at the Cordillera Blanca, fed by moist air from the Amazon basin (Garreaud et al., 2009). However, during El Niño the westward flow of moist air is hampered by stronger westerlies (Vuille, 2013), which increases the precipitation gradient across the Cordillera Blanca.

Long-term glaciological mass balance measurements are only available for two glaciers (WGMS, 2018). For the Artesonraju and Yanamarey glaciers, both located in subregion R1, average mass budgets of −911 and −1.491 km m⁻² a⁻¹ (2013–2016) are revealed from field measurements, whereas we observed $-\mathrm{962}\pm \mathrm{228}$ and $-\mathrm{3167}\pm \mathrm{2317}$ kg m⁻² a⁻¹, respectively. Both methods show higher loss rates for the Yanamarey Glacier compared to the Artesonraju Glacier. The deviation of ∼50 % between both approaches at Yanamarey Glacier (overlap within the error range) can be partly attributed to different observation intervals and basin delineations. Moreover, glaciological mass balance measurements are typically based on some stake measurements in the lower, often most accessible, ablation zone and only a few measurements in the accumulation region. At the Zongo Glacier, Bolivia, Soruco et al. (2009b) performed a comparison of different mass balance methods (glaciological, hydrological, and geodetic) and revealed a strong offset in the glaciological mass balance estimates. The authors suggest that the limited number of field measurement points is not representative of the whole glacier and that the interpolation between the measurement sites can bias glacier-wide specific mass balance information. On the other hand, to obtain glacier-wide mass balance estimates using the geodetic method, there is interpolation over data voids in the elevation change fields needed (in our case we used the hypsometric elevation change distribution). Thus, we attribute the offset between the results from both methods to uncertainties due to interpolation, spatial sampling and different observation periods.

The correlation of the annual glaciological measurements with the average ONI of the respective observation periods indicates a trend towards increased mass loss and higher ELA during El Niño conditions (Fig. S30). These tendencies fit with the observations by Silvero and Jaquet (2017) and Morizawa et al. (2013) (see Sect. 6.1) and support our revealed correlation between the increased glacier wastage in the period 2013–2016 and the strong El Niño event during this period.

Since the highest density of glaciological point measurements are collected at the lowest section at both glaciers and considering the observation of Soruco et al. (2009b) at Zongo Glacier, Bolivia, that the terminus region strongly dominates the mass balance variations in a tropical glacier, we did an analysis of the glaciological mass balance at regions below 4900 and 5000 m a.s.l. at Yanamarey and Artesonraju glaciers, respectively. A trend toward increased mass losses after 2011 is obvious (Fig. S31, unfortunately the data at Yanamarey Glacier for the period 2009–2011 were incomplete), fitting with the revealed more negative geodetic mass balances after 2013. Moreover, the most negative values are derived for the period 2015–2016, which also supports our suggestion that ENSO forced the increased glacier wastage after 2013. The comparison between the geodetic and the average glaciological measurements at the terminus regions in the period 2013–2016 revealed a good agreement of both methods at Yanamarey Glacier (average glaciological: −3152 kg m⁻² a⁻¹; geodetic: −3164 kg m⁻² a⁻¹), whereas at Artesonraju Glacier, the geodetic measurements indicate higher mass loss (average glaciological: −4206 kg m⁻² a⁻¹; geodetic: −2696 kg m⁻² a⁻¹). This difference can be partly attributed to the different glacier basin definitions and slightly different observation intervals, but also to limitations of the individual methods, as discussed above.

In subregion R1, Huh et al. (2017) derived volume losses ranging from −0.019 to −0.150 km³ at six glaciers in the Cordillera Blanca from various elevation datasets in the period 1962–2008. This ice depletion corresponds to surface lowering rates of −0.20 to −1.4 m a⁻¹. Mark and Seltzer (2005) reported lowering rates ranging from −0.14 to −0.59 m a⁻¹ for three glaciers at the Nevado Queshque (Cordillera Blanca, subregion R1) in the interval 1962–1999. The findings of both studies are in the range of our revealed elevation change rates for subregion R1 of $-\mathrm{0.208}\pm \mathrm{0.065}$ and $-\mathrm{1.067}\pm \mathrm{0.273}$ m a⁻¹ in the periods 2000–2016 and 2013–2016, respectively.

In subregion R2, Salzmann et al. (2013) estimated an ice volume loss in the Cordillera Vilcanota of ∼40 %–45 % between 1962 and 2006 based on ice thickness derived from glacier inventory parameters and thickness–volume scaling. The authors pointed out that nearly all ice loss occurred after 1985, leading to a lowering rate of about −0.39 m a⁻¹ for the period 1985–2006. This value is similar to our measured average surface lowering in subregion R3 of $-\mathrm{0.440}\pm \mathrm{0.069}$ m a⁻¹ in the interval 2000–2013.

At the Coropuna volcano in subregion R3, Racoviteanu et al. (2007) measured an average glacier surface lowering of −5 m (−0.1 m a⁻¹) based on a SRTM DEM and digitizing of a topographic map from 1955. Peduzzi et al. (2010) calculated a mean surface lowering of $-\mathrm{0.2}\pm \mathrm{0.3}$ m a⁻¹ for the period 1955–2000 and 1955–2002. Our revealed average surface lowering of $-\mathrm{0.44}\pm \mathrm{0.045}$ m a⁻¹ at the ice cap of the Coropuna volcano indicates an increased glacier wastage after 2000, which correlates with the increase in glacier retreat and SLA uplift observed by Veettil et al. (2016) since 2000.

The continent-wide surface mass balance simulation by Mernild et al. (2017) revealed an average mass change of $-\mathrm{1230}\pm \mathrm{690}$ kg m² a⁻¹ for the Andes north of 27^∘ S in the epoch 2000–2014, which is much higher than our countrywide average of $-\mathrm{169}\pm \mathrm{43}$ kg m² a⁻¹ for the period 2000–2013. This large offset can only be partly explained by the different study domains, time intervals, and glacier inventories applied and could be caused by limitations in the applied statistical downscaling of global circulation data.

On the countrywide scale, we observed ∼19 % less mass loss than Braun et al. (2019) (Region 02-04) in the period 2000–2012/13 (first density scenario). However, their analysis also includes the glacier areas in Bolivia, like the Cordillera Real, where significant glacier wastage is reported (e.g. Soruco et al., 2009a), leading to higher average mass loss rates. The change rates of subregions R1 and R2 in the period 2000–2012/13 and the results of Braun et al. (2019) in Regions 02 and 03 show good agreement, even though they are based on different glacier inventories. However, we computed more negative mass budgets in subregion R3 compared to Braun et al. (2019) (Region 04). This can be partly explained by the differences in the region delineation. Braun et al. (2019) included the few ice-covered volcanoes in south-west Bolivia. However, a more considerable reason for this offset is the fact that we found the highest number of misclassified ice-covered areas in RGI 6.0 in this subregion (Sect. 6.1). This explains the bias toward lower mass loss rates in the results of Braun et al. (2019), who measured surface elevation changes based on RGI 6.0. In a continent-wide analysis of geodetic mass budgets, like that of Braun et al. (2019), it is beyond the scope of the studies to map glacier outlines fitting to the whole elevation database as well. However, the revealed offset suggests the inaccuracy caused by imprecise glacier outlines and highlights the need for large-scale temporally consistent glacier outlines.

Gardner et al. (2013) carried out a comprehensive worldwide estimate of the glacier contribution to sea level rise. They computed a mass budget of $-\mathrm{1080}\pm \mathrm{360}$ kg m⁻² a⁻¹ (2003–2009) at the low latitudes (RGI region definition) by means of the extrapolation of glaciological measurements. A nearly similar value of $-\mathrm{1.03}\pm \mathrm{0.83}$m w.e. a⁻¹ ($-\mathrm{1030}\pm \mathrm{830}$ kg m² a⁻¹) is revealed by Zemp et al. (2019) by upscaling of glaciological and geodetic mass balances for the period 2006–2016. Both studies reported ice loss rates that are about 3 times higher than our countrywide average of $-\mathrm{296}\pm \mathrm{41}$ kg m⁻² a⁻¹ (2000–2016). However, they are based on very few measurements (e.g. < 2 % spatial coverage for Zemp et al., 2019) in this region, resulting in the large uncertainties. Moreover, the representativeness of the settings of the sampled glaciers, which is not assessed, can also strongly influence region-wide estimates as discussed in Sect. 6.1. A comparison with GRACE measurements by Jacob et al. (2012) is also difficult since their spatial domain covers the whole of South America, only excluding Patagonia. Thus, the mass balance of the Peruvian glaciers cannot be disentangled from their results. These factors depict the limitations of using the present global mass balance estimates at the country or even mountain range level.

The above-mentioned global or countrywide estimates do not cover multiple periods and thus do not provide any information on temporal variability of the mass balances. Our analysis reveals strong temporal variations in the glacier changes that correlate with changing climate conditions and specific climatic events. These findings underline that mono-temporal analysis, especially when using short time intervals, can be biased by short-term climate anomalies like El Niño and highlights the need for further monitoring of the proceeding glacier recession in the tropical Andes.