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 15 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 year BP to 1900 CE), constituting an extension of the previous record by 7,000 years. The database incorporates new-generation ice-core aerosol records with sub-annual temporal resolution and 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 20 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 TgS, of which 329 (39%) are located in the low latitudes with bipolar sulfate deposition, 426 (50%) are located in the Northern Hemisphere (NH) extratropics and 88 (10%) are located in the Southern Hemisphere (SH) extratropics. The spatial distribution of reconstructed eruption locations is in agreement with prior reconstructions for the past 2,500 years. In total, these eruptions injected 7410 25 teragram of sulfur (TgS) 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 30 of volcanic eruptions (+16%) and cumulative VSSI (+86%) compared to the past 2,500 years. This increase coincides with the rapid retreat of ice sheets during deglaciation, providing context for potential future increases of volcanic activity in regions a of volcanic stratospheric sulfur injection (VSSI) from volcanic eruptions extracted from a of four ice-core sulfate and sulfur records from Greenland as well as Antarctica covering the Holocene (i.e., from 11.5 9500 With a data coverage of 95% in Greenland and 100% in Antarctica we this reconstruction to be virtually complete for all eruptions that injected at least 5 TgS into the stratosphere, about half as much as the of Pinatubo in 1991. The timing of estimated volcanic eruptions is based on the high-precision 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., 10 Be, 14 C) confirms that the absolute dating accuracy is on average better than years during the Early-Mid Holocene and better than years for the past 5,000 years. Using the latitudinal mean distribution of volcanic eruptions from other geological records (i.e., GVP) as a guide, we 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 number of eruptions in the lower latitudes (40%) and extra-tropical eruptions in the Northern (50%) and Southern (10%) hemisphere, respectively, is very similar to the previously reconstructed structure from a larger network of ice cores over the past 2,500 years. distribution closely resembles the of landmasses and distribution of global subaerial volcanic activity. VSSI estimates from HolVol the Author contributions. conceived this study, performed ice-core analyses, developed age-models and analyzed data; MT 820 performed calculations, analyzed data and led data curation; JRM, MSe and JC-D performed ice-core analyses; MSi led the manuscript writing with input from coauthors.

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 these limitations have so far hampered the identification of stratospheric tropical eruptions over the Holocene. 70 Reconstruction of volcanic forcing requires that all individual ice-core records are synchronized to a common timescale achieved using records of volcanic fallout (e.g., acidity, sulfur, sulfate) archived in the ice-sheets (Parrenin et al., 2012;Seierstad et al., 2014;Severi et al., 2007;. Aligning the records is possible because volcanic aerosols from eruptions is well mixed in the stratosphere and quickly dispersed often on a hemispheric to global scale (Robock, 2000;Toohey et al., 2013). As reference chronology, we use the annual-layer counted 'WD2014' timescale from the WAIS-Divide (WD) 75 ice core in Antarctica providing the highest absolute dating accuracy for 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 abbreviations used in this paper is included (Table S1, Supplement). The exceptional high resolution of the WD sulfur and sulfate records (Cole-Dai et al., 2021) in tandem with the increased dating precision (Sigl et al., 2016) enable us to conduct a firm bi-polar synchronization with ice cores from Greenland over the Holocene as was previously demonstrated for the 80 Common Era (Plummer et al., 2012;Sigl et al., 2013;Sigl et al., 2015). Large volcanic eruptions from the lower latitudes resulting in 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 85 layer, allowing an independent test of the association of sulfate peaks and stratospheric eruptions with global sulfate distribution (Baroni et al., 2008;Baroni et al., 2007;Burke et al., 2019;Gautier et al., 2019). Combining information from Antarctica and Greenland enables us therefore to disentangle likely source regions of volcanic eruptions (i.e., northern hemisphere extratropics (NHET), 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 90 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 likely to occur also 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 13 th centurymarking the inception of the Little Ice Ageand the 1 st century CE 95 (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 21 st and 22 nd centuries with projected changes of the cryosphere under global warming (Schmidt et al., 2013;Tuffen, 2010).
Understanding of how future volcanic activity may affect climate is strongly dependent on understanding the statistical nature 100 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).

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 105 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 -2 yr -1 ) 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) at 75.00°S, 00.07°E, 2882 m a.s.l., has with 68 kg m -2 yr -1 a 2-3 times higher accumulation rate than the one at Dome  (Seierstad et al., 2014). Figure S1 (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 115 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.  (Wolff et al., 1997) 2.2 Ice-core measurements on WD Sulfur concentrations between 1300 m and 2003 m depth (4060-10000 BCE, 6010-11950 year BP) covering the early-to mid-Holocene were analysed 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., 125 2017;McConnell et al., 2018). Depth resolution for sulfur achieved with this system is 1 cm in ice, allowing to achieve nominal monthly time resolution over the entire Holocene. Sulfate concentrations between 577 m and 1300 m depth  covering the mid-late Holocene and the brittle zone (Neff, 2014) of the WD ice core were analysed using ion chromatography in discrete and continuous flow analysis mode (Cole-Dai et al., 2021;Cole-Dai et al., 2013;Cole-Dai et al., 2006) at the South Dakota State University, USA. Depth resolution for sulfate was 2 cm. High sampling resolution throughout the Holocene 130 permitted detection of annual cycles in impurity data, allowing for 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 0.
For key events and time periods we also report ages as years Before Present (BP, years before 1950) a notation used frequently in archaeology, geology and other scientific disciplines. 135

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 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;Svensson et al., 2020). For Greenland and the Arctic, many ice cores (e.g., NGRIP, GRIP, GISP2) had been synchronized 140 Seierstad et al., 2014) on the GICC05 chronology (Rasmussen et al., 2006;Svensson et al., 2008;Vinther et al., 2006), whereas for Antarctica the WD2014 chronology (Buizert et al., 2015;Sigl et al., 2016) serves as the reference chronology (Buizert et al., 2021;Buizert et al., 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 145 signals in Greenland ice-core from high-latitude volcanic eruptions (e.g., Iceland, Alaska, Kamchatka) 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 .
The accuracy of stratigraphic matches further depends on volcanic signal (e.g., temporal resolution) and ice-core site-specific 150 properties (e.g., accumulation rate variability) of both 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 17 th century BCE, 2910 BCE, 45-43 th century BCE, 67-63 th century BCE) were synchronized (major tie-points). Confidence in these match points derived from the combination of (1) a temporally sequence of distinctive signals, (2) comparable magnitudes, (3) a uniform evolution of layer thickness between stratigraphic 155 tie-points, (4) a distinctive shape of the common signals in some cases and finally (5) independent age constraints from 10 Be reflecting 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, Supplement). Relative accumulation rates in Antarctica calculated for longer time periods usually show low variability, narrowing the window for potential stratigraphic 160 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  and dielectric profiling (DEP) from GRIP (Wolff et al., 1997). Confidence in the match points derived from the 165 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 10 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). 170   Sporadic volcanic sulfate deposition at the ice cores sites is superimposed on background deposition of marine sulfate and other sulfuric species (e.g., methanesulfonic acid). This background is seasonally variable, but 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 non-volcanic 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 195 al., 2007;Sigl et al., 2013; 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 for Antarctica, for which three continuous ice-core records are available, a stacking approach similar to the one used for the past 2,000 years  and in accordance with previous work (Gao et al., 2008;Crowely & Unterman 2013). 200

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, 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 -1 ) by 5.1 and 6.7, respectively, with the scale factors determined by matching average Holocene sulfate deposition with sulfur concentrations at WD. The scale 205 factors thus account for differences in 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. The resulting sulfur time series of EDML and EDC can thus be interpreted as the equivalent sulfur concentration at the WD site allowing the construction of an annual-resolved sulfur concentration stack by averaging the three ice-core records (ANT12k, N=3, Fig. 3). We use this timeseries to assess the plausibility of the age synchronization based on the relative agreement of peak amplitudes as an additional diagnostic criterion 210 (see Figure S3, Supplement). The choice of alternative methods such as standardization or normalization) has no significant influence on the results of this reconstruction (see Figure S4, Supplement). Since 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 non-215 volcanic 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 comparable to previous work in Antarctica Gao et al., 2008). A year was classified to contain volcanic sulfur if the annual sulfur concentration exceeded this threshold. After 220 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. Annual volcanic sulfate concentration is calculated as the difference between the total sulfur concentrations of that year and the RRM of the non-volcanic sulfate of that year. The cumulative sulfate mass deposition (kg km -2 ) by an eruption often referred to as (cumulative) "volcanic sulfate flux" f(volc-SO4 2-) is the 225 sum of annual volcanic sulfate concentrations in the years when volcanic deposition occurred multiplied by the mean annual accumulation rate at WD (210 kg m -2 yr -1 ). 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) using the relation f(volc-SO4 2-)ANT12k =1.769 x f(volc-SO4 2-)AVS2k (R 2 =0.93, N=105) to estimate the ice sheet average sulfate fluxes for Antarctica, henceforth f(volc-SO4 2-)AVS12k. As start date of the volcanic eruption 230 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.

Volcanic sulfate deposition in Greenland
For Greenland, we followed a similar approach with some adjustments owing to the different properties of the available 250 volcanic proxy records. With GISP2, only a single ice core with continuous sulfate measurements exists with 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, NEEM-2011-S1) and electrical records (e.g., DEP from GRIP), which we attributed to long-lasting volcanic episodes from Iceland (Fig. 5). We, therefore, tagged all GISP2 volcanic sulfate values exceeding for a minimum of 10 consecutive 255 years the volcano detection threshold as 'prolonged eruption' (Table 3) and applied an additional correction to estimate sulfur injection (see section 2.5). The non-volcanic background sulfate concentration was initially approximated in GISP2 with a 121-points (window) RM fit to the biannual sulfate data. This is equivalent to a 240-year median and thus better suited to detect decadal-to-multidecadal volcanic sulfate variability than with shorter window lengths. Similar to the approach used for Antarctica, we detected volcanic 270 events that exceeded a threshold of RM+1.5×MAD and samples were deemed to contain volcanic fallout if the sulfate concentration exceeds this threshold. After 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. Annual volcanic sulfate concentration is calculated as the difference between the total sulfate concentration of that sample and the RRM of the non-volcanic sulfate 275 of that sample. Finally, the cumulative sulfate mass deposition flux is the sum of volcanic sulfate concentrations in the years when volcanic deposition occurred multiplied by the mean annual accumulation rate at GISP (210 kg m -2 yr -1 ). We tested the performance of our detection criteria during the pre-industrial 19 th century and found that volcanic sulfate is detected during the same periods for which volcanic eruptions had been previously detected by other higher-resolved sulfate records from Greenland ( Fig. S5, Supplement). No false positive events are reconstructed from GISP2 in this test, and volcanic signals 280 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) with 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. a: tephra in ice cores from multiple ice cores from Greenland indicates the Icelandic eruption of Grimsvötn (Saksunarvatn Ash) as a potential 290 source contributing to the ice-core sulfate (Gronvold et al., 1995) ; b: tephra in an ice core from TUNU2013 Greenland indicates Katla (Iceland) as a potential source contributing to the ice-core sulfate Plunkett et al., 2020). See Table S1 (Supplement) for a list of lava shield and fissure eruptions >2 km 3 from Iceland following Hjartarson, (2003) and Sinton et al., (2005).

Ice Core Uncertainties
The timing of volcanic eruptions from ice-core records is uncertain due to interpretation uncertainties during the construction 295 of the annual-layer dating. Based on the comparison of WD2014 with accurately-dated tree-ring records Sigl et al., 2016) we estimate that absolute age uncertainties in the ice-core records used in HolVol are better than ±1 to 5 years on average over the last 2,500 years and better than ±5 to 15 years for the rest of the Holocene. A 5,000 year-long tree-ring record with strong sensitivity for abrupt post-volcanic cooling (Salzer and Hughes, 2007) allows to further assess the absolute age accuracy of WD2014 following some of the largest late Holocene eruptions (see section 2.8). Another source of uncertainty 300 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% for estimating Greenland ice-sheet-wide average flux from mean sulfate flux from the single GISP2 record (Toohey and Sigl of the AVS12k composite stack including three ice cores. The estimated total error of the mean for Antarctica are thus slightly 305 above 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).

Injection locations and dates
Over the last 2,500 years, the localities and the timing of several (N=31, Toohey and Sigl 2017) stratospheric sulfur injections 310 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 1477 CE (Abbott et al., 2021a), Samalas 1257 CE (Lavigne et al., 2013), Changbaishan 946 CE Sun et al., 2014), Eldgjá 315 939-940 CE (Oppenheimer et al., 2018;Zielinski, 1995), Mt. Churchill 853 CE (Jensen et al., 2014), Katla 822 CE Plunkett et al., 2020) and Ilopango 431 CE . 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). Only a handful of ice-core sulfate peaks in the Holocene have to date been linked geochemically 320 to known eruptions, including the caldera-forming 43 BCE Okmok II eruption in Alaska (McConnell et al., 2020a), the c. 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 in Kamchatka (Cook et al., 2018) and the c. 10 ka Grímsvotn "Saksunarvatn Ash" eruption series from Iceland (Oladottir et al., 2020). The vast majority of large eruptions such as the caldera-forming Bronze Age 325 Thera/Santorini eruption, or the c. 6440 ±25 BCE caldera-forming Kurile Lake eruption in Kamchatka remained so far 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 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 330 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 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 eruptions in the extratropics of the Northern and Southern Hemisphere, and the tropics, respectively. We further attribute all volcanic events for which volcanic sulfate deposition to Greenland persisted for more than 10 years to 335 Icelandic source eruptions most likely from the Katla, Bárðarbunga and Grímsvotn volcanic systems or from shield volcanoes in the Western Volcanic Zone (Hjartarson, 2003;Sinton et al., 2005; and assign 64°N as the default latitude for these prolonged episodes (Fig. 5, Table 3; Table S1, Supplement). We are aware that the use of default latitudes paints a crude picture of the geographical distribution of past global volcanic activity (see Figure 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 340 of aerosol optical properties and radiative forcing using the EVA forcing generator (see section 2.7), only the broad region of the eruption site is important (i.e., tropics, NHET, SHET), 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 345 Sigl (2017), an eruption date of 1 January is assigned to unidentified eruptions.

Stratospheric sulfate injection estimation
Stratospheric sulfate injections are estimated from the ice sheets sulfate flux composites using a method described in detail by Toohey and Sigl (2017). Briefly, ice sheet average sulfate fluxes for Antarctica and Greenland are related to injected sulfur mass following Eq. (1): 350 where and are transfer functions accounting for the spatial distribution of sulfate deposition over each hemisphere. Based on analysis of the spread and deposition of nuclides from nuclear bomb test, sulfate from prior volcanic eruptions and atmospheric model simulations (Gao et al., 2007), the transfer functions and are estimated to be 1 x10 9 km 2 for tropical eruptions and 0.57 x 10 9 km 2 for extratropical eruptions. The method described above is based on the assumption that the ice-355 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 (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 with both stratospheric and tropospheric transport pathways deposited in Greenland following the large VEI=6 eruption of Katmai/Novarupta in Alaska, upwind of the Greenland ice sheet (Burke et 360 al., 2019). Of particular importance are long-lasting, effusive (i.e. non-explosive) 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-40 and Laki 1783-84 AD) are the most prominent examples and the extent to which sulfate was injected during these eruptions into the lower stratosphere is 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 available of 365 key eruption source parameters (e.g., plume height, SO2 dispersion height, duration, season), which determine how much sulfur gets injected into the stratosphere. Detailed volcanological fieldwork could delineate 10-11 distinctive eruptive episodes during the Laki 1783-84 event , allowing the development of detailed SO2 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 370 (Thordarson and Larsen, 2007;. 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-84 and derived a ratio of stratospheric to total sulfate deposition of 0.15. Due to a lack of data on the stratospheric versus tropospheric distribution of injected sulfur for other major Icelandic eruptions of the Holocene we adopt this approach and used a transfer function of 0.10 x 10 9 km 2 for those long-lasting deposition signals we assume to be from 375 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 8 subsequent injections), and durations 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. 380 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 which aim to estimate uncertainties from the first two terms. Uncertainty related to the limited number of ice cores and related sampling of the ice sheets has been 385 estimated (see section 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% uncertainty to the NH and SH transfer function (LNH and LSH), respectively, but these estimates may be model-dependent and recent work points to potentially larger values (Marshall et al., 2019, Marshall et al., 2018. It remains difficult to quantify errors arising from a potentially wrong 390 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 Northern hemisphere extra-tropics may not be recognized as separate eruptions in the biannual-resolution GISP2 record and thus erroneously attributed as a prolonged eruptive period when sulfate levels remain increased for >10 395 years. Under such a scenario, the true VSSI would be underestimated by up to 80%. To account for this later case, reported VSSI uncertainties for prolonged eruptions have been inflated to include magnitudes which 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 to explosive extratropical eruptions. We note that only specific eruptions may be subject to errors caused by wrong attributions which can be subsequently assessed and 400 corrected in future updates of this database if independent constraints for source locations from crypto-tephra, 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;McConnell et al., 2020a). A primary source of systematic error in the VSSI estimates is likely to originate from the uncertainty in the transfer functions and used to estimate hemispheric stratospheric sulfate aerosol burdens. These are originally derived using ice-sheet sulfate fluxes in Antarctica and observed from stratospheric sulfate burden from the 405 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).

Aerosol optical depth estimation 410
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 of 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 spatio-temporal structure of the EVA output fields is based on a simple three-box model of stratospheric transport which is optimized to produce agreement with observations of the aerosol cloud 415 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 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 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 420 removal time has been implemented in EVA based on ECHAM5-HAM aerosol model simulations of eruptions with a wide range of magnitude (Metzner et al., 2014). Based on the model results, stratospheric aerosol removal timescale is varied between its nominal value of 11 months and a minimum of 6 months as global stratospheric sulfate burden rises above 10 TgS.
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 425 HolVol used as an input to EVA. SAOD time series are shown as either monthly, annual or centennial averages. Peak SAOD values can therefore differ significantly depending on the time resolution of the time series.

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 10 Be in ice 430 cores  and 14 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 BP) is estimated to be better than ±20 years. During the early Holocene at 9500 BCE (11.5 ka BP), the WD2014 age uncertainty was estimated at ±66 years.
Matching the common production signal in cosmogenic isotopes ( 14 C, 10 Be) has further allowed us to assess the WD2014 ages relative to the radiocarbon calibration curve, which is during the Holocene based on dendrochronology and thus has virtually 435 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 that WD2014 ages are slightly too young during the early Holocene and are slightly too old between 7000 BCE to 1 CE (Sigl et al., 2016). No 10 Be measurements from WD were previously available and thus no assessment of the WD2014 timescale was possible for the time period between 3500 BCE and 500 BCE. To fill 440 this gap, we employ a multi-millennial compilation of the occurrence of ring-width minima and frost-rings in a bristlecone pine chronology from SW USA covering the past 5,000 years (Salzer and Hughes, 2007). A strong association between frostring 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 445 from the compilation by Salzer and Hughes (2007) all marker events (N=10) in which at least two consecutive ring-width minima corresponded with a frost damaged ring within an error margin of ±1 year (Table 4). Additional marker years in which a frost-ring corresponds with a ring-width minima within ±1 year are provided in Supplement (Table S2).   a: data gap in GISP2; VSSI is only based on Antarctica assuming a SHET source eruption; VSSI may be underestimated if a comparable large sulfate anomaly is detected in Greenland ice core records. b: period with long lasting reductions of ring-width and frequent frost-ring appearance following a large tropical eruption in 682 CE (Table  S2, Supplement). 460

Volcanic stratospheric sulfur injections
The HolVol v.1.0 VSSI time series is shown in Figure 6 and Figure 7. With 100% data coverage in Antarctica and 95% data coverage from Greenland we consider this record to be virtually complete for all volcanic eruptions with a strong climate impact potential (i.e. VSSI comparable or larger to the 1991 Pinatubo eruption). The ability to detect and quantify smaller events is primarily limited by data gaps (equivalent to 560 years) and coarse (bi-annual) temporal resolution of the GISP2 465 record from Greenland. Therefore, there is a potential of under-recording smaller eruptions situated in the NHET. With only 88% data coverage between 3000-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. into the stratosphere, of which 75% was emitted in the tropics, 4% in Iceland, 18% in the NH extra-tropics and 4% in the SH, respectively (Fig. 8). figure   490 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).

Figure 8: Number and cumulative volcanic stratospheric sulphur 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.
Number of eruptions and cumulative centennial VSSI varied within the Holocene (Fig. 6, Fig. 7). In general, the number of eruptions and the cumulative VSSI were enhanced during the early-to-mid-Holocene (10 th to 5 th millennium BCE, on average 500 76 TgS per century). Between the 4 th millennium BCE to the present, both the average number of eruptions and the cumulative VSSI (45 TgS per century) were lower by 21% and 41%, respectively (Fig. 6, Fig. S4, Fig. S7, Supplement). The period 9500 BCE to 7000 BCE, when glacial ice-sheets were retreating rapidly and 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 TgS injected in the stratosphere per century this period which we term 'Deglaciation Active Period' is 43% above of the 505 Holocene mean VSSI rate of 63 TgS. An increased VSSI rate is noted for eruptions in the NHET as well as in the tropics, but is absent for the SHET (see Figure S6, Supplement). The majority of the events we attributed to prolonged eruptive episodes (379 years of cumulative duration) also falls into this time period (see Fig. 6, Table 3).
The window 4000 BCE to 1000 BCE has the lowest frequency of eruptions and smallest VSSI rates in the Holocene. With 21% less eruptions and 36% smaller VSSI rates than the Holocene mean we term this period as the 'Holocene Quiet Period' 510 in analogy to other time periods with reduced volcanic activity such the Medieval Quiet Period (700-1100 CE) and the Roman Quiet Period . Throughout the Holocene the longest subsequent time period without an eruption with VSSI >1 TgS is 77 years (ending in 3206 BCE), that without VSSI >10 Tg is 317 years (ending in 2258 BCE). These volcanically quiet periods are slightly longer in comparison to the same metrics for the Medieval Quiet Period (55 and 217 years) and the Roman Quiet Period (71 and 212 years), respectively. 515 Eight of the ten largest VSSI injections are recorded between 6700 BCE 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 Mid-Holocene Active Period.

Comparison with other Holocene volcanic reconstructions 520
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 BCE and 50 BCE (Hammer et al., 1980). A direct comparison on an event-base is not possible owing to the different chronology compared to this study.
Age uncertainties in the Camp Century record of ±170 years 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 recently 525 has 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 all acids were from H2SO4) from Camp Century was 40 TgS, slightly below the estimate of 56 TgS in HolVol. The highest global volcanic fallout estimates in Camp Century were 85 TgS using latitudinal correction functions that were assuming a high-latitude eruption source for the vast majority of the ice-core signals. These are significantly smaller compared to the highest estimates in HolVol of up to 190 TgS, which were 530 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 ka), 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 535 ±150 years some 8 ka before present (Meese et al., 1997). Zielinski et al., (1994) did not estimate sulfur injection, or 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 are under-recording eruptions from the SHET and prone to systematically overestimate the forcing 540 from Icelandic eruptions and many other NHET eruptions by at least 43% and up to a factor of 10 for specific events.

Comparison with eVolv2k
The eVolv2k volcanic eruption catalogue (500 BCE -1900 CE) was reconstructed from bipolar ice-core records using a similar methodology as 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 545 for transient climate model simulations covering the past millennium or the past 2 ka including experiments (Jungclaus et al., 2017) within the Paleoclimate Model 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 explaining 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 TgS from 180 eruptions with >1 Tg S injection between 500 BCE and 1900 CE which is only 10% above the value as estimated 550 based on eVolv2k (Fig. S8, 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., above 10 TgS) while there is a larger scatter for eruptions with smaller VSSI (Fig. S8, Supplement). In order to perform a seamless Holocene-long simulation with climate models we recommend to merge the VSSI or SAOD reconstructions from HolVol v.1.0 with those from eVolv2k (see Fig. S9, Supplement) at the year 500 BCE or 1 CE. 555

Stratospheric aerosol optical depth and radiative forcing
Global mean SAOD from the EVA(HolVol) reconstruction is shown in Fig. 9. SAOD follows closely the spatio-temporal structure of VSSI in the Holocene, albeit with relatively less pronounced peaks for the largest eruptions due to the nonlinear parameterizations used in EVA. Global mean SAOD over the Holocene was 0.0153, SAOD over the Northern hemisphere was with 0.0182 almost 50% higher than over the Southern hemisphere (0.0124). Global mean SAOD between 9500 BCE and 560 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 Northern hemisphere (0-90°N) and less pronounced (+21%) when integrating over SHET (30-90°S). The largest annual global 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 global annual SAOD of 0.50; the largest non-tropical eruption of the Common Era in 565 536 CE produced NHET SAOD of 0.43. We stress that for such large eruptions, significantly larger than any eruption observed in the instrumental era, uncertainties in the SAOD should be understood to be large. The EVA(HolVol) reconstruction is compared to that of Kobashi et al., (2017) in Fig. 10. Since the Kobashi et al (2017) reconstruction contains estimates of radiative forcing (RF, in units of Wm -2 ), the EVA(HolVol) SAOD values are converted to RF using the linear scaling (RF=-25*SAOD) of Hansen et al., (2005) used in prior IPCC reports (Myhre et al., 2013). We 575 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 of the multi ice-core HolVol reconstruction and the single ice-core (GISP2) reconstruction from Greenland are the smaller magnitudes (minima of -21 W m -2 ) of RF for large volcanic eruptions in HolVol compared to values as strong as -50 W m -2 as reconstructed by Kobashi et al., (2017) which we attribute to applying a nonlinear scaling to HolVol. Furthermore, HolVol RF 580 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, therefore, 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 585 Greenland are underestimated. Some of the difference is also due to the additional inclusion of a non-zero 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).

Dating accuracy and precision
for details). The previous assessment based on correlating multi-decadal to centennial-scale production rates in the cosmogenic radionuclides 10 Be and 14 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 (through overcounting of 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 600 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 cores and tree-ring indicated cooling events, had previously been proposed as strong evidence for an age-scale bias in GICC05 over the late Holocene (Baillie, 2008(Baillie, , 2010McAneney and Baillie, 2019), including some of the very same tree-ring marker years (e.g. 1627 BCE, 43 BCE, 542 CE) we now correlate against the WD ice-core record. Before 1 CE, the age differences show 605 only subtle variations within very narrow margins of 3 to 5 years, with WD2014 ages being on average a few years too old.
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 ka, and probably also over the full Holocene given that WD ice-core quality and data resolution improving again, below the brittle zone (i.e. before 4000 BCE) of the WD ice core (Sigl et al., 2016).

Deglaciation, volcanism and potential climate feedbacks 620
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 regionallyincreased frequencies of subaerial eruptions in volcanically active areas, due to increased mantle melting rates (Huybers and 625 Langmuir, 2009). While quantifying the increase in some of these areas (e.g., Southern Andes, Cascades, Kamchatka) 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 peaking between 9800 and 8300 BCE coinciding with an increase of volcanic eruption rates (mass discharge per time) which were 30-50 times higher than at present day. These high eruption rates persisted for over 1000 years 630 after the deglaciation in each area investigated (Maclennan et al., 2002). When reconstructing a 110 ka 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 which they tentatively linked to crustal responses following the deglaciation. A similarly long-lasting sulfate signal dated to 3160 BCE in HolVol wasin the absence of known volcanic eruptions at the time or confirmative data from other ice coresattributed by 635 Zielinski et al. (1994) to anomalous marine biogenic emissions. In light of 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 that volcanic sulfate emissions can be sustained for centuries (McConnell et al., 2017) we here 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 spatio-640 temporal structure of volcanic emissions in HolVol, with 75% higher SAOD during the early Holocene (9500-7500 BCE) than between 4000 BCE to 1000 CE and these 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 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, since these areas were not glaciated during the glacial period. To 645 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 (Supplement Figure S10). We find that the mean asymmetry ratio for attributed tropical eruptions between 9500-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 -1900 CE (N=218). This difference arises due to a larger frequency of events 650 with a large (>0.75) asymmetry ratio. We interpret this result as an indication that the former group (during deglaciation) contains more eruptions that occurred further north than the latter group. The mean asymmetry ratio of both these groups of bipolar eruptions are significantly (p<0.01) different (i.e. again indicating a stronger asymmetry of 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 655 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, Kamchatka), but contemporaneous volcanic sulfate detected in Antarctica had caused us to attribute incorrectly the eruption source to the tropics. A consequence of that is we would have overestimated the true frequency for tropical eruptions and underestimated the true frequency of NHET eruptions (some of which plausibly originated from volcanic areas 660 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., 2020;Pearson et al., 2022). Consequently, in recently published studies "bipolarity" is no longer a definitive criterion for tropical volcanic sources (Abbott et al., 2021b;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 665 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; but 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 alignment of climate proxy and volcanic records in the North-Atlantic region (Oladottir et al., 2020). Though difficult to date precisely, large volume (≥2 670 km 3 ) fissure eruptions and lava shield eruptions from the Western Iceland Zone (e.g., Hallmundarhraun, Leitahraun, Skjaldbreiður, Þingvallahraun) overlap within age uncertainties with some prominent prolonged sulfate signals in the Greenland ice cores ( Fig. 5; Table S1, Supplement). Since Holuhraun, a comparable small fissure eruption producing slightly above 1 km 3 lava in 2014-15 (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 675 be linked to the aforementioned long-lived eruptions of larger volume with their characteristically low average effusion rates (Sinton et al., 2005).
Whereas carbon emissions from volcanoes are dwarfed by human emissions Le Quere et al., 2018) several studies have suggested that post-deglaciation increases of 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., 680 2013). Observations further showed that individual volcanoes such as Katla in Iceland can act as large point sources of CO2 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. Since estimating of the global volcanic CO2 flux to the atmosphere often involves relating volcanic fluxes of sulfur dioxide with measured or estimated C/S molar ratios 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.

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. Observation of volcanic SO2 emissions using remote sensing show 690 that a large proportion of total volcanic sulfur emissions remains within the troposphere. Between 1979 and 2018, a total of 44 Tg SO2 were emitted globally from effusive eruptions , of which only 5% were stratospheric emissions. The mean plume height of these effusive eruptions was < 7 km. Over the same time period 54 Tg SO2 were emitted from explosive eruptions, of which >80% were stratospheric emissions with a mean plume height of 16 km .
For the HolVol reconstruction, as in earlier reconstructions, stratospheric sulfur injection is estimated from the ice core sulfate 695 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 TgS was emitted up to 1-3 km high during the 6-month long fissure eruption of Holuhraun starting in September 2014 . An increase of volcanic sulfate and fluoride dated to late 2014 in a NE Greenland snow-pit sample indicates that volcanic fallout from this effusive eruption is preserved on the Greenland ice-sheet (Du et al.,700 2019), but the deep ice-cores used in this study to estimate volcanic fallout over the Holocene have not been updated to the present. Future work linking recent observed eruptions to ice cores may help to constrain how much of Icelandic tropospheric emissions from effusive eruptions can get deposited over Greenland and improve our interpretation of ice cores 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 705 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 TgS  and of Eldgjá with 110 TgS (Thordarson et al., 2001) quantified using the petrological method exceeded that of most stratospheric eruptions in the Common Era . For comparison, these sulfur emissions are twice as much as present day global annual sulfur emissions from fossil fuel burning and industrial processes (Lamarque et al., 2010). Such extreme 710 eruptions can thus be seen as natural analogues for the massive tropospheric sulfur emissions that occurred during the 1970-80s 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 715 development and application of new research approaches and more specific diagnostic ice-core proxies for discrimination.

Constraining eruption source parameters
With the eruption year and the VSSI, detailed information about the two primary eruptions source parameters that define the eruptions' climatic impact are provided in HolVol v.1.0. However, observations and climate modelling suggest that additional eruption source parameters are important (Aubry et al., 2020;Marshall et al., 2019). Specifically, the location, the season of 720 the eruption and the plume height of the SO2 emission are important since they have an effect on the specific climate footprint of a given eruption. Refining these secondaryand often interlinkederuption parameters and quantifying their effect constitute a challenge that is ideally addressed in a multidisciplinary approach using evidence from classical proximal deposits (geologic records), distal fallout (ice cores) and from climate models. The detection of volcanic ash (i.e. cryptotephra) preceding the sulfate deposition in Greenland in 43 BCE, for example, allowed to geochemically pinpoint the eruption to the 725 caldera-forming Okmok II in Alaska eruption (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 to sulfate. The known location (53°N) can in turn be used to evaluate the performance of climate-aerosol models (and statistical emulators) used to project radiative effects of volcanic eruptions over a wide 730 covarying range of SO2 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 26.5 km plume height and 50 TgS (Marshall et al., 2019;Marshall et al., 2021). This is only slightly south to the real location of the eruption at 53°N. Besides from cryptotephra, information about the potential volcanic sources may also been drawn on the basis of halogen content or the 735 trace element chemistry as several case studies have demonstrated Kellerhals et al., 2010;McConnell et al., 2017). In addition, sulfur isotopes ( 33 S, 34 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 are now applied to ice-core records on a sub-annual time resolution (Burke et al., 2019). Further elucidating eruption parameters by using novel methodology such as targeted crypto-tephra analyses and high-resolution S-740 isotope and trace element measurements thus holds great potential to substantially refine and improve HolVol v.1.0 and reduce existing uncertainties in the future.

Conclusion
The Holocene is the latest interglacial period, characterized by rather warm and fairly stable climate conditions compared to glacial periods. It is therefore a critical baseline period for present and future climate change caused by anthropogenic 745 emissions of greenhouse gases. Despite the apparent stability, 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 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) were mainly based on a single ice-core record from Greenland, lacking critical information about the timing, magnitude and 750 potential locations of past eruptions, due to limited chemical measurement, temporal resolution and dating accuracy.
We present here 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 as well as Antarctica covering the Holocene (i.e., from 11.5 ka BP 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 TgS into the stratosphere, about half as much as 755 the eruption of Pinatubo in 1991. The timing of 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., 10 Be, 14 C) confirms that the absolute dating accuracy is on average better than ±15 years during the Early-Mid Holocene and better than ±1-5 years for the past 5,000 years. Using the latitudinal mean distribution of volcanic eruptions from other geological records (i.e., GVP) as a guide, we 760 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 number of eruptions in the lower latitudes (40%) and extra-tropical eruptions in the Northern (50%) and Southern (10%) hemisphere, respectively, that is very similar to the previously reconstructed structure from a larger network of ice cores over the past 2,500 years. The distribution closely resembles the situation of landmasses and distribution of global subaerial volcanic 765 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 out of the four ice cores from HolVol were also included in the eVolv2k database.
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 are the higher (and reproducible) decadal to multi-decadal 770 variations of volcanic sulfate concentrations, in particular during the early Holocene. Based on tephra geochemistry available from ice cores linking some of these signals to Icelandic eruptions, and based on the known surge of post-glacial volcanic activity in Iceland at this time, we interpret these records as evidence for 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 for these eruptive episodes are currently poorly constrained 775 resulting in large uncertainties in the VSSI estimates. Our results give further support to a strong causal connection between glaciation and volcanic activity 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 size of volcanic emissions were recorded in the Southern Hemisphere where large ice sheets were comparably small.
Besides the mentioned increase of volcanic activity during the early Holocene, the most notable time periods of increased 780 volcanic activity and emissions were in the mid-Holocene (6700-4300 BCE), whereas the 3rd millennium BCE in contrast was the most "quiet" in a Holocene context comparable to the Medieval or Roman Quiet periods, 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 8 largest eruptions of the Holocene which were strongly clustered in the early mid-Holocene in agreement with the age estimates available for the few known eruptions (e.g., Crater Lake, USA; Kikai, Japan, Kurile Lake, Kamchatka) 785 with a Volcanic Explosivity Index (VEI) of 7.
We further used the timing, location and sulfur mass injection to estimate the changes of stratospheric aerosol properties deriving a time-and space resolved continuous reconstruction of SAOD. The HolVol reconstruction thus 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. 790 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) to 1900 CE. Between 500 BCE and 1 CE both reconstructions are based on 4 individual ice cores as original input data. After 1 CE, eVolv2k is based on up to 16 ice-core 795 records, reconstructed using a very similar methodology, with only subtle differences such as the default latitude used for unknown eruptions (45°N/0°/45°S versus 48°N/5°N/37°S in HolVol).
Future work should be targeted to reduce the existing uncertainties which are largest for time periods with increased volcanic background sulfate which also hamper the correct identification of bipolar (i.e. likely tropical) eruptions. Currently the attribution of these periods is based solely on time-duration of sulfate deposition as the discriminating factor. More objective, 800 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 non-explosive prolonged eruptions similar to those of Laki in 1783-84 and Holurauhn 2014-15. Besides employing present-day observations, remote sensing and aerosol modeling, ice-core records need to be extended in time to the present. 805

Data availability and data versioning
An archival version of the dataset  is stored on the website of the World Data Center PANGAEA (https://doi.pangaea.de/10.1594/PANGAEA.928646). Since 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 810 versioning scheme is proposed to track changes, assigning 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, C2 is a counter updated every time a modification (latitude, VSSI value, time) is made to the data or metadata for an individual eruption. The volcanic forcing published here is thus 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 815 version, will be posted at PANGAEA. We recommend to use this dataset for all applications focusing on the entire Holocene.
Author contributions. MSi conceived this study, performed ice-core analyses, developed age-models and analyzed data; MT 820 performed calculations, analyzed data and led data curation; JRM, MSe and JC-D performed ice-core analyses; MSi led the manuscript writing with input from all coauthors.
Competing interests. The authors declare that they have no conflict of interest. which in turn received support from the Swiss Academy of Sciences and the Chinese Academy of Sciences. We thank the students and staff at South Dakota State University, the Desert Research Institute and the University of Firenze for contributing to the ice-core chemical analysis of the WD and EPICA ice cores. M.Si. also thanks Eric Wolff for sharing of ice-core data.