Mapping Three Decades of Changes in the Tropical Andean Glaciers Using Landsat Data Processed in the Earth Engine

Abstract

The fast retreat of the tropical Andean glaciers (TAGs) is considered an important indicator of climate change impact on the tropics, since the TAGs provide resources to highly vulnerable mountain populations. This study aims to reconstruct the glacier coverage of the TAGs, using Landsat time-series images from 1985 to 2020, by digitally processing and classifying satellite images in the Google Earth Engine platform. We used annual reductions of the Normalized Difference Snow Index (NDSI) and spectral bands to capture the pixels with minimum snow cover. We also implemented temporal and spatial filters to have comparable maps at a multitemporal level and reduce noise and temporal inconsistencies. The results of the multitemporal analysis of this study confirm the recent and dramatic recession of the TAGs in the last three decades, in base to physical and statistical significance. The TAGs reduced from 2429.38 km² to 1409.11 km² between 1990 and 2020, representing a loss of 42% of the total glacier area. In addition, the time-series analysis showed more significant losses at altitudes below 5000 masl, and differentiated changes by slope, latitude, and longitude. We found a more significant percentage loss of glacier areas in countries with less coverage. The multiannual validation showed accuracy values of 92.81%, 96.32%, 90.32%, 97.56%, and 88.54% for the metrics F1 score, accuracy, kappa, precision, and recall, respectively. The results are an essential contribution to understanding the TAGs and guiding policies to mitigate climate change and the potential negative impact of freshwater shortage on the inhabitants and food production in the Andean region.

1. Introduction

Glacier meltwater plays a vital role in providing water resources to populations downstream of tropical glacier mountain ranges. However, glaciers are decreasing in extent and mass during the dry season at an accelerated rate [1,2,3], due to the increase in global temperature caused by anthropogenic and natural causes [4,5,6]. As a result, glacier retreat is an essential indicator of climate change [2] and climate variability [7,8]. The tropical Andean glaciers (TAGs) are undergoing a rapid coverage reduction, which can economically impact populations in the tropical Andes, with effects on agriculture, drinking water, electricity generation, ecosystem integrity [9,10,11,12], cultural assets, and tourism activities [13,14,15], due to their cultural and tourism value. In addition, the accelerated retreat of glaciers can generate an increase in glacial lake outburst flood (GLOF) risks [7,16].

Most high mountains with glaciers have experienced glacier shrinkage [17]. Understanding the current status of glacier retreat is vital for better management of water resources [18] and adaptation to climate change. Many tropical glaciers, such as those in Peru and Bolivia, are critical buffers against reduced precipitation during the dry season. Ninety-nine percent of tropical glaciers are located in the Andes of South America, extending over Venezuela, Colombia, Ecuador, Peru, Bolivia, Chile, and Argentina [19].

The direct mapping of glaciers is difficult and dangerous due to the remoteness and inaccessibility of the terrain and challenges for conducting extensive fieldwork [20]. The application of remote sensing in glaciers began in the 1970s, initiating a new phase in glacier monitoring and solving to some extent the limitations and difficulties of glacier mapping [19]. The data obtained by remote sensors provide valuable, constantly updated, and almost real-time information on glaciers and associated geographic features [19,20,21].

At global and regional scales, the extent of glaciers can be derived through automated classification algorithms applied to multispectral satellite data [20], such as those from Landsat missions with more than three decades of freely available global information (http://landsat.usgs.gov/ (accessed on 1 March 2021)). Automated methods for glacier classification from multispectral data have been developed and tested in the last decade, and have proven to be suitable for detecting clean or slightly dirty glacier ice and fresh snow [20,22,23]. However, the efforts of mapping glacier coverage do not have a unique record of glacier area. They contain information from different dates on each glacier [24], making their temporal comparison difficult, as in the Global Land Ice Measurements from Space (GLIMS) and World Glacier Monitoring Service (WGMS) initiatives [19,20]. Therefore, it is necessary to have a comparable temporal record of glacier areas to observe their response to climate change and variability [21].

In recent years, the availability, accessibility, and temporal and spatial continuity of satellite data (such as Landsat and Sentinel) [25,26,27,28] and the development of cloud computing platforms, such as Google Earth Engine (GEE) [29], have created an opportunity to map the historical evolution of land cover change, such as water bodies and glaciers [16,18,21,26,30,31,32,33,34,35,36]. Despite the progress made in recent years in the assessment of glacier areas, there are still many cases in which the automatic delimitation of glacier boundaries in satellite images is complex, error prone, or time consuming [13,20,37]. In most cases, glacier cover change estimates cover 5- or 10-year periods [10,35,38].

This research estimates the annual extensions of the entire tropical Andean glaciers in the last three decades, based on a semi-automated classification method and the complete Landsat satellite images historical archives using the GEE platform. Unlike previous studies, this work’s spatial scope covers the entire tropical Andes, temporally comprising a 36-year time series by using an innovative methodological approach to map exposed glaciers based on cloud computing technologies to retrieve and process large quantities of data and the application of post-classification filters. Based on its yearly results, a detailed analysis of all the annual changes of the TAGs is presented, considering site morphometric characteristics (altitude, aspect, and slope) and geographic location.

2. Materials and Methods

2.1. Study Area

The study area refers to the tropical Andean glaciers (TAGs), geographically located between the Tropic of Cancer and the Tropic of Capricorn (between latitudes 23°26′13.3″N and 23°26′13.3″S) within the Intertropical Convergence Zone (ITCZ) [2,19]. According to the Randolph Glaciers Inventory (RGI), the TAGs have a surface area of 2337.55 km² representing approximately 99% of all tropical glaciers [19,39]. According to the classification of glaciers from a climatic approach [40,41], the TAGs contain four climatic zones: Northern Wet Outer Tropics, Southern Wet Outer Tropics, Dry Outer Tropics, and Inner Tropics (Figure 1).

For the delimitation of the TAGs study area, we used the RGI (Version 6.0) [39], records to which we applied a catchment area or buffer of 1.5 km around each glacier to include the surrounding areas. We also visually identified and included glacier areas that were not mapped by the RGI in southern Peru, western Bolivia, and northern Chile, totaling an area addition of approximately 15%. The scope of our research precisely concerns the clean or marginally contaminated glacier surfaces. According to Herreid and Pellicciotti (2020) [42], debris-covered glaciers have a very low representativeness, accounting for approximately 5.4% of the total tropical glacial areas. Hence, the debris-covered glacial areas have not been dealt with in this re-search, for they also require specific methods for their mapping, as reported by Shukla et al. (2017) [43]. Finally, we ought to mention as well that glacier surges are not approached in this work. Contrary to other Andean glaciers [44], where these events are more commonly observed, glacier gains in the TAGs, usually resulting from surges, are sporadic [8] and occur to a minor extent.

2.2. Landsat Data

As the present study is an analysis of glacier coverage change at a multiannual level, the following data corresponding to the period 1985–2020 [26] were obtained from the image collections of the Landsat sensors: Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM +), and Operational Land Imager (OLI), on board the Landsat 5, Landsat 7, and Landsat 8 satellites, respectively [45]. We processed Landsat Collection 1 Level 1 and Tier 1 surface reflectance (SR) products [46], already orthorectified and with a spatial resolution of 30 m in the multispectral bands. The Landsat image collections used were accessed and processed in the Google Earth Engine (GEE) platform [29] (https://earthengine.google.com/ (accessed on 20 June 2021)).

2.3. Landsat Data Processing Workflow

Figure 2 summarizes the methodological workflow of our approach for the TAGs’ change detection. It comprises six stages: selection and masking of clouds and cloud shadows in the satellite images, calculation of the Normalized Difference Snow Index (NDSI), creation of annual Landsat mosaics, annual glacier classification, post-classification, and finally, accuracy assessment.

2.4. Images Selection and Cloud Mask

To select multiple images for each analyzed year, we carried out a sequence of procedures proposed by [20,47]. All Landsat images with cloud cover equal to or less than 70% acquired during the three-decade analyzed period were filtered from the Landsat archive and selected (Table S1 and Figure S1). After selecting the images based on the procedures above, we masked the clouds and cloud shadows. For cloud masking, we used the CloudScore technique [48,49]; in the case of cloud shadow, we used the Temporal Dark Outlier Mask (TDOM) algorithm [49,50] and the Band Quality Assessment (BQA) information available in the Landsat Collection. The BQA was not used for cloud masking, since clouds are overestimated in the glacier edge areas. We implemented these techniques in the GEE platform, which have been thoroughly described by the MapBiomas Project [51].

2.5. Normalized Difference Snow Index (NDSI)

The Normalized Difference Snow Index (NDSI) [52] is a widely used index to separate snow from other coverages [49], and therefore it is also useful to identify glacier coverage. We applied the NDSI index to each selected and masked image in the previous steps, based on Equation (1).

$$\mathrm{NDSI}=\left(\frac{\mathrm{Green}-\mathrm{SWIR}}{\mathrm{Green}+\mathrm{SWIR}}\right)$$

(1)

2.6. The Compositing and Mosaicking of the Annual Images

We created annual image composites for each year of the time series (36 years) using GEE’s image-reducing techniques, having as an input all the cloud-masked Landsat scenes available in GEE’s Landsat collection with less than 70% cloud cover. We applied the image reducers to the NDSI, RED, and NIR bands, as these variables are most suitable for glacier classification [53]. The NDSI values in each pixel of each year were reduced to a band of $NDS{I}_{min}$ (annual minimum NDSI), and the annual reduction in the RED and NIR bands were based on the 25th NDSI percentile. This percentile value was defined on a trial and error basis, aiming to reduce cloud cover and eliminate noise at the same time. In brief, the time-series data dimensionality was reduced by being condensed into a three-band artificial image for each year with the minimum snow cover, composed of manifold pixels from different dates throughout each of the analyzed years [51]. These three bands are the $NDS{I}_{min}$, $RE{D}_{minsnow},$ and $NI{R}_{minsnow}$, and the two latter ones were calculated using Equations (2) and (3), respectively.

$$RE{D}_{minsnow}=\mathrm{median}\left(\mathrm{RED} \mathrm{with} \mathrm{the} 25\mathrm{th} \mathrm{NDSI} \mathrm{percentile}\right)$$

(2)

$$NI{R}_{minsnow}=\mathrm{median}\left(\mathrm{NIR} \mathrm{with} \mathrm{the} 25\mathrm{th} \mathrm{NDSI} \mathrm{percentile}\right)$$

(3)

The $RE{D}_{minsnow}$ reduction is interpreted as the RED band with minimum snow in the year, and the $NI{R}_{minsnow}$ as the NIR band with minimum snow. We used these criteria to reduce glacier confusion with temporary snow [20]. A similar work about exploiting a time-series stack of images and its corresponding reduction was reported by Winsvold et al. (2016) [54].

2.7. Annual Glacier Classification

Annual glacier classification for the 1985–2020 period was obtained using an empirical decision tree (EDT) constructed using samples of $NDS{I}_{min}$, $RE{D}_{minsnow}$, and $NI{R}_{minsnow}$ variables. The sampling strategy was performed for both Landsat 5 and 7 in 2000, and for Landsat 8 in 2020, respectively, and a visual sample collection was performed for the classes: glacier, water, and other coverages. Samples were collected in the entire study area considering the classes spectral variation (e.g., pixels in different slopes) and prioritizing pixels where there is greater confusion among the three classes (e.g., shaded pixels). In total, 9800 samples were collected for each sensor. Then, we conducted a statistical analysis of the sample set for each year using a boxplot violin diagram (Figure 3 and Figure 4) in order to find the optimal thresholds.

The decision tree used is shown in Figure 5. We empirically defined the decision tree nodes (1–1, 2–2, and 3–2) with thresholds α, β, and γ defined for each TAG classification region, classification year, and sensor type—TM, ETM+, and OLI (see Figure S2 for the classification regions details). EDT nodes were constructed with the objective of detecting glacier ($\mathrm{variables} NDS{I}_{min}$) and eliminating confusion with water bodies and other coverages ($RE{D}_{minsnow}$ and $NI{R}_{minsnow}$ variables) [53].

The average value and standard deviation of the threshold α are 0.20 and 0.082, respectively, while the β and γ thresholds were set as constant at 0.11 and 0.10, respectively, which are very similar values to those used by Macander et al. (2015) [53]. The thresholds were adjusted using visual inspection in each classification region and year.

2.8. Post-Classification

To remove noise and information gaps due to cloud persistence, we designed and applied the following sequence of five post-classification filters:

• Gap-fill filter to replace information gaps due to clouds with consistent temporal classification values;
• Temporal filter to identify spurious class transitions (e.g., change in short periods from glacial to non-glacial);
• Frequency filter to identify stable pixels [51];
• Temporal permanence filter;
• Spatial filter to remove salt-and-pepper classification noise without changing the boundary of the classified patches.

These filters were applied to reconstruct land cover and land use classes over large time-series, as in this study [51]. Below, we provide the details of each filter’s utility.

The filter sequence begins with the gap-fill technique (Figure 6). Severely cloud-affected regions commonly present information gaps in the annual TAGs classifications in long time-series, as a result of the cloud masked regions. The gap-fill tool reduces this effect by replacing ‘no data’ pixels with the temporally closest class value starting from the near past and moving forward in time (example Figure 6). The use of this filter allows for the achievement of comparable glacier areas in the time series [46].

Next, a temporal filter was applied [51]. The temporal filter uses sequential rankings in a one-way moving window of 3, 4, or 5 years’ kernel sizes to identify the temporally inconsistent transitions (e.g., glacier/non-glacier expansion or contraction when the historical record shows the continuity of either class) based on the intermediate-, first-, and last-year rules (Figure 7). The first- and last year rules were initially applied to respectively eliminate time inconsistencies in the years 1985 and 2020 (Figure 7a,b), followed by the intermediate years’ rules (i.e., 1986 to 2019) in the following order: the three-year rule (Figure 7c), the four-year rule (Figure 7d), and finally the five-year rule (Figure 7e).

Next, the frequency filter was applied. This filter improves the classification according to the historical occurrence of the class over the analyzed 36-year period. For this purpose, it considers the frequency percentage of the glacier class in the time series, from which it updates the ‘non glacier’ class occurrences, which can be safely considered classification errors in most cases. This study used a 70% frequency threshold; in other words, occurrences smaller than or equal to 30% of ‘non-glacier’ values in the time series were replaced by glacier values. This filter was applied only to the first and intermediate years; years at the end of the time series were excluded to prevent the loss of information of transitions occurring at the end of the time series (e.g., glacier loss in recent years) (Figure 8).

Based on the GLIMS definition, glaciers are defined as permanent ice or snow-mass [20]. We proposed applying a temporal permanence filter to separate temporary snow from glacier coverage (Figure 9).

Finally, to eliminate salt-and-pepper classification errors [37], a spatial filter based on GEE’s ‘connectedPixelCount’ function was applied. The spatial filter identifies sets of pixels that share the same pixel value based on their neighborhood. Consequently, only pixels that do not meet a predefined connection criterion (i.e., minimum number of identical-value connected pixels) are defined as isolated pixels and reclassified [51]. We considered a minimum of 5 grouped pixels for glaciers, representing a minimum area of approximately 0.5 ha.

2.9. Accuracy Assessment

Thematic accuracy analysis is the primary way to evaluate the quality of maps [55]. The accuracy assessment of annual maps in this study used a random sampling of 7299 points. This sample size guarantees a confidence level of 95%, with a normal distribution (Z = 1.96), a level of significance of 5%, and a maximum tolerable error below 1% (e = 0.5%). A minimum of 200 sample points per class and country was assured. Each point was interpreted and assigned a class using the QGIS Earth Engine plugin (qgis-earthengine-plugin) [56] and the thematic map accuracy evaluation plugin AcATaMa [57]. The reference dataset was obtained by means of visual inspection of the sample points, conducted by independent expert image interpreters, who relied on an RGB image composite consisting of the following bands: SWIR with the 25th NDSI percentile, NIR with the 25th NDSI percentile, and RED with the 25th NDSI percentile, for each analyzed year. This allowed us to calculate the error matrix and other accuracy metrics. We also conducted accuracy analyses according to the work of Olofsson et al. (2014), on good practices for estimating area and assessing accuracy of land change [58].

The accuracy of the classifications was tested using a set of binary metrics. The possible outcomes are true positive (TP), true negative (TN), false positive (FP), and false-negative (FN). The metrics tested are given in the following Equations (4)–(7) [59,60]:

$$\mathrm{Precision}: \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}$$

(4)

$$\mathrm{Recall}: \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$$

(5)

$$\mathrm{Accuracy}: \frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+ \mathrm{FN}}$$

(6)

$$\mathrm{F}1\mathrm{score}=2\ast \frac{\mathrm{Precision}\ast \mathrm{Recall}}{\mathrm{Precision}+\mathrm{Recall}}=\frac{2\ast \mathrm{TP}}{2\ast \mathrm{TP}+\mathrm{FP}+\mathrm{FN}}$$

(7)

3. Results

3.1. The Area Statistics of the TAGs

We estimated the annual area of the TAGs from 1985 and 2020, comprising 36 years of annual mapping. This annual TAGs time series is the largest one currently available, to the best of our knowledge. In order to have a good quality base year, we analyzed the results from the year 1990 onwards, using ten-year intervals until 2020 (

Table 1. Glacier-mapping summary statistics estimated for decennial periods.

Period (# Year)	Minimum (km2)	Maximum (km2)	Mean	Standard Deviation (s)
			(km2/Year)
1991–2000	2089.20	2364.80	2220.35	95.58
2001–2010	1815.57	2069.41	1928.37	89.03
2011–2020	1409.12	1807.83	1670.35	123.14

). The smallest TAGs (minimum) for the period 2011–2020 had a total surface of 1409.12 km² (Table 1).

We analyzed the results of the TAGs coverage at the country level using the Global Administrative Unit Layers (GAUL) country boundary available in GEE (https://developers.google.com/earth-engine/datasets/catalog/FAO_GAUL_2015_level2 (accessed on 14 January 2021)). We found that the most extensive glacier area in 2020 was found in Peru and Bolivia, with 72.76% and 20.35% of the total TAGs surface, respectively. Ecuador, Colombia, Chile, and Argentina combined represented 6.89% of the TAGs coverage in the Andes Cordillera (3.89%, 2.18%, 0.78%, and 0.04%, respectively). Venezuela has less TAGs coverage, with a percentage of less than 0.01%, accounting for nearly 0.03 km² (

Table 2. Glacier-mapping summary statistics estimated for decennial periods per country.

Country	Glacier Area (km2)
	1990	2000	2010	2020
Peru	1743.45	1503.67	1320.85	1025.29
Bolivia	499.58	434.75	364.02	286.69
Ecuador	86.10	73.70	68.52	54.78
Colombia	77.21	60.18	49.21	30.73
Chile	20.77	15.73	12.15	10.96
Argentina	1.15	0.73	0.68	0.63
Venezuela	1.12	0.44	0.14	0.03
Total	2429.38	2089.20	1815.57	1409.12

and Figure 10a).

The TAGs occur at different altitudes (Figure 10b). We verified that, in 2020, most of the glaciers are distributed between 4801 and 5600 masl, representing 82.63% of the total TAGs, while approximately only 1.72% at altitudes less than or equal to 4800 masl, and 15.64% are found above 5600 masl.

The results show that, in 2020, 19.47% of the TAGs face north, 21.56% east, 31.30% south, and 27.67% face west (Figure 10c). Most TAGs (73.31%) occur in areas with slope ranges from 11° to 40° (Figure 10d). Flat and low slopes terrains (i.e., from 0° to 10°) contain 11.57% of the TAGs, while areas with pronounced slopes (above 41°) concentrate 15.11% of the TAGs (Figure 10d).

The results show that, in the year 2020, TAGs are located between latitudes −24° and 11° (Figure 11a). The largest glacier areas occur in the range between −14° to −10°, which contains 44.75% of the total TAGs; the Peruvian Cordilleras (mountain ranges), Blanca, Vilcanota, Vilcabamba, and Urubamba, are located within these latitudes. The longitude range of the TAGs varies from −78° to −67° (Figure 11b).

3.2. The Spatial and Temporal Patterns of the TAGs

The total area of the TAGs considerably reduced between 1990 and 2020, following a linear decline pattern. We found an average glacier area reduction of 28.42 km²/year with an R² of 0.98. In 1990, the mapped area of the TAGs reached 2429.38 km², reducing by 42% by 2020 to an area of 1409.11 km² (Figure 12a). Of the last three decades, the decade of 2011–2020 witnessed the greatest loss of glacier areas (associated with the decade’s minimum surface value) (Figure 12b). Figure 13 presents a map for the Vilcanota Cordillera (mountain range) in Peru. The retreat buffer width varies along the TAGs edges and over the years. The smallest TAGs (i.e., those equal to or less than 10 ha) lost more than 90% of their area in the period from 1990 to 2020 (Figure 13).

For the analysis of glacier loss or retreat, 1990 was adopted as the base year, and annual maps of losses were created (Table S2). The area losses of the TAGs by country are shown in percentages in Figure 14, where, in the year 2020, Peru, Bolivia, Ecuador, Colombia, Chile, Argentina, and Venezuela presented retreats of 41.19%, 42.61%, 36.37%, 60.19%, 47.24%, 45.47%, and 96.93%, respectively, with respect to the base year of 1990.

There are three distinct patterns of TAGs retreat when area loss is disaggregated by country (Figure 14, Table S3). First, Venezuela, Argentina, and Chile showed a rapid decline in the TAGs between 1990 and 1995, with approximately 20% of glacier area loss during that specific period (Figure 14), although this can be due to the low satellite image quality in those initial years. By 2020, Venezuela’s TAGs losses exceeded 90%, but in the case of Chile and Argentina, the losses reached 45.47% and 47.24%, respectively, since the year 1990. A second pattern is found in Colombia’s TAGs, which showed an acceleration from 1995 onwards, reaching a loss of 60.19% by 2020. The third pattern of gradual linear TAGs loss occurred in Peru, Bolivia, and Ecuador (Figure 14). Ecuador showed a lower rate of retreat.

We analyzed the glacier area loss by different latitudes and longitudes, with respect to the average loss of all TAGs (Tables S4 and S5). We found that the latitudes farther south (−24° to −20°) have higher losses with regard to the total mean loss, with more than 62.54% of the glacier loss observed during the analyzed time series (Figure S3), however these glaciers represent less than 1% of the total area of the TAGs in 2020. The TAGs in the southern latitudes, from −9° to 0°, had a reduction of 31.61%, lower than the total average (42.00%); these latitudes account for 33.6% of the total area of the TAGs in 2020. At the same time, we verified that glaciers with latitudes ranging from 1° to 11° north had an accelerated glacier retreat in recent years (since 1996), with a 54.27% loss from 1990 to 2020 (Table S4). We also observed differentiated change patterns when analyzing glaciers at different geographic longitudes. The glaciers between longitudes −74° and −73° experienced changes of 59.94%, well above the total average (i.e., mean total retreat TAGs over the period 1990–2020), while longitudes −78° and −77° had losses of 31.68%, below the total average (Table S5 and Figure S3).

The area loss of the TAGs with regard to aspect (i.e., different glacier slope orientations) was lower in the south and east (i.e., 41.59% and 40.63%, respectively) and higher in the north and west orientations (42.86%) (Table S6 and Figure S4). Our analyses also estimated that glaciers with a larger slope than 80° rapidly reduced by more than 30% in 1996, with regard to the 1990 TAGs state. The retreat then slowly decelerated in the following years, until it reached a total loss of 41.56% by 2020. On the other hand, the TAGs with slopes ranging from 61° to 70° had a greater acceleration of loss from 2015 onwards, representing a loss of 46.37% by 2020. Nevertheless, in other cases, we could not find a slope range that controls the retreat of the TAGs area (Table S7 and Figure S5). In summary, we also observed that the TAGs areas below 4600 masl lost an average of 92.41% of their extension in 1990, while at altitudes above 6000 masl, the losses have been lower, 9.73% on average (Table S8 and Figure S6).

Due to the importance of the Amazon basin in South America, we verified that the glaciers in this region have an approximate area of 869.59 km² in 2020, representing 61.71% of the TAGs total area. We calculated and compared the evolution of the TAGs area loss in the period from 1990 to 2020 in the Amazon basin, and concluded that it has lost 43.07% of its glacier area, unlike other non-Amazonian basins, which have lost 40.18% along the analyzed time series. We found a more significant acceleration in recent years, approximately from 1995 onwards, where the loss of the Amazon basin exceeds that of other basins (Table S9 and Figure S7).

3.3. Accuracy of the TAGs Mapping

We conducted the accuracy assessment of the TAGs annual maps over the entire time series of the reference dataset generated for the study area based on 7299 random points and visual inspection, as previously described. The average annual accuracy of the original classification was 95.29% (Table S10). After the application of the filters, we observed a gradual improvement in map accuracy to 95.30%, 96.31%, 96.37%, 96.27%, and 96.32% for the classifications relying on the cumulative sequential use of the gap-fill, temporal, frequency, temporal permanence, and spatial filters, respectively (Figure 10, Table S11). Although there was a slight decrease in the average accuracy, the standard deviation of the time-series classification accuracy was reduced from 1.24 to 0.44 with the filtering techniques, hence exemplifying the utility of the post-classification filtering routines. The marginal gain of the filtering techniques (i.e., 1.03 ± 0.44) is justified with the better spatial and temporal consistency in the classification of glacier coverage for multiannual comparisons. The average classification accuracy of the final filtered product attained 92.81%, 96.32%, 90.32%, 97.56%, and 88.54% for the validation metrics F1 score, accuracy, kappa, precision, and recall, respectively. The kappa index and the accuracy showed the highest values, while recall was the lowest (Figure S8).

This indicates an omission error of 11.46% and a commission error of 2.44%. The omission errors mostly occurred in shaded areas and in small glaciers, while the commission errors took place at the glaciers’ borders. The highest omission and commission errors were found in the first and last years (1985 and 2020), caused by the low amount of Landsat images in the early years (Figure S3), since the filters still have limitations in the early and late years (Figure 15a). The accuracy reduction valleys in given moments throughout the time series (Figure 15a) are related to the mosaic quality in those years (e.g., 2002), due to the limited availability of cloud-free images, which had an important effect on the accuracy of those years’. The first five years of the time series were also impacted by the limited image availability.

The accuracy analyses considering the area proportions according to [58] are presented in Table S12, which also shows the annual accuracy metrics disregarding these proportions. The difference between these two types of mean accuracy is only 0.42%. This minor discrepancy is probably due to the fact that we are dealing with a binary classification.

4. Discussion

The thematic accuracy of the obtained results present acceptable values for the multitemporal maps, with an average accuracy of 96.32%. Similar studies in the literature reported inferior accuracies in general. Albert (2002) [22] conducted several semi-automatic experiments of ice-area classification in the Quelccaya Ice Cap, Peru, and compared their results to a hand-digitization. The employed methods comprised pixel-based supervised and unsupervised techniques, including non-parametric approaches, such as fuzzy logic and empirical decision trees (thresholds). Although the hand digitization reached an accuracy of 99.8%, the semi-automated results obtained accuracies varying from 73.7% (fuzzy technique) to 96.7% (Spectral Angle Mapper—SAM). The mean accuracy was 84.8%. Specifically regarding the empirical results for the decision trees, the accuracies ranged from 93.0% to 93.9%. In their turn, Rittger et al. (2013) [60] compared the accuracies of three MODIS global products, namely, the binary and fractional snow product (MOD10A1), based on the NDSI, and another fractional snow product (MODSCAG) based on spectral mixture analysis. The study areas relate to the mountain ranges located in the USA (Sierra Nevada, Rocky Mountains, the Upper Rio Grande) and in the Nepal Himalaya, southern Asia, which are not included in tropical regions. The attained precision values ranged from 71.2% to 99.9%, and the mean precision values for the four study areas were 95,8% for the binary MOD10A1, 98.4% for the fractional MOD10A1, and 99.2% for the fractional MODSCAG. We obtained a mean precision (considering all TAGs and all years of the time series) of nearly 97.6%.

Post-classification filters specially adapted for the detection of glacial cover according to the GLIMS definition were also applied [20]. The application of filters helped filling the gap of non-observed areas and correct temporal and spatial inconsistencies, improving the classification accuracy by up to 6.91% in some years. Nevertheless, the gap-filling filters did not improve the map accuracy in the first and last years of the time series due to the smaller size of the filter’s temporal kernel and the limited number of images in the first years of the time series [37]. The first years of the time series were not consistent with the remaining time series due to the limited number of images. As a result, we used the year 1990 as the baseline for the gain–loss analysis and removed the mapping between the years 1985 and 1989. The accuracy assessment, particularly in the intermediate years, indicated that the proposed methodology generates comparable maps, presenting spatial and temporal coherence, useful for glacier change analysis.

As mentioned in [8] and other studies about the mass balance analysis of the tropical Andean glaciers, although there have been some sporadic gains in the TAGs, the general trend is negative, i.e., glaciers are retreating in the Andes in the past decades. Our study corroborates these findings. Detecting clouds and seasonal snow is critical to obtain consistent estimates of glaciers and measure their change over time with satellite imagery [42]. To overcome this obstacle, we refined cloud masking methods (i.e., CloudScore and TDOM) and applied temporal and spatial filters to fill the unobserved gaps due to cloud (using scenes with cloud cover <70%).

The annual reduction in the NIR and RED bands based on the 25th percentile of the NDSI values allowed seasonal snow detection. Our algorithm excluded it as soon as the non-glacial ground became exposed. This procedure enabled us to select the pixels with less snow coverage in the image collection, improving the Minimum NDSI composite used by Huang [49]. The method proposed by Huang [49], based on the minimum reducer in the NIR and RED bands, generates noise when there are pixels where clouds and cloud shadow have not been properly masked out, leading to underestimation of the glacier cover in the classification. Moreover, an appropriate selection of images and spectral bands [8] was also crucial for mapping glaciers to avoid confusion with temporary and seasonal snow cover. Finally, using regional thresholds in the empirical decision tree led to better classification results than a similar study for this region that exclusively relied on a single threshold-based decision tree [7].

When referring to the distribution of tropical glaciers, a worldwide reference is the RGI of the GLIMS initiative [19,39]. However, despite being widely used in many studies, this reference has a significant limitation to calculating the distribution of tropical glaciers and making comparisons of multitemporal glacier loss, since the RGI data are from different dates in each place as it is targeted at producing a glacier inventory. Therefore, in our study, we reconstructed the distribution of the tropical Andean glaciers in the last three decades and we updated the distribution of glaciers to the year 2020. Our new data set allows comparing glacier change across different countries. Peru and Bolivia comprise 72.76% and 20.35% of the total TAGs, while Ecuador, Colombia, Chile, and Argentina account for 3.89%, 2.18%, 0.78%, and 0.04%, respectively. In turn, Venezuela is the country with the smallest glacier extent. The countries with less glacier coverage in 1990 are the ones that have had the most significant glaciers losses in percentage by 2020.

Overall, the TAGs are experiencing a steady retreating trend, and their extent reduced by an average of 42% between 1990 to 2020, which indicates a loss of nearly half of the tropical Andean glaciers extension in the last three decades. This fast loss rate of the TAGs is an unprecedented change in tropical glaciers [8], potentially linked with the ongoing warming of the planet as a consequence of climate change and climate variability [19,38]. The increase in global temperature [61] as well as the variability of sea surface temperatures (SSTs) of the Pacific Ocean can be the potential drivers of the TAGs retreat [8] we presented in this study. Therefore, the retreat of the TAGs constitutes an early case of the need to adapt to climate change to reduce the economic and social impacts associated with the loss of glacier cover and the freshwater they produce [11].

TAGs dynamics are also influenced by altitude, latitude, longitude, slope, aspect, and country. We realized that the tropical Andean glaciers below 5000 masl are the ones that changed the most in the last three decades. This result is consistent with previous studies, such as a glacier in Bolivia and Peru [10,35]. The TAGs located below 5000 masl lost almost 80.25% of their initial area in the year 1990. These lower-height glaciers retreated faster than higher ones, which can be explained by the high negative impact of the net short-wave radiation linked with cloud cover and surface albedo [8,35]. The albedo feedback is suggested as the primary cause of accelerated retreat [8]. It is worth mentioning that the slope variable also influences the glacier loss, since areas with a slope above 80° have lost more glacier cover than average. Therefore, if the current environmental conditions remain, it is very likely that in the future the glacier surface will disappear at low elevations below 5000 masl in the TAGs, accounting for approximately 171.26 km², representing 12.15% of the TAGs in the year 2020.

We also found an accelerated retreat of a glacier in the hillsides belonging to the Amazon basin. Climatic and non-climatic factors could explain this significant change in the Amazon in recent years. One of the non-climatic factors could be the increase in forest burning in recent years in the Amazon, which generates black carbon [62,63] that can accelerate glacial retreat when entering the glacier surface [19,62]. However, it can also be explained by solar radiation in the tropical zone, which is more or less constant throughout the year [8] and lately intensified by climatic changes. Even though the hydrological contribution of glaciers to the Amazon basin is minor, this may influence the rivers fed by the Andean glaciers in the Amazon basin [35], as is the case of the Madeira River in Brazil [64], fed by the glaciers in Peru and Bolivia.

The countries with less TAG coverage, such as Venezuela, Argentina, Colombia, and Chile, are the ones that have lost their glacier coverage more rapidly in the last three decades. Furthermore, this may have a possible relationship with the geographic latitude and longitude of these TAGs. We observed that glaciers with latitudes farther from the equator are the ones that have lost more glacier coverage with respect to the total TAGs mean. In addition, we found that glaciers with latitudes ranging from −9° to 0° lost less area than the total average, and we noted that they are glaciers with more extensive areas in the tropical zone. The glaciers with latitudes closer to the Pacific Ocean, where the Blanca Cordillera (mountain range) in Peru is located, had lower losses than the total average, which could be explained by a possible association of glacier cover loss in the TAGs with the climatic behavior of the Pacific Ocean [8].

Several studies have shown the high (negative) correlation between annual air temperature and glacier mass balance in tropical regions [65]. Notably, the freezing level heights (FLHs), i.e., the elevation where the 0 °C air temperature is found, has been used as a proxy for the TAGs mass balance [15,66,67]. Higher FLH values correspond to glacier mass loss (or more negative mass balance). The continuous tropical FLH increase is basically linked with the Pacific ocean’s sea surface warming trend [67,68]. In some years during the ENSO (El Niño and La Niña events, with recurrence cycles from 2 to 7 years) occurrence, the Pacific ocean’s warming conditions could produce significant changes in the glaciers´ surface extension [67,68], but with different impacts on each country and region. Those sudden and considerable changes of the TAGs surface detected in the previous analyses can be due, in several cases, to ENSO effects, but this is out of the scope of the present research.

As limitations of the proposed methodology, we ought to mention that the data dimensionality reduction depends on the quantity and quality of the input images. The adopted filters do not add information into the classification process, but merely correct the temporal and spatial inconsistencies. The assessment of the glaciers´ retreat was conducted on a 2D basis, being restricted to the observable surface in the images. It would be highly recommendable to advance this research into a volumetric estimate of glacier loss, and this initiative would probably rely on sensors dedicated to measure the ice-sheet thickness (e.g., IceSat-2/Atlas). Finally, the classification methodology is not fully automatic, depending on expert image interpreters and heuristic routines for setting the decision tree thresholds.

5. Conclusions

In this research, we reconstructed the annual extensions of the TAGs for the last three decades, based on a semi-automated method and the historical Landsat series in the GEE platform. This approach built a new and updated dataset of the TAGs based on remote-sensing time-series data. The change detection analysis presented in this study presented a new understanding of annual glacier coverage changes across different zones of the tropical Andean Mountain Ranges. As a result, we provided the first multiannual and spatially detailed analysis covering the total extent of the TAGs. The observed changes highlight the worrisome retreat of glaciers throughout the tropical Andes, which will likely lead to considerable environmental, economic, and cultural problems in this region. In addition, the results represent an important source of information for developing water management strategies capable of meeting the challenges of climate change. The results shown in the present work are a source of data for understanding the relationship between glacier recession and different factors that determine either the progressive loss or permanence of tropical glaciers. Future steps in our research project are to provide the TAG dataset through a dashboard to support better decision making and water resource planning and management, climate change adaptation, and improving glacier monitoring.