Bipolar volcanic synchronization of abrupt climate change in Greenland and Antarctic ice cores during the last glacial period

The last glacial period is characterized by a number of millennial climate events that have been identified in both Greenland and Antarctic ice cores and that are abrupt in Greenland climate records. The mechanisms governing this climate variability remain a puzzle that requires a precise synchronization of ice cores from the two hemispheres to be resolved. Previously, Greenland and Antarctic ice cores have been synchronized primarily via their common records of gas concentrations or isotopes from the trapped air and via cosmogenic isotopes measured on the ice. In this work, we apply ice core volcanic proxies and annual layer counting to identify large volcanic eruptions that have left a signature in both Greenland and Antarctica. Generally, no tephra is associated with those eruptions in the ice cores, so the source of the eruptions cannot be identified. Instead, we identify and match sequences of volcanic eruptions with bipolar distribution of sulfate, i.e. unique patterns of volcanic events separated by the same number of years at the two poles. Using this approach, we pinpoint 82 large bipolar volcanic eruptions throughout the second half of the last glacial period (12–60 ka). This improved ice core synchronization is applied to determine the bipolar phasing of abrupt climate change events at decadal-scale precision. In response to Greenland abrupt climatic transitions, we find a response in the Antarctic water isotope signals (δ¹⁸O and deuterium excess) that is both more immediate and more abrupt than that found with previous gas-based interpolar synchronizations, providing additional support for our volcanic framework. On average, the Antarctic bipolar seesaw climate response lags the midpoint of Greenland abrupt δ¹⁸O transitions by 122±24 years. The time difference between Antarctic signals in deuterium excess and δ¹⁸O, which likewise informs the time needed to propagate the signal as described by the theory of the bipolar seesaw but is less sensitive to synchronization errors, suggests an Antarctic δ¹⁸O lag behind Greenland of 152±37 years. These estimates are shorter than the 200 years suggested by earlier gas-based synchronizations. As before, we find variations in the timing and duration between the response at different sites and for different events suggesting an interaction of oceanic and atmospheric teleconnection patterns as well as internal climate variability.

1 Introduction

Greenland and Antarctic ice cores provide high-resolution records of abrupt climate events occurring throughout the last glacial period (11.7–115 ka). In Greenland ice cores, Dansgaard–Oeschger (DO) events describe a series of characteristic climate events (Dansgaard et al., 1993; North Greenland Ice Core Project members, 2004) that involve warming transitions of up to 16.5 ^∘C (Kindler et al., 2014) occurring within decades (Erhardt et al., 2019). Each DO event consists of a relatively mild climatic period, referred to as a Greenland Interstadial (GI), that is followed by a cold climatic period, known as a Greenland Stadial (GS). The duration of GIs and GSs range from centuries to millennia. Detailed investigation of the stratigraphy of the 25 major DO events originally identified has revealed that some of the events are composed of several separate warming and cooling events, leading to a total of 31–33 abrupt warming events during the last glacial period depending on the definition employed (Rasmussen et al., 2014). Whereas the onset of a GI event is abrupt and occurs in less than a century, the cooling transitions from GI to GS are more gradual and typically occur over several centuries. DO events are believed to originate in the North Atlantic but have a global climatic impact that is documented in a wide range of paleoclimate archives across the Northern (Voelker and workshop participants, 2002) and Southern Hemispheres (Pedro et al., 2018). In Antarctic ice cores, the corresponding Antarctic Isotopic Maxima (AIM) are characteristic warm events that are more gradual and of smaller amplitude than the Greenland events (EPICA community members, 2006). The AIMs are believed to be related to the DO events through the so-called bipolar seesaw mechanism (Bender et al., 1994; Stocker and Johnsen, 2003), but the detailed mechanism is a matter of debate (Landais et al., 2015; Pedro et al., 2018). Knowledge of the exact phasing of climate in the two hemispheres is crucial for deciphering the driving mechanism of the abrupt climate variability of the last glacial period and the climatic teleconnection patterns that connect the two hemispheres.

Three different techniques have been applied to progressively improve the synchronization of Greenland and Antarctic ice cores: globally well-mixed atmospheric gases – in particular the methane concentration (Blunier et al., 1998; Lemieux-Dudon et al., 2010; Rhodes et al., 2015; WAIS Divide Project Members, 2015) and the isotopic compositions of O₂ – δ¹⁸O_atm (Bender et al., 1994; Capron et al., 2010), cosmogenic isotopes such as ¹⁰Be (Raisbeck et al., 2017; Steinhilber et al., 2012), and the identification of large volcanic eruptions with bipolar sulfate deposition (Sigl et al., 2013; Svensson et al., 2013). A strength of the bipolar methane matching approach is that atmospheric methane concentrations change almost in phase with abrupt Greenland climate change allowing for those events to be synchronized in ice cores. A weakness of the bipolar gas matching approach is the dependency on a precise determination of the so-called Δage that refers to the offset in age between the ice and the air enclosed in an ice core at a given depth (Blunier et al., 2007; Schwander and Stauffer, 1984). Modelling past Δage requires an understanding of the physical processes taking place in the firn as well as knowledge or assumptions about past accumulation and temperature variations, introducing substantial age uncertainties associated with the synchronization.

Cosmogenic isotope production rates are modulated by the Earth's magnetic field and by solar variability, and they therefore carry a global signal that is shared by Greenland and Antarctic cosmogenic ice core records. Bipolar ice core synchronization using cosmogenic isotopes has mostly been done in the Holocene (Mekhaldi et al., 2015; Sigl et al., 2015; Steinhilber et al., 2012) and around the geomagnetic Laschamp event that occurred some 41 kyr ago (Raisbeck et al., 2017, 2007). Furthermore, the ice core cosmogenic signal enables the comparison with ¹⁴C records of other archives, such as dendrochronologies (Adolphi and Muscheler, 2016; Sigl et al., 2016) and stalagmites (Adolphi et al., 2018). Weaknesses of this technique include the sparsity of significant events, climatic influences on radionuclide transport and deposition masking the cosmogenic signal, and the very costly and time-consuming analyses that limit the possibility of obtaining continuous high-resolution records. Furthermore, archive noise in the ice core records hampers unambiguous peak detection and synchronization.

This study focuses on volcanic bipolar synchronization of ice cores in the second half of the last glacial period (12–60 ka). The volcanic record of the last glacial period in Greenland ice cores includes more than 100 confirmed Icelandic and high-latitude eruptions that have left predominantly cryptotephra (invisible to the naked eye) deposits in the ice (Abbott and Davies, 2012; Bourne et al., 2015; Cook et al., 2020); and there are presumably many more Icelandic eruptions that have not been identified as such. In addition to those, a large number of more distant eruptions have left an acidity signature in the ice cores but no tephra (Zielinski et al., 1997). In fact, during the last glacial only tephra from mid- and high-latitude eruptions has been identified in Greenland, whereas, to date, there is no evidence of tropical, low-latitude, or even continental European tephra in Greenland. Whether the lower-latitude tephra never makes it to Greenland or whether it is too small to be identified by conventional optical microscopy techniques and thus masked by more abundant, similar-sized background dust of continental origin is an open question.

In Antarctica, there are many visible tephra layers of Antarctic origin as well as a large number of acidity spikes associated with more distant eruptions (Narcisi et al., 2017; Severi et al., 2007). It has been proposed that tephra of tropical origin is present in the WAIS Divide ice core, but the evidence is solely based on dust size distributions and not on geochemical fingerprinting (Koffman et al., 2013). A pioneering study has suggested to have identified a bipolar tephra at around 1257 CE (Palais et al., 1992), and there is recent support for that conclusion suggesting that the source is the Indonesian Samalas volcano (Lavigne et al., 2013).

Although the ice core records lack bipolar tephra layers they do hold evidence of volcanic eruptions that are powerful enough to leave sulfuric acid in the stratosphere, from where it may be distributed to both Greenland and Antarctica. More than 80 such events have been identified for the last 2500 years (Sigl et al., 2015). For the earlier part of the Holocene, which has been less intensely studied, some 75 bipolar events have been found (Veres et al., 2013), and from around the time of the Indonesian Toba eruption occurring in Sumatra some 74 kyr ago, a handful of bipolar volcanic events have been identified (Svensson et al., 2013). Recently, there has been progress in identifying stratospheric volcanic peaks in ice cores based on their sulfur isotopic fingerprints (Burke et al., 2019; Gautier et al., 2019). In this work, we expand the bipolar volcanic matching approach systematically throughout the 12–60 ka time interval.

2 Methods

The approach taken to synchronize Greenland and Antarctic ice cores is to identify large volcanic eruptions with a bipolar acidity or sulfur or sulfate signature. Because such events generally do not leave tephra in ice cores from both polar regions, individual eruptions cannot be matched geochemically between the two hemispheres, as there is no way to verify that they have the same source. What can be matched up, however, are sequences of eruptions that show the same relative timing in both Greenland and Antarctica. To determine the time interval between eruptions, and thereby the relative timing of events, annual layer counting is carried out over the volcanic sequence in both Greenland and Antarctic ice cores. When identical volcanic peak patterns are identified in north and south, it is seen as a strong indication of a bipolar link. There is, however, always a risk of making an incorrect link, because an assumed volcanic sequence could consist of regional (non-bipolar) eruptions with coincidental similar temporal spacing.

The peak heights of the recorded eruption intensities in a bipolar volcanic ice core sequence cannot be expected to be similar at the two poles, because the strength of the recorded signal depends on the geographical location of the eruption, the atmospheric circulation at the time of the eruption, and the variability in deposition of acids at the ice coring sites (Gautier et al., 2016). Furthermore, the annual layer counting comes with an uncertainty that adds to the possibility of making an incorrect bipolar match (Rasmussen et al., 2006). On the other hand, the bipolar timing of Greenland and Antarctic ice core records is already well constrained by existing gas- and ¹⁰Be-based bipolar synchronizations. At the onset of each GI, the uncertainty in the bipolar methane matching between the Greenland NGRIP and the Antarctic WDC ice cores is around a century (WAIS Divide Project Members, 2015), which constrains the time windows for matching of bipolar volcanic sequences. Similarly, the cosmogenic isotope link around the time of the Laschamp event firmly constrains the volcanic matching in that time period (Raisbeck et al., 2017). Therefore, the risk of making false bipolar volcanic matches is strongly reduced by the existing bipolar synchronization.

We perform annual layer counting in sections of the Greenland NGRIP (North Greenland Ice Core Project members, 2004) and the Antarctic EPICA (European Project for Ice Coring in Antarctica) Dronning Maud Land (EDML) (EPICA community members, 2006) ice cores using high-resolution records of chemical impurities (Bigler, 2004; Ruth et al., 2008), dust (Ruth et al., 2003; Wegner et al., 2015), and visual grey-scale intensity (Faria et al., 2018; Svensson et al., 2005). The approach is the same as that applied for the glacial section of the Greenland Ice core Chronology 2005 (GICC05) (Andersen et al., 2006; Svensson et al., 2008). For this study, most sections of the NGRIP ice core have been recounted, and the Greenland timescale has been slightly modified as the bipolar matching allows for obtaining an improved precision from annual counting in both NGRIP and EDML. For most of the glacial period, the EDML ice core has not previously been layer-counted, but for the investigated time interval the annual layer thicknesses are comparable to those of NGRIP (Veres et al., 2013) and layer counting can be done in a similar way. Examples of annual layer counting in NGRIP and EDML across four intervals applied to match up patterns of bipolar volcanic eruptions are shown in Figs. S17–S20 in the Supplement. The WAIS Divide ice core chronology 2014 (WD2014) timescale has been applied for guidance, and for most of the intervals within the layer-counted section of WD2014 there is agreement within error estimates between the EDML and WDC interval durations. The bipolar layer counting is not continuous but is focused on periods of abrupt climate variability or high volcanic activity. In order to allow for comparison to published records, ages in all tables and figures have been converted to GICC05 ages using the year 2000 CE as a datum (referred to as “b2k”). We note that ages published on the Antarctic Ice Core Chronology 2012 (AICC2012) or the WD2014 timescale use the year 1950 CE as a datum, and so ages reported relative to b2k will be 50 years greater than those in AICC2012 and WD2014.

In order to obtain a robust identification of volcanic sequences in the ice cores, all available acidity records from the Greenland GRIP (Wolff et al., 1997), GISP2 (Mayewski et al., 1997; Taylor et al., 1997), NGRIP (Bigler, 2004), and NEEM (Schüpbach et al., 2018) ice cores have been included. For Antarctica, records from the EDML ice core from the Atlantic sector (EPICA community members, 2006), the EDC core from the East Antarctic Plateau (EPICA community members, 2004), and the West Antarctic Ice Sheet Divide (WDC) (Fudge et al., 2016; Sigl et al., 2016; WAIS Divide Project Members, 2013) ice core are used. Records of sulfur, sulfate, chloride, and electrical conductivity measurements (ECM) of the ice (Hammer et al., 1980), dielectric profiling (DEP) (Moore et al., 1989; Wilhelms et al., 1998), and the liquid conductivity of meltwater are also employed as good indicators of volcanic signals in ice cores. With the inclusion of those records, the ability of distinguishing large global volcanic events from more regional eruptions is improved, in particular for Antarctica.

The Greenland ice cores used here previously have been synchronized by volcanic events (Rasmussen et al., 2013; Seierstad et al., 2014). Likewise, the Antarctic ice cores have been linked internally by volcanic matching (Buizert et al., 2018; Ruth et al., 2007). In addition to the published volcanic match points made for Antarctica, some 25 additional Antarctic match points have been identified in the present study to strengthen the synchronization in the neighbourhood of Greenland abrupt climate change events. The non-bipolar or “local” volcanic matching applied here is in agreement with the published synchronizations for Greenland and Antarctica, respectively.

To investigate the bipolar climate signal, we employ stable water isotopes (δ¹⁸O) from the GRIP (Dansgaard et al., 1993; Johnsen et al., 2001), GISP2 (Grootes et al., 1993; Stuiver and Grootes, 2000), NGRIP (Gkinis et al., 2014; North Greenland Ice Core Project members, 2004), and NEEM (Schüpbach et al., 2018) ice cores as well as δ¹⁸O and deuterium excess from the EDML (EPICA community members, 2006; Stenni et al., 2010b), EDC (EPICA community members, 2004; Stenni et al., 2010b), WDC (Markle et al., 2016; WAIS Divide Project Members, 2013), Dome Fuji (DF) (Kawamura et al., 2007; Watanabe et al., 2003), and Talos Dome (TAL) (Landais et al., 2015; Stenni et al., 2010a) ice cores. The sources of the employed datasets are listed in Table S1.

3 Results

The bipolar volcanic match points identified in the 12–60 ka interval are shown in Fig. 1 and listed in Table S2. Of the 87 bipolar match points listed, five are previously published cosmogenic match points associated with the Laschamp geomagnetic excursion occurring in the 40.5–42.0 ka interval (Raisbeck et al., 2017). For the interval 16.5–24.5 ka, roughly corresponding to the Last Glacial Maximum (LGM), the ice cores are notoriously difficult to match up, and no bipolar match points are reported. We note that most of the identified bipolar match points fall within Greenland interstadial periods and rather few are located in stadials. The main reasons for this are the elevated dust concentrations in the colder periods that mute the ice conductivity signal as well as the elevated sulfuric background signal of colder periods that obscures the volcanic signal of the ice (Seierstad et al., 2014). Besides this, precise annual layer counting is also more difficult in the colder periods where accumulation is lower (Andersen et al., 2006). Wind-scouring is another factor affecting low-accumulation Antarctic sites particularly during colder periods.

Figure 1 Greenland (NGRIP) and Antarctic (EDML, WDC, and EDC) climate records (δ¹⁸O) throughout the 10–60 ka time period based on volcanic matching. Ages are on the GICC05 timescale relative to the year 2000 CE (“b2k”). Grey vertical lines show the position of bipolar volcanic match points identified in this study (Table S2) together with five match points based on cosmogenic isotopes around 41 ka (Raisbeck et al., 2017). Blue-shaded intervals indicate the Greenland Interstadial (GI) periods according to the definition of Rasmussen et al. (2014). The bipolar synchronization for the 16.5–24.5 ka interval is tentative as there are no bipolar match points in that interval.

All of the bipolar volcanic match points are identified in the Greenland NGRIP and the Antarctic EDML ice cores, and most of them are also identified in the Antarctic WDC and EDC ice cores. About half of the bipolar match points are also identified in the Greenland GRIP, GISP2, and NEEM ice cores, but the lower resolution of available sulfate and conductivity records for those cores is often insufficient for precise identification of weaker volcanic events. The Greenland ice cores are however internally synchronized throughout the last glacial period (Rasmussen et al., 2013; Seierstad et al., 2014) by northern hemispheric eruptions, typically Icelandic in origin, that leave a stronger fingerprint in Greenland. Therefore, for most of the bipolar match points all of the Greenland ice cores are precisely matched by interpolation between Greenland match points.

Table 1 Bipolar onsets and terminations. Depths of the Greenland interstadial onsets (a) and terminations (b) in the NGRIP, NEEM, GRIP, GISP2, EDML, EDC, and WDC ice cores based on volcanic matching. The NGRIP onsets are defined as the midpoints of the δ¹⁸O transitions and are identical to those applied in WAIS Divide Project Members (2015) except for the onset of GI-1. The Greenland match points are from Seierstad et al. (2014), Antarctic match points are from Buizert et al. (2018), and the bipolar match points are from Table S2. Corresponding GICC05 ages are provided with reference to year 2000 CE. “Distance” refers to the temporal distance to the nearest bipolar match point in Table S2 with negative values signifying the match point occurring before the onset or termination. The relative uncertainty of the bipolar matching is stated as 10 % of the “distance” assuming a 5 % maximum counting error of the annual layer counting in both Greenland and Antarctica. “YD-PB” refers to the Younger Dryas–Preboreal transition (the onset of the Holocene) and “BA-YD” refers to the Bølling–Allerød–Younger Dryas transition (the onset of GS-1). All EDC depths are on the EDC99 depth scale.

The bipolar match points are unevenly distributed over the 12–60 ka interval, but bipolar volcanic events have been identified within a range of 500 years of all major onsets and terminations of Greenland interstadials (GIs) with the exception of GI-2 (Fig. 1). All of the Greenland abrupt climate change events in that period (except for GI-2) are thus synchronized with the Antarctic climate at high precision. The derived depths and ages for the onsets and terminations of GI events are shown in Table 1. Based on the 5 % average counting uncertainty during the last glacial period (Svensson et al., 2008), the relative uncertainty of the bipolar linking related to the GI events is taken as 10 % of the distance to the nearest bipolar match point (Table 1). On average, this relative uncertainty is less than 15 years and reaches a maximum of 50 years. The definition of GI onsets and terminations applied in this study are the midpoints of the NGRIP isotopic transitions as identified in WAIS Divide Project Members (2015), except for the onset of GI-1 (the Bølling–Allerød), which is taken from Steffensen et al. (2008).

In the following sections, we provide examples of the bipolar synchronization from selected time intervals. In the Supplement detailed figures are provided for all the DO events (Figs. S1–S14).

3.1 The termination of GI-1/onset of the Younger Dryas

The onset of the Younger Dryas (GS-1) is synchronized between the two hemispheres by four large acidity spikes clustered around 13 ka and spanning 110 years (Fig. 2). The two outermost spikes are most significant, but all four spikes are present in all investigated cores. All four volcanic eruptions are interpreted as bipolar and they are therefore most likely associated with low-latitude eruptions. This is in conflict with the hypothesis that one of them should be related to the German Laacher See eruption (Baldini et al., 2018) that is believed to have a primarily northern hemispheric fingerprint (Graf and Timmreck, 2001). A tephra layer in the NGRIP ice core occurring close to the oldest of the four spikes has previously been tentatively associated with a Hekla eruption (Mortensen et al., 2005), but with the clear bipolar signature there is likely a temporal overlap between this Icelandic and an additional lower-latitude eruption. We notice that the eruption associated with the very significant North Atlantic Vedde Ash layer (Lane et al., 2012; Mortensen et al., 2005) located at 12.17 ka in the Greenland ice cores (Rasmussen et al., 2006) has potentially left a weak acidic signal in Antarctica.

Figure 2 Bipolar volcanic synchronization of the investigated ice cores across the transition from GI-1/Bølling–Allerød (BA) to GS-1/Younger Dryas (YD) (a) and the onset of GI-8 (b). Grey vertical lines are bipolar volcanic match points (Table S2). The records have different units; some are uncalibrated and peak heights are not comparable on an absolute scale, which is the reason why no scales are provided.

Figure 3 Synchronized climate records of the investigated ice cores across the GS-1 onset (a) and the GI-8 onset (b) applying the volcanic synchronization shown in Fig. 2. The acidity records in the bottom of the figure are for reference and are also shown in Fig. 2. All other records are δ¹⁸O in per mille (see Fig. 1 for scales). Grey vertical lines are bipolar volcanic match points (Table S2).

Our bipolar volcanic linking allows for synchronizing the climate signal of the investigated cores (Fig. 3). The four Greenland cores show quite variable climate patterns for the termination of GI-1, making it difficult to define the duration of the transition, but they all have the most significant drop in isotopic values in the interval constrained by the two bipolar events at 12.75 and 12.92 ka, respectively.

It has been proposed that the Younger Dryas period/GS-1 was initiated by a cosmic impact, for which there is indirect and debated evidence in a large number of sites in the NH surrounding the North Atlantic region (Kennett et al., 2015). A very significant Platinum (Pt) spike has been identified in the GISP2 ice core (Petaev et al., 2013) and at several North American sites (Moore et al., 2017) that potentially originate from the same impact event. The Pt spike occurs about 45 years after the volcanic quadruplet, i.e. after the Greenland cooling has initiated but before it has reached its minimum (Fig. S1b). The hypothesis of the YD initiation by a cosmic impact is currently debated (Holliday et al., 2020).

The onset of Greenland Interstadial 1/Bølling–Allerød

The very steep Greenland onset of the GI-1/Bølling–Allerød period is preceded by a 1.8 kyr long period of strong global volcanic activity (Fig. S2a). The bipolar phasing is well constrained by several significant eruptions leading up to the onset, and the bipolar volcanic matching pattern is easily recognized. In agreement with Steffensen et al. (2008), we note that NGRIP appears to be the Greenland ice core with the steepest δ¹⁸O transition at the GI-1 onset (Fig. S2b). The Antarctic ice cores all peak close to 200 years after the Greenland mid-transition in agreement with the methane matching of NGRIP and WDC (WAIS Divide Project Members, 2015). The very strong volcanic double spike in NGRIP close to 15.68 ka is associated with tephra from the explosive caldera-forming Towada-H eruption (Bourne et al., 2016) located in present-day Japan close to 40^∘ N. This eruption appears to have no significant Antarctic imprint.

3.2 Greenland Stadial 3 (GS-3)

In the late GS-3, at 24.67 ka, there is a characteristic volcanic triplet spike that constitutes a strong bipolar link (Fig. S3a). The three spikes are separated by 20±1 and then 10±1 years, making the match point unique (Fig. S3a inset). At 25.46 ka (b2k GICC05 age) we hypothesize a record of traces from the Oruanui eruption from the Taupo volcano in present-day New Zealand in the Greenland record. Tephra of this eruption has been previously identified and dated to 25.37 ka (b2k WD2014 age) in the WDC ice core (Dunbar et al., 2017). There appears to be a major Greenland acidity spike associated with this eruption despite its latitude being close to 40^∘ S. Unfortunately, there are no adjacent bipolar eruptions within several hundreds of years, making the bipolar Oruanui link somewhat uncertain. The nearest pair of bipolar eruptions is found at 25.76 and 25.94 ka, leaving enough room in the layer counting uncertainty for the Greenland acidity spike to potentially be a “false match” offset by up to 30 years from the Antarctic Oruanui spike. The link needs to be investigated by the bipolar sulfur isotopes method to rule out a coeval local source for the NGRIP event (Burke et al., 2019). Besides being relevant for comparing Greenland and Antarctic ice core timescales (Sigl et al., 2016), the Oruanui eruption constitutes an important comparison point for ¹⁴C ages and ice core chronologies as tephra from the eruption is widely distributed (Muscheler et al., 2020).

3.3 The onset of Greenland Interstadial 8

The very prominent onset of GI-8 is associated with four significant bipolar eruptions within a 400-year period (Fig. 2). The eruption occurring at 38.13 ka, close to a century after the onset, shows a very significant signal in all acidity proxies of all investigated ice cores, and it appears to be one of the largest eruptions of the last glacial period. We thus see it as a potential candidate for the H1 horizon identified in Antarctica by radio-echo sounding (Winter et al., 2019). Other prominent bipolar events are situated at 37.97, 38.23, and 38.37 ka. The 38.23 ka event occurs right at the initiation of the GI-8 onset.

In the climate records across the GI-8 onset, the Greenland records behave quite similarly, whereas the Antarctic records show rather distinct patterns (Fig. 3). EDC expresses a prominent warming coinciding with the Greenland warming; WDC also shows a warming although much less significant. Both EDC and WDC exhibit a cooling trend initiating some decades after the Greenland onset. In contrast, EDML shows a century-long cooling period starting right at the Greenland onset. Being rather noisy, it is hard to separate signal from noise in the Antarctic records based on just one warming event. However, the stacking exercise across several warming events discussed below reveals that the pattern expressed by the Antarctic records at the GI-8 onset is an archetypical expression of an Antarctic response to a major GI onset.

3.4 Greenland interstadials 9 and 10

For the period 40–43 ka that covers GI-9 and GI-10, the bipolar matching is very well constrained by 16 bipolar match points (Fig. S8a), five of which are independent cosmogenic match points related to the Laschamp geomagnetic excursion (Raisbeck et al., 2017). The volcanic match presented here is in agreement within uncertainties with the cosmogenic matching, and it replaces the existing bipolar volcanic matching in this region (Svensson et al., 2013) that has been shifted by some 30 years.

The onsets of GI-9 and GI-10 provide examples of less prominent GI events where the Antarctic isotopic response pattern is different from that of the larger events, such as GI-8 and GI-12 (Figs. 2 and S10b). For the GI-9 onset, it is practically impossible to distinguish a bipolar seesaw response in the adjacent periods in the Antarctic climate records. For the GI-10 onset, EDC has a small spike and EDML has a dip, but none of them stand out from the background. WDC appears to enter GI-10 without any climatic response. For the smaller or weaker GI events the response pattern in Antarctica is likewise smaller, making it harder to identify it within the isotopic background variability (Fig. 1).

4 Bipolar phasing of abrupt climate change

The characteristics of the climate record of the individual GI onsets may vary from event to event, from ice core to ice core, and from proxy to proxy both in Greenland and in Antarctica. In Greenland, the main pattern associated with a GI onset is similar for all deep ice cores, but transition durations and the relative phasing of individual parameters, such as water isotopes and impurity concentrations, vary between events and to a lesser degree between cores (Erhardt et al., 2019; Rasmussen et al., 2014; Steffensen et al., 2008). In Antarctica, the investigated cores are located further apart, they are exposed to the climatic influence from different ocean basins, and they do in general show greater variability in relation to the Greenland GI onsets than is the case for the more closely located Greenland coring sites. In addition to the noise in water isotope records caused by wind-driven redeposition at low-accumulation sites, millennial-scale climate variability in Antarctica has a profoundly lower signal-to-noise ratio than that in Greenland, contributing to the difficulty of interpreting individual events.

Keeping in mind this variability among individual GI events, we find it useful to extract a general pattern across the GI onsets and terminations by aligning the individual events and stacking (or compositing) them, as has been done using bipolar methane synchronization (Buizert et al., 2018; WAIS Divide Project Members, 2015). The stacking provides us with an overall phasing relation that may be helpful in unravelling the governing mechanisms of the abrupt climate change, for example by comparison to model experiments that typically do not capture the details of individual events (Buizert et al., 2018). We stress that the associated underlying assumption of the stacking approach, that the complete temporal progression of all events is the result of the same underlying process, may not be fully justified.

Figure 4 Stacks of isotopic records across GI onsets (a, c, e) and terminations (b, d, f) for the events listed in Table 1a and b, applying the bipolar volcanic synchronization. The time “t=0” refers to the midpoint of the NGRIP onset for each GI event as defined in Table 1. (a, b) Stack of NGRIP δ¹⁸O (blue; left axis) and WDC CH₄ (green; right axis). (c, d) Stack of Antarctic δ¹⁸O at the indicated locations and the average curve (mean). (e, f) Stack of Antarctic d_ln at the indicated locations and the mean. The figure is modified from Buizert et al. (2018). See Fig. S15 for Antarctic core site locations.

In Fig. 4, we show the δ¹⁸O records of NGRIP together with five Antarctic ice cores stacked across all of the GI onsets and terminations in the 12–60 ka interval (except for GI-2), centred at the Greenland transition midpoint as defined in Table 1. Besides the EDML, EDC, and WDC ice cores applied for the bipolar synchronization in this study, we expand the geographical coverage by including the Antarctic Dome Fuji (DF)(Kawamura et al., 2007) and Talos Dome (TAL) (Stenni et al., 2010a) records applying the existing Antarctic volcanic synchronization (Buizert et al., 2018; Fujita et al., 2015; Severi et al., 2012). Figure 5 shows the stacking of the five Antarctic ice cores from Fig. 4, thus resulting in a stack of 21 events in five Antarctic ice cores totaling 105 events. We note that the Greenland onsets are aligned according to the midpoint of the warming transition (set to t=0), implying that the initiation of the Greenland event occurs earlier than the alignment point. A recent study suggests that the abrupt NGRIP δ¹⁸O and Calcium (Ca) transition onsets on average precede the δ¹⁸O transition midpoint by 25±7 and 33±15 years, respectively (Erhardt et al., 2019); in Figs. 4 and 5 this places the Greenland δ¹⁸O and Ca event onsets at $t=-\mathrm{25}$ and $t=-\mathrm{33}$ years, respectively.

Figure 5 (a) Stack of Greenland NGRIP isotopes and CH₄ for the onsets of GI events listed in Table 1a. (b) Antarctic five-core mean δ¹⁸O (stack of 5×21 warming events). Black curve is the same as “mean” in Fig. 4, centre; grey curve applies the bipolar methane synchronization of WAIS Divide Project Members (2015). (c) Antarctic five-core mean deuterium excess (d_ln). Black curve is the same as “mean” in Fig. 4e, f; grey curve applies the bipolar methane synchronization of WAIS Divide Project Members (2015). Orange curves are fitting functions using the change-point analysis applied in WAIS Divide Project Members (2015). The stated 2σ uncertainty estimates are obtained from a Monte Carlo sampling (see main text). When the same uncertainty estimate is made for the GI terminations (Fig. 4) the δ¹⁸O timing is at 101±29 years, the d_ln transition onset is at $-\mathrm{59}\pm \mathrm{58}$ years, and the d_ln transition end is at 95±34 years; the duration between the d_ln onset and the δ¹⁸O change point is 160±65 years. The earlier d_ln onset for the GI termination probably reflects the GI terminations being generally more gradual than the GI onsets.

For Antarctica, the stacked EDC, WDC, and EDML records show distinct δ¹⁸O patterns similar to those identified for the onset of GI-8 (Fig. 3). EDC, WDC, and TAL show a peak of accelerated warming, the onset of which is synchronous with Greenland warming and that lasts for close to a century. DF likewise shows an accelerated warming, albeit somewhat later than the aforementioned cores. This direct Antarctic warming response to the Greenland warming is likely to be associated with fast atmospheric changes on a global scale (Markle et al., 2016). In particular, it has been proposed that a northward shift in the SH westerlies in response to NH warming (Lee et al., 2011; Pedro et al., 2018) may drive a warming anomaly in most of the Antarctic continent through enhanced zonal heat transport in the atmosphere (Buizert et al., 2018; Marshall and Thompson, 2016).

Another process that is likely to contribute to the alignment of the water isotope records at the GI onsets is the local cycle of sublimation–condensation in summer on the Greenland and Antarctic ice sheets that is currently under investigation (Kopec et al., 2019; Pang et al., 2019). Sublimation affects the isotope concentration and the deuterium excess of snow through kinetic fractionation. Snow sublimation requires large amounts of energy, and it is controlled by the relative humidity, which in turn is linked to the large-scale atmospheric circulation. Sublimation effects are poorly constrained on the East Antarctic Plateau.

EDML, however, shows an immediate cooling response that is distinct among the cores investigated (Figs. 4 and S15), perhaps reflecting regional effects such as wind-driven changes to the Weddell Sea stratification, gyre circulation, sea-ice extent, or polynya activity.

Based on the volcanic bipolar synchronization, the general Antarctic response time (Fig. 5) to a Greenland warming event is shorter than that obtained from bipolar methane linking (WAIS Divide Project Members, 2015). Instead of the 200-year lag found in the methane-based synchronization, we find an average response time of 122±24 years (2σ uncertainty) using the same fitting routine as used in WAIS Divide Project Members (2015) (see discussion of uncertainty estimates below). This difference is mostly due to uncertainty in the WDC Δage calculation; the new synchronization suggests that the glacial Δage was too small by around 70 years on average (Fig. S16).

Besides δ¹⁸O, we also stack records of Antarctic deuterium excess using the logarithmic definition (d_ln) introduced by Uemura et al. (2012). Markle et al. (2016) showed that in the WAIS Divide ice core, d_ln abruptly increases in synchrony with the onset of GIs; at the termination of GIs, d_ln abruptly decreases. Markle et al. (2016) used a climate model simulation with moisture tagging to show that this relationship could be explained by north–south shifts in the location of moisture sources associated with changes in the shifts of the Southern Hemisphere (SH) subpolar jet, and westerly winds. This is consistent with work of Schmidt et al. (2007), who had previously shown with climate model simulations that the deuterium excess should be inversely correlated with the Southern Annular Mode (SAM) index. Masson-Delmotte et al. (2010) made a similar argument on the basis of the EDC core, but without sufficient dating precision to demonstrate the close relationship found by Markle et al. (2016). These findings were later extended to multiple Antarctic sites by Buizert et al. (2018).

Stacks using a methane-based synchronization show a d_ln transition that takes ∼220 years, followed by a broad peak (Fig. 5); in contrast, our volcanic synchronization suggests a shorter transition (152±37 years) and a much sharper d_ln transition. Any chronological errors in the bipolar synchronization will misalign the events being composited, thereby broadening the climatic features in the stacked record. The fact that our volcanic synchronization yields sharper features is thus indirect evidence that it is more accurate than existing gas-based synchronizations. The duration of the d_ln transition we observe in our stack is still an upper bound on its true duration, given that our event alignment includes uncertainties due to annual layer counting to the nearest volcanic tie point as well as potentially incorrectly identified volcanic tie points. It was suggested that the gradual d_ln trends before and after the transition follow the gradual source-water sea-surface–temperature trends of the SH via the bipolar seesaw (Markle et al., 2016).

Our volcanic bipolar synchronization also shifts the onset of the d_ln transition towards older ages, placing it at $t=-\mathrm{30}\pm \mathrm{29}$ years (2σ uncertainty) relative to the Greenland δ¹⁸O transition midpoint. Such an early onset may seem surprising but is actually in very good agreement with other Greenland proxies that suggest that low-latitude changes precede the Greenland δ¹⁸O signals. In particular, changes in Greenland dust or Ca concentrations appear to lead the Greenland δ¹⁸O by a decade at the onset of the transitions (Erhardt et al., 2019), which has been attributed to early changes in the Inter-Tropical Convergence Zone (ITCZ) position and atmospheric circulation (Steffensen et al., 2008). Meridional shifts in the SH eddy-driven jet and westerlies are suggested to be dynamically linked to the ITCZ position (Ceppi et al., 2013). The onset of the Greenland Ca transition (presumably reflecting NH atmospheric circulation shifts) precedes the Greenland δ¹⁸O transition midpoint (which we set as t=0) by 33±15 years on average (Erhardt et al., 2019), in good agreement with the 30±29 years we find for the onset of SH atmospheric circulation changes.

The Antarctic d_ln and δ¹⁸O signals (Fig. 5) are recorded in the same physical ice cores, and therefore the uncertainty in their relative phasing is small and only related to the stacking of the Antarctic cores and the change-point determinations. Errors in the bipolar synchronization will blur the abruptness of their transitions but should not alter their relative phasing; by contrast, the phasing relative to Greenland proxies is very sensitive to bipolar synchronization errors. The 152±37-year duration between the onset of the d_ln response and the breakpoint in the δ¹⁸O curve therefore represents a robust estimate of the climatic lag of the mean Antarctic temperature response behind the first atmospheric manifestation of the GI event in the Southern Hemisphere high latitudes.

At the terminations of GI events, the stacked NGRIP δ¹⁸O shows a less prominent but still sharp transition over a ∼100-year interval (Fig. 4). For δ¹⁸O, all of the stacked Antarctic cores reach a minimum in the interval 100–150 years following the Greenland termination, with the strongest response seen for TAL. For d_ln, EDC and DF show a significant response related to the Greenland terminations, whereas the other Antarctic cores express a less coherent signal. Note that EDML has its coldest temperatures near t=0, suggesting again a fast response at this site of opposite sign to that in the DO warming case.

Uncertainty estimate of stacked records

To estimate the uncertainty in the change-point analysis of the stacked records (Fig. 5) we use a Monte Carlo scheme with 1000 iterations (WAIS Divide Project Members, 2015). In each iteration, the alignment of the individual events is shifted randomly prior to the stacking, and the change point is identified in the new stack using an automated algorithm. The applied time shifts are randomly drawn from normal distribution with widths corresponding to the following event-specific uncertainties: (1) the uncertainty in the NGRIP event midpoint detection (WAIS Divide Project Members, 2015); (2) the uncertainty in the Antarctic layer count from the bipolar eruption to the event (Table 1a); (3) the uncertainty in the Antarctic volcanic synchronization (Buizert et al., 2018).  In each iteration, the user-specified parameters of the fitting algorithm (such as the time interval used in the fitting) are likewise perturbed randomly. The stated 2σ uncertainty values therefore reflect uncertainty in the bipolar synchronization, stacking procedure, and change-point detection.

The Antarctic delay times and uncertainties identified for δ¹⁸O and d_ln, respectively, are valid for the stacked (averaged) transitions (Fig. 5) but are not representative of the variation among individual transitions. When performing an event-by-event fitting of the Antarctic five-core average we find a much greater range of delay times. For the δ¹⁸O change point, the mean and standard deviation of the individual event timings is $t=\mathrm{138}\pm \mathrm{89}$ years; for the d_ln it is $-\mathrm{6}\pm \mathrm{78}$ and 116±80 years for transition beginning and end, respectively. This larger variability reflects both differences in timing between individual events as well as the much smaller signal-to-noise ratio when fitting individual events.

5 Conclusions and outlook

Overall, our new bipolar volcanic synchronization confirms the centennial-scale delay of the mean Antarctic bipolar seesaw temperature response behind abrupt Greenland DO variability (WAIS Divide Project Members, 2015); however, the improved age control offered by volcanic synchronization significantly reduces the estimated duration of this delay compared to previous work based on CH₄ synchronization. Our reduced estimates are more in line with, but still larger than, results from climate model simulations (Pedro et al., 2018; Vettoretti and Peltier, 2015) that typically give an oceanic response on multi-decadal timescales.

WAIS Divide Project Members (2015) interpreted the 208±96-year delay they observed as characteristic of an oceanic teleconnection. The reduced delay timescale we infer here (122±24 years) is still consistent with the original interpretation. However, the new observations urge some caution. Despite our best efforts, the new bipolar volcanic framework likely contains some incorrect matches, and as the bipolar synchronization continues to be refined, the inferred Antarctic delay may be reduced further. The divergent climate response at various sites (Figs. 4, S15 and WAIS Divide Project Members, 2015) as well as the relatively gradual transition in d_ln over ∼145 years (suspiciously similar to the updated timescale Antarctic temperature delay presented here) perhaps suggest a more complex interplay of atmospheric and oceanic processes that is currently very poorly understood (see, e.g., Kostov et al., 2017). We suggest future work along two parallel lines of inquiry.

First, further refinement and confirmation of our bipolar synchronization is called for. Analysis of sulfur mass-independent isotopic fractionation (Burke et al., 2019) is needed for the proposed bipolar volcanic events. For low-latitude eruptions to show up in both the Arctic and Antarctic almost certainly requires the injection of materials into the stratosphere, which is reflected in Δ³³S. High-resolution records of ¹⁰Be can further refine bipolar matching (Adolphi et al., 2018) in critical intervals where the volcanic record is ambiguous. Second, climate modelling studies are needed to better understand the interaction between atmospheric and oceanic changes during the Dansgaard–Oeschger cycle. In particular the anomalous response of the EDML site during both Heinrich (Landais et al., 2015) and Dansgaard–Oeschger (Buizert et al., 2018) abrupt climate change calls for detailed investigation.

The volcanic bipolar synchronization has a wide range of potential applications that go beyond the objectives of this paper. These include the development of consistent bipolar ice core timescales, constraining ice core delta-gas ages, the investigation of impacts of volcanism on abrupt climate change, the quantification of the last glacial global volcanic eruption record, and the discussion of solar variability through the synchronization of cosmogenic isotopes. Furthermore, the precise bipolar synchronization should allow for an improved understanding of the mechanisms governing the glacial climate through a comparison to model studies and non-ice core records.