Introduction
Cryogenic ecosystems by definition depend on the dynamics of snow and ice. They are undergoing widespread and intense transformations because of climate change, especially in the tropics (Huss et al., Reference Huss, Bookhagen, Huggel, Jacobsen, Bradley and Clague2017; Masiokas et al., Reference Masiokas, Rabatel, Rivera, Ruiz, Pitte and Ceballos2020). Tropical glacier ecosystems are restricted to the highest elevations of mountain ranges (typically above 4,700 m) where snow and ice have accumulated over years, and the ice, melting water and rocky substrates sustain microbiological communities and other meso- and macro-biota (Hotaling et al., Reference Hotaling, Hood and Hamilton2017). They are naturally isolated and increasingly exposed to rising temperatures and changing precipitation patterns that create an imbalance in ice mass, with more ice lost because of melting and ablation than the amount gained by snowfall and rainfall (Sagredo & Lowell, Reference Sagredo and Lowell2012; Rabatel et al., Reference Rabatel, Francou, Soruco, Gomez, Cáceres and Ceballos2013; Veettil & Kamp, Reference Veettil and Kamp2019). Historical losses of glacial ice mass have been documented in several countries, with > 60% having been lost in Colombia and > 50% in Peru over the last 60 years (Masiokas et al., Reference Masiokas, Rabatel, Rivera, Ruiz, Pitte and Ceballos2020). The tropical Andes are considered a global hotspot for cryogenic change (Vuille et al., Reference Vuille, Yarleque, Elison Timm, Carey, Huggel and Salzmann2018), having lost over 1 gigatonne of ice during 2000–2018 (Dussaillant et al., Reference Dussaillant, Berthier, Brun, Masiokas, Hugonnet and Favier2019), and the whole tropics region is expected to lose at least 60% of its glacier ice mass by 2100 if current global climate change trends are maintained (Rounce et al., Reference Rounce, Hock, Maussion, Hugonnet, Kochtitzky and Huss2023).
The functions of glacier ice as a water reservoir, a regulator of local temperature and a substrate for biotic processes are being disrupted by continuing losses of mass (Buytaert et al., Reference Buytaert, Moulds, Acosta, De Bièvre, Olmos and Villacis2017; Stibal et al., Reference Stibal, Bradley, Edwards, Hotaling, Zawierucha and Rosvold2020). The changes in biota and ecosystem processes associated with glacial retreat will have cascading effects as novel ecosystems emerge in recently exposed areas (Llambí et al., Reference Llambí, Melfo, Gámez, Pelayo, Cárdenas and Rojas2021; Rosero et al., Reference Rosero, Crespo-Pérez, Espinosa, Andino, Barragán and Moret2021; Anthelme et al., Reference Anthelme, Carrasquer, Ceballos and Peyre2022) and changes in downstream aquatic habitats, biota and water quality affect human livelihoods in the surrounding landscape (Cauvy-Fraunié et al., Reference Cauvy-Fraunié, Andino, Espinosa, Calvez, Jacobsen and Dangles2016; Huss et al., Reference Huss, Bookhagen, Huggel, Jacobsen, Bradley and Clague2017; Llambí et al., Reference Llambí, Becerra, Peralvo, Avella, Baruffol and Díaz2020). Long-term monitoring in the region indicates that plant diversity in the tropical Andes is particularly vulnerable given the high proportion of endemic species that are dependent on low temperatures and have narrow thermal niches (Cuesta et al., Reference Cuesta, Tovar, Llambí, Gosling, Halloy and Carilla2020, Reference Cuesta, Carilla, Llambí, Muriel, Lencinas and Meneses2023). Andean landscapes are being transformed, including via the contraction of puna and paramo and degradation/loss of high Andean wetlands (Cuesta et al., Reference Cuesta, Llambí, Huggel, Drenkhan, Gosling and Muriel2019).
Monitoring the status and trends of ecosystems is fundamental to informing biodiversity conservation and sustaining the benefits from nature (Bland et al., Reference Bland, Nicholson, Miller, Andrade, Carré and Etter2019). The IUCN Red List of Ecosystems provides a standard protocol for monitoring through a rigorous assessment of the risk of ecosystem collapse. Ecosystem collapse is defined as a large transformational change involving loss of biodiversity and/or major degradation and loss of ecosystem function. The risk status of an ecosystem is represented by the assigned Red List of Ecosystems category of risk that summarizes the probability of collapse based on the diagnosis of qualitative and quantitative indicators (Keith et al., Reference Keith, Rodríguez, Rodríguez-Clark, Nicholson, Aapala and Alonso2013).
Previous Red List of Ecosystems assessments in the tropical Andes have focused on vegetated ecosystem types including tropical alpine vegetation, but have not considered glacial ecosystems (Oliveira-Miranda et al., Reference Oliveira-Miranda, Huber, Rodríguez, Rojas-Suarez, De Oliveira-Miranda, Zambrano-Martinez, Hernandez, Rodriguez, Rojas-Suárez and Hernández2010; Etter et al., Reference Etter, Andrade, Saavedra, Cortés, Moreno, Rueda and Andrade2018). Globally there are few examples of such assessments of cryogenic ecosystems, and those that exist are primarily focused on seasonal snow (Williams et al., Reference Williams, Stott, Camac, White, Burns and Harris2015; Kontula & Raunio, Reference Kontula and Raunio2019), freeze–thaw glacial lakes or permanent snowfields (Murray et al., Reference Murray, Keith, Duncan, Tizard, Ferrer-Paris and Worthington2020).
Here we contribute to filling this gap by presenting a precursory example of a global assessment of all tropical glacier ecosystems. We focus on the glaciers of the Cordillera de Mérida in Venezuela, which have a documented history of decline in number and extent (Ramírez et al., Reference Ramírez, Melfo, Resler and Llambí2020). We use information from local, regional and global studies to apply the Red List of Ecosystems criteria, and we provide a comprehensive ecosystem description supported by knowledge from related ecosystem types and information from local studies, a diagnosis of ecosystem functioning, threats and state of collapse and the ecosystem risk status based on spatial and functional symptoms of collapse.
Study area
The tropical glacier ecosystem of the Cordillera de Mérida is a distinct ecosystem type that formerly occupied several peaks above 4,600 m in Sierra Nevada National Park in Venezuela (Fig. 1; Braun & Bezada, Reference Braun and Bezada2013; Ramírez et al., Reference Ramírez, Melfo, Resler and Llambí2020).
The study area in the Cordillera de Mérida, Venezuela. This assessment focuses on three glaciated peaks: Bolívar, La Concha and Humboldt. These peaks and other historically glaciated areas are encompassed in a single 10 × 10 km cell (highlighted square, area of occupancy = 1). The isolated Mucuñuque Peak at 4,609 m lost its ice before 1930.
Following the IUCN Global Ecosystem Typology, all tropical glacier ecosystems belong to the Ecosystem Functional Group T6.1 Ice sheets, glaciers and perennial snowfields in the Cryogenic functional biome (T6; Keith et al., Reference Keith, Ferrer-Paris, Nicholson, Bishop, Polidoro and Ramirez-Llodra2022). Tropical glaciers share key features with other glaciers, such as the icy substrate, atmospheric deposition of nutrients, a biota dominated by cold-adapted microorganisms that inhabit several habitat compartments on (supraglacial), in (englacial) or below (subglacial) the icy substrate, truncated trophic networks and low productivity and diversity (Anesio et al., Reference Anesio, Lutz, Chrismas and Benning2017; Hotaling et al., Reference Hotaling, Hood and Hamilton2017; Keith et al., Reference Keith, Ferrer-Paris, Nicholson, Bishop, Polidoro and Ramirez-Llodra2022). Tropical glaciers differ from others in their exposure to daily variations in temperature that far exceed the differences in monthly mean temperatures across the year. Precipitation can occur throughout the year as rainfall or snowfall, although it tends to concentrate during the rainy or wet season, which typically extends for > 8 months (Sagredo & Lowell, Reference Sagredo and Lowell2012). Annual precipitation at the highest elevations in the Cordillera de Mérida is estimated to be 1,000–1,200 mm, with 10–17% of the total falling during December–March, and the mean annual temperature is −0.4 °C at the highest station where records are available in the country (4,766 m), with a ratio of diurnal variation to annual variation of 2.15–2.36 (Pulwarty et al., Reference Pulwarty, Barry, Hurst, Sellinger and Mogollon1998; Andressen, Reference Andressen and Grau2007).
The long-term mass balance of the icy substrate of tropical glaciers is therefore probably dominated by interannual fluctuations in precipitation and exposure to solar radiation, but this process can be accelerated by scale and edge effects (Ceballos et al., Reference Ceballos, Euscátegui, Ramírez, Cañon, Huggel, Haeberli and Machguth2006; Andressen, Reference Andressen and Grau2007; Vuille et al., Reference Vuille, Yarleque, Elison Timm, Carey, Huggel and Salzmann2018). The Cordillera de Mérida has lost several glaciated areas in the last 150 years (Braun & Bezada, Reference Braun and Bezada2013), and the glacier retreat rates increased after 1998, exceeding 7% area loss per year, which is higher than the rates reported for the last 30 years in other glaciers of the tropical Andes (Ramírez et al., Reference Ramírez, Melfo, Resler and Llambí2020).
Atmospheric or aeolian deposition (windfall) provide key nutrients such as carbon, nitrate and ammonium to the micro- and mesobiota of the supraglacial zone (Edwards, Reference Edwards1987; Hotaling et al., Reference Hotaling, Hood and Hamilton2017), but this is also a source of light-absorbing particles that increase melting processes (Gilardoni et al., Reference Gilardoni, Di Mauro and Bonasoni2022). Nutrients and meltwater can be transported through interglacial cracks and crevasses to reach the subglacial zone (contact zone between rock and ice substrate), where they combine with small particles produced by rock comminution (Hotaling et al., Reference Hotaling, Hood and Hamilton2017). Greater concentrations of black carbon (sub-micron particle volume increasing from 0.19 to 1.4 μm3/cm3) at high elevations of the Cordillera de Mérida have been linked to biomass burning in Venezuelan savannah, with greater fire activity and higher concentrations observed following El Niño years (Hamburger et al., Reference Hamburger, Matisāns, Tunved, Ström, Calderon and Hoffmann2013).
The ice substrate, proglacial waters and glacier forefield share part of their microbiota because of dispersal through melting water (Hotaling et al., Reference Hotaling, Hood and Hamilton2017) and are here considered compartments of the same ecosystem. Glacier retreat associated with climate change and subsequent soil development triggers multitrophic changes from simple food webs dominated by microbes and heterotrophic organisms to more complex networks as deglaciated areas are colonized by lichen, mosses and vascular plants, marking the transition to a novel ecosystem (Ficetola et al., Reference Ficetola, Marta, Guerrieri, Gobbi, Ambrosini and Fontaneto2021; Khedim et al., Reference Khedim, Cécillon, Poulenard, Barré, Baudin and Marta2021; Llambí et al., Reference Llambí, Melfo, Gámez, Pelayo, Cárdenas and Rojas2021; Rosero et al., Reference Rosero, Crespo-Pérez, Espinosa, Andino, Barragán and Moret2021; Anthelme et al., Reference Anthelme, Carrasquer, Ceballos and Peyre2022).
Prospective samples from englacial and subglacial microbiota were collected from Bolívar and Humboldt peaks (Fig. 1, Plates 1 & 2; Ball et al., Reference Ball, Gómez, Magallanes, Rosales, Melfo and Yarzábal2014; Balcazar et al., Reference Balcazar, Rondón, Rengifo, Ball, Melfo, Gómez and Yarzábal2015; Rondón et al., Reference Rondón, Gómez, Ball, Melfo, Rengifo and Balcázar2016). These include samples taken prior to the complete disappearance of ice from Bolívar peak, but they do not allow for a full assessment of the role of the microbiota in the ecosystem. However, a more detailed study of the microbiota in the glacier forefield at Humboldt peak is underway (B. Huber, pers. comm., 2023). In the forefield area, a diverse array of pioneer lichen and bryophyte species (including many new records for Venezuela and several endemic species) appears to facilitate the establishment of pioneer wind-dispersed and wind-pollinated vascular plants during the first 6 decades of primary succession after glacier ice retreat (Llambí et al., Reference Llambí, Melfo, Gámez, Pelayo, Cárdenas and Rojas2021).
View of the glacier at Humboldt Peak in the Cordillera de Mérida, Venezuela, in August 2022. The photograph was taken from the peak, looking west towards the direction of (1) Bolívar Peak and (2) La Concha Peak. Photo: J.A. González.
Mosaic of satellite images (Sentinel-2 MultiSpectral Instrument, Level 2A; December 2021–March 2022; clouds and shadows removed) over the three peaks in the Cordillera de Mérida, Venezuela: (1) Bolívar, (2) La Concha and (3) Humboldt. The square represents the same 10 × 10 km grid cell shown in Fig. 1.
Methods
We applied the IUCN Red List of Ecosystems protocol (Bland et al., Reference Bland, Keith, Miller, Murray and Rodríguez2017) that requires assessment of five criteria: A, declining distribution; B, restricted distribution and exposure to threats; C, degradation of abiotic environment; D, disruption of biotic processes and E, quantitative risk analysis. Three of the five criteria (A, C and D) include sub-criteria applied to three different timeframes to capture the effects of historical, current and future threats. Criterion B evaluates the current exposure to threats using various spatial metrics, and criterion E is evaluated over 50 and 100 years into the future. The complete assessment is provided in the Supplementary Material.
We captured the main ecosystem processes and the threats of climate change and air pollution in a conceptual ecosystem model (Fig. 2) and identified spatiotemporal indicators for the assessment of each criterion and sub-criterion. We defined a threshold of collapse for each indicator and interpolated or extrapolated time series data to infer rates of change over the relevant timeframes.
Conceptual ecosystem model for the tropical glacier ecosystem of the Cordillera de Mérida, Venezuela. White ovals and boxes: those with a dashed outline represent the main characteristic biota and biotic processes, those with a solid outline the main elements of the abiotic environment and abiotic processes, respectively. Dark grey boxes: threatening processes. Black lines: relationships between elements (arrow: increases; circle: reduces). Light grey boxes: indicator variables used in the assessment and indirectly linked to ecosystem processes or components by grey dashed lines. Some elements and processes are excluded for clarity of visualization.
We defined the collapsed state of the ecosystem as the complete loss of the icy substrate that can sustain a cold-resistant microbiota (zero ice extent for spatial criteria A and B or zero ice mass for criterion E). We used three indicators to assess climatic suitability for glacial function under criterion C: the equilibrium-line altitude (the altitude at which rates of ice accumulation and ablation are equal); the atmospheric freezing level height (the altitude of the 0 °C isotherm); and a suitability threshold based on a correlative bioclimatic suitability model. An increase in equilibrium-line altitude or freezing level height reduces the available area for long-term glacier persistence, with thresholds of collapse for both indicators set to the altitude of the mountain summit (Polissar et al., Reference Polissar, Abbott, Wolfe, Bezada, Rull and Bradley2006; Braun & Bezada, Reference Braun and Bezada2013). We based the threshold of bioclimatic suitability on the optimal classification of the current occurrence of glacier outlines. There is currently no suitable indicator of collapse for criterion D because the microbiota of the supraglacial zone and its relationships with other ecological zones of the glacier are still poorly understood (Ball et al., Reference Ball, Gómez, Magallanes, Rosales, Melfo and Yarzábal2014; Rondón et al., Reference Rondón, Gómez, Ball, Melfo, Rengifo and Balcázar2016). Data collected at Humboldt Peak in 2019 and 2021 (including supraglacial, englacial and subglacial samples as well as soil from the glacial forefront) may shed light on this issue. Preliminary results indicate marked changes in microbiota composition and function and a sizeable role of the microbiota on ecosystem processes in four chronosequence sites deglaciated during 1910–2009 (B. Huber, pers. comm., 2023); similar findings have been reported for a tropical glacier in Ecuador (Díaz et al., Reference Díaz, Monfort-Lanzas, Quiroz-Moreno, Rivadeneira, Castillejo and Arnau2023).
We used available measures of glacier ice extent to calculate best estimates of past decline and to extrapolate rates of future declines in ice extent (sub-criteria A1, A2b and A3). We obtained measurements and their standard errors from published, ground-validated reconstructions based on cartographic data from 1910, aerial photographs from 1952 and 1998, satellite images from 2009, 2015 and 2016 and field measurements using GPS and drone imagery from 2019 (Ramírez et al., Reference Ramírez, Melfo, Resler and Llambí2020). Indirect and less accurate evidence of the extent of glacier ice prior to 1910 is available from a few other sources and is provided for context (Polissar et al., Reference Polissar, Abbott, Wolfe, Bezada, Rull and Bradley2006; Jomelli et al., Reference Jomelli, Favier, Rabatel, Brunstein, Hoffmann and Francou2009; Braun & Bezada, Reference Braun and Bezada2013). We calculated past declines based on pairs of measurements and considering error propagation. For future declines we calculated first proportional rates of decline (as percentage area loss per year) and then used this rate to project the expected magnitude of decline over a 50-year period (1948–1998). Given the evidence suggesting accelerating rates of decline, we compared the results of using different timeframes for calculation of proportional rates of decline.
We used two spatial measures of exposure to threats (criterion B): the extent of occurrence was calculated as the area of the convex hull around the glacier outlines from the RGI 6.0 database (Randolph Glacier Inventory Consortium, 2017), and the area of occupancy was calculated as the minimum number of 10 × 10 km cells that include all of these occurrences.
For criterion C1 we used a time series of freezing level height in metres calculated from climate reanalysis data for 1948–2011 (Braun & Bezada, Reference Braun and Bezada2013) and fitted a local polynomial regression (loess with gaussian distribution, span = 0.75 and degree = 2, equivalent number of parameters = 4.35) to smooth the temporal trend. We calculated relative severity of degradation (RS) using the formula: RS = OD/MD, where OD (observed decline) is the difference between the initial and final values and MD (maximum decline) is the difference between the initial value and the collapse threshold (Keith et al., Reference Keith, Rodríguez, Rodríguez-Clark, Nicholson, Aapala and Alonso2013). We used the initial and final values of the smoothed trend in freezing level height and two collapse thresholds of 4,920 and 4,970 m, which span the range of plausible estimates. We took estimates of historical changes in equilibrium-line altitude from Polissar et al. (Reference Polissar, Abbott, Wolfe, Bezada, Rull and Bradley2006), and these are based on reconstructed palaeo-glacier topography near Mucuñuque Peak (located in nearby Sierra de Santo Domingo; Fig. 1) and cumulative elevation profiles of modern glaciers at Bolívar Peak. For the calculation of relative severity of degradation we used change in equilibrium-line altitude as a direct measurement of observed decline and the maximum possible change as maximum decline (sub-criterion C3).
We also used local results of correlative and mechanistic global models to project trends in climatic suitability (sub-criterion C2a) and the probability of persistence of tropical glaciers into the future (criterion E). These models incorporated different global circulation models and scenarios of projected greenhouse gas emissions based on alternative socio-economic developments and climate policies (shared socio-economic pathways; Karger et al., Reference Karger, Conrad, Böhner, Kawohl, Kreft and Soria-Auza2017).
The correlative model of environmental suitability focused exclusively on tropical regions and compares the bioclimatic conditions of high-elevation areas with and without glaciers using the Gradient Boosting Machine algorithm (Ferrer-Paris & Keith, Reference Ferrer-Paris and Keith2024). We selected occurrence records using stratified random sampling from glacier outlines in tropical areas (from the tropical Andes, Mexico, Africa and Asia) and non-glaciated areas with a minimum of 3,500 m elevation and a maximum distance of 25 km from the nearest glacier outlines. We withheld occurrence records from the Cordillera de Mérida for final evaluation of prediction performance of the model; we divided all other records into calibration (80%) and test partitions (20%) and used them for tuning of the model parameters (number of trees, interaction depth, shrinkage and minimum number of observations per node) and final model fitting using cross-validation. We compared the predicted probability or suitability index with the known occurrences in the test partition and in the Cordillera de Mérida to calculate cut-off values of equal sensitivity and specificity (balancing errors of omission and commission) and maximum accuracy (minimizing misclassification errors). We then used the model to predict future suitability under 15 combinations of global circulation models and socio-economic pathways, assuming the alternative cut-off values represented thresholds of ecosystem collapse (sub-criterion C2a).
The quantitative analysis of collapse for criterion E was based on a global hybrid model (Rounce et al., Reference Rounce, Hock, Maussion, Hugonnet, Kochtitzky and Huss2023) that combined a mass balance module and a glacier dynamics module to model all glaciers in the world independently using globally available datasets of glacier outlines (Randolph Glacier Inventory Consortium, 2017), glacier-wide geodetic mass balance data and regional ice volume estimates for calibration (Farinotti et al., Reference Farinotti, Huss, Fürst, Landmann, Machguth, Maussion and Pandit2019; Hugonnet et al., Reference Hugonnet, McNabb, Berthier, Menounos, Nuth and Girod2021). In the case of small regions with no direct measurements (such as the Cordillera de Mérida) the model used initial estimates of ice volume based on digital elevation models and most likely overestimated initial mass, making the projections conservative (Rounce et al., Reference Rounce, Hock, Maussion, Hugonnet, Kochtitzky and Huss2023). The model was run with monthly and annual time steps from 2000 to 2100 for various ensembles of global circulation models and socio-economic pathways. We used the ice mass projections of the model for the glacier outlines of the Cordillera de Mérida (Rounce et al., Reference Rounce, Hock and Maussion2022) to estimate the year of collapse (first year when ice mass reaches zero) for 48 combinations of global circulation models and socio-economic pathways. We used the empirical cumulative distribution function of the year of collapse to calculate the proportion of models indicating collapse for each year and estimate probability of collapse.
We performed all statistical analysis using R 4.3.1 (R Core Team, 2023). For the suitability model we used functions in the R package caret (Kuhn, Reference Kuhn2008).
Results
There are no direct field measurements of decline in extent of glacier ice for the last 50 years, but cartographic estimates of glacier extent suggest a decline of at least 89% over 21 years (proportional rate of decline = 5.173 ± SE 0.130%/year) and nearly 99% if we consider the last 67 years (proportional rate of decline = 10.220 ± SE 0.473%/year). These rates of decline would translate to future declines of at least 94% by the year 2048 (Table 1).
Estimated magnitude of past and future tropical glacier declines in the Cordillera de Mérida, Venezuela (Fig. 1), based on previous measurements (Ramírez et al., Reference Ramírez, Melfo, Resler and Llambí2020) and estimated proportional rates of decline.
Given the rapid decline and disappearance of the glacier at Bolívar Peak during 1998–2017 and that the size of the Humboldt glacier at the time of the last measurement (0.045 km2 in 2019) was comparable to that of the Bolívar glacier in 1998, it is reasonable to expect that the Humboldt glacier will disappear completely over the next 20 years. The area of glacier ice in Humboldt Peak was last measured at 0.046 ± 0.004 km2 in July 2019 (Ramírez et al., Reference Ramírez, Melfo, Resler and Llambí2020) and was estimated to be c. 0.01 km2 in August 2022 (Plate 1; N. Ramírez, pers. comm., 2023). This last field observation would confirm the recent acceleration of the rate of decline of the Humboldt glacier and increases the likelihood that it will disappear sooner than predicted, probably within the next 5 years.
Reconstruction of the historical evolution of glaciers in South America suggests that the most recent maximum glacial extent in Venezuela occurred c. 1730 (during the Little Ice Age), and glaciers retreated continuously in the following centuries, with only minor readvances c. 1760, 1820 and 1880 (Polissar et al., Reference Polissar, Abbott, Wolfe, Bezada, Rull and Bradley2006; Jomelli et al., Reference Jomelli, Favier, Rabatel, Brunstein, Hoffmann and Francou2009). Thus, we can assume that the observed decline between 1910 and the present represents a lower bound of the total decline between 1750 and the present (Table 1).
The extent of occurrence for the Humboldt glacier and surrounding areas is 0.892 km2. If the proglacial waters and glacier forefields of recently collapsed glaciers are included, the extent of occurrence is at most 5.957 km2. The main occurrences of known glaciers (extant and recently collapsed) occupy a single 10 × 10 km cell (Fig. 1, Plate 2). This cell probably contains all occurrences of associated habitats with connected microbiota (proglacial waters and glacier forefields).
The relationship between recent temporal changes in climatic conditions and the disappearance of glacier ice is best illustrated with the case of La Concha Peak. The estimated maximum elevation of ice was 4,840 m in the year 1952, and it disappeared before 1998. The mean freezing level height was almost 10 m below the maximum elevation in 1948 and rose to more than 83 m above the maximum elevation in 2010, and the smoothed freezing level height mean surpassed the 4,840 m threshold during 1972–1981 (Fig. 3). Using the thresholds of maximum elevation of glacier ice in the two other peaks, we calculated the relative severity of change for this indicator to be 67.2 ± SE 18.5% (Bolívar Peak) and 104.4 ± 28.5% (Humboldt Peak). The large standard error is because of the wide variation in the time series.
Local polynomial regression of freeze level height in the Cordillera de Mérida, Venezuela, for the period 1948–2011, with 95% confidence intervals (based on data from Braun & Bezada, Reference Braun and Bezada2013). Horizontal lines represent the last recorded elevation of glacier ice at the different peaks, with years given in parentheses.
The best-performing bioclimatic suitability model had high predictive performance during model training (area under the curve = 0.967 ± SE 0.008, sensitivity = 0.621 ± 0.058 and specificity = 0.981 ± 0.004, with Gradient Boosting Machine model parameters: 200 trees, interaction depth = 5, shrinkage = 0.1 and a minimum of 12 observations per node) and when evaluating model predictions under current climatic conditions in the Cordillera de Mérida (area under the curve = 0.990; sensitivity = 0.667 and specificity = 0.996; full model details in Supplementary Material). The best estimate of mean relative severity for this indicator across all combinations of five global circulation models, three pathways, two timeframes and two collapse thresholds is 97%, with a 90% confidence interval of 63–100%.
There are three plausible values of difference in equilibrium-line altitude: the historical reconstruction of equilibrium-line altitude between the maximum glacier extent (pre-1820) and 1972 suggests a value of −300 m, whereas comparison with the historically and recently collapsed Mucuñuque and Bolívar Peaks suggests values between −400 and −500 m (Polissar et al., Reference Polissar, Abbott, Wolfe, Bezada, Rull and Bradley2006). Using a value of −550 m as the maximum decline, we calculated the relative severity of degradation for this indicator and timeframe to be 54–90%, with a median value of 72%.
The hybrid model predicted a high probability of collapse in the future. Focusing on the 50-year period 2020–2070, we found that 79% of the models end in collapse (Fig. 4). Uncertainty in mass estimates (mean absolute deviation) does not have a significant effect on the estimated year of collapse. Considering each scenario separately, the proportion of models that predict collapse by 2070 is higher than 50% in all cases, except for the sustainable development scenario (socio-economic pathway 1-2.6).
Empirical cumulative distribution function of year of collapse from a global hybrid model of mass balance and glacier dynamics based on future climate predictions from 48 combinations of global circulation models and shared socio-economic pathways (SSPs; thick solid line) and disaggregated by the four different pathway (labelled lines): SSP1-2.6 is a sustainable development scenario, SSP2-4.5 is intermediate, SSP3-7.0 prioritizes national development and SSP5-8.5 is fossil-fuelled development.
These quantitative results support the assessment of the risk of collapse using the IUCN Red List of Ecosystems categories and criteria (Table 2).
Summary of the IUCN Red List of Ecosystems assessment of the tropical glacier ecosystem of the Cordillera de Mérida, Venezuela. The first column indicates the sub-criteria assessed under each criterion. A more detailed report is available in the Supplementary Material. The categories of ecosystem risk are: Collapsed (CO), Critically Endangered (CR), Endangered (EN), Vulnerable (VU), Near Threatened (NT), Least Concern (LC) and Data Deficient (DD). A threat defined location is ‘a geographically or ecologically distinct area in which a single threatening event can rapidly affect all occurrences of an ecosystem type' (Bland et al., Reference Bland, Keith, Miller, Murray and Rodríguez2017).
Discussion
Based on the analyses presented here, the current Red List status of Cordillera de Mérida glacier is Critically Endangered, with plausible evidence that it may already be Collapsed: CR (CR-CO). Almost all of the criteria and sub-criteria assessed indicate the CR status, implying similar magnitudes of historical, present and future threats (Table 2). For the Cordillera de Mérida the evidence from recent field visits suggests that its glacier ice will probably disappear before 2048 (sub-criterion A2b). Other lines of evidence suggest similar conclusions, although with more uncertainty because of the use of indirect indicators and global instead of local data. During the last 100 years, the freezing level height and equilibrium-line altitude have increased close to the point of collapse, and permanent ice has completely disappeared from all previously glaciated peaks except Humboldt (sub-criteria A1, A3 and C1). The current distribution of glaciated areas is much reduced and exposed (sub-criteria B1 and B2), is undergoing rapid decline at a rate similar to the recently disappeared Bolívar glacier and is unlikely to recover mass before its complete collapse (sub-criterion C2a and criterion E).
This is the first Red List of Ecosystems assessment of a tropical glacier ecosystem. Cryogenic ecosystems are often under-represented in national assessments that focus on vegetation units, but despite differences in processes and functions they can be assessed using the same criteria and assigned to comparable categories of risk (Keith et al., Reference Keith, Rodríguez, Rodríguez-Clark, Nicholson, Aapala and Alonso2013). Risk assessments and long-term monitoring of all ecosystems in tropical mountains can provide valuable information for comprehensive conservation planning and ecosystem management in these complex landscapes (Cuesta et al., Reference Cuesta, Llambí, Huggel, Drenkhan, Gosling and Muriel2019; Llambí et al., Reference Llambí, Becerra, Peralvo, Avella, Baruffol and Díaz2020).
This quantitative assessment of several symptoms of collapse illustrates how the Red List of Ecosystems protocol can be applied to other tropical glaciers (and cryogenic ecosystems outside tropical areas) by combining all available evidence within a common framework. We anticipate that similarly isolated and exposed tropical glaciers in Indonesia, East Africa, Mexico and some Andean regions (e.g. Santa Isabel in Colombia, Carihuairazo in Ecuador) will also be at extreme risk of collapse, and the larger regional units in the tropical Andes may experience widespread decline and degradation (Ceballos et al., Reference Ceballos, Euscátegui, Ramírez, Cañon, Huggel, Haeberli and Machguth2006; Rabatel et al., Reference Rabatel, Ceballos, Micheletti, Jordan, Braitmeier and González2018; Veettil & Kamp, Reference Veettil and Kamp2019; Ferrer-Paris & Keith, Reference Ferrer-Paris and Keith2024). The Red List of Ecosystems protocol provides a valuable framework for consistent and comparable assessments across all continents that synthesizes local and regional data from geology, climatology, palaeontology/palynology, microbiology, entomology and soil and vegetation ecology to describe ecosystems and diagnose the threats that they face.
The Red List of Ecosystems assessment explores different pathways towards collapse, and the associated uncertainty of methods and projections is consolidated within the categories of risk and their plausible bounds. Here we fully exploited this versatility, assessing spatial and functional symptoms, undertaking quantitative risk analysis and considering alternative interpretations of collapse thresholds. Our study shows that the multiple lines of evidence offer a broad picture of threatening processes and pathways towards collapse and that the assessment protocol handles several sources of uncertainty (Keith et al., Reference Keith, Rodríguez, Rodríguez-Clark, Nicholson, Aapala and Alonso2013; Bland et al., Reference Bland, Keith, Miller, Murray and Rodríguez2017). Sufficient data are available for the comprehensive risk assessment of tropical glacier ecosystems in different regions of the world, and global datasets and analyses can be used as informative indicators of the climate and cryosphere dynamics that have not been studied locally, with the caveat of the limited ground validation and potential bias of these sources (Sagredo & Lowell, Reference Sagredo and Lowell2012; Rounce et al., Reference Rounce, Hock, Maussion, Hugonnet, Kochtitzky and Huss2023).
Because of the rapid decline of the icy substrate in recent decades, the small remnants of ice in Venezuela and on other tropical mountains are often considered static or extinct glaciers because dynamic processes of ice accumulation no longer operate and ice loss is accelerated by scale and edge effects (Ceballos et al., Reference Ceballos, Euscátegui, Ramírez, Cañon, Huggel, Haeberli and Machguth2006; Braun & Bezada, Reference Braun and Bezada2013; Rabatel et al., Reference Rabatel, Ceballos, Micheletti, Jordan, Braitmeier and González2018; Ramírez et al., Reference Ramírez, Melfo, Resler and Llambí2020). Despite rapid contraction, the little remaining ice substrate could sustain much of its original microbiota (Ball et al., Reference Ball, Gómez, Magallanes, Rosales, Melfo and Yarzábal2014; Balcazar et al., Reference Balcazar, Rondón, Rengifo, Ball, Melfo, Gómez and Yarzábal2015; Rondón et al., Reference Rondón, Gómez, Ball, Melfo, Rengifo and Balcázar2016). Studies of the glacier forefield at Humboldt Peak provide insights into the postglacial chronosequence: pioneer lichens, mosses and vascular vegetation are already present 10 years after the retreat of the glacier. However, vascular plants only increase slowly in cover, species and functional diversity during the first 100 years of primary succession, and soil properties (e.g. soil organic matter, total nitrogen and exchangeable bases) change significantly after 21 and 60 years (Llambí et al., Reference Llambí, Melfo, Gámez, Pelayo, Cárdenas and Rojas2021).
Although collapse seems imminent and is probably unavoidable, we can still learn much from inventories of native biota and monitoring of ecosystem transitions. Multiple changes in climate, soil and ecological communities operate at different temporal scales in high Andean summits and the recently deglaciated lands (Ficetola et al., Reference Ficetola, Marta, Guerrieri, Gobbi, Ambrosini and Fontaneto2021; Khedim et al., Reference Khedim, Cécillon, Poulenard, Barré, Baudin and Marta2021; Llambí et al., Reference Llambí, Melfo, Gámez, Pelayo, Cárdenas and Rojas2021; Anthelme et al., Reference Anthelme, Carrasquer, Ceballos and Peyre2022; Cuesta et al., Reference Cuesta, Carilla, Llambí, Muriel, Lencinas and Meneses2023). A key research question is whether some biotic components of glacial ecosystems will be lost or whether they will persist in cold, ice-free refuges after glacial collapse (Stibal et al., Reference Stibal, Bradley, Edwards, Hotaling, Zawierucha and Rosvold2020). However, we have limited time to study the remainder of this glacier microbiota in the novel ecosystems that form after the complete retreat of ice.
Author contributions
Study design: JRF-P, DAK; data review: JRF-P, LDL, AM; data analysis: JRF-P; writing: all authors.
Acknowledgements
We thank Nerio Ramírez for his invaluable work documenting glacier retreat and Barbara Huber for kindly sharing preliminary results on glacier microbiota in the Venezuelan Andes. The National Geographic Society provided financial support for LDL and AM (Grant No. NGS-55170R-19). We also thank the Mucumbarila cable car system and INPARQUES for their support in accessing the study sites near Humboldt Peak.
Conflicts of interest
None.
Ethical standards
This research abided by the Oryx guidelines on ethical standards.
Data availability
The data and code that support the findings of this study are available at osf.io/y3279.
Open access
