Development of a spatial projection map of glacial retreat based on vulnerability maps in the Central Cordillera, Peru

Tropical glaciers are particularly sensitive to climate change. Therefore, their behavior and response to climate variability is differentiated (Hastenrath 1994; IPCC 2018). Unlike in higher latitudes, the processes of accumulation (increase of glacier mass) and ablation (loss of glacier mass) occur simultaneously; the highest precipitation occurs in the months of greatest radiative contribution, so that a snow cover that blocks melting does not form and the glacier is more prone to melting processes (Vuille et al. 2018).

At the beginning of the last decade, some studies agreed that the average temperature had increased by about 0.7 °C between 1939 and 2006 in the Tropical Andes (Bradley et al. 2009; Buytaert et al. 2017). In addition, between 2000 and 2016, 29% of the Tropical Peruvian Glaciers (which represent 71% of the world's total tropical glaciers; Kaser & Osmaston 2002) were lost (Buytaert & De Bièvre 2012). This occurred in a global context where the average temperature over the continents in the period 2011–2020 was 1.59 °C higher than in the period 1850–1900 (pre-industrial levels) and 0.99 °C (for the period 2001–2020) higher than the equivalent global average temperature change (IPCC 2018). Under these scenarios, and also since the behavior of glaciers in the context of climate change is variable, it is convenient to presume a relationship between the increase in average global temperature and an accelerated retreat of South American tropical glaciers (Kaser & Osmaston 2002; Herzog et al. 2011; Schauwecker et al. 2017; INAIGEM 2018).

Glacial recession has significant impacts on the ecosystem and downstream Andean communities. While glacial retreat directly results in increased meltwater discharge and thus increased water supply to streams, it also results in the formation of new ponds and flood zones that pose a danger to downstream communities. In addition, the water supply from glacial runoff is only positive up to a certain point; beyond the critical threshold, glacial input is of little significance to the streams. This represents a major problem with respect to water availability for the development of economic activities such as agricultural practices, livestock, hydroelectric power, and mining, among others. Given the importance of glaciers for maintaining the balance of mountain ecosystems and managing water resources across entire basins, research on glacier dynamics, particularly those associated with populated areas, is crucial (Rojas et al. 2021).

Traditionally, the studies on glacier surface changes focus on the evaluation of hydrological and meteorological variables and/or satellite images (Salzmann et al. 2013; Baraer et al. 2015; Veettil et al. 2016; Seehaus et al. 2019). These studies enable the quantification of glacier loss for the evaluated period. Furthermore, under a climate change perspective, they correlate the trend of glacier loss with the influence of climatic factors.

Zhang et al. (2016) developed a temperature index-based model to project glacier mass loss in the Altai Mountains under different climate scenarios (RCP4.5 and RCP8.5), representing a type of research with a different focus. They found that under a high emission scenario, glacier loss is more pronounced with an estimated disappearance of 53% of glaciers by 2100 compared to 2005. However, Hongkai et al. (2021) also developed a model to project changes in glaciers with future climate scenarios. However, their model integrated temperature and precipitation data from eight climate models of the CMIP6 project under two emission scenarios (RCP2.5 and RCP8.5) to evaluate changes in glacier ice thickness, area, volume, and flow for the period from 2030 to 2100. In contrast to the work by Zhang et al. (2016), this study not only quantified glacier area loss through 2100 but also provided a detailed visual representation of glacial retreat over time.

Yalcin (2020), on the other hand, approaches the study of glacial retreat from a perspective of their topography and geometry. A bivariate statistical analysis is used to determine the areas of a glacier that are most susceptible to melting based on past information and established criteria. The methodology is an adaptation of Vlcko et al. (1980), which mixes concepts related to Shannon entropy and frequency index (FR), originally used to evaluate slope stability on a regional scale (i.e., the degree of vulnerability of a terrain zone to landslides). The use of these methods to generate landslide vulnerability maps (VMs) is common (Bednarik et al. 2009; Constantin et al. 2011). It is important to remember that the concept of vulnerability can be defined as the degree of susceptibility or inability of a system to respond to the adverse effects of climate change, including climate variability and extreme events. In other words, it is a function of the character, magnitude, and rate of climate variation to which a system is exposed, its sensitivity, and its adaptive capacity (IPCC 2018). In the present study, the system is the glacier focus zone (GFZ) of the Central Cordillera of the Peruvian Andes.

Building upon the methodology proposed by Yalcin (2020) and adapting it to glaciers located in the Central Cordillera of Peru, this study distinguishes itself by not only quantifying glacier area loss but also providing a detailed visual representation of glacial retreat over time. While most glacier research has predominantly concentrated on understanding temporal dynamics of glacier mass loss, documenting glacier retreat processes, and identifying factors contributing to its accelerated pace in recent years, the novel contribution of this study lies in its holistic approach to transcend historical observations and prioritize modeling and projection of future glacier conditions. By integrating statistical analysis and satellite image processing, this research offers a comprehensive understanding of the regions most affected by glacier disappearance, thus providing valuable perceptions for effective decision-making in climate adaptation and water resource management. This emphasis on both quantitative analysis and visual representation enhances the utility of the findings for policymakers, resource managers, and researchers striving to address the challenges posed by climate change in glacierized regions.

All glaciers that occupy the Central Cordillera experiment a process of glacial retreat that dates back decades. In 1975, through the reconstruction of glacier surfaces with Landsat satellite images, a glacier coverage of 117.20 km² was estimated for the Central Cordillera. However, this area has been considerably reduced, covering only 79.32 km² in the second half of the 1990s, 51.91 km² in 2014, and 42.44 km² in 2016 (INAIGEM 2018).

For the present study, the GFZs were subdivided into two: ‘northern GFZ’ and ‘southern GFZ’ (Figure 1). This subdivision was enforced to comply with the methodology as the distance between glacier groups within the Central Cordillera is about ∼5 km² from the northern end of the southern GFZ (glacier name) to the southern end of the northern GFZ (glacier name), which makes it difficult to design parametric maps, especially the orientation map. This map should take the highest elevation point and divide the glacier group into cardinal zones from there. If the two groups of glaciers were considered one, it would be impossible to create this map.

The elevation in both glacial zones is between 4,600 and 5,800 m above sea level. The average annual temperature varies between 1 and 6 °C and the precipitation between 200 and 1,000 mm (SENAMHI 2021a). The predominant climate in this zone is semi-frigid, rainy in summer and deficient in winter.

In this first phase, a study area and a suitable time frame were selected based on the problems described in Section 1, the response to seasonality in the Central Peruvian Andes, and a literature review with respect to a minimum 30-year period as a reference over which climate variations can be evaluated. Subsequently, freely available satellite images were obtained, considering season, cloud cover, and El Niño influence. The selected images were preprocessed and delineated to the study area.

The choice of the study area responds to the problems defined in Section 1. Also for practical reasons, this area was chosen for the current study because it is a mountain range that responds well to the annual seasonality marked in Peru. Specifically, it is an area that is normally rainy during the southern hemisphere summer and dry in winter. For this reason, it is more likely to find the glacial areas free of clouds, unlike the glaciers located on the eastern slope of the Andes. The choice of time period is based on the recommendation found in the literature regarding a minimum period of 30 years as a reference over which climate variations can be evaluated (WMO 2017). A study involving the climatic analysis of a region should consider a period of at least 20–30 years. For this reason, the most recent year for which Landsat satellite images are available is 2021, and the first year of study would be 1990.

Initially, 26 Landsat TM5 and OLI 8 satellite images with a 30-m spatial resolution were obtained from the earthexplorer.usgs.gov portal and analyzed for the identification and subsequent discrimination of glacier areas (see Table 1). Of the 26 images, those that met the following criteria were selected: (a) cloud cover is less than 15% and, if possible, no cloud cover overlaps the areas of interest; (b) belong to the dry season, between May and September. At this time of the year, there is less ablation, and the snowpack is at its annual minimum, so it is expected that precipitation would not affect glacier volume (Medina & Mejía 2014). In Peru, in terms of precipitation, there are two distinct periods: the wet and dry seasons. Accumulation occurs during the wet season. Ablation occurs throughout the year, although it is higher during the wet season due to humidity and temperature conditions (Peruvian summer). These ablation processes are lesser during the dry season, when water demand is lower and the contribution of glaciers in response to water supply is reduced (Rivera et al. 2017; INAIGEM 2018). Finally, (c) the dates of the images should correspond to periods of El Niño or La Niña with neutral or weak intensity according to (A) the Coastal El Niño Index (ICEN, for its acronym in Spanish) proposed by the Peruvian Geophysical Institute – IGP (Takahashi et al. 2014), and (B) neutral to moderate intensity in the Southern Oscillation Index (SOI) proposed by the National Oceanic and Atmospheric Administration – NOAA (Rasmusson & Carpenter 1982) (see Table 1 and Tables S1 and S2 in the Supplementary material). Peruvian glaciers, as discussed in previous sections, are extremely sensitive to changes in climate. For this reason, Peruvian glaciers may be affected by climatic events that are a result of anomalies in global temperature patterns and regional ocean circulation, such as El Niño or La Niña (Fyffe et al. 2021).

Table 1

Landsat TM5 and OLI 8 image detail

View Large

Finally, seven images met the criteria and were used for this study (Table 1), corresponding to the dates August 1990, August 1999, July 2003, July 2005, July 2011, August 2013, and August 2021. Although there are images from other years that also meet the above selection requirements, they have not been considered because selecting a larger number of images is equivalent to investing more effort in preprocessing and processing them, when the main objective of the study is not to analyze a detailed change in the glacier surface, but to have some references as accurate as possible about this change over the 31-year period.

The seven satellite images were first reprojected to UTM 18 S coordinates. The images were then corrected using the semiautomatic classification plugin available in the GIS software (i.e., conversion of digital radiance values to values in W m⁻² r⁻¹ μm⁻¹, as recommended by Hantson et al. (2011) and Chavez (1996), among others). With these radiance values, the atmospheric correction procedure can be performed. Finally, the band set of the corrected images was cut on a mask of the Central Cordillera to focus on the study area.

To determine the glacier coverage in each selected year, satellite images were processed by applying spectral indices.

The Normalized Difference Snow Index (NDSI) introduced by Dozier (1989) is used here to identify and to map glacier cover. This index is based on the theory that snow mainly reflects radiation in the visible spectrum and absorbs radiation in the infrared wavelengths. In this sense, it is possible to distinguish snow from vegetation, soil, and clouds using the following relationship. ⁠. The NDSI is the normalized difference between the green and near infrared spectral bands. For Landsat 4–7, ⁠; and for Landsat 8–9, ⁠.

A limitation of the NDSI method is the confusion between glacier cover with light blue water bodies and shaded areas within the glacier cover. For this reason, we decided to follow the methodology proposed by Cea et al. (2005) to apply a cloud and shadow mask to eliminate the confusion between snow cover and cloud cover. In addition, to eliminate water bodies in the study area, the near infrared channel corresponding to band 4 in Landsat 5 was used. Clear water has a higher reflectivity in the blue band and decreases toward the near and mid-infrared, so water is excluded for values greater than 0.11 in band 4 (Salcedo 2011).

Although the NDSI distinguishes glacier cover from other types, it does not help differentiate between snow and ice. For this purpose, there is a methodology proposed by Keshri et al. (2009) and further explored by Herrera-Ossandón (2016) and Gutiérrez (2021). The process consists of applying various spectral indices to differentiate ‘supraglacial’ canopies. This methodology was used in this study, but it was concluded that the proportions of snow and ice found through the years are similar. Also, the percentage of ice mixed with debris agrees with the description made by INAIGEM for the Central Cordillera (INAIGEM 2018). This suggests that it may not be as necessary to discriminate snow from ice when using images from the dry season and when the analysis shows that the variation in the snow–ice ratio percentage over time is not significant. In addition, while using this methodology is an alternative to discriminating snow from ice, it is also somewhat arbitrary in relation to the threshold that defines the boundary between snow and ice for NDSII-2.

Two considerations were considered in the final generation of the glacier coverage maps. First, the minimum mappable area had to be greater than 0.005 km² as recommended by INAIGEM (2018); second, glacier areas ∼10 km away from the core of each GFZ were not included. The final areas found for the different periods are shown in Table 3.

As mentioned by Bueno et al. (2023), manual delineation of glaciers is a laborious and time-consuming activity. Therefore, due to the large amount of glacier area analyzed, the techniques mentioned in Section 2.3.1 as NDSI (Dozier 1989) apply cloud and shadow masks (Cea et al. 2005); near infrared channel to eliminate water bodies (Salcedo 2011) and histograms to approximate thresholds (Keshri et al. 2009) are used.

To validate the NDSI thresholds it was necessary to manually contrast with the 5-3-2 (for 1990, 1999, 2003, 2005, and 2011 Landsat TM5 images) or 6-4-3 (for 2013 and 2021 Landsat OLI 8 images) band set of the original satellite image to distinguish the glacier surface from other similar shapes, such as shiny ground, water, and clouds. In addition, for most images (1990, 1999, 2003, 2005, 2011, and 2013), the found glacier area was compared to high-quality imagery from Google Earth satellite imagery.

Finally, the thresholds of this study (0.31–0.51) were compared with other thresholds suggested for the Central Cordillera, Huaytapallana Cordillera, Vilcanota Cordillera, Vilcabamba Cordillera, Ampato Cordillera, and Blanca Cordillera. In this sense, INAIGEM (2018) uses thresholds of 0.4, as well as Soto et al. (2021), Calizaya et al. (2021), López-Moreno et al. (2014), and Zubieta & Lagos (2010); Muñoz et al. (2021) uses thresholds around 0.55 and 0.65; and Veettil et al. (2016) around 0.45 and 0.55. Thus, it has been verified that there is concordance between the thresholds used and those recommended. The differences between this study and the references are explained by the necessary adjustments made to the thresholds for the region and the specific date of image.

To validate the glacier coverage, the First National Glacier Inventory – The Glacier Ranges of Peru (INAIGEM 2018) was used as a reference. It was sought that there is some concordance and similarity between our results and those presented in this document. For example, the inventory mentions that for 1997, an average glacier area of 79.32 km² was found, while this study shows a glacier area of 77.52 and 80.72 km² for the years 1990 and 1999, respectively. The inventory also indicates that in 2014 the central mountain range had 51.91 km² of glacier area, while this study found 52.57 km² of glacier area in 2013. The values are equivalent, and the small differences can be explained by the specific NDSI threshold used and by the climatic peculiarity of the year in question (for example, there was a weak coastal El Niño in 1994–1995, a moderate La Niña in 1996, an extraordinary coastal El Niño in 1997–1998, and a weak El Niño in 2014).

It is important to note that the main objective of this study is not to introduce/implement/develop a glacier identification methodology but rather to perform commonly used techniques to generate inputs for the development of VMs and melt projections. For this reason, a more technical or in-depth validation of the NDSI thresholds or glacier coverage (statistical validation, field verification, etc.) was not performed.

The glacier coverage maps were used to generate VMs for the two GFZs. During this phase, the VM was also validated by comparing the actual estimated glacier retreat with the glacier coverage maps generated in the second phase. The purpose of this validation was to verify the accuracy of the vulnerability model.

The VM is a representation, based on certain criteria, of the areas within the GFZ that are most susceptible to glacier retreat. The criteria are those abiotic components of the environmental system that, individually and collectively, can provide clues to surface changes in glacier dynamics. They are defined by understanding climate as a function of elevation, aspect, topography, orientation, and other topographic parameters of glaciers. In fact, it is these parameters that drive changes in glacier dynamics (Davies et al. 2011; Li & Li 2014; Yang et al. 2015).

For the determination of the criteria used in the present work, the recommendations of Yalcin (2020), Davies et al. (2011), Li & Li (2014), and Yalcin (2019) were followed, in addition to a map of precipitation and flow direction, according to the definition of the criteria revealed in the previous paragraph.

Since landslides or glacial melting are complex systems of material and energy exchange with the environment, these events can be measured and described by the information entropy method (Yang & Qiao 2009). It is also widely used as a complement to the FR model (Demir et al. 2014; Pirasteh & Li 2017; Thapa & Bhandari 2019; Park et al. 2012), where p(x_i) is the probabilistic density of landslide occurrence for feature i, i.e., p(x_i) = FR/FR_total for a given criterion or causal factor. Landslide entropy can be understood as a measure of the influence of causal factors (criteria) on landslide occurrence or, in any case, it refers to the extent to which different factors influence the development of a landslide.

Two versions of the VM have been developed. The first one considers elevation, orientation, aspect, slope, flow direction, precipitation, and atmospheric temperature criteria; and the second one considers the same maps except the flow direction criteria and replaces the atmospheric temperature map by the surface temperature. This was done to achieve the highest possible percentage of validation. For representation purposes, the VM was classified into five classes using the natural breaks method (Jenks) (Yalcin 2020).

The validation process of the VM with respect to the actual glacier retreat estimated from the glacier cover maps was performed by simulating glacier melt areas in the VM corresponding to the periods 1990–1999, 1990–2003, 1990–2005, 1990–2011, 1990–2013, and 1990–2021. The glacier melt areas were estimated as the difference between the most recent glacier cover and the 1990 glacier cover. That is, the average melt (i.e., difference between the total areas of each GFZ) is not used, as this value considers the new areas of the most recent image. The areas used for validation are described in Table 4.

In general, the procedure consists of obtaining, for each period, a simulated melting map based on the VM and comparing this simulated melting with the one obtained by the difference in glacier coverage found with NDSI of the last year of the period (1999, 2003, 2005, 2011, 2013, 2021, respectively) with respect to the first year of the period (1990 = 24.55 km² for the southern GFZ and 56.18 km² for the northern GFZ).

Once the map was validated, areas of future melting were simulated according to the glacial retreat rate indicated by INAIGEM (2018) for the Central Cordillera equivalent to −1.39 km²/year. The first step was to cut the final VM according to the edge of glacier cover in 2021. From there, a melt projection map was created for the periods 2021–2025, 2025–2030, 2030–2035, 2035–2040, 2040–2045, 2045–2050, and 2050–2055.

It is worth mentioning that the retreat rate indicated by INAIGEM applies to the entire glacier area of the Central Cordillera. In this sense, to estimate the retreat of each GFZ, two things are assumed: first, that the sum of the southern and northern GFZs represents the total glacier area of the Central Cordillera; second, that the proportion of each GFZ area is maintained over time. The latter is supported by Table 3, where the percentages of each GFZ with respect to the total have a standard deviation of 1.82% for the southern GFZ and 1.78% for the northern GFZ. Therefore, to estimate the future retreat of each zone, it is assumed that the proportion will continue to be ∼30% for the southern GFZ and ∼70% for the northern GFZ with respect to the total area of the Central Cordillera. Finally, it is assumed that to maintain the proportion of glacier cover, the annual rate of glacier retreat must also maintain the same proportion. Then, based on the annual retreat ratio of 1.39 km²/year calculated by INAIGEM for the entire Central Cordillera, a ratio of 0.417 km²/year is estimated for the southern GFZ (30%) and 0.973 km²/year for the northern GFZ (70%). These values allow the estimation of future glacier cover in the following section.

The Tables 5 and Supplementary Tables S3 and S4 summarize the process for calculating the FR value and the weights for each criterion. In addition, the surface temperature weighting calculated from the 1990 satellite image is shown in Tables S5 and S6 in the Supplementary material.

The results corresponding to the areas of overlap between the glacier melt simulated from the VMs and the glacier melt obtained from the analysis of satellite imagery are shown in Tables 6 and 7, for the first and second versions of the VM, respectively. Note that the total glacier area for the first year (August 1990) is 23.83 km² for the southern GFZ and 53.69 km² for the northern GFZ.

Table 8 shows the results corresponding to the areas of glacier cover estimated for the years 2021–2055, the year in which it is estimated that all glacier cover in the Cordillera could disappear (see Table 8).

Four glacier VMs were generated, including two versions for the southern GFZ and two versions for the northern GFZ, based on (1) satellite imagery, (2) the selection of eight topographic and meteorological criteria, and (3) the use of statistical models such as the FR and Shannon's entropy method. Both versions of the maps agree on the trend of increasing validity over increasingly longer time periods (see overlap columns in Tables 6 and 7). These results are mainly due to the methodology used to calculate the VM and to the definition of the FR in this study (see Section 2). By using the year 2021 as a reference for the estimation of the glacier retreat area ratio, it is logical that the highest percentage of validation occurs for the period 1999–2021. The calculation method may also partly explain the relatively low validity values for the first two periods (1990–1999 and 1990–2003), since the glacier melt for these periods was different from that for 2021, both quantitatively and spatially.

When analyzing the overlap values (see Tables 6 and 7), it is important to consider the nonuniform distribution of melting values for each zone during different periods. Although glacier cover generally tends to decrease over time, the melting values for each zone do not necessarily follow this trend. In the case of southern GFZ, this assumption would only be fulfilled by removing the year 2005, which had the greatest melting for the study period. However, it is necessary to include a year between 2003 and 2011 to avoid a significant gap between the periods. Only images from the years 2004 and 2005 were available because images from other years did not meet the requirements mentioned in Section 2.2.2. These requirements were related to the El Niño phenomenon from August 2006 to January 2007, from January to September 2008, from May to September 2009, and La Niña phenomenon from April to December 2007 and from August to November 2010. The low glacier coverage during the selected period may be attributed to the low precipitation values recorded at surrounding stations during the March–April–May quarter of 2005. The precipitation range of 72–115 mm reached was lower than those recorded in the years before 2004 and after 2006. The overlapping results for this case, however, confirm the hypothesis that the higher the melting value, the greater the accuracy of the VM, regardless of the year.

These differences, and in general the results of the validation, are mainly due to the characteristics of the parametric maps. The choice of morphometric criteria (elevation, aspect, slope, flow direction, and orientation) is based on the premise that these criteria reflect the climatic and hydrological conditions, as well as the actual landscape composition of the study area. This choice is also supported by various studies that implement the FR model and Shannon entropy model to evaluate landslides (Chung & Fabbri 1999; Sharma et al. 2010, 2014; Pourghasemi et al. 2012; Demir et al. 2014; Jaafari et al. 2014; Lee & Pradhan 2016; Pirasteh & Li 2017; Thapa & Bhandari 2019; Panchal & Shrivastava 2021; Park et al. 2012). According to the validation of the VMs (see overlap columns in Tables 6 and 7), it has been shown that the decision to include or suppress certain parametric maps in the calculation of the final VM directly affects its accuracy. In fact, in version 2 of the VM, the parametric flow direction map is suppressed and the SENAMHI annual mean temperature map is replaced by a surface temperature map based on the infrared band of the 1990 satellite image. In this case, the validity increases for both the southern and northern GFZ.

One of the purposes of flow direction maps is to suggest areas where surface runoff may flow and to determine the spatial contribution to water accumulation zones. It can be assumed that in the case of glacial melting, a map of this type could be closely related to the location of new ponds, new channels, or new infiltration zones due to melting ice and snow. However, the weight assigned to it according to the statistical model used is the lowest of all (0.0032 and 0.0027 for southern and northern GFZ respectively). Since the classes are fixed by the nature of the criterion (i.e., the criterion is divided into eight classes according to the angle through which the water flows), the maximum Shannon entropy value H_jmax is also fixed and the final weighting responds mostly to the p_ij values. When calculating the p_ij values of the flow direction criterion, the distribution of the probabilities is quite uniform and there is not much variation between them (0.078, 0.119, 0.081, 0.103, 0.136, 0.168, 0.208, and 0.109, respectively, for each of the eight classes described in Figure 3(e)) compared to the probability distribution of the other criteria. This is confirmed by the Shannon entropy value H_j, which is known to be high since it is close to the H_jmax and suggests that classes in this case have little causal relationship with glacier melt, at least in terms of surface cover change. That is, melting between 1999 and 2021 is not characteristic of any of the flow direction map classes, or in other words the classes do not provide relevant information about glacier melting for the study area.

There are two main issues with the annual mean temperature map. First, this map does not represent the temperature of the glacier itself, and second, the original raster map has a resolution of 100 m; thus, to perform the corresponding calculations in this study, it had to be resampled to 30 m pixel size using the nearest neighbor method (Rukundo & Cao 2012). The new pixel values then contain a degree of error inherent in this interpolation method. This bias, added to the original error of the map according to the methodology used by SENAMHI (2021b) means to be working with inaccurate values of atmospheric temperature in the study area. High altitudinal temperatures and their interpolations require local corrections (Rau et al. 2013; Rau et al. 2018). Therefore, it is reasonable that the VM loses accuracy.

It is not possible to compare the results obtained here with other reference results, since this type of methodology has not been applied before in the study of Peruvian tropical glaciers. For this reason, the present study represents a novelty in the evaluation of glacier zones at risk. In any case, it is recommended to perform a prior validation of the methodology used in the present study, in particular the choice of number and the distribution of classes, such as that of Panchal & Shrivastava (2021). They examine the validity of the FR ratios and the weighting of each criterion by calculating the area under the ‘receiver operating curve’. This gives an idea of how accurate the obtained FR values and weights are or, in other words, how well the criteria and classes of each criterion have been chosen. Also, it would be interesting if an intermediate VM were produced for each period and then calculate a final VM as an average of all intermediate maps. In this way, the final map does not depend only on the conditions of the last year of the study but combines the characteristics of each melting period. The increased validity with respect to increasingly longer time periods and with respect to increasingly larger melt areas also suggests that the map for future melt areas (periods greater than 2021) would be consistently accurate, or at least more accurate than 2021.

Finally, in this study, satellite image processing and glacier cover identification represents intermediate results that serve as input for the preparation of parametric maps, the preparation of VMs, and their validation. These results show an average retreat of 0.92 km²/year for the total area, i.e., the sum of the areas of the southern GFZ and the areas of the northern GFZ (see Table 3). This retreat ratio differs from the ratio of 1.39 km²/year proposed by INAIGEM for several reasons not discussed in this study, but related to the number of images analyzed, periodicity, mappable areas, and thresholds for NDSI.

This study presents an analysis in the glaciated Central Cordillera of the Peruvian Andes that integrates remote sensing techniques based on spectral indices, specifically the NDSI, together with a statistical evaluation to project future glacier retreat from Landsat images. The statistical analysis developed was aimed at developing VMs and was based on glacier melt indices (FR) and on the numerical importance of certain factors that a priori influence glacier dynamics (weighting of criteria according to Shannon's entropy model), to quantify the vulnerability of the glaciers evaluated. Two versions of the VM were developed for each GFZ based on eight criteria (aspect, flow direction, altitude, orientation, slope, precipitation, atmospheric temperature, and surface temperature), one using all eight criteria but surface temperature and the other using all eight criteria but atmospheric temperature and flow direction, resulting in a total of four VMs. The weight of these criteria in the occurrence of the glacier retreat process was estimated. For both zones evaluated (north and south), it was found that the factors that most influence the dynamics of glaciers in the Central Cordillera are surface temperature (0.0272 for southern GFZ and 0.0426 for northern GFZ), atmospheric temperature (0.0271 for southern GFZ and 0.0253 for northern GFZ), elevation (0.0249 for southern GFZ and 0.0234 for northern GFZ), and precipitation (0.0236 for southern GFZ and 0.0212 for northern GFZ). Based on the vulnerability values, the pixels of the maps were classified into five groups (very low, low, medium, high, and very high). In general, for all maps generated, loose glacier zones and regions located at lower elevations are identified as the most vulnerable.

It was also found that the second VM is more accurate than the first, with an average overlap of 68.83% for the southern GFZ and 69.63% for the northern GFZ compared to the first (64.63% for the southern GFZ and 62.86% for the northern GFZ). Differences between the two versions of the VMs are explained by the choice of the parametric maps (especially those of temperature, flow direction, and precipitation), the choice of the number and distribution of classes for the other criteria, the methodology for calculating the VMs, the FR values, and the final weighting of each criterion. It is concluded that it is necessary to iterate on the classification of the criteria to find the maximum possible validity.

Finally, with the validated VM and based on the annual glacier retreat rate (1.89 km²/year) in the Central Cordillera, the loss of glacier coverage was estimated for the years 2025, 2030, 2035, 2045, 2050, and 2055. It was observed that during the evaluated period, the surface area occupied by glaciers decreased from 33.94 and 14.32 km² to 0.86 and 0.15 km², in the North and South zones, respectively. Furthermore, considering the retreat rate, total disappearance is expected by the year 2056. In addition, regarding the regions that will disappear first, it was observed in the southern GFZ that the highest altitude region of the main glacier body would be the last to vanish. In contrast, in the northern GFZ, it was observed that for the last estimated glacier retreat period (2050–2055), scattered remnants of the different glacier bodies will remain, without being concentrated in a specific region. Moreover, it is noted that the eastern margin of both GFZs is diminishing its glacier coverage at a faster rate than the western margin.

This study represents the first attempt to assess the vulnerability of glaciers in the Central Cordillera through statistical analysis and GIS, utilizing the methodology developed by Yalcin (2020). The findings are anticipated to provide valuable insights for government agencies involved in water resource management and mountain ecosystem preservation. By leveraging the data obtained, these entities now have a starting point for investigating potential hazards to populations residing near the identified vulnerable regions within the Central Cordillera and for establishing preventive measures to manage them. Furthermore, it is hoped that this work will inspire the replication of the methodology in other tropical glaciers of the region in future studies.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.