Climate Projections Very Likely Underestimate Future Volcanic Forcing and Its Climatic Effects

Standard climate projections represent future volcanic eruptions by a constant forcing inferred from 1850 to 2014 volcanic forcing. Using the latest ice‐core and satellite records to design stochastic eruption scenarios, we show that there is a 95% probability that explosive eruptions could emit more sulfur dioxide (SO2) into the stratosphere over 2015–2100 than current standard climate projections (i.e., ScenarioMIP). Our simulations using the UK Earth System Model with interactive stratospheric aerosols show that for a median future eruption scenario, the 2015–2100 average global‐mean stratospheric aerosol optical depth (SAOD) is double that used in ScenarioMIP, with small‐magnitude eruptions (<3 Tg of SO2) contributing 50% to SAOD perturbations. We show that volcanic effects on large‐scale climate indicators, including global surface temperature, sea level and sea ice extent, are underestimated in ScenarioMIP because current climate projections do not fully account for the recurrent frequency of volcanic eruptions of different magnitudes.

• There is a 95% chance that the time-averaged 2015-2100 volcanic SO 2 flux from explosive eruptions exceeds the time-averaged 1850-2014 flux • Standard climate projections very likely underestimate the 2015-2100 stratospheric aerosol optical depth and volcanic climate effects • Small-magnitude eruptions (<3 Tg SO 2 ) contribute 30%-50% of the volcanic climate effects in a median future eruption scenario

Supporting Information:
Supporting Information may be found in the online version of this article.
Up-to-date ice-core and satellite volcanic sulfur emission datasets enable us to account for the occurrence of (a) eruptions larger in magnitude than those that occurred between 1850 and 2014, which injected as much as 300 Tg of SO 2 into the atmosphere, and (b) small-magnitude eruptions below the detection threshold of ice-core datasets (Figure 1a), which can contribute a significant fraction to stratospheric aerosol optical depth (SAOD) (Santer et al., 2014;Schmidt et al., 2018).
In addition, whether they apply a constant volcanic forcing (e.g., CMIP6 ScenarioMIP) or use stochastic eruption scenarios (Bethke et al., 2017), existing climate projections use prescribed volcanic aerosol optical properties derived from simplified volcanic aerosol models. Climate models with interactive stratospheric aerosols (Timmreck et al., 2018) showed a better agreement between the simulated surface temperature responses and tree-ring surface temperature reconstructions for the 1257 Mount Samalas and 1815 Mount Tambora eruptions (Stoffel et al., 2015) and the 1783-1784 Laki eruption (Pausata et al., 2015;Zambri et al., 2019). Furthermore, the prescribed aerosol approach cannot account for the impacts of global warming on the life cycle of volcanic sulfate aerosols (Aubry et al., 2021), including the impact of changing atmospheric stratification on volcanic plume height (Aubry et al., 2019). Such climate-volcano feedbacks might amplify the peak global-mean radiative forcing associated with large-magnitude tropical eruptions by 30% (Aubry et al., 2021).
Our study aims to improve our understanding of future volcanic impacts on climate. To this end, we perform model simulations from 2015 to 2100 with two innovations: (a) a stochastic resampling approach using the latest ice-core and satellite datasets to generate improved future volcanic eruption scenarios; and (b) a plume-aerosol-chemistry-climate modeling framework (named UKESM-VPLUME), which combines a volcanic plume model and an Earth System Model with interactive stratospheric aerosols to simulate volcanic climate effects while accounting for climatic controls on plume-rise height.

Stochastic Future Eruption Scenarios
We generate 1,000 stochastic future eruption scenarios for 2015-2100 by resampling SO 2 mass from volcanic emission inventories from a bipolar ice-core array covering the past 11,500 years (Holvol; Sigl et al., 2022) and a multi-satellite record from 1979 to 2021 (S. Carn, 2022;S. A. Carn et al., 2016) (Figure 1a and Figure S1 in Supporting Information S1). Before resampling, we filter out: (a) effusive eruptions; (b) in the satellite record, eruptions with eruptive plume heights more than 3 km below the thermal tropopause (obtained from NCEP/ NCAR Reanalysis 1; Kalnay et al., 1996); we assume that aerosol lofting could result in stratospheric injections for tropospheric plumes less than 3 km below the tropopause. By examining the eruption frequency-magnitude (i.e., in this study, SO 2 mass) distribution of both ice-core and satellite records (Figure 1a), we identify 3 Tg of SO 2 as a threshold: (a) below which ice-core records underestimate eruption frequency due to under-recording; and (b) above which the short duration of the satellite record precludes it from capturing the true frequency of eruptions with higher magnitude. Accordingly, we use a 3 Tg of SO 2 threshold to define "small-magnitude" and "large-magnitude" eruptions. We resample small-magnitude eruptions from the satellite record only, and large-magnitude ones from the combined ice-core and satellite record. Details of the resampling procedures are discussed in Supporting Information S1.

UKESM-VPLUME
Atmospheric stratification, wind and humidity affect volcanic plume dynamics and SO 2 injection height (e.g., Mastin, 2014), but SO 2 height is commonly prescribed in modeling studies of volcanic forcing (e.g., Timmreck et al., 2018). To account for meteorological controls on plume dynamics, we have developed UKESM-VPLUME, which couples the UK Earth System Model (UKESM; Mulcahy et al., 2023) with Plumeria (1-D eruptive plume model; Mastin, 2007Mastin, , 2014 (details in Supporting Information S1). We use version 1.1 of UKESM with fully-coupled atmosphere-land-ocean and interactive stratospheric aerosols. In brief, for each time step of the UKESM atmospheric model during an eruption, UKESM-VPLUME interactively passes the atmospheric conditions simulated at the eruption location to Plumeria. Plumeria then computes the neutral buoyancy height of the volcanic plume based on atmospheric conditions and the mass eruption rate generated for each eruption in the stochastic scenarios. Volcanic SO 2 is injected into UKESM at the neutral buoyancy height calculated in Plumeria using a Gaussian profile with a width of 10% of the plume height (consistent with large-eddy simulations of volcanic plumes, Aubry et al., 2019). This approach ensures that plume heights of volcanic eruptions are consistent with the meteorological conditions simulated by UKESM.

Experimental Design
We perform simulations using the UKESM-VPLUME framework for four stochastic future eruption scenarios at the 2.5th, 50.0th, 50.5th and 98.0th percentiles (termed VOLC2.5, VOLC50-1, VOLC50-2, VOLC98) of the distribution of the 2015-2100 average SO 2 flux across the 1,000 future eruption scenarios ( Figure 1b). We choose scenarios close (within 0.5 percentile) to the 2.5th, 50th and 97.5th to sample the median and 95% confidence interval of the future volcanic stratospheric SO 2 injections. To test future climate trajectory sensitivity to the temporal and spatial distribution of eruptions, we run two scenarios near the 50th percentile. For instance, VOLC50-2 has more large-magnitude eruptions than VOLC50-1 in the early 21st century (Figure 1c). We also performed the VOLC50 runs with small-magnitude eruptions only (VOLC50-1S and VOLC50-2S) to isolate their contribution to the overall climate effects caused by eruptions of all magnitudes. We compare the results from  Figure 2 shows the global monthly-mean SAOD at 550 nm and the time-averaged values over 2015-2100. The time-averaged ensemble-mean SAOD ranges from 0.015 ± 0.0004 (VOLC2.5, one standard deviation uncertainty) to 0.062 ± 0.0018 (VOLC98), with an average value of 0.024 ± 0.0012 for the two median future eruption scenarios (VOLC50), while the SAOD in CONST, which followed the ScenarioMIP design, is 0.012 ± 0.0018. Small-magnitude eruptions contribute 0.010-0.013 ± 0.0002 to the time-averaged SAOD in the VOLC50 scenarios, that is, about 50% of the total SAOD. Comparing VOLC2.5 to CONST and assuming that the rank for the 2015-2100-year mean volcanic SO 2 flux and SAOD are the same, it is thus very likely (i.e., >90% probability following IPCC guidance note; Mastrandrea et al., 2010) that the actual global 2015-2100 mean SAOD will be higher than that prescribed in ScenarioMIP, with the median (VOLC50) SAOD value being double that used in ScenarioMIP. The result is consistent with Figure 1, given that the 1850-2014 time-averaged SO 2 flux used to define the ScenarioMIP volcanic forcing is close to the 2.5th percentile of the future volcanic SO 2 flux distribution. Beyond the time-averaged SAOD value, owing to the sporadic nature of volcanic eruptions, the global monthly-mean SAOD values in VOLC scenarios can be up to a factor of 60 greater than that in ScenarioMIP ( Figure 2 and Figure S2 in Supporting Information S1). Figure 3a shows the global annual-mean surface air temperature at 1.5 m (GMST) relative to the 1850-1900 period. Large-magnitude volcanic eruptions lead to a short-term drop in the annual-mean GMST for at least 1 year and up to 6-7 years for the largest eruptions. In the VOLC98 scenario where clusters of large-magnitude eruptions occur, they can induce multi-decadal global cooling. The 2015-2100 time-averaged GMST relative to detrended NOVOLC ensemble mean (Figure 3b) ranges between −0.16°C (VOLC2.5) and −0.56°C (VOLC98), with CONST lying outside this range at −0.12°C. Volcanic cooling for median eruption scenarios (VOLC50-1 and VOLC50-2) is 0.20-0.24°C, double that of CONST, and 0.09-0.10°C of cooling is attributable to small-magnitude eruptions (Table S1 in Supporting Information S1).

Results
The IPCC defines global warming as an increase, relative to 1850-1900, in the global mean surface air and sea surface temperatures over a period of 30 years (IPCC, 2021). Using this definition, we examine the year of crossing of 1.5, 2, and 3°C warming thresholds for VOLC and CONST runs (Figure 3c). Volcanic eruptions delay the time of crossing 1.5°C by about 1.6-3.2 years when compared to NOVOLC (Table S2 in Supporting Information S1), consistent with Bethke et al. (2017). Compared to CONST, times of temperature threshold crossings are significantly delayed by 1.8-2.5 years in VOLC50-2, but unaffected in VOLC50-1. This highlights the sensitivity of the time of crossing to the temporal distribution of large-magnitude eruptions. The occurrence of volcanic clusters in VOLC98 causes an extended cooling period between 2034 and 2060 ( Figure 3a) which delays the crossing of 2 and 3°C by 7 and 14 years, respectively.
In Figure 4, we examine volcanic effects on large-scale climate indicators other than GMST. The 2015-2100 time-averaged global annual-mean precipitation fluxes in all VOLC runs show a greater reduction than CONST, with a range between −0.014 mm/day (VOLC2.5) to −0.052 mm/day (VOLC98), and −0.010 mm/day for CONST (Figure 4a). In VOLC50 scenarios, the global annual-mean precipitation flux is reduced by 0.019 mm/ day with small-magnitude eruptions alone contributing between 0.008 and 0.009 mm/day, comparable to the effects of the volcanic forcing implemented in ScenarioMIP. It is thus very likely that the reduction of global mean precipitation due to volcanic effects is underestimated in ScenarioMIP.
Volcanic-induced surface cooling penetrates into the deep ocean layer and decreases the global ocean heat content (Figure 4b and Figure S3 in Supporting Information S1), which in turn leads to less thermal expansion in seawater and a reduction in thermosteric sea level (Figure 4c). Volcanic forcing in VOLC50 reduces global ocean heat content and thermosteric sea level by 6%-7% compared to NOVOLC by 2100, whereby about half is attributed to small-magnitude eruptions ( Figure S3 in Supporting Information S1). Although volcanic forcing can cause considerable impacts on large-scale ocean metrics, it does not offset the anthropogenic-induced ocean warming trends even for the upper-end volcanic emission scenario VOLC98 (Figures 4b, 4c, and 4f).
Depending on the eruption magnitude and location, the Arctic and Antarctic sea ice extents show an immediate increase for 1-2 years after large-magnitude eruptions (Figures 4d and 4e). The time-averaged global  sea ice extent in VOLC runs over 2015-2100 increases by 0.43 million km 2 (VOLC2.5) to 1.53 million km 2 (VOLC98) as compared to 0.20 million km 2 for CONST. Comparing VOLC2.5 to CONST suggests that for similar time-averaged SAOD, the use of a constant forcing instead of a stochastic eruption distribution halves the magnitude of the sea ice response.
The time-averaged Atlantic Meridional Overturning Circulation (AMOC) at 26°N over 2015-2100 is strengthened by between 0.26 Sv (VOLC2.5) and 0.93 Sv (VOLC98) as compared to NOVOLC, with all VOLC scenarios exhibiting an increased decadal mean AMOC strength (Figure 4f). Small-magnitude eruptions alone can increase the time-averaged AMOC strength by 0.36-0.38 Sv (VOLC50-1S and VOLC50-2S), which is greater than CONST at 0.28 Sv, and contribute to over 77% of the AMOC response in the median future scenarios (Table  S1 in Supporting Information S1). One of the median future scenarios (VOLC50-1) has a weaker time-averaged AMOC response than the same run with small-magnitude eruptions only (VOLC50-1S) due to an extended period of weakened AMOC after the occurrence of large-magnitude eruptions with aerosols dispersed to the Southern Hemisphere ( Figures S2 and S4 in Supporting Information S1). Our results suggest that AMOC may have different responses caused by differences in the asymmetries of the stratospheric aerosol burdens from large-magnitude eruptions, which warrants further investigation.

Discussion
Small-magnitude eruptions (<3 Tg of SO 2 ) contribute a considerable fraction (between 33% and 40%) of the total upper atmospheric volcanic SO 2 emissions in VOLC50, and in turn, are responsible for 30%-50% of the volcanic impact on selected large-scale climate indicators and over 77% of the AMOC response ( Figure 5 and Table S1 in Supporting Information S1). For future eruption scenarios with fewer eruptions than VOLC50, the contribution from small-magnitude eruptions is expected to be even greater because the total mass injected by small-magnitude eruptions is relatively similar across all scenarios. Despite the importance of volcanic forcing from small-magnitude eruptions, they are mostly unaccounted for in historical simulations before satellite measurements are available. In the pre-satellite historical period , the Neely and Schmidt (2016) and Sigl et al. (2022)  Our stochastic scenarios imply that CMIP6 ScenarioMIP very likely (95 ± 2.5%) underestimates the 2015-2100 volcanic SO 2 flux from explosive eruptions and, in turn, forcing ( Figure 1b). Figure 1b shows the cumulative probability against the annual SO 2 flux obtained by resampling ice-core record of volcanic SO 2 injection only (i.e., Holvol; Sigl et al., 2022) and both ice-core and satellite (S. Carn, 2022) records as in our stochastic scenarios. CMIP6 ScenarioMIP uses a constant volcanic forcing inferred from the 1850-2014 period during which the mean volcanic SO 2 flux recorded in emission inventories was 0.7 ± 0.06 Tg per year. However, we find a 95% confidence interval for the 2015-2100 mean volcanic SO 2 flux between 0.64 and 5.28 Tg per year in our eruption scenarios (Figure 1b). Our stochastic approach, which represents better the frequency-magnitude distribution of small-magnitude eruptions, results in a higher annual SO 2 flux than resampling from the ice-core record only (e.g., Bethke et al., 2017).
Our future volcanic eruption scenarios greatly enhance the variability of large-scale climate indicators as compared to the ScenarioMIP forcing ( Figure 5). Future volcanic emissions in our scenarios cause a 3.5% (VOLC2.5) to 15.0% (VOLC98) decrease in the 2015-2100 time-averaged global net radiative forcing at the top-of-the-atmosphere relative to the anthropogenic contribution ( Figure 5 and Figure S5 in Supporting Information S1). The time-averaged climate responses of our selected climate indicators scale with the magnitude of volcanic forcing except for the Antarctic sea ice extent and AMOC, which may depend on the latitudinal distribution of eruptions.
We also find that the magnitude of volcanic effects on climate indicators are comparable between CONST and VOLC2.5, which is a scenario with only one Pinatubo-like eruption over 2015-2100. Our results suggest that due to the low volcanic forcing used in ScenarioMIP, it is very likely (97.5%) that ScenarioMIP underestimates the climate effects of the large-scale climate indicators examined in this study.
Our simulation results show that for the SSP3-7.0 scenario, volcanic forcing can offset 2.1%-18.2% of the anthropogenic effects to large-scale climate indicators depending on the future eruption scenarios ( Figure 5). In a future scenario with low-end anthropogenic emission, we would expect the relative effect between future volcanism and anthropogenic forcing to be much greater. Our work highlights how the high level of uncertainty on volcanic forcing affects climate projections. For the same future eruption scenario, the volcanic effects on climate will also vary between SSP scenarios owing to climate-volcano feedbacks (e.g., Aubry et al., 2022;Fasullo et al., 2018;Hopcroft et al., 2018), which need to be quantified.

Conclusion
We performed climate model simulations from 2015 to 2100 with stochastic future eruption scenarios using UKESM-VPLUME (a plume-aerosol-chemistry-climate model framework that accounts for climate-volcano feedbacks) to examine how the uncertainties on volcanic forcing affect climate projections. Using the latest ice-core and satellite datasets, we show that the 2015-2100 volcanic SO 2 flux from explosive eruptions has a 95% probability to exceed the 1850-2014 flux, which was used to derive volcanic forcing in CMIP6 ScenarioMIP. Our simulations suggest that the time-averaged SAOD in a median future scenario is 0.024 (95% uncertainty: 0.015-0.062), which is double that in ScenarioMIP, and that ScenarioMIP very likely underestimates the future volcanic effects on climate. Our study emphasizes the importance of the climate effects of future volcanic eruptions relative to the anthropogenic contribution, which even for an upper end anthropogenic forcing scenario (SSP3-7.0) can range between 2.1% and 18.2% for large-scale climate indicators. We also highlight the climate-relevance of small-magnitude eruptions, which are responsible for 30%-50% of the volcanic effects on selected climate indicators. Future CMIP experiment designs should update historical forcings for the pre-satellite historical period  to include volcanic sulfur emissions from small-magnitude eruptions in this period, and ideally use interactive stratospheric aerosols that can explicitly simulate climate-volcano feedbacks. Future climate projection studies should either use stochastic eruption scenarios generated using state-of-the-art volcanic emission inventories or time-averaged constant forcing derived from such scenarios that better represents long-term volcanic activity. VPLUME framework. M.M. Chim is supported by the Croucher Foundation and The Cambridge Commonwealth, European & International Trust through a Croucher Cambridge International Scholarship. Thomas J. Aubry acknowledges support from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 835939. Jane Mulcahy and Jeremy Walton were supported by the Met Office Hadley Centre Climate Programme, funded by BEIS. Anja Schmidt acknowledges support from NERC grants NE/ S000887/1 (VOL-CLIM) and NE/ S00436X/1 (V-PLUS). This work used Monsoon2, a collaborative high-performance-computing facility funded by the Met Office and the Natural Environment Research Council. This work used JASMIN, the UK collaborative data analysis facility.