Articles | Volume 20, issue 2
Research article
09 Feb 2024
Research article |  | 09 Feb 2024

State-dependent impact of major volcanic eruptions observed in ice-core records of the last glacial period

Johannes Lohmann, Jiamei Lin, Bo M. Vinther, Sune O. Rasmussen, and Anders Svensson

Recently, a record of large, mostly unknown volcanic eruptions occurring during the younger half of the last glacial period (12–60 ka) has been compiled from ice-core records. In both Greenland and Antarctica these eruptions led to significant deposition of sulfate aerosols, which were likely transported in the stratosphere, thereby inducing a climate response. Here we report the first attempt to identify the climatic impact of volcanic eruptions in the last glacial period from ice cores. Average negative anomalies in high-resolution Greenland and Antarctic oxygen isotope records suggest a multi-annual volcanic cooling. Due to internal climate variability, glaciological noise, and uncertainties in the eruption age, the high-frequency noise level often exceeds the cooling induced by individual eruptions. Thus, cooling estimates for individual eruptions cannot be determined reliably. The average isotopic anomaly at the time of deposition also remains uncertain, since the signal degrades over time as a result of layer thinning and diffusion, which act to lower the resolution of both the oxygen isotope and sulfur records.

Regardless of these quantitative uncertainties, there is a clear relationship of the magnitude of isotopic anomaly and sulfur deposition. Further, the isotopic signal during the cold stadial periods is larger in Greenland and smaller in Antarctica than during the milder interstadial periods for eruptions of equal sulfur deposition magnitude. In contrast, the largest reductions in snow accumulation associated with the eruptions occur during the interstadial periods. This may be the result of a state-dependent climate sensitivity, but we cannot rule out the possibility that changes in the sensitivity of the isotope thermometer or in the radiative forcing of eruptions of a given sulfur ejection may play a role as well.

1 Introduction

Several studies on ice-core and tree-ring records, as well as climate models, show that volcanism plays a major role in generating the climate variability observed in the Common Era (PAGES 2k Consortium2019). During this period, all of the most pronounced episodes of reduced tree growth in composite tree-ring records can be associated with large volcanic eruptions and their tropospheric cooling effect due to the ejection of sulfur aerosols (Sigl et al.2015). This suggests that volcanic eruptions are responsible for the strongest multi-annual summer temperature decreases in middle- to high-latitude regions of the Northern Hemisphere. On longer timescales, clusters of large eruptions coincide with centennial cold periods during the Holocene similar to the Little Ice Age, as shown in tree-ring (Helama et al.2021) and ice-core records (Kobashi et al.2017). In climate model simulations of the past millennium, the temperature variability due to volcanic forcing exceeds the variability due to solar forcing (Schurer et al.2014), as well as the internal multi-decadal variability (Mann et al.2021).

Large future eruptions are unpredictable hazardous perturbations, which may compound the increasing climate extremes that stress ecosystems and societies and which may increase risks of potential tipping points (Lenton et al.2008). However, the impact of very large eruptions on the climate is not understood in detail, and it may change with the average state of the climate. In particular, the climatic impact may differ from glacial to interglacial conditions or the warmer world of the next centuries. Modeling studies investigating eruptions under future warming scenarios have reported both an enhanced (Fasullo et al.2017) and a reduced (Hopfcroft et al.2018) volcanic cooling or a change in cooling that depends on the eruption magnitude (Aubry et al.2021). Another modeling study found no evidence for a difference in the global temperature response during the Last Glacial Maximum and present-day conditions Ellerhoff et al. (2022). These contrasting results may be due to different biases in feedbacks or missing physics that could be responsible for a potential real-world state dependency. Detailed and direct observations are needed to complement modeling results. But even the largest eruptions of the satellite era are not large compared to eruptions that eventually occur over time spans of a hundred years or more. Thus, the impact of such eruptions needs to be reconstructed by paleoclimate proxy records that go beyond the observational period. It is challenging to obtain such records with sufficient temporal resolution and accurate dating. Ice cores arguably provide the most detailed records covering timescales of years up to several hundred millennia. This is because the temporal resolution of the material is large compared to other common stratigraphic archives, which often allows for a layer-counted timescale.

The ejection of sulfate aerosols into the stratosphere by large volcanic eruptions leads to a sharp peak in polar ice-core sulfate records with a delay of roughly 1–2 years (Burke et al.2019). Based on the integrated sulfate concentration in Greenland and Antarctic ice cores, continuous records of volcanic eruptions along with rough estimates of the magnitude of the eruptions can be constructed (Zielinski et al.1997; Castellano et al.2004; Gao et al.2007; Sigl et al.2015, 2022). Here we use two recently compiled datasets: first, a record of volcanic eruptions in the period 12–60 ka with sulfate peaks detected simultaneously in Greenland and Antarctica (Svensson et al.2020), and second, continuous records of volcanic eruptions detected in either Greenland or Antarctic ice cores (Lin et al.2022). The former represents significant volcanic eruptions, which most likely distributed sulfate aerosols globally in the stratosphere and are thus expected to have a global climatic impact. The latter is a much larger set that also includes eruptions with more regional aerosol distribution.

By analyzing eruptions during the time interval 12–60 ka and comparing them to large historic eruptions, we provide a first attempt of using ice-core data to quantify the cooling effect of very large eruptions with return periods of hundreds of years and more. To this end, sulfate-derived records of volcanic eruptions are combined with high-resolution δ18O records from the same ice cores. δ18O is a widely used proxy for surface temperature at the accumulation site, which can be measured with up to sub-annual time resolution. The variability at such short timescales may not represent reliable climatic information, however, because the original temperature signal is altered by post-depositional processes (Münch et al.2016). These result in high-frequency noise, referred to as stratigraphic or glaciological noise, as well as a smoothing of short-term anomalies. It is unknown how much climatic information remains at sub-decadal timescales in the glacial ice-core record (Vinther et al.2010). Here we compare the average short-term cooling signal of a large number of volcanic eruptions to the non-volcanic proxy variability. This provides insights into the signal preservation of the δ18O proxy that are useful for future studies on ice-core data of increasingly high resolution. There are large quantitative uncertainties in the calibration of the δ18O temperature proxy in the glacial period, which make it difficult to estimate the volcanic cooling in absolute terms. Thus, we complement our analysis with direct ice-core observations of changes in (annual) snow–water accumulation following the detected eruptions, which are not subjected to an unknown calibration. Snow accumulation is known as a climate-sensitive parameter on large ice sheets, and reductions in precipitation are expected after large volcanic eruptions (Robock and Liu1994; Bala et al.2008).

The glacial volcanic record also allows us to assess a potential state dependency of the climate response, since it features the so-called Dansgaard–Oeschger (DO) cycles. These are abrupt regime shifts between quasi-stable colder and milder Northern Hemisphere climate conditions, known as Greenland stadials (GSs) and Greenland interstadials (GIs). The glacial climate resides in these quasi-stable states for centuries up to several millennia. Using different subsets of eruptions, we investigate how the volcanic δ18O anomaly depends on the climate background state, as well as the sulfate deposition magnitude of the eruptions. This yields observational evidence that complements ongoing investigations into the state dependency of climate sensitivity (Caballero and Huber2013; Köhler et al.2015; von der Heydt et al.2016; Ashwin and von der Heydt2020).

2 Methods and materials

2.1 Records of volcanism

We investigate two records of volcanic eruptions. First, we study the 82 volcanic eruptions identified simultaneously in Greenland and Antarctic ice cores by Svensson et al. (2020) in the period 12–60 ka. These are referred to as bipolar eruptions hereafter. Due to difficulties in matching Greenland and Antarctic ice cores around the time of the Last Glacial Maximum, this dataset has a gap from 16.5–24.5 ka. Further, in Svensson et al. (2020) not all eruptions could be identified with a sulfate spike in all ice cores under consideration. It is unknown whether in these cases the eruption did not yield any sulfate deposition at the ice-core site, whether the sulfate deposition was wiped away by snow redistribution, or whether the missing eruption is due to limitations in data resolution and the synchronization procedure. It is likely that highly localized phenomena strongly influence the amount of sulfate deposition that is preserved. A previous study showed that even large events such as the 1815 CE Tambora eruption can be entirely missing in several of a handful of very nearby replicate cores (Gautier et al.2016). Since for our study a precise alignment of the δ18O records to the sulfate spikes is crucial, only those ice cores enter our analysis where a given eruption has been identified.

The second dataset is a record of volcanic sulfate depositions in either Greenland or Antarctic ice cores in the period 9–60 ka compiled by Lin et al. (2022), which we restrict to the glacial period 11.7–60 ka. This dataset consists of the depth of several hundred eruptions in the NGRIP (N=780), NEEM (N=311), GISP2 (N=282), EDC (N=211), WAIS (N=470), and EDML (N=470) ice cores, along with estimated magnitudes derived from the integrated sulfate deposition in the respective cores. Due to differences in resolution and quality of the underlying sulfate records (see Sects. 2.2 and S1), some ice cores allow for the identification of a larger number of eruptions with smaller average sulfate peaks compared to other cores (Fig. S3). Thus, for instance, the NGRIP dataset contains many more small eruptions compared to NEEM.

There are also large differences in the estimated sulfate deposition of different cores for the same eruptions. Reasons for this include actual differences in deposition quantity due to the different locations on the ice sheet, local relative differences in wet and dry deposition, differences in the sulfate measurement method and resolution, post-depositional snow redistribution, and potential biases in the thinning function used for the different ice cores. Lin et al. (2022) give one composite volcanic record for Greenland and one for Antarctica with 1019 and 691 eruptions, respectively. Therein, the Greenland Ice Core Chronology 2005 (GICC05) age has been derived from the sulfate spikes in one or more ice cores, and a large subset of eruptions has been matched within cores of the same hemisphere. In this dataset, the sulfate deposition estimate is given individually for all cores where the eruption has been identified and as an average over the cores.

In the subset of eruptions that were matched in two or more cores of one hemisphere the scatter of deposition values for individual eruptions in different cores is large (Fig. S4), and the mean deposition values can differ significantly. For instance, on average, the same eruptions in NEEM and GISP2 have a larger estimated sulfate deposition compared to NGRIP. As a result, when cores with larger sulfate deposition only occasionally contribute to the calculation of the average deposition values given by Lin et al. (2022), a certain degree of statistical noise is introduced. But since the differences in the mean deposition of cores are much smaller than the scatter of deposition values among the same eruptions in different cores (Fig. S4), we use the average deposition in our analysis, unless noted otherwise. This is supported by the abovementioned observation that even large eruptions can be entirely missing in individual cores, which underlines the fact that the deposition values of individual cores are often not reliable.

Most of the eruptions in the dataset of Lin et al. (2022) are not matched across Greenland and Antarctica. But the dataset does include the bipolar eruptions previously identified by Svensson et al. (2020). Importantly, this dataset will be referred to as unipolar hereafter, even though the eruptions from Svensson et al. (2020) are still included.

2.2 Fine tuning and calibration of the eruption ages

The depths of the eruptions are not known with arbitrary precision, especially in ice cores where the underlying sulfur datasets are of low resolution and/or very noisy. Here we use the nominal depths reported in Lin et al. (2022) when investigating the unipolar dataset and the nominal depths from Svensson et al. (2020) when analyzing the bipolar eruptions. These depths are then transferred to the common age scale (see next section), followed by a slight recalibration of the eruption ages, as explained in the following. First, there are slight systematic average offsets of the nominal depths compared to the sulfate maxima. This is a result of the detection of individual eruptions from noisy data combined with a slight asymmetry of the sulfate peaks, as well as the usage of multiple proxies in Svensson et al. (2020). In Fig. S1, the average sulfate peaks over bipolar and unipolar eruptions in all cores are shown, and one can see slight offsets of up to 2 years with respect to the nominal ages. Here we choose to correct these offsets and shift the eruption ages such that in each core the sulfate peaks on the age scales are aligned on average (see Sect. S1 for more details).

Second, we further shift the ages slightly by a fixed amount to account for the fact that the maximum sulfate peak in the ice core is delayed with respect to the eruption age. For large historic eruptions comparable in size to the bipolar eruptions investigated here, this delay is estimated to be around 1.5 years (Burke et al.2019). We shift all eruption ages back in time by 1.5 years relative to the time of maximum sulfate deposition. Ideally, one would do this individually for each eruption by determining the start depth of the sulfate peak as an estimate of the actual starting year of the eruption. However, an individual age adjustment would only increase the jitter along the time axis, since the sulfate records are noisy due to intermittent deposition and snow redistribution and since the peaks of volcanic origin are subjected to smoothing by diffusion and different measurement techniques and resolution, which leads to peak widths that vary greatly across cores and time periods (Figs. S1 and S2). Thus, when interpreting our reported δ18O anomalies averaged over different eruptions, it should be kept in mind that the events are aligned using the maximum sulfur deposition shifted by 1.5 years toward older ages and not using the unknown, true time of the eruption start. In the plots where we report the time before eruption along the horizontal axis, the year 0 indicates our estimate of the starting time of the eruptions as described here.

Table 1Median time resolution of δ18O records (in years). The WAIS data were measured by a different technique, leading to a very high sample resolution. The true time resolution of the δ18O signal is lower (see main text).

Download Print Version | Download XLSX

2.3 High-resolution oxygen isotopes

To quantify the climatic impact of the eruptions, we use high-resolution δ18O records from four Greenland ice cores (NGRIP – NGRIP Members2004; Gkinis et al.2014; GRIP – Johnsen et al.1997; GISP2 – Stuiver and Grootes2000; NEEM – Rasmussen et al.2013) on the annual layer-counted GICC05 timescale (Svensson et al.2006, 2008; Rasmussen et al.2013; Seierstad et al.2014), as well as from three Antarctic ice cores (WAIS – Buizert et al.2015; Jones et al.2018; EDC – Jouzel et al.2007; EDML – EPICA Community Members2006) that have been matched to GICC05 at the bipolar volcanic eruptions (Svensson et al.2020). All records cover the period from 11 700 years b2k (years before 2000 AD) to 60 000 years b2k.

Since the records were measured at different depth resolutions and were taken at sites with different accumulation and thinning rates, their time resolution varies (see Table 1). The WAIS record was measured with continuous-flow analysis (CFA), yielding a dataset with high depth resolution of 5 mm, which results in the sub-annual time resolution given in Table 1. Note, however, that the effective measurement resolution is slightly lower due to mixing of material within the CFA apparatus (Jones et al.2017). Perhaps apart from the low-resolution EDML dataset, in all records the true δ18O time resolution, i.e., the degree of preservation of the temporal isotopic variability that was originally deposited on the ice sheet, is lower than given in Table 1. This is because of diffusion of the water molecules in firn and ice, which highly attenuates any variability below multi-annual timescales for ice of the last glacial period, also in the high-resolution WAIS record (Cuffey and Steig1998; Jones et al.2017, 2018).

For all time series used here, each data point consists of the average δ18O value of bulk material in contiguous depth intervals. The data are thus not point samples but averages over contiguous intervals. For our study the δ18O records are processed in the following way. The midpoints of the depth intervals are interpolated linearly to the GICC05 time–depth scale, yielding an unequally spaced time series. Next, this series is oversampled to a 1-year equidistant grid using nearest-neighbor interpolation. An exception is the sub-annual WAIS record, which is first averaged to 1-year resolution. Like this, the nature of the measurements as contiguous depth averages and the original measurement values are preserved, and all records are placed on the same equidistant time grid. We furthermore construct a stacked Greenland record in time slices around the bipolar volcanic eruptions. For a given eruption all individual cores where a depth has been recorded in Svensson et al. (2020) are centered around the eruption depth and averaged.

For comparison with known historic eruptions we consider high-resolution Holocene δ18O records from four different Greenland ice cores (NGRIP, GRIP, GISP2, and Dye-3; Vinther et al.2006), covering the last 2000 years. The time resolution in this period varies from monthly (Dye-3, GRIP) to biennial (GISP2). All four records have annual or higher resolution from 0–1.2 ka, and for the period >1.2 ka this is still the case for all cores except GISP2. The measured data on the depth scale are processed by interpolating the midpoints of the depth intervals linearly to the GICC05 time–depth scale, yielding an unequally spaced time series. Then, this series is oversampled to a monthly equidistant grid using nearest-neighbor interpolation. Only the Dye-3 record features a seasonal cycle, which is removed by a running yearly average. Subsequently, the records are stacked and a time series without trends and centennial or millennial variability is obtained via high-pass-filtering the record by removing a Gaussian kernel smoother with 150-year standard deviation.

Figure 1(a) Average NGRIP δ18O anomaly centered at the bipolar eruptions, defined with respect to the mean of the period 10–50 years prior to the eruption. The average signal is shown in blue, and the gray bands are the 16th to 84th percentiles of detrended time slices covering individual eruptions. In orange (green) we show the mean signal of eruptions during GI (GS). (b) Same for the Greenland stack, where detrended slices of all cores for every eruption are averaged, using only cores where a depth is identified – 49 eruptions with four cores, 14 (10) with three (two) cores, and nine represented by NGRIP only. (c) Distribution of 6-year average anomalies of the Greenland δ18O stack around the eruptions (blue) compared to a bootstrap of randomly chosen 6-year anomalies from the stack using all four cores on the GICC05 synchronization (gray). The dashed black line is the 16th percentile of the bootstrap distribution, and the blue line is the mean of the volcanic anomalies. Dashed orange (green) lines show individual eruptions during GI (GS). Red lines are eruptions preceding the onsets of Dansgaard–Oeschger events within less than 50 years, as identified in Lohmann and Svensson (2022). (d) Anomalies of the Holocene Greenland stack (see “Methods and materials”). Shown is the cooling amplitude of several major Common Era eruptions, as well as the bootstrap distribution of random segments from the past 2 kyr. The historic eruptions are 1815 CE Tambora, 1258 CE Samalas, and 43 BCE Okmok, as well as the 536–540 CE doublet. For the latter we chose the age of 536. In most δ18O records the doublet is merged due to diffusion. GICC05 ages are taken from McConnell et al. (2020) and Sigl et al. (2015).


2.4 Records of layer thickness

The NGRIP and WAIS cores have been layer-counted up to a certain depth. Subsequent depths of counted layers comprise an annual record of the layer thickness, which we use to study post-eruptive changes in accumulation rate. In NGRIP the layer counting was performed until 60.2 ka BP, and thus the resulting record of annual layer thickness (Rasmussen et al.2023) covers the entire investigated period. The counting includes certain and uncertain layers. For the certain layers, the depth increment corresponds to a 1-year time increment. In uncertain layers, which make up 10.1 % of all layers, subsequent depths are defined as a half-year time increment (Andersen et al.2006). To obtain the layer thickness record, we first convert the depth–age pairs of the GICC05 chronology to thickness–age pairs by taking the increments of subsequent depths. Then, to homogenize the record of full and half-years, we linearly interpolate the record to a 0.1-year grid. The WAIS core was layer-counted until 31.2 ka BP (Sigl et al.2016), thus only covering the younger part of the glacial. Here we use the layer-counted WD2014 chronology, which does not include half-years for uncertain layers. Otherwise, it is processed in the same way as for NGRIP.

3 Results

3.1 Volcanic cooling observed in Greenland after bipolar eruptions

We first consider the Greenland δ18O signal following the bipolar eruptions. For each eruption, the δ18O records are centered around the estimated time of eruption, and a segment of 50 years before and after the eruption is chosen and detrended linearly. To obtain anomalies we subtract the mean value of the detrended signal in the interval 10 to 50 years before the eruption. Finally, we extract the mean cooling anomaly from the non-volcanic variability by averaging these 100-year anomaly time slices over all eruptions. In Fig. 1a, the results are shown for the NGRIP core, which has the best temporal resolution. A negative multi-annual anomaly is seen, which clearly exceeds the variability in the mean signal leading up to the eruption. However, the mean anomaly is only approximately half the size of the high-frequency isotope variability around individual eruptions (gray bands). The other Greenland ice cores show the same qualitative behavior, but the signals are less sharp due to the lower resolution (Fig. S5).

We attempt to remove non-climatic noise by averaging across all Greenland cores, as shown in Fig. 1b. Here the average isotopic cooling anomaly begins significantly prior to the estimated eruption age. This is due to diffusion of water molecules in firn and ice, as well as the averaging introduced by the isotope measurement on bulk material at multi-annual to decadal resolution. In addition, since the eruptions are aligned at the sulfate maxima and a constant 1.5-year shift towards older ages was applied to estimate the true eruption ages (see Sect. 2.2), for many eruptions we can expect the true age to be older than our estimate. This especially holds for larger eruptions with longer-lasting sulfate deposition, as well as for records with poor resolution and wider sulfate peaks (Fig. S1).

To quantify the isotopic cooling, we use the average signal (Fig. 1b) to define a time period corresponding to the most pronounced anomaly. This period consists of the estimated year of the eruption as well as the following 5 years (green shading in Fig. 1b). The average anomaly over this time period gives a scalar estimate of the volcanic cooling for each eruption, which we call the cooling amplitude hereafter. There is a rather weak correlation of this scalar cooling estimate of the individual eruptions among the Greenland cores (Fig. S6), suggesting large non-climatic noise in the high-resolution records. From the distribution of amplitudes in individual eruptions (Fig. 1c) it is clear that a non-negligible number of eruptions are followed by a positive δ18O anomaly, i.e., a potential warming associated with the eruption. For the Greenland stack, 23 out of the 82 bipolar eruptions feature a positive δ18O anomaly. It is unclear whether these eruptions indeed induced no volcanic cooling in Greenland or whether it is masked by positive anomalies in the non-climatic noise and multi-annual climate variability. We construct a bootstrap distribution of 6-year anomalies of randomly chosen segments from the Greenland δ18O stack for the entire 12–60 ka period (gray distribution in Fig. 1c), which shows that natural fluctuations in δ18O anomalies of up to 1 ‰ are common. Using present-day calibrations of the δ18O thermometer of 0.69 ‰ K−1 to 0.8 ‰ K−1 (Sjolte et al.2011; Buizert et al.2014), this would correspond to 6-year mean temperature anomalies of up to 1.25–1.45 K, which would have the potential to mask the volcanic cooling of even very significant eruptions. Due to this large uncertainty, one cannot interpret the amplitude of individual eruptions as a quantitative estimate of the volcanic cooling. Nevertheless, the distribution of amplitudes is clearly shifted towards negative values, unlike the bootstrap distribution, which is symmetric and centered at 0.

We define a signal-to-noise ratio (SNR) of the volcanic signal by dividing the mean volcanic anomaly (blue line in Fig. 1c) by the 16th percentile of the bootstrap distribution (as a measure of standard deviation, dashed black line in Fig. 1c). This is not an SNR of the record as a temperature proxy in absolute terms. We consider only the volcanic cooling anomaly to be a signal, while internal temperature variability not caused by volcanism is considered noise, alongside actual proxy noise from intermittent precipitation and post-depositional processes. Thus, what we define as SNR measures the strength of the multi-annual volcanic cooling signal relative to the multi-annual climatic and non-climatic proxy variability. For the Greenland stack this yields SNR=0.66, and for NGRIP we find SNR=0.48. Thus, stacking improves the SNR, but the average anomaly still does not exceed the non-volcanic variability. While noise in the vertical axis of Fig. 1b is reduced when stacking different cores, additional noise is introduced in the horizontal axis since not all cores have an equally good alignment of the isotope record relative to the true eruption age. This is because the precise eruption depth is less certain in some cores due to low resolution of the underlying sulfate records (Figs. S1 and S2). Further, there are small systematic offsets in the depth scale of δ18O and sulfate measurements of the same ice core, as they are not obtained from the same samples.

The average 6-year cooling amplitude is 0.48 ‰ in the stack and 0.63 ‰ in NGRIP. This may be compared to the largest eruptions in the Common Era. These are much better constrained since most of them have an identified source, a well-quantified magnitude, and a precise date, allowing them to be matched to other paleoclimate proxies, such as from tree rings (Sigl et al.2015). In Fig. 1d we show a distribution of anomalies from randomly chosen segments of a Greenland δ18O stack covering the last 2000 years (“Methods and materials”), together with the cooling amplitude of four major historic eruptions that have been estimated to be among the five largest eruptions during this time interval (Sigl et al.2015). These feature an average negative δ18O anomaly of 0.73 ‰. The average isotopic anomaly of the bipolar eruptions during the glacial is thus slightly weaker by comparison. However, the calculated glacial δ18O anomalies are likely underestimated compared to the Common Era eruptions due to several factors discussed in Sect. 3.3.

Figure 2Same as Fig. 1b–c, but for the Antarctic cores WAIS (a, c) and EDC (b, d). For EDC, only 78 out of 82 bipolar eruptions were detected in Svensson et al. (2020).


3.2 Volcanic cooling observed in Antarctica after bipolar eruptions

The average volcanic isotopic anomaly in Antarctica is more subdued, which may be expected as Antarctica is climatically relatively isolated and more volcanos are located in the Northern Hemisphere. The WAIS record is a priori best suited to show a clear volcanic cooling signal due to its high accumulation rate and measurement resolution. A roughly 4-year-long average negative δ18O anomaly is found (Fig. 2a), but it is only marginally significant as it is not much larger than the variations in the mean anomaly before and after the eruption. The average δ18O cooling anomaly in EDC is much less sharp (Fig. 2b) because of the low accumulation rate. This yields an average δ18O resolution of almost 10 years, along with pronounced diffusion and non-climatic noise (Münch et al.2016). The EDML core does not show any cooling signal for the bipolar dataset (Fig. S5d), partly because its isotopic resolution is too low. In addition, δ18O records close to coastal regions have been found to capture only very little local temperature variability on short timescales (Vega et al.2016; Goursaud et al.2019). Nevertheless, by averaging over many eruptions from the unipolar dataset a slight cooling anomaly can be discerned (not shown here).

Figure 3(a, b) Average δ18O anomalies in the Greenland stack (a) and the WAIS record (b) aligned to the bipolar eruptions. The gray shading and blue curve are the same as in Figs. 1b and 2a. The red (yellow) curves correspond to the average δ18O anomaly of the younger (older) half of the bipolar dataset. (c) Average δ18O anomalies in the NGRIP record aligned to the unipolar eruptions. The average signal is shown in black, and the gray bands are the 16th to 84th percentiles of detrended time slices covering individual eruptions. The colored curves are the average signals of four equally sized subsets of eruptions divided according to age.


Figure 4(a) Average integrated δ18O anomaly in NGRIP for the unipolar dataset, where the eruptions are separated according to age in consecutive bins with an equal number of eruptions (10 bins with 78 eruptions each). The δ18O time slice around each eruption in a bin is detrended as described in the main text. The detrended slices of eruptions in a bin are then averaged to yield the mean δ18O anomaly, from which the integrated anomaly is defined as the area above the consistently negative signal around the time of the eruption (yellow shading in the inset). The red square corresponds to the value for the youngest bin of eruptions, occurring from 15 685 to 12 022 years BP, and it is a thus far unexplained deviation from the roughly linear decrease in the anomaly over time, as indicated by the linear fit. (b) Signal-to-noise ratio in the NGRIP record estimated in a 12 kyr moving window. The method, explained in Sect. 3.1, is applied here to the unipolar dataset using 4-year average anomalies starting with the year of the eruptions. Shown are curves for all eruptions, as well as for the GI and GS subsets. The GI curve is interrupted from 20–28 ka, as there are too few eruptions (fewer than 20 per 12 kyr) for a robust SNR estimation.


We again define a period of most pronounced cooling based on the average anomaly curves. For WAIS this corresponds to the estimated eruption year, as well as the year before and the 2 after. Figure 2c shows that the cooling amplitudes associated with bipolar eruptions, which average −0.20 ‰, are only shifted slightly towards negative values compared to randomly selected periods of the record. For EDC we choose an almost symmetric period with 7 years before and 6 years after the estimated eruption year. This also yields an average anomaly of −0.20 ‰ and a slight negative shift of the distribution of individual anomalies (Fig. 2d). Since the EDC record has a much lower sample resolution and thus more pronounced smoothing due to averaging, the original peak anomaly would be clearly larger in absolute terms compared to WAIS. Still, compared to the proxy background variability, the average volcanic signal is similar for the two cores. The SNR derived from the distributions in Fig. 2c and d is SNR=0.30 for WAIS and SNR=0.28 for EDC. These low values highlight the fact that the volcanic cooling signals are weaker in the Antarctic cores compared to Greenland, which may be in part due to a more muted or variable Antarctic climate response, but also due to poorer performance of the δ18O proxy. Indeed, cooling amplitudes of individual eruptions in WAIS and EDC are not significantly correlated, and the amplitudes of both Antarctic cores are also not correlated with the amplitudes from the Greenland stack (Fig. S7).

3.3 Preservation of the cooling signal in the isotope record

The above estimates of the average multi-annual isotopic cooling anomaly are a compound of young eruptions and older ones, for which the signal is degraded due to several effects. First, multi-annual δ18O anomalies are smoothed out by diffusion of water molecules in the ice. The older the ice the more time has elapsed for the diffusion to act. Additionally, deeper annual layers become thinner due to ice flow, which leads to increasing ice diffusion length (in years) with depth (and thus age). Second, there is additional smoothing due to the measurement of δ18O on contiguous pieces of the ice core at constant depth intervals. For thinning annual layers with age, this smoothing by averaging is more pronounced the older the eruption. Third, the effective temporal resolution of the underlying sulfate records typically decreases with age (Table S1). As a result the eruption age can be determined less accurately for older eruptions (Fig. S2a), which again leads to a smearing-out of the average cooling anomaly.

Figure 5(a, b) Average δ18O anomalies in the Greenland stack (a) and the EDC record (b) aligned to the bipolar eruptions. The black curve is the average signal over all bipolar eruptions, and the red (blue) curves correspond to the average δ18O anomaly of half the bipolar dataset with a larger (smaller) total sulfur aerosol loading (Lin et al.2022) compared to the median (127 Tg). (c) Mean δ18O anomalies (per mill) in the WAIS core in response to the unipolar eruptions. Black bars show the mean signal, and the gray bands are the 16th to 84th percentiles of detrended time slices covering individual eruptions. The red curve corresponds to the mean anomaly in the largest 10 % of eruptions in terms of their WAIS sulfate deposition magnitude.


Figure 6(a) Correlation of the integrated NGRIP δ18O anomalies (eruptions from the unipolar dataset) and the associated sulfur deposition in the NGRIP core (Lin et al.2022). Individual dots represent the average integrated anomaly of the eruptions divided into bins according to the 15th, 30th, 45th, 60th, 75th, and 90th percentiles of the NGRIP sulfur deposition. The integrated anomaly is defined by the sum of averaged δ18O anomaly values in the years around the estimated year of eruption where the anomaly is negative, as shown with the shaded area in the inset. We also show a linear regression with 95 % confidence intervals, which has a nonzero intercept of −4.3 ‰. (b) Average NGRIP δ18O anomaly as a function of the baseline δ18O level, defined by the mean δ18O value from 50 up until 3 years before the unipolar eruptions. The data are averaged in equally sized bins according to the δ18O baseline values. Shown is the minimum of the average δ18O anomaly (black), as well as the integrated average δ18O anomaly, as was used in panel (a) (dashed green).


Consequently, while the average magnitude of the eruptions measured by their sulfate deposition does not appear to change over the course of the glacial (see Figs. S8 and S9, as well as Lin et al.2022), younger eruptions show a more pronounced cooling anomaly compared to older ones (Figs. 3a, b, and S10 for all other cores). In the Greenland stack, the younger half of eruptions show a minimum anomaly of −0.75 ‰ in the year after the eruption. With the abovementioned present-day calibrations of the δ18O thermometer of 0.69 ‰ K−1 to 0.8 ‰ K−1 (Sjolte et al.2011; Buizert et al.2014), this yields a peak cooling of 0.94–1.09 K, which comes close to the 1.24 K summer NH cooling estimated from tree rings for the four largest eruptions of the Common Era (Sigl et al.2015).

The evolution over time of the δ18O cooling anomaly can be investigated more precisely using the NGRIP core in isolation, which has the highest δ18O and sulfate resolution, as well as the best dating. Here, in the younger half of the bipolar eruptions the 6-year mean isotope amplitude is −0.77 ‰ and the peak cooling amplitude is −0.90 ‰. Using the unipolar eruption record, which features many more eruptions and thus a less noisy mean signal, we find that eruptions occurring in the period 12–32 ka feature a minimum anomaly of −1.7 ‰ 2 years after the estimated eruption age (Fig. 3c). Thus, despite a return period of only 65 years, the cooling anomaly of the youngest glacial eruptions clearly exceeds the anomaly after the largest eruptions in the Common Era. For the older eruptions, the minimum anomaly is attenuated by roughly a factor of 2.

The amplitude of the isotopic cooling signal at the time of deposition is expected to be even larger because (a) the δ18O records are not perfectly aligned to the true eruption year, and (b) the smoothing effect of diffusion has not been accounted for in our study. There are techniques to achieve the latter if the diffusion length in ice and firn is known (Johnsen et al.2000). Here we refrain from doing so because the variations over time of the cooling anomalies do not appear to follow a simple diffusion process. While the peak cooling anomalies in Fig. 3c decrease over time, the anomaly does not get visibly smeared out further in time. The area under the curve corresponding to the negative anomalies does not stay constant, as expected for a simple diffusion of temperature fluctuations over time, but decreases over time (Fig. 4a). Further, in contrast to a constant SNR due to a roughly equal diffusive attenuation of the noise background and volcanic signal, the SNR decreases over time (Fig. 4b). This may reflect the additional attenuation effect on the volcanic signal of the decreasing precision of the eruption alignment going further back in time.

Figure 7Mean δ18O anomalies (per mill) in response to the bipolar eruptions in the Greenland stack (a) and EDC (b) records. The black curve and gray shading are as in Figs. 1b and 2b, and the colored curves are the average anomalies of the subsets of eruptions classified by Lin et al. (2022) as Northern Hemisphere (NH, above 40 N) and low-latitude Southern Hemisphere (LL/SH, below 40 N) based on the relative sulfur deposition measured in Greenland versus Antarctica.


3.4 Correlation of cooling signal with volcanic magnitude and hemispheric sulfur deposition

We next investigate whether eruptions that were large in terms of their sulfate deposition also led to a more pronounced isotopic cooling. Considering the total sulfur aerosol loading, which is a combined metric of the Greenland and Antarctic sulfate deposition (Lin et al.2022), the bipolar dataset is divided into two, depending on whether an eruption was larger or smaller compared to the median. The larger eruptions feature a much more pronounced δ18O anomaly, as shown for the Greenland stack and EDC in Fig. 5a and b. For individual eruptions there is only a weak correlation of δ18O anomaly and aerosol loading (Fig. S11). This is because of the high noise levels in the records, which, together with the relatively small size of the bipolar sample, limit our ability to quantify the dependence of the δ18O anomaly on the sulfate deposition. For this purpose the larger unipolar dataset is used in the following, with the caveat that only the unipolar deposition is available and not the total aerosol loading.

Also, in the unipolar dataset eruptions with a larger sulfate deposition show a more pronounced average δ18O anomaly. The WAIS core, for instance, which showed only a weak average anomaly after bipolar eruptions, features a much stronger cooling signal for the eruptions with largest unipolar sulfate deposition (Fig. 5c). We next define a scalar metric for the cooling associated with eruptions of a certain size category. To this end we separate the eruptions from the unipolar dataset of a given core into bins according to their unipolar sulfate deposition. Here we consider the two cases: (a) the deposition values are only taken from one core, or (b) the deposition values are averaged over all cores where a given eruption was identified (see Sect. 2.1). We then calculate the average δ18O anomalies of all eruptions in one bin and integrate the resulting mean signal over the years around the eruption where a negative anomaly is found (yellow area in inset of Fig. 6a, see also Fig. S12). For all cores, this integrated isotopic anomaly shows a clear relation with the unipolar deposition magnitude (Figs. 6a and S13). A linear fit to the data seems justified in most cases, but we cannot rule out a nonlinear relation. For most cores the linear fit indicates that there are still significant negative isotopic cooling anomalies when extrapolating to zero sulfate deposition. We speculate this could be because (a) the linear relationship breaks down for the smallest eruptions that still have a global cooling effect but no polar sulfate deposition or because (b) the largest sulfate deposition values are inflated due to a significant proportion of local or regional eruptions with a large tropospheric sulfate transport and polar deposition. There is generally a better correlation of the anomaly with the deposition in the individual core and not the deposition averaged over multiple cores (Fig. S13). This may be surprising since the latter should be a more reliable estimate for the magnitude (Sect. 2.1). A reason for this may be that the eruptions with a large depositional sulfate peak in the respective core feature a more precise depth estimate, leading to a better average alignment of the isotopic cooling anomaly to the true age of the eruption.

Based on the relative deposition of bipolar eruptions in Greenland and Antarctica, the source latitudes have been classified in binary categories with Northern Hemisphere (NH above 40 N) or Southern Hemisphere and low-latitude (SH/LL) eruptions (Lin et al.2022). Since there is a correlation of the isotopic anomaly with the unipolar deposition magnitude in all cores, we see an according stronger Greenland (Antarctic) isotopic response for NH (SH/LL) eruptions (Fig. 7). For eruptions with a larger Greenland sulfate deposition (classified as NH) there is no significant EDC δ18O cooling anomaly. It may be that a non-negligible number of these eruptions even feature a warming anomaly, which might be reflected in the positive δ18O excursion in the confidence bands for the lower-resolution Antarctic cores (Figs. 2b, S5d and 7b) and in the slight indication of a bimodal distribution in Fig. 2d. Note, however, that a certain widening of the confidence bands is expected in the low-resolution records due to the detrending and nudging of the anomaly to the period prior to the eruption. Moreover, since there are relatively more NH classified eruptions in GS compared to GI, we cannot clearly separate the effect of the estimated eruption latitude on the δ18O anomaly from the even more pronounced GI–GS contrast (see next section). Larger bipolar datasets would be required to resolve this, and at this stage we believe that neither the determination of the eruption latitude nor the inferred volcanic cooling from the δ18O proxy is precise enough to warrant much speculation on the dependence of the climate response as a function of the eruption site.

Figure 8Average δ18O anomaly in the NGRIP (a–c) and EDC (d–f) core, obtained by aligning the records at the volcanic eruptions from the (a, b) bipolar and (c, d) unipolar datasets, as well as a subset of the unipolar dataset representing a time period with an equal proportion of GS and GI conditions (e, f). The eruptions are separated into subsets occurring during GI (GS), and the average δ18O anomaly is shown in black (red). The dashed blue lines show the anomaly curves obtained when resampling the GS subset in the NGRIP and the GI subset in EDC such that the distribution of the associated magnitudes of sulfur deposition matches the distribution of the eruptions in the corresponding GI subset for the NGRIP and the GS subset for EDC (see main text for more detail).


3.5 State dependency of the climate response

In Sect. 3.3 we found that the younger, best-preserved glacial eruptions in NGRIP feature a significantly stronger isotopic cooling compared to the largest eruptions during the Common Era. This indicates a state dependency of the volcanic δ18O anomaly, which could reflect a state dependency of climate sensitivity, i.e., of the response of the global average surface temperature to the radiative forcing of the volcanic eruption. Instead of global climate sensitivity, there could also be a state dependency in regional climate sensitivity, i.e., a difference in the spatial pattern of the temperature change. Here we use the term climate sensitivity to describe the temporary temperature change following the relatively short-lived volcanic perturbation. This differs from the concept of equilibrium climate sensitivity, which describes the long-term increase in global temperature as a response to an instantaneous, permanent doubling of atmospheric CO2. The response to short-term volcanic forcing involves fast components of the climate sensitivity, but longer-term climate feedbacks and changes in deep-ocean heat content are limited. Accordingly, modeling studies have found difficulties in using observations of volcanic cooling to constrain equilibrium climate sensitivity (Boer et al.2007; Bender et al.2010; Merlis et al.2014).

It is also possible that the temperature response is essentially independent of the background climate state and that instead there is a state dependency of the proxy sensitivity, i.e., of the magnitude of δ18O change for a given volcanic temperature change. Compared to the present-day proxy sensitivity, previous work suggests that the δ18O proxy in Greenland reacts more sensitively to the temperature changes across glacial regime shifts, such as the last deglaciation and DO events (Guillevic et al.2013; Buizert et al.2014), while the opposite is the case for Antarctic δ18O (Uemara et al.2012; Buizert et al.2021). This is due to a combination of changes in accumulation seasonality, moisture source, and ice sheet topography associated with the regime shifts. But the sensitivity of the proxy to short-term temperature changes without major regime shifts is unknown, and its dependence on the background climate state remains an active subject of research (Liu et al.2023; Cauquoin et al.2023). Comparing the mean volcanic δ18O anomaly to the baseline δ18O values at which the corresponding eruptions occurred, we find a nonlinear dependence of the anomaly on the background state (Fig. 6b). This could be interpreted as a state dependency of the proxy or climate response, but it partly reflects the better signal preservation for the predominantly young eruptions occurring at low δ18O baseline values as a result of the gradual decrease in δ18O values throughout the glacial.

Figure 9(a) Frequency of eruptions in the Greenland and Antarctic unipolar datasets in different magnitude categories (defined as the logarithm of the unipolar sulfate deposition in kg km−2), which was derived as averages over the deposition in the individual cores where the eruptions could be detected (Lin et al.2022). The datasets are further split up into eruptions occurring during GS and GI. The error bars on the frequency estimates are given as black lines and represent the 10th to 90th percentile computed analytically assuming a Poisson distribution for the number of eruptions occurring in the respective time intervals. (b–e) Distributions of the magnitude of the eruptions, given by the logarithm of the unipolar sulfate deposition, for the Greenland (b–c) and Antarctic (d–e) datasets, which are divided into eruptions occurring during GS and GI.


A more conclusive picture of the state dependency for both Greenland and Antarctica can be obtained by dividing the datasets into eruptions occurring during the cold (GS) and mild (GI) periods of DO cycles. While changes in Antarctic climate over DO cycles are much weaker compared to Greenland, the DO cycles are nevertheless the most pronounced large-scale climate regime shifts of the last glacial. Thus, dividing the data into GI and GS periods indicates which part of the DO cycle the global climate state occupies. This seems a reasonable target to test the climate state dependency for both Greenland and Antarctica. Eruptions occurring during GS show a more pronounced isotope anomaly in Greenland compared to eruptions during GI, while the opposite is the case for the Antarctic EDC core (colored lines in Figs. 1a, b and 2b). The response pattern in WAIS seems similar to Greenland (Fig. 2a), but it is inconclusive since the signals are not larger than the variability before the eruptions. Figure 8a and b show the mean anomaly signals in NGRIP and EDC in more detail. The stronger Greenland GS response is surprising because we would a priori expect a sharper volcanic cooling response in GI due to the higher accumulation rate resulting in a higher resolution and less pronounced non-climatic noise. Further, the higher accumulation rate also leads to a higher-resolution sulfate record and thus a sharper estimate of the eruption depth.

Figure 10(a, b) Average layer thickness change after the volcanic eruptions of the unipolar dataset during GI and GS in the NGRIP (a) and WAIS (b) core. The anomalies are defined with respect to the 50-year period before the individual eruptions. (c) Same as (a), but in addition the average signal is shown in blue, and the gray bands show the 16th to 84th percentiles of detrended time slices covering individual eruptions.


Using the much larger unipolar datasets, the difference in response is also seen clearly (Fig. 8c, d). However, this dataset (unlike the bipolar one) contains the Last Glacial Maximum, which features almost exclusively stadial conditions and which occurs during the younger part of the glacial where the signal preservation is better (Sect. 3.3). For a more fair comparison, we choose the interval 32–47.5 ka in the middle of our time period, which features an equal number of years with GI and GS conditions (Fig. S14). Even though reduced in NGRIP, the contrasting isotopic response is still significant in this interval (Fig. 8e, f). This difference in GI versus GS could be due to several factors.

  1. There is a different global or regional climate sensitivity to identical radiative forcing.

  2. The effective radiative forcing of (identical) sulfur-rich eruptions is different.

  3. The global or regional volcanic activity was different in GS versus GI.

  4. The dependence of δ18O on annual mean surface temperature in Greenland and Antarctica varied for GS and GI.

  5. The influence of factors other than annual mean temperature on δ18O anomalies is different in GS and GI.

Since the SNR in GI and GS is similar for most parts of the record (Fig. 4b), the increase in inferred volcanic cooling during GS compared to GI equals the increase in non-volcanic proxy variability, which is consistent with a state dependency of both climate and proxy sensitivity. A state dependence of climate sensitivity (point 1) would be an intriguing finding, but it is hard to rule out the confounding factors (points 2–5). Since the state dependency of the δ18O anomaly is opposite in Greenland and Antarctica, the potential state dependency is more likely in the regional climate sensitivity, i.e., the spatial pattern of the volcanic cooling, rather than in global average temperature. In the next section, we analyze differences in the volcanic forcing between GS and GI (point 3 and to some extent 2). In the section thereafter we employ records of relative snow accumulation rate in an attempt to gather more evidence for state dependency of the δ18O–temperature relationship (point 4), as well as for influences of relative accumulation rate changes on the δ18O signal (as part of point 5).

3.6 State dependency of volcanic forcing

There is generally a higher frequency of eruptions detected in Greenland during GS (Fig. 9a and see also Lin et al.2022). To some degree, this may be an artifact of the automatic detection of eruptions because the estimated eruption magnitudes could depend on the background noise level in the sulfate records, which is very different in GS and GI (Lin et al.2022). But for the average climate impact of eruptions only the relative distribution of the magnitudes counts and not the absolute frequency of the eruptions. The distribution of sulfate deposition in Greenland is skewed towards larger values during GS (Fig. 9b, c), whereas in Antarctica the distribution of GI eruptions is skewed to larger values (Fig. 9d, e).

This gives a consistent pattern with larger eruptions and more pronounced isotopic cooling in Greenland during GS and conversely larger eruptions and more cooling during GI in Antarctica. But in the following we show by resampling that the differences in sulfur deposition magnitude cannot explain the contrasting δ18O response. In particular, we resample the subset of eruptions with a larger average deposition (i.e., the GS eruptions for Greenland and the GI eruptions for Antarctica) with replacement such that they match the deposition magnitude distribution of the other subset with lower average deposition. From this resampled set of eruptions we calculate the average δ18O response and compare it to the subsets before resampling. The resampling method is explained in the Appendix of Lohmann and Svensson (2022), and it is similar to established Monte Carlo methods such as importance sampling, which aim to generate samples from a particular distribution, when only having observed samples from another distribution. The sulfur deposition distributions before and after resampling are shown in Fig. S15, and the resulting resampled average δ18O anomaly is shown by the dashed lines in Fig. 8a–b and e–f. In Greenland, the average isotopic anomaly is still more pronounced for the GS eruptions, and the same holds true for the GI eruptions in Antarctica. The results also hold for the other Greenland cores (see Fig. S16 for NEEM). Thus, the contrasting isotopic response in GI versus GS cannot be explained by the observed differences in the distribution of sulfate depositions. It may be that the latter is due to differences in sulfur transport and not the amount of sulfur ejected. There could be GI–GS differences in wind speeds and circulation patterns, which may or may not influence the aerosol climate forcing. Differences in atmospheric moisture may also modulate the lifetime and climate forcing of sulfur aerosols. A longer sulfate lifetime in the drier GS may be visible in broader Greenland sulfate peaks (Fig. S2c), but we cannot distinguish this from a broadening of the peaks due to the lower resolution during GS.

Figure 11(a) Scatter plot of the layer thickness change and δ18O anomaly in NGRIP for the unipolar dataset. The blue dots show the eruptions where both a negative δ18O anomaly and a layer thickness reduction are found. The red line and associated 95 % confidence interval are given by exponentiating the linear Deming regression line of Δlog λ and δ18O. A corresponding exponential CC relationship assuming α=0.8 is shown in green (see main text for more details). (b) Δlog λ and δ18O anomalies averaged in bins, which are given by every 5th percentile of the δ18O data. A linear regression without (with) intercept is shown in blue (red).


3.7 State-dependent volcanic impact on accumulation rate

Due to the unknown and potentially varying sensitivity α of the δ18O proxy to temperature changes, the implicated state dependency of the volcanic cooling may be spurious. To get additional evidence, we reconstruct changes in precipitation after the eruptions. Precipitation changes are expected to follow radiatively induced changes in temperature, since the atmospheric moisture capacity varies exponentially with temperature (Clausius–Clapeyron relation, CC). Indeed, volcanic cooling leads to a reduction in precipitation due to a weakened hydrological cycle (Robock and Liu1994; Bala et al.2008). In the polar regions, short-term relative changes in snow accumulation rates λ can be almost directly monitored in layer-counted ice cores by comparing the average layer thickness (implied by the depths of the counted annual layers) of nearby time intervals. Unlike δ18O, this is a direct measurement and its annual resolution is only slightly blurred by the imprecisions of the layer identification. We follow CC by assuming a change in λ with temperature,

(1) λ T = γ λ ,

and obtain λeγT, where γ is the accumulation sensitivity. Thus, the logarithm of the ratio of λ before and after a temperature change ΔT=T-T0 is linearly related to ΔT and by extension to the measured δ18O change for a given isotope sensitivity α:

(2) Δ log λ log λ ( T ) λ ( T 0 ) = γ Δ T = γ α Δ δ 18 O .

γCC=0.073 would be found when deriving Eq. (1) from a linearized Clausius–Clapeyron equation, assuming that the total precipitation amount is proportional to water vapor pressure. But the true value of γ for the climate system is lower and varies with location and T0 (Allen and Ingram2002).

We now consider anomalies with respect to the average state (T0, λ(T0)≡λ0) during the 50 years prior to the eruption. Figure 10a and b show the percentage change anomalies of λ, defined as (λλ0-1)100, for the NGRIP and WAIS unipolar datasets. Indeed, there are reduced accumulation rates in NGRIP associated with the eruptions in both GI and GS. The reduction is clearly more pronounced in GI. For WAIS, while the reduction in GS is not significant compared to the variability of the mean, there seems to be a more pronounced reduction in GI. However, the signal-to-noise ratio is low because the layer thickness is strongly affected by surface snow redistribution, and thus an average over a large number of eruptions is needed to extract the signal (compare mean signal and gray band in Fig. 10c). The seemingly delayed peak reduction in WAIS for GI eruptions might thus only be a random feature of the variability in the mean due to the small sample size of GI eruptions.

Focusing on NGRIP, the maximum layer thickness change averaged over all eruptions is 5.6 ± 0.9 %. The error is estimated by the standard deviation of the mean before the eruptions (fluctuations of the mean curve in Fig. 10c). This is larger than the 3 % global precipitation reduction inferred via sea level changes after five major eruptions during the last century (Grinsted et al.2007) or the modeled reductions of 1 %–2 % for the same eruptions (Iles and Hegerl2014). This could be due to an amplified polar response, but it also reflects the shorter return time of the eruptions in question compared to our unipolar dataset (20 vs. 65 years). The corresponding maximum δ18O anomaly is 1.54 ± 0.09 ‰, derived from the youngest quarter of eruptions (Fig. 3c) to minimize the effect of diffusion. This yields γα=0.037 [0.030, 0.046] (confidence band via the above standard deviations). Assuming the present-day value of α=0.8, this would correspond to an accumulation sensitivity of 2.9 [2.3, 3.6] % K−1.

An alternative estimate is obtained by regression of the anomalies of individual eruptions. This is shown in Fig. 11a, where a 4-year average accumulation anomaly (period of significant volcanic anomaly as determined from Fig. 10c) was used, along with a 10-year average isotopic cooling anomaly (period of significant anomaly determined from Fig. 3c). The exponential CC relationship assuming α=0.8 is shown in green. Clearly, the large data scatter may permit a variety of functional relationships. Nevertheless, we perform a fit to Eq. (2), where we take into account noise in both variables. The noise levels can be estimated via the SNR, as explained in Sect. 3.1. We find SNR=0.24 and SNR=0.44 for Δlog λ and Δδ18O, respectively. With the ratio of SNRs we can perform so-called Deming regression on the normalized data, which avoids underestimating the slope as in regular linear regression (attenuation bias). This yields γ/α=0.029±0.004, and assuming α=0.8 the accumulation sensitivity is 2.4 [2.0, 2.7] % K−1, in agreement with a model-derived global sensitivity of 2.4 % K−1 (Bala et al.2008). Note, however, that our given confidence interval does not reflect the significant freedom of choice in defining the average anomalies and performing the linear regression.

This accumulation sensitivity derived for the whole glacial ignores the clearly different accumulation reductions in GI and GS, where GI eruptions lead to a more pronounced reduction (Fig. 10a, b). In contrast, if our δ18O analysis reflects a genuine state dependency of the regional temperature response, we expect the stronger Greenland cooling in GS to yield a larger accumulation reduction. In NGRIP, the peak accumulation reduction is 4.3 ± 1.0 % in GS and 8.5 ± 1.5 % in GI, while the peak δ18O anomaly is 1.8 ± 0.1 ‰ and 0.78 ± 0.23 ‰, respectively (Fig. S17). If λ were perfectly proportional to ΔT at constant γ, this would imply a GS–GI contrast of the isotopic sensitivity of αGSαGI=4.5 [2.2, 9.8]. This large difference may be unrealistic, indicating that the accumulation sensitivity may also not be constant over time, as suggested by a previous analysis of the WAIS core (Fudge et al.2016).

In all above estimates of the sensitivity we assumed that a vanishing δ18O anomaly is accompanied by a vanishing λ anomaly (as in Eq. 2). However, when reducing the noise level by averaging the data in δ18O bins we can see that the linear relationship does not pass through the origin (Fig. 11b). Thus, the response of either of the proxies includes one or more processes that are not directly dependent on the underlying temperature change. This underlines the fact that the true values and state dependencies of γ and α cannot be revealed here. Nevertheless, on a qualitative level, the state dependency of Δλ in GI versus GS, which is opposite to the state dependency of Δδ18O, strongly suggests a state-dependent climate response to volcanic eruptions, albeit of a more complicated nature than simple differences in the annual (regional or global) temperature response. The state-dependent accumulation rate reduction also makes it plausible that the seasonality of precipitation after volcanic eruptions may be altered in a different way for GI and GS, which could in turn partly explain the differences in annual mean δ18O anomalies.

4 Discussion and conclusions

Here we attempt for the first time to quantify the volcanic cooling following large eruptions during the last glacial period from ice-core data. This is done by precise alignment of two recent datasets of volcanism to high-resolution δ18O records from the same ice cores. Going back in time far beyond the observational and historical periods enables us to investigate the impact of eruptions with very large return times. We find that the volcanic cooling signal is preserved in the ice-core records (Sect. 3.1 and 3.2), highlighting their potential to constrain the climatic impact of past volcanic eruptions in addition to tree-ring and lake sediment records (Sigl et al.2015; Tejedor et al.2021). However, the preservation depends critically on a high measurement resolution of the δ18O records, a high accumulation rate at the ice-core site, and a moderate thinning of the annual layers (Sect. 3.3). Further, detecting a sharp multi-annual cooling relies on high-resolution sulfur and conductivity records, which are used to define the precise depths of the volcanic eruptions in the ice cores. Not all cores used here fulfill these criteria. Given these limitations, we find that the observed isotopic anomaly after individual eruptions with centennial return periods is smaller than the high-frequency variability of the proxies (Figs. 1 and 2). The latter comprises multi-annual internal climate variability and non-climatic noise. As a result, we cannot give reliable estimates for the isotopic anomaly associated with individual eruptions and therefore also cannot estimate the cooling effect. Even the average anomaly at the time of deposition cannot be fully reconstructed since the signal degrades over time in a way that is not well-understood (Figs. 3, 4a and S10).

With this caveat in mind, the amplitude of the Greenland isotopic response to bipolar eruptions during the younger half of the investigated time interval is consistent with their observed return period of 500 years, since the four largest eruptions of the last 2000 years show roughly the same isotopic anomaly (Sect. 3.3). On the other hand, as the glacial records feature more diffusion, lower resolution, and less accurate alignment to the eruptions, the true isotopic anomaly should be larger. Indeed, the youngest glacial eruptions in the larger unipolar dataset show a clearly larger isotopic signal compared to the largest Common Era eruptions, despite a much lower return time of approximately 65 years (Fig. 3c). To better understand these differences, an in-depth comparison to eruptions of similar sulfate deposition covering the entire Holocene (Sigl et al.2022) may be helpful in a future study.

Eruptions with larger sulfate deposition magnitude also lead to increased δ18O cooling anomalies (Figs. 5 and 6a). Due to the large noise levels, we cannot determine with confidence whether this relationship is linear. Future studies with larger datasets covering longer periods should be able to reveal whether eruptions with increasing return times simply have a linearly increasing amplitude and/or duration of volcanic cooling or whether this relationship could be nonlinear and potentially have effects beyond a short-term cooling by compounding climatic regime shifts (tipping points). To do this it may be necessary to complement our methodology with idealized modeling of the proxy degradation over time.

By separating the data into eruptions occurring during the cold GS and milder GI periods, we find that the Greenland δ18O anomaly is larger during GS, while, on the other hand, the anomaly in the Antarctic EDC core is larger in GI (Fig. 8). This suggests a state-dependent climate response with more pronounced Greenland (Antarctic) cooling following eruptions during GS (GI). If the state dependency is indeed robust, the pronounced Greenland cooling during GS eruptions may play a role in the apparent influence of bipolar eruptions on the transitions from GS to GI (Lohmann and Svensson2022). But our results are also compatible with a more complicated difference in the climate response that is encoded in different sensitivities of the δ18O proxy to the volcanic cooling. In addition, the volcanic forcing could be state-dependent as a result of differences in atmospheric moisture and circulation or of a modulation of the volcanic activity by the climate state (Cooper et al.2018; Swindles et al.2018; Farquharson and Amelung2022). We indeed find slightly larger sulfur deposition estimates in Greenland (Antarctica) during GS (GI) (Fig. 9). However, this cannot explain the state-dependent δ18O anomalies, as shown in Sect. 3.6 by resampling the data such that the samples of eruptions during GI and GS have an equivalent distribution of sulfur deposition magnitudes.

It remains possible that the differences in δ18O arise despite an identical climate response in GI and GS, for instance due to a fixed seasonality of the volcanic cooling (Robock2000) in combination with different seasonalities of precipitation for GS and GI (Steig et al.1994; Werner et al.2000; Li et al.2005; Andersen et al.2006). In particular, if there is less winter precipitation in GS compared to GI, a less pronounced volcanic cooling (or even a warming) in winter compared to summer (equally in GI and GS) would give a more depleted annual mean δ18O signal in GS relative to GI. A similar situation could arise if there are different average precipitation source areas in GI and GS, for instance due to the differences in sea ice extent. If there is a latitudinal gradient in the volcanic cooling (Pausata et al.2020), this could mean that the change in temperature gradient from source to sink after an eruption would be higher in GS, which also results in more depleted δ18O.

Due to the shortcoming of unknown glacial δ18O sensitivity we also analyzed changes in accumulation rate after the eruptions (Sect. 3.7). While precipitation is generally believed to decrease proportionally to atmospheric cooling, we find that accumulation decreases in WAIS and NGRIP are clearly larger during GI eruptions, in contrast to the larger GS cooling suggested by Greenland δ18O. This reinforces the idea that there is a kind of state dependency, but the opposing tendencies cast doubt on whether the larger GS δ18O anomaly reflects more pronounced Greenland volcanic cooling. Since a vanishing volcanic δ18O anomaly does not coincide with a vanishing accumulation anomaly (Fig. 11b), it is clear that at least one of the two does not depend on temperature in a simple way. Just like δ18O, the local accumulation rate can be influenced by many factors apart from local temperature. Our analysis cannot reveal these factors, leaving the sensitivities of the δ18O proxy and of the accumulation rate to temperature unknown. An extension of our analysis to other ice-core proxies may give further insights into the climate response. Besides the response, the actual climate forcing of large volcanic eruptions can be much more varied compared to the simple surface cooling and drying assumed here, as evidenced by the recent Hunga–Tonga Hunga eruption (Millán et al.2022).

Nevertheless, we provide a proof of concept to use ice-core proxy records in assessing the multi-annual climate response to volcanic eruptions, as well as its change with time and climate background state. The provided observational evidence of a state-dependent response of δ18O and accumulation rate may be tested in studies with comprehensive climate models. Previous modeling argues both for and against a state dependency of the global climate response to volcanic eruptions (Zanchettin et al.2013; Berdahl and Robock2013; Muthers et al.2015; Ellerhoff et al.2022), as well as for state dependencies with opposite sign in a future warm climate (Fasullo et al.2017; Hopfcroft et al.2018). A state dependency with respect to DO cycles has not been investigated yet. Studies considering volcanic eruptions in models that can simulate glacial DO-like switches between GI and GS states (Vettoretti and Peltier2016; Klockmann et al.2020; Zhang et al.2021; Kuniyoshi et al.2022; Armstrong et al.2023) would be helpful. To test our results, direct comparisons with oxygen-isotope-enabled climate models (Schmidt et al.2007; LeGrande and Schmidt2009) should be performed, especially by extending previous studies on volcanic eruptions (Colose et al.2016) to other climate background states. Understanding the spatial patterns of the precipitation and temperature response to volcanic perturbations at different latitudes (Colose et al.2016), as well as the seasonality of the response (Zambri et al.2017), may help explain our results. Further, comparisons of our results with data from non-ice-core archives are needed. Tree-ring records, for instance, are beginning to reach into the last glacial (Reinig et al.2018; Pauly et al.2020; Reinig et al.2021) and may thus be used to assess the state dependency of the climate response to volcanic eruptions in future studies.

The presented methodology may also foster studies on climate variability and signal preservation in proxy records. Together with constraints on the strength of volcanic forcing, variability in climate records could be calibrated by the average volcanic climate response signal. Our preliminary analysis based on the signal-to-noise ratio suggests that the increase in the volcanic Greenland δ18O response during GS compared to GI is roughly the same as the increase in the non-volcanic proxy variability (Fig. 4b). Assuming equal volcanic forcing, one might thus speculate that the much-discussed state dependency of climate variability inferred from Greenland ice cores (Ditlevsen et al.1996; Rehfeld et al.2018) is due to a state-dependent proxy sensitivity. But more detailed modeling of the proxy evolution over time is required to make a fair comparison between GI and GS states, as well as glacial and interglacial periods. Specifically, it would be insightful to model the post-depositional alteration and subsequent diffusion of an idealized volcanic cooling signal and compare this to the observed average signals reported here.

In summary, we show that multi-annual cooling after major volcanic eruptions is preserved in high-resolution δ18O records of polar ice cores. The inferred average δ18O anomaly remains smaller than the proxy variability, however. This may suggest that volcanism is not the main driver of multi-annual and decadal temperature variability during the last glacial, as opposed to what has been found from tree-ring records during the Common Era (Sigl et al.2015). However, the temperature change at the time of eruption is uncertain due to attenuation of the volcanic δ18O signal over time and an unknown sensitivity of the proxy. In addition, the glacial δ18O variability is inflated by non-climatic noise resulting from low accumulation rates. The Greenland δ18O cooling anomaly during the cold GS periods is larger than during the milder GI. The opposite holds for Antarctica. This may indicate that the climate response to the radiative cooling of the eruptions is state-dependent. But due to other effects, such as precipitation seasonality, it may also be the sensitivity of annual mean δ18O to the volcanic cooling that is state-dependent. Post-eruptive cooling is accompanied by a reduction in ice-core accumulation rates. In contrast to the pattern observed in δ18O, GI periods feature a larger volcanic accumulation reduction than GS. The mechanisms behind this complicated state dependency of the post-eruptive ice-core signals cannot be revealed here. This may be achieved in future studies that test our observations in climate models or analyze the volcanic signals in additional proxies. Further usage of the volcanic cooling signal to understand the climate variability implied by the δ18O proxy may also be fruitful, especially as larger volcanic datasets become available.

Code and data availability

The bipolar volcanic record is available in the Supplement of Svensson et al. (2020), and the unipolar records are available in the Supplement of Lin et al. (2022). The high-resolution oxygen ice-core records of the individual cores are publicly available in the following online resources. NGRIP: (NGRIP Members2004); GISP2: (White2009); NEEM: (Gkinis et al.2020); EDML: (EPICA Community Members2010); EDC: (Jouzel and Masson-Delmotte2007); WAIS: (White et al.2019). The GRIP record is available upon request from the corresponding author. The high-resolution sulfate records shown in the Supplement are available in the following online resources. NGRIP: Supplement of Lin et al. (2022); WAIS: (McConnell2017); NEEM: Supplement of Schüpbach et al. (2018); GISP2: (Mayewski1999); EDC: (Severi2020). Code created for the statistical analyses can be found at (last access: 31 January 2024; DOI:, Lohmann2024).


The supplement related to this article is available online at:

Author contributions

JLo designed and performed the research. JLi and AS analyzed the volcanic sulfur peaks. The paper was written by JLo with input from all co-authors. All authors discussed and interpreted the results.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Views and opinions expressed are those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


We thank Vasileios Gkinis for help with the Antarctic high-resolution δ18O records.

Financial support

This research has been supported by the European Union's Horizon 2020 research and innovation program under grant agreement no. 820970 (TiPES) and by Danmarks Frie Forskningsfond (grant no. 2032-00346B). Anders Svensson was funded by the European Union (ERC, Green2Ice, 101072180).

Review statement

This paper was edited by Amaelle Landais and reviewed by two anonymous referees.


Allen, M. R. and Ingram, W. J.: Constraints on future changes in climate and the hydrologic cycle, Nature, 419, 224–232, 2002. a

Andersen, K. K., Svensson, A., Johnsen, S. J., Rasmussen, S. O., Bigler, M., Röthlisberger, R., Ruth, U., Siggaard-Andersen, M.-L., Steffensen, J. P., Dahl-Jensen, D., Vinther, B. M., and Clausen, H. B.: The Greenland Ice Core Chronology 2005, 15–42 ka. Part 1: constructing the time scale, Quaternary Sci. Rev., 25, 3246–3257, 2006. a, b

Armstrong, E., Izumi, K., and Valdes, P.: Identifying the mechanisms of DO-scale oscillations in a GCM: a salt oscillator triggered by the Laurentide ice sheet, Clim. Dynam., 60, 3983–4001,, 2023. a

Ashwin, P. and von der Heydt, A. S.: Extreme Sensitivity and Climate Tipping Points, J. Stat. Phys., 179, 1531–1552, 2020. a

Aubry, T. J., Staunton-Sykes, J., Marshall, L. R., Haywood, J., Abraham, N. L., and Schmidt, A.: Climate change modulates the stratospheric volcanic sulfate aerosol lifecycle and radiative forcing from tropical eruptions, Nat. Commun., 12, 4708,, 2021. a

Bala, G., Duffy, P. B., and Taylor, K. E.: Impact of geoengineering schemes on the global hydrological cycle, P. Natl. Acad. Sci. USA, 105, 7664–7669, 2008. a, b, c

Bender, F. A.-M., Ekman, A. M. L., and Rodhe, H.: Response to the eruption of Mount Pinatubo in relation to climate sensitivity in the CMIP3 models, Clim. Dynam., 35, 875–886, 2010. a

Berdahl, M. and Robock, A.: Northern Hemispheric cryosphere response to volcanic eruptions in the Paleoclimate Modeling Intercomparison Project 3 last millennium simulations, J. Geophys. Res., 118, 12359–12370, 2013. a

Boer, G. J., Stowasser, M., and Hamilton, K.: Inferring climate sensitivity from volcanic events, Clim. Dynam., 28, 481–502, 2007. a

Buizert, C., Gkinis, V., Severinghaus, J. P., He, F., Lecavalier, B. S., Kindler, P., Leuenberger, M., Carlson, A. E., Vinther, B., Masson-Delmotte, V., White, J. W. C., Liu, Z., Otto-Bliesner, B., and Brook, E. J.: Greenland temperature response to climate forcing during the last deglaciation, Science, 345, 1177–1180, 2014. a, b, c

Buizert, C., Cuffey, K. M., Severinghaus, J. P., Baggenstos, D., Fudge, T. J., Steig, E. J., Markle, B. R., Winstrup, M., Rhodes, R. H., Brook, E. J., Sowers, T. A., Clow, G. D., Cheng, H., Edwards, R. L., Sigl, M., McConnell, J. R., and Taylor, K. C.: The WAIS Divide deep ice core WD2014 chronology – Part 1: Methane synchronization (68–31 ka BP) and the gas age–ice age difference, Clim. Past, 11, 153–173,, 2015. a

Buizert, C., Fudge, T. J., Roberts, W. H. G. et al.: Antarctic surface temperature and elevation during the Last Glacial Maximum, Science, 372, 1097–1101, 2021. a

Burke, A., Moore, K. A., Sigl, M., Nita, D. C., McConnell, J. R., and Adkins, J. F.: Stratospheric eruptions from tropical and extra-tropical volcanoes constrained using high-resolution sulfur isotopes in ice cores, Earth Planet. Sc. Lett., 521, 113–119, 2019. a, b

Caballero, R. and Huber, M.: State-dependent climate sensitivity in past warm climates and its implications for future climate projections, P. Natl. Acad. Sci. USA, 110, 14162–14167, 2013. a

Castellano, E., Becagli, S., Jouzel, J., Migliori, A., Severi, M., Steffensen, J. P., Traversi, R., and Udisti, R.: Volcanic eruption frequency over last 45 ky as recorded in Epica-Dome C ice core (East Antarctica) and its relationship with climatic changes, Global Planet. Change, 42, 195–205, 2004. a

Cauquoin, A., Abe-Ouchi, A., Obase, T., Chan, W.-L., Paul, A., and Werner, M.: Effects of Last Glacial Maximum (LGM) sea surface temperature and sea ice extent on the isotope–temperature slope at polar ice core sites, Clim. Past, 19, 1275–1294,, 2023. a

Colose, C. M., LeGrande, A. N., and Vuille, M.: Hemispherically asymmetric volcanic forcing of tropical hydroclimate during the last millennium, Earth Syst. Dynam., 7, 681–696,, 2016. a, b

Cooper, C. L., Swindles, G. T., Savov, I. P., Schmidt, A., and Bacon, K. L.: Evaluating the relationship between climate change and volcanism, Earth Sci. Rev., 177, 238–247, 2018. a

Cuffey, K. M. and Steig, E. J.: Isotopic diffusion in polar firn: implications for interpretation of seasonal climate paratneters in ice-core records, with emphasis on central Greenland, J. Glaciol., 44, 273–284, 1998. a

Ditlevsen, P. D., Svensmark, H., and Johnsen, S.: Contrasting atmospheric and climate dynamics of the last-glacial and Holocene periods, Nature, 379, 810,, 1996. a

Ellerhoff, B., Kirschner, M. J., Ziegler, E., Holloway, M. D., Sime, L., and Rehfeld, K.: Contrasting State-Dependent Effects of Natural Forcing on Global and Local Climate Variability, Geophys. Res. Lett., 49, e2022GL098335,, 2022. a, b

EPICA Community Members: One-to-one coupling of glacial climate variability in Greenland and Antarctica, Nature, 444, 195–198, 2006. a

EPICA Community Members: Stable oxygen isotopes of ice core EDML, PANGAEA [data set],, 2010. a

Farquharson, J. I. and Amelung, F.: Volcanic hazard exacerbated by future global warming-driven increase in heavy rainfall, R. Soc. Open Sci., 9, 220275,, 2022. a

Fasullo, J. T., Tomas, R., Stevenson, S., Otto-Bliesner, B., Brady, E., and Wahl, E.: The amplifying influence of increased ocean stratification on a future year without a summer, Nature Comm., 8, 1236,, 2017. a, b

Fudge, T. J., Markle, B. R., Cuffey, K. M., Buizert, C., Taylor, K. C., Steig, E. J., Waddington, E. D., Conway, H., and Koutnik, M.: Variable relationship between accumulation and temperature in West Antarctica for the past 31,000 years, Geophys. Res. Lett., 43, 3795–3803, 2016. a

Gao, C., Oman, L., Robock, A., and Stenchikov, G. L.: Atmospheric volcanic loading derived from bipolar ice cores: Accounting for the spatial distribution of volcanic deposition, J. Geophys. Res., 112, D09109,, 2007. a

Gautier, E., Savarino, J., Erbland, J., Lanciki, A., and Possenti, P.: Variability of sulfate signal in ice core records based on five replicate cores, Clim. Past, 12, 103–113,, 2016. a

Gkinis, V., Simonsen, S. B., Buchardt, S. L., White, J. W. C., and Vinther, B. M.: Water isotope diffusion rates from the NorthGRIP ice core for the last 16,000 years – Glaciological and paleoclimatic implications, Earth Planet. Sc. Lett., 405, 132–141, 2014. a

Gkinis, V., Vinther, B. M., Quistgaard, T., et al.: NEEM ice core High Resolution (0.05 m) Water Isotope Ratios (18O/16O, 2H/1H) covering 8–129 ky b2k, PANGAEA [data set],, 2020. a

Goursaud, S., Masson-Delmotte, V., Favier, V., Preunkert, S., Legrand, M., Minster, B., and Werner, M.: Challenges associated with the climatic interpretation of water stable isotope records from a highly resolved firn core from Adélie Land, coastal Antarctica, The Cryosphere, 13, 1297–1324,, 2019. a

Grinsted, A., Moore, J. C., and Jevrejeva, S.: Observational evidence for volcanic impact on sea level and the global water cycle, P. Natl. Acad. Sci. USA, 104, 19730–19734, 2007. a

Guillevic, M., Bazin, L., Landais, A., Kindler, P., Orsi, A., Masson-Delmotte, V., Blunier, T., Buchardt, S. L., Capron, E., Leuenberger, M., Martinerie, P., Prié, F., and Vinther, B. M.: Spatial gradients of temperature, accumulation and δ18O-ice in Greenland over a series of Dansgaard–Oeschger events, Clim. Past, 9, 1029–1051,, 2013. a

Helama, S., Stoffel, M., Hall, R. J., Jones, P. D., Arppe, L., Matskovsky, V. V., Timonen, M., Nöjd, P., Mielikäinen, K., and Oinonen, M.: Recurrent transitions to Little Ice Age-like climatic regimes over the Holocene, Clim. Dynam., 56, 3817–3833, 2021. a

Hopfcroft, P. O., Kandlbauer, J., Valdes, P. J., and Sparks, R. S. J.: Reduced cooling following future volcanic eruptions, Clim. Dynam., 51, 1449–1463, 2018. a, b

Iles, C. E. and Hegerl, G. C.: The global precipitation response to volcanic eruptions in the CMIP5 models, Environ. Res. Lett., 9, 104012,, 2014. a

Johnsen, S. J., Clausen, H. B., Dansgaard, W., Gundestrup, N., Hammer, C. U., Andersen, U., Andersen, K. K., Hvidberg, C. S., Dahl-Jensen, D., Steffensen, J. P., Shoji, H., Sveinbjörnsdóttir, A. E., White, J., Jouzel, J., and Fisher, D.: The δ18O record along the Greenland Ice Core Project deep ice core and the problem of possible Eemian climatic instability, J. Geoph. Research, 102, 26397–26410, 1997. a

Johnsen, S. J., Clausen, H. B., Cuffey, K. M., Hoffmann, G., Schwander, J., and Creyts, T.: Diffusion of stable isotopes in polar firn and ice: the isotope effect in firn diffusion, in: Physics of Ice Core Records, edited by: Hondoh, T., Hokkaido University Press, 2000, Sapporo, 121–140,, 2000. a

Jones, T. R., White, J. W. C., Steig, E. J., Vaughn, B. H., Morris, V., Gkinis, V., Markle, B. R., and Schoenemann, S. W.: Improved methodologies for continuous-flow analysis of stable water isotopes in ice cores, Atmos. Meas. Tech., 10, 617–632,, 2017. a, b

Jones, T. R., Roberts, W. H. G., Steig, E. J., Cuffey, K. M., Markle, B. R., and White, J. W. C.: Southern Hemisphere climate variability forced by Northern Hemisphere ice-sheet topography, Nature, 554, 351–355, 2018. a, b

Jouzel, J. and Masson-Delmotte, V.: EPICA Dome C Ice Core 800 KYr deuterium data and temperature estimates, PANGAEA [data set],, 2007. a

Jouzel, J., Masson-Delmotte, V., Cattani, O., Dreyfus, G., Falourd, S., Hoffmann, G., Minster, B., Nouet, J., Barnola, J. M., Chappellaz, J., Fischer, H., Gallet, J. C., Johnsen, S., Leuenberger, M., Loulergue, L., Luethi, D., Oerter, H., Parrenin, F., Raisbeck, G., Raynaud, D., Schilt, A., Schwander, J., Selmo, E., Souchez, R., Spahni, R., Stauffer, B., Steffensen, J. P., Stenni, B., Stocker, T. F., Tison, J. L., Werner, M., and Wolff, E. W.: Orbital and Millennial Antarctic Climate Variability over the Past 800,000 Years, Science, 317, 793–796, 2007. a

Klockmann, M., Mikolajewicz, U., Kleppin, H., and Marotzke, J.: Coupling of the Subpolar Gyre and the Overturning Circulation During Abrupt Glacial Climate Transitions, Geophys. Res. Lett., 47, e2020GL090361,, 2020. a

Kobashi, T., Menviel, L., Jeltsch-Thömmes, A., Vinther, B. M., Box, J. E., Muscheler, R., Nakaegawa, T., Pfister, P. L., Döring, M., Leuenberger, M., Wanner, H., and Ohmura, A.: Volcanic influence on centennial to millennial Holocene Greenland temperature change, Sci. Rep., 7, 1441,, 2017. a

Köhler, P., de Boer, B., von der Heydt, A. S., Stap, L. B., and van de Wal, R. S. W.: On the state dependency of the equilibrium climate sensitivity during the last 5 million years, Clim. Past, 11, 1801–1823,, 2015. a

Kuniyoshi, Y., Abe-Ouchi, A., Sherriff-Tadano, S., Chan, W.-L., and Saito, F.: Effect of climatic precession on Dansgaard-Oeschger-like oscillations, Geophys. Res. Lett., 49, e2021GL095695,, 2022. a

LeGrande, A. N. and Schmidt, G. A.: Sources of Holocene variability of oxygen isotopes in paleoclimate archives, Clim. Past, 5, 441–455,, 2009. a

Lenton, T. M., Held, H., Kriegler, E., Hall, J. W., Lucht, W., Rahmstorf, S., and Schellnhuber, H. J.: Tipping elements in the Earth’s climate system, P. Natl. Acad. Sci USA, 105, 1786–1793, 2008. a

Li, C., Battisti, D. S., Schrag, D. P., and Tziperman, E.: Abrupt climate shifts in Greenland due to displacements of the sea ice edge, Geophys. Res. Lett., 32, L19702,, 2005. a

Lin, J., Svensson, A., Hvidberg, C. S., Lohmann, J., Kristiansen, S., Dahl-Jensen, D., Steffensen, J. P., Rasmussen, S. O., Cook, E., Kjær, H. A., Vinther, B. M., Fischer, H., Stocker, T., Sigl, M., Bigler, M., Severi, M., Traversi, R., and Mulvaney, R.: Magnitude, frequency and climate forcing of global volcanism during the last glacial period as seen in Greenland and Antarctic ice cores (60–9 ka), Clim. Past, 18, 485–506,, 2022. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q

Liu, Z., He, C., Yan, M., Buizert, C., Otto-Bliesner, B. L., Lu, F., and Zeng, C.: Reconstruction of Past Antarctic Temperature Using Present Seasonal δ18O–Inversion Layer Temperature: Unified Slope Equations and Applications, J. Climate, 36, 2933–2957, 2023. a

Lohmann, J.: Python code used in the analysis of the publication Lohmann et al. “State-dependent impact of major volcanic eruptions observed in ice-core records of the last glacial period”, Zenodo [code],, 2024. a

Lohmann, J. and Svensson, A.: Ice core evidence for major volcanic eruptions at the onset of Dansgaard–Oeschger warming events, Clim. Past, 18, 2021–2043,, 2022. a, b, c

Mann, M. E., Steinman, B. A., Brouillette, D. J., and Miller, S. K.: Multidecadal climate oscillations during the past millennium driven by volcanic forcing, Science, 371, 1014–1019, 2021. a

Mayewski, P. A.: GISP2 Ions: Deep (D) Core, PANGAEA [data set],, 1999. a

McConnell, J.: WAIS Divide Ice-Core Aerosol Records from 1300 to 3404 m, U.S. Antarctic Program (USAP) Data Center [data set],, 2017. a

McConnell, J. R., Sigl, M., Plunkett, G., Burke, A., Mi Kim, W., Raible, C. C., Wilson, A. I., Manning, J. G., Ludlow, F., Chellman, N. J., Innes, H. M., Yang, Z., Larsen, J. F., Schaefer, J. R., Kipfstuhl, S., Mojtabavi, S., Wilhelms, F., Opel, T., Meyer, H., and Steffensen, J. P.: Extreme climate after massive eruption of Alaska’s Okmok volcano in 43 BCE and effects on the late Roman Republic and Ptolemaic Kingdom, P. Natl. Acad. Sci. USA, 117, 15443–15449, 2020. a

Merlis, T. M., Held, I. M., Stenchikov, G. L., Zeng, F., and Horowitz, L. W.: Constraining Transient Climate Sensitivity Using Coupled Climate Model Simulations of Volcanic Eruptions, J. Climate, 27, 7781–7795, 2014. a

Millán, L., Santee, M. L., Lambert, A., Livesey, N. J., Werner, F., Schwartz, M. J., Pumphrey, H. C., Manney, G. L., Wang, Y., Su., H., Wu, L., Read, W. G., and Froidevaux, L.: The Hunga Tonga-Hunga Ha'apai Hydration of the Stratosphere, Geophys. Res. Lett., 49, e2022GL099381,, 2022. a

Münch, T., Kipfstuhl, S., Freitag, J., Meyer, H., and Laepple, T.: Regional climate signal vs. local noise: a two-dimensional view of water isotopes in Antarctic firn at Kohnen Station, Dronning Maud Land, Clim. Past, 12, 1565–1581,, 2016. a, b

Muthers, S., Arfeuille, F., Raible, C. C., and Rozanov, E.: The impacts of volcanic aerosol on stratospheric ozone and the Northern Hemisphere polar vortex: separating radiative-dynamical changes from direct effects due to enhanced aerosol heterogeneous chemistry, Atmos. Chem. Phys., 15, 11461–11476,, 2015. a

NGRIP Members: High-resolution record of Northern Hemisphere climate extending into the last interglacial period, Nature, 431, 147–151,, 2004 (data available at:, last access: 31 January 2024). a, b

PAGES 2k Consortium: Consistent multidecadal variability in global temperature reconstructions and simulations over the Common Era, Nat. Geosci., 12, 643–649, 2019. a

Pauly, M., Helle, G., Büntgen, U., Wacker, L., Treydte, K., Reinig, F., Turney, C., Nievergelt, D., Kromer, B., Friedrich, M., Sookdeo, A., Heinrich, I., Riedel, F., Balting, D., and Brauer, A.: An annual-resolution stable isotope record from Swiss subfossil pine trees growing in the late Glacial, Quaternary Sci. Rev., 247, 106550,, 2020. a

Pausata, F. S. R., Zanchettin, D., Karamperidou, C., Caballero, R., and Battisti, D. S.: ITCZ shift and extratropical teleconnections drive ENSO response to volcanic eruptions, Sci. Adv., 6, eaaz5006,, 2020. a

Rasmussen, S. O., Abbott, P. M., Blunier, T., Bourne, A. J., Brook, E., Buchardt, S. L., Buizert, C., Chappellaz, J., Clausen, H. B., Cook, E., Dahl-Jensen, D., Davies, S. M., Guillevic, M., Kipfstuhl, S., Laepple, T., Seierstad, I. K., Severinghaus, J. P., Steffensen, J. P., Stowasser, C., Svensson, A., Vallelonga, P., Vinther, B. M., Wilhelms, F., and Winstrup, M.: A first chronology for the North Greenland Eemian Ice Drilling (NEEM) ice core, Clim. Past, 9, 2713–2730,, 2013. a, b

Rasmussen, S. O., Dahl-Jensen, D., Fischer, H., Fuhrer, K., Hansen, S. B., Hansson, M., Hvidberg, C. S., Jonsell, U., Kipfstuhl, S., Ruth, U., Schwander, J., Siggaard-Andersen, M.-L., Sinnl, G., Steffensen, J. P., Svensson, A. M., and Vinther, B. M.: Ice-core data used for the construction of the Greenland Ice-Core Chronology 2005 and 2021 (GICC05 and GICC21), Earth Syst. Sci. Data, 15, 3351–3364,, 2023. a

Rehfeld, K., Münch, T., Ho, S. L., and Laepple, T.: Global patterns of declining temperature variability from the Last Glacial Maximum to the Holocene, Nature, 554, 356–359, 2018. a

Reinig, F., Nievergelt, D., Esper, J., Friedrich, M., Helle, G., Hellmann, L., Kromer, B., Morganti, S., Pauly, M., Sookdeo, A., Tegel, W., Treydte, K., Verstege, A., Wacker, L., and Büntgen, U.: New tree-ring evidence for the Late Glacial period from the northern pre-Alps in eastern Switzerland, Quaternary Sci. Rev., 186, 215–224, 2018. a

Reinig, F., Wacker, L., Jöris, O., Oppenheimer, C., Guidobaldi, G., Nievergelt, D., Adolphi, F., Cherubini, P., Engels, S., Esper, J., Land, A., Lane, C., Pfanz, H., Remmele, S., Sigl, M., Sookdeo, A., and Büntgen, U.: Precise date for the Laacher See eruption synchronizes the Younger Dryas, Nature, 595, 66–69, 2021. a

Robock, A.: Volcanic eruptions and climate, Rev. Geophys., 38, 191–219, 2000. a

Robock, A. and Liu, Y.: The Volcanic Signal in Goddard Institute for Space Studies Three-Dimensional Model Simulations, J. Climate, 7, 44–55, 1994. a, b

Schmidt, G. A., LeGrande, A. N., and Hoffmann, G.: Water isotope expressions of intrinsic and forced variability in a coupled ocean-atmosphere model, J. Geophys. Res., 112, D10103,, 2007. a

Schüpbach, S., Fischer, H., Bigler, M. et al.: Greenland records of aerosol source and atmospheric lifetime changes from the Eemian to the Holocene, Nat. Commun., 9, 1476,, 2018. a

Schurer, A. P., Tett, S. F. B., and Hegerl, G. C.: Small influence of solar variability on climate over the past millennium, Nat. Geosci., 7, 104–108, 2014. a

Seierstad, I. K., Abbott, P. M., Bigler, M., Blunier, T., Bourne, A. J., Brook, E., Buchardt, S. L., Buizert, C., Clausen, H. B., Cook, E., Dahl-Jensen, D., Davies, S. M., Guillevic, M., Johnsen, S. J., Pedersen, D. S., Popp, T. J., Rasmussen, S. O., Severinghaus, J. P., Svensson, A., and Vinther, B. M.: Consistently dated records from the Greenland GRIP, GISP2 and NGRIP ice cores for the past 104 ka reveal regional millennial-scale δ18O gradients with possible Heinrich event imprint, Quaternary Sci. Rev., 106, 29–46, 2014. a

Severi, M., Udisti, R., Castellano, E., Wolff, E. W.: NOAA/WDS Paleoclimatology – EPICA Dome C 203,000 Year High-Resolution FIC Sulfate Data, NOAA National Centers for Environmental Information [data set],, 2020. a

Sigl, M., Winstrup, M., McConnell, J. R. et al.: Timing and climate forcing of volcanic eruptions for the past 2,500 years, Nature, 523, 543–549, 2015. a, b, c, d, e, f, g, h

Sigl, M., Fudge, T. J., Winstrup, M., Cole-Dai, J., Ferris, D., McConnell, J. R., Taylor, K. C., Welten, K. C., Woodruff, T. E., Adolphi, F., Bisiaux, M., Brook, E. J., Buizert, C., Caffee, M. W., Dunbar, N. W., Edwards, R., Geng, L., Iverson, N., Koffman, B., Layman, L., Maselli, O. J., McGwire, K., Muscheler, R., Nishiizumi, K., Pasteris, D. R., Rhodes, R. H., and Sowers, T. A.: The WAIS Divide deep ice core WD2014 chronology – Part 2: Annual-layer counting (0–31 ka BP), Clim. Past, 12, 769–786,, 2016. a

Sigl, M., Toohey, M., McConnell, J. R., Cole-Dai, J., and Severi, M.: Volcanic stratospheric sulfur injections and aerosol optical depth during the Holocene (past 11 500 years) from a bipolar ice-core array, Earth Syst. Sci. Data, 14, 3167–3196,, 2022. a, b

Sjolte, J., Hoffmann, G., Johnsen, S. J., Vinther, B. M., Masson-Delmotte, V., and Sturm, C.: Modeling the water isotopes in Greenland precipitation 1959–2001 with the meso-scale model REMO-iso, J. Geophys. Res., 116, D18105,, 2011. a, b

Steig, E. J., Grootes, P. M., and Stuiver, M.: Seasonal precipitation timing and ice core records, Science, 266, 1885–1887, 1994. a

Stuiver, M. and Grootes, P. M.: GISP2 Oxygen Isotope Ratios, Quaternary Res., 53, 277–284, 2000. a

Svensson, A., Andersen, K. K., Bigler, M., Clausen, H. B., Dahl-Jensen, D., Davies, S. M., Johnsen, S. J., Muscheler, R., Rasmussen, S. O., Röthlisberger, R., Steffensen, J. P., and Vinther, B. M.: The Greenland Ice Core Chronology 2005, 15–42 ka. Part 2: comparison to other records, Quaternary Sci. Rev., 25, 3258–3267, 2006. a

Svensson, A., Andersen, K. K., Bigler, M., Clausen, H. B., Dahl-Jensen, D., Davies, S. M., Johnsen, S. J., Muscheler, R., Parrenin, F., Rasmussen, S. O., Röthlisberger, R., Seierstad, I., Steffensen, J. P., and Vinther, B. M.: A 60 000 year Greenland stratigraphic ice core chronology, Clim. Past, 4, 47–57,, 2008. a

Svensson, A., Dahl-Jensen, D., Steffensen, J. P., Blunier, T., Rasmussen, S. O., Vinther, B. M., Vallelonga, P., Capron, E., Gkinis, V., Cook, E., Kjær, H. A., Muscheler, R., Kipfstuhl, S., Wilhelms, F., Stocker, T. F., Fischer, H., Adolphi, F., Erhardt, T., Sigl, M., Landais, A., Parrenin, F., Buizert, C., McConnell, J. R., Severi, M., Mulvaney, R., and Bigler, M.: Bipolar volcanic synchronization of abrupt climate change in Greenland and Antarctic ice cores during the last glacial period, Clim. Past, 16, 1565–1580,, 2020. a, b, c, d, e, f, g, h, i, j, k

Swindles, G. T., Watson, E. J., Savov, I. P., Lawson, I. T., Schmidt, A., Hooper, A., Cooper, C. L., Connor, C. B., Gloor, M., and Carrivick, J. L.: Climatic control on Icelandic volcanic activity during the mid-Holocene, Geology, 47, 47–50, 2018. a

Tejedor, E., Steiger, N. J., Smerdon, J. E., Serrano-Notivoli, R., and Vuille, M.: Global hydroclimatic response to tropical volcanic eruptions over the last millennium, P. Natl. Acad. Sci. USA, 118, e2019145118,, 2021. a

Uemura, R., Masson-Delmotte, V., Jouzel, J., Landais, A., Motoyama, H., and Stenni, B.: Ranges of moisture-source temperature estimated from Antarctic ice cores stable isotope records over glacial–interglacial cycles, Clim. Past, 8, 1109–1125,, 2012. a

Vega, C. P., Schlosser, E., Divine, D. V., Kohler, J., Martma, T., Eichler, A., Schwikowski, M., and Isaksson, E.: Surface mass balance and water stable isotopes derived from firn cores on three ice rises, Fimbul Ice Shelf, Antarctica, The Cryosphere, 10, 2763–2777,, 2016. a

Vettoretti, G. and Peltier, W. R.: Thermohaline instability and the formation of glacial North Atlantic super polynyas at the onset of Dansgaard-Oeschger warming events, Geophys. Res. Lett., 43, 5336–5344, 2016. a

Vinther, B. M., Clausen, H. B., Johnsen, S. J., Rasmussen, S. O., Andersen, K. K., Buchardt, S. L., Dahl-Jensen, D., Seierstad, I. K., Siggaard-Andersen, M.-L., Steffensen, J. P., Svensson, A., Olsen, J., and Heinemeier, J.: A synchronized dating of three Greenland ice cores throughout the Holocene, J. Geophys. Res., 111, D13102,, 2006.  a

Vinther, B. M., Jones, P. D., Briffa, K. R., Clausen, H. B., Andersen, K. K., Dahl-Jensen, D., and Johnsen, S. J.: Climatic signals in multiple highly resolved stable isotope records from Greenland, Quaternary Sci. Rev., 29, 522–538, 2010. a

von der Heydt, A. S., Dijkstra, H. A., van de Wal, R. S. W., Caballero, R., Crucifix, M., Foster, G. L., Huber, M., Köhler, P., Rohling, E., Valdes, P. J., Ashwin, P., Bathiany, S., Berends, T., van Bree, L. G. J., Ditlevsen, P., Ghil, M., Haywood, A. M., Katzav, J., Lohmann, G., Lohmann, J., Lucarini, V., Marzocchi, A., Pälike, H., Baroni, I. R., Simon, D., Sluijs, A., Stap, L. B., Tantet, A., Viebahn, J., and Ziegler, M.: Lessons on Climate Sensitivity From Past Climate Changes, Curr. Clim. Change Rep., 2, 148–158, 2016. a

Werner, M., Mikolajewicz, U., Heimann, M., and Hoffmann, G.: Borehole versus isotope temperatures on Greenland: Seasonality does matter, Geophys. Res. Lett, 27, 723–726, 2000. a

White, J.: GISP2 Stable Isotopes, Version 1.0, UCAR/NCAR – Earth Observing Laboratory [data set],, 2009. a

White, J., Bradley, E., Garland, J., Jones, T. R., Morris, V., Price, M., and Vaughn, B.: Stable Isotopes of Ice in the Transition and Glacial Sections of the WAIS Divide Deep Ice Core, U.S. Antarctic Program (USAP) Data Center [data set],, 2019. a

Zambri, B., LeGrande, A. N., Robock, A., and Slawinska, J.: Northern Hemisphere winter warming and summer monsoon reduction after volcanic eruptions over the last millennium, J. Geophys. Res.-Atmos., 122, 7971–7989, 2017. a

Zanchettin, D., Bothe, O., Graf, H. F., Lorenz, S. J., Luterbach, J., Timmreck, C., and Jungclaus, J. H.: Background conditions influence the decadal climate response to strong volcanic eruptions, J. Geophys. Res., 118, 4090–4106, 2013. a

Zhang, X., Barker, S., Knorr, G., Lohmann, G., Drysdale, R., Sun, Y., Hodell, D., and Chen, F.: Direct astronomical influence on abrupt climate variability, Nat. Geosci., 14, 819–839, 2021. a

Zielinski, G. A., Mayewski, P. A., Meeker, L. D., Grönvold, K., Germani, M. S., Whitlow, S., Twickler, M. S., and Taylor, K.: Volcanic aerosol records and tephrochronology of the Summit, Greenland, ice cores, J. Geophys. Res., 102, 26625,, 1997. a

Short summary
We present the first attempt to constrain the climatic impact of volcanic eruptions with return periods of hundreds of years by the oxygen isotope records of Greenland and Antarctic ice cores covering the last glacial period. A clear multi-annual volcanic cooling signal is seen, but its absolute magnitude is subject to the unknown glacial sensitivity of the proxy. Different proxy signals after eruptions during cooler versus warmer glacial stages may reflect a state-dependent climate response.