Volcanic stratospheric sulfur injections and aerosol optical depth during the Holocene (past 11 500 years) from a bipolar ice-core array

The injection of sulfur into the stratosphere by volcanic eruptions is the dominant driver of natural climate variability on interannual to multidecadal timescales. Based on a set of continuous sulfate and sulfur records from a suite of ice cores from Greenland and Antarctica, the HolVol v.1.0 database includes estimates of the magnitudes and approximate source latitudes of major volcanic stratospheric sulfur injection (VSSI) events for the Holocene (from 9500 BCE or 11 500 years BP to 1900 CE), constituting an extension of the previous record by 7000 years. The database incorporates new-generation ice-core aerosol records with a sub-annual temporal resolution and a demonstrated sub-decadal dating accuracy and precision. By tightly aligning and stacking the ice-core records on the WD2014 chronology from Antarctica, we resolve long-standing inconsistencies in the dating of ancient volcanic eruptions that arise from biased (i.e., dated too old) ice-core chronologies over the Holocene for Greenland. We reconstruct a total of 850 volcanic eruptions with injections in excess of 1 teragram of sulfur (Tg S); of these eruptions, 329 (39 %) are located in the low latitudes with bipolar sulfate deposition, 426 (50 %) are located in the Northern Hemisphere extratropics (NHET) and 88 (10 %) are located in the Southern Hemisphere extratropics (SHET). The spatial distribution of the reconstructed eruption locations is in agreement with prior reconstructions for the past 2500 years. In total, these eruptions injected 7410 Tg S into the stratosphere: 70 % from tropical eruptions and 25 % from NH extratropical eruptions. A long-term latitudinally and monthly resolved stratospheric aerosol optical depth (SAOD) time series is reconstructed from the HolVol VSSI estimates, representing the first Holocene-scale reconstruction constrained by Greenland and Antarctica ice cores. These new long-term reconstructions of past VSSI and SAOD variability confirm evidence from regional volcanic eruption chronologies (e.g., from Iceland) in showing that the Early Holocene (9500–7000 BCE) experienced a higher number of volcanic eruptions (+16 %) and cumulative VSSI (+86 %) compared with the past 2500 years. This increase coincides with the rapid retreat of ice sheets during deglaciation, providing context for potential future increases in volcanic activity in regions under projected glacier melting in the 21st century. The reconstructed VSSI and SAOD data are available at https://doi.org/10.1594/PANGAEA.928646 (Sigl et al., 2021).

1 Introduction

Volcanoes impose various hazards on our climate, societal and economic systems. By injecting large amounts of sulfur into the atmosphere and thereby reducing the insolation reaching the Earth's surface (Robock, 2000), volcanic eruptions have been identified as the main drivers of natural climate variability on interannual to decadal timescales. They have been responsible for numerous cooling extremes in the past 2500 years (Anchukaitis et al., 2012; Guillet et al., 2017; Luterbacher et al., 2016; McConnell et al., 2020a; Sigl et al., 2015; Stoffel et al., 2015; Tejedor et al., 2021; Toohey et al., 2019, 2016a), often promoting crop failures and famines (Büntgen et al., 2020, 2016; Helama et al., 2018; Huhtamaa and Helama, 2017; Gao et al., 2021; Luterbacher and Pfister, 2015; McConnell et al., 2020b; Raible et al., 2016). Earth's largest volcanic eruptions (Crosweller et al., 2012; Mason et al., 2004) since the emergence of human civilization were more than an order of magnitude larger than the largest eruptions (e.g., Pinatubo 1991) for which we have direct observational evidence of the resulting atmospheric radiative perturbations and associated climatic effects (Douglass and Knox, 2005; Graf et al., 1993; Kremser et al., 2016).

Although attribution studies have affirmed a pivotal role of volcanism in driving climate variability during the past 1000 years (Owens et al., 2017; Schurer et al., 2014), volcanic forcing estimates have rarely (to date) been included in comprehensive climate model experiments aiming at simulating the climate evolution over the Holocene (Braconnot et al., 2012; Harrison et al., 2014). Therefore, the extent to which changes in global (Huybers and Langmuir, 2009) or regional volcanic activity (Maclennan et al., 2002) and volcanic extreme events (Zdanowicz et al., 1999) during the Holocene influenced climate evolution on various timescales is unknown. Abrupt large-magnitude changes in temperature (Mayewski et al., 2004; Wanner et al., 2011) and hydroclimate (Bond et al., 1997; Donges et al., 2015) frequently occurred throughout the Holocene and cannot be fully explained by the mix of external forcing and feedbacks considered at present (Liu et al., 2014; Wanner et al., 2008). Climate model simulations that include volcanic forcing, however, produce hemisphere-wide centennial-scale to millennial-scale temperature variability in better agreement with proxy evidence (Dallmeyer et al., 2021; Kobashi et al., 2017).

The ice sheets of Antarctica and Greenland contain invaluable information regarding the role that volcanic eruptions have played in driving past variations in the Earth's climate. Thus, ice cores obtained from polar ice sheets are the primary archives for the reconstruction of volcanic activity and its associated atmospheric aerosol loading (Gao et al., 2008; Sigl et al., 2014; Zielinski, 1995). To date, robust reconstruction of the timing and sulfate injection of explosive (and effusive) volcanism based on multiple ice cores exists only for the period of the past 2500 years (Sigl et al., 2015). While several individual ice-core histories have been developed for Greenland (Zielinski et al., 1994) and Antarctica (Castellano et al., 2004; Hammer et al., 1997), their use for reconstructing global volcanic forcing has been limited until recently (Cole-Dai et al., 2021), owing to previously poorly constrained age models (Parrenin et al., 2007; Plunkett et al., 2022; Torbenson et al., 2015) and post-depositional processes (e.g., wind erosion) at some low-snow-accumulation sites in East Antarctica that are able to disturb the original deposition record (Gautier et al., 2016). Eruptive histories from Greenland, on the other hand, are strongly dominated by events from proximal volcanic activity, in particular from nearby Icelandic volcanism (Abbott and Davies, 2012; Clausen et al., 1997; Coulter et al., 2012; Sigl et al., 2013; Thordarson and Hoskuldsson, 2008; Thordarson and Larsen, 2007). Both of these limitations have, thus far, hampered the identification of stratospheric tropical eruptions over the Holocene.

Reconstruction of volcanic forcing requires that all individual ice-core records are synchronized to a common timescale, which is achieved using records of volcanic fallout (e.g., acidity, sulfur and sulfate) archived in the ice sheets (Parrenin et al., 2012; Seierstad et al., 2014; Severi et al., 2007; Sigl et al., 2014). Aligning the records is possible because volcanic aerosols from eruptions are well mixed in the stratosphere and quickly dispersed, often on a hemispheric to global scale (Robock, 2000; Toohey et al., 2013). As the reference chronology, we use the “WD2014” timescale based on annual-layer counting from the WAIS Divide (WD) ice core in Antarctica, which provides the highest absolute dating accuracy of the ice chronologies currently available (Sigl et al., 2016). Abbreviated names are used for the numerous ice cores utilized to reconstruct the volcanic chronology. A list of all of the abbreviations used in this paper is included in the Supplement (Table S1). The exceptional high resolution of the WD sulfur and sulfate records (Cole-Dai et al., 2021) and the increased dating precision (Sigl et al., 2016) enable us to conduct a firm bipolar synchronization with ice cores from Greenland over the Holocene, as has previously been demonstrated for the Common Era (Plummer et al., 2012; Sigl et al., 2013, 2015). Large volcanic eruptions from the lower latitudes resulting in the global distribution of sulfate over both hemispheres are recognized by synchronous sulfate deposition in Greenland and Antarctica, employing constraints provided by the high relative age precision of the two layer-counted chronologies in both hemispheres between subsequent volcanic marker events (Sigl et al., 2016; Vinther et al., 2006). The isotopic composition of sulfur in volcanic sulfate also contains information on whether aerosol formation occurred at altitudes above the ozone layer, allowing an independent test of the association of sulfate peaks and stratospheric eruptions with global sulfate distribution (Baroni et al., 2008, 2007; Burke et al., 2019; Gautier et al., 2019). Combining information from Antarctica and Greenland therefore enables us to disentangle likely source regions of volcanic eruptions (i.e., the Northern Hemisphere extratropics, NHET; the Southern Hemisphere extratropics, SHET; and low latitudes) which are important to analyze volcanic activity through time and to estimate their radiative forcing on past climate (Crowley and Unterman, 2013; Gao et al., 2008; Toohey et al., 2016b). The boundaries between these conceptualized likely source regions are understood to be permeable, with interhemispheric mixing of stratospheric sulfate aerosols also likely to occur after large eruptions in the extratropics (Aubry et al., 2020; Marshall et al., 2019; Toohey et al., 2013; Wu et al., 2017).

Volcanic stratospheric sulfur injections (VSSI) from global volcanic activity, summed over centuries, have varied by an order of magnitude between the highly active 13th century – marking the inception of the Little Ice Age – and the 1st century CE (Toohey and Sigl, 2017). Even larger variations have likely occurred during the warm Early Holocene, when the rapid melting of large ice sheets during deglaciation (Clark et al., 2012) regionally triggered a strong acceleration in volcanic activity (Maclennan et al., 2002; Sigmundsson et al., 2010; Watt et al., 2013) through feedback chains that may also operate during the 21st and 22nd centuries with projected changes in the cryosphere under global warming (Schmidt et al., 2013; Tuffen, 2010). Understanding how future volcanic activity may affect climate is strongly dependent on understanding the statistical nature of volcanic activity: its variability and the degree of temporal clustering of eruptions (Bethke et al., 2017; Man et al., 2021; Tuel et al., 2017).

2 Method

2.1 Ice-core sites

The drilling site for the WD ice core (79.48^∘ S, 112.11^∘ W; 1766 m a.s.l.) was selected to obtain a precisely dated, high-time-resolution ice-core record that would be the Southern Hemisphere equivalent of the deep Greenland ice cores (WAIS Divide Project Members, 2013, 2015). The 3404 m long WD ice core was collected from a cold (mean annual temperature −31 ^∘C), high-snowfall (200 kg m⁻² yr⁻¹) West Antarctic site. Within the European Project for Ice Coring in Antarctica (EPICA), more deep ice cores were drilled in Antarctica (EPICA-Community-Members, 2004, 2006). The ice core from Dronning Maud Land (EDML; 75.00^∘ S, 00.07^∘ E; 2882 m a.s.l.) has, at 68 kg m⁻² yr⁻¹, a 2–3 times higher accumulation rate than the ice core at Dome C (EDC; 75.10^∘ S, 123.35^∘ E; 3233 m a.s.l.). Multiple deep ice cores have also been drilled in Greenland, including the Greenland Ice Core Project (GRIP; 72.60^∘ N, 35.80^∘ W; 3232 m a.s.l.), the Greenland Ice Sheet Project 2 (GISP2; 72.58^∘ N, 38.47^∘ W; 3053 m a.s.l.) and the North Greenland Ice Core Project (NGRIP or NorthGRIP) ice cores, providing continuous records of atmospheric impurities over the Holocene (Seierstad et al., 2014). Figure S1 (in the Supplement) summarizes the depth–age relation for these deep ice cores on the common, annual-layer-counted WD2014 chronology (Sigl et al., 2016) after applying volcanic synchronization during the glacial (Buizert et al., 2018) and Holocene (this study). The specific datasets used for aligning these chronologies are shown in Table 1.

Table 1 Ice-core records and analytical methods used (ICPMS: inductively coupled plasma mass spectrometry; IC: ion chromatography; FIC: fast ion chromatography; DEP: dielectric profiling)

2.2 Ice-core measurements on WD

Sulfur concentrations between 1300 and 2003 m depth (4060–10 000 BCE, 6010–11 950 years BP) covering the Early to Middle Holocene were analyzed using trace element continuous flow analysis (TE-CFA) at the Desert Research Institute (DRI) in Reno, USA. The DRI Ultra-Trace Chemistry Laboratory used a method that allowed continuous, simultaneous measurement of a large number of trace elements at very high depth resolution (Cole-Dai et al., 2021; McConnell, 2002; McConnell et al., 2017, 2018). The depth resolution for sulfur achieved with this system is 1 cm in ice, allowing for the achievement of a nominal monthly time resolution over the entire Holocene. Sulfate concentrations between 577 and 1300 m depth (396–4060 BCE) covering the Middle to Late Holocene and the brittle zone (Neff, 2014) of the WD ice core were analyzed using ion chromatography in discrete and continuous flow analysis mode (Cole-Dai et al., 2021, 2013, 2006) at the South Dakota State University, USA. The depth resolution for sulfate was 2 cm. High sampling resolution throughout the Holocene permitted the detection of annual cycles in impurity data, allowing for the precise and accurate annual-layer dating of the ice-core records during the entire Holocene (Sigl et al., 2016). For consistency with Toohey and Sigl (2017), we report the calendar ages using the ISO 8601 international standard, which does (in contrast to the historical Gregorian calendar) include a year zero. For key events and time periods, we also report ages as years before present (BP, years before 1950) – a notation used frequently in archeology, geology and other scientific disciplines.

2.3 Volcanic synchronization

Synchronization is based on matching volcanic sulfate, sulfur, acidity or conductivity peaks of the dependent core to equivalent peaks in an independently dated reference core, and it is used to transfer or to synchronize ice-core timescales. It is widely used in the ice-core community to align ice-core chronologies on a common reference chronology (Langway et al., 1995; Sigl et al., 2014; Svensson et al., 2020). For Greenland and the Arctic, many ice cores (e.g., NGRIP, GRIP and GISP2) have been synchronized (Rasmussen et al., 2013; Seierstad et al., 2014) on the Greenland Ice Core Chronology 2005 – GICC05 (Rasmussen et al., 2006; Svensson et al., 2008; Vinther et al., 2006), whereas the WD2014 chronology (Buizert et al., 2015; Sigl et al., 2016) serves as the reference chronology for Antarctica (Buizert et al., 2021, 2018; Sigl et al., 2015; Winski et al., 2019). Ice cores have also been synchronized across the hemispheres (Langway et al., 1995), but the density and certainty of these match points have been much lower owing to hitherto low dating precision in ice cores from Antarctica and the large abundance of volcanic eruption signals in the Greenland ice core from high-latitude volcanic eruptions (e.g., Iceland, Alaska and the Kamchatka Peninsula) hampering reliable source attribution. In an attempt to synchronize ice cores from Greenland and Antarctica over the entire Holocene, a total of 74 match points have been suggested between the NGRIP and EDML ice cores (Veres et al., 2013), about as many as were identified between WD and EDML during the Common Era (Sigl et al., 2014).

The accuracy of stratigraphic matches further depends on the volcanic signal (e.g., temporal resolution) and ice-core site-specific properties (e.g., accumulation rate variability) of both the dependent and reference ice-core records through time. We synchronized ice-core records in this study using an iterative approach. First, volcanic signals with outstanding magnitudes and characteristic temporal spacing that are virtually certain (e.g., in the 17th century BCE, 2910 BCE, the 45th–43th century BCE and the 67th–63th century BCE) were synchronized (major tie points). Confidence in these match points was derived from the combination of (1) a temporal sequence of distinctive signals, (2) comparable magnitudes, (3) a uniform evolution of layer thickness between stratigraphic tie points, (4) a distinctive shape of the common signals in some cases and, finally, (5) independent age constraints from ¹⁰Be that reflected variations in cosmic ray flux (Adolphi and Muscheler, 2016; Sigl et al., 2016). Using linear interpolation of the derived initial mean annual layer thickness calculated between age markers, secondary stratigraphic links from moderate volcanic eruptions became obvious and were matched to WD2014 (see Fig. S2 in the Supplement). Relative accumulation rates in Antarctica calculated for longer time periods usually show low variability, narrowing the window for potential stratigraphic tie points between the two records. We applied volcanic synchronization against WD first to EDML and EDC and verified that the individual selected tie points were consistent with the previous volcanic synchronization between EDML and EDC (Severi et al., 2007). We performed several iterations allowing for 218 (EDML) and 148 (EDC) volcanic matches with WD (Fig. 1, Table 2). We repeated this approach for the two Greenland ice-core records of sulfate from GISP2 (Mayewski et al., 1997) and dielectric profiling (DEP) from GRIP (Wolff et al., 1997). Confidence in the match points was derived from the combination of (1) a distinctive sequence of common signals, (2) a uniform evolution of layer thickness between stratigraphic tie points, (3) sequential annual-layer counts between volcanic age markers and (4) constraints from ¹⁰Be matching. We performed several iterations allowing for 164 (GISP2) and 93 (GRIP) volcanic matches with WD (Fig. 2, Table 2). We verified that all major bipolar tie points identified in GRIP and GISP2 are consistent with the previous synchronization between GRIP and GISP2 (Seierstad et al., 2014).

Figure 1 Volcanic synchronization for Antarctica: panel (a) shows the changes in mean annual layer thickness ($\mathrm{d}z/\mathrm{d}t$) and age difference (WD2014 minus AICC2012 – Veres et al., 2013) calculated between volcanic tie points for the EDML and EDC ice cores, respectively; panel (b) shows the WD non-sea-salt sulfur (nssS) record as well as the EDML and EDC sulfate records for the 8600–8000 BCE time period and volcanic tie points; panel (c) shows the same records for the 4850–4250 BCE time period. All records are shown at an annual resolution on the annual-layer-counted WD2014 chronology (Sigl et al., 2016).

Table 2 Number of volcanic tie pints identified between different deep ice cores.

This study (bold). ^a Sigl et al. (2014). ^b Toohey and Sigl (2017). ^c Severi et al. (2007). ^d Veres et al. (2013). ^e Seierstad et al. (2014).

Figure 2 Bipolar volcanic synchronization: panel (a) shows the changes in mean annual layer thickness ($\mathrm{d}z/\mathrm{d}t$) and age differences (WD2014 minus GICC05 – Vinther et al., 2006; WD2014 minus GISP2 annual-layer-counted timescale – Meese et al., 1997) calculated between volcanic tie points for the GISP2 ice core; panel (b) shows the WD non-sea-salt sulfur (nssS) record as well as the GISP2 sulfate record for the 4700–4200 BCE time period and volcanic tie points used for the synchronization; panel (c) shows the same records for the 8600–8000 BCE time period. All records are shown on the WD2014 chronology (Sigl et al., 2016).

2.4 Volcanic signal detection and sulfate mass deposition

Sporadic volcanic sulfate deposition at the ice-core sites is superimposed on the background deposition of marine sulfate and other sulfuric species (e.g., methanesulfonic acid). This background is seasonally variable; however, without volcanic input, it has very limited variability between years (Cole-Dai, 2010). Therefore, a method to differentiate between volcanic sulfur or sulfate and the nonvolcanic background requires quantification of the background and its variability (Traufetter et al., 2004). To detect and quantify volcanic sulfate deposition, we used established methods (Cole-Dai, 2010; Cole-Dai et al., 2021; Gao et al., 2007; Sigl et al., 2013, 2014), which are summarized below. Sulfate deposition over both polar ice sheets varies systematically with snow accumulation and may be further modified randomly by site-specific post-depositional effects (such as redistribution of snow through wind). To account for these random and systematic differences, we selected a stacking approach for Antarctica (for which three continuous ice-core records are available) similar to the one used for the past 2000 years (Sigl et al., 2014) and in accordance with previous work (Gao et al., 2008; Crowley and Unterman, 2013).

Figure 3 Holocene sulfur records from Antarctica: sulfur concentration records from the (a)  WD, (b) EDML and (c) EDC ice cores as well as (d) a stack (“Antarctica” or ANT12k) of all three records for the Holocene (10 000 BCE–2000 CE). Measured as sulfate, the EDML and EDC records are synchronized on the WD2014 chronology (Sigl et al., 2016), annualized and scaled to the WD record. The upper part from EDML is based on the 200 m long B40 ice core drilled at the same site in 2012 (Sigl et al., 2015). Signals from two large historic eruptions, Tambora (1815) and Samalas (1257), are marked.

2.4.1 Volcanic sulfate deposition in Antarctica

We first resampled and annualized the sulfate and sulfur concentration records by averaging all samples within a calendar year (WD and EDML) or by interpolation (EDC). To compare the relative magnitudes of sulfur deposition at the three ice-core sites over the past 11 500 years, we scaled the EDML and EDC sulfate concentrations (in ng g⁻¹) by 5.1 and 6.7, respectively, with the scale factors determined by matching average Holocene sulfate deposition with sulfur concentrations at WD. Thus, the scale factors account for differences in the molar masses of sulfur (32 g mol⁻¹) and sulfate (96 g mol⁻¹) as well as differences in accumulation rates and emission sensitivities between the ice-core sites. Therefore, the resulting sulfur time series of EDML and EDC can be interpreted as the equivalent sulfur concentration at the WD site, allowing for the construction of an annually resolved sulfur concentration stack by averaging the three ice-core records (ANT12k, N=3; Fig. 3). We use this time series to assess the plausibility of the age synchronization based on the relative agreement of peak amplitudes as an additional diagnostic criterion (see Fig. S3 in the Supplement). The choice of alternative methods (such as standardization or normalization) has no significant influence on the results of this reconstruction (see Fig. S4 in the Supplement). As we do not know a priori which ice core best represents the stratospheric sulfate burden after volcanic eruptions, we use an unweighted average of all three ice cores.

We also employed this stack (in addition to the individual WD sulfur record) to synchronize the Greenland GISP2 sulfate record to the WD2014 chronology by identifying synchronous sulfate deposition in Antarctica and Greenland. The nonvolcanic background sulfur concentration was first estimated in ANT12k using the 101-year running median (RM) of the annually averaged sulfur data. The mean absolute deviation (MAD) from the RM was then determined for each 101-year window, which is a robust measure of background variability in the presence of outliers. To detect volcanic events against the variable background, a threshold of RM + 2 × MAD was set, which is comparable to previous work in Antarctica (Sigl et al., 2014; Gao et al., 2008). A year was classified to contain volcanic sulfur if the annual sulfur concentration exceeded this threshold. After removing all years with concentrations above this threshold, the reduced running mean (RRM) was calculated for the remaining years in the 101-year window of the time series. The duration of the volcanic event is defined as the length of time in which the sulfate concentrations exceeded RM + 1.5 × MAD. The annual volcanic sulfate concentration is calculated as the difference between the total sulfur concentrations of that year and the RRM of the nonvolcanic sulfate of that year. The cumulative sulfate mass deposition (kg km⁻²) from an eruption, often referred to as the (cumulative) “volcanic sulfate flux” f(volc-SO${}_{\mathrm{4}}^{\mathrm{2}-}$), is the sum of the annual volcanic sulfate concentrations in the years when volcanic deposition occurred multiplied by the mean annual accumulation rate at WD (210 kg m⁻² yr⁻¹). Finally, we scaled the cumulative sulfate flux from ANT12k against a corresponding area-weighted composite sulfate deposition rate obtained from a more comprehensive “AVS2k” stack including more than 10 ice cores (see Fig. 4; Sigl et al., 2014) using the relation f(volc-SO${}_{\mathrm{4}}^{\mathrm{2}-}$)${}_{\mathrm{ANT}\mathrm{12}\mathrm{k}}=\mathrm{1.769}\times f$(volc-SO${}_{\mathrm{4}}^{\mathrm{2}-}$)_AVS2k (R²=0.93, N=105) to estimate the ice-sheet average sulfate fluxes for Antarctica, which is henceforth referred to as f(volc-SO${}_{\mathrm{4}}^{\mathrm{2}-}$)_AVS12k. As the start date of the volcanic eruption, we use the initial [nssS] increase from WD which provides the highest temporal resolution and the lowest degree of peak broadening (through wind drift and snow mixing) typical for low-accumulation ice-core sites.

Figure 4 Representativeness of individual and stacked ice-core records in Antarctica: panel (a) displays the mean annual sulfur concentrations from the “Antarctica” stack (ANT12k, N=3) over the Holocene compared to the “AVS2k” stack from, on average, 14 ice cores from Antarctica over the Common Era (Sigl et al., 2014); the time period from 1550 to 2000 CE (purple shading) is displayed in panel (b), with known large volcanic eruptions from the tropics highlighted; panel (c) presents a map of the ice-core sites from Antarctica used in this and a previous study; and panel (d) displays a scatterplot of cumulative volcanic sulfate mass deposition (“flux”) for individual ice cores (WD, black; EDML, red; EDC dark red), the “Antarctica” stack composite record (light blue) and the “AVS2k” stack. Included in the analysis are the 30 largest sulfate deposition events in “AVS2k”.

Figure 5 Prolonged eruptive episodes in Greenland ice cores. Panel (a) shows the GISP2 total sulfate and volcanic sulfate records and the GRIP DEP record between 8400 and 8200 BCE; panel (b) shows the same records between 7600 and 7400 BCE; panel (c) shows the same records between 3250 and 3050 BCE; panel (d) shows the same records between 800 and 850 CE as well as mean sulfate concentrations from a stack of three synchronized ice cores from Greenland (NGRIP, NEEM-2011-S1 and TUNU2013) on the NS1-2011 chronology (Sigl et al., 2015). Shading indicates the time periods and duration of prolonged volcanic activity (Hjartarson, 2003; Sinton et al., 2005); stars mark tephra from Icelandic sources identified in ice cores from NGRIP, GRIP, NEEM, GISP2 and TUNU2013. Panel (e) shows the GISP2 sulfate record for the Holocene with a 15-year running median.

2.4.2 Volcanic sulfate deposition in Greenland

For Greenland, we followed a similar approach with some adjustments owing to the different properties of the available volcanic proxy records. For GISP2, only a single ice core with continuous sulfate measurements exists with a biannual temporal resolution (Zielinski et al., 1994), hampering the detection of smaller and short-lived volcanic perturbations (Toohey and Sigl, 2017). Stronger decadal to multidecadal background variations are observed, reproduced by shorter ice-core records (e.g., NGRIP and NEEM-2011-S1) and electrical records (e.g., DEP from GRIP), which we attributed to long-lasting volcanic episodes from Iceland (Fig. 5). Therefore, we tagged all GISP2 volcanic sulfate values exceeding the volcano detection threshold for a minimum of 10 consecutive years as a “prolonged eruption” (Table 3) and applied an additional correction to estimate sulfur injection (see Sect. 2.5).

The nonvolcanic background sulfate concentration was initially approximated in GISP2 with a 121-point (window) RM fit to the biannual sulfate data. This is equivalent to a 240-year median and, thus, better suited for detecting decadal to multidecadal volcanic sulfate variability than shorter window lengths. Similar to the approach used for Antarctica, we detected volcanic events that exceeded a threshold of RM + 1.5 × MAD, and samples were deemed to contain volcanic fallout if the sulfate concentration exceeded this threshold. After the removal of all years with concentrations above this threshold, the RRM was computed for the years that remained in the moving 121-point window of the time series. The duration of the volcanic event is defined as the length of time in which the sulfate concentrations exceed RM + MAD. The annual volcanic sulfate concentration is calculated as the difference between the total sulfate concentration of that sample and the RRM of the nonvolcanic sulfate of that sample. Finally, the cumulative sulfate mass deposition flux is the sum of volcanic sulfate concentrations in the years in which volcanic deposition occurred multiplied by the mean annual accumulation rate at GISP (210 kg m⁻² yr⁻¹). We tested the performance of our detection criteria during the preindustrial 19th century and found that volcanic sulfate is detected during the same periods in which volcanic eruptions had previously been detected by other higher-resolution sulfate records from Greenland (see Fig. S5 in the Supplement). No false positive events were reconstructed from GISP2 in this test, and volcanic signals reconstructed from other ice cores and not detected in GISP2 were of small amplitude and duration. Based on this comparison, we conclude that volcanic eruptions comparable in strength (with respect to sulfur injection) to the Icelandic eruptions of Katla (1755, 1.2 Tg VSSI) or Hekla (1766, 2.5 Tg VSSI) are detectable in GISP2, providing a lower bound of the detection limit for Icelandic eruptions.

Table 3 Prolonged eruptions. All volcanic sulfate deposition signals lasting >20 years sorted by duration as well as the most recent signal persisting >10 years from the GISP2 on the WD2014 timescale. Start and end dates are shown in the before present (1950 CE) notation as well as in the Common Era (CE) and Before the Common Era (BCE) notation using the ISO 8601 standard, which includes a year zero.

^a The tephra in ice cores from multiple Greenland ice cores indicates the Icelandic eruption of Grímsvötn (Saksunarvatn Ash) as a potential source contributing to the ice-core sulfate (Gronvold et al., 1995); ^b The tephra in an ice core from TUNU2013 Greenland indicates Katla (Iceland) as a potential source contributing to the ice-core sulfate (Büntgen et al., 2017; Plunkett et al., 2020). See Table S1 (in the Supplement) for a list of lava shield and fissure eruptions >2 km³ from Iceland following Hjartarson (2003) and Sinton et al. (2005).

2.4.3 Ice-core uncertainties

The timing of volcanic eruptions from ice-core records is uncertain due to interpretation uncertainties during the construction of the annual-layer dating. Based on the comparison of WD2014 with accurately dated tree-ring records (Sigl et al., 2015, 2016), we estimate that absolute age uncertainties in the ice-core records used in HolVol database are better than ±1–5 years on average over the last 2500 years and better than ±5–15 years for the rest of the Holocene. A 5000-year-long tree-ring record with strong sensitivity for abrupt post-volcanic cooling (Salzer and Hughes, 2007) allows for further assessment of the absolute age accuracy of WD2014 following some of the largest Late Holocene eruptions (see Sect. 2.8). Another source of uncertainty arises from the limited number of ice-core locations (i.e., one from Greenland and three from Antarctica) available to estimate the mean ice-sheet deposition and, ultimately, the hemispheric sulfate burden. We have previously estimated 1σ errors of 33 % with respect to estimating the Greenland ice-sheet-wide average flux from the mean sulfate flux of the single GISP2 record (Toohey and Sigl, 2017). We further assume 1σ errors of 20 % with respect to estimating the Antarctica ice-sheet-wide average flux from the mean sulfate flux of the AVS12k composite stack including three ice cores. The estimated total error of the mean for Antarctica is, therefore, slightly above the typical (root-mean-square) uncertainties of approximately 13 % for a larger (AVS2k, N=14) Antarctic ice-core composite but below a constant uncertainty value of 26 % based on regression analysis between AVS2k and the composite of WD and B40 over the 1–2000 CE period (see Sigl et al., 2015; Toohey and Sigl, 2017).

2.5 Injection locations and dates

Over the last 2500 years, the localities and the timing of several (N=31; Toohey and Sigl, 2017) stratospheric sulfur injections reconstructed from ice cores could be assigned based on matching the ice-core inventory with observed historical eruptions using the Volcanoes of the World online database (Global Volcanism Program, 2013) and other sources of information. There is some degree of uncertainty and subjectivity associated with such matchings. However, for certain cases, geochemical analysis of tephra from ice cores has been used to establish or strengthen the matches, including Veiðivötn in 1477 CE (Abbott et al., 2021b), Samalas in 1257 CE (Lavigne et al., 2013), Changbaishan in 946 CE (Oppenheimer et al., 2017; Sun et al., 2014), Eldgjá in 939–940 CE (Oppenheimer et al., 2018; Zielinski, 1995), Mt Churchill in 853 CE (Jensen et al., 2014), Katla in 822 CE (Büntgen et al., 2017; Plunkett et al., 2020) and Ilopango in 431 CE (Smith et al., 2020). Attributing locations to ice-core eruption signals over the full Holocene is even more difficult due to the increasing incompleteness and decreasing dating precision (often based on radiocarbon dating) over time of the volcanic eruption inventory derived from proximal geological evidence (Brown et al., 2014; Crosweller et al., 2012). To date, only a handful of ice-core sulfate peaks in the Holocene have been linked geochemically to known eruptions, including the caldera-forming 43 BCE Okmok II eruption in Alaska (McConnell et al., 2020a), the ca. 1628 BCE Aniakchak II eruption in Alaska (McAneney and Baillie, 2019; Pearce et al., 2004; Pearson et al., 2022; Plunkett and Pilcher, 2018), the caldera-forming 5677±150 BCE Crater Lake “Mazama” eruption in Oregon (Zdanowicz et al., 1999), the 5922±50 BCE Khangar eruption on the Kamchatka Peninsula (Cook et al., 2018) and the ca. 10 ka Grímsvötn “Saksunarvatn Ash” eruption series from Iceland (Oladottir et al., 2020). The vast majority of large eruptions, such as the caldera-forming Bronze Age Thera/Santorini eruption or the ca. 6440±25 BCE caldera-forming Kurile Lake eruption on the Kamchatka Peninsula, has remained unidentified in the ice-core record. We assigned approximate eruption latitudes for most sulfate signals in ice cores that cannot be immediately attributed to a known eruption, based on the presence or absence of simultaneous signals in the Greenland and Antarctic ice cores. Volcanic sulfate deposition identified synchronously (within small possible dating errors) in both Greenland and the Antarctic composites are attributed to eruptions in the tropics, whereas signals that occur in only one hemisphere are assumed to be of extratropical origin, as described in Sigl et al. (2015). Characteristic latitudes for unidentified eruptions are inferred from the latitudinal distribution of known eruptions. Using the mean distribution of all VEI ≥4 (where VEI refers to the volcanic explosivity index) eruptions from a Holocene eruption database (Global Volcanism Program, 2013), we assigned average latitudes of 48^∘ N, 37^∘ S and 5^∘ N to all unidentified extratropical eruptions in the Northern and Southern hemispheres and in the tropics, respectively. We further attributed all volcanic events for which volcanic sulfate deposition to Greenland persisted for more than 10 years to Icelandic source eruptions, most likely from the Katla, Bárðarbunga and Grímsvötn volcanic systems or from shield volcanoes in the Western Volcanic Zone (Hjartarson, 2003; Sinton et al., 2005; Thordarson et al., 2003), and assigned 64^∘ N as the default latitude for these prolonged episodes (Fig. 5; Tables 3, S1). We are aware that the use of default latitudes paints a crude picture of the geographical distribution of past global volcanic activity (see Fig. 6), but we currently lack the necessary knowledge to more precisely assign individual volcanoes to ice-core signals. We also note that, in the construction of aerosol optical properties and radiative forcing using the Easy Volcanic Aerosol (EVA) forcing generator (see Sect. 2.7), only the broad region of the eruption site is important (i.e., tropics, NHET and SHET), as the exact latitude has no impact on the generated aerosol properties. When sulfur emissions are directly used in aerosol–climate models, differences in aerosol evolution depending on the latitude of the eruption within these broad regions may be relevant (see Toohey et al., 2019; Marshall et al., 2021), and our choice of the mean latitudes helps to minimize any potential bias in the long-term mean radiative forcing. Consistent with Toohey and Sigl (2017), an eruption date of 1 January is assigned to unidentified eruptions.

Figure 6 Spatiotemporal distribution of volcanic stratospheric sulfur injections from volcanic eruptions since 9500 BCE from HolVol 1.0 based on known and assigned locations (Iceland, 64^∘ N; NHET, 48^∘ N; tropics, 5^∘ N; SHET, 37^∘ S). Only eruptions with VSSI >1 Tg S are included; prominent historic and prehistoric eruptions are marked; source attributions for Aniakchak and for Crater Lake are based on the geochemistry of cryptotephra from Greenland ice cores (Coulter et al., 2012; Zdanowicz et al., 1999).

2.6 Stratospheric sulfate injection estimation

Stratospheric sulfate injections are estimated from the ice-sheet sulfate flux composites using a method described in detail by Toohey and Sigl (2017). Briefly, ice-sheet average sulfate fluxes for Antarctica f_A and Greenland f_G are related to injected sulfur mass M_S following Eq. (1):

$$\begin{array}{}\text{(1)}& M_{\mathrm{S}}={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\left[L_{\mathrm{G}}{f}_{\mathrm{G}}+L_{\mathrm{A}}{f}_{\mathrm{A}}\right],\end{array}$$

where L_G and L_A are transfer functions accounting for the spatial distribution of sulfate deposition over each hemisphere. Based on the analysis of the spread and deposition of nuclides from a nuclear bomb test, sulfate from prior volcanic eruptions and atmospheric model simulations (Gao et al., 2007), the transfer functions L_G and L_A are estimated to be 1×10⁹ km² for tropical eruptions and 0.57×10⁹ km² for extratropical eruptions. The method described above is based on the assumption that the ice-core sulfate deposition is proportional to the stratospheric sulfur emission. In fact, some of the sulfate deposited may originate from volcanic sulfur emissions into the troposphere, especially when volcanic eruptions are situated upwind of (e.g., in Alaska) or in close proximity to the Greenland ice sheet (e.g., in Iceland). Recently, sulfur isotopes from Greenland ice-core records have been used to detect the presence of sulfate deposited in Greenland via both the stratospheric and tropospheric transport pathways following the large VEI = 6 eruption of Katmai/Novarupta in Alaska, upwind of the Greenland ice sheet (Burke et al., 2019). Of particular importance are long-lasting, effusive (i.e., nonexplosive) eruptions from Iceland, which may produce significant sulfate deposition over Greenland, even when the stratospheric injection is minimal. The two largest fissure eruptions in Iceland in historical times (Eldgjá 939–940 and Laki 1783–1784 AD) are the most prominent examples, and the extent to which sulfate was injected into the lower stratosphere during these eruptions is the subject of ongoing research (Lanciki et al., 2012; Schmidt et al., 2012; Zambri et al., 2019). Only for very recent volcanic eruptions are direct observations of key eruption source parameters (e.g., plume height, SO₂ dispersion height, duration, season) available, which determine how much sulfur gets injected into the stratosphere. Detailed volcanological fieldwork could delineate 10–11 distinctive eruptive episodes during the Laki 1783–1784 event (Thordarson and Self, 2003), allowing for the development of detailed SO₂ emission scenarios for modeling the climatic impact of this episode (Schmidt et al., 2010; Zambri et al., 2019). However, such detailed information is not available for other large fissure eruptions in Iceland, of which at least 14 are known to have occurred over the Holocene (Thordarson and Larsen, 2007; Thordarson et al., 2003). To correct for a significant proportion of tropospheric sulfate when estimating stratospheric sulfur emissions, Crowley and Unterman (2013) adjusted Greenlandic sulfate depositions following the Laki eruption in 1783–1784 and derived a ratio of stratospheric to total sulfate deposition of 0.15. Due to a lack of data on the stratospheric vs. tropospheric distribution of injected sulfur for other major Icelandic eruptions of the Holocene, we adopted this approach and used a transfer function of 0.10×10⁹ km² for those long-lasting deposition signals that we assume to be from prolonged volcanic eruptions in Iceland. We implemented the sulfur injections within these eruptive episodes using the biannual resolution of the GISP2 ice-core record (i.e., a 16-year-long episode is implemented as eight subsequent injections), and the durations are reported so that injection can be spread uniformly over time in simulations. We stress that additional objective criteria to detect proximal eruption signals, correctly attribute these to specific source eruptions and subsequently correct the VSSI estimates are urgently needed.

Estimates of VSSI have significant uncertainty due to three major sources of potential errors: (1) random errors in the ice-core flux measurements, (2) uncertainties in the transfer functions used to translate the ice-core sulfate data to estimates of VSSI and (3) potential errors in the estimation of the latitudinal position (and explosivity) of the eruption (i.e., tropical vs. extratropical explosive vs. extratropical effusive). VSSI uncertainties are included in the HolVol dataset, with the aim of estimating the uncertainties from the first two terms. Uncertainty related to the limited number of ice cores and related sampling of the ice sheets has been estimated (see Sect. 2.4.3). As in Toohey and Sigl (2017), this uncertainty is added in quadrature to an estimate of the uncertainty related to using ice-sheet deposition to estimate hemispheric deposition. Based on an ensemble of aerosol model simulations (Toohey et al., 2013), this term is estimated to contribute ∼16 % and 9 % of uncertainty to the Northern Hemisphere (NH) and Southern Hemisphere (SH) transfer functions (L_NH and L_SH), respectively, but these estimates may be model-dependent, and recent work has pointed to potentially larger values (Marshall et al., 2019, 2018). It remains difficult to quantify errors arising from a potentially incorrect attribution of the source location for individual eruptions. VSSI from an eruption erroneously attributed to a tropical source, which in reality may have been from two different eruptions in the high latitudes of both hemispheres, will be overestimated by 43 %. As another example, sulfate deposition in Greenland resulting from a potential cluster of several subsequent volcanic eruptions in the NHET may not be recognized as separate eruptions in the biannual-resolution GISP2 record; thus, it may be erroneously attributed to a prolonged eruptive period when sulfate levels remain increased for >10 years. Under such a scenario, the true VSSI would be underestimated by up to 80 %. To account for the latter case, reported VSSI uncertainties for prolonged eruptions have been inflated to include magnitudes that would be calculated if the eruption was not prolonged, which results in uncertainties of over 100 %. This large error also signifies a relatively low confidence in the adjustment to the transfer function used for prolonged eruptions compared with explosive extratropical eruptions. We note that only specific eruptions may be subject to errors caused by incorrect attributions which can be subsequently assessed and corrected in future updates of this database if independent constraints for source locations from cryptotephra, sulfur isotope and trace metal analyses of archived and new ice cores become available (Burke et al., 2019; Gautier et al., 2019; McConnell et al., 2017, 2020a). A primary source of systematic error in the VSSI estimates is likely to originate from the uncertainty in the transfer functions (L_G and L_A) used to estimate hemispheric stratospheric sulfate aerosol burdens. These are originally derived using ice-sheet sulfate fluxes in Antarctica and observed from the stratospheric sulfate burden of the Pinatubo 1991 eruption as well as deposited nuclear bomb test fallout in Greenland and climate model simulations (Gao et al., 2007). Continued efforts to constrain observational uncertainties with aerosol model simulations have been unsuccessful due to significant inter-model differences – for example in the simulation of aerosol spread and deposition after the Tambora 1815 eruption (Marshall et al., 2018).

2.7 Aerosol optical depth estimation

The Easy Volcanic Aerosol (EVA) version 1.2 forcing generator (Toohey et al., 2016b) is employed to convert sulfur emissions into optical properties of volcanic aerosols. We specifically consider the variation in the stratospheric aerosol optical depth (SAOD) at 550 nm. Using a time series of VSSI and eruption latitudes as input, EVA generates aerosol optical properties as required for use in climate model simulations. The spatiotemporal structure of the EVA output fields is based on a simple three-box model of stratospheric transport that is optimized to produce agreement with observations of the aerosol cloud from the eruption of Pinatubo in Indonesia in 1991. Internally, EVA first computes the transport of sulfate mass and then scales the sulfate mass to the SAOD. While this scaling is linear for most eruptions, following Crowley and Unterman (2013), a nonlinear scaling between mass and SAOD is adopted for very large eruptions (i.e., eruptions with a VSSI in excess of that of Tambora in 1815). Furthermore, to account for the self-limiting effect of aerosol growth on the stratospheric lifetime of aerosol after large eruptions implied in model studies (Pinto et al., 1989; Timmreck et al., 2009), a simple parameterization of variable removal time has been implemented in EVA based on ECHAM5-HAM aerosol model simulations of eruptions with a wide range of magnitudes (Metzner et al., 2014). Based on the model results, the stratospheric aerosol removal timescale is varied between its nominal value of 11 months and a minimum of 6 months as the global stratospheric sulfate burden rises above 10 Tg S. We refer to the SAOD results presented below that were generated by the EVA forcing generator using the HolVol VSSI database as “EVA(HolVol)”. This naming convention emphasizes the two-stage procedure of the SAOD reconstruction, with HolVol used as an input to EVA. SAOD time series are shown as monthly, annual or centennial averages. Therefore, peak SAOD values can differ significantly depending on the time resolution of the time series.

2.8 Assessment of dating accuracy and precision

Nominal age uncertainty for the WD2014 chronology due to ambiguities in the interpretation of annual layering has been estimated to linearly increase with age over most of the Holocene (Sigl et al., 2016). Constrained at 775 CE using ¹⁰Be in ice cores (Mekhaldi et al., 2015) and ¹⁴C in tree rings (Büntgen et al., 2018) to detect the distinctive 774/775 CE solar proton event (Miyake et al., 2012), the age error from annual-layer interpretation down to 3000 BCE (5 ka) is estimated to be better than ±20 years. During the Early Holocene, at 9500 BCE (11.5 ka), the WD2014 age uncertainty was estimated at ±66 years. Matching the common production signal in cosmogenic isotopes (¹⁴C, ¹⁰Be) has further allowed us to assess the WD2014 ages relative to the radiocarbon calibration curve, which is based on dendrochronology during the Holocene and, thus, has virtually no age uncertainty (Sigl et al., 2016). The best fit necessary to align both chronologies had been found to vary by small margins of less than ±15 years throughout the Holocene, suggesting that the cumulative error estimated from the annual-layer counting of WD2014 is very conservative. There is a tendency for WD2014 ages to be slightly too young during the Early Holocene and slightly too old between 7000 BCE and 1 CE (Sigl et al., 2016). No ¹⁰Be measurements from WD were previously available; thus, no assessment of the WD2014 timescale was possible for the time period between 3500 and 500 BCE. To fill this gap, we employ a multi-millennial compilation of the occurrence of ring-width minima and frost rings in a bristlecone pine chronology from the southwestern USA covering the past 5000 years (Salzer and Hughes, 2007). A strong association between frost-ring formation and climatically effective volcanic eruptions has been previously noted (Baillie, 2010; Lamarche and Hirschboeck, 1984; McAneney and Baillie, 2019; Salzer and Hughes, 2007; Sigl et al., 2015). To assess the temporal relation between major volcanic eruptions reconstructed with HolVol and cooling extremes indicated by the tree-ring series, we extract all marker events from the compilation by Salzer and Hughes (2007) (N=10) in which at least two consecutive ring-width minima corresponded to a frost-damaged ring within an error margin of ±1 year (Table 4). Additional marker years in which a frost ring corresponds to a ring-width minima within ±1 year are provided in the Supplement (Table S2).

Table 4 Dating assessment using tree rings. Marker events are shown in which a ring-width minima (Salzer et al., 2014) corresponded to a frost-damaged ring within an error margin of ±1 year (Salzer and Hughes, 2007) in relation to reconstructed volcanic deposition events over the Late Holocene (this study) and the past 2500 years (Toohey and Sigl, 2017). WD2014 ages are provided for bipolar eruption signals (Sigl et al., 2016). Ages from attributed Northern Hemisphere extratropical (NHET) eruptions are on the NS1-2011 chronology (Sigl et al., 2015). Eruptions with VSSI >10 Tg (comparable to Krakatau 1883) within ±5 years of the cooling start are highlighted in bold. All ages are reported using the ISO 8601 international standard, which does (in contrast to the historical Gregorian calendar) include a year zero.

^a Data gap in GISP2; VSSI is only based on Antarctica assuming a SHET source eruption and may be underestimated if a comparable large sulfate anomaly is detected in Greenland ice-core records. ^b A period with long-lasting reductions in ring width and frequent frost-ring appearance following a large tropical eruption in 682 CE (Table S2). n/a: not applicable.

Figure 7 Holocene volcanic stratospheric sulfur injection (VSSI) from explosive eruptions: panel (a) shows the reconstructed VSSI for single eruptions over the Holocene, panel (b) shows the number of eruptions per century and panel (c) shows the total VSSI per century. A version of this figure with VSSI shown separately for the three major source regions – NHET (30–90^∘ N), SHET (30–90^∘ S) and tropical (30^∘ N–30^∘ S) eruptions – is provided in the Supplement (Fig. S6).

3 Results

3.1 Volcanic stratospheric sulfur injections

The HolVol v.1.0 VSSI time series is shown in Figs. 6 and 7. With 100 % data coverage for Antarctica and 95 % data coverage for Greenland, we consider this record to be virtually complete for all volcanic eruptions with a strong climate impact potential (i.e., VSSI comparable to or larger than the 1991 Pinatubo eruption). The ability to detect and quantify smaller events is primarily limited by data gaps (equivalent to 560 years) and the coarse (biannual) temporal resolution of the GISP2 record from Greenland. Therefore, there is the possibility that smaller eruptions situated in the NHET are under-recorded. With only 88 % data coverage between 3000 and 1000 BCE obtained within the “brittle zone” of the GISP2 record, under-recording and ambiguities in matching and correctly attributing source latitudes pose some limitations at the moment.

Figure 8 Number of eruptions and cumulative volcanic stratospheric sulfur injection (VSSI) from eVolv2k (Toohey and Sigl, 2017) and from HolVol for the period of overlap (500 BCE–1900 CE) and the full Holocene reconstruction. Only eruptions with VSSI >1 Tg S are included.

HolVol v.1.0 contains a total of 1189 volcanic eruptions, resulting in a cumulative VSSI of 7412 Tg S between 9500 BCE and 1900 CE. On average, a detected eruption occurred once every 10 years. Of these eruptions, 850 injected at least 1 Tg S into the stratosphere, the equivalent of the eruptions of Nabro (Eritrea 2011) and Kasatochi (Alaska 2008), which were implicated to have contributed to the slowdown of warming in the 21st century (Carn et al., 2016; Ridley et al., 2014; Santer et al., 2014). Figures 6 and S6 show the latitudinal distribution for each reconstructed VSSI. Figure 8 summarizes their mean distribution over the Holocene. A total of 40 % of the eruptions (with VSSI >1 Tg S) are attributed to tropical eruptions (30^∘ N–30^∘ S), 5 % of eruptions are attributed to effusive Icelandic events (>63^∘ N), 45 % of eruptions are attributed to other NHET (30–60^∘ N) events and 10 % of eruptions are attributed to SHET (30–90^∘ S) events. The mean frequency of reconstructed volcanic eruptions >1 Tg S is 0.074 yr⁻¹ (i.e., an eruption every 14 years on average). These 850 eruptions injected a total of 7260 Tg S into the stratosphere, of which 75 % was emitted in the tropics, 4 % was emitted in Iceland, 18 % was emitted in the NHET and 4 % was emitted in the SH (Fig. 8).

The number of eruptions and the cumulative centennial VSSI varied within the Holocene (Figs. 6, 7). In general, the number of eruptions and the cumulative VSSI were enhanced during the Early to Middle Holocene (10th to 5th millennium BCE; 76 Tg S per century on average). Between the 4th millennium BCE and the present, both the average number of eruptions and the cumulative VSSI (45 Tg S per century) were 21 % and 41 % lower, respectively (Figs. 6, S4, S7). The period from 9500 to 7000 BCE, when glacial ice sheets were retreating rapidly and were widespread (Carlson and Clark, 2012), is characterized by the highest frequency of eruptions as well as the largest cumulative VSSI over the entire Holocene. With an average of 90 Tg S injected in the stratosphere per century, this period – which we term the “Deglaciation Active Period” – is 43 % above the Holocene mean VSSI rate of 63 Tg S. An increased VSSI rate is noted for eruptions in the NHET as well as in the tropics, but this feature is absent for the SHET (see Fig. S6 in the Supplement). The majority of the events that we attributed to prolonged eruptive episodes (379 years of cumulative duration) also falls into this time period (Fig. 6, Table 3).

The window from 4000 to 1000 BCE has the lowest frequency of eruptions and the smallest VSSI rates in the Holocene. With 21 % less eruptions and 36 % smaller VSSI rates than the Holocene mean, we term this period the “Holocene Quiet Period”, in analogy to other time periods with reduced volcanic activity such the “Medieval Quiet Period” (700–1100 CE) and the “Roman Quiet Period” (40 BCE–200 CE). Throughout the Holocene, the longest subsequent time period without an eruption with VSSI >1 Tg S is 77 years (ending in 3206 BCE); the longest period without VSSI >10 Tg is 317 years (ending in 2258 BCE). These volcanically quiet periods are slightly longer compared with the same metrics for the Medieval Quiet Period (55 and 217 years, respectively) and the Roman Quiet Period (71 and 212 years, respectively).

A total of 8 of the 10 largest VSSI injections are recorded between 6700 and 4300 BCE, all exceeding the VSSI of the largest known volcanic eruptions of the Common Era except for Samalas (Lavigne et al., 2013; Vidal et al., 2016) in 1257 CE (ranked ninth). The highest recorded VSSI reach values >150 Tg S. With a VSSI rate 22 % above the Holocene mean, we term this period the “Mid-Holocene Active Period”.

3.1.1 Comparison with other Holocene volcanic reconstructions

A limited number of previous reconstructions of volcanic sulfate injections exist for the Holocene. Based on the Camp Century ice core in Greenland, global acid fallout was estimated from a total of 18 eruptions between 8000 and 50 BCE (Hammer et al., 1980). A direct comparison on an event basis is not possible owing to the different chronology compared with this study. Age uncertainties of ±170 years in the Camp Century record are an order of magnitude larger than our estimates for HolVol. The only unambiguous match with HolVol is the large sulfate anomaly dated in Camp Century to 50±30 BCE, which has recently been pinned to the caldera-forming Okmok II eruption in 43 BCE in Alaska using tephra in the GISP2 ice core (McConnell et al., 2020a). The estimated equivalent global sulfur fallout (assuming that all acids were from H₂SO₄) from Camp Century was 40 Tg S, which is slightly below the estimate of 56 Tg S in HolVol. The highest global volcanic fallout estimates in Camp Century were 85 Tg S using latitudinal correction functions that assumed a high-latitude eruption source for the vast majority of the ice-core signals. These are significantly smaller than the highest estimates in HolVol of up to 190 Tg S, which were eruptions with bipolar sulfate deposition. A more comprehensive reconstruction of volcanic sulfate deposition was performed using the GISP2 ice-core record since 7000 BCE (Zielinski et al., 1994). In the GISP2 record, a total of 298 eruptions were detected in the residual volcanic sulfate. Using a less conservative volcano detection threshold (aided by a larger number of now available ice-core records during the past 2 kyr), we detect (for the same time period and in the same GISP2 sulfate dataset) a total of 555 eruptions. Age uncertainties in the GISP2 ice core were previously estimated at ±2 % of the age or approximately ±150 years some 8 ka before present (Meese et al., 1997). Zielinski et al. (1994) did not estimate sulfur injection nor changes in SAOD from the GISP2 record, but recent studies have estimated volcanic forcing from the GISP2 record (Bader et al., 2020; Brovkin et al., 2019; Kobashi et al., 2017). In the absence of a well-synchronized ice-core sulfate record from Antarctica at the time, these studies have assumed that all sulfate signals in Greenland were from eruptions located in the low latitudes. As a result, these reconstructions under-record eruptions from the SHET, and they are prone to systematically overestimate the forcing from Icelandic eruptions and many other NHET eruptions by at least 43 % and up to a factor of 10 for specific events.

3.1.2 Comparison with eVolv2k

The eVolv2k volcanic eruption catalogue (500 BCE–1900 CE) was reconstructed from bipolar ice-core records using a similar methodology to that used for HolVol v.1.0, but it was based on a larger number of sulfur (and sulfate proxy) records, including 3 from Greenland and up to 14 from Antarctica (Toohey and Sigl, 2017). Thus, eVolv2k remains the recommended volcanic forcing for transient climate model simulations covering the past millennium or the past 2 kyr, including experiments (Jungclaus et al., 2017) within the “Paleoclimate Modelling Intercomparison Project” (PMIP) contributing to the fourth phase of the PMIP (PMIP4). Ice-core records from the same sites employed by HolVol were also used in eVolv2k, which explains the strong similarity in the underlying Antarctica sulfur stacks (Fig. 4). Using HolVol, we estimate (from the four ice cores) a cumulative VSSI of 1278 Tg S from 180 eruptions with >1 Tg S injection between 500 BCE and 1900 CE, which is only 10 % above the value estimated based on eVolv2k (see Fig. S8 in the Supplement). The source distribution of the eruptions is also virtually identical between the different reconstructions during the period of overlap. On an event basis, the agreement between HolVol and evolv2k is strongest for larger eruptions (i.e., >10 Tg S), whereas there is a larger scatter for eruptions with smaller VSSI (see Fig. S8 in the Supplement). In order to perform a seamless Holocene-long simulation with climate models, we recommend merging the VSSI or SAOD reconstructions from HolVol v.1.0 with those from eVolv2k (see Fig. S9 in the Supplement) at the year 500 BCE or 1 CE.

Figure 9 Global annual mean stratospheric aerosol optical depth (SAOD) from the EVA(HolVol) reconstruction. Years are shown using the ISO 8601 standard, which includes a year zero.

3.2 Stratospheric aerosol optical depth and radiative forcing

The global mean SAOD from the EVA(HolVol) reconstruction is shown in Fig. 9. The SAOD closely follows the spatiotemporal structure of VSSI in the Holocene, albeit with relatively less pronounced peaks for the largest eruptions due to the nonlinear parameterizations used in EVA. The global mean SAOD over the Holocene was 0.0153; SAOD over the NH (0.0182) was almost 50 % higher than that over the SH (0.0124). The global mean SAOD between 9500 and 4500 BCE was 48 % higher than between 4500 BCE and 1900 CE. The difference in SAOD between these two time windows was stronger (+57 %) when integrating over the NH (0–90^∘ N), whereas it was less pronounced (+21 %) when integrating over the SHET (30–90^∘ S). The largest global annual SAOD reached 0.85; the largest SAOD over the NHET reached 1.45 following the Crater Lake eruption (Oregon, USA). For comparison, the largest eruption during the Common Era, Samalas 1257 CE, is estimated in eVolv2k to have produced a global annual SAOD of 0.50; the largest nontropical eruption of the Common Era in 536 CE produced a NHET SAOD of 0.43. We stress that for such large eruptions, which are significantly larger than any eruption observed in the instrumental era, uncertainties in the SAOD should be understood to be large.

Figure 10 (a) Global mean radiative forcing (RF) from the EVA(HolVol) reconstruction (inverted axis) and from a reconstruction based on the GISP2 ice core (Kobashi et al., 2017). (b) Centennial mean RF for the two reconstructions.

The EVA(HolVol) reconstruction is compared to that of Kobashi et al. (2017) in Fig. 10. As the Kobashi et al. (2017) reconstruction contains estimates of radiative forcing (RF, in units of W m⁻²), the EVA(HolVol) SAOD values are converted to RF using the linear scaling (RF = −25 ⋅ SAOD) of Hansen et al. (2005), which has been employed in prior Intergovernmental Panel on Climate Change (IPCC) reports (Myhre et al., 2013). We note that several recent studies have suggested that consideration of rapid adjustments (e.g., in cloud formation) leads to a reduction in the scaling factor in the order of 20 % (Marshall et al., 2020; Schmidt et al., 2018). The major difference between the multi-ice-core HolVol reconstruction and the single-ice-core (GISP2) reconstruction from Greenland is the smaller magnitudes (minima of −21 W m⁻²) of RF for large volcanic eruptions in HolVol compared with values as strong as −50 W m⁻², as reconstructed by Kobashi et al. (2017), which we attribute to applying a nonlinear scaling to HolVol. Furthermore, HolVol RF values are constrained by ice-core records from Antarctica, whereas Kobashi et al. (2017) assumed that the GISP2 sulfate record from Greenland is representative of the global volcanic sulfate burden, thereby inevitably overestimating RF for all eruptions with unipolar sulfate distribution (e.g., eruptions from Iceland) or eruptions with a strong asymmetry of the sulfate burden in the NH. While the negative radiative forcing from large events is very likely overestimated by Kobashi et al. (2017), the negative RF from smaller and moderate eruptions (that are often not detected in the single-ice-core reconstruction from Greenland) are underestimated. Some of the difference is also due to the additional inclusion of a nonzero background SAOD in the EVA(HolVol) reconstruction. The effect of these methodological differences is that HolVol depicts smaller variability than the previous reconstruction of global RF (Kobashi et al., 2017).

Figure 11 Estimated age uncertainty, and a comparison of WD2014 and independent chronologies based on dendrochronology: panel (a) shows the filtered WD ¹⁰Be (blue) and ¹⁴C (red) tree-ring data on their respective timescales, and panel (b) shows the most likely time shift (red line, 2 kyr sliding window) for the highly significant correlations along with the 2σ uncertainty range (see Sigl et al., 2016, for details). Green circles mark the age difference between major reconstructed eruptions (VSSI >10 Tg) and cooling anomalies (i.e., co-occurrence of frost-ring and ring-width minima within ±1 year) in a 5 kyr bristlecone pine chronology from the southwestern USA. A complete list of the selected event years is given in Table 4 and in the Supplement (Table S2) and was extracted from the compilation by Salzer and Hughes (2007).

3.3 Dating accuracy and precision

The results of the assessment of dating accuracy and precision are summarized in Fig. 11 (see Table 4 and Table S1 in the Supplement for details). The previous assessment based on correlating multidecadal- to centennial-scale production rates in the cosmogenic radionuclides ¹⁰Be and ¹⁴C showed only slightly varying age-scale differences over most of the Holocene (Sigl et al., 2016). This indicates a high accuracy and precision of the WD2014 timescale, in contrast to other annual-layer-dated chronologies (e.g., GICC05) that consistently overestimated (via overcounting annual layers) ages throughout most of the Holocene (Adolphi and Muscheler, 2016; Muscheler et al., 2014). Age differences between ice-core-indicated sulfate deposition and tree-ring-indicated summer cooling extremes (i.e., frost-ring formation co-occurring with reduced ring width) are calculated between 3000 BCE and 1640 CE for 14 major volcanic sulfate signals. The characteristic spacing of sulfate peaks in ice-core-indicated and tree-ring-indicated cooling events has previously been proposed as strong evidence of an age-scale bias in GICC05 over the Late Holocene (Baillie, 2008, 2010; McAneney and Baillie, 2019), including some of the very same tree-ring marker years (e.g., 1627 BCE, 43 BCE, 542 CE) that we now correlate against the WD ice-core record. Before 1 CE, the age differences show only subtle variations within very narrow margins of 3–5 years, with WD2014 ages being a few years too old on average. This indicates that the WD2014 ice-core timescale and, thus, the HolVol v.1.0 eruption database are highly accurate as well as precise for at least the past 5 kyr (and probably also over the full Holocene, given that WD ice-core quality and data resolution improve again) below the brittle zone (i.e., before 4000 BCE) of the WD ice core (Sigl et al., 2016).

4 Discussion

4.1 Deglaciation, volcanism and potential climate feedbacks

The effect of deglaciation on mantle melting beneath Iceland leading to increased eruption frequencies has been recognized since the 1990s (Hardarson and Fitton, 1991; Jull and McKenzie, 1996). Globally, ice sheets reached their maximum extent during the Last Glacial Maximum (LGM), before retreating rapidly during the deglaciation until the Early Holocene (Clark et al., 2012). It is generally thought that the postglacial ice retreat and mass unloading after the LGM resulted in regionally increased frequencies of subaerial eruptions in volcanically active areas, due to increased mantle melting rates (Huybers and Langmuir, 2009). While quantifying the increase in some of these areas (e.g., southern Andes, Cascades and the Kamchatka Peninsula) remains difficult due to the incomplete nature of the geologic eruption record (Watt et al., 2013), the evidence is particularly strong for Iceland (Maclennan et al., 2002; Sinton et al., 2005). Coming out of the glacial, the final deglaciation of Iceland was dominated by rapid ice unloading that peaked between 9800 and 8300 BCE and coincided with an increase in volcanic eruption rates (mass discharge per time) which were 30–50 times higher than the present day. These high eruption rates persisted for over 1000 years after the deglaciation in each area investigated (Maclennan et al., 2002).

When reconstructing a 110 kyr volcanic eruption record from the GISP2 ice core, Zielinski et al. (1996) and Lin et al. (2022) noted a strong increase in the number, magnitude and duration of volcanic sulfate peaks during the termination of the Last Glacial Period and the Early Holocene; they tentatively linked this to crustal responses following the deglaciation. A similarly long-lasting sulfate signal dated to 3160 BCE in HolVol was – in the absence of known volcanic eruptions at the time or confirmative data from other ice cores – attributed by Zielinski et al. (1994) to anomalous marine biogenic emissions. In light of the independent verification of long-lasting acid deposition (see Fig. 5) in 3160 BCE from the GRIP core (Wolff et al., 1997) and new ice-core records that corroborate the fact that volcanic sulfate emissions can be sustained for centuries (McConnell et al., 2017), we interpret the increased frequency and duration of volcanic sulfate deposition in Greenland as additional evidence of the increased volcanic activity, predominantly from Iceland, following the large-scale warming during the deglaciation (Geirsdottir et al., 2009). The spatiotemporal structure of volcanic emissions in HolVol, with 75 % higher SAOD during the Early Holocene (9500–7500 BCE) than between 4000 BCE and 1000 CE and the increased SAOD concentrating in the NH, is consistent with a causal coupling of subaerial volcanism and rapid deglaciation in formerly glaciated volcanic regions.

The increased VSSI rates during the Early Holocene, also from eruptions that we have attributed to the tropics, seem to be at odds with the idea of a coupling between volcanic activity and rapid ice unloading during the deglaciation, as these areas were not glaciated during the glacial period. To investigate this further, we compared the relative distribution of sulfur deposition between Greenland and Antarctica (defined as the asymmetry ratio) for all eruptions in HolVol v.1.0 with those of known tropical eruptions (Fig. S10). We find that the mean asymmetry ratio for attributed tropical eruptions between 9500 and 7000 BCE (N=98) is significantly (p<0.05) different (i.e., indicating a stronger asymmetry of sulfate burden towards Greenland) from the mean asymmetry ratio for attributed tropical eruptions between 7000 BCE and 1900 CE (N=218). This difference arises due to a larger frequency of events with a large (>0.75) asymmetry ratio. We interpret this result to be an indication that the former group (during deglaciation) contains more eruptions that occurred further north than the latter group. The mean asymmetry ratios of both of these groups of bipolar eruptions are significantly (p<0.01) different (i.e., again indicating a stronger asymmetry of the sulfate burden towards Greenland) from the mean asymmetry ratios calculated for the limited number (N=11) of known tropical volcanic eruptions in ice-core records in both HolVol v.1.0 and eVolv2k. Based on this we hypothesize that the apparent increased volcanic activity in the tropics is likely an artifact from our applied and traditionally used source attribution (i.e., bipolar signal equals tropical source). We further hypothesize that some of the increased activity actually took place in glaciated NHET regions (e.g., Iceland, Alaska and the Kamchatka Peninsula), but contemporaneous volcanic sulfate detected in Antarctica had caused us to incorrectly attribute the eruption source to the tropics.

A consequence of this is that we would have overestimated the true frequency of tropical eruptions and underestimated the true frequency of NHET eruptions (some of which plausibly originated from volcanic areas experiencing ice unloading). This hypothesis is gaining increasing support from new tephra identifications for bipolar volcanic events that can be geochemically assigned to eruptions in Iceland (Lin et al., 2022; Svensson et al., 2020) or Alaska (McConnell et al., 2020a; Pearson et al., 2022). Consequently, in recently published studies, “bipolarity” is no longer a definitive criterion for tropical volcanic sources (Abbott et al., 2021a; Lin et al., 2022; Pearson et al., 2022). Linking specific volcanic eruptions with ice-core-indicated sulfate signals in HolVol, however, remains difficult due to the often low age precision of proximal volcanic deposits and the scarcity of tephra to fingerprint and identify eruptions in ice cores during the Holocene (Abbott and Davies, 2012). The 10 ka Grímsvötn tephra series (i.e., the Saksunarvatn Ash, found in numerous Greenland ice cores; see Fig. 5) is one of the few exceptions; however, these Grímsvötn tephra layers are increasingly considered to represent a time interval marker spanning approximately 500 years, rather than a sharp marker horizon, further complicating the alignment of climate proxy and volcanic records in the North Atlantic region (Oladottir et al., 2020). Although difficult to date precisely, large volume (≥2 km³) fissure eruptions and lava shield eruptions from the Western Iceland Zone (e.g., Hallmundarhraun, Leitahraun, Skjaldbreiður and Þingvallahraun) overlap within age uncertainties with some prominent prolonged sulfate signals in the Greenland ice cores (Fig. 5, Table S1). Since Holuhraun, a comparable small fissure eruption producing slightly above 1 km³ of lava in 2014–2015 (Bonny et al., 2018) is readily detectable in snow samples from Greenland; despite increased sulfate background from industrial emissions (Du et al., 2019), the prolonged sulfate signals over the Holocene may plausibly be linked to the aforementioned long-lived eruptions of larger volume with their characteristically low average effusion rates (Sinton et al., 2005).

Although carbon emissions from volcanoes are dwarfed by human emissions (Fischer et al., 2019; Le Quéré et al., 2018), several studies have suggested that post-deglaciation increases in subaerial volcanism evoked potential feedbacks in the climate system primarily through the co-emission of greenhouse gases (e.g., carbon dioxide (Huybers and Langmuir, 2009; Kutterolf et al., 2013). Observations have further shown that individual volcanoes such as Katla in Iceland can act as large point sources of CO₂, emitting 12–24 kt d⁻¹, even during quiescent time periods (Ilyinskaya et al., 2018). Potential gas emissions from prolonged volcanic eruptions lasting for decades are likely several orders of magnitude larger. As estimating the global volcanic CO₂ flux to the atmosphere often involves relating volcanic fluxes of sulfur dioxide with measured or estimated $\mathrm{C}/\mathrm{S}$ molar ratios (Fischer et al., 2019; Werner et al., 2019), our comprehensive HolVol reconstruction of volcanic sulfate now provides a basis for further research to advance our understanding of the coupling between climate and volcanism during the last deglaciation.

4.2 Tropospheric volcanic sulfur emissions

Reconstructions of volcanic aerosol forcing commonly assume that the vast majority of the volcanic sulfate deposited on the polar ice sheets derives from fallout from the stratosphere. Observations of volcanic SO₂ emissions using remote sensing show that a large proportion of the total volcanic sulfur emissions remains within the troposphere. Between 1979 and 2018, a total of 44 Tg SO₂ was emitted globally from effusive eruptions (Carn et al., 2016), of which only 5 % was stratospheric emissions. The mean plume height of these effusive eruptions was <7 km. Over the same time period, 54 Tg SO₂ was emitted from explosive eruptions, of which >80 % was stratospheric emissions with a mean plume height of 16 km (Carn et al., 2016).

For the HolVol reconstruction, as in earlier reconstructions, stratospheric sulfur injection is estimated from the ice-core sulfate fluxes using simple scaling factors, which are assumed to be unbiased in an average sense but uncertain for individual eruptions. Indeed, some sulfur spikes recorded in ice cores may be the result of purely tropospheric emissions. Remote sensing suggests that 1–2 Tg S was emitted up to 1–3 km high during the 6-month-long fissure eruption of Holuhraun starting in September 2014 (Schmidt et al., 2015). An increase in volcanic sulfate and fluoride dated to late 2014 in a northeastern Greenland snow-pit sample indicates that volcanic fallout from this effusive eruption is preserved on the Greenland ice sheet (Du et al., 2019), but the deep ice cores used to estimate volcanic fallout over the Holocene in this study have not been updated to the present. Future work linking recent observed eruptions to ice cores may help to constrain how much of the Icelandic tropospheric emissions from effusive eruptions can be deposited over Greenland and to improve our interpretation of ice-core sulfate records in terms of stratospheric vs. tropospheric content.

While tropospheric sulfur emissions from eruptions lead to aerosols with shorter atmospheric lifetimes, the climate impacts of such emissions may not be negligible, especially for the largest such eruptions. In terms of atmospheric sulfur mass injection, the prolonged fissure eruptions of Lakigar, with 61 Tg S (Thordarson and Self, 2003), and of Eldgjá, with 110 Tg S (Thordarson et al., 2001), quantified using the petrological method, exceeded that of most stratospheric eruptions in the Common Era (Thordarson and Self, 2003; Thordarson et al., 2003). For comparison, these sulfur emissions are twice as high as present-day global annual sulfur emissions from fossil fuel burning and industrial processes (Lamarque et al., 2010). Thus, such extreme eruptions can be seen as natural analogues for the massive tropospheric sulfur emissions that occurred during the 1970s–1980s in industrialized Europe and North America, leading to increased SAOD and reduced solar irradiance at the surface known as “global dimming” (Wild, 2009). In particular, during the Early Holocene, when these types of emissions were more frequent and persistent, as reconstructed from proximal geologic records (Maclennan et al., 2002; Sinton et al., 2005) and from HolVol, our understanding of stratospheric and tropospheric volcanic sulfur emissions remains fragmentary, calling for the development and application of new research approaches and more specific diagnostic ice-core proxies for discrimination.

4.3 Constraining eruption source parameters

Comprising the eruption year and the VSSI, detailed information about the two primary eruption source parameters that define the eruptions' climatic impact are provided in HolVol v.1.0. However, observations and climate modeling suggest that additional eruption source parameters are important (Aubry et al., 2020; Marshall et al., 2019). Specifically, the location, the season of the eruption and the plume height of the SO₂ emission are important, as they have an effect on the specific climate footprint of a given eruption. Refining these secondary – and often interlinked – eruption parameters and quantifying their effect constitute a challenge that is ideally addressed using a multidisciplinary approach with evidence from classical proximal deposits (geologic records), distal fallout (ice cores) and climate models. The detection of volcanic ash (i.e., cryptotephra) preceding the sulfate deposition in Greenland in 43 BCE, for example, allowed for the eruption to be geochemically pinpointed to the caldera-forming Okmok II event in Alaska (McConnell et al., 2020a). The known location not only gave access to key eruption source parameters, such as the magnitude or the volcanic plume height from the proximal deposits (Burgisser, 2005; Crosweller et al., 2012), but also helped to narrow down the eruption date to the winter season, due to the shorter (< weeks) atmospheric lifetime of ash compared with sulfate. The known location (53^∘ N) can, in turn, be employed to evaluate the performance of aerosol–climate models (and statistical emulators) used to project the radiative effects of volcanic eruptions over a wide covarying range of SO₂ emission magnitudes, injection heights and eruption latitudes (e.g., (Marshall et al., 2019). In the case of Okmok II, model simulations constrained by ice-core deposition values from Greenland and Antarctica would imply a most likely eruption latitude between 44^∘ N and 9^∘ S for an eruption in January with a 26.5 km plume height and a 50 Tg S injection (Marshall et al., 2019, 2021). This is only slightly south of the actual location of the eruption at 53^∘ N. Besides information from cryptotephra, information about the potential volcanic sources may also been drawn from the halogen content or the trace element chemistry, as several case studies have demonstrated (Clausen et al., 1997; Kellerhals et al., 2010; McConnell et al., 2017). In addition, sulfur isotopes (³³S and ³⁴S) have been established as a powerful tool (Baroni et al., 2007; Gautier et al., 2019) to differentiate sulfate produced above the ozone layer (i.e., stratospheric) from sulfate forming below the ozone layer (i.e., tropospheric), and they are now applied to ice-core records at a sub-annual time resolution (Burke et al., 2019). Thus, further elucidating eruption parameters using novel methodology, such as targeted cryptotephra analyses and high-resolution S-isotope and trace element measurements, holds great potential with respect to substantially refining and improving HolVol v.1.0 as well as reducing existing uncertainties in the future.

5 Data availability

An archived version of the dataset used in this work (Sigl et al., 2021) is stored on the website of the World Data Center PANGAEA (https://doi.org/10.1594/PANGAEA.928646). As this reconstruction is expected to be updated as new ice-core records become available, or as existing records are revised or reprocessed and new attributions are made, a systematic versioning scheme is proposed to track changes that assigns a unique identifier to each version. The versioning scheme proposed is as follows: the version number for a data compilation is of the form C1.C2, where C1 is a counter associated with the publication of a set of sulfate ice-core records, and C2 is a counter updated every time a modification (latitude, VSSI value, or time) is made to the data or metadata for an individual eruption. Thus, the volcanic forcing published here is v.1.0 of the HolVol dataset. Future versions of the dataset, along with a change log that specifies the modifications associated with each new version, will be posted on PANGAEA. We recommend the use of this dataset for all applications focusing on the entire Holocene. For shorter time periods, we endorse the use of the recommended eVolv2k database (500 BCE–1900 CE) or the “historical” volcanic forcing (1850–present) recommended by CMIP6 archived at https://cera-www.dkrz.de/WDCC/ui/cerasearch/ (last access: 6 July 2022).

6 Conclusion

The Holocene is the latest interglacial period, characterized by rather warm and fairly stable climatic conditions compared with glacial periods. Therefore, it is a critical baseline period for present and future climate change caused by anthropogenic emissions of greenhouse gases. Despite its apparent stability, the Holocene climate shows variability on different scales, including rapid cooling events (Mayewski et al., 2004; Wanner et al., 2011). These changes are induced by internal processes of the coupled climate system and by changes in the external natural forcing, including changes in volcanic, solar, orbital and greenhouse gas forcing. Among these, estimates of volcanic forcing over the Holocene (e.g., Zielinski et al., 1994; Kobashi et al., 2017) have mainly been based on a single-ice-core record from Greenland, which lacks critical information about the timing, magnitude and potential locations of past eruptions, due to limited chemical measurement, temporal resolution and dating accuracy.

Here, we present a reconstruction of volcanic stratospheric sulfur injection (VSSI) from volcanic eruptions extracted from a network of four ice-core sulfate and sulfur records from Greenland and Antarctica that cover the Holocene (i.e., from 11.5 ka or 9500 BCE onwards). With a data coverage of 95 % in Greenland and 100 % in Antarctica, we consider this reconstruction to be virtually complete for all eruptions that injected at least 5 Tg S into the stratosphere, which is about half as much as the eruption of Pinatubo in 1991. The timing of the estimated volcanic eruptions is based on the high-precision WD2014 chronology from the WD ice core in Antarctica, and cross-comparison with absolutely dated tree-ring chronologies throughout the Holocene using frost-ring occurrences and cosmogenic radionuclides (i.e., ¹⁰Be and ¹⁴C) confirms that the absolute dating accuracy is better than ±15 years on average during the Early to Middle Holocene and better than ±1–5 years on average for the past 5000 years. Using the latitudinal mean distribution of volcanic eruptions from other geological records (i.e., GVP) as a guide, we estimate the likely latitudinal position of past volcanic eruptions for which the source volcano is generally unknown or unconfirmed, unless tephra was identified in ice cores and correlated to a known eruption. We find, on average, a distribution of the number of eruptions in the lower latitudes (40 %) and extratropical eruptions in the respective NH (50 %) and SH (10 %) that is very similar to the previously reconstructed structure from a larger network of ice cores over the past 2500 years. The distribution closely resembles the placement of landmasses and the distribution of global subaerial volcanic activity. VSSI estimates from HolVol v.1.0 and from eVolv2k over the period of overlap (500 BCE to 1900 CE) agree well, as should be expected given that three of the four ice cores from HolVol were also included in the eVolv2k database.

The eruption frequency and cumulative sulfur injection was elevated during the Early Holocene (9500–7000 BCE), most notably in the NHET, which we attribute to increased emissions from formerly glaciated volcanic regions such as Iceland. The most notable difference in the character of the Greenland ice-core proxy records is the higher (and reproducible) decadal to multidecadal variations in the volcanic sulfate concentrations, in particular during the Early Holocene. Based on tephra geochemistry available from ice cores that links some of these signals to Icelandic eruptions and based on the known surge in postglacial volcanic activity in Iceland at this time, we interpret these records as evidence of prolonged episodes of volcanic sulfate emissions from Icelandic shield volcanoes, lava floods and fissure eruptions. Dominated by a basaltic composition and effusive character, the plume heights and, ultimately, the climate impact potential of these eruptive episodes are currently poorly constrained, resulting in large uncertainties in the VSSI estimates. Our results further support a strong causal connection between glaciation and volcanic activity that is commonly explained through changes in mantle melting following rapid mass unloading of the retreating glacial ice sheets. No increases in the number of events or the size of volcanic emissions were recorded in the SH where large ice sheets were comparably small.

Besides the mentioned increase in volcanic activity during the Early Holocene, the most notable time periods of increased volcanic activity and emissions were in the Middle Holocene (6700–4300 BCE); in contrast, the 3rd millennium BCE was the most “quiet” period in a Holocene context, comparable to the Medieval Quiet Period or Roman Quiet Period, respectively, but of longer duration. The sulfur injections of the largest known eruptions of the Common Era (Samalas 1257 and Tambora 1815) do not rank among the eight largest eruptions of the Holocene, which were strongly clustered in the early Middle Holocene, in agreement with the age estimates available for the few known eruptions (e.g., Crater Lake, USA; Kikai, Japan; Kurile Lake, the Kamchatka Peninsula) with a volcanic explosivity index (VEI) of 7.

We further used the timing, location and sulfur mass injection to estimate the changes in stratospheric aerosol properties deriving a temporally and spatially resolved continuous reconstruction of SAOD. Thus, the HolVol reconstruction provides the necessary input data for climate model simulations aiming to include volcanic climate forcing in climate model experiments going as far back in time as 9500 BCE. Reconstructed VSSI can be directly incorporated into dedicated aerosol–climate models. As an alternative, the EVA forcing generator (Toohey et al., 2016b) can be used to determine the optical properties of the stratospheric aerosol (i.e., SAOD) on the basis of the VSSI data. For model experiments aiming to perform seamless simulations of climate from 9500 BCE to the present, we recommend using HolVol v.1.0 until 500 BCE (or 1 CE) and the eVolv2k database (the recommended forcing for PMIP4 past2k simulations) until 1900 CE. Between 500 BCE and 1 CE, both reconstructions are based on four individual ice cores as original input data. After 1 CE, eVolv2k is based on up to 16 ice-core records, reconstructed using a very similar methodology, with only subtle differences such as the default latitude used for unknown eruptions (45^∘ N, 0^∘ and 45^∘ S vs. 48^∘ N, 5^∘ N and 37^∘ S in HolVol).

Future work should be targeted at reducing the existing uncertainties, which are largest for time periods with increased volcanic background sulfate that also hampers the correct identification of bipolar (i.e., likely tropical) eruptions. Currently, the attribution of these periods is based solely on the duration of sulfate deposition as the discriminating factor. Moreover, objective geochemical tools are urgently needed for better identification of the source volcanoes, including cryptotephra, halogen content or trace element composition. An equally important goal in the future must be to reduce uncertainty in the transfer functions used to estimate atmospheric sulfate from ice-core sulfate fluxes, in particular for nonexplosive prolonged eruptions similar to those of Laki in 1783–1784 and Holuhraun in 2014–2015. Besides employing present-day observations, remote sensing and aerosol modeling, ice-core records need to be extended in time to the present.