PalVol v1: a proxy-based semi-stochastic ensemble reconstruction of volcanic stratospheric sulfur injection for the last glacial cycle (140 000–50 BP)

Perturbations in stratospheric aerosol due to explosive volcanic eruptions are a primary contributor to natural climate variability. Observations of stratospheric aerosol are available for the past decades, and information from ice cores has been used to derive estimates of stratospheric sulfur injections and aerosol optical depth over the Holocene (approximately 10 000 BP to present) and into the last glacial period, extending back to 60 000 BP. Tephra records of past volcanism, compared to ice cores, are less complete but extend much further into the past. To support model studies of the potential impacts of explosive volcanism on climate variability across timescales, we present here an ensemble reconstruction of volcanic stratospheric sulfur injection (VSSI) over the last 140 000 years that is based primarily on terrestrial and marine tephra records. VSSI values are computed as a simple function of eruption magnitude based on VSSI estimates from ice cores and satellite observations for identified eruptions. To correct for the incompleteness of the tephra record, we include stochastically generated synthetic eruptions assuming a constant background eruption frequency from the ice core Holocene record. While the reconstruction often differs from ice core estimates for specific eruptions due to uncertainties in the data used and reconstruction method, it shows good agreement with an ice-core-based VSSI reconstruction in terms of millennial-scale cumulative VSSI variations over the Holocene. The PalVol reconstruction provides a new basis to test the contributions of forced vs. unforced natural variability to the spectrum of climate and the mechanisms leading to abrupt transitions in the palaeoclimate record with low- to high-complexity climate models. The PalVol volcanic forcing reconstruction is available at https://doi.org/10.26050/WDCC/PalVolv1 (Toohey and Schindlbeck-Belo, 2023).

1 Introduction

Explosive volcanic eruptions transfer massive quantities of material from the solid Earth into the atmosphere. Eruptive plumes contain large amounts of solid material, as well as gaseous compounds including water vapor, carbon dioxide, and sulfur-containing species (mostly SO₂), often being combined into aerosols. The solid volcanic fragments are fragmented magma that is ejected by an eruption and are called tephra (mostly ash at ∅<2 mm and lapilli at ∅2 to 64 mm). They generally fall back to the surface of the Earth rather quickly, producing a tephra layer with decreasing thickness and grain size with increasing distance to the volcano. Such deposits persist as a record of past volcanic eruptions that can be seen in outcrops or in sediment cores extracted from marine and lacustrine environments (e.g., Kutterolf et al., 2016; Schindlbeck et al., 2016).

Gaseous emissions from eruptions can persist in the atmosphere much longer than solid emissions. While the emission of H₂O and CO₂ from a single eruption is generally insignificant compared to the atmosphere burden of these species, volcanic emissions of sulfur-containing species can produce significant increases in atmospheric sulfur (Robock, 2000). Under atmospheric conditions, sulfur-containing species produce sulfate aerosol particles, small droplets of sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$), and water in liquid form that are dispersed in a gaseous matrix (Robock, 2000). Sulfur emissions to the troposphere (from volcanic eruptions or anthropogenic activities) produce sulfate aerosols that have an atmospheric lifetime of days to weeks as they grow in size and are eventually rained out (Hamill et al., 1997). Sulfate aerosols in the cold and dry lower stratosphere do not generally grow as large as those in the troposphere and as a result persist in the stratosphere for months to years, over which time they are transported around the globe. These sulfate aerosols have important impacts on atmospheric radiative transfer by scattering solar radiation and absorbing longwave radiation, with the net effect of a decrease in downwelling net radiation at the Earth's surface, which leads to a cooling of surface temperatures (Robock, 2000).

Post-volcanic large-scale (global or hemispheric) cooling has been observed after recent eruptions (e.g., El Chichón in 1982, Mt. Pinatubo in 1991) and is apparent in millennial-scale climate reconstructions, for example the 1815 Tambora eruption in Indonesia (e.g., Rampino and Self, 1982, 1993). It has been shown that negative radiative forcing from volcanic aerosol perturbations is the primary driver of natural climate variability over the past 1000 to 2000 years (Sigl et al., 2015; Schurer et al., 2013). Additionally, Kobashi et al. (2017) showed that, during the Holocene, volcanic eruptions also played an important role in centennial to millennial temperature variability in Greenland. Representing the intermittent natural forcing in climate model experiments for the late Holocene and the Glacial period improves modeled variability on timescales from decades to centuries (Ellerhoff et al., 2022). Abbott et al. (2021) provide a 400-year reconstruction of volcanic forcing from 13.2 to 12.8 ka, which supports model simulations for the Younger Dryas inception.

The role of volcanic eruptions in longer-term climate variability has long been speculated about but remains poorly understood. Strong volcanic eruptions have been linked to multi-decadal periods of cooler-than-usual surface temperatures, for example during the Little Ice Age (Owens et al., 2017; Miller et al., 2012; Zanchettin et al., 2013; Schurer et al., 2014; Timmreck et al., 2021) and the so-called Late Antiquity Little Ice Age (LALIA, Büntgen et al., 2016; Toohey et al., 2016). The manifestation of volcanic radiative forcing as multi-decadal temperature anomalies has been suggested to be related to the thermal inertia of the Earth's oceans, which dampen the initial temperature response but also prolong it through the accumulation of negative energy input in the deep ocean (Gupta and Marshall, 2018). Indeed, ocean sea surface temperature changes over the 801–1800 CE period are well explained by the time series of volcanic eruptions (McGregor et al., 2015). There are climate-modeling studies that suggest that strong volcanic radiative forcing can perturb ocean circulation modes (Swingedouw et al., 2017), which may produce long-term perturbations to surface climate (Miller et al., 2014; Schleussner and Feulner, 2013; Schleussner et al., 2015; Otterå et al., 2010; Zhong et al., 2011). Indeed, at the global scale, clusters of eruptions have been linked to global mean temperature variations over the Common Era (Pages2k-Consortium, 2019). But other work has suggested a potential match in the timings of large eruptions or clusters of eruptions with the sudden transitions in climate between stadials and interstadials (Baldini et al., 2015; Bay et al., 2006; Lohmann and Svensson, 2022), although the robustness of this connection is limited by uncertainties in eruption magnitudes and timings.

Ice core sulfur (Huybers and Langmuir, 2009) and tephra chronologies (e.g., Praetorius et al., 2016; Sigvaldason et al., 1992) both attest to a marked increase in eruption frequencies during and after the last glaciation, especially in the Northern Hemisphere (NH) middle to high latitudes. On longer timescales, variations in eruption frequency from deep marine cores suggest periodicities on Milankovic timescales (Kutterolf et al., 2013, 2019; Schindlbeck et al., 2018b), suggesting a connection between the climate changes brought about by variations in orbital parameters and volcanic eruption frequencies. A leading theory is that mass transfer between ice sheets and the ocean due to changing temperatures leads to changes in the pressure on the Earth's crust, which can modulate crustal stress fields and enhance the possibility of associated eruptive events by providing pathways for the magma to rise (Kutterolf et al., 2013).

Recent ice core-based reconstructions of volcanic sulfur emissions confirm the increase in explosive eruptive frequency after the last deglaciation and its likely impact on stratospheric aerosol levels (Sigl et al., 2022). Ice core data presented by Lin et al. (2022) extend into the last glacial period and suggest that, during glacial conditions, the frequency of large eruptions with significant volcanic stratospheric sulfur injection (VSSI) was similar to that after the deglaciation.

The emissions of subaerial eruptions are eventually deposited to the Earth's surface and, in some cases, are preserved, providing us with records of past volcanic activity. Tephra layers from past eruptions are often discernable in terrestrial outcrops, but the records are often incomplete since younger deposits overlay the stratigraphy, or the deposits are already heavily weathered, eroded, or covered by vegetation. Ash and tephra, either transported by fallout or density flow processes, are well preserved in the marine sediments since wide areas of the seafloor are relatively unaffected by erosion or bioturbation (Freundt et al., 2021). This makes marine sediments outstanding archives for previous eruptions. By piston or gravity coring (∼ up to 20 m depth) or drilling (several hundreds of meters) the marine sediment records can be recovered, and the marine ash layers can provide a stratigraphically controlled record. While drilling is time consuming and expensive, and while the lateral coverage is thus limited to a few sites, multiple shorter piston or gravity cores can be taken in a certain region, providing good coverage, which enables detection of lateral stratigraphic changes or local erosion. Due to their length, however, they are often limited to the last two glacial cycles. In general, marine tephra records from sediment cores can cover several millions of years depending on the used coring or drilling technique. There are geographical gaps or regions with sparse coverage of cores, and, by far, not all cores have been studied for their tephra inventory in detail.

Figure 1 Global overview map created using GeoMapApp (http://www.geomapapp.org, last access: March 2023; GMRT-Global Multi-Resolution Topography) (Ryan et al., 2009). Red triangles mark global distribution of volcanoes from the LaMEVE database (Crosweller et al., 2012). Black squares mark regions with marine sediment cores that have been used in this publication. Red circles show positions of ice cores in Greenland and Antarctica used by Cole-Dai et al. (2021) and Sigl et al. (2015, 2022, and references therein).

Ice cores (Fig. 1) provide a good archive for past volcanic eruptions as the volcanic sulfate aerosols and, in some cases, ash particles are deposited to the surface and incorporated into the glacial ice. Chemical analyses of ice cores therefore provide time series, which are especially valuable when the dating of the ice cores is of high quality (Hammer, 1977; Gao et al., 2008; Sigl et al., 2015). Past volcanic eruptions are evidenced by strong increases in the sulfur or sulfate content of the ice or by the electrical conductivity of the ice as the analysis traverses the ice core and thus reaches backwards in time. By synchronizing multiple ice cores from Greenland and Antarctica, recent efforts have produced estimates of the VSSI from volcanic eruptions covering the past 2500 years (Toohey and Sigl, 2017), the Holocene (Cole-Dai et al., 2021; Sigl et al., 2015), and the late glacial period and deglaciation (Lin et al., 2022).

Figure 2 Scatter plot showing the distribution of eruptions by magnitude over the last 140 000 years. Data comprises all eruptions from the LaMEVE database, as well as additional eruptions from marine sediment cores, as described in Sect. 2.1.

Due to limitations with ice cores' absolute length, thinning, and synchronization, there is a limit as to the length of time ice cores can be used to reconstruct past volcanism. Tephra records from sediment cores, on the other hand, extend much further back in time. However, they have their own limitations: incompleteness, dating uncertainty, and not being as direct a measure of the stratospheric sulfate aerosol load as ice core sulfate records are. There is a strong temporal trend in the data set, which is described by a decreasing number of detected events going back in time (Fig. 2). Brown et al. (2014), for example, emphasize in the LaMEVE database, which covers the last 1.8 Ma, about 40 % of the detected eruptions included occurred during the Holocene (the past 11 kyr), and they conclude that the decrease going back in time is mainly due to under-recording of eruptions. The time trend in underreporting is found to depend on the magnitude of eruption, with the frequency of smaller eruptions ($\mathrm{4.0}\le M<\mathrm{5.0}$) falling off much faster than that for larger eruptions (M>6) (Fig. 2). Despite these limitations, far-reaching tephra records do provide valuable information about specific events (e.g., Pinatubo, Toba) and changes in eruption frequency with time. For example, analysis of tephra records long enough to cover several glacial cycles shows that eruption frequencies vary by periods representative of Milankovitch cycles, supporting the claims that the Earth's climate influences eruption frequencies on long timescales through changes in sea level and associated crustal stresses (Paterne et al., 1990; Rampino et al., 1979; Kutterolf et al., 2013, 2019; Schindlbeck et al., 2018b).

Here we present a new time series of volcanic stratospheric sulfur injection called PalVol (Toohey and Schindlbeck-Belo, 2023), covering the last 140 kyr, that is based on terrestrial and marine tephra records and includes stochastically generated synthetic eruptions to correct for the incompleteness of tephra records, in particular regarding small to medium eruptions (M<6). We provide an ensemble of time series, with each realization providing different timings of events found in the tephra record according to the uncertainties in the dating of the events, as well as different timings of the synthetic events. Our aim is not to provide an accurate reconstruction of the actual timing and magnitude of eruptions over the last glacial cycle. This is potentially even impossible. Rather, our aim is to provide a plausible set of such time series, each of which might approximate a possible true history given the information available and some basic assumptions. While we do not guarantee the accuracy of the timing and magnitude of specific eruptions, we do aim to produce a time series which represents our best estimate of the stochastic forcing provided by volcanic eruptions, as well as some accuracy in terms of millennial-scale variability in volcanic forcing.

The paper is organized as follows: in Sect. 2, the data and methods used to produce the ensemble VSSI time series product are introduced. In Sect. 3 the ensemble VSSI product is compared to existing ice-core-based reconstructions. A discussion and conclusions are included in Sect. 5.

2 Data and methods

2.1 Tephra data

The majority of tephra data used in this study are extracted from the LaMEVE database (Crosweller et al., 2012; LaMEVE Version 3), a global compilation of large-magnitude eruptions with VEI ≥4 (volcanic explosivity index) and/or magnitude 4 that are known from terrestrial deposits and outcrops. The LaMEVE database is publicly available and summarizes information regarding (1) erupted mass or volume and therefore magnitude or VEI; (2) eruption dates, as well as the applied dating techniques and uncertainties; (3) the source volcano; (4) eruption parameters (e.g., column height); and (5) rock types. We focused on the time interval from 140 000 years BP to 2014 AD (comprising 1565 events).

The eruption data taken from LaMEVE are complemented and corrected by records from a suite of marine cores and recent studies (Fig. 1; Tables A1, A2 in the Appendix), including cores and samples from Central America, offshore the Izu-Bonin volcanic arc, offshore the Cape Verdes, offshore the Kamchatka peninsula, and from the Hellenic Arc. Eruptive volumes and magnitudes for eruptions with no published values have been calculated by applying a single isopach approach following the methods described in Schindlbeck et al. (2018a).

The tephra record is characterized by increasing incompleteness back in time (Papale, 2018; Kiyosugi et al., 2015; Brown et al., 2014) (Fig. 2). The reasons for missing eruptions in terrestrial tephra records are mostly erosion and alteration or burial by younger deposits (e.g., Lavigne et al., 2004; Pollard et al., 2003). However, generally, it is thought that the incompleteness is more apparent for smaller-magnitude eruptions (e.g., M=4) compared to larger magnitudes (e.g., M≥5) (Brown et al., 2014) since preservation in marine sediments and on land is strongly increasing with the thicknesses of deposits (Wetzel, 2009; Freundt et al., 2021) (Fig. 2). In the marine environment, increasing time for bioturbation and diagenesis (the time in which bio-organisms and chemical interaction with formation waters can actively modify the ash layers physically and chemically), as well as the depositional environment, plays a role in the incompleteness of the records (Hopkins et al., 2020; Freundt et al., 2021). Additionally, for long records, the plate motions (except for ocean island volcanoes) might be responsible for a decreasing number of recorded eruptions in the past since most of the eruptive records are coupled with convergent margins and associated arc volcanism (e.g., Schindlbeck et al., 2015). The deeper and older the ash layers in the marine sediments are, the greater the distance between the volcanic source and deposition location at the time of the eruption due to plate motion (e.g., Schindlbeck et al., 2015). This mechanism could therefore reduce the possibility that smaller eruptions are recorded in the marine sediment archive from convergent margins.

Figure 3 Dating uncertainties over time for eruptions included in the LaMEVE database over the X–Y period, color-coded according to the different dating techniques as provided by Crosweller et al. (2012) and Cisneros de León et al. (2021a, b). Dashed gray lines indicate the 1 % and 10 % error range.

Marine or lacustrine ash layers, as well as terrestrial deposits, can be either directly or indirectly dated. The techniques comprise radiometric techniques (e.g., radiocarbon (¹⁴C), ⁴⁰Ar ${/}^{\mathrm{39}}$Ar mineral dating, zircon dating) and orbital tuning of oxygen isotope curves and/or sedimentation rates (Drexler et al., 1980; Kutterolf et al., 2008), which lead to age uncertainties of approximately 1 %–10 % of the estimated age (Fig. 3). The most precise dates are obtained for tephra that can be associated with historical observations of an eruption or with ice core or tree ring signals. The current literature has updated several eruption ages, which have not been included in PalVol v1 but which will be included in the next version. These comprise, in particular, ages obtained by ice core or dendrochronological studies, including Bárdarbunga, 877 CE (Plunkett et al., 2023); White River Ash, 853 CE (Mackay et al., 2022); Ilopango, 431 CE (Smith et al., 2020); Okmok II, 43 BCE (McConnell et al., 2020); Aniakchak II, 1628 BCE (Pearson et al., 2022); and Laacher See, 13 006 BP (Reinig et al., 2021).

Magnitude (M) is the preferred measure of eruption size used in the LaMEVE database (Crosweller et al., 2012) and is the quantity used here to estimate VSSI (Newhall and Self, 1982). Magnitude is a function of erupted mass, which is typically derived from estimates of erupted volumes, classically calculated by drawing isopach maps (areas of equal thickness of the ash layer; e.g., Bonadonna and Houghton, 2005; Fierstein and Nathenson, 1992; Pyle, 1995) and integrating over area. For eruptions used in this study which are not included in LaMEVE, magnitudes are taken from the respective publications. For the volcanic eruptions described in Derkachev et al. (2020), no volume estimates were available, and we therefore calculated minimum volumes by applying a single-isopach approach (Schindlbeck et al., 2018a) based on the given thickness data. The estimated uncertainty in mass estimates (and thus magnitude) is of the order of 10 % to 20 % (e.g., Kutterolf et al., 2021; Klawonn et al., 2014) and depends on the sample size (available outcrops on land and core density in the sea) and uncertainties in rock densities used for converting volumes into masses.

2.2 Ice-core-based VSSI

VSSI time series derived from polar ice cores are used in this work, both as a basis for the PalVol VSSI reconstruction methods (see Sect. 2.3, 2.4) and to validate the new reconstruction through statistical comparisons (Sect. 3).

For the −500 to 1900 CE period (256 eruptions), VSSI is taken from the eVolv2k reconstruction (Toohey and Sigl, 2017). These estimates are based on a set of continuous sulfate records from a suite of ice cores from Greenland and Antarctica. As in earlier ice-core-based reconstructions (e.g., Gao et al., 2008), tropical eruptions are identified by simultaneous deposition in both Antarctica and Greenland, and signals present in only one hemisphere are assumed to result from extratropical eruptions. VSSI is estimated by applying empirically derived transfer functions to the ice sheet average sulfate flux values (Toohey and Sigl, 2017; Gao et al., 2008).

Volcanic stratospheric sulfur injection estimates for the Holocene (from 9500 BCE or 11 500 years BP to 1900 CE, comprising 1496 eruptions) are available from the HolVol v1.1 reconstruction (Sigl et al., 2022). The method of VSSI reconstruction is very similar to that of eVolv2k, although the number of ice cores used is necessarily smaller since fewer cores cover the full Holocene. Despite the lower number of samples, the HolVol reconstruction shows good agreement with the 2500-year-long eVolv2k record, strengthening confidence in its accuracy (Sigl et al., 2022). Uncertainty in the timing of eruptions in HolVol is estimated to be ±1 to 5 years on average over the last 2500 years and better than ±5 to 15 years for the rest of the Holocene. Reported uncertainties in the VSSI values in HolVol for explosive (i.e., non-effusive) eruptions are typically between approximately 20 % and 40 % (Sigl et al., 2021).

Ice-core-derived VSSI estimates have been recently reconstructed for the 60–9 ka period, which covers the late glacial period, as well as the Early Holocene (Lin et al., 2022). Due to the thinning of ice sheets with age, the reconstruction of Lin et al. (2022) is limited to eruptions which produced strong deposition to both Greenland and Antarctica. The reconstruction provides estimates of stratospheric sulfate loading for 85 eruptions with bipolar deposition. To convert the sulfate aerosol mass estimates of Lin et al. (2022) to units of mass sulfur, we divide by a factor of 3 to account for the ratio of molar masses for sulfur (32 g mol⁻¹) to sulfate (96 g mol⁻¹).

2.3 Deriving VSSI from tephra

VSSI is estimated assuming a linear relationship between VSSI and erupted volume, i.e., a power-law relationship between VSSI and eruption magnitude, as used by Pyle et al. (1996) and Metzner et al. (2014). Here, we derive a fit of VSSI to magnitude using ice-core-derived VSSI (Sect. 2.2) and tephra-based magnitudes (Sect. 2.1) for identified eruptions, as well as recent eruptions for which estimates of sulfur emission are available from satellite instruments (Carn, 2022) (Fig. 4).

Figure 4 Volcanic stratospheric sulfur injection (VSSI) derived from ice cores (Toohey and Sigl, 2017; Sigl et al., 2021) and satellite observations (Carn, 2022) as a function of eruption magnitude from the LaMEVE database (Crosweller et al., 2012). A least-squares power-law best fit is shown in black, compared to similar fits from Pyle et al. (1996) and Metzner et al. (2014). A 1σ uncertainty range to the fit is shown as gray shading.

Table A1 lists eruption data which are used to derive a relationship between VSSI and eruption magnitude. For the period 1980 to 2014, we use sulfur emission estimates from satellite instruments compiled by Carn (2022). Emitted sulfur is matched to the eruptions listed in LaMEVE for this period. For the period before the satellite era, we rely on VSSI estimates and eruption attributions included in the eVolv2k reconstruction (Toohey and Sigl, 2017), supplemented with a few events from the HolVol v1 reconstruction (Sigl et al., 2022). Compared to prior investigations of the VSSI–magnitude relationship (Pyle et al., 1996; Metzner et al., 2014), this data compilation includes more satellite-based estimates, as well as ice-core-based estimates, extending the coverage to larger-magnitude eruptions, particularly due to the inclusion of the recently attributed ice core signals of Okmok (McConnel et al., 2020) and Crater Lake (Sigl et al., 2022) (Fig. 4).

2.4 Semi-stochastic VSSI time series generation

In order to correct for the incompleteness of the tephra-based eruption time series, we manufacture a synthetic eruption time series with the same statistical characteristics as an input data set, with timing and magnitudes of eruptions randomized. To construct our supplementary synthetic eruption time series, we draw on the efficient methodology of Bethke et al. (2017) but using as a basis the recent HolVol ice-core-based reconstruction over the years −4000 to 1900 CE. The algorithm of Bethke et al. (2017) stochastically produces a new eruption time series based on the eruption magnitude–frequency distribution of the input time series. Our time range is chosen so as to include a large period to enhance the statistical basis of the frequency distribution while excluding the early to mid-Holocene (11 to 6 ka), for which the eruption frequency is amplified compared to the middle to late Holocene, likely due, at least in part, to changes in ice loading at the high latitudes. When an eruption is randomly generated, characteristics of the eruption from HolVol, including VSSI, eruption region (tropical, NH extratropical, or Southern Hemisphere (SH) extratropical), and month, are copied into the constructed synthetic time series. For eruptions that are unidentified in the HolVol base data (which is the case for the majority of events), eruption parameters of month and precise latitude (within three latitude ranges) are unknown and set to default values. In our synthetic eruptions, the eruption month is randomized uniformly across the calendar year, and the eruption latitude is randomized within the identified tropical, NH, and SH bands using the probability density of the LaMEVE data set between 10 ka and the present. We repeat this process 100 times to produce 100 synthetic eruption time series.

To produce an ensemble of final VSSI time series therefore requires merging each synthetic time series – based on the statistics of the HolVol ice core data – with the evidence from the tephra record. Merging the two eliminates the decreasing eruption frequency backwards in time in the tephra record, assuming that this characteristic of the tephra record is a product of incompleteness. It also assumes that the true eruption frequency distribution is approximately static with time. For each of the 100 ensemble members, the synthetic and tephra records are merged so that each event in the tephra record is inserted into the synthetic record while also removing from the synthetic record a synthetic event with the closest matching magnitude within a window of 500 years. To represent the dating uncertainty in the tephra-based events in the ensemble of forcing time series and also to avoid clumping of eruptions around intervals of thousands of years due to the limited resolution of the reported dates of the tephra-based events, thus potentially creating an artificial millennial-scale periodicity in the radiative forcing, we add a random, normally distributed perturbation to the reported date of the tephra-based eruptions based on the estimated dating uncertainty (see Sect. 2.1 and Fig. 3). Thus, the dates of tephra-based eruptions in our reconstruction match the original dates for an eruption within the reported uncertainty. The dating perturbation is performed separately for each of the 100 ensemble members so that the date assigned differs in each realization of the data product and so that the spread in dates between realizations depends on the dating uncertainty of the tephra event.

2.5 Aerosol optical properties

The impact of volcanic aerosol radiative forcing is implemented in climate models at different levels of complexity. The simplest models take as input variations in the top-of-atmosphere radiative flux anomalies (W m⁻²), which represent the net effect of scattering and absorption of radiation by stratospheric aerosol (e.g., the energy balance model used in Pages-2k-Consortium, 2019). Comprehensive climate models, on the other hand, may require detailed optical properties of the aerosol as a function of latitude, height, wavelength, and time, to be used in the radiative calculations of the model (as, e.g., for the Max Planck Institute Earth System Model used in Bader et al., 2020).

To produce time series of the radiative impacts of past eruptions, including the detailed optical properties required by comprehensive models, we use the Easy Volcanic Aerosol (EVA) forcing generator (Toohey et al., 2016). This simple model takes as input the eruption timing, VSSI estimates, and eruption latitude and produces the aerosol extinction, single scattering albedo, and scattering asymmetry factor as a function of latitude, altitude, wavelength, and time. These variables are the result of a simple three-box model of stratospheric transport; scaling approximations between aerosol mass and AOD at 0.55 µm; and Mie theory, which describes the scattering of radiation by spheres. The overall impact of stratospheric aerosol on the Earth's energy balance is roughly proportional to the aerosol optical depth (AOD) in the visible part of the electromagnetic spectrum so it is common to illustrate the volcanic forcing simply as the AOD at 0.55 µm, either as a function of latitude or as a global (area weighted) mean.

2.6 A note on date formatting

Throughout the paper, we use two conventions with regard to dates. For periods extending no further than around 2500 years into the past, we use a variation of the ISO8601 format, which is very similar to the usual Common Era (CE) system, though it differs in the sense that the ISO8601 system includes a year 0, while the Common Era system does not. The two systems therefore differ by 1 year for years before 1 CE; for example, 44 BCE would correspond to the year −43 in the ISO8601 system.

For dates further in the past, we use the widely used Before Present (BP) system, which indicates the number of years before 1950.

3 Results

3.1 VSSI–eruption-magnitude relationship

The relationship between eruption magnitude and VSSI as observed by satellite instruments over the most recent decades and estimated from ice cores over the last 2500 years and Holocene is shown in Fig. 4. Satellite observations predominantly offer information on the VSSI from M=4 to M=5 eruptions: the only eruptions larger than M=5 that have been observed directly are the 1991 Pinatubo (M=6.1) and Cerro Hudson (M=5.8) (Crossweller et al., 2012, and references therein). For M=4 eruptions, the VSSI observed by satellites covers a range from 0 to approximately 1 Tg S. The largest values of VSSI observed from satellites are associated with the eruptions of Pinatubo (M=6.1, VSSI = 7.6 Tg S) and El Chichón (M=5.1, VSSI = 4.0 Tg S). Eruptions from the ice core reconstructions extend the range of eruption magnitudes included in the analysis; here, we use 24 eruptions, of which 21 have magnitudes greater than or equal to 5.0.

The VSSI values from satellite instruments and ice core reconstructions taken together show a proportionality with eruption magnitude: larger eruption magnitudes lead generally to larger VSSI values. Following Pyle et al. (1996), we fit a power-law relation to the data to obtain a best-fit relationship for VSSI (in Tg S): VSSI $=(\mathrm{1.67}\times {\mathrm{10}}^{-\mathrm{5}})\times {\mathrm{6.27}}^M$.

This fit is very similar to fits produced by Pyle et al. (1996) based solely on satellite data from the period 1979–1993 and is also somewhat similar to the fit presented by Metzner et al. (2014), who based their fit on petrologically obtained sulfur emission values from Central American eruptions (see Fig. 4).

There is clearly a significant amount of scatter around the line of best fit. This scatter can be the result of uncertainties in the eruption magnitude estimates, as well as the VSSI values, but is likely to be dominated by actual variability in sulfur emissions between different eruptions with similar magnitudes (e.g., Andres et al., 1993; Sigurdsson, 1990) and of the proportion of emitted sulfur that is injected into the stratosphere. The strong apparent scatter for M=4 eruptions in Fig. 4 represents a modest scatter in term of absolute numbers, with VSSI ranging from approximately 0 to 1 Tg S, reflecting variations in the sulfur output, as well as the modulating effect of plume heights, which, in many cases, are not high enough to bring sulfur into the stratosphere. For M≥5, the proportionality between VSSI and magnitude is more compact, with the majority of events falling within approximately 1 order of magnitude around the best-fit line. Notable outliers for which the VSSI is larger than the best-fit-line relationship include Laki (1783), for which the Greenland VSSI estimate has been suggested to potentially include a significant amount of tropospheric aerosol (Lanciki et al., 2012), and Hekla (1766), another Icelandic eruption for which the same scenario may hold. In contrast, a VSSI much smaller than expected, based on the best-fit relationship, is seen for the Millennium eruption of Changbaishan (940 CE), which has been discussed previously and may be due to some combination of sulfur-poor magma (Horn and Schmincke, 2000) and short stratospheric lifetime due to injection height or seasonal atmospheric dynamics (Iacovino et al., 2016). Accordingly, VSSI estimated from eruption magnitudes should be understood to have significant uncertainty for any individual eruption.

3.2 Common Era

During the Common Era (0–2000 CE) the eruption rate as recorded in the tephra records is larger than that of the ice core records during the chosen base period (−4000 to 1400 CE). Therefore, we can base the PalVol reconstruction purely on the tephra records without any synthetic events that could otherwise be included to compensate for undersampling by the tephra (Fig. 5).

Figure 5 VSSI estimates for the Common Era from (a) tephra and (b) ice cores (eVolv2k, Toohey and Sigl, 2017). (c) Centennial total VSSI for both data sets, with the correlation of the two centennial time series.

The time series of VSSI calculated from the tephra magnitude estimates shows reasonable agreement with the ice-core-based eVolv2k reconstruction (Fig. 5). The Rinjani (Samalas) eruption of 1257 CE, which produced the largest VSSI in eVolv2k (59.4 Tg S), is well reproduced in the tephra time series, with an estimated VSSI of 69.4 Tg S. The large eruptions of Tambora (1815) and Krakatau (1883) are apparent in both time series, albeit they are larger by a factor of approximately 2 in the tephra-based estimates compared to eVolv2k. The tephra record includes a large VSSI associated with the Changbaishan eruption of 942, which is not found in ice cores, which is likely due to the sulfur-poor content of the erupted magma (Horn and Schmincke, 2000). Other prominent events in the first millennium CE of the tephra VSSI time series, including TBJ (450±30 CE, VSSI = 38.9 Tg S), Taupo (230±16 CE, VSSI = 57.2 Tg S), and Ambrym (190±135, VSSI = 47.2 Tg S), are likely to be responsible for ice core signals within the dating uncertainty of the tephra events.

Given the large uncertainties in VSSI estimated from tephra (and the non-negligible uncertainties in the same value from ice cores), we would not and should not expect the VSSI values for individual eruptions to agree to a precision better than 1 order of magnitude. When multiple eruptions are averaged together, we expect the errors to average out to some degree and the cumulative VSSI values to become more reliable. Centennial total VSSI is shown for the tephra reconstruction and eVolv2k in Fig. 5c. The correlation between the two centennial time series is notable, with an overall correlation coefficient of R=0.50. The correlation is especially good in the second millennium, with both reconstructions showing elevated mean VSSI amounts for the 13th and 19th centuries, as well as more modest elevated values for the 15th and 17th centuries. While the tephra time series misses the elevated VSSI values of the 6th century, there is some agreement in the elevated VSSIs during the 2nd and 3rd centuries.

Figure 6 Comparing tephra and ice-core-based VSSI time series over the Holocene. (a) VSSI time series derived from tephra records. (b) HolVol ice-core-based VSSI time series (Sigl et al., 2022). Gray shading shows the period used to base eruption statistics for the synthetic time series.

3.3 Holocene

VSSI reconstructed from tephra is shown in Fig. 6 for the Holocene period (roughly 11 to 0 ka) and compared to the HolVol ice-core-based reconstruction. Both data sets show an increase in large-magnitude VSSI in the early Holocene, around 11 to 7 ka. Particularly, both data sets include four quite large events between 9 and 7 ka, with the tephra-based values reaching approximately 100 Tg S, while the ice core values reach up to over 150 Tg S.

Figure 7 Millennial-scale variations in eruptive characteristics over the Holocene in ice core (HolVol), tephra, and PalVol time series. (a) Number of volcanic events per millennium in each data set, (b) cumulative VSSI per millennium, (c) average VSSI per event, (d) histogram of correlation coefficients calculated between the ensemble of semi-synthetic millennial cumulative VSSI time series with the HolVol time series. Each bar indicates the fraction of all ensemble members with correlation coefficients between the values defining the edges of the bar along the horizontal axis. In (d), the correlation between the tephra time series and HolVol is indicated by the orange cross. In panels (a)–(c), gray lines indicate values for each of the 100 stochastic realizations, and the black line indicates the ensemble mean.

In a next step we compare the statistical characteristics of the tephra data, the PalVol reconstruction ensemble, and the HolVol ice-core-based VSSI reconstruction over the period of the Holocene on millennial timescales. Comparing first the number of events in the tephra and ice core data sets (Fig. 7a) we see that the number of tephra events drops rapidly with increasing age, from ∼275 events per kyr from 0 to 1 kyr to less than 100 events per kyr before 2 kyr and then to less than 50 per kyr in the early Holocene, which is consistent with prior analyses of the LaMEVE database (Brown et al., 2014).

From roughly 0 to 9 ka, the ice core number of events was held approximately steady at around 100 events per kyr, which increases during the deglaciation period, reaching a maximum of ∼160 events per kyr around 11 ka. In the most recent 2 ka, the frequency of tephra events is larger than the frequency of ice core events. This is likely due in part to overcounting of tephra events due to overestimation of the magnitude of events (increasing the absolute number of VEI ≥4 events included in LaMEVE). We speculate that another source of discrepancy may be undercounting of VEI ≥4 eruptions by ice cores since not all eruptions may leave traces in ice cores if not explosive enough or sulfur poor. Before 0 CE, the incompleteness of the tephra record is clear compared to the ice cores, and the degree of incompleteness increases roughly linearly until around 6 ka. There may be a weak local maximum in tephra events around 9 to 10 ka in agreement with the ice core increase here, but the difference in tephra events between this period and the local minimum at 7 to 8 ka is rather small.

Next, we compare the cumulative VSSI per millennium derived directly from the tephra data compared to the HolVol VSSI database (Fig. 7b). The tephra-derived VSSI per kyr is, in all millennia before 2 ka, smaller than the HolVol values, which is perhaps unsurprising given the much smaller number of events in the tephra database compared to the ice cores. In contrast, over the Common Era (0–2 ka), the tephra VSSI per kyr is larger or roughly equal to the ice core VSSI per kyr.

Over the Holocene, the average VSSI per millennium from HolVol is 638 Tg S, while that from the pure tephra data is 231 Tg S per millennium.

Despite clear differences between the tephra and ice core VSSI reconstructions for specific events and a bias with tephra showing lower values for most millennia of the Holocene, there is q correlation between the two data sets in terms of millennial VSSI over the Holocene, with both showing larger VSSI per kyr values in the early Holocene (11 to 7 ka) compared to the mid-Holocene (2 to 7 ka). The correlation coefficient between the millennial VSSI totals for the tephra and ice-core-based data sets over the 10 to 2 ka period is R=0.46. This correlation comes about despite the strong undercounting of the tephra database in the early Holocene. Evidently, although the number of events captured by the tephra data sets is small, the events that are counted tend to be the larger eruptions, which contribute the most to the VSSI millennial sums. Figure 7c shows the average VSSI per event as a function of millennium, which shows a similar structure for both the tephra and ice core data, with larger average VSSI per event in the early Holocene compared to the middle and late Holocene. This implies that the increase in VSSI per kyr in the early Holocene arises from an increase in the frequency of large eruptions, which is evident in both the ice core and tephra data.

Now we compare the semi-stochastic PalVol VSSI reconstruction, constructed by adding stochastically generated events to the tephra data, to the statistics of the HolVol reconstruction. By construction, the PalVol reconstruction includes a relatively constant number of events per millennium, around 100 events per kyr (Fig. 7a), in agreement with the frequency of events in the HolVol reconstruction over the chosen base period of −4000 to 1900 CE (roughly 6 to 0 ka). The number of events per kyr for each individual ensemble member of the reconstruction will vary around this number based on the stochastic event generation, and the inter-millennial variance matches well with the variance of the HolVol reconstruction over the base period. The PalVol reconstruction clearly does not include the increase in events per kyr in the early Holocene that is seen in HolVol (Fig. 7a). This is by construction since we do not adjust the eruption frequency probability with time.

Nonetheless, in Fig. 7b, we see that after adding the stochastically generated events to create the PalVol forcing reconstruction, the agreement of the millennial distribution of VSSI per kyr with the HolVol reconstruction improves compared to the pure tephra time series. First, the overall bias is reduced (but not eliminated), with an average VSSI per kyr of 505 Tg S for PalVol compared to 639 Tg S for HolVol. Secondly, since more stochastically generated events are added to the early Holocene compared to the late Holocene to make up for the larger bias in event frequency (see Fig. 7a), the VSSI per kyr is boosted more in the early Holocene. This improves the correlation: the mean PalVol VSSI per kyr time series shows a correlation of R=0.73 with HolVol. A histogram of the correlation coefficient between each individual PalVol VSSI time series with HolVol (Fig. 7d) shows that the largest proportion of ensemble members have a correlation of 0.7–0.8. We conclude that, over the Holocene, the addition of stochastic events to the tephra data improves resulting time series in comparison to ice-core-derived reconstructions, both by reducing the low bias in the tephra data and improving the inter-millennial variability of VSSI.

Figure 8 Volcanic stratospheric sulfur injection for the 60–9 ka period from (a) PalVol and (b) the ice core-based reconstruction of Lin et al. (2022). Stochastically generated PalVol VSSI (in gray) are from a single ensemble member, while tephra-based estimates (red) are shown at the reported dates (not randomized).

3.4 Last glacial cycle

The PalVol VSSI reconstruction is compared to the ice-core-based VSSI reconstruction of Lin et al. (2022) in Fig. 8. The Lin et al. (2022) events used here constitute 85 eruptions with matched bipolar signals with a deposition of >20 kg km⁻² in Antarctica and >10 kg km⁻² in Greenland between 60–9 ka. Tephra events in this period have a median dating uncertainty of 2300 years (900–4600 years, 0.25–0.75 interquartile range); therefore, we do not expect clear temporal matches to the ice core events for the reported dates of tephra events. Despite this, we find a decent agreement in the VSSI frequency distribution of the largest events in both time series for the period overall. For example, we find 11 events with VSSI >100 Tg S in the Lin et al. (2022) time series, while in the PalVol record for the same 60–9 ka BP period, we have 10 such events. The largest VSSI signal in the PalVol reconstruction over this period is associated with the 27 ka Taupo Oruanui eruption with M=8.1, leading to an estimated VSSI of 480 Tg S, almost 5 times greater than the ice-core-derived value of 127 Tg S. This represents an overestimate compared to the ice-core-based reconstruction for this event, which is not overly surprising given the uncertainty in the tephra-based reconstruction method.

Figure 9 Tephra-based volcanic stratospheric sulfur injection (VSSI) reconstruction. (a) VSSI based on the pure tephra record. (b) A semi-synthetic VSSI time series based on merging a stochastic synthetic record with the tephra reconstruction. (c) Millennial VSSI rates for the tephra, PalVol, and HolVol reconstructions.

Reconstructed VSSI is shown in Fig. 9 for the full PalVol reconstruction period. Individual eruptions are shown in panels (a) and (b), for which the decreasing number of detected eruptions with age and the resulting increasing number of synthetic eruptions included are apparent. The four largest VSSI estimates in PalVol exceed the range plotted in Fig. 9: the corresponding VSSI values are listed in Table 1 as part of the top 20 VSSI estimates. The eruption of Toba is the largest eruption of the past 140 000 years and results in an estimated 3000 Tg S VSSI using our method, with a wide estimated uncertainty range spanning 310 to 29 000 Tg S. Prior estimates of Toba's sulfate emission fall within a very broad range; multiple studies were compiled by Oppenheimer (2002) to define a range of 35 to 3300 Tg S. More recently, Costa et al. (2014) estimated Toba's sulfur emission as 850–1750 Tg S, which falls within our uncertainty range, while Crick et al. (2021) estimated a range of 72–233 Tg S, which falls outside our uncertainty range. It must be stressed that our estimates of VSSI for the strongest eruptions are based on an extrapolation of the VSSI–magnitude relationship beyond what has been observed, and non-linearities in the physical processes (e.g., plume collapse) are not considered here.

The Los Chocoyos eruption of Atitlan, ∼75 ka, has the same magnitude as the Taupo eruption and therefore the same estimated VSSI of 480 Tg S, with a range of 57 to 4000 Tg S. This value is very similar to the 343 Tg S estimated by Metzner et al. (2014), although the method used there was very similar to that used here.

Our estimate for the Changbaishan eruption (946 CE) is 130 Tg S, while Horn and Schmincke (2000) estimated a release of 5.7 Tg S, and Iacovino et al. (2016) estimated a release of 45 Tg S. A lack of strong ice core sulfate signals around the documented time of the eruption has been used as evidence that the sulfur emission from Changbaishan must have been quite minor (Oppenheimer et al., 2017) and likely much lower than the estimate here based only on the eruption magnitude.

Table 1 Estimated VSSI for the top 20 largest eruptions of the past 140 000 years based on the LaMEVE database.

Figure 9c shows the millennial average VSSI for the PalVol reconstruction, as well as the pure tephra data and the HolVol ice core reconstruction. Each PalVol ensemble member is shown in gray: these reconstructions generally show millennial average VSSI ranging from approximately 0.3 to 0.7 Tg S yr⁻¹, consistent with the mean value and variability from HolVol. Exceptions occur in millennia which contain the strongest eruptions, for which the millennial average VSSI can increase by a factor of 2 or more. The ensemble average PalVol millennial VSSI is approximately constant with age except for the perturbations due to the largest eruptions and some small increases due to enhanced numbers of detected strong eruptions, for example in the 55–40 ka period.

Figure 10 The PalVol global VSSI product, with the coordinates of each circle representing the time and latitude of an eruption and the circle size representing the VSSI in Tg S. Circle colors signify the type of event, either detected events compiled in LaMEVE (red) or additional records from marine sediment cores (purple) or the stochastic synthetic events used here to fill the record (gray).

One realization of PalVol VSSI time series is visualized in Fig. 10, with VSSI magnitudes plotted as bubbles as a function of ka and latitude. Many characteristics of the LaMEVE database are apparent: the decreasing eruption sampling with increasing age and the dependence of this sampling on latitude, with a generally more complete sampling in the NH middle latitudes. This figure illustrates the distribution of the tephra-based and synthetic eruptions in time and with latitude and specifically the latitude distribution of the synthetic events based on the latitudinal distribution of the LaMEVE eruptions.

4 Data availability

The PalVol volcanic stratospheric sulfur injection stratospheric aerosol optical depth data sets described herein are available through the World Data Center for Climate in netCDF format (https://doi.org/10.26050/WDCC/PalVolv1; Toohey and Schindlbeck-Belo, 2023).

5 Conclusions and discussion

We have produced PalVol VSSI reconstructions that cover the last glacial cycle (the last 140 kyr). The ensembles incorporate tephra data of past eruptions, as well as stochastically generated events to attempt to correct for the undersampling of events with increasing age. VSSI values are derived from the eruption magnitude, which is itself estimated based on the thickness of tephra layers in the surrounding of a volcano. Importantly, the VSSI values for individual eruptions have significant uncertainties as we (and prior works; e.g., Andres et al., 1993; Sigurdsson, 1990) showed that the amount of sulfur released by an eruption can vary by orders of magnitude for any given eruption magnitude (Fig. 4). Especially for M=4 eruptions, the VSSI varies significantly since some reach the stratosphere and some do not (Fig. 4), which depends on their actual plume height but also on latitude. Uncertainties in VSSI are included in our reconstruction as upper and lower bounds of the VSSI for each eruption based on the uncertainty in our derived magnitude–VSSI relationship. We note as well that this relationship may not be constant in time, and, indeed, changes may occur, e.g., to changes in the height of the tropopause or in stratospheric circulation from glacial conditions to interglacials. This would be similar to predicted changes in future volcanic radiative forcing due to changes in the tropopause height due to climate change (Aubry et al., 2016). This is something we did not take into account in our reconstructions. Nonetheless, although the VSSI for any specific eruption is quite uncertain, we expect that, due to the observed relationship between VSSI and magnitude, this will be compensated for when averaged over a sufficiently large set of eruptions, and therefore the tephra data may contain information on variations in VSSI on long timescales.

Furthermore, the PalVol reconstruction shows good agreement with the ice-core-based HolVol VSSI time series in terms of millennial variations in cumulative VSSI but not particularly in terms of the timing and magnitude of individual events.

This agreement arises from the shared increase in cumulative VSSI in the early Holocene compared to the middle and late Holocene. In both data sets, this arises from the increased frequency of relatively large eruptions. Evidently, while the tephra time series is incomplete during this period, the sampling of large events is good enough to detect the increase in the amount of large eruptions, which translates into increases in millennial cumulative VSSI.

We assumed a constant eruption frequency distribution, but there is evidence that variations in eruption frequency are driven by changes in the mass distribution of ice sheets and respective sea level changes. The PalVol reconstruction includes a small increase in VSSI per ka in the early Holocene, which is qualitatively consistent with observations from the ice cores but with an amplitude which is actually smaller. Therefore, we propose that future iterations of stochastic forcing could take into account such variations, for example by using sea level reconstructions as a basis to estimate variations in eruption frequency.

The PalVol reconstruction takes the form of an ensemble of realizations, where each realization differs in terms of the timing and size of the stochastically generated events and the timing of the events taken from the tephra corresponding to the uncertainty in the tephra dating. While any single realization is unlikely to be a very accurate reconstruction of the true history of VSSI, it is quite possible that some realizations include aspects that are realistic. Therefore, the ensemble of VSSI time series represents a probability distribution of the probable forcing from volcanic eruptions over this period given the information we currently have from proxies.

We emphasize that users wanting the most accurate reconstruction of VSSI over the last glacial cycle could consider using a merged product, for example by concatenating the HolVol ice core time series with PalVol for the period, which occurs before the beginning of HolVol. Future improvements to the PalVol reconstruction will be possible with the addition of new tephra data, improvements in the dating, and magnitude estimates of volcanic events.

Appendix A

Table A1 Compilation of matched eruption magnitudes from LaMEVE with volcanic stratospheric sulfur injection (VSSI) estimates from satellite measurements or ice-core-based reconstructions. Volcano, eruption year, and magnitude taken from the LaMEVE database. VSSI (in Tg S) is taken from three sources. Indicated by last column: 1 – from satellite observations compiled by Carn (2022), 2 – the eVolv2k ice core reconstruction (Toohey and Sigl, 2017), 3 – the HolVol ice core reconstruction (Sigl et al., 2022).

Table A2 Additional data from a suite of marine cores and recent studies including cores samples from Central America, offshore the Izu-Bonin volcanic arc, offshore the Cape Verdes, and offshore the Kamchatka peninsula.

* Volumes are calculated for a 30° single isopach following the methods described in Schindlbeck et al. (2018b).