# Radiative forcing by super-volcano eruptions

###### Abstract

We investigate the climatic effects of volcanic eruptions spanning from Mt.
Pinatubo-sized events to super-volcanoes. The study is based on ensemble
simulations in the Community Earth System Model Version 2 (CESM2) climate
model using the Whole Atmosphere Community Climate Model Version 6 (WACCM6)
atmosphere model. Our analysis focuses on the impact of different SO2-amount
injections on stratospheric aerosol optical depth (AOD), effective radiative
forcing (RF), and global temperature anomalies. Unlike the traditional linear
models used for smaller eruptions, our results reveal a non-linear
relationship between RF and AOD for larger eruptions. We also uncover a
notable time-dependent decrease in aerosol forcing efficiency across all
eruption magnitudes during the first post-eruption year. In addition, the
study reveals that larger as compared to medium-sized eruption events produce
a delayed and sharper peak in AOD, and a longer-lasting temperature response
while the time evolution of RF remains similar between the two eruption types.
When including the results of previous studies, we find that relating SO2 to
any other parameter is inconsistent across models compared to the
relationships between AOD, RF, and temperature anomaly. Thus, we expect the
largest uncertainty in model codes to relate to the chemistry and physics of
SO2 evolution. Finally, we find that the peak RF approaches a limiting value,
and that the peak temperature response follows linearly, effectively bounding
the temperature anomaly to at most $\sim$-12\text{\,}\mathrm{K}$$.

JGR: Atmospheres

UiT The Arctic University of Norway, Tromsø, Norway

Eirik R<EMAIL_ADDRESS>
The linear RF dependence on AOD breaks down for eruptions larger than Mt.
Pinatubo

The RF to AOD ratio has a time-after-eruption dependence on eruption latitude

Temperature and RF peak values has a linear dependence and reaches an upper
limit

## Plain Language Summary

Volcanic eruptions can have a significant impact on the Earth’s climate.
Eruptions large enough that the aerosols they emit reach the stratosphere
cause a cooling effect by reflecting sunlight. Typically, an eruption is
measured by its impact on the opacity of the stratosphere and the change in
the energy balance at the top of the atmosphere. The two measures are often
assumed to be linearly related, but the linearity is tested only against
eruptions seen in the last two millennia. We use a coupled climate model to
simulate the impact of eruptions of sizes up to the largest known eruptions.
The smallest eruptions we simulate are still large enough to cause global
climate effects. We find a clear non-linear relationship for eruptions larger
than the ones seen in the past two millennia. Our simulations and supporting
data shows that the eruption latitude significantly influences the development
of the relationship between energy imbalance and stratospheric opacity with
time after the eruption. Additionally, we find evidence that the peak energy
imbalance reaches a limit, and that the peak temperature response follows
linearly with the peak energy imbalance, also reaching a limiting value.

##  1 Introduction

Effective radiative forcing (RF) and stratospheric aerosol optical depth (AOD)
are crucial metrics representing the energy imbalance at top-of-the-atmosphere
(TOA) and the stratospheric opacity due to aerosol scattering, respectively.
They are extensively used to quantify the impact of major volcanic eruptions.
The assumption of a linear dependency of RF on AOD is commonly adopted [Myhre
. (2013), Andersson . (2015)], and applying such a linear relationship has
yielded reasonably accurate estimates in climate model simulations of volcanic
eruptions [Mills . (2017), Hansen . (2005), Gregory . (2016), Marshall .
(2020), Pitari . (2016)]. Yet, a wide spread in the estimated aerosol forcing
efficiencies (RF normalised by AOD) exists among studies, spanning
approximately from $\sim$-15\text{\,}\mathrm{W}\mathrm{m}^{-2}\ce{AOD}^{-1}$$
[Pitari . (2016)] to
$\sim$-25\text{\,}\mathrm{W}\mathrm{m}^{-2}\ce{AOD}^{-1}$$ [Myhre . (2013)].
Additionally, these estimates are predominantly based on small eruptions with
AOD values up to at most $\sim 0.7$.

Although H2O, N2, and CO2 are the most abundant gases emitted by volcanoes
[Robock (2000)], sulphur species such as SO2 provide a greater influence due
to the comparatively high background concentrations of the former gases in the
atmosphere. The transformation of SO2 molecules through reactions with OH and
H2O leads to the formation of sulphuric acid (H2SO4) [Robock (2000)], which
scatters sunlight thereby elevating planetary albedo and reducing the RF. As
the conversion from SO2 to H2SO4 occurs over weeks [Robock (2000)], the peak
RF experiences a slight delay from the eruption’s peak SO2 injection. The
lifetime of the H2SO4 aerosols in the stratosphere depends on various factors,
including latitude [Marshall . (2019), Toohey . (2019)], volcanic plume height
[Marshall . (2019)], aerosol size [Marshall . (2019)], the quasi-biennial
oscillation phase [Pitari . (2016)] and the season of the year (determining to
which hemisphere aerosols are transported) [Toohey . (2011), Toohey . (2019)].
In the case of tropical eruptions, aerosols are typically transported poleward
in the stratosphere and descend back to mid-latitude troposphere within one to
two years [Robock (2000)]. Upon descending below the tropopause, these
aerosols are readily removed by wet deposition [Liu . (2012)].

Before the current era of significant anthropogenic climate forcing, volcanic
eruptions were the primary forcing mechanism dictating Earth’s climate
variability during the Holocene period [Sigl . (2022)]. Despite this
substantial impact, few climate-model experiments have included volcanic
forcing when simulating climate evolution during the Holocene [Sigl . (2022)],
likely implying an exaggerated positive forcing [Gregory . (2016), Solomon .
(2011)]. This absence of persistent cooling is one of several factors that
have been suggested to contribute to the common disparity between simulated
and observed global warming [Andersson . (2015)]. Despite extensive attention
on understanding the way volcanic eruptions influence climate, questions
regarding aerosol particle processes—such as growth and creation rates when OH
is scarce—remain unanswered [<]e.g.
¿[]robock2000,zanchettin2019,marshall2020,marshall2022. These processes impact
aerosol scattering efficiency and potentially the RF to AOD relationship.
marshall2020 observe higher aerosol forcing efficiency in post-eruption years
$2$ and $3$ compared to year 1, and attribute this post-eruption increase in
aerosol forcing efficiency to strong spatial concentration in the initial year
and subsequent distribution of aerosols over a larger area. This spatial
redistribution increases the albedo per global mean AOD thereby causing a
stronger RF to AOD ratio [Marshall . (2020)].

Previous studies of both Mt. Pinatubo [Mills . (2017), Hansen . (2005)] and
volcanoes within the instrumental era [Gregory . (2016)] have been used to
estimate the relationship between the RF energy imbalance and change in AOD
caused by volcanic eruptions. While myhre2013 employ a formula scaling RF by
AOD to obtain $-25\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{AOD}^{-1}$, recent
literature reports estimates down to
$-19.0(5)\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{AOD}^{-1}$ [Gregory .
(2016)] and $-18.3(10)\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{AOD}^{-1}$
[Mills . (2017)]. Synthetic volcano simulations in marshall2020 yield a
scaling factor of
$-20.5(2)\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{AOD}^{-1}$ across an
ensemble of $82$ simulations featuring varying injection heights and latitudes
of volcanic emissions, with injected SO2 ranging from $10$ to
$100\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$.

A similar simulation setup, albeit with notable differences, was conducted by
niemeier2015, involving an ensemble of $14$ levels of injected sulphur
spanning between
$1\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{S}\mathrm{)}\mathrm{yr}^{-1}$
($2\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}\mathrm{yr}^{-1}$)
and
$100\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{S}\mathrm{)}\mathrm{yr}^{-1}$
($200\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}\mathrm{yr}^{-1}$).
These geoengineering simulations maintained continuous sulphur injections,
running until a steady sulphur level was achieved. Results indicated an
inverse exponential relationship between RF and injected SO2 rate, converging
to $-65\text{\,}\mathrm{W}\mathrm{m}^{-2}$ (Eq. 1). Even the $100\times$ Mt.
Pinatubo super-volcano simulation by jones2005, which obtained a peak RF of
$-60\text{\,}\mathrm{W}\mathrm{m}^{-2}$, is below the suggested limit of
$-65\text{\,}\mathrm{W}\mathrm{m}^{-2}$. Moreover, timmreck2010 find a peak RF
anomaly of $-18\text{\,}\mathrm{W}\mathrm{m}^{-2}$ from a
$1700\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$ eruption
simulation, which corresponds well with the function estimated by niemeier2015
at the given SO2 level. Several studies have demonstrated a linear
relationship of approximately
$-$20\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{AOD}^{-1}$$ between RF and AOD,
although substantial variability exists in the slope among studies [Mills .
(2017), Hansen . (2005), Gregory . (2016), Marshall . (2020), Pitari .
(2016)]. Moreover, a time-after-eruption dependence on the RF to AOD ratio is
found in marshall2020, whereas niemeier2015 revealed a non-linear relationship
between RF and injected SO2 rate. Thus, a consensus on the relationship
between injected SO2, AOD, and RF has yet to be established.

One avenue that has garnered considerable attention is comparing the magnitude
of volcanic or volcano-like forcings to increased CO2 levels. Several studies
explore the connection between volcanic forcing and the climate sensitivity to
a doubling of CO2 [Boer . (2007), Marvel . (2016), Merlis . (2014), Ollila
(2016), Richardson . (2019), Salvi . (2022), Wigley . (2005)]. The comparison
of forcing from volcanoes and CO2 aims to mitigate the large uncertainty in
estimates of the sensitivity of the real climate system. Inferring climate
sensitivity from volcanic eruption events has been attempted as a way to
constrain the sensitivity [Boer . (2007)] by assuming that volcanic and CO2
forcings produce similar feedbacks [Pauling . (2023)]. Earlier studies suggest
the potential for constraining equilibrium cilmate sensitivity (ECS) using
volcanoes [Bender . (2010)], provided that ECS is constrained by effective
radiative forcing (ERF) rather than instantaneous radiative forcing (IRF), as
ERF accounts for rapid atmospheric adjustments in contrast to IRF [Richardson
. (2019)]. However, other studies refute this approach, pointing out that
different sensitivities of volcanic forcing and CO2 doubling seem to exist
[Douglass . (2006)], or that constraining the ECS by ERF lacks accuracy due to
the precision of climate simulations [Boer . (2007), Salvi . (2022)]. Although
ERF offers a more suitable indicator of forcing than IRF [Marvel . (2016),
Richardson . (2019)], more recent studies conclude that ECS cannot be
constrained from volcanic eruption events [Pauling . (2023)].

Employing eruptions in the medium to super-volcano size enhances the signal-
to-noise ratio without necessitating an extensive and computationally
expensive ensemble, and as such, is a tempting way to mimic a large ensemble
of smaller volcanic eruptions. However, the AOD, RF, and temperature
signatures are not necessarily a simple scaling of that of smaller volcanic
eruptions. Previous studies have simulated super-volcanoes using AOD as the
input forcing, where the AOD was that of Mt. Pinatubo scaled by a factor of
one hundred [Jones . (2005)]. This approach may yield incorrect results, both
because the peak of the AOD may be too small or too big, but also because the
evolution of the AOD could be inappropriate. Likewise, a different AOD
evolution may be found from similar size eruptions, but at different
latitudes. To investigate this issue, our simulations are based on four levels
of injected SO2 covering three orders of magnitude and the inclusion of one
high latitude eruption of the second largest injected SO2 case.

We conducted ensemble simulations of volcanic eruptions in the Community Earth
System Model Version 2 (CESM2) coupled with the Whole Atmosphere Community
Climate Model Version 6 (WACCM6). The ensembles span four different levels of
injected SO2: $26\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$,
$400\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$,
$1629\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$ and
$3000\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$. Details
regarding the experimental setup are provided in section 2. Our findings
reveal non-linear RF to AOD dependencies for medium to super-volcano size
eruptions. Additionally, we observe a time-dependent variation in the RF to
AOD ratio, detailed in section 3 and discussed in section 4. Furthermore, our
data, along with insights from previous studies, suggest that the RF
dependency on injected SO2 identified by niemeier2015 acts as a lower
boundary. Our conclusions are presented in section 5.

##  2 Method

###  2.1 Model

We use the CESM2 [Danabasoglu . (2020)] in conjunction with the WACCM6
[Gettelman . (2019)] and the fully dynamical ocean component Parallel Ocean
Program version 2 (POP2) [R. Smith . (2010), Danabasoglu . (2020)]. The
atmosphere model was run at a nominal $2\text{\,}\mathrm{\SIUnitSymbolDegree}$
resolution with $70$ vertical levels in the middle atmosphere (MA)
configuration.

The WACCM6 version employed in the MA configuration uses the three mode
version of the Modal Aerosol Module (MAM3) [Gettelman . (2019)], a simplified
and computationally efficient default setting within the Community Atmosphere
Model version 5 (CAM5) [Liu . (2016)], as described in liu2012. The MAM3 was
developed from MAM7 and features the modes Aitken, accumulation, and coarse
[Liu . (2016)].

###  2.2 Simulations

Appendix A provides a description of the simulation setup and utilised output
variables. Table 1 summarises the simulations, encompassing four SO2 injection
magnitudes and up to four seasons: 15 February, 15 May, 15 August, and 15
November. The magnitudes vary over three orders of magnitude:
$26\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$,
$400\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$,
$1629\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$, and
$3000\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$.

The smallest eruption case, C2W$\downarrow$, is similar in magnitude as
compared to events like Mt. Pinatubo [<]$\sim
10$–$20\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$;¿[]timmreck2018
and Mt. Tambora
[<]$\sim$56.2\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$$;¿[]zanchettin2016.
The intermediate case, C2W$-$, resembles the magnitude of the Samalas eruption
in 1257
[<]$\sim{144}$–$170\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$;¿[]vidal2016,
while the second largest and largest eruption cases, C2W$\uparrow$ and
C2W$\uparrow\uparrow$, is in the likely range of the Young Toba Tuff (YTT)
eruption occurring about $72\,000\text{\,}\mathrm{yr}$ ago
[<]$100$–$10\,000\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$;¿[]jones2005.
All eruptions were situated at the equator
($0\text{\,}\mathrm{\SIUnitSymbolDegree}\mathrm{N}$,
$1\text{\,}\mathrm{\SIUnitSymbolDegree}\mathrm{E}$) with SO2 injected between
$18\text{\,}\mathrm{km}$ and $20\text{\,}\mathrm{km}$ altitude. Collectively,
the four tropical eruption cases C2W$\downarrow$, C2W$-$, C2W$\uparrow$, and
C2W$\uparrow\uparrow$ are referred to as C2WTrop. An additional high-latitude
eruption ensemble, labelled C2WN$\uparrow$, of the same injected SO2 magnitude
as C2W$\uparrow$ was simulated at
$56\text{\,}\mathrm{\SIUnitSymbolDegree}\mathrm{N}$,
$287.7\text{\,}\mathrm{\SIUnitSymbolDegree}\mathrm{E}$ with a six-month
separation (15 February and 15 August).

Table 1: Simulations done with the CESM2a Ensemble name | $\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$ | Lat [$\mathrm{\SIUnitSymbolDegree}\mathrm{N}$] | Lon [$\mathrm{\SIUnitSymbolDegree}\mathrm{E}$] | Alt [$\mathrm{km}$] | Eruption months
---|---|---|---|---|---
C2W$\uparrow\uparrow$ | $3000$ | $\hphantom{1}0$ | $\hphantom{28}1\hphantom{.7}$ |  $18$–$20$ |  Feb,May,Aug,Nov
C2WN$\uparrow$ | $1629$ | $56$ | $287.7$ |  $18$–$20$ | Feb,May,Aug,Nov
C2W$\uparrow$ | $1629$ | $\hphantom{1}0$ | $\hphantom{28}1\hphantom{.7}$ |  $18$–$20$ | Feb,May,Aug,Nov
C2W$-$ | $\hphantom{1}400$ | $\hphantom{1}0$ | $\hphantom{28}1\hphantom{.7}$ |  $18$–$20$ | Feb,May,Aug,Nov
C2W$\downarrow$ | $\hphantom{14}26$ | $\hphantom{1}0$ | $\hphantom{28}1\hphantom{.7}$ |  $18$–$20$ | Feb,May,Aug,Nov
aThe ensembles C2WN$\uparrow$ and C2W$\uparrow$ have the same eruption
magnitude, but while C2W$\uparrow$ is located at the equator, C2WN$\uparrow$
is located at a high northern latitude. C2W$\uparrow\uparrow$, C2W$-$ and
C2W$\downarrow$ are located at the equator, but with different magnitudes
compared to C2W$\uparrow$. The three smallest tropical ensembles have four
members, indicated by the number of eruption months, while the northern
latitude and the extra large super-volcano ensemble consists of two members.

##  3 Results

###  3.1 Analysis of the time series

Figure 1 presents time series of global mean AOD, RF, and surface air
temperature. The black lines represent the medians across the ensembles, while
shading indicates the 5th to 95th percentiles. The four distinct forcing
magnitudes (C2W$\downarrow$, C2W$-$, C2W$\uparrow$, and C2W$\uparrow\uparrow$)
outlined in table 1 have been used. The time series in Fig. 1 are normalised
by setting the peak value to unity, defined based on the peak of a fit from a
Savitzky-Golay filter of 3rd order and a one-year window length [Savitzky
Golay (1964)].

A notable feature across the subfigures of Fig. 1 is the peak occurrence of
the C2W$\downarrow$ case compared to the larger eruption cases. The peak of
C2W$\downarrow$ arrives earlier for both AOD (Fig. 1a) and temperature (Fig.
1c), while the RF time series in Fig. 1b are all indistinguishable. Cases
C2W$-$, C2W$\uparrow$, and C2W$\uparrow\uparrow$ are indistinguishable in
their temperature development, and while C2W$\downarrow$ peaks at an earlier
time, it decays similarly to the other cases. Interestingly, the same
development between C2W$-$ and C2W$\uparrow$ is not found in the AOD time
series. C2W$\downarrow$ peaks at an earlier time, but also spends more time
around the peak and as such decays at a later time post-eruption. Likewise,
C2W$-$ has a faster rise and slower decay compared to C2W$\uparrow$, but where
both peak at a similar time. C2W$\uparrow$ and C2W$\uparrow\uparrow$ have
similar AOD developments, but where C2W$\uparrow\uparrow$ show a slightly
faster decay from the peak.

The timescale of the perturbation of AOD and RF is shorter than that of the
temperature. While the AOD and RF time series return to their equilibrium
state within roughly three years, the temperature time series remain heavily
perturbed three years post-eruption. Even when running the simulations for 20
years post-eruption, the temperature time series are still decaying.

Figure 1: AOD (a), RF (b) and temperature response (c) time series to the four
tropical volcanic eruption cases, C2W$\downarrow$, C2W$-$, C2W$\uparrow$, and
C2W$\uparrow\uparrow$. The time series have been normalised to have peak
values at unity, where $C$ is the normalisation constant. Black lines indicate
the median across the ensembles, while shading marks the 5th and 95th
percentiles.

###  3.2 RF dependency on AOD

We next focus on the development of the AOD and RF time series relative to
each other. Similar comparisons were conducted in ¡¿[their Fig. 4]gregory2016
and ¡¿[their Fig. 1]marshall2020, with RF plotted against AOD. Figure 2
displays annual mean values from the five simulation cases in table 1; the
small eruption case (C2W$\downarrow$) as blue downward-pointing triangles, the
intermediate eruption case (C2W$-$) as orange thick diamonds, the large
tropical eruption case (C2W$\uparrow$) as green upward-pointing triangles, the
extra large eruption case (C2W$\uparrow\uparrow$) as small pink upward-
pointing carets, and the large northern hemisphere eruption case
(C2WN$\uparrow$) as brown upward-pointing three-branched twigs. Also shown are
the data from ¡¿[Fig. 4, black crosses from HadCM3 sstPiHistVol]gregory2016 as
grey crosses labelled G16 (described in Appendix B, section B.3).
Additionally, the estimated peak values from the Mt. Pinatubo and Mt. Tambora
eruptions are plotted as a black star and plus, while the peak from the
jones2005 simulation is shown as a pink square labelled J05. Finally, red
circles represent the peak values obtained from the C2W eruption cases. The
straight lines are the same as shown by gregory2016. The full data range is
shown in Fig. 2a while Fig. 2b highlights a narrow range, focusing on the
C2W$\downarrow$ case.

The annual mean data from the Pinatubo-like C2W$\downarrow$ case in Fig. 2b
have RF values as a function of AOD that follow almost the same constant slope
as the G16 data. However, in Fig. 2a we observe that the stronger eruptions
lead to dissimilar responses in AOD and RF, where C2W$-$ seems to follow close
to a $-10$ slope and C2W$\uparrow$ is closer to a $-5$ slope. The peak values
(red circles) suggest a non-linear dependence, while within each eruption
strength (same colour) the annual mean values fall relatively close to a
straight line.

Figure 2: RF as a function of AOD, yearly means. Data from the five
simulations listed in table 1 (C2W$\downarrow$, C2W$-$, C2W$\uparrow$,
C2WN$\uparrow$, and C2W$\uparrow\uparrow$) are shown along with the data from
the HadCM3 sstPiHistVol simulation by gregory2016 (grey crosses, G16). Also
shown are the estimated peak values of the Mt. Pinatubo (black star) and Mt.
Tambora (black plus) eruptions. The peak values from the C2W simulations are
shown as red circles. Additionally in (a) the simulated super-volcano of
jones2005 (pink square) is shown. All peak values (as opposed to annual means)
have an asterisk ($\ast{}$) in their label. The grey lines are the same
regression fits as in ¡¿[Fig. 4]gregory2016, where the solid line is the fit
to G16. (b): Zooming in on the smallest AOD values.

To investigate the time dependence of the ratio between RF and AOD, we present
seasonal means of this ratio in Fig. 3. The plot shows the eruption cases
given in table 1, as well as the tropical eruptions from marshall2020dataset
($6$ of $82$ eruptions), labelled M20 and described in Appendix B, section
B.2. The C2W$\uparrow$ case is similar to C2W$\uparrow\uparrow$ as indicated
in table 2, but is not shown in the plot to better highlight C2WN$\uparrow$.
In Fig. 3a, lines are linear regression fits to the seasonal means across all
ensemble members, summarised in table 2. Shaded regions are the standard
deviation around the seasonal means. A similar shading is plotted in Fig. 3b,
but where the regression fits have been omitted for clarity. As the AOD and RF
time series start from zero, the ratio from the first season is not included.
Likewise, after three years both time series are almost fully equilibrated
(Fig. 1a,b). The data is further divided into two periods; a pre-peak period
where the peak of both the AOD and the RF is included (consisting of the first
post-eruption year), and a post-peak period for the decaying part (consisting
of the second and third post-eruption years).

Although the ratio changes across the eruption magnitudes, we find that all
the tropical cases follow a positive slope during the pre-peak period, as seen
in Fig. 3a and described in table 2. The northern latitude case in
C2WN$\uparrow$ shows a much flatter slope compared to C2WTrop and M20. The
distinction between the slopes from the tropical and non-tropical cases is
perhaps more clear in Fig. 3b and corresponding rows in table 2. Again,
C2WN$\uparrow$ shows an almost flat slope compared to the tropical cases.
During the post-peak period, more noise is introduced, but a weak tendency of
negative slopes is found among the tropical cases, as well as in the
C2WN$\uparrow$ case up to the last season where the noise is also the largest.

Table 2: Slope and standard deviation for the data in Fig. 3a Figure | Ensemble name | Pre-peak | Post-peak
---|---|---|---
| C2WN$\uparrow$ | $0.45\pm 1.15$ | $1.51\pm 1.45$
| C2W$\uparrow\uparrow$ | $3.38\pm 0.97$ | $-2.74\pm 0.77$
3a | C2W$\uparrow$ | $3.85\pm 0.52$ | $-3.29\pm 0.60$
C2W$-$ | $4.36\pm 0.82$ | $-3.37\pm 0.59$
| C2W$\downarrow$ | $3.64\pm 2.41$ | $-1.41\pm 3.25$
| M20 | $6.34\pm 1.77$ | $-0.36\pm 1.33$
| C2WN$\uparrow$ | $0.08\pm 0.20$ | $0.27\pm 0.26$
| C2W$\uparrow\uparrow$ | $0.86\pm 0.25$ | $-0.70\pm 0.19$
3b | C2W$\uparrow$ | $0.75\pm 0.10$ | $-0.64\pm 0.12$
C2W$-$ | $0.43\pm 0.08$ | $-0.34\pm 0.06$
| C2W$\downarrow$ | $0.18\pm 0.12$ | $-0.07\pm 0.16$
| M20 | $0.33\pm 0.07$ | $-0.02\pm 0.08$
aThe regression fits in the top half of the table are for Fig. 3a, while the
bottom half is for Fig. 3b. The columns “pre-peak” and “post-peak” refer to
the two periods as shown in Fig. 3. The ensembles are the same as those given
in table 1, in addition to the $6$ tropical eruptions from the $82$ member
ensemble in marshall2020.

¡¿[their Fig. 1c,d]marshall2020 present results that demonstrate a time-
dependent relationship in the conversion between AOD and RF. They obtain an RF
to AOD ratio with a negative slope when comparing the first post-eruption year
to the second and third. As such, marshall2020 find that, on average, the
aerosol forcing efficiency increases during the first two to three post-
eruption years. This phenomenon is explained by marshall2020 as the aerosols
initially being spatially confined to the hemisphere where the eruption
occurred. Subsequently, during the second and third years, they spread
globally, resulting in a higher global-mean albedo per AOD and consequently a
stronger RF per AOD ratio with time. However, as noted above, a decrease in
aerosol forcing efficiency is found when analysing the M20 data with seasonal
resolution during the pre-peak period (first year post-eruption) while
constraining the ensemble to only include eruptions within $-10$ to
$10\text{\,}\mathrm{\SIUnitSymbolDegree}\mathrm{N}$. The post-peak period
shows an increasing aerosol forcing efficiency, and during the full first
three post-eruption years (pre-peak and post-peak), both the tropical subset
and the full M20 data yield an increasing efficiency, as expected. Likewise,
the first three post-eruption years of the C2W$-$, C2W$\uparrow\uparrow$, and
C2WN$\uparrow$ cases show a weak negative slope and thus an increasing
efficiency, while C2W$\downarrow$ shows an elevated post-peak ratio as seen in
Fig. 3b.

We also note that while the aerosol forcing efficiency is decreasing for
tropical M20 data in the pre-peak period, the full dataset shows increasing
efficiency. This is in line with what we find from C2WN$\uparrow$, which is
the only eruption case that does not show a clear aerosol forcing efficiency
decrease during the pre-peak period.

Figure 3: (a): The ratio of RF to AOD, with time-after-eruption on the
horizontal axis. Straight lines indicate linear regression fits and are
described in table 2, while shaded regions are the standard deviation across
the ensembles for each season. Regression fits and shadings are made for the
pre-peak and post-peak periods. (b): Same as in (a), but where the underlying
AOD and RF time series have been scaled to have peak values at unity. Shown
are data from table 1 along with tropical eruptions from M20.

###  3.3 Parameter scan

In Fig. 4, we compare the peak values of all investigated CESM2 output
parameters against each other as well as to injected SO2. For our tropical
cases (C2WTrop), we observe in Fig. 4a an almost linear yet notably weakening
relationship between AOD peak values and injected SO2. The latitude also plays
a role in the magnitude of the AOD perturbation, evident from C2WN$\uparrow$.
This weak yet notable latitude dependence aligns with findings by
marshall2019, indicating that $72\text{\,}\mathrm{\char 37\relax}$ of the AOD
variance can be attributed to injected SO2, while latitude accounts for only
$16\text{\,}\mathrm{\char 37\relax}$ of the variance. Peak values from their
data (82 simulations) plotted as red thin diamonds display a similar pattern,
with AOD exhibiting close to linear dependence on injected SO2, but with
latitude introducing a spread in AOD. Peak values from Mt. Pinatubo (P) and
Mt. Tambora (T) are shown for reference, along with peak values from jones2005
labelled J05 and timmreck2010 labelled T10. The J05 is a simulation of a
super-volcano based on a $100$ times scaling of the AOD from Mt. Pinatubo,
while T10 is a simulation of the YTT eruption based on SO2 injections.

In Fig. 4b, RF plotted against injected SO2 (with the absolute value of RF on
the $y$-axis) indicates a substantial damping effect on RF as injected SO2
increases for the C2W data, in agreement with results from ottobliesner2016,
labelled OB16. The OB16 data come from a $2500$ year long simulation using
historic volcanoes as the only external forcing. The analysis details of OB16
can be found in Appendix B, section B.1. Despite the model complexity
difference, ottobliesner2016’s simulations using Community Earth System Model
version 1 (CESM1) with a low-top atmosphere (CAM5) produce RFs comparable to
our findings.

Figure 4: (a) AOD, (b) RF, and (c) temperature anomaly as a function of
injected SO2. (d) RF and (e) temperature anomaly as a function of AOD. (f)
Temperature anomaly as a function of RF. Blue diamonds labelled C2WTrop
represent tropical cases (C2W$\downarrow$, C2W$-$, C2W$\uparrow$,
C2W$\uparrow\uparrow$), the brown three-branched twig signifies the
C2WN$\uparrow$ case, and green downward triangles denote OB16 data from
ottobliesner2016. The red thin diamonds labelled M20 display the
marshall2020dataset data. Black star and plus indicate Mt. Pinatubo and Mt.
Tambora estimates based on observations. The pink square labelled J05 refers
to the one-hundred times Mt. Pinatubo super-volcano from jones2005, and the
pink disk labelled T10 represents the YTT super-volcano from timmreck2010. The
pink dashed line labelled N15 is from niemeier2015, indicating the function in
Eq. 1.

niemeier2015 conducted simulations of continuous sulphur injections up to
$200\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}\mathrm{yr}^{-1}$
in the ECHAM5’s middle atmosphere version [Giorgetta . (2006)] with aerosol
microphysics from HAM [Stier . (2005)]. They observed an RF dependence on SO2
injection rate following an inverse exponential, which converges to
$-65\text{\,}\mathrm{W}\mathrm{m}^{-2}$, depicted in Fig. 4b as the stippled
pink line labelled N15 and given as;

$\Delta
R_{\mathrm{TOA}}=-$65\text{\,}\mathrm{W}\mathrm{m}^{-2}$\mathrm{e}^{-{\left(\frac{$2246\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\mathrm{S}\mathrm{)}\mathrm{y}\mathrm{r}^{-1}$}{x}\right)}^{0.23}}.$
(1)

Both our simulations and OB16 exhibit a notably faster increase than this
exponential relationship. The results by N15, on which Eq. 1 is based, are all
averages over at least three years of steady sulphur burdens, substantially
longer than the time it takes for RF to reach peak values after an eruption.
Combined with their lack of a full chemistry model [Niemeier Timmreck (2015)],
a direct comparison between Eq. 1 to peak RF values (occurring about one year
post-eruption) may not reflect the same chemical and physical processes. In
Eq. 1, $x$ represents S, while the axis shows values of SO2, thus halving of
the SO2 values on the axis gives the appropriate shape of Eq. 1 as a function
of S.

With these caveats in mind, we observe that T10’s results closely align with
the function described in Eq. 1. Starting with an initial input of
$850\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{S}\mathrm{)}$ (equivalent to
$1700\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$, representing
the YTT eruption), their estimated AOD led to a peak RF of
$-18\text{\,}\mathrm{W}\mathrm{m}^{-2}$, depicted as a pink filled circle in
Fig. 4b. The results from T10 came from a simulation using the MPI-ESM climate
model, driven by AOD data from the HAM aerosol model. This alignment likely
stems from the utilization of the same aerosol microphysical model in both
timmreck2010 and niemeier2015, as well as the application of similar climate
models, MPI-ESM and ECHAM5, respectively. The relationship between climate
model families and their implications are further described in Appendix C.
Notably, the peak values from M20 fit well within an upper boundary defined by
C2WTrop and OB16, and a lower boundary defined by Eq. 1. Eruptions closer to
the equator within M20 align with data points near the upper boundary, whereas
eruptions at more extreme latitudes tend to yield weaker peak RF values,
closer to the lower boundary. Importantly, none of the eruption simulations
shown in Fig. 4b exceeded the upper threshold of
$-65\text{\,}\mathrm{W}\mathrm{m}^{-2}$ as suggested in Eq. 1.

Figure 4c illustrates the response of temperature against injected SO2. The
increase in temperature response with injected SO2 decreases for higher
injected SO2, showing a similar relationship between C2WTrop, C2WN$\uparrow$,
and OB16. Notably, T10 and J05 exhibit respectively much weaker and much
stronger temperature responses to injected SO2 than C2WTrop. T10 has a maximum
temperature anomaly of only $-3.5\text{\,}\mathrm{K}$ for their
$1700\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$ eruption,
while J05 records a substantially larger maximum temperature anomaly of
$-10.7\text{\,}\mathrm{K}$. Since the M20 experiment was conducted with
prescribed sea-surface temperatures [Marshall . (2020)], preventing the
temperature from being fully perturbed, we do not focus on the M20 data in the
temperature plots but include them for completeness.

In Fig. 4d, we revisit the relationship between RF and AOD, focusing on peak
values rather than annual and seasonal averages. As previously discussed, the
RF to AOD ratio displays weaker slopes than previous studies, with the C2W
peak values not conforming to a linear trend. The relationship between RF and
AOD suggests potential substantial dependencies on the model and its input
parameters, such as latitude, but most notably to an inherent non-linear RF
dependence on AOD. Both the G16 data in Fig. 2 and the J05 data originate from
the same climate model. Similarly to what we find from the C2W data, the ratio
is much stronger for small eruptions in the industrial era (G16) compared to
the super-volcano eruption (J05).

In Fig. 4e, we again find that the response of the C2WTrop data decreases with
injected SO2, this time in temperature anomaly. Additionally, both the
C2WN$\uparrow$ and the J05 cases align well with C2WTrop, with the T10 case
following a similar dependence.

Finally, in Fig. 4f, we compare the temperature and RF responses. Both C2WTrop
and OB16 show a near-linear relationship between temperature and RF. The
C2WTrop data indicate a steeper slope, implying stronger temperature
perturbations as compared to OB16. However, potential biases exist in the
values from the analysis of OB16, as outlined in Appendix B, section B.1.
This, along with considerable noise, results in the analysis of OB16
temperature anomalies being less reliable. As in Fig. 4e, the C2WN$\uparrow$
case along with both the T10 and J05 cases closely follow the temperature to
RF dependence of C2WTrop.

The almost linear relationship between AOD and injected SO2 for the C2WTrop
data in Fig. 4a suggests a comparable trend for RF versus injected SO2 in Fig.
4b, as seen for RF versus AOD in Fig. 4d. For the same reason, we expect Fig.
4e to show a similar pattern for C2WTrop as observed in Fig. 4c.

This relationship, along with the functional relationships between all other
parameters shown in Fig. 4, are illustrated in Fig. 5. There, we show that
from assuming a linear dependency of AOD on injected SO2 ($ax+b$), and of
temperature on RF ($cx+d$), we must have that $f$, $g$, $h$, and $k$ all have
the same functional form, where $f:\ce{SO2}\to\mathrm{RF}$,
$g:\mathrm{AOD}\to\mathrm{T}$, $h:\ce{SO2}\to\mathrm{T}$, and
$k:\mathrm{AOD}\to\mathrm{RF}$. From this, we deduce that $f(x)=k(ax+b)$ and
$h(x)=f(cx+d)=g(ax+b)$, and finally that $h(x)=k(acx+ad+b)$, concluding that
$f$, $g$, $h$, and $k$ have the same functional form.

AODTRFSO2$g$$cx+d$$f$$ax+b$$h$$k$ Figure 5: Diagram describing the functional
relationships of the parameters shown in Fig. 4.

###  3.4 Climate sensitivity estimate

As previously mentioned, the J05 experiment is similar to C2W$\uparrow$
concerning RF values, yet differ in both AOD and temperature. At the same time
J05 is similar to C2W$\uparrow\uparrow$ in AOD and RF. To investigate this
discrepancy, we here conduct a comparison between the J05 climate feedback
parameter $\alpha$ (where $s=1/\alpha$ is the climate sensitivity parameter)
with our climate resistance, denoted as $\rho$, and the transient climate
response parameter (TCRP) $1/\rho$ (where
$\mathrm{TCS}=F_{2\times\ce{CO2}}\times\mathrm{TCRP}$ is the transient climate
sensitivity and $F_{2\times\ce{CO2}}$ is the forcing due to a doubling of pre-
industrial CO2 concentration). As the forcing of volcanic eruptions typically
lasts for about a year, a duration too short for $F=\rho T$ to remain valid,
an alternative approach using a time-integral form introduced by merlis2014 is
applied:

$\int_{0}^{\tau}F\mathrm{d}t=\rho\int_{0}^{\tau}T\mathrm{d}t$ (2)
$\rho=\frac{\int_{0}^{\tau}F\mathrm{d}t}{\int_{0}^{\tau}T\mathrm{d}t}.$ (3)

If the upper bound of the integral, $\tau$, is sufficiently large, so that the
upper ocean heat content is the same at $t=0$ and $t=\tau$ (merlis2014 used
$\tau=$15\text{\,}\mathrm{yr}$$), this approach agrees with $F=\rho T$ for
long-term forcing [Gregory . (2016)]. Additionally, we note that the climate
resistance and the climate feedback parameter are associated with the ocean
heat uptake efficiency ($\kappa$) through $\rho=\alpha+\kappa$ [Gregory .
(2016)].

The climate feedback parameter estimated by jones2005 is
$\alpha\simeq$4\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{K}^{-1}$$, exceeding
twice the value obtained by gregory2016 in their simulations of Mt. Pinatubo
using the same HadCM3 climate model. We determine the climate resistance using
the integral-form computation outlined in Eq. 3 and adopting
$\tau=$20\text{\,}\mathrm{yr}$$. The estimated climate resistance from the
three tropical simulation cases (with four in each ensemble) converges to
$\rho=$2.5(4)\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{K}^{-1}$$, and TCRP
values of $1/\rho=$0.41(5)\text{\,}\mathrm{K}\mathrm{W}^{-1}\mathrm{m}^{2}$$,
as reported in table 3, and is therefore assumed to be a good estimate of
$\alpha$.

Importantly, our estimate agrees well with G16, while the J05 estimate of
$\alpha\simeq$4\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{K}^{-1}$$ is still
notably higher. Since the temperature perturbation obtained by J05 was larger
than in any of our CESM2 cases, it indicates that the forcing used by J05 must
be stronger. The peak value of the J05 RF is similar to the C2W$\uparrow$
case, and as such, the overall stronger forcing must originate from the
development of the forcing time series rather than the peak value.

Table 3: Estimated climate resistance and TCRPa Simulation type | $\rho[$\mathrm{W}\mathrm{m}^{-2}\mathrm{K}^{-1}$]$ | $1/\rho$
---|---|---
C2W$\uparrow$ | $2.21(5)$ | $0.45(1)$
C2W$-$ | $2.51(6)$ | $0.40(1)$
C2W$\downarrow$ | $2.9(6)$ | $0.36(7)$
Total | $2.5(4)$ | $0.41(5)$
aEstimates are based on ensembles with four members and
$\tau=$20\text{\,}\mathrm{yr}$$ using Eq. 3.

##  4 Discussion

Figures 2, 3, and 4d demonstrate that as the AOD exceeds approximately $1.0$,
the linear RF dependence of approximately
$-20\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{AOD}^{-1}$ no longer holds. The
sublinear increase in RF with injected SO2 in Fig. 4b for large eruptions is
consistent with previous results from simulations using similar climate models
of smaller historic eruptions (G16) and of super-volcanoes (J05). Such a
change in ratio has been attributed to larger eruptions, injecting more SO2,
leading to larger aerosols, and hence less effective radiation scattering,
thereby reducing the RF for the same injected SO2 [English . (2013), Timmreck
. (2010), Timmreck . (2018)].

The non-linear relationship between peak RF and AOD values is a strong
signature in both Figs. 2 and 3. Across eruptions of the same strength, the
ratio stays relatively constant, leading to a close to
$-10\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{AOD}^{-1}$ slope for C2W$-$ and
a $-5\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{AOD}^{-1}$ slope for
C2W$\uparrow$ and C2W$\uparrow\uparrow$. Still, a non-linear development in
the RF to AOD ratio is found across all tropical eruptions. Similar to the
results of marshall2020, we find in C2W$-$, C2W$\uparrow$, and
C2W$\uparrow\uparrow$ that the post-peak period (second and third post-
eruption years) has a stronger aerosol forcing efficiency compared to the pre-
peak period (first post-eruption year). The post-peak period of
C2W$\downarrow$ is elevated as compared to the pre-peak period, resulting in a
decreasing aerosol forcing efficiency from the first to the second and third
post-eruption years, in contrast to the other tropical eruptions.

Focusing on the pre-peak period, we find tropical eruptions to differ from
eruptions at high latitudes. During the pre-peak period, all tropical
eruptions show a decreasing aerosol forcing efficiency, while no significant
change in the RF to AOD ratio is found from the C2WN$\uparrow$ case. The full
M20 dataset indicates an increasing aerosol forcing efficiency also during the
pre-peak period, contrasting the decreasing efficiency found from their
tropical eruptions and supporting the latitudinal dependence we find with
C2WN$\uparrow$. While we find a linear relationship to be a useful
approximation of RF dependence on AOD for eruptions similar to or smaller than
Mt. Pinatubo, additional factors must be considered for larger eruptions.
These factors, such as OH scarcity and aerosol growth, influence reflectance
and their gravitational pull, substantially impacting both AOD and RF
evolution, is highlighted by timmreck2010. The large difference in ratio found
when comparing eruption magnitudes suggests that injected SO2 is crucial when
estimating the time-average of the RF to AOD ratio. However, latitude and, in
particular, aerosol dispersion are more influential in determining the post-
eruption evolution of the ratio, particularly during the pre-peak period.

We find that the suggested upper threshold from Eq. 1 is not violated by any
eruption simulation, and most notably that the temperature peak value follow
the RF trend in reaching a limiting value. The C2WTrop cases follow a close to
linear temperature dependency on RF, with the J05, T10, and C2WN$\uparrow$
aligning close to the same slope. The linear relationship between temperature
and RF is the strongest dependence found between the parameters in Fig. 4, and
a strong signature across both eruption magnitudes and latitudes, but also
across highly different climate models. Thus, from a maximum RF of
$-65\text{\,}\mathrm{W}\mathrm{m}^{-2}$, we expect temperature anomalies to
reach at most $\sim$-12\text{\,}\mathrm{K}$$.

The biggest spread in the data shown in Fig. 4 is found when relating injected
SO2 to any of the three output parameters. As the amount of injected SO2
increases, both AOD, RF, and temperature across models have a big spread. The
AOD to injected SO2 relationship is consistent within similar models, even
when comparing simulations of volcanic eruptions [Timmreck . (2010)] and
continuous injection of SO2 [Niemeier Timmreck (2015)], but has a wide spread
at high values of injected SO2 across model families (Figs. 4a,b,c).
Comparatively, the RF (Fig. 4d) and temperature (Fig. 4e) as a function of
AOD, as well as temperature as a function of RF (Fig. 4f), demonstrate a
smaller spread across models. marshall2019, marshall2020, marshall2021 use a
code with seven log-normal modes to simulate aerosol mass and number
concentrations, along with an atmosphere-only configuration of the UM-UKCA
with prescribed sea-surface temperatures and sea-ice extent [Marshall .
(2019)]. This approach is in contrast with CESM2, operating as an Earth System
Model, but with a simpler aerosol chemistry model in the MAM3. The family of
models to which M20 is based is different from that of C2W and OB16, and also
different from the T10 and N15, as described in Appendix C. Based on Fig. 4,
we find the model family to be pivotal in determining the estimated AOD and RF
magnitudes from injected SO2, whereas the various models generally demonstrate
more consistency in representing RF from AOD.

timmreck2010 highlights that for sufficiently large eruptions, OH radicals are
too scarce, which limits SO2 oxidation. The AOD peak in the YTT simulation of
T10 occurs six months after Mt. Pinatubo’s peak. This aligns with our results,
as illustrated in Fig. 1a, where C2W$\downarrow$ shows an earlier AOD peak
compared to C2W$-$, C2W$\uparrow$, and C2W$\uparrow\uparrow$. While the peak
RF value of T10 occurs $7$–$8$ months post-eruption, similar to C2W, the J05
peak anomaly occurs one year post-eruption. Additionally, as jones2005 obtains
a climate feedback parameter larger than both what gregory2016 found for the
same climate model and larger than the climate resistance obtained here from
C2W, we conclude that such a simple approach of scaling the AOD of smaller
eruptions to represent larger eruptions is insufficient. Moreover, having a
small ensemble of large eruptions to represent smaller eruptions is also
insufficient when simulating from injected SO2, as both AOD and temperature
evolution are found to develop differently.

##  5 Summary and conclusions

We consider five medium to super-volcano sized eruption ensembles and compare
them to previously reported results. We find the commonly adopted RF
dependence on AOD of
$\sim$-20\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{AOD}^{-1}$$ to be
representative for Mt. Pinatubo-sized eruptions. Larger eruptions, with one to
two orders of magnitude larger injections of SO2, are found to have an RF
dependence on AOD closer to
$\sim$-10\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{AOD}^{-1}$$ and
$\sim$-5\text{\,}\mathrm{W}\mathrm{m}^{-2}\mathrm{AOD}^{-1}$$. A shallower
slope for larger eruptions is also consistent with peak values from previous
studies of super-volcanoes.

The time-after-eruption dependence of the ratio between RF and AOD is found to
weaken with time, resulting in a decreasing aerosol forcing efficiency in the
pre-peak period. The effect is found across all eruption sizes, but only the
tropical cases show a clear trend. The high-latitude case displays an almost
constant efficiency with time. These results agree with a reanalysis of the
tropical data in marshall2020dataset. Thus, these findings provide strong
supporting evidence that latitude is generally significant in determining the
aerosol forcing efficiency, particularly as a function of time-after-eruption.
These findings emphasise the complexity of volcanic impacts on climate,
demonstrating significant differences in climatic response depending on
eruption magnitude and latitude.

We find that the AOD peak arrives later for larger eruptions than for smaller
ones, and also that larger eruptions produce a sharper peak in the AOD time
series. The RF time series are similar across all eruption sizes, and while
the smallest eruption experiences a faster temperature decay, the larger
eruptions produce time series indistinguishable in development for both RF and
temperature. Thus, a simple scaling of the AOD or temperature time series from
a smaller eruption is insufficient in representing that of larger volcanic
eruptions.

Considering injected SO2 and the peak values of AOD and RF, a large spread is
found across model families in Fig. 4. Improving the consistency between model
families in how the chemistry and physics of SO2 and H2SO4 are represented is
an important step in enhancing the accuracy of simulated volcanic eruptions’
influence on climate by models. More simulations of larger volcanic eruptions
with injected SO2 greater than
$200\text{\,}\mathrm{T}\mathrm{g}\mathrm{(}\mathrm{SO2}\mathrm{)}$ would
provide useful information for a more precise determination of the RF to AOD
ratio in the non-linear regime. This would also serve as a useful test to
check if a comparison between SO2 injection events and continuous SO2
injection is reasonable. Introducing a spread in latitude similar to the
marshall2020dataset dataset would allow for better comparison between
eruptions across all latitudes and the suggested lower limit following Eq. 1,
describing a situation of aerosol saturation.

##  Appendix A Simulation set up and output

Input files used in the simulations were created by modifying the file
available at
http://svn.code.sf.net/p/codescripts/code/trunk/ncl/emission/createVolcEruptV3.ncl,
using a Python package available on GitHub at
https://github.com/engeir/volcano-cooking or through the Python Package Index
(PyPI). The package is available both as a library and a Command Line
Interface (CLI), and is used to create volcanic eruptions with a specified
amount of SO2 that is injected over six hours at a given latitude, longitude,
and altitude. All volcanic SO2 files are created from a shell script by
setting the eruption details in a JSON file that is read by the volcano-
cooking CLI at a fixed version, ensuring a reproducible experiment setup.

We are using the coupled model version BWma1850 component setup to run the
CESM2, and an accompanying fixed sea-surface temperature version, fSST1850, to
obtain estimates of the RF. The applied fSST1850 is not from a standardised
component setup but is instead explicitly specified as
1850_CAM60%WCCM_CLM50%BGC-
CROP_CICE%PRES_DOCN%DOM_MOSART_CISM2%NOEVOLVE_SWAV_TEST%. The component setup
BWma1850 and fSST1850 differ in that the latter uses a prescribed sea-ice
(CICE -> CICE%PRES), a prescribed data ocean (POP2%ECO%DEP -> DOCN%DOM) and a
stub wave component instead of the full Wave Watch version 3 (WW3 -> SWAV).

The important input data used in the model simulations are injected SO2 in
units of teragrams ($\mathrm{T}\mathrm{g}\mathrm{(}\ce{SO2}\mathrm{)}$), used
to simulate volcanic eruptions. RF is calculated as the combined (short wave
and long wave) all-sky TOA energy imbalance, where the CESM2 provide the
output variables “net solar flux at the top of the model” (FSNT) and “net
longwave flux at the top of the model” (FLNT). Thus,
$\mathrm{RF_{*}}=\mathrm{FSNT}-\mathrm{FLNT}$, and taking the difference
between volcanic forcing simulations and a control simulation gives the final
estimate of RF ($\mathrm{RF}=\mathrm{RF_{VOLC}}-\mathrm{RF_{CONTROL}}$)
[Marshall . (2020)]. The RF calculation is based on fSST1850, hence this
outline specifically describes how to calculate ERF as opposed to IRF, which
instead is the difference between the ERF and the sum of all rapid atmospheric
adjustments [Marshall . (2020), CJ. Smith . (2018)]. The AOD is obtained from
the output variable “stratospheric aerosol optical depth 550 nm day nigth”
(AODVISstdn), while global temperature is saved by CESM2 to the variable
“reference height temperature” (TREFHT). The analysis of this work is
performed using these four variables.

During analysis, one outlier was found in the ensemble representing
C2W$\downarrow$, specifically in the temperature time series. This ensemble
member was the February 15, 1850, eruption, which was changed in favor of a
February 15, 1851, eruption in the C2W$\downarrow$, C2W$-$, and C2W$\uparrow$
ensembles. For completeness, the February 15, 1850, eruption is still included
in the online archive.

##  Appendix B External data

###  B.1 Otto-Bliesner data analysis

Data from ottobliesner2016 are the original input data of injected SO2 as used
in their model simulations, along with RF and temperature output data. The
injected SO2 can be found at https://www.cesm.ucar.edu/working-
groups/paleo/simulations/ccsm4-lm. Only the peak values of the SO2 dataset
were used in the analysis. Output variables are available at
www2.cesm.ucar.edu/models/experiments/LME.

Since the OB16 dataset contains a five-member ensemble, the final RF and
temperature time series used were ensemble means. A single control simulation
time series is used to remove seasonal dependence from the temperature, where
the control simulation is averaged into a climatology mean. Further, a drift
in the temperature is removed by subtracting a linear regression fit. RF has
seasonality removed in the Fourier domain.

The time of an eruption is found based on a best attempt at aligning the SO2
time series with both the RF time series and the temperature time series. The
RF and temperature peak values are taken as the value of the time series at
the time of an eruption according to the SO2 time series. Missing the true
peak means the found peaks will be biased towards lower values. However,
instances where eruptions occur close in time will contribute a bias to higher
values. These biases contribute to a greater uncertainty related to OB16 in
Figs. 4b,c,f.

###  B.2 Marshall data analysis

Data used to compute the M20 values were from marshall2020dataset, available
at https://doi.org/10.5285/232164e8b1444978a41f2acf8bbbfe91. As each file
includes a single eruption, peak values of AOD, RF, and temperature were found
by applying a Savitzky-Golay filter of third order and one-year window length,
and choosing the maximum value [Savitzky Golay (1964)].

###  B.3 Gregory data analysis

Data used to compute G16 values were kindly provided by Jonathan Gregory
(personal communication). The full 160-year-long time series were further
analysed by computing annual means.

##  Appendix C Model families

The model used here was the CESM2 with the WACCM6 atmosphere in the MA
configuration. The MA configuration uses the MAM3 [Gettelman . (2019)], a
simplified and computationally efficient default setting within the CAM5 [Liu
. (2016)], as described in liu2012. The MAM3 was developed from MAM7,
consisting of the seven modes Aitken, accumulation, primary carbon, fine dust,
fine sea salt, coarse dust, and coarse sea salt. Instantaneous internal mixing
of primary carbonaceous aerosols with secondary aerosols and instantaneous
ageing of primary carbonaceous particles are assumed by emitting primary
carbon in the accumulation mode [Liu . (2016)]. As dust absorbs water
efficiently and is expected to be removed by wet deposition similarly to sea
salt, fine dust is merged with fine sea salt into the accumulation mode and
coarse dust is merged with coarse sea salt into a coarse mode. The coarse mode
will quickly revert to its background state below the tropopause [Liu .
(2012)]. Consequently, MAM3 features the three modes Aitken, accumulation, and
coarse [Liu . (2016)].

The CESM2 is an ancestor of CESM1 used by OB16. They belong to a different
model family than both the HadCM3 (J05 and G16) and the UM-UKCA (M20), which
is an extended version of HadGEM3 [Dhomse . (2014)], and an ancestor of
HadCM3. A third model family is represented through ECHAM5 (N15) and MPI-ESM
(T10), where the latter is related to the former via the ECHAM6. A summary of
the model code genealogy is in table 4, based on the model code genealogy map
created by kuma2023.

Table 4: Model code family relationsa Family relation | Model name | Data name
---|---|---
CESM1 $\rightarrow$ CESM1-CAM5 $\rightarrow$ CESM2 | CESM1 | OB16
CESM2 | _This contribution_
| HadCM3 | J05, G16
HadCM3 $\rightarrow$ HadGEM1 $\rightarrow$ HadGEM2 $\rightarrow$ HadGEM3 $\rightarrow$ UM-UKCA | UM-UKCA | M20
ECHAM5 $\rightarrow$ ECHAM6 $\rightarrow$ MPI-ESM | ECHAM5 | N15
MPI-ESM | T10
aOverview of various model codes grouped into families according to the model
code genealogy map by kuma2023, with each table entry also indicating the
specific model code used in the referenced papers of this study.

AODVISstdn “stratospheric aerosol optical depth 550 nm day night” AOD
stratospheric aerosol optical depth CAM5 Community Atmosphere Model Version 5
CESM1 Community Earth System Model Version 1 CESM2 Community Earth System
Model Version 2 ECS equilibrium climate sensitivity ERF effective radiative
forcing FLNT “net longwave flux at the top of the model” FSNT “net solar flux
at the top of the model” IRF instantaneous radiative forcing MAM3 three mode
version of the Modal Aerosol Module MA middle atmosphere POP2 Parallel Ocean
Program Version 2 RF effective radiative forcing TCRP transient climate
response parameter TOA top-of-the-atmosphere TREFHT “reference height
temperature” WACCM6 Whole Atmosphere Community Climate Model Version 6 YTT
Young Toba Tuff

## Open Research Section

The direct output data of CESM2 are too large to be easily archived and
transferred. Instead, data generated directly from output fields of CESM2 are
made available in a NIRD Research Data Archive [Enger (20242)], and were
generated using scripts available at https://github.com/engeir/cesm-data-
aggregator. Analysis scripts are available at GitHub
(https://github.com/engeir/code-to-radiative-forcing-by-super-volcano-
eruptions) and is published to Zenodo [Enger (20241)]. Source code used to
generate CESM2 input files are available at
https://github.com/engeir/cesm2-volcano-setup.

###### Acknowledgements.

The simulations were performed on resources provided by Sigma2 — the National
Infrastructure for High Performance Computing and Data Storage in Norway. This
work was supported by the Tromsø Research Foundation under Grant Number
19_SG_AT. Thanks to both Maria Rugenstein and Martin Rypdal for valuable
discussions. We would also like to thank the authors of gregory2016,
ottobliesner2016, and marshall2020dataset for making their data available.

## References

  * Andersson . (2015) andersson2015Andersson, SM., Martinsson, BG., Vernier, JP., Friberg, J., Brenninkmeijer, CAM., Hermann, M.Zahn, A.  2015\.  Significant radiative impact of volcanic aerosol in the lowermost stratosphere Significant radiative impact of volcanic aerosol in the lowermost stratosphere.  Nature Communications67692-.  https://doi.org/10.1038/ncomms8692 10.1038/ncomms8692
  * Bender . (2010) bender2010Bender, FAM., Ekman, AML. Rodhe, H.  201010\.  Response to the eruption of Mount Pinatubo in relation to climate sensitivity in the CMIP3 models Response to the eruption of Mount Pinatubo in relation to climate sensitivity in the CMIP3 models.  Climate Dynamics355875–886.  http://link.springer.com/10.1007/s00382-010-0777-3 10.1007/s00382-010-0777-3
  * Boer . (2007) boer2007Boer, GJ., Stowasser, M. Hamilton, K.  200702\.  Inferring climate sensitivity from volcanic events Inferring climate sensitivity from volcanic events.  Climate Dynamics285481–502.  http://link.springer.com/10.1007/s00382-006-0193-x 10.1007/s00382-006-0193-x
  * Danabasoglu . (2020) danabasoglu2020Danabasoglu, G., Lamarque, JF., Bacmeister, J., Bailey, DA., DuVivier, AK., Edwards, J.Strand, WG.  2020\.  The Community Earth System Model Version 2 (CESM2) The community earth system model version 2 (CESM2).  Journal of Advances in Modeling Earth Systems122e2019MS001916.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019MS001916 e2019MS001916 2019MS001916  10.1029/2019MS001916
  * Dhomse . (2014) dhomse2014Dhomse, SS., Emmerson, KM., Mann, GW., Bellouin, N., Carslaw, KS., Chipperfield, MP.Thomason, LW.  2014\.  Aerosol microphysics simulations of the Mt. Pinatubo eruption with the UM-UKCA composition-climate model Aerosol microphysics simulations of the mt. pinatubo eruption with the um-ukca composition-climate model.  Atmospheric Chemistry and Physics142011221–11246.  https://acp.copernicus.org/articles/14/11221/2014/ 10.5194/acp-14-11221-2014
  * Douglass . (2006) douglass2006Douglass, DH., Knox, RS., Pearson, BD. Jr., AC.  2006\.  Thermocline flux exchange during the Pinatubo event Thermocline flux exchange during the pinatubo event.  Geophysical Research Letters3319L19711.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2006GL026355 10.1029/2006GL026355
  * Enger (20241) enger2024analysisEnger, ER.  2024102\.  Accompanying code to ’Radiative forcing by super- volcano eruptions’ Accompanying code to ’Radiative forcing by super- volcano eruptions’ [Software].  Zenodo.  https://doi.org/10.5281/zenodo.10722164 10.5281/zenodo.10722164
  * Enger (20242) enger2024Enger, ER.  20242\.  CESM2(WACCM6) single super-volcano simulations CESM2(WACCM6) single super-volcano simulations [Dataset].  Norstore.  https://doi.org/10.11582/2024.00025 10.11582/2024.00025
  * English . (2013) english2013English, JM., Toon, OB. Mills, MJ.  2013\.  Microphysical simulations of large volcanic eruptions: Pinatubo and Toba Microphysical simulations of large volcanic eruptions: Pinatubo and toba.  Journal of Geophysical Research: Atmospheres11841880–1895.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/jgrd.50196 https://doi.org/10.1002/jgrd.50196
  * Gettelman . (2019) gettleman2019Gettelman, A., Mills, MJ., Kinnison, DE., Garcia, RR., Smith, AK., Marsh, DR.Randel, WJ.  2019\.  The Whole Atmosphere Community Climate Model Version 6 (WACCM6) The whole atmosphere community climate model version 6 (WACCM6).  Journal of Geophysical Research: Atmospheres1242312380–12403.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JD030943 10.1029/2019JD030943
  * Giorgetta . (2006) giorgetta2006Giorgetta, MA., Manzini, E., Roeckner, E., Esch, M. Bengtsson, L.  2006\.  Climatology and Forcing of the Quasi-Biennial Oscillation in the MAECHAM5 Model Climatology and forcing of the quasi-biennial oscillation in the maecham5 model.  Journal of Climate19163882 - 3901.  https://journals.ametsoc.org/view/journals/clim/19/16/jcli3830.1.xml 10.1175/JCLI3830.1
  * Gregory . (2016) gregory2016Gregory, JM., Andrews, T., Good, P., Mauritsen, T. Forster, PM.  20161201\.  Small global-mean cooling due to volcanic radiative forcing Small global-mean cooling due to volcanic radiative forcing.  Climate Dynamics473979–3991.  https://doi.org/10.1007/s00382-016-3055-1 10.1007/s00382-016-3055-1
  * Hansen . (2005) hansen2005Hansen, J., Sato, M., Ruedy, R., Nazarenko, L., Lacis, A., Schmidt, GA.Zhang, S.  2005\.  Efficacy of climate forcings Efficacy of climate forcings.  Journal of Geophysical Research: Atmospheres110D18.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2005JD005776 10.1029/2005JD005776
  * Jones . (2005) jones2005Jones, GS., Gregory, JM., Stott, PA., Tett, SFB. Thorpe, RB.  20051201\.  An AOGCM simulation of the climate response to a volcanic super-eruption An AOGCM simulation of the climate response to a volcanic super-eruption .  Climate Dynamics257725–738.  https://doi.org/10.1007/s00382-005-0066-8 10.1007/s00382-005-0066-8
  * Kuma . (2023) kuma2023Kuma, P., Bender, FAM. Jönsson, AR.  2023\.  Climate Model Code Genealogy and Its Relation to Climate Feedbacks and Sensitivity Climate model code genealogy and its relation to climate feedbacks and sensitivity.  Journal of Advances in Modeling Earth Systems157e2022MS003588.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022MS003588 e2022MS003588 2022MS003588  https://doi.org/10.1029/2022MS003588
  * Liu . (2012) liu2012Liu, X., Easter, RC., Ghan, SJ., Zaveri, R., Rasch, P., Shi, X.Mitchell, D.  2012\.  Toward a minimal representation of aerosols in climate models: description and evaluation in the Community Atmosphere Model CAM5 Toward a minimal representation of aerosols in climate models: description and evaluation in the community atmosphere model CAM5.  Geoscientific Model Development53709–739.  https://gmd.copernicus.org/articles/5/709/2012/ 10.5194/gmd-5-709-2012
  * Liu . (2016) liu2016Liu, X., Ma, PL., Wang, H., Tilmes, S., Singh, B., Easter, RC.Rasch, PJ.  2016\.  Description and evaluation of a new four-mode version of the Modal Aerosol Module (MAM4) within version 5.3 of the Community Atmosphere Model Description and evaluation of a new four-mode version of the modal aerosol module (MAM4) within version 5.3 of the community atmosphere model.  Geoscientific Model Development92505–522.  https://gmd.copernicus.org/articles/9/505/2016/ 10.5194/gmd-9-505-2016
  * Marshall . (2019) marshall2019Marshall, LR., Johnson, JS., Mann, GW., Lee, L., Dhomse, SS., Regayre, L.Schmidt, A.  2019\.  Exploring How Eruption Source Parameters Affect Volcanic Radiative Forcing Using Statistical Emulation Exploring how eruption source parameters affect volcanic radiative forcing using statistical emulation.  Journal of Geophysical Research: Atmospheres1242964–985.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018JD028675 10.1029/2018JD028675
  * Marshall . (2022) marshall2022Marshall, LR., Maters, EC., Schmidt, A., Timmreck, C., Robock, A. Toohey, M.  20220504\.  Volcanic effects on climate: recent advances and future avenues Volcanic effects on climate: recent advances and future avenues.  Bulletin of Volcanology84554.  https://doi.org/10.1007/s00445-022-01559-3 10.1007/s00445-022-01559-3
  * Marshall . (2021) marshall2021Marshall, LR., Schmidt, A., Johnson, JS., Mann, GW., Lee, LA., Rigby, R. Carslaw, KS.  2021\.  Unknown Eruption Source Parameters Cause Large Uncertainty in Historical Volcanic Radiative Forcing Reconstructions Unknown eruption source parameters cause large uncertainty in historical volcanic radiative forcing reconstructions.  Journal of Geophysical Research: Atmospheres12613e2020JD033578.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020JD033578 e2020JD033578 2020JD033578  10.1029/2020JD033578
  * Marshall Smith (2020) marshall2020datasetMarshall, LR. Smith, CJ.  2020September25.  Vol-Clim: UM-UKCA interactive stratospheric aerosol model summary data for perturbed parameter ensemble of volcanic eruptions Vol-Clim: UM-UKCA interactive stratospheric aerosol model summary data for perturbed parameter ensemble of volcanic eruptions [dataset].  Centre for Environmental Data Analysis.  https://dx.doi.org/10.5285/232164e8b1444978a41f2acf8bbbfe91 10.5285/232164e8b1444978a41f2acf8bbbfe91
  * Marshall . (2020) marshall2020Marshall, LR., Smith, CJ., Forster, PM., Aubry, TJ., Andrews, T. Schmidt, A.  2020\.  Large Variations in Volcanic Aerosol Forcing Efficiency Due to Eruption Source Parameters and Rapid Adjustments Large variations in volcanic aerosol forcing efficiency due to eruption source parameters and rapid adjustments.  Geophysical Research Letters4719e2020GL090241.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020GL090241 e2020GL090241 2020GL090241  10.1029/2020GL090241
  * Marvel . (2016) marvel2016Marvel, K., Schmidt, GA., Miller, RL. Nazarenko, LS.  20160401\.  Implications for climate sensitivity from the response to individual forcings Implications for climate sensitivity from the response to individual forcings.  Nature Climate Change64386–389.  https://doi.org/10.1038/nclimate2888 10.1038/nclimate2888
  * Merlis . (2014) merlis2014Merlis, TM., Held, IM., Stenchikov, GL., Zeng, F. Horowitz, LW.  2014\.  Constraining Transient Climate Sensitivity Using Coupled Climate Model Simulations of Volcanic Eruptions Constraining transient climate sensitivity using coupled climate model simulations of volcanic eruptions.  Journal of Climate27207781–7795.  https://journals.ametsoc.org/view/journals/clim/27/20/jcli-d-14-00214.1.xml 10.1175/JCLI-D-14-00214.1
  * Mills . (2017) mills2017Mills, MJ., Richter, JH., Tilmes, S., Kravitz, B., MacMartin, DG., Glanville, AA.Kinnison, DE.  2017\.  Radiative and Chemical Response to Interactive Stratospheric Sulfate Aerosols in Fully Coupled CESM1(WACCM) Radiative and chemical response to interactive stratospheric sulfate aerosols in fully coupled CESM1(WACCM).  Journal of Geophysical Research: Atmospheres1222313,061–13,078.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017JD027006 10.1002/2017JD027006
  * Myhre . (2013) myhre2013Myhre, G., Shindell, D., Bréon, FM., Collins, W., Fuglestvedt, J., Huang, J.Zhang, H.  2013\.  Anthropogenic and natural radiative forcing Anthropogenic and natural radiative forcing.  TF. Stocker . (), Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change Climate change 2013: The physical science basis. contribution of working group i to the fifth assessment report of the intergovernmental panel on climate change ( 659–740).  Cambridge, UKCambridge University Press.  10.1017/CBO9781107415324.018
  * Niemeier Timmreck (2015) niemeier2015Niemeier, U. Timmreck, C.  201508\.  What is the limit of climate engineering by stratospheric injection of SO2 ? What is the limit of climate engineering by stratospheric injection of SO2 ?  Atmospheric Chemistry and Physics15169129–9141.  https://doi.org/10.5194%2Facp-15-9129-2015 10.5194/acp-15-9129-2015
  * Ollila (2016) ollila2016Ollila, A.  201602\.  Climate Sensitivity Parameter in the Test of the Mount Pinatubo Eruption Climate sensitivity parameter in the test of the mount pinatubo eruption.  Physical Science International Journal,91–14.  10.9734/PSIJ/2016/23242
  * Otto-Bliesner . (2016) ottobliesner2016Otto-Bliesner, BL., Brady, EC., Fasullo, J., Jahn, A., Landrum, L., Stevenson, S.Strand, G.  2016\.  Climate Variability and Change since 850 CE: An Ensemble Approach with the Community Earth System Model Climate variability and change since 850 CE: An ensemble approach with the community earth system model.  Bulletin of the American Meteorological Society975735–754.  https://journals.ametsoc.org/view/journals/bams/97/5/bams-d-14-00233.1.xml 10.1175/BAMS-D-14-00233.1
  * Pauling . (2023) pauling2023Pauling, AG., Bitz, CM. Armour, KC.  2023\.  The Climate Response to the Mt. Pinatubo Eruption Does Not Constrain Climate Sensitivity The climate response to the Mt. Pinatubo eruption does not constrain climate sensitivity.  Geophysical Research Letters507e2023GL102946.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2023GL102946 e2023GL102946 2023GL102946  10.1029/2023GL102946
  * Pitari . (2016) pitari2016bPitari, G., Genova, GD., Mancini, E., Visioni, D., Gandolfi, I. Cionni, I.  2016\.  Stratospheric Aerosols from Major Volcanic Eruptions: A Composition-Climate Model Study of the Aerosol Cloud Dispersal and e-folding Time Stratospheric aerosols from major volcanic eruptions: A composition-climate model study of the aerosol cloud dispersal and e-folding time.  Atmosphere76.  https://www.mdpi.com/2073-4433/7/6/75 10.3390/atmos7060075
  * Richardson . (2019) richardson2019Richardson, TB., Forster, PM., Smith, CJ., Maycock, AC., Wood, T., Andrews, T.Watson-Parris, D.  2019\.  Efficacy of Climate Forcings in PDRMIP Models Efficacy of climate forcings in PDRMIP models.  Journal of Geophysical Research: Atmospheres1242312824–12844.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JD030581 10.1029/2019JD030581
  * Robock (2000) robock2000Robock, A.  2000\.  Volcanic eruptions and climate Volcanic eruptions and climate.  Reviews of Geophysics382191–219.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/1998RG000054 10.1029/1998RG000054
  * Salvi . (2022) salvi2022Salvi, P., Ceppi, P. Gregory, JM.  2022\.  Interpreting Differences in Radiative Feedbacks From Aerosols Versus Greenhouse Gases Interpreting differences in radiative feedbacks from aerosols versus greenhouse gases.  Geophysical Research Letters498e2022GL097766.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022GL097766 e2022GL097766 2022GL097766  10.1029/2022GL097766
  * Savitzky Golay (1964) savitzky1964Savitzky, A. Golay, MJE.  1964\.  Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Smoothing and differentiation of data by simplified least squares procedures.  Analytical Chemistry3681627-1639.  https://doi.org/10.1021/ac60214a047 10.1021/ac60214a047
  * Sigl . (2022) sigl2022Sigl, M., Toohey, M., McConnell, JR., Cole-Dai, J. Severi, M.  2022\.  Volcanic stratospheric sulfur injections and aerosol optical depth during the Holocene (past 11 500 years) from a bipolar ice-core array Volcanic stratospheric sulfur injections and aerosol optical depth during the holocene (past 11 500 years) from a bipolar ice-core array.  Earth System Science Data1473167–3196.  https://essd.copernicus.org/articles/14/3167/2022/ 10.5194/essd-14-3167-2022
  * CJ. Smith . (2018) smith2018Smith, CJ., Kramer, RJ., Myhre, G., Forster, PM., Soden, BJ., Andrews, T.Watson-Parris, D.  2018\.  Understanding Rapid Adjustments to Diverse Forcing Agents Understanding rapid adjustments to diverse forcing agents.  Geophysical Research Letters452112,023-12,031.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018GL079826 https://doi.org/10.1029/2018GL079826
  * R. Smith . (2010) smith2010Smith, R., Jones, P., Briegleb, B., Bryan, F., Danabasoglu, G., Dennis, J.Yeager, S.  20100323\.  The Parallel Ocean Program (POP) reference manual The parallel ocean program (POP) reference manual.  LAUR-10-01853.  https://www.cesm.ucar.edu/models/cesm1.0/pop2/doc/sci/POPRefManual.pdf
  * Solomon . (2011) solomon2011Solomon, S., Daniel, JS., Neely, RR., Vernier, JP., Dutton, EG. Thomason, LW.  2011\.  The Persistently Variable “Background” Stratospheric Aerosol Layer and Global Climate Change The persistently variable “background” stratospheric aerosol layer and global climate change.  Science3336044866-870.  https://www.science.org/doi/abs/10.1126/science.1206027 10.1126/science.1206027
  * Stier . (2005) stier2005Stier, P., Feichter, J., Kinne, S., Kloster, S., Vignati, E., Wilson, J.Petzold, A.  2005\.  The aerosol-climate model ECHAM5-HAM The aerosol-climate model echam5-ham.  Atmospheric Chemistry and Physics541125–1156.  https://acp.copernicus.org/articles/5/1125/2005/ 10.5194/acp-5-1125-2005
  * Timmreck . (2010) timmreck2010Timmreck, C., Graf, HF., Lorenz, SJ., Niemeier, U., Zanchettin, D., Matei, D.Crowley, TJ.  2010\.  Aerosol size confines climate response to volcanic super-eruptions Aerosol size confines climate response to volcanic super-eruptions.  Geophysical Research Letters3724.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2010GL045464 10.1029/2010GL045464
  * Timmreck . (2018) timmreck2018Timmreck, C., Mann, GW., Aquila, V., Hommel, R., Lee, LA., Schmidt, A.Weisenstein, D.  2018\.  The Interactive Stratospheric Aerosol Model Intercomparison Project (ISA- MIP): motivation and experimental design The interactive stratospheric aerosol model intercomparison project (ISA- MIP): motivation and experimental design.  Geoscientific Model Development1172581–2608.  https://gmd.copernicus.org/articles/11/2581/2018/ 10.5194/gmd-11-2581-2018
  * Toohey . (2011) toohey2011Toohey, M., Krüger, K., Niemeier, U. Timmreck, C.  2011\.  The influence of eruption season on the global aerosol evolution and radiative impact of tropical volcanic eruptions The influence of eruption season on the global aerosol evolution and radiative impact of tropical volcanic eruptions.  Atmospheric Chemistry and Physics112312351–12367.  https://acp.copernicus.org/articles/11/12351/2011/ 10.5194/acp-11-12351-2011
  * Toohey . (2019) toohey2019Toohey, M., Krüger, K., Schmidt, H., Timmreck, C., Sigl, M., Stoffel, M. Wilson, R.  20190201\.  Disproportionately strong climate forcing from extratropical explosive volcanic eruptions Disproportionately strong climate forcing from extratropical explosive volcanic eruptions.  Nature Geoscience122100–107.  https://doi.org/10.1038/s41561-018-0286-2 10.1038/s41561-018-0286-2
  * Vidal . (2016) vidal2016Vidal, CM., Métrich, N., Komorowski, JC., Pratomo, I., Michel, A., Kartadinata, N.Lavigne, F.  201610\.  The 1257 Samalas eruption (Lombok, Indonesia): the single greatest stratospheric gas release of the Common Era The 1257 samalas eruption (lombok, indonesia): the single greatest stratospheric gas release of the common era.  Scientific Reports61.  http://dx.doi.org/10.1038/srep34868 10.1038/srep34868
  * Wigley . (2005) wigley2005Wigley, TML., Ammann, CM., Santer, BD. Raper, SCB.  2005\.  Effect of climate sensitivity on the response to volcanic forcing Effect of climate sensitivity on the response to volcanic forcing.  Journal of Geophysical Research: Atmospheres110D9.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2004JD005557 10.1029/2004JD005557
  * Zanchettin . (2016) zanchettin2016Zanchettin, D., Khodri, M., Timmreck, C., Toohey, M., Schmidt, A., Gerber, EP.Tummon, F.  2016\.  The Model Intercomparison Project on the climatic response to Volcanic forcing (VolMIP): experimental design and forcing input data for CMIP6 The model intercomparison project on the climatic response to volcanic forcing (volMIP): experimental design and forcing input data for CMIP6.  Geoscientific Model Development982701–2719.  https://gmd.copernicus.org/articles/9/2701/2016/ 10.5194/gmd-9-2701-2016
  * Zanchettin . (2019) zanchettin2019Zanchettin, D., Timmreck, C., Toohey, M., Jungclaus, JH., Bittner, M., Lorenz, SJ. Rubino, A.  2019\.  Clarifying the Relative Role of Forcing Uncertainties and Initial-Condition Unknowns in Spreading the Climate Response to Volcanic Eruptions Clarifying the relative role of forcing uncertainties and initial-condition unknowns in spreading the climate response to volcanic eruptions.  Geophysical Research Letters4631602–1611.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018GL081018 10.1029/2018GL081018