Increase in average global temperature over the last few decades has caused an accelerated retreat of tropical glaciers. Andean populations live in strict dependence on the water services provided by mountains and glaciers. The present study aims to generate a glacier melt projection map in the Peruvian Central Cordillera based on vulnerability maps over the 1990–2021 period. Seven satellite images were selected to determine the changes in glacier coverage based on normalized indexes. Subsequently, seven parametric maps consisting of terrain and climate characteristics were assimilated into a vulnerability analysis based on the frequency index and the Shannon entropy index model, allowing one to identify areas most susceptible to glacial retreat. The results show that the most important criteria for the southern and northern glacial study areas are surface temperature, elevation, precipitation, aspect, orientation, and slope. The validation results revealed the most accurate set of parameters from the vulnerability map in terms of projecting melting areas and were used to produce a spatial projection map for the period 2021–2055. From 2021, a glacier loss in the range 84–98% would be reached by the 2050s.

  • A versatile method is proposed to study vulnerabilities in complex terrains and climatic uncertainty such as glacial context.

  • Updated glacier retreat maps over a long time period from 1990 to 2021 based on remote sensing in Central Andes have been used in this study.

  • A methodology to project glacier retreat in a context of scarce or limited information is used in this study.

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.

Study area

The study area is called the GFZ within Peru's Central Cordillera. This area is a group of different glaciers that extend throughout the departments of Lima and Junín, specifically in the provinces of Huarochirí, Yauyos, Jauja, and Yauli (see Figure 1). The study area is located between southern latitudes 12°20′41.8″ and 11°43′50.93″, and western longitudes 76°8′54.87″ and 75°51′49.63″.
Figure 1

Study area. (a) Map of Peru showing the departments of Lima and Junín. (b) Central Cordillera within the districts of Huarochirí, Yauli, Jauja, and Yauli. (c) Northern GFZ in red and Southern GFZ in blue within the Central Cordillera.

Figure 1

Study area. (a) Map of Peru showing the departments of Lima and Junín. (b) Central Cordillera within the districts of Huarochirí, Yauli, Jauja, and Yauli. (c) Northern GFZ in red and Southern GFZ in blue within the Central Cordillera.

Close modal

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 km2 was estimated for the Central Cordillera. However, this area has been considerably reduced, covering only 79.32 km2 in the second half of the 1990s, 51.91 km2 in 2014, and 42.44 km2 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 km2 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.

Collect and preprocess data

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.

Study area and time frame

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.

Landsat images

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

 
 

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.

Image preprocessing

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−2 r−1 μm−1, 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.

Determination of glacier cover

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

Glacier cover identification

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

To select an appropriate NDSI threshold, we followed the recommendation of Keshri et al. (2009), who suggested the use of histograms. For an NDSI image, the histogram usually reveals a bimodal distribution, and the threshold can be understood as the range of minimum values that separates one distribution from the other. To approximate a threshold, the NDSI image was cropped to a representative glacier in the study area so that negative values (debris, bare ground, shadows, gaps, etc.) would be minimized and positive values would have more weight in the histogram. Figure 2 shows the NDSI of the Acopallca, Cullec, and Ticlla glaciers in the southern GFZ (which formed all three a single glacier body in 1990) and the histogram. In this case, the threshold is around the value 0.4, which coincides with the thresholds recommended by Dozier (1989) and Hall et al. (1995). The same procedure was used to determine the NDSI threshold for the other years (detailed in Table 2). It is important to mention that it should not matter which glacier is chosen to perform this analysis; the NDSI range and threshold are almost invariant between one and another glacier of the same image because a threshold varies depending on the sensor, seasons, and dates, unless a bad cut is made to the representative glacier or the chosen glacier is covered by clouds or contains a high percentage of shadows (Dozier 1989; Hall et al. 1995).
Table 2

Thresholds for glacier cover classification

ID Landsat sceneDateNDSI threshold
LT50070681990216CUB00 4 August 1990 0.34 
LT50070681999225COA00 13 August 1999 0.31 
LT50070682003188CUB00 7 July 2003 0.41 
LT50070682005209CUB00 28 July 2005 0.49 
LT50070682011210CPE00 29 July 2011 0.36 
LT80070682013231LNG01 19 August 2013 0.51 
LT80070682021237LGN00 25 August 2021 0.41 
ID Landsat sceneDateNDSI threshold
LT50070681990216CUB00 4 August 1990 0.34 
LT50070681999225COA00 13 August 1999 0.31 
LT50070682003188CUB00 7 July 2003 0.41 
LT50070682005209CUB00 28 July 2005 0.49 
LT50070682011210CPE00 29 July 2011 0.36 
LT80070682013231LNG01 19 August 2013 0.51 
LT80070682021237LGN00 25 August 2021 0.41 
Figure 2

(a) Acopallca, Cullec, and Ticlla glaciers NDSI; (b) histogram of NDSI map.

Figure 2

(a) Acopallca, Cullec, and Ticlla glaciers NDSI; (b) histogram of NDSI map.

Close modal

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.

Quantification of glacier coverage

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

Table 3

Thresholds and glacier areas for the period 1990–2021

ID Landsat sceneDateNDSI thresholdTotal area (km2)Southern GFZ
Northern GFZ
km2% of total areakm2% of total area
LT50070681990216CUB00 04 August 1990 0.34 77.52 23.83 30.74% 53.69 69.26% 
LT50070681999225COA00 13 August 1999 0.31 80.72 24.55 30.41% 56.18 69.60% 
LT50070682003188CUB00 07 July 2003 0.41 73.40 23.37 31.84% 50.04 68.17% 
LT50070682005209CUB00 28 July 2005 0.49 50.19 13.41 26.72% 36.68 73.08% 
LT50070682011210CPE00 29 July 2011 0.36 83.71 23.57 28.16% 60.13 71.83% 
LT80070682013231LNG01 19 August 2013 0.51 52.57 14.17 26.95% 38.40 73.05% 
LT80070682021237LGN00 25 August 2021 0.41 48.26 14.32 29.67% 33.94 70.33% 
ID Landsat sceneDateNDSI thresholdTotal area (km2)Southern GFZ
Northern GFZ
km2% of total areakm2% of total area
LT50070681990216CUB00 04 August 1990 0.34 77.52 23.83 30.74% 53.69 69.26% 
LT50070681999225COA00 13 August 1999 0.31 80.72 24.55 30.41% 56.18 69.60% 
LT50070682003188CUB00 07 July 2003 0.41 73.40 23.37 31.84% 50.04 68.17% 
LT50070682005209CUB00 28 July 2005 0.49 50.19 13.41 26.72% 36.68 73.08% 
LT50070682011210CPE00 29 July 2011 0.36 83.71 23.57 28.16% 60.13 71.83% 
LT80070682013231LNG01 19 August 2013 0.51 52.57 14.17 26.95% 38.40 73.05% 
LT80070682021237LGN00 25 August 2021 0.41 48.26 14.32 29.67% 33.94 70.33% 

Validation of NDSI thresholds and glacier coverage

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 km2 was found, while this study shows a glacier area of 77.52 and 80.72 km2 for the years 1990 and 1999, respectively. The inventory also indicates that in 2014 the central mountain range had 51.91 km2 of glacier area, while this study found 52.57 km2 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.

Glacier VMs

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.

Criteria definition

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.

Parametric maps

A total of eight parametric maps were created for each GFZ based on (1) elevation, (2) orientation, (3) aspect, (4) slope, (5) flow direction, (6) precipitation, (7) atmospheric temperature, and (8) surface temperature.
  • – For maps of elevation, aspect, slope, and flow direction, the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER GDEM v2) was used.

  • – For the orientation map, the Lines around Points plugin and the intersection and dissolution functions were used to obtain the cardinal orientation maps (north, northeast, east, southeast, south, southwest, west, northwest) based on the highest point of each GFZ (see Figure 3). The division of the glaciated area of the Central Cordillera into two zones is in part a response to the development of this map, since taking the highest point is a difficulty when considering the entire Central Cordillera as a single glacial zone.

  • – For the precipitation map, the 100 m resolution raster product of the Annual Precipitation Map (1981–2010) (SENAMHI 2021b) was used.

  • – For the atmospheric temperature, the 100 m resolution raster products of the Annual Minimum Temperature Map (1981–2010) and the Annual Maximum Temperature Map (1981–2010) (SENAMHI 2021b) were used.

  • – The surface temperature map is an alternative to the average temperature map obtained from SENAMHI. The construction of this map consists of using the pixel values of the shortwave infrared band (band 6 in Landsat 5 TM) and converting them to absolute radiance values following the methodology from Giannini et al. (2015).

Figure 3

Parametric maps. (a) Elevation maps, (b) orientation maps, (c) aspect maps, (d) slope maps, (e) flow direction maps, (f) annual precipitation maps, (g) annual temperature maps, and (h) surface temperature maps.

Figure 3

Parametric maps. (a) Elevation maps, (b) orientation maps, (c) aspect maps, (d) slope maps, (e) flow direction maps, (f) annual precipitation maps, (g) annual temperature maps, and (h) surface temperature maps.

Close modal

Vulnerability maps

The physical sense of entropy refers to the dynamic thermal disorder of the microstructure of a system; disorder that gives rise to temperature, heat capacity, and heat energy (Kostic 2014). Shannon (1948), in his Mathematical Theory of Communication, used this idea to refer to entropy as a measure of the uncertainty associated with a probabilistic distribution within any system. In other words, Shannon's entropy takes into account the inputs of a system to indicate the amount of information required to obtain (guess or correctly guess) the output of the system. The average uncertainty, denoted by the function H(X) for an event or probabilistic distribution, is given by the following expression:
formula
where H(X) is the Shannon entropy measured in bits; p is the probability of occurrences of xi; log2 is the bases chosen to obtain entropy in units of (bits).

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(xi) is the probabilistic density of landslide occurrence for feature i, i.e., p(xi) = FR/FRtotal 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.

In this study, the VM is generated from the values of the FR of each class belonging to each criterion and the resulting weight of each criterion in the form:
formula
where VM is each pixel value of the VM, Wj is the resulting weight calculated from the Shannon entropy method and the maps are those reclassified (RECL) according to the FR assigned to each class (the set of pixels belonging to a class gets that class's FR value by means of a raster reclassification).
The procedure for calculating the FR and the weight is to calculate the areas based on the following equations:
formula
where b is the melting area ratio; a is the glacial area ratio; i (class) = 1, 2, 3, …, n.; and j (criteria) = 1, 2, 3, …, m.
formula
formula
formula
where Sj is the total number of classes of the criterion j.
formula
formula
where
  • FR is the ratio of the area of glacial melt to the total study area;

  • pij represents the relative probability, which compares how much each class contributes to the total melting of the criterion;

  • pij (mean) is the average of all the relative probabilities pij;

  • Hj is the Shannon entropy value in units of ‘bits’. This value represents how evenly the probabilities pij are distributed across all classes for a given criterion j. A high entropy means that the probabilities are more dispersed or evenly distributed across all classes i of criterion j. This would indicate that glacier melt as a function of criterion j is less predictable. A low entropy, however, means that the probabilities pij are more concentrated or less uniformly distributed in the classes i for criterion j. This indicates that glacier melt is more predictable as a function of criterion j;

  • Hjmax is the maximum Shannon entropy value in units of ‘bits’. This would be the state of maximum uncertainty. That is, the entropy value when melting is equally likely in all classes i of a factor j;

  • Ij is the information coefficient, which is the additional information provided by criterion j relative to maximum uncertainty Hjmax. For example, if Hj is equal to Hjmax, it means that melting is equally likely for all i classes of a criterion j, and Ij would be 0 since no additional information has been provided and the uncertainty is still maximum. If Hj is less than Hjmax, it would indicate that the distribution is less uniform and Ij would approach 1. Here, criterion j provides valuable information or reduces the uncertainty in the glacier melt projection in terms of class distribution;

  • Wj represents the relative importance of each criterion in melting the GFZ. This value is a way of weighting the informative importance of criterion j (Ij) with its average distribution across all classes i (pij (mean)). It can be thought of as a composite measure of how criterion j contributes to glacier melt vulnerability in terms of information and average distribution.

Note that the FR is fitted as follows:
formula
where
formula
and
formula

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

Validation of VMs

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.

Table 4

Areas used to validate the VMs

iPeriodSouthern GFZ
Northern GFZ
Last year of period area (km2)Melted area (km2) (Ai)Average melted area (km2)Last year of period area (km2)Melted area (km2) (Ai)Average melted area (km2)
90–99 24.55 4.16 −0.71 56.18 10.23 −2.50 
90–03 23.37 4.97 0.47 50.04 14.89 3.65 
90–05 13.41 10.74 10.42 36.68 23.20 16.92 
90–11 23.57 6.30 0.26 60.13 12.38 −6.44 
90–13 14.17 10.07 9.66 38.40 21.86 15.29 
90–21 14.32 10.24 9.51 33.94 24.37 19.75 
iPeriodSouthern GFZ
Northern GFZ
Last year of period area (km2)Melted area (km2) (Ai)Average melted area (km2)Last year of period area (km2)Melted area (km2) (Ai)Average melted area (km2)
90–99 24.55 4.16 −0.71 56.18 10.23 −2.50 
90–03 23.37 4.97 0.47 50.04 14.89 3.65 
90–05 13.41 10.74 10.42 36.68 23.20 16.92 
90–11 23.57 6.30 0.26 60.13 12.38 −6.44 
90–13 14.17 10.07 9.66 38.40 21.86 15.29 
90–21 14.32 10.24 9.51 33.94 24.37 19.75 

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 km2 for the southern GFZ and 56.18 km2 for the northern GFZ).

The first step consists of taking a matrix of pixels from the VM, where each pixel has an ID, an area value in km2 (0.0009 km2 per pixel, since these are 30 m resolution maps), and a vulnerability value. The next step is to sort the matrix in descending order of vulnerability and accumulate the area value until the melted area value Ai is reached. This is because the pixels with the highest vulnerability value should be melted first. The pixels that were selected are grouped into a separate map that represents the simulated melting for that period. This map is shown as orange areas in Figure 4. An overlap (green areas in Figure 4) analysis was then performed to determine the percentage of overlap area between the actual melting zones described in Table 4 and the melting zones according to the VMs. The results of this analysis are shown as green areas in Figure 4. This figure provides an example of the validation for the period 1990–2021.
Figure 4

Validation for the period 1990–2021. (a) Simulated melting areas of period 1990–2021, total area of 4.16 km2 calculated from the most vulnerable pixels of the VM v1. (b) Actual melting areas (1990–2021), total area of 4.16 km2 calculated as the difference between glacier cover of 1990 and glacier cover of 2021. (c) Intersection between simulated melting areas and actual melting areas; difference areas representing those areas of simulated melt that do not overlap with actual melt.

Figure 4

Validation for the period 1990–2021. (a) Simulated melting areas of period 1990–2021, total area of 4.16 km2 calculated from the most vulnerable pixels of the VM v1. (b) Actual melting areas (1990–2021), total area of 4.16 km2 calculated as the difference between glacier cover of 1990 and glacier cover of 2021. (c) Intersection between simulated melting areas and actual melting areas; difference areas representing those areas of simulated melt that do not overlap with actual melt.

Close modal

Estimation of future glacier retreat scenarios

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 km2/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 km2/year calculated by INAIGEM for the entire Central Cordillera, a ratio of 0.417 km2/year is estimated for the southern GFZ (30%) and 0.973 km2/year for the northern GFZ (70%). These values allow the estimation of future glacier cover in the following section.

Criteria weighting

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.

Table 5

Final criteria weighting

Southern GFZ
Northern GFZ
CriteriaWeightCriteriaWeight
Surface temperature 0.0272 Surface temperature 0.0426 
Atmospheric temperature (annual average) 0.0271 Atmospheric temperature (annual average) 0.0253 
Elevation 0.0249 Elevation 0.0234 
Precipitation 0.0236 Precipitation 0.0212 
Aspect 0.0206 Aspect 0.0070 
Orientation 0.0068 Orientation 0.0055 
Slope 0.0054 Slope 0.0027 
Flow direction 0.0032 Flow direction 0.0027 
Southern GFZ
Northern GFZ
CriteriaWeightCriteriaWeight
Surface temperature 0.0272 Surface temperature 0.0426 
Atmospheric temperature (annual average) 0.0271 Atmospheric temperature (annual average) 0.0253 
Elevation 0.0249 Elevation 0.0234 
Precipitation 0.0236 Precipitation 0.0212 
Aspect 0.0206 Aspect 0.0070 
Orientation 0.0068 Orientation 0.0055 
Slope 0.0054 Slope 0.0027 
Flow direction 0.0032 Flow direction 0.0027 

Vulnerability maps

The first versions of the hazard maps are shown in Figure 5(a). The second versions of the VMs (those that include surface temperature from the satellite and do not include the flow direction map, precipitation and average temperature) are shown in Figure 5(b).
Figure 5

VMs (a) Version 1, elaboration includes (1) elevation, (2) orientation, (3) aspect, (4) slope, (5) flow direction, (6) precipitation, and (7) average temperature. (b) Version 2, elaboration includes (1) elevation, (2) orientation, (3) aspect, (4) slope, and (5) surface temperature. The blue pixels in the legend indicate areas of very low vulnerability (lower probability of melting in the near future), while the red ones indicate areas of very high vulnerability (higher probability of melting in the near future). Light blue (low vulnerability), yellow (medium vulnerability), and orange (high vulnerability) are intermediate vulnerability zones.

Figure 5

VMs (a) Version 1, elaboration includes (1) elevation, (2) orientation, (3) aspect, (4) slope, (5) flow direction, (6) precipitation, and (7) average temperature. (b) Version 2, elaboration includes (1) elevation, (2) orientation, (3) aspect, (4) slope, and (5) surface temperature. The blue pixels in the legend indicate areas of very low vulnerability (lower probability of melting in the near future), while the red ones indicate areas of very high vulnerability (higher probability of melting in the near future). Light blue (low vulnerability), yellow (medium vulnerability), and orange (high vulnerability) are intermediate vulnerability zones.

Close modal

VMs validation

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 km2 for the southern GFZ and 53.69 km2 for the northern GFZ.

Table 6

Validation of the first version of the VM

VM 1st version
PeriodSouthern GFZ
Northern GFZ
Melted area (km2)Overlap
Melted area (km2)Overlap
Area (km2)%Area (km2)%
90–99 4.2 2.4 58.3 10.2 5.7 55.2 
90–03 5.0 2.4 48.3 14.9 8.8 59.3 
90–05 10.7 7.3 68.1 23.2 14.6 63.0 
90–11 6.3 4.3 68.7 12.4 7.6 61.4 
90–13 10.1 7.1 70.8 25.7 14.7 67.3 
90–21 10.2 7.5 73.6 24.4 17.3 71.0 
VM 1st version
PeriodSouthern GFZ
Northern GFZ
Melted area (km2)Overlap
Melted area (km2)Overlap
Area (km2)%Area (km2)%
90–99 4.2 2.4 58.3 10.2 5.7 55.2 
90–03 5.0 2.4 48.3 14.9 8.8 59.3 
90–05 10.7 7.3 68.1 23.2 14.6 63.0 
90–11 6.3 4.3 68.7 12.4 7.6 61.4 
90–13 10.1 7.1 70.8 25.7 14.7 67.3 
90–21 10.2 7.5 73.6 24.4 17.3 71.0 
Table 7

Validation of the second version of the VM

VM 2nd version
PeriodSouthern GFZ
Northern GFZ
Melted area (km2)Overlap
Melted area (km2)Overlap
Area (km2)%Area (km2)%
90–99 4.2 2.5 58.8 10.2 6.3 61.3 
90–03 5.0 2.8 56.3 14.9 10.2 68.4 
90–05 10.7 7.8 72.2 23.2 17.0 73.3 
90–11 6.3 4.5 71.8 12.4 7.7 62.4 
90–13 10.1 7.6 75.6 25.7 16.5 75.6 
90–21 10.2 8.0 78.3 24.4 18.7 76.8 
VM 2nd version
PeriodSouthern GFZ
Northern GFZ
Melted area (km2)Overlap
Melted area (km2)Overlap
Area (km2)%Area (km2)%
90–99 4.2 2.5 58.8 10.2 6.3 61.3 
90–03 5.0 2.8 56.3 14.9 10.2 68.4 
90–05 10.7 7.8 72.2 23.2 17.0 73.3 
90–11 6.3 4.5 71.8 12.4 7.7 62.4 
90–13 10.1 7.6 75.6 25.7 16.5 75.6 
90–21 10.2 8.0 78.3 24.4 18.7 76.8 

Future scenarios

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

Table 8

Estimation of the area of future glacier cover in the Central Cordillera

YearGlacier coverage of the central cordillera (km2)Northern GFZ (km2)Southern GFZ (km2)
2021 48.26 33.94 14.32 
2025 42.71 30.05 12.66 
2030 35.75 25.18 10.57 
2035 28.81 20.32 8.49 
2040 21.85 15.45 6.40 
2045 14.91 10.59 4.32 
2050 7.95 5.72 2.23 
2055 1.01 0.86 0.15 
YearGlacier coverage of the central cordillera (km2)Northern GFZ (km2)Southern GFZ (km2)
2021 48.26 33.94 14.32 
2025 42.71 30.05 12.66 
2030 35.75 25.18 10.57 
2035 28.81 20.32 8.49 
2040 21.85 15.45 6.40 
2045 14.91 10.59 4.32 
2050 7.95 5.72 2.23 
2055 1.01 0.86 0.15 

Northern GFZ and Southern GFZ.

Based on the estimated glacier surface loss and the generated VMs, it is also possible to generate a melt projection map for the periods 2021–2025, 2025–2030, 2030–2035, 2035–2040, 2040–2045, 2045–2050, and 2050–2055 (see Figures 6 and 7). From 2021, a glacier loss in the range of 84–98% would be reached by the 2050s over the entire Central Cordillera.
Figure 6

Melt projection map for northern GFZ.

Figure 6

Melt projection map for northern GFZ.

Close modal
Figure 7

Melt projection map for southern GFZ.

Figure 7

Melt projection map for southern GFZ.

Close modal

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.

The validation results also show that the second versions of the VMs are more accurate than the first. In fact, the second version shows an increase in the percentage of overlap for all periods in both GFZs. The maximum values of overlap obtained in this study correspond to 78.3 and 76.8% for the southern and northern GFZ, respectively. The difference between the two versions lies mainly in the fact that the first version overestimates the melting areas in loose zones (as in the northern part of the southern GFZ) and in zones of intermediate elevation (5,000–5,300 m). In addition, there are areas that the first version of the VM simply does not consider to be highly vulnerable, while the second version does. For example, Figure 8 shows the actual melt area for the period 1999–2021 and the melt areas simulated with each version of the VMs. There are some clear differences, such as the fact that the first version does not consider the northern and eastern zones of Cullec (framed in the red ellipses from Figure 8) as melt areas for this period.
Figure 8

Comparison of real and simulated melted areas in the Southern GFZ.

Figure 8

Comparison of real and simulated melted areas in the Southern GFZ.

Close modal

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 Hjmax is also fixed and the final weighting responds mostly to the pij values. When calculating the pij 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 Hj, which is known to be high since it is close to the Hjmax 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 km2/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 km2/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 km2/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 km2 to 0.86 and 0.15 km2, 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.

Baraer
M.
,
McKenzie
J.
,
Mark
B.
,
Gordon
R.
,
Bury
J.
,
Condom
T.
,
Gomez
J.
,
Knox
S.
&
Fortner
S.
2015
Contribution of groundwater to the outflow from ungauged glacierized catchments: A multi-site study in the Tropical Cordillera Blanca, Peru
.
Hydrological Processes
29
,
2561
2581
.
Bednarik
M.
,
Magulová
B.
,
Matys
M.
&
Marschalko
M.
2009
Landslide susceptibility assessment of the Kraľovany–Liptovský Mikuláš railway case study
.
Physics and Chemistry of the Earth
3–5
,
162
171
.
Buytaert
W.
&
De Bièvre
B.
2012
Water for cities: The impact of climate change and demographic growth in the tropical Andes
.
Water Resources Research
48
(
8
),
8503
2012
.
Buytaert
W.
,
Moulds
S.
,
Acosta
L.
,
De Bièvre
B.
,
Olmos
C.
,
Villacis
M.
,
Tovar
C.
&
Verbist
K.
2017
Glacial melt content of water use in the tropical Andes
.
Environmental Research Letters
12
,
114014
.
Calizaya
E.
,
Mejía
A.
,
Barboza
E.
,
Calizaya
F.
,
Corroto
F.
,
Salas
R.
,
Vásquez
H.
&
Turpo
E.
2021
Modelling snowmelt runoff from tropical Andean glaciers under climate change scenarios in the Santa River sub-basin (Peru)
.
Water
13
,
3535
.
Cea
C.
,
Cristóbal
J.
,
Serra
P.
&
Pons
X.
2005
Mejoras en la detección semiautomática de nubes y sombras en imágenes Landsat
. In
XI Congreso Nacional de Teledetección
,
Puerto de la Cruz, Tenerife
.
Chavez
P.
1996
Image-based atmospheric corrections-revisited and improved
.
Photogrammetric Engineering and Remote Sensing
62
(
9
),
1025
1035
.
Chung
C.
&
Fabbri
A.
1999
Probabilistic prediction models for landslide hazard mapping
.
Photogrammetric Engineering and Remote Sensing
65
(
12
),
1388
1399
.
Constantin
M.
,
Bednarik
M.
,
Jurchescu
M.
&
Vlaicu
M.
2011
Landslide susceptibility assessment using the bivariate statistical analysis and the index of entropy in the Sibiciu Basin (Romania)
.
Environmental Earth Sciences
2
(
63
),
397
406
.
Davies
B.
,
Carrivick
J.
,
Glasser
N.
,
Hambrey
M.
&
Smellie
J.
2011
A new glacier inventory for 2009 reveals spatial and temporal variability in glacier response to atmospheric warming in the Northern Antarctic Peninsula, 1988–2009
.
The Cryosphere Discussions
5
,
3541
3595
.
Fyffe
C. L.
,
Potter
E.
,
Fugger
S.
,
Orr
A.
,
Fatichi
S.
,
Loarte
E.
,
Medina
K.
,
Hellström
R.
,
Bernat
M.
,
Aubry-Wake
C.
,
Gurgiser
W.
,
Perry
B.
,
Suarez
W.
,
Quincey
D.
&
Pellicciotti
F.
2021
The energy and mass balance of Peruvian glaciers
.
Journal of Geophysical Research: Atmospheres
126
,
e2021JD034911
.
Giannini
M.
,
Belfiore
O.
,
Parente
C.
&
Santamaria
R.
2015
Land surface temperature from Landsat 5 TM images: Comparison of different methods using airborne thermal data
.
Journal of Engineering Science and Technology Review
8
(
3
),
83
90
.
Gutiérrez
F.
2021
Aplicación de índices Espectrales Para la Identificación de Cubiertas Glaciares en el Glaciar Universidad, Región del Libertador Bernardo O'Higgins
.
Thesis
, Universidad de Chile,
Santiago, Chile
.
Hall
D. K.
,
Riggs
G. A.
&
Salomonson
V. V.
1995
Development of methods for mapping global snow-cover using moderate resolution spectroradiometer data
.
Remote Sensing of Environment
54
,
127
140
.
Hantson
S.
,
Chuvieco
E.
,
Pons
X.
,
Domingo
C.
,
Cea
C.
,
Moré
G.
,
Cristobal
J.
,
Peces
J. J.
&
Tejeiro
J.
2011
Cadena de pre-procesamiento estándar para las imágenes Landsat del Plan Nacional de Teledetección
.
Revista de Teledetección
36
,
51
61
.
Hastenrath
S.
1994
Recession of tropical glaciers
.
Science
265
(
5180
),
1790
1791
.
Herrera-Ossandón
M.
2016
Estimación de las Altitudes de las Líneas de Equilibrio en Glaciares de Montaña Para el último Ciclo Glacial-Interglacial en los Andes de Santiago, Chile Central
.
Thesis
, Universidad de Chile,
Santiago, Chile
.
Herzog
S.
,
Martinez
R.
,
Jørgensen
P.
&
Tiessen
H.
2011
Climate Change and Biodiversity in the Tropical Andes
.
Inter-American Institute for Global Change Research, SCOPE, Montevideo
.
Hongkai
G.
,
Zijing
F.
,
Tong
Z.
,
Yuzhe
W.
,
Xiaobo
H.
,
Li
H.
,
Xicai
P.
,
Ze
R.
,
Xi
C.
,
Wenxin
Z.
&
Zheng
D.
2021
Assessing glacier retreat and its impact on water resources in a headwater of Yangtze River based on CMIP6 projections
.
Science of The Total Environment
765
,
142774
.
Instituto Nacional de Investigación en Glaciares y Ecosistemas de Montaña [INAIGEM]
2018
Inventario Nacional de Glaciares. Las Cordilleras Glaciares del Perú
.
Huaraz
,
Peru
.
Intergovernmental Panel on Climate Change [IPCC]
2018
Global Warming of 1.5°C. An IPCC Special Report on the Impacts of Global Warming of 1.5°C above pre-industrial levels and related Global Greenhouse Gas Emission pathways, in the context of strengthening the global response to the threat of climate change, sustainable development, and efforts to eradicate poverty (V. Masson-Delmotte, P. Zhai, H.-O. Pörtner, D. Roberts, J. Skea, P.R. Shukla, A. Pirani, W. Moufouma-Okia, C. Péan, R. Pidcock, S. Connors, J. B. R. Matthews, Y. Chen, X. Zhou, M.I. Gomis, E. Lonnoy, T. Maycock, M. & T. Waterfield, eds.). Report IPCC SR1.5, Cambridge University Press, Cambridge, UK and NY, USA
.
Jaafari
A.
,
Najafi
A.
,
Pourghasemi
H. R.
,
Rezaeian
J.
&
Sattarian
A.
2014
GIS-based frequency ratio and index of entropy models for landslide susceptibility assessment in the Caspian Forest, northern Iran
.
International Journal of Environmental Science and Technology
11
(
4
),
909
926
.
Kaser
G.
&
Osmaston
H.
2002
The Nature of Tropical Glaciers. Tropical Glaciers
.
Cambridge University Press
,
Cambridge
.
Keshri
A.
,
Shukla
A.
&
Gupta
R.
2009
ASTER ratio indices for supraglacial terrain mapping
.
International Journal of Remote Sensing
30
(
2
),
519
524
.
López-Moreno
J.
,
Fontaneda
S.
,
Bazo
J.
,
Revuelto
J.
,
Azorin-Molina
C.
,
Valero-Garcés
B.
,
Morán-Tejada
E.
,
Vicente-Serrano
S.
,
Zubieta
R.
&
Alejo-Cochachín
A.
2014
Recent glacier retreat and climate trends in Cordillera Huaytapallana, Peru
.
Global and Planetary Change
112
,
1
11
.
Muñoz
R.
,
Huggel
C.
,
Drenkhan
F.
,
Vis
M.
&
Viviroli
D.
2021
Comparing model complexity for glacio-hydrological simulation in the data-scarce Peruvian Andes
.
Journal of Hydrology: Regional Studies
37
,
100932
.
National Oceanic and Atmospheric Administration [NOAA]
n.d.
Southern Oscillation Index (SOI)
.
Available from: https://www.ncei.noaa.gov/access/monitoring/enso/soi (accessed November 2022)
.
Rau
P.
,
Condom
T.
&
Lavado
W.
2013
Spatio-temporal analysis of monthly temperature in the mountainous regions of Peru. An approach for NCEP NCAR Reanalysis data correction
.
Proceedings of the IAHR World Congress
12
,
10602
10612
.
Rau
P.
,
Bourrel
L.
,
Labat
D.
,
Frappart
F.
,
Ruelland
D.
,
Lavado
W.
,
Dewitte
B.
&
Felipe
O.
2018
Hydroclimatic change disparity of Peruvian Pacific drainage catchments
.
Theoretical and Applied Climatology
134
(
1–2
),
139
153
.
Rivera
A.
,
Bown González
F.
,
Napoleoni
F.
,
Muñoz
C.
&
Vuille
M.
2017
Manual Balance de Masa Glaciar
.
Centro de Estudios Científicos, CECs
.
Rojas
T.
,
Quincey
D.
,
Rau
P.
,
Horna
D.
&
Abad
J.
2021
Adapting to receding glaciers in the tropical Andes
.
Eos
102
.
Rukundo
O.
&
Cao
H.
2012
Nearest neighbour value interpolation
.
International Journal of Advanced Computer Science and Applications
3
(
4
),
1
6
.
Salcedo
A.
2011
Estimación de área Cubierta de Nieve en Cuencas con Elevado Aporte de Fusión Utilizando Datos ERS-2
.
Magister Thesis
, Universidad Nacional de Córdoba,
Argentina
.
Salzmann
N.
,
Huggel
C.
,
Rohrer
M.
,
Silverio
W.
,
Mark
B.
,
Burns
P.
&
Portocarrero
C.
2013
Glacier changes and climate trends derived from multiple sources in the data scarce Cordillera vilcanota region, southern Peruvian Andes
.
The Cryosphere
7
,
103
118
.
Schauwecker
S.
,
Rohrer
M.
,
Huggel
C.
,
Endries
J.
,
Montoya
N.
,
Neukon
R.
,
Perry
B.
,
Salzmann
N.
,
Schwarb
M.
&
Suarez
W.
2017
The freezing level in the tropical Andes, Peru: An indicator for present and future glacier extents
.
Journal of Geophysical Research: Atmospheres
122
,
10
.
Seehaus
T.
,
Malz
P.
,
Sommer
C.
,
Lippl
S.
,
Cochachin
A.
&
Braun
M.
2019
Changes of the tropical glaciers throughout Peru between 2000 and 2016 – Mass balance and area fluctuations
.
The Cryosphere
13
(
10
),
2537
2556
.
SENAMHI
2021a
Climas del Perú. Mapa de clasificación climática del Perú. Servicio Nacional de Meteorología e Hidrología del Perú, Lima, Perú. Available from: http://idesep.senamhi.gob.pe/geonetwork/srv/api/records/9f18b911-64af-4e6b-bbef-272bb20195e4/attachments/Resumen%20ejecutivo%20Climas%20del%20Peru%CC%81.pdf (accessed November 2022)
.
SENAMHI
2021b
Catalogo de Metadatos Cartograficos IDESEP. Mapa Anual de Precipitación. Mapa Anual de Temperatura Máxima. Mapa Anual de Temperatura Mínima (1981–2010). Available from: https://idesep.senamhi.gob.pe/geonetwork/srv/spa/catalog.search#/metadata/d4d2a77b-cbc5-41d0-8874-b32848acef57 (accessed November 2022)
.
Shannon
C. A.
1948
Mathematical theory of communication
.
Bell System Technical Journal
27
(
3
),
379
423
.
Soto
C.
,
Zúñiga
J.
,
Paucar
J.
,
Jiménes
M.
,
Ibarra
M.
,
Narváez
A.
&
Paucar
S.
2021
Multi-temporal analysis of the glacier retreat using Landsat satellite images in the Nevado of the Ampay National Sanctuary, Peru
.
Journal of Sustainable Development of Energy, Water and Environment Systems
10
(
1
),
1080380
.
Takahashi
K.
,
Mosquera
K.
&
Reupo
J.
2014
El Índice Costero El Niño (ICEN): Historia Y Actualización
.
Instituto Geofísico del Perú
,
Lima, Perú
. .
Vlcko
J.
,
Wagner
P.
&
Rychlíková
Z.
1980
Evaluation of regional slope stability
.
Mineralia Slovaca
3
(
12
),
275
283
.
Vuille
M.
,
Carey
M.
,
Huggel
C.
,
Buytaert
W.
,
Rabatel
A.
,
Jacobsen
D.
,
Soruco
A.
,
Villacis
M.
,
Yarleque
C.
,
Elisson
O.
,
Condom
T.
,
Salzmann
N.
&
Sicart
J.-E.
2018
Rapid decline of snow and ice in the tropical Andes – Impacts, uncertainties and challenges ahead
.
Earth-Science Reviews
176
,
195
213
.
World Meteorological Organization [WMO]
2017
WMO Guidelines on the Calculation of Climate Normals
.
World Health Organization
,
Geneva
.
Yang
Z.
&
Qiao
J.
2009
Entropy-based hazard degree assessment for typical landslides in the Three Gorges Area, China
. In:
Environmental Science and Engineering
(Wang, F. & Li, T., eds.)
.
Springer
,
Berlin
, pp.
519
529
.
Yang
J.
,
Ding
Y.
,
Liu
S.
&
Tan
C.
2015
Vulnerability of mountain glaciers in China to climate change
.
Advances in Climate Change Research
6
(
3–4
),
171
180
.
Zhang
Y.
,
Enomoto
H.
,
Ohata
T.
,
Kitabata
H.
,
Kadota
T.
&
Hirabayashi
Y.
2016
Projections of glacier change in the Altai Mountains under twenty-first century climate scenarios
.
Climate Dynamics
47
,
2935
2953
.
Zubieta
R.
&
Lagos
P.
2010
Estudio de la dinámica superficial de glaciares en la Cordillera Huaytapallana por sensoramiento remoto: Período 1976-2006 Huancayo – Perú
.
Sociedad Geológica del Perú
9
,
351
356
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).

Supplementary data