We present a continentwide study of 600 glaciers located on and near 37 ice-clad volcanoes in South America. Results demonstrate glacier sensitivity to volcanic heat. We distinguished between “volcanic glaciers” (≤1 km from volcanic centers; n = 74), and “proximal glaciers” (1–15 km; n = 526) and calculated their equilibrium line altitudes (ELAs). For each ice-clad volcano, we compared the ELAs of its volcanic glaciers to those of its proximal glaciers, which showed that the ELAs of the former are higher than the ELAs of the latter. ΔELAmean, defined as the offset between the mean ELA of the volcanic glaciers compared with that of the proximal glaciers, was calculated for each ice-clad volcano. ΔELAmean was positive for 92% of the 37 volcanoes, and a quantitative relationship between ΔELAmean and volcanic thermal anomaly was established. Results highlight the impact of volcanic heat on glacier elevation; emphasize the need to exclude glaciers on, or near, volcanoes from glacier-climate investigations; and demonstrate the first-order potential for glaciers as “volcanic thermometers.” Volcanic-glacier monitoring could contribute to our understanding of magmatic and thermal activity, with changes in glacier geometries potentially reflecting long-term fluctuations in volcanic heat and unrest.

Volcanic eruptions are natural socioeconomic hazards with devastating consequences, including the displacement of communities; damage to businesses and infrastructure; disruption of air traffic; and loss of human life (Loughlin et al., 2015). A major challenge in the management of such hazards is the identification and monitoring of precursors to forthcoming volcanic eruptions. The measurement of thermal anomalies is one such technique, with some volcanoes exhibiting signs of thermal unrest for several years prior to an eruptive event (Reath et al., 2019; Girona et al., 2021). Although thermal anomalies may be detected using remote-sensing methods, glaciers on volcanoes may mask the thermal anomalies, impacting monitoring efforts.

While glaciers on volcanoes are considered a major risk (Tuffen, 2010; Edwards et al., 2020) and a hindrance to obtaining accurate temperature measurements, they are likely to be affected by volcanic heat (Barr et al., 2018). If the impact of volcanic heat on glaciers can be demonstrated and quantified, this could help to improve magmatic system dynamics models and be used as a novel tool to monitor long-term changes in the thermal state of ice-covered volcanoes, which may otherwise be obscured from most conventional remote-monitoring systems. Mapping of volcanic glaciers using geospatial tools can be used to analyze the interplay between glacier geometries, glacier equilibrium line altitudes (ELAs), and volcanic activity (Rivera et al., 2006; Rivera and Bown, 2013; Reinthaler et al., 2019). This is the first large continental-scale analysis and quantitative assessment of the potential link between glacier geometries and volcanic heat.

This study focused on the Andes (Fig. 1), where many volcanoes have glaciers within 1 km (volcanic glaciers) and between 1 and 15 km (proximal glaciers) from their center, and, importantly, maximum thermal anomaly measurements are available for some of these volcanoes (Reath et al., 2019). We assumed that a volcanic glacier (located on a volcano) will likely experience a basal melt rate significantly higher than a proximal glacier (not on a volcano; Fig. 2). To assess this effect, we could look at metrics such as the minimum, median, or average glacier elevation for volcanic versus proximal glaciers. However, these values can be impacted significantly by the local topography, so we calculated a different metric, the glacier ELA, using the area-altitude balance ratio (AABR) method. This metric accounts for glacier geometry via the hypsometry (i.e., the distribution of the surface area with altitude) and recognizes that the surface accumulation and ablation gradients differ (Rea, 2009). The ELA is the point on the glacier where the surface mass balance, measured over 1 yr, is zero; i.e., accumulation (snowfall) equals ablation (snowmelt/sublimation). The ELA can be measured in the field via repeated (time-consuming and logistically challenging) observations. Calculating ELAs using the AABR method is a good proxy for measured ELA (Oien et al., 2021), provided glaciers are clean (no debris cover) and terrestrially terminating (not in water), and basal melt contributes a negligible component of the overall mass balance. For volcanic glaciers, where the basal melt rate may be significant in the overall mass balance, the calculated ELAs (hereafter ELAs) are not a good proxy for the measured ELA, but they nonetheless represent an ideal metric with which to characterize the elevation of the glaciers, taking into account their hypsometry. By calculating the AABR ELAs for volcanic versus proximal glaciers, we attempted to identify the potential impact of volcanic heat on glacier geometries; i.e., we expected the ELA for the volcanic glaciers to be higher than the ELA for the proximal glaciers. We also made the reasonable assumptions that the ELAs of the proximal glaciers are comparable within a restricted geographic area (i.e., a 15 km radius), as they will experience a similar climate (Sagredo et al., 2014), and that the ELA of the volcanic glaciers is a function of both climate and the volcanic heat that drives additional ice loss (Jóhannesson et al., 2020). In this continent-scale study, we combined data from worldwide glacier (Randolph Glacier Inventory [RGI] 6.0) and volcano (Glopal Volcanism Program [GVP], 2013) inventories to identify 37 Holocene Andean volcanoes that host glaciers on (volcanic) and near them (proximal). Water-terminating and debris-covered glaciers were excluded, as were glaciers <0.1 km2, to limit the effect of complex glacier dynamics and niche microclimates. We calculated ELAs for 74 volcanic glaciers and 526 proximal glaciers (Table S1 in the Supplemental Material1) distributed latitudinally from 5°N to 41°S along the Andes (Fig. 1). For each selected volcano, we assessed how ELA varied with distance from the volcano and calculated the difference in mean ELA between the volcanic glaciers and the proximal glaciers (i.e., ΔELAmean; Fig. 2). Glacier ΔELAmean was then compared with Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER)–based volcano temperature anomalies (Reath et al., 2019), acquired by the National Aeronautics and Space Administration (NASA) Terra satellite, for 13 of the 37 ice-clad volcanoes (Table S2). We also compared ΔELAmean with climate data from WorldClimVersion 2 (Fick and Hijmans, 2017). Full methodological details are provided in the Supplemental Material.

ΔELAmean: Assessing the Impact of Volcanic Heat on Glaciers

Results highlight that 469 (89%) of the 526 proximal glaciers are characterized by an ELA lower than the mean ELA of the nearby volcanic glaciers. For 50% of these, a statistically significant correlation (R2 > 0.50 at p < 0.05) was established between glacier ELA and distance from the volcano center (Fig. 2). For example, proximal glacier ELAs gradually decreased away from the volcano by as much as 655 m for Copahue volcano (Chile-Argentina), from an ELA of 2807 m on the volcanic glacier to 2152 m on a proximal glacier located 8.46 km away (R2 = 0.87; Figs. 2C and 2D). Weaker relationships between glacier ELA and distance from the volcano (Table S1) were rarely found and might be due to local microclimate (e.g., aspect/orientation of slope [Evans, 2006] or shading) and/or volcanogenic (e.g., offset magma reservoirs; Lerner et al., 2020) factors.

We used the ΔELAmean to investigate how glaciers are affected by the volcano heat, i.e., measured thermal anomalies. For 92% of the ice-clad volcanoes (Fig. 2A; Table S1), the ΔELAmean was positive (i.e., the mean ELA of volcanic glaciers was higher than that for the proximal glaciers), with a mean ΔELAmean of 229 m and a median of 187 m. Given the relatively short distances considered (<15 km) and the large number (n = 600) of glaciers analyzed (comprising different slope aspects, etc.), local climate variations cannot be invoked to explain our results. Instead, we take this as a strong indication that the offset in ELA between proximal and volcanic glaciers is controlled primarily by the volcanic heat source.

How Does Volcanic Heat Affect Glacier Elevation?

To demonstrate a quantitative, empirical relationship between glaciers and volcanoes, continuous volcanic heat measurements covering a temporal interval longer than the glacier response times would have been ideal, but these are not available. Instead, Reath et al. (2019) provided direct observations of volcano maximum thermal anomalies, recorded between 2000 and 2018, which were obtained from Terra satellite data for 88 volcanoes in Central and South America, including some of the volcanoes analyzed here. For 13 (Table S2) of the original 37 Holocene ice-clad volcanoes, it was possible to analyze the correlation between the mean volcano maximum thermal anomaly (mean δTmax) and the ΔELAmean to establish a first-order quantitative assessment of glacier elevation sensitivity to volcanic thermal state. Our results demonstrated a strong, positive relationship (R2 = 0.72, p < 0.001) between mean δTmax and ΔELAmean (Fig. 3), with δTmax = 0.22 × ΔELAmean – 25.90 and a slope uncertainty of ±0.04 (at 95% confidence level). For example, on the quiescent Falso Azufre (Chile-Argentina), Parinacota (Chile-Bolivia), and Pular (Chile) volcanoes, the ΔELAmean was relatively low (133–156 m), as was mean δTmax (3.7–5.7 °C). However, on the presently active Copahue (Chile-Argentina) and Villarrica (Chile) volcanoes, both the ΔELAmean and mean δTmax were much higher (Fig. 3), with values ranging 252–270 m and 35.9–44.7 °C, respectively. This provides confidence that volcanic-glacier geometries respond to increased volcanic heat.

The response time of a glacier to a surface mass balance perturbation related to climate can be approximated as a function of the ice thickness and ablation rate at the terminus (Paterson, 1994), such that, for a 50–100-m-thick glacier with an ablation rate at the terminus of 5 m yr–1 (reasonable values for the glaciers under consideration here), the response time of the glacier is ~10–20 yr. If the ablation rate is increased to 10 m yr–1, due to enhanced basal melt, it would likely reduce the response time to 5–10 yr.

Although the exact response time of glaciers to volcanic-induced enhanced basal melt is not known, the strong relationship between mean δTmax and ΔELAmean is evidence that volcanic heat does enhance glacier basal melt, resulting in volcanic glaciers located at higher elevations, with concomitantly higher ELAs than for their proximal glaciers. This is particularly encouraging given that we used the maximum thermal anomaly and not a continuous measurement. These results indicate that the ΔELAmean could be used as a first-order approximation, over reasonable time scales, e.g., 5–10 yr (Girona et al., 2021), to identify changes in volcanic heat output and contribute to monitoring ice-clad volcanoes, particularly where the presence of glacial ice may otherwise obstruct or complicate the volcano thermal signature.

Variation of ΔELAmean with Climatic Region

The Andes are characterized by three distinct megaclimatic zones (Garreaud et al., 2009; Sagredo and Lowell, 2012), corresponding well with the Northern, Central, and Southern volcanic zones (Fig. 1; Tilling, 2009), and the climate influence on glacier ELAs along the Andes is well documented (Vuille et al., 2008; Rabatel et al., 2013; Braun et al., 2019). In principle, it is possible that climate also affects ΔELAmean and could limit its applicability as a proxy for volcano thermal activity. However, the correlation between ΔELAmean and total annual precipitation (Ptot) is very weak (R2 = 0.05, p = 0.164), as is that for mean annual air temperature (Tmean; R2 = 0.08, p = 0.06; Fig. 4).

Given the lack of correlation between ΔELAmean and climate, we concluded that, while the ELAs of volcanic and proximal glaciers are, respectively, in part and fully controlled by climate, the ΔELAmean is little impacted by variations therein. For example, the substantial decrease in volcanic-glacier ELAs (from an average of 6021 m to 2983 m) when migrating from the dry subtropical Andes of Peru and northern Chile (Ptot of 115–757 mm, Tmean of −9.48 °C to −2.38 °C) in the Central volcanic zone (15.52°S–27.20°S) to the warmer and wetter semiarid regions along the border of Chile and Argentina (Ptot of 495–1652 mm, Tmean of −7.98° to 3.99 °C) in the Southern volcanic zone (34.16°S–40.97°S) is likely due to climate. Significantly, the ΔELAmean remains consistent throughout these volcanic zones (an average of 209 m from Peru and northern Chile to an average of 182 m across the Chile-Argentine border).

Uncertainties

ELAs have a computational accuracy of 5 m (Pellitero et al., 2015). An ~5% gross geometry error for the Southern Andes (region 17) for RGI 6.0 glacier outlines, due to the erroneous inclusion of seasonal glacier-peripheral snow and transient ice, was reported by Pfeffer et al. (2014). Our exclusion of glaciers <0.1 km2 will have reduced the likelihood of including some erroneously mapped snow patches and seasonal ice cover. While checking the mapping of all 600 glaciers would have been unfeasible, we remapped outlines for the 13 volcanoes with a record of δTmax, also to align the temporal observation of δTmax with that of the glacier extent (Supplemental Material).

The algorithm for calculating ASTER-derived temperatures is accurate to ± 1–2 °C (Abrams, 2000). Magma vent and tectonic structures could be complex (García et al., 2019), and hence GVP volcano center points can underestimate the extent and location of volcanic activity (and thus glacio-volcanic interactions). However, a volcano by volcano (field-based) analysis of the magmatic geometry and geothermal heat flux was beyond the scope of this project.

In this study, we analyzed 74 volcanic glaciers and 526 proximal glaciers. For most locations, the ELA of proximal glaciers was lower than that of the volcanic glaciers, with a tendency for proximal glacier ELAs to decrease with distance from the volcanic center. For 92% of the 37 ice-clad volcanoes, the difference in mean ELA between the volcanic and proximal glaciers (i.e., the ΔELAmean) was positive. For a subset (13) of these 37 volcanoes, a strong, positive correlation was identified between ΔELAmean and observed volcano maximum temperature anomalies (i.e., mean δTmax).

These results indicate that volcanic heat alters glacier geometries, which we have highlighted using calculated ELAs, through what is assumed to be enhanced basal melting. Volcanic glaciers tend to be confined to higher elevations and so have higher ELAs relative to their proximal glacier neighbors. For this reason, glaciers located on, or near, Holocene volcanoes should be excluded from studies assessing the impact of recent or ongoing climate forcing on glacier dynamics and elevation. Conversely, and importantly, this study shows that ΔELAmean can be used as a first-order approximation for volcanic thermal anomalies; i.e., high ΔELAmean means high volcanic heat. Monitoring ΔELAmean for glacio-volcanic complexes may help to identify changes in the thermal state of a volcano and could provide a long-term (e.g., 5–10 yr) indication of increased/renewed activity that can be used to improve our understanding of magma dynamics, identify volcanoes of concern, and help assess future periods of volcanic unrest.

1Supplemental Material. Further methodological details and tabular data of analyzed volcanoes and glaciers. Please visit https://doi.org/10.1130/GEOL.S.24008367 to access the supplemental material, and contact editing@geosociety.org with any questions.

This project was supported by the Natural Environment Research Council (NERC) Global Partnerships Seedcorn grant NE/W003724/1 and the Leverhulme Trust Research Project RPG-2019–093. We thank Kevin A. Reath for support in analyzing volcanic thermal anomalies. A. Gomez-Patron and M.E. Pritchard were partly supported by the National Aeronautics and Space Administration Science Mission Directorate Earth Surface and Interior grant 80NSSC21K0842. J. Jaszewski was supported by the NERC Quadrat Doctoral Training Partnership. We express our gratitude to John Smellie and an anonymous reviewer for fruitful feedback, which greatly improved the manuscript.

Gold Open Access: This paper is published under the terms of the CC-BY license.