A New Multi‐Method Assessment of Stratospheric Sulfur Load From the Okmok II Caldera‐Forming Eruption of 43 BCE

The 43 BCE eruption of Okmok Volcano has been proposed to have had a significant climate cooling impact in the Northern Hemisphere. In this study, we quantify the climate cooling potential of the Okmok II eruption by measuring sulfur concentration in melt inclusions (up to 1,606 ppm) and matrix glasses and estimate a total of 62 ± 16 Tg S released. The proportion reaching the stratosphere (2.5%–25%, i.e., 1.5–15.5 Tg S) was constrained by physical modeling of the caldera‐collapse eruption. Using the NASA Goddard Institute for Space Studies E2.2 climate model we found a linear response between cooling and stratospheric sulfur load (0.05–0.08°C/Tg S). Thus, the 1–2°C of cooling derived from proxy records would require 16–32 Tg sulfur injection. This study underscores the importance of combining approaches to estimate stratospheric S load. For Okmok II, we find all methods are consistent with a range of 15–16 Tg S.

The primary goal of this study is to provide the first petrological measurements of the total amount of sulfur released by Okmok II magma.The other goal is to estimate the fraction of this load injected to the stratosphere, critical to determining its climate impact.This quantity is challenging to determine by any one method, and the novelty of our study is that we employ three independent methods and data sets for estimating stratospheric sulfur: (a) ice sheet deposition and sulfur isotopes, (b) global climate model and climate proxy data, and (c) total petrologic sulfur load and eruption dynamics modeling for stratospheric injection.

Materials and Methods
Okmok II PDC and andesite fall samples were selected for melt inclusion study (Figure 1a, Table S3 in Supporting Information S1).Polished crystals and matrix glasses from four samples were measured by electron microprobe for major element composition and sulfur concentration.Selected melt inclusions and matrix glasses were then measured for water and sulfur content by ion microprobe (Figure 1d; Texts S3 and S4 in Supporting Information S1).
To quantify the climate impacts of different stratospheric sulfur injection amplitudes, we simulated the effect of the Okmok II eruption on the global climate system using the NASA Goddard Institute for Space Studies (GISS) Model E2.2 (Orbe et al., 2020;Rind et al., 2020).It was configured following DallaSanta and Polvani (2022) and is a state of the art atmosphere-ocean-land-sea ice coupled model which participated in the Coupled Model Intercomparison Project, Phase 6 (CMIP6) (Text S7 in Supporting Information S1).

Melt Inclusion Composition
The sulfur concentrations in 19 melt inclusions from Okmok II range from 152 to 1,606 ppm S (with 10% relative uncertainty in the ionprobe analysis) in melt compositions that span from basaltic to rhyodacitic (51-65 wt% SiO 2 , concentrations expressed volatile-free).To estimate sulfur load, it is critical to determine the melt inclusion compositions that best reflects the ∼50 km 3 of bulk PDC magma, as both the whole rocks and matrix glasses are degassed in sulfur.The more evolved composition of the GISP2 tephra (54.7-59.9wt% SiO 2 ) and PDC matrix glass (57.1-60.8wt% SiO 2 ) is inferred to be the result of degassing-driven crystallization of microlites upon ascent as described by Larsen et al. (2007), and thus does not represent the bulk PDC composition (Figure 1b; Text S5 in Supporting Information S1).The MI compositions in AOK147-plag have SiO 2 (55.2-56.4wt%) that overlap best with the PDC whole rocks (54.6-55.7 wt% SiO 2 ), which we take to best represent the bulk composition of the PDC magma (Figure 1b).These inclusions have the highest sulfur concentrations and are in equilibrium with their plagioclase hosts and the bulk PDC composition, suggesting they represent the original volatile content of the magma before degassing, with sulfide saturation inferred from the Cu-Th trend in Okmok whole rocks (1,263-1,606 ppm S, Figure 1c, Text S6 in Supporting Information S1).Selection of these inclusions is critical in accounting for the sulfur concentration in the parental PDC magma body at depth, which would otherwise be underestimated by ∼1,000 ppm S if the petrologic method was applied to melt inclusions similar in composition to the matrix glass (Su et al., 2016;Figure 1c).We thus take the pre-eruptive concentration for the Okmok II PDC magma as the maximum sulfur measured in these inclusions: 1,606 ± 160 ppm.Total sulfur concentration in the parental PDC magma could be even higher due to minor sulfide fractionation (see Cu-Th systematics in Figure S4 in Supporting Information S1) or prior degassing.

Total Mass Erupted and Total Sulfur Load
Calculating the total sulfur load from Okmok II requires an estimate of the erupted mass.PDCs account for 98% of the volume erupted in Okmok II, and thus represent the most important volume component for estimating total sulfur load.Burgisser (2005) estimated a total of 24 km 3 of PDC deposits on Umnak and Unalaska Island, but recognized a large uncertainty in the volume of underwater deposits.To address this, we assume two areas of deposition around Umnak Island that extend the thickness of the deposits measured on land to sea (see Figure S3 in Supporting Information S1).This produces an additional 25 km 3 of pyroclastic flows deposited underwater for a total volume during Okmok II PDC stage of 50 ± 10 km 3 , consistent with the volume of the caldera (Burgisser, 2005).We calculate a total erupted magma mass of 4.29 ± 1.01 × 10 13 kg, based on average deposit density and corrections for pre-existing lithic material and deposit void fraction (Burgisser, 2005; Text S1 in Supporting Information S1).The sulfur load is then derived from Equation 1 of Devine et al. (1984): 10.1029/2023GL103334 3 of 10 where M s is the mass of sulfur erupted, W x is the mass fraction of crystals, C i is the concentration of sulfur as measured in melt inclusions, and C m is the concentration of sulfur as measured in degassed matrix glass.W x is taken as 1.7% ± 1.8% (Burgisser, 2005).C i is 1,606 ± 160 ppm (as above) and C m is 127 ± 23 ppm S. Thus, ∼1,500 ppm of sulfur was degassed during this eruption, producing a total sulfur load of 62 ± 16 Tg S (Equation 1, Table S1 in Supporting Information S1).This is likely a minimum estimate for total Okmok II sulfur emission, given our assumptions of zero pre-exsolved gas phase.Further, seawater vapourization as a result of contact with hot PDC material could have led to additional S release (Rowell et al., 2022), though it's unlikely that this source would reach a significant altitude (Dufek et al., 2007).

Climate Model Scenarios
We performed model simulations to study the effect of three different stratospheric S loads: 62 Tg (total petrologic S load), 38 Tg (estimated by McConnell et al., 2020 from ice core load), and 10 Tg (lower bound) and averaged response over 43 and 42 BCE, when the signal is largest (Figure 2d).In the Mediterranean region, the  2020), who used a different climate model than ours, but an identical aerosol forcing model.Over the Northern Hemisphere (Figures S1a-S1d in Supporting Information S1), we observe highly significant surface temperature responses.The milder responses over Northern Eurasia, as compared to North America, are associated with a positive phase of the Northern Annular Mode in boreal winter, which emerges for sufficiently large eruptions (DallaSanta & Polvani, 2022).The lack of response in the Southern Hemisphere is due to the Easy Volcanic Aerosol (EVA) parameterization of hemispheric asymmetry, though a small volcanic sulfate peak is observed in Antarctic ice cores around the same time (Figure S1 in Supporting Information S1; Toohey et al., 2016).Notably, the climate response, as defined by the ensemble averaged departure from control run, shows roughly linear cooling as a function of S load (0.05-0.08°C cooling/Tg S) (Figure 4).
Overall, we find reduced precipitation over Northern Europe and increased precipitation over Southern Europe following the eruption, with magnitudes up to 12 mm/month.Our dipole response for 38 Tg is slightly shifted in latitude from that in McConnell et al. (2020), pointing to intermodel uncertainties.As with temperature, the precipitation response strengthens approximately linearly with stratospheric sulfur injection mass, with significance emerging at 38 Tg S (Figure 2f) and strengthening at 62 Tg S (Figure 2g).However, the month-to-month spread across ensemble members is high (Figure 2h).Thus, the Roman Empire likely experienced overall precipitation changes after the eruption, but any given month could have been wetter or drier.The precipitation changes over Africa reflect heightened interannual variability due to contraction of the Hadley Cell, which could have led to drying at the Nile headwaters (Figures 2e-2i, Figures S1e-S1i in Supporting Information S1), consistent with historical records (McConnell et al., 2020).On a global scale (Figures S1e-S1i in Supporting Information S1), the precipitation response is more robust, associated with changes to the zonally symmetric circulation.
Lastly, we find that the month of eruption does not appear to be critical for the temperature and precipitation impacts (Figure S2 in Supporting Information S1).This provides confidence that our results are meaningful for historical interpretation, as the actual eruption month of the Okmok II eruption is not known.

Among the Largest Sulfur Loads in History
A petrologic sulfur load of 62 Tg S for Okmok II is on par with historical eruptions of similar size (Laki, Samalas, Paektu, Minoan), and ranks as one the largest volcanic sulfur loads in the last 2,500 years (Figure 3).Comparable eruptions have been linked to climate instability and major historical milestones: the 1783 flood-lava eruption  S2 in Supporting Information S1) and (b) ice core-based estimates of stratospheric sulfur loading (b).Lines in (a) show sulfur emissions for magmas that undergo 100, 500, 1,000, 2,000, and 3,000 ppm of sulfur degassing, for reference.Dashed lines in panel (b) show the offset between total sulfur emission from petrological studies and the stratospheric S estimated from ice core studies.The difference between the two likely reflects variable stratospheric injection of the total sulfur load during any given eruption.Eruptions are colored by dominant composition, with eruption years noted.Dates in panel (b) extend from 500 BCE to 1990 CE, reflecting the coverage of the ice core record presented in the eVolv2k v3 database (Toohey & Sigl, 2017).Eruption name codes: (AG: Agung, EL:E ldgja, HP: Huayanaputina, LK: Laki, KR: Krakatau, KW: Kuwae, MN: Minoan/ Thera, PK: Paektu, PT: Pinatubo, SA: Samalas, SM: Santa Maria, TBJ: Ilopango TBJ, and TM: Tambora).
of Laki emitted around 60 Tg of sulfur and drove significant cooling in the Northern Hemisphere (Schmidt et al., 2011;Thordarson & Self, 2003;Zambri et al., 2019), while the 1257 eruption of Mt.Samalas generated ∼1-2°C global cooling and has been connected with the initiation of the Black Death in Europe (Büntgen et al., 2022;Fell et al., 2020;Guillet et al., 2017).Notably, petrologic evidence suggests very little halogen degassing from Okmok II (∼0 Tg for Cl and F), with matrix glass concentrations of Cl and F near or greater than that observed in melt inclusions.
Only S injected into the stratosphere exerts a climate cooling influence (e.g., Robock, 2000).Sulfate concentrations in ice cores have been used to estimate stratospheric sulfur loading, relying on the transport and deposition of volcanic sulfate aerosols to ice sheets 1-2 years post-eruption (Toohey & Sigl, 2017).While this strategy is useful for identifying the timing of large eruptions and their relative magnitudes, it does not provide source constraints nor does it distinguish stratospheric fall-out versus sulfur transported through the troposphere (Marshall et al., 2021).In Figure 3b, we compare petrological estimates of sulfur load and ice core estimates over the past 2,500 years, and find that petrologic estimates are systematically higher than the modeled ice core deposition.This offset likely reflects variable stratospheric injection due to the physical dynamics of any given eruption and illustrates the importance of studies of past volcanic sulfur loads.
To address partial stratospheric injection of the total S load from the Okmok II eruption, we compare below three independent methods: (a) ice core deposition, (b) eruption dynamics, and (c) climate proxy data.

Stratospheric Injection From Ice Cores
Synchronized estimates made by Toohey and Sigl (2017) predict stratospheric sulfur injection of 38.6 ± 11.3 Tg S from Okmok II based on its source location and ice core sulfate concentrations.McConnell et al. (2020) assumed all of this was stratospheric, and used this as input to their climate model.A recent study (Pearson et al., 2022) Figure 4. Land-averaged cooling response following modeled stratospheric sulfur injections (10, 38, 62 Tg S) from Okmok II over three representative areas: Northern Europe (green line, 40°-70°N, 10°W, 30°E), Asia (orange line, 30°-55°N, 40°, 140°E), and the Northern Hemisphere (blue line).The three injection scenarios show a roughly linear cooling response for all regions.Proxy cooling reconstruction is denoted by the shaded areas: tree ring temperature reconstructions from Northern Europe predict −1.52 ± 0.42°C of cooling in 43 BCE (Luterbacher et al., 2016), while a Chinese speleothem record predicts −1.76 ± 0.2°C of cooling within a 5-year dating uncertainty (Tan et al., 2003).Combined with climate model runs, overlapping proxy evidence suggests a stratospheric sulfur load of 18-22 Tg S. Stratospheric sulfur injection and observed change in global mean surface temperature from the 1991 eruption of Pinatubo is noted for reference (Guo et al., 2004;Pauling et al., 2023).
proposes an even larger stratospheric sulfur load of 48 ± 15 Tg S, which would represent ∼77% of the total erupted sulfur given our petrologic estimate.This same study, however, makes the important point that some ice-deposited S will be tropospheric, especially for volcanoes like Okmok which are relatively proximal to Greenland.
A closer examination of mass independent fractionation in sulfur isotopes (S-MIF) may be used to resolve the stratospheric versus tropospheric ice deposition of sulfur (Burke et al., 2019).A study of the Aniakchak II eruption (1628 BCE) uses mass balance, along with triple isotope measurements of sulfate deposited in Greenland to estimate the proportion of ice-deposited sulfate that reached the stratosphere (Pearson et al., 2022).The Aniakchak II eruption, which had a similar ice deposition (52 Tg) to Okmok II (48 Tg) and is at similar latitude in Alaska (56.9°N) to Okmok (53.4°N), deposited a sulfur peak on Greenland that was calculated to be 61% stratospheric and 39% tropospheric.Given that the magnitude of the S-MIF peak measured for Okmok was much smaller than that for Aniakchak (∆ 33 S max = 0.20‰ for Okmok II vs. ∆ 33 S max = 1.30‰ for Aniakchak II), the stratospheric proportion for the total volcanic S deposited following the Okmok II eruption will be <61%.As the measurements in McConnell et al. (2020) did not include key sulfur concentrations or background isotope ratios, it is not currently possible to carry out a meaningful sulfur mass balance for Okmok.Furthermore, S-MIF may not occur in the lower stratosphere at high latitudes, and thus relying solely on ∆ 33 S may result in an underestimation of stratospheric sulfur injection and resulting climate cooling (Thordarson & Self, 2003).Nonetheless, the existing measurements suggest a smaller stratospheric proportion for Okmok II versus Aniakchak II, and thus we can conclude that <29 Tg (61% of 48 Tg) of Okmok II sulfur load measured in ice cores was stratospheric.This would constitute <47% of the total S erupted.

Stratospheric Injection From Eruption Dynamics Models
The presence of basaltic andesite tephra and S-MIF in Greenland ice cores indicate that the Okmok II PDC material reached across the Northern Hemisphere and had a stratospheric component (McConnell et al., 2020).This presents a conundrum, as the Okmok II PDC deposits are not associated with an airfall deposit.In order to address how a predominantly pyroclastic flow-generating eruption could also produce stratospheric sulfur and hemispheric ash without producing a Plinian air-fall, a companion study uses a numerical model to explore the conditions that could lead to this eruptive behavior.Burgisser et al. (2023) apply the MFIX-TFM two-phase flow model to simulate the dynamics of PDC emplacement.Model results are constrained by field stratigraphic observations of asymmetric deposition, maximum flow runout, a lack of airfall deposits and total mass erupted, with the stratospheric injection of gas during the eruption estimated via mass balance based on the fraction of stratospheric solids.
The results of this physical model indicate that 2.5%-25% of the solids and gases from the PDC phase of the Okmok II eruption would be injected into the stratosphere (Burgisser et al., 2023).The low proportion of stratospheric injection reflects the lack of an airfall deposit from the central plume, with the majority of stratospheric injection occurring via successive pulses of co-ignimbrite material issued from the main pyroclastic flows.We note that complex eruption dynamics may also affect the transport of S into the stratosphere, such as scavenging in a water-rich plume (Thordarson et al., 1996).Nevertheless, if only 2.5%-25% of the total Okmok II sulfur is stratospheric, that would result in a stratospheric injection of 1.55-15.5Tg, well below the 29 Tg maximum outlined above based on ice deposition.

Constraints From Climate Models and Proxies
Temperature proxy evidence from tree ring and speleothem records (McConnell et al., 2020;Sigl et al., 2015) indicate >2°C cooling in the Northern Hemisphere around the time of the Okmok II eruption.As discussed above, our climate modeling indicates a roughly linear relationship between cooling and mass of S injected into the stratosphere (Figure 4).Thus, the amount of cooling recorded by temperature proxies can constrain the amount of stratospheric S needed to create such a signal.

Multi-Method Convergence of Stratospheric S Injection
The upper limit from the eruption model (15.5 Tg) approaches the lower limit of the climate model (16 Tg, as discussed above) and both are permissible given the upper bound from ice deposition (29 Tg).Thus, 15-16 Tg of stratospheric S from Okmok II is consistent with all data, models and approaches.
Despite this coherence, it is worth considering the uncertainties due to the sparsity of data and number of assumptions.It is possible that additional gaseous sulfur could have been injected during the rhyodacitic Plinian phase that preceded the PDCs and caldera collapse.However, the lack of any glass >60% SiO 2 in the Greenland ice cores is inconsistent with the earlier Plinian phases producing a major stratospheric load.Alternatively, it is possible that our petrologic estimate is too low.In order to overlap with the proxy cooling-based estimate for Okmok II (18-22 Tg S), and assuming a 25% stratospheric portion, the total petrologic load would have to be increased to 80 Tg S. A higher petrologic load could be achieved with additional MI measurements that could extend to S concentrations higher than 1,600 ppm, although we note that Zimmer et al. (2010) measured less than 1,300 ppm in all 20 MI from several more recent mafic eruptions of Okmok (Figure 1e).A less conservative estimate of total eruptive volume, or the addition of exsolved sulfur from an unerupted primitive parent could also supply an additional sulfur load.
Conversely, the observed climate cooling could be overestimated.Proxy estimates for cooling on interannual timescales are relatively sparse due to the lack of ancient tree ring records.It is also questionable how well a regionally averaged ensemble of climate models might relate to location-specific (and seasonally specific) data.The GISS climate model results are also very likely to depend on the injection height and aerosol size, and the EVA model used here may not be a good representation of the Okmok II eruption.
Regardless, all evidence suggests less than half of the total Okmok II petrologic sulfur load was injected into the stratosphere, where it produced a significant climate cooling effect.We see the relative agreement in multiple reconstructions as a success that demonstrates the utility of this three-pronged approach to evaluate other eruptions with more robust climate proxy data.

Conclusions
The Okmok II eruption released a total sulfur load of 62 ± 16 Tg S, one of the largest estimated sulfur loads in recorded history.However, eruption dynamics models and ice core S-MIF data suggest that much of the sulfur released during Okmok II did not reach the stratosphere.Thus, only a fraction of this total load was available to impact surface climate.Comparison of model-generated climate anomalies with coincident cooling recorded by tree ring and speleothem proxies suggests that a stratospheric sulfur load of 18-22 Tg is consistent with the observed climate effect.This constitutes ∼30% of the total petrologic sulfur load, in good agreement with the upper end of stratospheric injection predicted by physical modeling (25%) and permissible given the ice core S-MIF (<61%).
Moreover, we find a systematic mismatch between petrological and ice core estimates of volcanic sulfur loads, suggesting that partial stratospheric injection should be considered when examining other large eruptions.Okmok II serves as a case study for an interdisciplinary approach to reconstructing major historical eruptions and their related impacts-one that combines petrological source constraints, physical modeling, and climate models to interpret the ice core record.

Figure 1 .
Figure 1.Stratigraphic sequence of Okmok II units (a, left) and map of sample locations across Umnak and Unalaska Islands ((a), right).Units are labeled by type and composition (RD = rhyodacite fall, AN = andesite fall, and PDC = pyroclastic density current), with bulk rock SiO 2 wt% in brackets.Fall deposits are noted by units (A, B, C) identified in Burgisser (2005), and are scaled relative to maximum deposit thickness.Pyroclastic flow facies have an average thickness of 30-60 m depending on distance from the caldera.Panels (b, c) show sulfur and major element concentrations in melt inclusions, matrix glasses, and whole-rocks from Okmok II.CaO and SiO 2 are reported volatile-free, normalized to 100%.Matrix glasses are plotted together in panel (b), and separated in panel (c).They largely overlap in our data set.Photomicrographs (d) are of representative olivine-hosted melt inclusions from the JLOK6H andesite fall (left) and plagioclase-hosted melt inclusions from pyroclastic density current deposit sample AOK147 (which contains the maximum S concentrations).The H 2 O and S concentrations in panel (e) were measured by ionprobe, and include melt inclusion measurements from the 1997 eruption of Okmok from Zimmer et al. (2010).Panel (a) made with GeoMapApp (Ryanet al., 2009).

Figure 2 .
Figure 2. Impact of Okmok II eruption on Mediterranean surface temperature (column 1) and precipitation (column 2).The responses are defined as ensemble-averaged departures from the corresponding control run, averaged over 43 and 42 BCE.Stippling indicates an inability to reject the null hypothesis at 95% confidence.The three sulfate amplitudes are labeled for rows 1-3.Row 4 shows the time evolution of the responses, averaged over the specified region.Dots indicate where all 10 members are less than zero.

Figure 3 .
Figure 3. Okmok II sulfur load compared with (a) sulfur emission estimates of other historical eruptions (TableS2in Supporting Information S1) and (b) ice core-based estimates of stratospheric sulfur loading (b).Lines in (a) show sulfur emissions for magmas that undergo 100, 500, 1,000, 2,000, and 3,000 ppm of sulfur degassing, for reference.Dashed lines in panel (b) show the offset between total sulfur emission from petrological studies and the stratospheric S estimated from ice core studies.The difference between the two likely reflects variable stratospheric injection of the total sulfur load during any given eruption.Eruptions are colored by dominant composition, with eruption years noted.Dates in panel (b) extend from 500 BCE to 1990 CE, reflecting the coverage of the ice core record presented in the eVolv2k v3 database(Toohey & Sigl, 2017).Eruption name codes: (AG: Agung, EL:E ldgja, HP: Huayanaputina, LK: Laki, KR: Krakatau, KW: Kuwae, MN: Minoan/ Thera, PK: Paektu, PT: Pinatubo, SA: Samalas, SM: Santa Maria, TBJ: Ilopango TBJ, and TM: Tambora).