Icy thermometers: Quantifying the impact of volcanic heat on glacier elevation

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 km², 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 Material¹) 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., ΔELA_mean; Fig. 2). Glacier ΔELA_mean 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 ΔELA_mean with climate data from WorldClimVersion 2 (Fick and Hijmans, 2017). Full methodological details are provided in the Supplemental Material.

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 (R² > 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 (R² = 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 ΔELA_mean 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 ΔELA_mean was positive (i.e., the mean ELA of volcanic glaciers was higher than that for the proximal glaciers), with a mean ΔELA_mean 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.

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 δT_max) and the ΔELA_mean to establish a first-order quantitative assessment of glacier elevation sensitivity to volcanic thermal state. Our results demonstrated a strong, positive relationship (R² = 0.72, p < 0.001) between mean δT_max and ΔELA_mean (Fig. 3), with δT_max = 0.22 × ΔELA_mean – 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 ΔELA_mean was relatively low (133–156 m), as was mean δT_max (3.7–5.7 °C). However, on the presently active Copahue (Chile-Argentina) and Villarrica (Chile) volcanoes, both the ΔELA_mean and mean δT_max 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 δT_max and ΔELA_mean 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 ΔELA_mean 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.

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 ΔELA_mean and could limit its applicability as a proxy for volcano thermal activity. However, the correlation between ΔELA_mean and total annual precipitation (P_tot) is very weak (R² = 0.05, p = 0.164), as is that for mean annual air temperature (T_mean; R² = 0.08, p = 0.06; Fig. 4).

Given the lack of correlation between ΔELA_mean and climate, we concluded that, while the ELAs of volcanic and proximal glaciers are, respectively, in part and fully controlled by climate, the ΔELA_mean 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 (P_tot of 115–757 mm, T_mean 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 (P_tot of 495–1652 mm, T_mean 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 ΔELA_mean 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).

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 km² 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 δT_max, also to align the temporal observation of δT_max 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 ΔELA_mean) was positive. For a subset (13) of these 37 volcanoes, a strong, positive correlation was identified between ΔELA_mean and observed volcano maximum temperature anomalies (i.e., mean δT_max).

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 ΔELA_mean can be used as a first-order approximation for volcanic thermal anomalies; i.e., high ΔELA_mean means high volcanic heat. Monitoring ΔELA_mean 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.

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.