Dwindling Relevance of Large Volcanic Eruptions for Global Glacier Changes in the Anthropocene

Large volcanic eruptions impact climate through the injection of ash and sulfur‐containing gases into the atmosphere. While the ash particles fall out rapidly, the gases are converted to sulfate aerosols that reflect solar radiation in the stratosphere and cause a lowering of the global mean surface temperature. Earlier studies have suggested that major volcanic eruptions resulted in positive mass balances and advances of glaciers. Here, we perform a multivariate analysis to decompose global glacier mass changes from 1961 to 2005 into components associated with anthropogenic influences, volcanic and solar activities, and the El Niño‐Southern Oscillation. We find that the global glacier mass loss was mainly driven by the anthropogenic forcing, interrupted by a few years of intermittent mass gains following large volcanic eruptions. The relative impact of volcanic eruptions has dwindled due to strongly increasing greenhouse gas concentrations since the mid‐20th century.

In this study, we quantify the impact of stratospheric aerosols from large volcanic eruptions on glacier changes at global scale. We follow an approach by Lean and Rind (2008) to perform multiple linear regression analyses to decompose a data set of observed global glacier mass changes (Zemp et al., 2019b) from the second half of the 20th century into components associated with anthropogenic influences, volcanic and solar activities, and El Niño-Southern Oscillation (ENSO); we also extend the analysis back to 1850. The empirical model allows us to quantify the impact of large volcanic eruptions as compared to other forcings and to discuss related changes since Little Ice Age (LIA) maxima. Being aware of the complexity of atmosphere-ice interactions, we deliberately chose this simplified approach to establish a zero-order estimate of the magnitude of the effects on the global scale. In a controlled model run, we repeat the 1991 eruption of Pinatubo to demonstrate its relative impact as compared to other forcings over time. Finally, we discuss the limitations of our approach and if glacier front variations can be used to extend the analysis further back in time.

Glacier Distribution and Fluctuations
Glaciers distinct from the Greenland and Antarctic ice sheets cover an area of ∼706,000 km 2 globally (RGI, 2017), with an estimated total volume of 160,000 cubic kilometers, or 0.32 ± 0.08 m of potential sea-level rise (Farinotti et al., 2019). They form where snow accumulation exceeds the melt of snow, firn, and ice over periods of several consecutive years. A climatic change directly influences the energy and mass balance at the glacier surface. Over longer time periods, changes in mass balance cause volume and thickness changes, which in turn affect the flow of ice via internal deformation and basal sliding (Cuffey & Paterson, 2010). This dynamic reaction eventually leads to glacier length changes, that is, the advance or retreat of glacier tongues. As a consequence, glacier front variations represent an indirect, delayed, and filtered response to changes in climate over glacier-specific reaction times of years to decades and response times of up to several decades or even centuries (Haeberli & Hoelzle, 1995;Jóhannesson et al., 1989;Oerlemans, 2007).
ZEMP AND MARZEION 10.1029/2021GL092964 2 of 11 Figure 1. The impact of volcanic eruptions on climate and glaciers. Schematic diagram of volcanic inputs to the atmosphere and their effects on climatic and glacier mass changes in relation to the explosivity (y-axis) of and spatial distance (x-axis) to the eruption. At local to regional scales, glaciers react with opposed mass-balance anomalies: 1. Negative due to ice removal from the blast and melting by hot lava (first mountain in red sector), 2. Positive due to insulation under thick layers of ash (second mountain in blue sector), and 3. Negative due to reduced snow and ice albedo caused by thin ash layers (third mountain in red sector). Eruptions reaching the stratosphere result in 4. (fourth mountain below blue sector): positive mass-balance anomalies due to reduced solar flux and net cooling at the land surface. Figure modified after Robock (2000).
Fluctuations of glaciers have been compiled and disseminated by the World Glacier Monitoring Service (WGMS) and its predecessor services for more than 125 years (Allison et al., 2019;WGMS, 1998). We used global glacier mass changes from a recent assessment (Zemp et al., 2019b) as dependent variable for our multiple regression model (cf., section 2.7). This data set combines the temporal variability from glaciological field measurements (from 450 glaciers) with glacier-specific change rates from air and spaceborne geodetic surveys (from 19,130 glaciers), as compiled by the WGMS (WGMS, 2020b;Zemp et al., 2015). These calibrated annual time series were regionally extrapolated from the observational to the full glacier sample (215,000 glaciers) by spatially interpolating the individual specific mass changes to all glacier locations in the region. The regional (area-weighted) mean specific mass changes (in meters water equivalent) were multiplied by regional glacier surface areas (in km 2 , RGI, 2017) to get regional mass changes in gigatonnes (1 Gt = 10 12 kg). Global sums provide glacier mass changes for the hydrological years from 1961/1962 to 2015/2016. Error bars of these observations correspond to 95% confidence bands. Full details on the methodology are found in Zemp et al. (2019b) and the global glacier mass-change estimates are available from the Zenodo repository (Zemp et al., 2019a). For the discussion of glacier front variations (cf., section 3.4), we used the global data set from direct observations (Hoelzle et al., 2003;Leclercq et al., 2014) and reconstructions (Zemp et al., 2011), consisting of more than 47,500 observations of 2,500 glaciers ( Figure S1; WGMS, 2020b; Zemp et al., 2015).

Climate Forcings
We used the time series of radiative climate forcings as prepared for the Coupled Model Intercomparison Project (CMIP5) (Taylor et al., 2012) for the extended historical period from 1850 to 2012. Global-average instantaneous forcings at the tropopause, as calculated for the GISS E2-R NINT ensemble (R. L. Miller et al., 2014), are available as anomalies relative to the net radiative flux in 1850 in Watt per square meter (Wm −2 ). We combined these individual forcings into sums of volcanic (stratospheric aerosols), solar (orbital variations, solar irradiance), and anthropogenic (well-mixed greenhouse gases, ozone, land use, snow albedo, direct and indirect tropospheric aerosols) drivers ( Figure S2). Full details on these forcings are found in R. L. Miller et al. (2014). To create the scenario for the repetition of the Pinatubo eruption, we set the forcing from stratospheric aerosols to zero but then repeated the negative forcing from the 5 years following the 1991 eruption at intervals of 45 years, that is, in 1856, 1901, 1946, 1991, and 2036. Further, we assumed a continued increase of the anthropogenic forcings by extending the linear trend from the 30-years period from 1983-2012 to 2040. For simplicity and in line with its limited influence (cf., Section 3.1), we kept the anomaly in solar forcings to zero after 2012.

Volcano Names and Stratospheric Aerosol Optical Depth
Volcano names and eruption years are taken from the Volcanos of the World database (from the field StartDate; Global Volcanism Program, 2013). Reconstructions of global mean stratospheric aerosol optical depth at 550 nm come from the EVA(eVolv2k) reconstruction of Toohey and Sigl (2017b) and from Sato et al. (1993Sato et al. ( , updated 2012 for the time before and after 1850, respectively. We note a discrepancy in the timing of the eruption of Sheveluch (Kamchatka) between the eruption year 1854, as given in the Volcanos of the World database (Global Volcanism Program, 2013), and the increase of volcanic aerosols (with its peak in 1857) in Sato et al. (1993Sato et al. ( , updated 2012.

Global Temperatures
We used global temperatures from 1961 to 2005 from the HadCRUT4 data set (Morice et al., 2012) as developed by the Climate Research Unit (CRU) of the University of East Anglia and the Hadley Center, UK, combining land (CRUTEM4) and marine (HadSST3) observations for the global mean temperature. For the analysis of temperature at the glacier locations, we used the CRU Time-Series (TS) version 4.01 of high-resolution gridded data of month-by-month variation in climate (Harris & Jones, 2017). From these data, expressed as anomalies from 1961 to 1990, we calculated mean annual global temperatures and mean annual temperatures at glacier locations, weighted by the glaciers' areas (RGI, 2017). Temperature increase was calculated from the slope of linear trends over corresponding time periods.

ENSO Information
The multivariate ENSO index (MEI) is a weighted average of the main ENSO features contained in sea-level pressure, surface wind, surface sea and air temperatures, and cloudiness (Wolter & Timlin, 1998). The values are normalized for each of 12 sliding bi-monthly seasons (December/January, January/February, …, November/December) so that the values over the full period  have an average of zero and a standard deviation of one. We used the extended version of MEI (Wolter & Timlin, 2011) with data from 1871 to 2005 and calculated annual means from the bi-monthly seasons. For the model runs over the full period, we assumed the index to be zero before 1871 and after 2005.

Superposed Epoch Analysis
We used a superposed epoch analysis (cf., Singh & Badruddin, 2006) to detect significant positive masschange anomalies in the years after large volcanic eruptions ( Figure S3). We first calculated mass-change anomalies with respect to the arithmetic average calculated over a time period covering the 5 years before and after each eruption in order to remove the general trend toward more negative mass balances (which is most pronounced over the Pinatubo period). From the superposed time series of these mass-balance anomalies, we then calculated the mean anomaly for each post-eruption year and used a t-test to detect significant positive mean anomalies with respect to the pre-eruption years.

Multiple Regression Model
Following the approach by Lean and Rind (2008), we reconstructed global annual glacier mass changes ΔM as: from time series of ENSO E, solar irradiance S, volcanic aerosols V, and anthropogenic influences A. The fitted coefficients c 0 for intercept and c for the independent variables, related uncertainties, and the correlation matrix were obtained by a multiple linear regression model using ordinary least squares estimation (Table S1). Linear correlation of the individual independent variables was tested for different lags (Δt, in years, Table S2). The introduction of lags can help to maximize the explained variance (Lean & Rind, 2008) but complicates the interpretation of the results (Muryshev et al., 2017). Consequently, we set all lags to zero for the main model run. In addition to global glacier mass changes (Figures 2, S4 and S5), we ran the model with the same independent forcings to reconstruct mean annual global temperatures (Figures S6a−S6e) and mean annual temperatures at glacier locations ( Figures S6f−S6j). The overall model uncertainty is expressed by the r-squared value (adjusted for multiple independent variables), which indicates the percentage of the variance in the observations that can be explained by the (multivariate) model, and by the probability of the related F-statistics. The significance of the independent variables is indicated by the probability of the related t-test. We provide error bars corresponding to 95% confidence bands of the overall regression model and of the individual coefficients of the independent components. For multi-annual periods, annual uncertainties are considered to be dependent and, hence, cumulated linearly. We note that this conservative assumption results in rather large multi-year uncertainties.

Attribution of Global Glacier Mass Changes to Anthropogenic and Natural Forcings From 1961 to 2005
Mass balance is the most direct reaction of a glacier to atmospheric conditions and, hence, most appropriate for quantifying the impact of climatic changes. Extrapolations of glaciological and geodetic observations (Zemp et al., 2019b) show that global glacier mass changes cumulated to −6,000 Gt, or 17 mm sea-level equivalent, over the observation period from 1961/62 to 2004/05. Over the same time period, annual air temperatures-which are strongly correlated with glacier mass changes (Braithwaite, 1981;Ohmura, 2001)-increased by 0.6 K for the global mean and by 1.2 K at glacier locations ( Figure 2a). From the individual components (Figures 2c−2f), anthropogenic influences and volcanic aerosols are the main drivers of global glacier mass changes (exceeding 99% confidence levels). The anthropogenic component shows a linear increase from about −75 ± 45 to −250 ± 150 Gt yr −1 over the observation period. The eruptions of Agung, El Chichón, and Pinatubo resulted in positive mass changes with peaks of 126 ± 75, 122 ± 72, and 203 ± 120 Gt yr −1 , respectively. Solar irradiance and ENSO vary with standard deviations of less than 10 Gt yr −1 around mean values close to zero and do not significantly contribute to the model. This lack of correlation is noteworthy since both solar irradiance and ENSO show significant contributions in a corresponding model applied to mean annual global temperatures (Figures S6a−S6e). When applying the model to mean annual temperatures at glacier locations, however, anthropogenic and volcanic forcings remain the only ones contributing to the explained variance, with solar irradiance and ENSO becoming insignificant again (Figures S6f−S6j).
The limited influence of solar irradiance at global scale can be explained by the small variability in instantaneous radiative forcing (±0.05 Wm −2 ) that is one to two orders of magnitude smaller than the changes in anthropogenic or volcanic forcings ( Figure S2). Also, the spatial pattern of the planetary energy balance shows that solar radiation dominates the influx of energy at low latitudes, while the energy balance at high latitudes-where most glaciers are located-depends more strongly on the meridional transport of energy within the Earth's climate system (Stephens & L'Ecuyer, 2015). The missing representation of meridional transport of energy (Liu et al., 2020) plausibly explains the underestimation of natural variability in our model. Similarly, the spatial patterns associated with ENSO show dominant impacts between ±30° latitudes but distinct meridional asymmetry in the Northern hemisphere (Lean & Rind, 2008), and probably in the Arctic as well. We also tested annual precipitation (both as global means and at glacier locations; not shown) as a potential driver representing internal climate variability. However, the corresponding model runs showed no significant contributions, which indicates that precipitation rather plays a role at local to regional and at seasonal scales (Fujita, 2008;Tejedor et al., 2021;Thibert et al., 2018).

Impact of Large Volcanic Eruptions on Glacier Mass Changes Since 1850
To assess the impact of major eruptions over time, we ran the multiple regression model for the entire period of available forcings back to 1850 ( Figure S4). This experiment shows that the major post-LIA volcanic eruptions resulted in positive anomalies of global glacier mass changes with distinct annual peaks and durations of a few years. As such, the 1883 eruption of Krakatau, Indonesia, shows a glacier mass gain of 440 ± 260 Gt accumulated over three years, with the annual maximum in 1884 of 250 ± 148 Gt ( Figure S4, Table S3). Pinatubo featured a 1-year delayed annual peak of glacier mass gain of about 200 ± 120 Gt, which is about twice as large as the corresponding peaks following the eruptions of Sheveluch (1854, Kamchatka), Santa Maria (1902, Guatemala), Agung and El Chichón. The anthropogenic component shows a steady increase of mass-loss rates from 0 ± 0 Gt yr −1 in 1850 to −275 ± 160 Gt yr −1 at present, with distinct acceleration after 1950. As a result, the large volcanic eruptions of the 19th and early 20th century led to distinct periods of a few years of global glacier mass gain, whereas the eruption of Agung just resulted in 1 year of a slightly positive mass balance of 30 ± 71 Gt yr −1 in 1964. However, Pinatubo only managed to bring glaciers close to steady-state conditions in 1992. Meanwhile, the continuously increasing impact of the anthropogenic forcing is so strong that a hypothetic repetition of the Pinatubo eruption in the coming decades would still result in a negative global glacier mass balance. We demonstrate this in a controlled model run setting the volcanic forcing to zero but repeating the negative forcing of the Pinatubo eruption at 45-years intervals in 1856,1901,1946,1991, and 2036 ( Figure S5, Table S3). This experiment clearly demonstrates the dwindling relevance of volcanic eruptions for the global net balance of glaciers. While the hypothetic eruptions in the 19th century resulted in distinct periods of global mass gain over a few years, a repetition of the Pinatubo eruption in 2036 would be nowhere near able to compensate for the ice loss due to anthropogenic forcings. For this future scenario, an event with two to three times the radiative effect of the Pinatubo eruption would be required for a global net mass gain of glaciers. Global glacier changes from the combined natural and anthropogenic forcings sum up to −10,600 ± 10,900 Gt for the period from 1850 to 2005 (Table S4), corresponding to a sea-level equivalent (SLE) of about 0.03 ± 0.03 m (incl. Antarctic periphery). This is only about half of earlier estimates which were based on scaling from observed front variations (0.058 ± 0.014 m SLE) (Leclercq et al., 2011;Marzeion et al., 2015) and on a numerical model driven by observed climate (0.063 ± 0.008 m SLE) (Marzeion et al., 2012(Marzeion et al., , 2015, both excluding glaciers in the Antarctic periphery. A direct comparison of these approaches is not straightforward and our lower values are not necessarily in contradiction to the other studies. As such, our model does not account for internal climate variability and glacier imbalance from before 1850. Corresponding sea-level contributions from these two components have recently been estimated from numerical modeling to be 0.030 ± 0.013 m SLE for the period from 1900 to 2005 (Slangen et al., 2016). This could explain the observed gap and why the majority of glaciers have been retreating since 1850 ( Figure S1c), while our reconstructed mass balances are only marginally negative from 1850 to 1900 (and-to lesser degree-up to 1950). Furthermore, our model suggests that 95 ± 52% of the explained variance in the global glacier mass changes since 1850 can be attributed to anthropogenic influences (Table S4). However, our statistical model explains only 43% of the variance in the observations  and does neither consider natural variability nor glacier imbalance from before 1850. Consequently, our results cannot directly be compared to attribution studies using numerical models that include such internal forcings. As such, Marzeion et al. (2014) attributed 25 ± 35% of post-LIA glacier changes to anthropogenic causes, while Roe et al. (2021) found that the anthropogenic component of the centennial mass loss is essentially 100%. Their differences are explained by a potential bias in the estimated glacier response times and different assumptions about how to account for the glacier disequilibrium with concurrent climatic conditions (Roe et al., 2021). Studies agree that anthropogenic causes are the dominant driver of the temperature increase in the industrial era (Allen et al., 2018) and of global glacier mass loss since at least about 1970 (Marzeion et al., 2014;Slangen et al., 2016), but it remains open if the same holds true for global glacier mass loss since 1850 (Roe et al., 2017).

Limitations of the Multivariate Analysis of Recent Mass-Balance Observations
Our multiple regression model is subject to limitations from several sources. First, our model is calibrated with global glacier mass changes estimated from available glaciological and geodetic samples that might not be fully representative in all regions (Zemp et al., 2019b). Second, our model is calibrated over the second half of the 20th century when volcanic eruptions were relatively small. The stratospheric aerosol optical depth after the eruption of Pinatubo reached only about half the values of the major eruptions of the past 2,500 years (Toohey & Sigl, 2017b). As a consequence, the volcanic impact on global glacier mass changes was probably more important for major eruptions in times of limited anthropogenic influences. Third, the global forcings might not fully cover regional effects. As such, the dimming of solar radiation from the 1950s until the 1980s and the subsequent brightening during the 1990s (Wild et al., 2005) have been shown to be in line with corresponding positive and negative mass-change anomalies, respectively, in the European Alps (Huss et al., 2009). The calibration of our model during the longer period of dimming, if not fully reflected in the global solar forcing (R. L. Miller et al., 2014), might have led to a positive bias in the resulting mass changes. Fourth, the (significant) forcings do not include internal climate variability (or "noise"), which could explain additional inter-annual to decadal variability of global glacier changes (Oerlemans, 2000;Roe et al., 2017). Fifth, our model relates instantaneous forcings to annual glacier mass changes neither considering potential lags in the forcings (Lean & Rind, 2008;Muryshev et al., 2017) nor committed mass changes due to the imbalance of glaciers with climatic conditions (Marzeion et al., 2018;Mernild et al., 2013;Slangen et al., 2016). Finally, the climatic impact of volcanic eruptions might be mediated through atmospheric or oceanic feedbacks which in turn depend on the seasonal timing of the eruption (Stevenson et al., 2017), the initial atmospheric conditions (Pausata et al., 2016), and on the hemispheric distribution of the aerosol optical depth (Colose et al., 2016).
Future studies might also look into the impact of seasonal timing and latitudinal location of the eruption on glaciers (cf., Gao et al., 2008;Sato et al., 1993). While strong eruptions are expected to have a relatively uniform global impact due to the larger amount and the longer lifetime of stratospheric aerosols, we would expect that global glacier mass changes are also sensitive to weaker eruptions taking place in the northern hemisphere, and before the ablation period. In addition, it might be worthwhile to include temperature at glacier locations-after removing the effect of all other forcings used (cf., Foster & Rahmstorf, 2011)-as an independent variable in the model in order to increase the explained inter-annual variability of glacier mass changes. Finally, it would be worthwhile to repeat the present study at regional scale, which would allow for testing the significance of the various forcings at the regional scale and for investigating potential ZEMP AND MARZEION 10.1029/2021GL092964 competing effects such as the occurrence of El Niño events during/after volcanic eruptions (Emile-Geay et al., 2008;Stevenson et al., 2016).

Can Glacier Front Variations Extend the Analysis Back in Time?
Observational estimates of global glacier mass changes are only available for the second half of the 20th century when volcanic eruptions were relatively small. However, glacier front variations can potentially be used to extend the analysis back in time. Corresponding observations have been collected from direct measurements since the 19th century and extended from reconstructions back into the LIA (Zemp et al., 2011(Zemp et al., , 2015. The global compilation of glacier front variations is illustrated by long reconstruction series at selected glaciers and as a ratio of advancing or retreating glaciers since 1550 ( Figure S1). It shows a general trend of frontal retreat at the centennial scale with intermittent decadal periods of re-advances for up to 40% of the glacier sample. Major volcanic eruptions, such as Huaynaputina (Peru) in 1600, Komagatake (Japan) and Parker (Philippines) in 1640, or Tambora (Indonesia) in 1815, coincide with periods of glacier advances. However, the increase in the number of advancing glaciers typically starts years before the eruptions, and other distinct periods of glacier advances (e.g., around 1975) are without preceding major eruptions. From process understanding and in view of the time needed for a glacier tongue to react to a climatic change (i.e., years to decades Haeberli & Hoelzle, 1995;Jóhannesson et al., 1989;Oerlemans, 2007), we do not expect that a short pulse of mass gain, lasting only one or a few years, will result in a distinct period of major global glacier advances. However, such advances are well possible when large volcanic eruptions coincide with other forcings, internal climate variability, and/or when they contribute to positive feedback mechanisms (e.g., related to sea ice, atmospheric or oceanic circulation patterns) with decadal-scale persistence (Berdahl & Robock, 2013;G. H. Miller et al., 2012;Zanchettin et al., 2012). Future modeling studies are needed to quantify the impacts of large volcanic eruptions on glacier mass balance and front variations at local, regional, and global scales.

Conclusions
This study is a first attempt to quantify the impact of stratospheric aerosols from large volcanic eruptions on glacier changes at global scale. We performed multiple linear regression analyses to decompose an observational data set of global glacier mass changes into components associated with anthropogenic influences, volcanic events, solar activities, and ENSO. Our empirical model suggests that the global glacier mass loss over the industrial era is dominantly driven by anthropogenic causes modified by short periods of positive mass-change anomalies following large volcanic eruptions. With increasing greenhouse gas concentrations, however, the relevance of volcanic eruptions is dwindling; future eruptions are unlikely to result in a net gain of glacier mass at global level. Furthermore, our study indicates that solar activity and ENSO have limited impacts on climate conditions at glacier locations, and that volcanic eruptions alone can hardly explain decadal periods of glacier advances documented since the 16th century.

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.

Data Availability Statement
The reconstructed global mass changes of the present study, including components from the different forcings, have been deposited in the Zenodo repository (https://doi.org/10.5281/zenodo.4706597). The full samples of front-variation and mass-balance observations for individual glaciers are publicly available from the WGMS (WGMS, 2020a). The global glacier mass-change estimates (Zemp et al., 2019b) for the period 1961-2016 are available from the Zenodo repository (Zemp et al., 2019a). Climate forcings are available from NASA/GISS (2014), stratospheric aerosol optical depth from the World Data Center for Climate (Toohey & Sigl, 2017a) and from NASA/GISS (2016), global temperatures from CRU (2020), and ENSO from NOAA (2011). The analytical scripts are available from the authors on request.