Magnitude, frequency and climate forcing of global volcanism during the last glacial period as seen in Greenland and Antarctic ice cores (60-9 ka)

. Large volcanic eruptions occurring in the last glacial period can be detected by their accompanying sulfuric acid deposition in continuous ice cores. Here we employ continuous sulfate and sulfur records from three Green-land and three Antarctic ice cores to estimate the emission strength, the frequency and the climatic forcing of large volcanic eruptions that occurred during the second half of the last glacial period and the early Holocene, 60–9 kyr before 2000 CE (b2k). Over most of the investigated interval the ice cores

in the last 2500 years and 69 eruptions are estimated to have larger sulfur emission strengths than the VEI-7 30 Tambora eruption that occurred in Indonesia in 1815 AD. The frequency of eruptions larger than the typical VEI-7 (VEI-8) eruption by the comparison of sulfur emission strength is found to be 5.3 (7) times higher than estimated from geological evidence. Throughout the investigated period, the frequency of volcanic eruptions is rather constant and comparable to that of recent times. During the deglacial period (16-9 ka b2k), however, there is a notable increase in the frequency of volcanic events recorded in Greenland and an obvious increase in the 35 fraction of very large eruptions. For Antarctica, the deglacial period cannot be distinguished from other periods.
These volcanoes documented in ice cores provide atmospheric sulfate burden and climate forcing for further research on climate impact and understanding the mechanism of the Earth system.

Introduction
The dispersal of gas, aerosols and ash particles by volcanic eruptions play a major role in the climate system 40 (Gao et al., 2007;Robock, 2000). Large volcanic eruptions injecting sulfuric gases into the stratosphere and forming sulfate aerosols have a global or hemispheric cooling effect of several degrees lasting for several years after the eruption (Sigl et al., 2015).
Estimations of volcanic stratospheric sulfur injections and of the timing and frequency of large volcanic eruptions is essential for the ability to understand and model past and future global climate conditions 45 (Timmreck et al., 2016). For the last 1200 to 2500 years, the ice-core based volcanic forcing records derived from Greenland and Antarctica (Crowley and Unterman, 2013;Gao et al., 2008;Toohey and Sigl, 2017) provide an essential record for climate model simulations (Jungclaus et al., 2017) supporting detection and attribution studies (Schurer et al., 2014) including those applied in the IPCC, but so far the global ice-core based volcanic record of the last glacial period is poorly documented.

Ice-core records of volcanic sulfate deposition
Several studies have reconstructed the volcanic sulfate deposition for part or all of the Holocene in Greenland (Cole-Dai et al., 2009;Gao et al., 2008;Sigl et al., 2013) or in Antarctica (Kurbatov et al., 2006;Castellano et al., 2004;Plummer et al., 2012;Nardin et al., 2020;Cole-Dai et al., 2021). Sigl et al. (2015) applied accurately dated ice cores synchronized between the two hemispheres to reconstruct global volcanism over the last 2500 eruptions, 5 of which had a sulfur emission strength large or similar to the Tambora eruption occurring in Indonesia in 1815 CE. Prior to the Holocene no bipolar volcanic sulfate deposition record is currently available from ice cores.

60
One conclusion drawn from historical eruptions is that there is a significant variability of the same volcanic event in the sulfate deposition as determined in different ice cores both on a regional and a local scale (Sigl et al., 2014;Gao et al., 2007). Part of this regional variability can be explained by the difference in sulfate deposition fluxes at different locations. For example, in Antarctica where geographical distances are large, the sulfate deposition at a specific site will be strongly dependent on the location of the eruption, governing wind sulfur deposition observed for recent eruptions, one would expect even larger variability for eruptions occurring during the last glacial period where snow accumulation was generally lower and the thinning of the annual 90 layers in the ice cores becomes more important.

Studies of the frequency of volcanic eruptions
The volcanic sulfate record of the Greenland GISP2 ice core has been investigated by Zielinski et al. (1997), who found that there has been increased volcanic activity during the deglacial period (22-8 ka b2k) as compared to the average activity of the last glacial cycle. This is interpreted as being related to the tectonic isostatic 95 response to the melting of the large ice sheets during that period. Based on the geological (non-ice-core) record, Huybers and Langmuir (2009) found that volcanism increased two to six times during the deglacial period, 12-7 ka b2k, as compared to the average level of eruptions during 40-0 ka b2k interval.
In Antarctica, Castellano et al. (2004) determined the frequency of volcanic eruptions over the last 45 ka based on the EDC ice core. They found a rather constant level of volcanic activity throughout that period except for 100 the most recent millennia, where the activity shows an increase. Kurbatov et al. (2006) detected volcanic signals during the last 12 ka in the Siple Dome A ice core from West Antarctica. They found that the number of volcanic sulfate signals is decreasing with age, possibly related to the relatively low sampling resolution in the deeper part of that core. Recently, Cole-Dai et al. (2021) applied the high-accumulation WAIS Divide ice core to determine a fairly constant Holocene eruption frequency with larger-than-Tambora (1815 AD) events 105 occurring approximately once per millennium.

Volcanic events identified in ice cores with tephra and sulfate peak synchronization
In historical times, the volcanic origin of an ice-core acidity spike may be pinpointed based on a precise dating of the ice core alone (Sigl et al., 2015). Further back in time, as the uncertainty of both the ice-core dating and the dating of the volcanic eruptions increases, the origin of a volcanic ice-core layer can only be determined if it 110 is associated with a volcanic ash (tephra) deposition in the ice (Gronvold et al., 1995). However, tephra layers are not always coinciding with sulfate peaks and not all tephra layers match volcanic sulfate signals (Davies et al., 2010). The ice-core volcanic source identification is important as it helps to constrain the magnitudeinterpreted here as sulfur emission strength rather than the mass of material erupted (Pyle, 2015) and the climate forcing of the eruption and, furthermore, it allows for a more detailed comparison to modelling studies.
In the last glacial period, many Greenland tephra deposits have been associated with Icelandic eruptions while around a dozen of identified tephra layers originate in North America and Eastern Asia Bourne et al., 2015;Davies et al., 2014;Cook et al., 2021 (in preparation)). In Antarctica, tephra layers have been identified and associated with eruptions occurring within Antarctica and in the Southern Hemisphere (Narcisi et al., 2005;Narcisi et al., 2010;Narcisi et al., 2012). Recently, Mcconnell et al. (2017) identified 120 tephra from the long-lasting and Halogen-rich Antarctic Mount Takahe eruption that occurred around 17.7 ka . Tephra of the Oruanui eruption from the Taupo volcano in present-day New Zealand has been identified and dated to 25.4 ka in the West Antarctic Ice Sheet Divide ice core (WDC) (Dunbar et al., 2017).
Volcanic eruptions generally do not deposit tephra in both Greenland and Antarctica, so the bipolar 125 synchronization of sulfur spikes in the ice cores is dependent on an alternative matching technique. Svensson et al. (2020) applied annual layer counting in both Greenland and Antarctic ice cores to match patterns of volcanic eruptions leading to the identification of some 80 bipolar eruptions in the 60-12 ka interval. For the Holocene, a bipolar synchronization of volcanic eruptions was released with the AICC2012 time scale (Veres et al., 2013). It has recently become possible to test if sulfate has indeed reached the stratosphere, which is a prerequisite for 130 being globally distributed, as the sulfate undergoes characteristic isotope fractionation in the stratosphere (Burke et al., 2019;Gautier et al., 2018;Crick et al., 2021;Baroni et al., 2008), but these analyses are still scarce for the Glacial.

135
Here we extend the ice-core record of sulfate deposition in Greenland and Antarctica by employing sulfate sulfate (or sulfur in one case; sulfate is mentioned for brevity) records from three Greenland and three Antarctic ice cores in the interval 60-9 ka. We investigate the sulfur emission strengths (i.e. defining the climate impact potential) and the frequency of volcanic eruptions detected in either Greenland or Antarctica. For eruptions identified in both hemispheres, we estimate the climate forcing against modern analogues (Pinatubo 1991AD) 140 and determine the occurrence of very large eruptions.
The sulfate records were obtained by different methods and with different temporal resolution as detailed in the 150 Table S1. For NGRIP the SO4 2concentration was obtained by a Continuous Flow Analysis system (CFA) (Bigler, Thesis, 2004). For WDC, sulfur was measured by ICP-MS coupled to a CFA system (Sigl et al., 2016).
For EDC and EDML, SO4 2was measured by Fast Ion Chromatography (FIC) coupled to a CFA system (Severi et al., 2015). For GISP2 and for a second NGRIP profile, SO4 2was obtained from discrete samples by Ion Chromatography (IC) Mayewski et al., 1993;Siggaard-Andersen, PhD thesis, 2004). The 155 temporal resolution of the sulfate records decreases with ice-core depth and ranges from sub-annual for the CFA profiles to decadal for the lower-resolution discrete profiles (Table S1 and Fig. S1). The high-resolution NGRIP and WDC records have been resampled to annual resolution, that we see as an upper limit for the effective resolution during the oldest part of the analysis in last glacial. As discussed in section 3.4, we are also resampling these records to lower resolution in order to obtain comparable resolution throughout the 160 investigated period and in order to investigate eruption occurrences among stadial and interstadial periods.
Minor data gaps were interpolated using linear interpolation and larger gaps are indicated as missing data. Sulfur and sulfate records were corrected for the sea-salt-sulfur contribution based on sodium concentrations as sea salt tracer by assuming 0.084 for the ratio of Na/S (Bowen et al., 1979), except for EDC. In all cases, the sea-salt correction is less than 15% of the measured sulfate background signal, giving a slightly change to the non-sea-165 salt-sulfur. Methanesulfonic acid (MSA) was subtracted from the Antarctic non-sea-salt-sulfate signal in a minor correction (Cole-Dai et al., 2021). In Greenland, the average concentration of MSA over the Holocene and last glacial period is low . As there are no continuous MSA records available for the investigated period, we do not make corrections by the MSA records to obtain the non-sea-salt sulfate.
For the 12-9 ka interval, the sulfate deposition is based on the NGRIP and GISP2 sulfate records of Greenland, and the WDC, EDML and EDC sulfate records of Antarctica. Besides the sulfate records, the depth assignment of the volcanic peaks was occasionally assisted by application of other high-resolution records that also are indicators of volcanic sulfate deposition. Those include the Electrical Conductivity Measurement (ECM) profile, the Di-electric Profiling (DEP) and the liquid conductivity record, that are also available for most of the applied cores (Table S1 for a complete list of records and 175 references).

Background signal determination and volcanic peak detection
To determine the biogenic sulfate background level of the ice-core sulfate records and detect volcanic peaks above this background, we apply methods similar to those applied for Holocene records (Fischer et al., 1998;Karlof et al., 2005;Gao et al., 2007;Sigl et al., 2013). A running median filter with a width of 50-180 years was 180 applied to estimate the non-volcanic background signal. Due to the different depth resolution of the individual records, and because of the highly variable and abruptly changing background levels of the Greenland sulfate levels across DO-events, it has been necessary to apply different filter widths for different records ( Fig. S10 (ay) and Table S1). A second iteration of a reduced running median (RRM) filter using the same filter widths was applied after removal of the spikes identified in the first iteration. For the coarse-resolution NGRIP IC sulfate -185 record, a continuous background determination was not possible and the background has been estimated manually for selected events only.
For all records, a volcanic detection threshold was estimated as RRM plus three times the running median of absolute deviations from the RRM (RMAD) using the same window widths as for the background determination ( Fig. S10 (a-y)). In Greenland, the RMAD value varies strongly across DO-events due to the much greater 190 background variability during stadial periods as compared to interstadials (Table S2). This implies that the volcanic detection may be compromised for the time windows covering the fast transitions from stadials and interstadials (lasting on the order of decades to up to about 150 years, Capron et al. (2021)) when the background variability changes in a short time. Thus, the volcano frequency found in these short intervals may be impacted, but the volcanic frequency outside of these short intervals is not affected. The duration of a 195 volcanic event is estimated from the fraction of the sulfate spike that is above the RRM and the sulfate deposition associated with the event is calculated by integrating the sulfate concentrations above the RRM across the duration of the event. The volcanic sulfate peak area S is calculated as follows: (1) https://doi.org/10.5194/cp-2021-100 Preprint. Discussion started: 9 August 2021 c Author(s) 2021. CC BY 4.0 License.
where y is the nss-sulfate, 0.917 is the ice-water density ratio, and the sulfate layer in the ice core is constrained 200 between the depths D1 and D2. At shallow depth, the temporal duration of the sulfate deposition can be determined precisely, but in the last glacial period peak broadening by diffusion is hampering such a determination.

Correction for ice flow / layer thinning
As the volcanic layer is buried in the ice sheet, the layer is being thinned by ice flow. In order to calculate the 205 amount of sulfate deposited on the ice sheet at the time of the eruption a correction for the thinning of the volcanic layer in the ice must be applied. (2) Where SF is the accumulated sulfate flux and T is the layer thinning. The thinning or the ratio between the original deposition thickness and the layer thickness in the ice core has been calculated by site-specific ice-flow 210 models that we applied here (Table S1 and Fig. S1). During the last glacial period, the thinning of the ice can be significant. For the Greenland cores, the thinning ranges from a 60% reduction of the layer thickness at 12 ka b2k to as much as 90% at 60 ka b2k (Table S1). For Antarctica, the thinning is most significant for the highaccumulation WDC core, where it approaches 95% at 60 ka b2k (Fudge et al., 2016;Buizert et al., 2015), whereas the EDC core exhibits less thinning of 30% at 60 ka b2k (Fig. S1). For the WDC core, we apply the 215 thinning function of Fudge et al. (2016) for the upper 2800 m, whereas below 2800m we apply a simple linear fit to the gas-based thinning function (Buizert et al., 2015) that has large and possibly unrealistic wiggles (Fig.   S1). As the volcanic sulfate deposition scales with the amount of thinning, an inaccurate thinning factor has large implications for the calculated sulfate depositions.
It should be noted that the volcanic sulfate flux calculation assumes that the sulfate deposition is mainly 220 occurring as dry deposition in the last glacial period. For high accumulation sites, such as those in Greenland and coastal Antarctica, it is quite likely that wet deposition is dominating at present day (Kreutz et al., 2000;Schupbach et al., 2018). We are not attempting to resolve this issue here, but we are aware that it may have a bias in our estimation of volcanic sulfate flux.

225
As the depth resolution of the sulfate records for the GISP2 and NEEM cores is relatively low (Fig. S1), adjacent acidity peaks may be merged into falsely large acidity spikes. We made a manual correction for this effect by comparing to the corresponding higher-resolution ECM and DEP records of the same core and removed or split falsely large peaks according to the associated ECM or DEP peaks.

230
We generate three lists of volcanic sulfate deposition for the 60-9 ka period: one for eruptions identified in the Greenland ice cores (Table S3), one for eruptions identified in the Antarctic ice cores (Table S4), and one for global (bipolar) eruptions identified in both Greenland and Antarctica (Table S5). The bipolar list contains the large global eruptions identified by Svensson et al. (2020) and includes two large bipolar eruptions of Veres et al. (2013) for the 12-9 ka period. Furthermore, two additional bipolar eruptions at 44.75 ka b2k and 44.76 ka 235 b2k have been included, applying the same methods as described in Svensson et al. (2020). For the bipolar list, there is no bipolar eruptions have been attributed in the 24-16 ka period, where ice core synchronization is difficult, so the list covers an effective period of 42 ka. The list of Northern Hemisphere (NH) volcanism contains all volcanic deposits larger than 20 kg km -2 in any of the applied Greenland ice cores, and the list containing the Southern Hemisphere includes all deposition events larger than 10 kg km -2 for the Antarctic ice 240 cores. This cut-off is necessary because of the highly variable sulfate background in the Greenland cold periods of the last glacial predominatly associated with mineral dust (e.g. gypsum (Svensson et al., 2000)) (Table S2). A sulfate deposition of 20 kg km -2 corresponds to half the Greenland deposition from the 1815 AD Tambora eruption, thus refers to quite large events in terms of total sulfur injections into the atmosphere. If the cut-off is set at a lower value, a large number of presumably non-volcanic spikes will be identified in stadial periods. A 245 similar cut-off applied to the sulfate deposition record of the last 2000 years  reduces the number of identified eruptions in Greenland from 138 to 49 events. Likewise, for Antarctica, we are only identifying events that are referred to as 'large' or 'very large' by Cole-Dai et al. (2021).
When eruptions are detected in several ice cores in Greenland or Antarctica, we provide the sulfate deposition value for each core, the range spanned by all the cores, and the calculated average sulfate deposition value. The 250 average Greenland sulfate deposition is calculated by the simple mean of the three cores due to their similar sulfate-deposition levels and proximity of the ice coring sites (Fig. S2). In Antarctica, however, the ice-core sites are far apart and the EDC core generally has the lowest sulfate deposition due to its low accumulation.
Therefore, when there is only a sulfate signal present in one or two cores, the average sulfate deposition of https://doi.org/10.5194/cp-2021-100 Preprint. Discussion started: 9 August 2021 c Author(s) 2021. CC BY 4.0 License.
Antarctica is calculated with a rescaling factor, similar to the method applied by Gao et al. (2008). The scaling 255 factor is based on the relative ratio of sulfate deposition from the 30 largest eruptions (in terms of the sulfate deposition) in the three cores in the period of 60-9 ka. For the 30 largest eruptions identified in Antarctica, the ratio of the sulfate deposition among the three Antarctic cores is EDC:EDML:WDC =0.72:0.87:1.4 (Fig. S2) .
The error estimation of volcanic sulfate deposition is according to the propagation of formula (1) and (2).
Considering the sulfate diffusion and layer thinning of ice cores, we derived the area of sulfate peak and then 260 estimated the uncertainty of the area of sulfate peak and layer thinning. We apply an estimated 21% uncertainty of the sulfate concentration levels independent of the analytical technique, an estimated 10% uncertainty of the applied thinning functions, and a 15% uncertainty related to the varying temporal resolution of the sulfate records. This leads to an error estimate of up to 26% for individual deposition events that we apply as error estimate for all the derived sulfate depositions.

Latitudinal band assignment of bipolar volcanic eruptions
The volcanic sulfate deposition in Greenland and Antarctica shows a distribution pattern (Table S5 and Table   S6), similar to that derived from the aerosol-climate modeling of volcanoes over past 2500 years (Martshall et al., 2020). To estimate the latitudinal band of bipolar volcanic eruptions of unknown origin, we applied the Support Vector Machine (SVM) classification model of Gapper et al. (2019) and Vapnik (1998). The 270 classification model is based on a kernel function generation and logistic regression (Fig. S3). For an optimal classification, a maximum-margin hyperplane was used to separate two classes. The kernel scale and box constraints were chosen for the model and a Bayesian optimization was used to optimize the above two parameters to yield the best classification model (Fig. S3). The model was trained using 21 eruptions for which the eruption site is known from tephra deposits in the ice (Table S6), and the bipolar eruptions of unknown 275 origin were predicted into two latitudinal bandsabove 40°N (NHHL) and below 40°N (LL or SH) (Table S5).
Due to the low number of known volcanoes erupted in the high latitudes of Southern Hemisphere, the method does not allow identification of eruptions potentially located in this region.

280
The sulfate deposition records derived from the NGRIP, NEEM and GISP2 ice cores are displayed in Figure 1 and Figures S9 (a-y). The background level of the sulfate signal in the Greenland cores is in the range of 50-200 https://doi.org/10.5194/cp-2021-100 Preprint. Discussion started: 9 August 2021 c Author(s) 2021. CC BY 4.0 License. ppb with both the absolute level and the signal variability being much higher in stadial periods than in interstadial periods (Table S2). At the abrupt climate transitions during the last glacial period the sulfate background changes abruptly, sometimes making it difficult to determine the background level and volcanic 285 detection limits (see above). The higher background level and variability during the stadial periods is partly due to the lower snow accumulation during these periods.
Using a 20 kg km -2 deposition threshold 1113 volcanic events are identified in Greenland in total (Table S3). Table S8 shows the number of volcanoes detected from one, two or three Greenland ice cores, respectively.
NGRIP has the highest event detection rate because of its superior depth resolution that becomes increasingly 290 critical in the cold climate of the investigated period where annual layers get down to a few centimeters of layer thickness (Fig. S1d).
The difference in sulfate deposition among the cores reflects to some degree a spatial variability in sulfate deposition within Greenland, but is also strongly influenced by the sulfate record type (IC, FIC, CFA), the data resolution, the thinning correction made for each site, and possibly the accumulation rates. A comparison of the 295 sulfate deposition between NGRIP, NEEM and GISP2 for the 30 largest volcanoes with identified signals in Antarctica ice cores shows that there is a very large spread in the signal from event to event ( Fig. S2 (a-d)). Due to the large variability of sulfate deposition values among the cores, we provide the full range of sulfate deposition values spanned by the cores as well as the average sulfate deposition value, when an eruption is identified in several ice cores. 300 Fig. S2 (e) shows a comparison between the 57 largest volcanic events identified in the NGRIP ice core, as derived from the high-resolution CFA SO4 2record and from the lower resolution IC SO4 2record, respectively (Table S3). This comparison illustrates the differences in calculated sulfate depositions due to measurement technique and record resolution only, as the sulfate deposition and the thinning factor uncertainties are eliminated. It is obvious from the comparison that the derived sulfate deposition of a single event in a single 305 core should not be taken at face value; there are very large uncertainties on top of the variability caused by the spatial distribution pattern. Still, being the first systematic compilation of large volcanic eruptions from the last glacial period, the estimated depositions give a good approximation about the order of sulfur emission strength and the frequency of large volcanic eruptions in Greenland.
Out of the 1113 volcanic events identified in Greenland we identify 10 very large events with sulfate deposition 310 estimated to be larger than 300 kg km -2where 9 have an Antarctic counterpart indicating their global impact.
We identify 87 volcanic events with sulfate deposition larger than 135 kg km -2 , thereby exceeding the Icelandic Laki eruption occurring in 1783 CE, which produced the largest Greenland sulfate deposition in the last 2500 https://doi.org/10.5194/cp-2021-100 Preprint. Discussion started: 9 August 2021 c Author(s) 2021. CC BY 4.0 License.
years (Sigl et al., 2015). The largest sulfate deposition in Greenland is by far the North Atlantic Ash Zone (NAAZ II) event occurring at 55.38 ka b2k (Rutledal et al., 2020;Austin et al., 2004). The sulfate deposition of 315 this event is in the range of 866-1436 kg km -2 , which is more than a factor of two higher than any other event occurring in the investigated period. The second largest event is an unknown eruption occurring at 45.55 ka b2k with a sulfate deposition in the range 200-849 kg km -2 . Those and other large events will be discussed in more detail in section 4.4.

320
The volcanic sulfate deposition extracted from the EDML, EDC and WDC ice cores are shown in Fig. 1 and Fig. S10 (a-y). The sulfate background level in Antarctica is comparable to that of Greenland (Table S2), but there is much less signal variability and an absence of abrupt shifts associated with Greenland DO transitions.
In total, 740 (350) volcanic sulfate deposition values estimated to be larger than 10 (20) kg km -2 are identified in the Antarctic ice cores (Table S4 and Table S8). WDC has the highest accumulation of the three records, but as 325 the deepest part of the applied WDC record is close to bedrock, layer thinning becomes very significant here (Fudge et al., 2016) (Table S1). Therefore, in the deepest section of the WDC core, fewer eruptions are detected as compared to EDML and EDC, and the derived WDC sulfate deposition has a very strong dependency on the applied thinning function. EDC has the lowest accumulation of the three and also the lowest temporal resolution of the sulfate record. Therefore, smaller events are generally not detected in the EDC core, and the low 330 accumulation most likely causes some eruptions to be partly or entirely absent from the record (Gautier et al., 2016).
The volcanic deposition of the 30 largest events in Antarctica are compared in Fig. S2. In general, there are large differences among the sulfate deposition in the different cores for the same events, primarily owing to the large gradient of accumulation rates over Antarctica (Sigl et al., 2014). The influence of different measurement 335 types (FIC or CFA-ICP-MS) and the temporal resolution of the records may also accounts for some of the variability, as well as post-depositional processes (Gautier et al., 2018).
In the 60-9 ka period, there are 32 large volcanic events with a sulfate deposition larger than 78.2 kg km -2 , which is the largest volcanic deposition in Antarctica over the most recent 2500 years (Sigl et al., 2015). Fourtynine volcanic events have a sulfate deposition larger than 63.6 kg km -2 corresponding to the Antarctic sulfate

The bipolar volcanic sulfate deposition record (60-9 ka)
For the 60-9 ka period, the 85 bipolar volcanic eruptions that have been identified in both Greenland and 345 Antarctic ice cores (Veres et al., 2013;Svensson et al., 2020) are listed separately in Table S5. We note that no bipolar eruptions have been identified in the interval 24.5-16.5 ka, because of the general difficulty to synchronize ice cores in this period (Seierstad et al., 2014). For the bipolar eruptions the volcanic forcing can be estimated using methods previously applied to Holocene eruptions (see section 4.3). Twenty-eight out of the 85 bipolar eruptions have been identified in all six ice cores. Two of the bipolar volcanic eruptions are known from 350 Greenland ice-core tephra deposits to be of Icelandic origin: The Vedde ash eruption (around 12.17 ka b2k) (Mortensen et al., 2005) and the NAAZ II eruption (around 55.38 ka b2k) (Rutledal et al., 2020). Furthermore, there are tephra deposits in Greenland from the Japanese Towada-H eruption (ca 15.68 ka b2k) (Bourne et al., 2016). In Antarctica there is tephra deposited from the New Zealand Oruanui, Taupo, eruption (ca 25.46 ka) (Dunbar et al., 2017) (Table S6).

Relationship between data resolution and the volcanic eruption record of the last glacial (60-9 ka)
To investigate the frequency of different size categories of volcanic eruptions over time, we need to consider the sample resolution of the underlying sulfate records (Fig. S1(a+b)). As we go back in time, the layer thinning becomes stronger which makes it increasingly difficult to detect smaller eruptions signals. Depending on the 360 ice-core sample resolution, the apparent number of smaller eruptions decreases relatively faster with increasing age (or depth) than the larger eruptions as the former will no longer exceed the detection threshold determined by the background signal due to smoothing of the records. The effect is seen when the frequency of detected eruptions is separated into fractions of sulfate deposition sizes. We separated the eruptions according to the 0.7 and 0.9 quantiles of the volcanic sulfate deposition distribution, which for Greenland are 68 kg km -2 and 140 kg 365 km -2 and for Antarctica are 25 kg km -2 and 50 kg km -2 (Fig. S4). When considering the frequency of volcanic eruptions across the investigated period, we should preferably have the same or at least a comparable temporal sulfate sample resolution throughout the period in order to minimize this age (or depth) bias.
For Greenland, only the NGRIP sulfate record has high (sub-annual) resolution for most of the investigated interval (Fig. S1b). Therefore, the number of eruptions detected in different periods will depend strongly on the 370 NGRIP sample resolution. The NGRIP sulfate record is measured in 1mm resolution, but the effective resolution is lower than that due to smoothing of the record in the ice andmost importantlyby signal dispersion during the measurement. Using an approach similar to that of Rasmussen et al. (2005), we performed a spectral analysis of the NGRIP sulfate record to determine a signal cut-off at around 2 cm, which is https://doi.org/10.5194/cp-2021-100 Preprint. Discussion started: 9 August 2021 c Author(s) 2021. CC BY 4.0 License.
comparable to that of other components analyzed with the same setup (Bigler et al., 2011) (Fig. S5). Two cm 375 corresponds to about 2 annual layers in NGRIP in the climatically coldest and oldest sections of the investigated interval. In Table S7 and Fig. S4 we investigate the effect of increased smoothing of the NGRIP sulfate record for the detection of sulfate spikes. As expected the number of detected eruptions decreases with increased smoothing, however, with 3 years smoothing, the number of larger sulfate spikes starts to increase as adjacent sulfate peaks are being merged.

380
As a compromise between compensating for the decreasing resolution of the record with depth and at the same time avoiding a strong influence of the peak merging effect, we smoothed the NGRIP sulfate record to 2-year resolution for the eruption frequency analysis, knowing that this is eliminating a fraction of the smaller eruptions in the younger section of the record. With the 2-year smoothing there is probably still a bias towards a higher detection rate in the stadial periods due to the high signal variability of the sulfate record in those periods (see 385 discussion in section 4.1). The GISP2 sulfate record has lower than 10-year resolution in the 60-14 ka period ( Fig. S1), so in this period we can only use the GISP2 sulfate record to estimate the sulfate deposition of the larger eruptions. The NEEM sulfate record has a constant resolution of 10 years and we apply this record throughout the investigated period. In Fig. 4 we show the records of detected Greenland eruptions per millennium with the sulfate deposition grouped into three size fractions.

390
For Antarctica, EDML and EDC have fairly constant thinning rates for the 60-9 ka period (Fig. S1c), whereas WDC has very strong thinning in the oldest part of the period, due to the high accumulation rate (Fig. S1c).
Spectral analysis of the Antarctic sulfate records shows the effective resolution (Table S1) and the signal cut-off is around 2 cm, 3 cm, and 1 cm for EDML, EDC and WDC, respectively (Fig. S5). In order to have comparable resolution throughout the investigated period, we smoothed the entire WDC sulfate record to 2-year resolution 395 to homogenize the record resolution across the investigated period. For EDML the correspondingly effective signal resolution are 1-2 year for the entire period, so we keep the original resolution of the EDML sulfate record for the entire period. The EDC sulfate record has almost constant resolution, so we apply this record throughout the investigated period. The result of the frequency analysis is discussed in section 4.1. al., 2018). In particular, the melting of the large ice sheets at the end of the last glacial period and the corresponding sea-level rise are thought to have created crustal stress imbalances and increased the volcanic 405 eruption frequency in the deglacial period (Watt et al., 2013;Huybers and Langmuir, 2009;Zielinski et al., 1996;Zielinski et al., 1997).

Millennium scale volcanic eruption variability
We investigate the variability in eruption frequency and sulfur emission strength with the DO cycles by separating the detected volcanic eruptions according to climate of 'cold' and 'milder' periods, which the onset 410 and termination of DO events defined by Rasmussen et al. (2014). 'Cold' periods constitute the stadial periods from 21 ka b2k to 60 ka b2k, covering a total of 21.8 ka period. 'Mild' periods are the remaining periods that include all of the interstadial periods from 21 ka b2k to 60 ka b2k, adding up to 17.2 ka years. We grouped the detected eruptions according to the size of deposited sulfate and also investigated the effect of smoothing the NGRIP sulfate record to 1, 2, 3, and 5 years resolution (Fig. S6). This observed difference leaves us two options: either there is a true millennium-scale volcanic eruption variability related to northern hemispheric climate variability, or the dependency is an artefact caused by the very different properties of the Greenland sulfate record in cold and milder periods. In Greenland, the variability of the sulfate background level is high in cold periods, stadials and the LGM, and low in milder periods such as 425 interstadials (see section 3.1). We try to take this difference into account by employing the variability-dependent RMAD to determine the eruption detection limit, but because of the very different character of the sulfate profiles in cold and milder periods there is still the possibility of having a climate related detection bias of sulfate spikes that is unrelated to the actual volcanic eruption record. For example, in cold periods, it could be that at strong seasonal sulfate spikes associated with marine or terrestrial (non-volcanic) sulfate production 430 could be detected as smaller volcanic eruptions. For recent times, Gao et al. (2007) has demonstrated that the volcanic sulfate deposition to some degree is accumulation dependent, with increased accumulation leading to higher sulfate deposition. Since, however, the highest accumulation occurs in the interstadial periods, this effect https://doi.org/10.5194/cp-2021-100 Preprint. Discussion started: 9 August 2021 c Author(s) 2021. CC BY 4.0 License. would lead to relatively higher sulfate deposition in interstadial periods, which is the opposite of what is observed.

435
We investigated the effect of a possible seasonal sulfate contributions of non-volcanic origin by a smoothing exercise where successively increased smoothing of the NGRIP sulfate record increasingly eliminates the detection of smaller short-lived (seasonal) sulfate spikes. Fig. S4, Fig. S6 and Table S7 show the number of sulfate spikes detected in cold and milder periods as a function of the degree of smoothing. It is seen that for cold periods, the number of detected small and short-lived sulfate spikes decreases significantly with increased 440 smoothing, whereas the effect is almost absent for milder periods. This suggests that the higher number of detected sulfate spikes in cold periods is mainly due to spikes of periodic duration of sulfate background which we believe are less likely to be of volcanic origin. For volcanic eruptions that are large enough to inject sulfate into the stratosphere modern observations show that the sulfate fallout typically will last several years. We thus conclude that the apparent increased volcanism detected in cold periods is most likely an artefact caused by non-445 volcanic seasonal contributions from high-sulfate marine or terrestrial sources. It could for example be contributions from gypsum that is known to have high deposition fluxes in Greenland during stadial periods (Legrand and Mayewski, 1997). Indeed, we do observe a higher correlation in stadial periods than interstadial periods between the volcanic sulfate deposition and both insoluble dust concentration and sodium concentration for bipolar events (Fig. S9). We cannot exclude, however, that other mechanisms such as different atmospheric 450 circulation patterns, or different sulfate deposition process in terms of dry versus wet deposition, are also affecting the Greenland volcanic detection rates in stadial versus interstadial periods. We note that because an increased smoothing of the sulfate record mostly decreases the number of smaller detected sulfate depositions, the effect of smoothing on the large accumulated sulfate deposition over time is less important (Fig. S4).

455
In order to estimate the long-term pattern of volcanic activity, the accumulated number of volcanic eruptions and resulting sulfate deposition in Greenland and Antarctica over the investigated period are compared in Fig. 3 (same data as applied in Figure 1). The cumulative volcanic number (sulfate deposition) in Greenland is larger than in Antarctica by a factor of 3.2 (4.7), considering only sulfate depositions larger than 20 kg km -2 at both poles. We attributed this mainly to the larger number of volcanoes situated in the NH, but it is possibly also 460 related to the proximity of Iceland to Greenland and other volcanic regions upwind of Greenland, and to the differing atmospheric circulation patterns between the two hemispheres (Toohey et al., 2019). In terms of cumulated volcanic number there is some variability that could be climate related, but overall the eruption https://doi.org/10.5194/cp-2021-100 Preprint. Discussion started: 9 August 2021 c Author(s) 2021. CC BY 4.0 License. frequency is more constant as compared to the geological record (Brown et al., 2014), most likely because the ice-core preservation and identification are much more homogeneous over time than that of radiometrically 465 dated eruptions.
The number of eruptions per millennium shows a fairly constant level both for Greenland and for Antarctica ( Fig. 4, Fig. S7 and Fig.S8). The Early Holocene increase in global volcanism as suggested by Watt et al. (2013) is not identified in our record. Neither is the periodicity of Campanian eruptions as identified by Kutterolf et al. (2019) but this pattern could be masked by eruptions from elsewhere.

470
Based on the EDC sulfate record, Castellano et al. (2004) determined the Antarctic volcanic sulfate deposition to be fairly constant over the last 45 ka, except for an increase during the last 2000 years. The eruption frequency level of 5-10 eruptions per millennium determined in that study is on the low side compared to the present study, in particular for the older part of the period (Fig. 4). This is likely because the EDC sulfate record has low resolution and could not pick up all of the smaller short-lived eruptions, also because the low 475 accumulation at the EDC site does not record some eruptions (see methods/results section). Over the 60-9 ka interval we have detected 224 volcanic eruptions in EDC as compared to 513 and 496 in the higher resolution WDC and EDML records, respectively (Table S4) Based on the same GISP2 sulfate dataset as applied in the present study and after introducing a rescaling to compensate for the decreasing sample resolution with depth, Zielinski et al. (1996Zielinski et al. ( , 1997 identified an enhanced 485 northern hemispheric volcanic eruption frequency in the deglacial period (17-6 ka) compared to the 60-0 ka average. In our 3-core Greenland sulfate deposition compilation, we also identified an increased volcanic eruption frequency over the deglacial period, which is 34% higher than the Holocene period (2-0) for the average volcanic number per millennium ( Fig. 2 and Fig. 4). In particular, there is an increased fraction of large sulfate deposition events in the deglacial period. In Antarctica, neither the detected eruption frequency nor the 490 total sulfate deposition over the deglacial period are higher than the average for the full 60-9 ka interval. We therefore conclude that the increased number of large eruptions occurring in the deglacial period are mostly of northern hemispheric origin. This hypothesis is supported by the elevated number of tephra deposits identified in Iceland records (Maclennan et al., 2002) and in Greenland during the deglacial period mostly of Icelandic https://doi.org/10.5194/cp-2021-100 Preprint. Discussion started: 9 August 2021 c Author(s) 2021. CC BY 4.0 License.
origin, a few from Japanese and unknow origins (Bourne et al., 2016;Zielinski et al., 1997;Mortensen et al., 495 2005;Gronvold et al., 1995). Van Vliet-Lanoe et al. (2020) found that the enhanced volcanic events are most likely related to stress unlocking during deglaciation events, which the occurrence of volcanism in Iceland is related to rifting events at the melting margin of the ice sheet. In the period of 14.6-13.1 ka in Southeast Alaska, the enhanced volcanic activity from the Mount Edgecumbe Volcanic Field has also been observed by Praetorius et al. (2016). The frequency of detected eruptions for the glacial period (60-21 ka) is similar to that of the last 500 2000 years, when comparing the same sulfate deposition fractions (Fig. 4). We note that for the last 2000 years only eruptions with sulfate depositions larger than 20 kg km -2 are included and the sulfate records have been smoothed with a two-year filter. The apparent increase in global volcanism over the last millennia in the geologic records (Brown et al., 2014) is therefore most likely due to an under-representation of older eruption identifications in that study (Papale, 2018;Deligne et al., 2010). When it comes to sulfur emission strengths of 505 eruptions and accumulated sulfate deposition, however, the last 2000 years appear under-represented in very large eruptions as compared to long periods of the last glacial period (Table 2).

Estimating the volcanic forcing of the last glacial 60-9 ka
To estimate the volcanic forcing from eruptions occurring in the last glacial and early Holocene we applied the approach of Gao et al. (2007), Crowley and Unterman. (2013) and Hansen et al. (2005). This method was 510 calibrated and evaluated using modern volcanic eruptions (Gao et al., 2008;Crowley and Unterman, 2013), and it is quite likely that the results of the glacial eruptions will be biased by the very different glacial conditions. Calibrated against Pinatubo 1991 (at 15°N), the method is best suited for bipolar eruptions, and it is necessary to separate NH high-latitude eruptions from other eruptions, as they require different scaling functions to relate ice-sheet wide sulfate deposition rates to stratospheric aerosol loading (Gao et al., 2007). We defined NH high 515 latitude eruption as eruptions that occurred at a latitude above 40° N. To identify the NH high latitude eruptions, we applied a Support Vector Machine learning classifier model (SVMsee methods section), that is being trained by the bipolar sulfate deposition of volcanic eruptions for which the eruption site is known. We applied 17 Holocene and 4 glacial volcanic eruptions of known origin (Table S6) to predict that 51 out of 85 bipolar eruptions of unknown origin are likely to have occurred in the NH high latitudes (Fig. 5 and Fig. 6). The Gao et 520 al. (2007) approach provides an estimate of the global stratospheric sulfate aerosol loading that is then linearly scaled to an atmospheric optical depth (Crowley and Unterman, 2013) and radiative forcing using the scaling factor of Hansen et al. (2005). The obtained radiative forcings are directly comparable to those obtained for the last 2500 years (Sigl et al., 2015). Having estimated the volcanic forcing of the glacial and early Holocene eruptions, we can compare to the well-525 constrained eruptions that occurred over the last 2500 years (Sigl et al., 2015). In Table 2 we compared the number of eruptions with a climate forcing larger than Oruanui, Taupo, 25.5 ka; the unidentified eruption of 426 BC and Tambora, 1815 AD for the periods of 2.5-0 ka, and 60-9 ka. Note that the investigated period of the last glacial and the early Holocene covers some 43 ka (60-9 ka minus the section of no identified bipolar eruptions 16.5 to 24.5 ka). Among the 85 bipolar events investigated in the last glacial and early Holocene, three eruptions 530 are found to be larger than Qruanui, Taupo, and 27 eruptions are found to be larger than the largest eruption of the past 2500 years.
Based on geological evidence, the occurrence of large volcanic eruptions was estimated in terms of their VEI (Volcanic explosivity index), which is based on the volume of ejected magma. The average occurrence rate of VEI-7 (VEI-8) volcanoes discovered in the last 125 ka BP (2600 ka BP) is 0.3 (0.01) per millennium, 535 respectively, based on databases of LaMEVE and GVP (Papale, 2018). From ice-core records, we cannot directly derive the VEI index, but the volcanic forcing estimated is likely to be another good indicator of volcanic eruption magnitude. The average rate of volcanoes in the 60-9 ka interval with a climate forcing larger than the 1815 AD Tambora eruption (VEI-7) is 1.60 per millennium (Table 2). In the same interval, we have 4 eruptions with a climate forcing larger or equal to the 25.47 ka Taupo Oruanui (VEI-8) eruption. Of those, the 540 Icelandic 55.38 ka NAAZ II volcanic event is likely to be overestimated due to its proximity to Iceland (see section 4.4.) and it may not be of VEI-8 magnitude. Thus, we identify three VEI-8 sized eruptions occurring at a rate of 0.07 eruptions per millennium. Those ice-core based eruption frequency estimates are thus about 5.3 (7) times larger than those based on geological evidence for VEI-7 (VEI-8).

545
The largest eruptions of the last glacial period and early Holocene are listed by the average climate forcing in Table 1. The highest ranking on the list is the Icelandic eruption associated with the NAAZ II occurring in GS-15.2 at 55.4 ka b2k (Gronvold et al., 1995;Zielinski et al., 1997). The sulfate deposition in Greenland from this eruption is enormous (almost 15 times higher than the deposition for Laki 1783 AD), but none or only very little sulfate is deposited in Antarctica. With its close proximity to Greenland, it is therefore likely that a large 550 fraction of the sulfate was transported to Greenland through the troposphere and thus the estimated climatic forcing of -(99.7-228.0) W m -2 is probably overestimated and should not be taken at face value. From the marine sediment records, the NAAZ II layer consists of different Icelandic events (Rutledal et al., 2020). For the https://doi.org/10.5194/cp-2021-100 Preprint. Discussion started: 9 August 2021 c Author(s) 2021. CC BY 4.0 License.
above reasons, the eruption should not be considered as the largest eruption of the last 60 ka in a global context.
The very wide range of sulfate deposition of 866.2-1435.6 kg km -2 among the different Greenland ice-core 555 records (Table S5) is probably due to a combination of 1) a true geographical variability in deposition in Greenland, 2) the quite significant and possibly inaccurate correction for layer thinning at great depth in the ice cores for this eruption, and 3) an unexplained difference related to the different analytical method applied (see results section).
The second largest event on the list occurred in GI-12 at 45.56 ka b2k and is identified in all Greenland and 560 Antarctic ice cores, but the deposition at EDC is very low, most likely due to its low accumulation and the low time-resolution of the sulfate record. The eruption is of unknown origin, but based on its relative Greenland and Antarctic sulfate deposition, it occurred in the NH extratropical area and its climate forcing is estimated to the range of -(71.8-176.5) W m -2 .
Number three on the list of large eruptions occurred at 38.13 ka b2k; 100 years after the onset of GI-8 and 11 565 years before the occurrence of the Faroe Marine Ash Zone III (FMAZ III) tephra in Greenland . The eruption is detected in all ice cores (but there is a data gap for the WDC sulfur in this interval), and with very similar Greenland and Antarctic sulfate depositions this is most likely a low-latitude eruption with an estimated average climate forcing of -82.8 W m -2 , about 4 times that of Tambora (1815 AD).

570
The forth largest bipolar ice-sheet sulfate deposits are from the New Zealand Taupo, Oruanui, eruption that occurred in GS-3 at 25.46 ka b2k, verified by the presence of a tephra deposit (Dunbar et al., 2017). The estimated climate forcing of the eruption is in the range of -(33.2-118.1) W m -2 , about 2-3 times higher than Samalas 1257AD (-32.8 W m -2 ) (Sigl et al., 2015), which caused large variability of atmospheric temperature of global climate for several years (Lim et al., 2016), prolonged the drought in North American and leading to the 575 famine (Herweijer et al., 2007). Despite its magnitude, the sulfate deposition of the eruption is undetected in the Greenland GISP2 ice cores with the lowest temporal sulfate-data resolution in this time period.
In GS-5.1 at 29.68 ka b2k there is a pair of volcanic eruptions separated by some 40 years, of which the younger event ranks 5 on the list of the largest eruptions. The estimated climatic forcing of this eruption is -65.5 W m -2 and it has a clear NH signature in the bipolar sulfate distribution. There is no associated tephra deposit in the ice volcanic eruption of that magnitude within a range of several millennia. The Y-3 eruption is independently Carbon-14 dated to 29.0 ka b2k (Albert et al., 2015).
Number 7 on the list of the largest eruptions occurred at 46.68 ka b2k in GI-12 and is the only eruption in this study apart from the 25.46 ka b2k Taupo, Oruanui, eruption that has a clear southern hemispheric origin. The

585
Antarctic sulfate deposition of the eruption is approximately twice that in Greenland (Fig. 6) Table 1. Both eruptions have bipolar sulfate distributions suggesting a NH 590 eruption. Because of their magnitude and their stratigraphic setting right at the onset of GS-9 these volcanic events are both possible candidates for the Italian Y-5 Campanian Ignimbrite eruption. There is no tephra evidence for this suggestion in the ice cores, but tephra from this eruption has been identified in the Black Sea in very similar stratigraphic setting and the eruption is independently dated by Ar/Ar to 39.9 ± 0.1 ka b2k (Giaccio et al., 2017).

595
In the late GI-1 right before the onset of GS-1 / Younger Dryas there is a quadruple of bipolar eruptions covering a period of 110 years (Svensson et al., 2020). The oldest (13,028 a b2k; -44.1 W m -2 ) and the youngest (12,917 a b2k; -40.2 W m -2 ) of those eruptions are number 17 and 19 on the list of large eruptions. Both eruptions are likely to have occurred in the NH at above 40°N. The 12,917 b2k eruption has been suggested as a candidate for the Laacher See eruption (LSE) that occurred in the East Eifel region in present-day Germany 600 (Brauer et al., 1999;Baldini et al., 2018), but no tephra has been identified in the ice cores from this eruption. A recent publication, however, suggested that the LSE was around 130 years older (Reinig, et al., 2021) and is synchronous instead (within age uncertainty) with the oldest event of this quadruple of bipolar eruptions, but also with some minor sulfate signals.
The Vedde Ash layer (number 36) originates from a significant euption from Katla, Iceland with widespread 605 tephra deposition in Greenland (Mortensen et al., 2005;Gronvold et al., 1995) and the North Atlantic region in the middle of GS-1 / Younger Dryas at 12.17 ka (Lane et al., 2012). Surprisingly, this Icelandic eruption not only deposited large amounts of sulfate in Greenland but also a much smaller amount in Antarctica, where it is identified in the WDC and EDML sulfate profiles suggesting that it had a large stratospheric injection. The volcanic forcing is estimated to about -30.9 W m -2 taken with the reservation that some of the sulfate could have 610 made it to Greenland through the troposphere. https://doi.org/10.5194/cp-2021-100 Preprint. Discussion started: 9 August 2021 c Author(s) 2021. CC BY 4.0 License.

Conclusion
We have employed three Greenland and three Antarctic sulfate/sulfur ice-core records to document the global volcanic activity over the period 60-9 ka b2k. Detection of volcanic signals in the last glacial is challenging due to extensive layer thinning of the ice cores with depth and a highly variable sulfate background level in the 615 Greenland ice cores across DO events, among other factors. Due to those challenges, we limited ourselves to identify Greenland (Antarctica) eruptions with sulfate deposition larger than 20 (10) kg km -2 , which for Greenland is about half the volcanic sulfate deposition of the Tambora 1815 AD and for Antarctica comparable to the Pinatubo 1991 AD eruption.
With those restrictions we identified 1113 volcanic eruptions in Greenland and 740 volcanic eruptions in 620 Antarctica verifying the Northern Hemisphere (NH) to be the most volcanic active. Of those, 85 volcanic eruptions had global impact (bipolar volcanos) as they were previously identified in Greenland and Antarctic ice-core records. Compared to the past 2500 years, the ratio of bipolar volcanoes to total identified volcanic events is very low in last glacial period, highlighting the challenges of confident synchronization over the glacial period. Based on the hemispheric partitioning of sulfate deposition following well-known historical eruptions, 625 we determined if the bipolar volcanoes are likely to be situated above or below 40 ° N latitude. We then estimated their climatic impact in terms of radiative forcing based on established methods. Throughout the investigated period, we find that 69 volcanoes are larger than the Tambora 1815 AD eruption, and one unknown volcano occurring at 45.56 ka b2k in the NH and one unknown volcano at 38.13 ka b2k in the low latitude or Southern Hemisphere are larger than the Taupo, Oruanui, eruption occurring at 25.46 ka b2k in present day New

630
Zealand. The Icelandic NAAZ II eruption (55.38 ka b2k) has by far left the largest sulfate deposition in Greenland, but due to its minor sulfate deposition in Antarctica, it is thought that only a fraction of the sulfate was injected to the stratosphere. In general, we observe significantly higher occurrences of very large eruptions (VEI 7 or larger) than estimated from the geological record, indicating that ice cores are more consistent than geological evidence for volcanic signals detecting.

635
Overall, the frequency of volcanic eruptions per millennium is rather constant throughout the investigated period and comparable to that of the most recent millennia. In agreement with previous studies, however, we find elevated levels of volcanic activity in the NH during the deglacial period (16-9 ka b2k). In particular, many very large eruptions occurred in this interval, quite likely associated with the redistribution of mass related to the melt down of the major glacial ice sheets in this period. An apparent increase of Northern Hemispheric volcanic activity in cold stadial periods as compared to milder interstadial periods is likely to be an artefact due to the highly variable nature of the Greenland sulfate signal of those periods.
Data availability: The high-resolution NGRIP CFA sulfate dataset is available as supplementary information for this publication. All other applied datasets are available elsewhere.   Table 1 The largest bipolar volcanic eruptions ranked by the volcanic aerosol forcing that is derived from the average Greenland and Antarctic sulfate depositions. The stratospheric aerosol loading was estimated with the same 1025 method as Gao et al. (2007) and the radiative volcanic aerosol forcing was established by scaling the atmospheric aerosol loading for Pinatubo 1991 (section 4.3). The 'Prediction of volcanic site' column is adapted from Table S5. Table 2 The number of volcanic bipolar eruptions of comparable magnitude to those of Oruanui, Taupo (25.4 ka b2k), 1030 the largest eruption in last 2500 years (426 BC) and the Tambora 1815 AD eruptions by the climate forcing.

Figure 1
The average volcanic sulfate deposition in Greenland and Antarctica (60-9 ka; 2.5-0 ka adapted from Sigl et al., 2013Sigl et al., , 2015. Blue line represents δ 18 O of NGRIP core. Volcanic sulfate deposition of Greenland is larger than 20 kg km -2 (a, cyan bars) and Antarctic volcanic sulfate deposition is larger than 10 kg km -2 (b, black bars).  As only the WDC ice core is included in the above period, a scaling factor is applied to obtain the average sulfate deposition in Antarctica. For the DG, G, GI, GS periods the number of volcanoes per millennium is detected in all Antarctic ice cores with WDC sulfate records smoothed to 2 years resolution. (b) The number of periods, NGRIP is applied in 2-year resolution, and GISP2 is smoothed in 5-year resolution in the period 14-7 ka and for the period older than 14 ka b2k only larger eruptions are included).  Greenland and Antarctica (60-9 ka). Volcanoes are classified into three groups: 20-68 kg km -2 , 68-140 kg km -2 , and >140 kg km -2 for Greenland and 10-25 kg km -2 , 25-50 kg km -2 , and >50 kg km -2 for Antarctica. smoothed sulfate record from GISP2 and NGRIP. For GISP2, the data was smooth with 5 years in the period of 14-9 ka and for the period older than 14 ka only large volcanic eruptions of GISP2 sulfate record is applied.
NGRIP sulfate data has been smoothed to 2-year resolution. From 12-9 ka, number and sulfate deposition of volcanoes per millennium is based on GISP2 of Greenland ice core and EDML, EDC and WDC of Antarctica ice core. For the latest period 911-1911 AD, 89-911 AD, volcanic sulfate deposition is defined from NEEM ice 1065 core of Greenland and WDC ice core of Antarctica ) and volcanic sulfate records of WDC and NEEM were detected with smoothed data of 2-year resolution. Volcanic list of Antarctica is detected with smoothed sulfate record of WDC in 2-year resolution.