Precise interpolar phasing of abrupt climate change during the last ice age Authors: WAIS Divide Project Members This is a postprint of an article that originaly appeared in Nature on April 2015. WAIS Divide Project Members. "Precise interpolar phasing of abrupt climate change during the last ice age." Nature 520, no. 7549 (April 2015): 661-665. DOI:10.1038/nature14401 Made available through Montana State University’s ScholarWorks scholarworks.montana.edu Precise interpolar phasing of abrupt climate change during the last ice age WAIS Divide Project Members* The last glacial period exhibited abrupt Dansgaard–Oeschger cli- matic oscilations, evidence of which is preserved in a variety of Northern Hemisphere palaeoclimate archives1. Ice cores show that Antarctica cooled during the warm phases of the Greenland Dansgaard–Oeschger cycle and vice versa2,3, suggesting an inter- hemispheric redistribution of heat through a mechanism caled the bipolar seesaw4–6. Variations in the Atlantic meridional overturn- ing circulation (AMOC) strength are thought to have been import- ant, but much uncertainty remains regarding the dynamics and trigger of these abrupt events7–9. Key information is contained in the relative phasing of hemispheric climate variations, yet the large, poorly constrained diference between gas age and ice age and the relatively low resolution of methane records from Antarctic ice cores have so far precluded methane-based synchron- ization at the required sub-centennial precision2,3,10. Here we use a recently driled high-accumulation Antarctic ice core to show that, on average, abrupt Greenland warming leads the corresponding Antarctic cooling onset by 218692 years (2s) for Dansgaard– Oeschger events, including the Bøling event; Greenland cooling leads the corresponding onset of Antarctic warming by 208696 years. Our results demonstrate a north-to-south direc- tionality of the abrupt climatic signal, which is propagated to the Southern Hemisphere high latitudes by oceanic rather than atmo- spheric processes. The similar interpolar phasing of warming and cooling transitions suggests that the transfer time of the climatic signal is independent of the AMOC background state. Our find- ings confirm a central role for ocean circulation in the bipolar seesaw and provide clear criteria for assessing hypotheses and model simulations of Dansgaard–Oeschger dynamics. Net heat transport by the Atlantic branch of the global ocean over- turning circulation is northwards at al latitudes, resulting in a heat flux from the Southern Hemisphere (SH) to the Northern Hemisphere (NH)4. Variations in this flux act to redistribute heat between the hemispheres, a mechanism commonly invoked to explain abrupt sub-orbital climatic variability and the asynchronous coupling of Greenland and Antarctic temperature variations on these timescales2,6. Milennial-scale AMOC variability is corroborated by North Atlantic proxies for deep-water ventilation and provenance that suggest decreased North Atlantic deep water production and the intrusion of southern-sourced water masses during stadial (that is, North Atlantic cold) periods11,12. The oceanic instabilities are accompanied by shifting atmospheric transport paterns. Proxy data and climate models consistently indicate northward and southward migrations of the intertropical convergence zone in response to abrupt NH warm- ing and cooling, respectively13,14. Abrupt NH events may also induce changes in the strength and position of the SH mid-latitude westerlies. Strengthened and/or southward-shifted westerlies during NH stadials have the potential to warm the Southern Ocean and Antarctica by enhancing the wind-driven upweling of relatively warm circumpolar deep waters, providing a direct atmospheric pathway for the bipolar seesaw to operate15. Climate model simulations further suggest that atmospheric readjustment in response to decreased North Atlantic deep water formation can induce inhomogeneous temperature changes over the Antarctic continent16,17, possibly by means of wind-driven changes in sea-ice distribution8. Atmospheric teleconnections operate on seasonal to decadal time- scales because of the fast response time of the atmosphere8,18,19.By contrast, oceanic teleconnections can operate on a wide range of time- scales from decadal to multi-milennial, depending on the processes and ocean basins involved19–21. Climatic signals can be rapidly propa- gated from the North Atlantic to the South Atlantic via Kelvin waves5, but models suggest a more gradual (centennial timescale) propagation from the South Atlantic to the SH high latitudes as a result of the absence of a zonal topographic boundary at the latitudes of the Antarctic circumpolar current (ACC)19,21. The timescale on which newly formed North Atlantic deep water is exported to the Southern Ocean is similarly around several centuries22. The bipolar seesaw rela- tionship observed between ice core records potentialy bears the imprint of both atmospheric and oceanic teleconnections; precise con- straints on the interhemispheric phasing can help distinguish which mode dominates, and can also identify leads and lags10,19,21. Completed in December 2011, the West Antarctic Ice Sheet (WAIS) Divide ice core (WDC)17was driled and recovered to a depth of 3,405 m. Here we present results from the deep part of the core (Fig. 1), extending back to 68 kyr before 1950 (before present;BP). The WDC climatic records show no stratigraphic disturbances (such as large-scale folds) that are commonly encountered near the base of ice cores, a fact we atribute to basal melting at the site that removes old ice before such disturbances can develop and that decreases shear stress by lubricating the bed. The Greenland NGRIP core similarly has strong basal melting and an undisturbed stratigraphy1. WDC d18O of ice, a proxy for local condensation temperature, exhibits clear milennial-scale variability as observed across the Antarctic continent (Fig. 1c, d). We evaluate the phasing of WDC milennial variability relative tod18O of the Greenland NGRIP core (Fig. 1a). For each Greenland Dansgaard–Oeschger (DO) event we can identify a corres- ponding Antarctic Isotopic Maximum (AIM) event3, although AIM 9 is only very weakly expressed in WDC (we adopt the naming conven- tion whereby AIMxis concomitant with DOx). To investigate the phasing of the bipolar seesaw, we have synchro- nized WDC to the Greenland NGRIP core by means of globaly wel- mixed atmospheric methane (CH4; Fig. 1b)2,23. WDC alows CH4synchronization at unprecedented, sub-centennial precision as a result of the smal diference between gas age and ice age (Dage) and continuous, centimetre-scale resolution CH4record (ExtendedData Figs 1 and 2). Because uncertainties in the relative phasing of CH4and Greenland climate24are smaler than uncertainties inGreenlandDage, we synchronize WDC CH4directly to NGRIPd18O (rather than to NGRIP CH4)23. As pointed out by several authors6,21, a strong anti-correlation exists between Greenlandd18O and the rate of change (that is, the first time derivative) of Antarcticd18O. This relationship is consistent with a simple thermodynamic model of the bipolar seesaw in which Antarctic temperature variations are moderated by the large thermal mass of the *Lists of participants and their afiliations appear at the end of the paper. Southern Ocean6. In evaluating this relationship at WAIS Divide (Fig. 2a, blue curve), we find the strongest anti-correlation at a NH lead of about 170 years (blue dot in Fig. 2a). This timing is distinct from the zero-year lead commonly assumed in mathematical descriptions of the bipolar seesaw. The analysis is repeated using an ensemble of 43103alternative WDC chronologies (Methods), showing that the centennial-scale timing is a robust result (Fig. 2d) with an estimated uncertainty of 69 years (2suncertainty bounds are used throughout this work). Evaluation of the cross-correlations using WDC CH4instead of NGRIPd18O gives almost identical results (Fig. 2a, green curve), in line with the notion that CH4is a good proxy for Greenlandclimate2. Next we analyse the DO–AIM coupling in more detail. The bipolar seesaw theory suggests that during Greenland stadial periods heat accumulates in the SH, causing gradual warming of Antarctica; during interstadial periods Antarctica is cooling. The abrupt stadial–intersta- dial transitions in Greenland are accompanied by a breakpoint in the WDC d18O record, where the warming trend changes to a cooling trend. Although the timing of abrupt transitions can be established unambiguously in both CH4and Greenlandd18O records, the smalersignal-to-noise ratio of Antarcticd18O time series complicates break- point identification for individual AIM events. To overcome this lim- itation, we stack and average the AIM events to detect their shared climatic signal (Methods and Extended Data Figs 3 and 4). We align WDC records at the midpoint of the DO CH4transitions, which is setto lag the midpoint in the NGRIPd18O transition by 56638 years24. The relative timing of the WDCd18OandCH4curves is determined byDage, which is the largest source of uncertainty (Methods). We find that Antarctic cooling is delayed relative to abrupt NH warming by 218692 years on average (Fig. 2b, e); similarly, Antarctic warming lags NH cooling by 208696 years (Fig. 2c, f). The robustness of this result is demonstrated with a Monte Carlo analysis in which we randomly perturb the relative alignment of the individual AIM events, in combination with an ensemble of 43103 alternative WDC chronologies. The timing of the WDCd18O break- point (orange dots in Fig. 2b, c) is determined by using an automated fiting algorithm (Methods and Extended Data Fig. 5). Performing the analysis separately on a stack of only the major AIM events (AIM 4, 8, 12, 14 and 17), only the minor AIM events or stacks of eight randomly selected AIM events gives nearly identical phasing (Extended Data Fig. 6), excluding the possibility that the result is dominated by the timing of a few prominent events. We repeat the stacking procedure for the WDC sea-salt sodium record (Extended Data Fig. 7), and find that this tracer, which has been interpreted as a proxy for sea-ice production25, changes almost synchronously with WDCd18O—possibly reflecting a common forcing by Southern Ocean temperatures6. At the onset of the NH Bøling warm period (DO 1; 14.6 kyr) we find a north–south phasing comparable to that of the glacial period, with a 2566133-year lag of the Antarcticd18O response (Fig. 3 and Extended Data Table 1). The interpolar phasing at the onset and termination of the Younger Dryas stadial is ambiguous in our record (Extended Data Table 1), mostly because of the relatively gradual nature of the Younger Dryas onset and the presence of two local maxima in the WAIS-Dd18O record around the Younger Dryas– Holocene transition (possibly reflecting regional climate). A recent study using a CH4-synchronized stack of five near-coastal Antarcticcores spanning the deglaciation suggested a near synchrony of Antarctic breakpoints and Greenland transitions within a dating uncertainty of 200 years10. Our result overlaps with this range. However, the much smaler dating uncertainty at WDC and our use of the ful sequence of glacial AIM events shows that a centennial-scale Antarctic lag is a systematic feature of DO–AIM coupling throughout the glacial period. The lead of NH climate revealed at WDC (Fig. 2) provides clues about DO climate dynamics. The simplest interpretation is that the abrupt phases of the DO cycle are initiated in the NH, presumably in the North Atlantic. This is consistent with several proposed DO mechanisms, such as North Atlantic sea-ice dynamics9, freshwater 253035404550556065 WD2014 age (kyr BP) 253035404550556065 2 SIM3 SIM 4 SIM –46 –44 –42 –40 –38 NG RIP δ 18 O ( ‰) 350 400 450 500 550 600 650 WD C C H 4 (p .p. b.) –42 –40 −38 −36 WD C δ18 O ( ‰) –10 –8 –6 –4 AT S Δ T ( °C) DO: 2345.15.26789101112131415161718 a b c d Figure 1|Records of glacial abrupt milennial-scale climatic variability. a, Greenland NGRIPd18O record1on GICC0531.0063 chronology (Methods).b, WDC discrete CH4record on the WD2014 chronology, which isbased on layer counting (0–31.2 kyr) and CH4synchronization to NGRIP (31.2–68 kyr)23.c, WDCd18O record.d, Antarctic temperature stack (ATS)30in degrees Celsius relative to the present day on AICC1231.0063 chronology. DO/AIM events are indicated with orange vertical bars, numbered at the botom of the figure. forcing8or ice shelf colapse7. We acknowledge, however, that the abrupt North Atlantic events could be the response to a remote process not visible in the ice core records; therefore, although we cannot ascer- tain the location of the elusive DO ‘trigger’ (if such a concept is even meaningful in a highly coupled dynamical system), our results clearly indicate a north-to-south directionality of the abrupt phases of the bipolar seesaw signal. The centennial timescale of the NH lead demonstrates the dom- inance of oceanic processes in propagating NH temperature anomalies to the SH high latitudes. Any atmospheric teleconnection would be manifested within at most a few decades after the abrupt NH event18; on this timescale the WDCd18O stacks (Fig. 2b, c) show no discernible temperature response above the noise. We estimate the noise level by subtracting the fited curves from thed18O stacks; the remaining signal has a 2svariability of 0.14%, or about 0.18uC assuming an isotope sensitivity of 0.8%K21. We therefore estimate an upper bound of 0.18uC on an atmosphericaly induced Antarctic temperature res- ponse (for comparison, milennial Antarctic temperature variations are on the order of 1–2uC). A readjustment of atmospheric transport may induce spatialy inhomogeneous temperature changes8,16that may be (partly) responsible for the heterogeneity in expression, or shape, of AIM events across Antarctica3. We find that on average the DO cooling signal is transmited as fast to Antarctica as the DO warming signal is (our sensitivity study sug- gests a diference in propagation time of 10689 years). This implies that the north-to-south propagation time is independent of the AMOC background state; that is, it is independent of whether the AMOC is in the weak or strong overturning state. Modeling work suggests that the meridional propagation time across the ACC latitudes depends prim- arily on ACC strength21; our inference of similar propagation times may thus reflect a stability of the ACC on milennial timescales. There is also a conspicuous synchronicity between the phasing of the bipolar seesaw and the duration of the abrupt increase in CO2at 14.6 kyr(Fig. 3d); if the former does indeed reflect the timescale of the oceanic response, this may hold clues about the (unidentified) source of this CO2(ref. 26). Proxy records of North Atlantic ventilation and overturning strength during Marine Isotope Stage 3 typicaly show the most prom- inent excursions during Heinrich stadials11,12, periods of extreme cold in the North Atlantic associated with layers of ice-rafted debris in ocean sediments that represent times of massive delivery of –500–250 0 250 500 750 –0.8 –0.6 –0.4 –0.2 0 0.2 SH leads NH leads Lag (years) Cor rel ati on 0 200 400 600 Fre qu en cy 0 50 100 150 200 250 300 169 ± 69 NH lead time (years) –500 –250 0 250 500 750 0 0.5 1.0 1.5 2.0 Nor mali ze d si gn al Time (years) 0 10 20 30 40 Fre qu en cy (× 10 3 ) 0 100 200 300 218 ± 92 NH lead time (years) –250 0 250 500 750 0 100 200 300 400 208 ± 96 ba c e fd Figure 2|Interhemispheric phasing of the bipolar seesaw. a, Lagged correlation between NGRIPd18O and WDC d(d18O)/dt(blue), and between WDC CH4and d(d18O)/dt(green). Milennial-scale variability is isolated byusing a fourth-order Buterworth bandpass filter (500–10,000-year window); the CH4-synchronized part of the records is used (31.2–68 kyr).b,DO3–18stack of NGRIPd18O (blue), WDC CH4(green) and WDCd18O (orange withcurve fit), aligned at the midpoint of the DO warming signal. Events are averaged with their original amplitudes and normalized after stacking for convenience of visualization.c,Asinb, but for NH abrupt cooling events (that is, the interstadial terminations).d–f, Histograms of NH lead time associated witha–c, respectively, generated by binning solutions from the sensitivity study. The total number of solutions is 43103ind, and 43105ineand f. Distribution mean and 2sprobability bounds are listed in the panels. Shaded vertical yelow bars (upper panels) show NH lead time; the error bar represents 2sas defined for the lower panels. 1112131415 1112131415 Holoc.YDB−A / ACROD –44 –42 –40 –38 –36 NG RIP δ 18 O ( ‰) 400 500 600 700 800 WD C C H 4 (p .p. b.) –39 –38 –37 –36 –35 –34 WD C δ18 O ( ‰) 230 240 250 260 270 WD C C O 2 (p .p. m.) a b c d WD2014 age (kyr BP) Figure 3|Timing of the last deglaciation. a, NGRIPd18O on GICC05 chronology1.b, WDC CH4.c, WDCd18O with breakpoints as orange dots anderror bars showing the 2suncertainty bounds (Extended Data Table 1). d, WDC CO2data (dots) with spline fit (solid line)26. Period abbreviations: OD,Oldest Dryas; B–A, Bøling–Alerød; ACR, Antarctic Cold Reversal; YD, Younger Dryas; Holoc., Holocene. Vertical orange lines correspond to the midpoints of the WDC CH4transitions. NGRIP and WDC chronologies areboth based on annual-layer counting, and are fuly independent. debris-laden icebergs (Heinrich events). This suggests that there may be a diference between climate dynamics during Heinrich and non- Heinrich stadials. We find no anomalous NH lead time for DO transi- tions directly after or before Heinrich stadials (Extended Data Fig. 6), and thus no evidence that Heinrich stadials are unusual from the perspective of the oceanic teleconnections that dominate the bipolar seesaw. Similarly, previous work has shown that an increase in Antarctic temperature is linearly related to Greenland stadial duration, irrespective of the occurrence of Heinrich events during these stadials3. Recent studies suggest that decreases in AMOC strength precede the Heinrich events, alowing the possibility that the later are a response to AMOC variations rather than their cause27,28. The main influence of Heinrich events on the bipolar seesaw may thus be to lengthen the stadial periods during which they occur by suppressing the AMOC by means of iceberg-delivered freshwater, alowing large amounts of heat to build up in the Southern Ocean. Although both the data5and the models19,21suggest a fast temper- ature response in the South Atlantic (decadal-scale in models), there is currently no consensus in existing model studies on the physical mechanisms and timescales of the processes that propagate temper- ature signals between the South Atlantic and the Southern Ocean; complicating factors to consider include the lack of a zonal boundary to support the propagation of Kelvin waves21, the large thermal inertia of the Southern Ocean6, the complexity of eddy heat transport across the Antarctic circumpolar current, and the coupling between the extent of Antarctic sea ice and the overturning circulation29. In this context, our precise phasing observations provide a new constraint with which to test future model simulations seeking to capture the dynamics of milennial-scale climate variability. Online ContentMethods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper. Received 11 September 2014; accepted 3 March 2015. 1. NGRIP Project Members. High-resolution record of Northern Hemisphere climate extending into the last interglacial period.Nature431,147–151 (2004). 2. Blunier, T. & Brook, E. J. Timing of milennial-scale climate change in Antarctica and Greenland during the last glacial period.Science291,109–112 (2001). 3. EPICA Community Members. One-to-one coupling of glacial climate variability in Greenland and Antarctica.Nature444,195–198 (2006). 4. Crowley, T. J. North Atlantic Deep Water cools the southern hemisphere. Paleoceanography7,489–497 (1992). 5. Barker, S.et al.Interhemispheric Atlantic seesaw response during the last deglaciation.Nature457,1097–1102 (2009). 6. Stocker, T. F. & Johnsen, S. J. A minimum thermodynamic model for the bipolar seesaw.Paleoceanography18,1087 (2003). 7. Petersen, S. V., Schrag, D. P. & Clark, P. U. A new mechanism for Dansgaard– Oeschger cycles.Paleoceanography28,24–30 (2013). 8. Rind, D.et al.Efects of glacial meltwater in the GISS coupled atmosphere–ocean model. 2. A bipolar seesaw in Atlantic Deep Water production.J. Geophys. Res. 106(D21) 27355–27365 (2001). 9. Dokken, T. M., Nisancioglu, K. H., Li, C., Batisti, D. S. & Kissel, C. Dansgaard– Oeschger cycles: interactions between ocean and sea ice intrinsic to the Nordic seas.Paleoceanography28,491–502 (2013). 10. Pedro, J. B.et al.The last deglaciation: timing the bipolar seesaw.Clim. Past.7, 671–683 (2011). 11. Piotrowski, A. M., Goldstein, S. L., Hemming, S. R., Fairbanks, R. G. & Zylberberg, D. R. Oscilating glacial northern and southern deep water formation from combined neodymium and carbon isotopes.Earth Planet. Sci. Let.272, 394–405 (2008). 12. Martrat, B.et al.Four climate cycles of recurring deep and surface water destabilizations on the Iberian margin.Science317,502–507 (2007). 13. Peterson, L. C., Haug, G. H., Hughen, K. A. & Ro¨hl, U. Rapid changes in the hydrologic cycle of the tropical Atlantic during the last glacial.Science290,1947–1951 (2000). 14. Chiang, J. C. H. & Bitz, C. M. Influence of high latitude ice cover on the marine Intertropical Convergence Zone.Clim. Dyn.25,477–496 (2005). 15. Toggweiler, J. R. & Lea, D. W. Temperature diferences between the hemispheres and ice age climate variability.Paleoceanography25,PA2212 (2010). 16. Timmermann, A.et al.Towards a quantitative understanding of milennial-scale Antarctic warming events.Quat. Sci. Rev.29,74–85 (2010). 17. WAIS Divide Project Members. Onset of deglacial warming in West Antarctica driven by local orbital forcing.Nature500,440–444 (2013). 18. Velinga, M. & Wood, R. A. Global climatic impacts of a colapse of the Atlantic thermohaline circulation.Clim. Change54,251–267 (2002). 19.Liu, Z. Y. & Alexander, M. Atmospheric bridge, oceanic tunnel, and global climatic teleconnections.Rev. Geophys.45,RG2005 (2007). 20. Masuda, S.et al.Simulated rapid warming of abyssal North Pacific waters.Science 329,319–322 (2010). 21. Schmitner, A., Saenko, O. A. & Weaver, A. J. Coupling of the hemispheres in observations and simulations of glacial climate change.Quat. Sci. Rev.22, 659–671 (2003). 22. England, M. H. The age of water and ventilation timescales in a global ocean model. J. Phys. Oceanogr.25,2756–2777 (1995). 23. Buizert, C.et al.The WAIS Divide deep ice core WD2014 chronology. Part 1. Methane synchronization (68–31 ka BP) and the gas age–ice age diference.Clim. Past.11,153–173 (2015). 24. Baumgartner, M.et al.NGRIP CH4 concentration from 120 to 10 kyr before present and its relation to ad15N temperature reconstruction from the same ice core.Clim. Past10,903–920 (2014). 25. Wolf, E. W., Rankin, A. M. & Ro¨thlisberger, R. An ice core indicator of Antarctic sea ice production?Geophys. Res. Let.30,2158 (2003). 26. Marcot, S. A.et al.Centennial-scale changes in the global carbon cycle during the last deglaciation.Nature514,616–619 (2014). 27. Gutjahr, M. & Lippold, J. Early arrival of Southern Source Water in the deep North Atlantic prior to Heinrich event 2.Paleoceanography26,PA2101 (2011). 28. Marcot, S. A.et al.Ice-shelf colapse from subsurface warming as a trigger for Heinrich events.Proc. Natl Acad. Sci. USA108,13415–13419 (2011). 29. Ferrari, R.et al.Antarctic sea ice control on ocean circulation in present and glacial climates.Proc. Natl Acad. Sci. USA111,8753–8758 (2014). 30. Parrenin, F.et al.Synchronous change of atmospheric CO2and Antarctictemperature during the last deglacial warming.Science339,1060–1063 (2013). Supplementary Informationis available in the online version of the paper. AcknowledgementsWe thank the WAIS Divide Driling Team (2006–2013) (E. Morton, P. Cassidy, M. Jayred, J. Robinson, S. Polishinski, J. Koehler, L. Albershardt, J. Goetz, B. Gross, R. Kulin, S. Haman, W. Neumeister, C. Zander, J. Kyne, L. Augustin, B. Folmer, S. B. Hansen, E. Alexander and J. Fowler) and the dozens of core handlers who processed the ice core in the field and at the National Ice Core Laboratory (NICL). This work is funded through the US National Science foundation grants 0944078, 0841308 (to M.A.), 1043528 (to R.B.A., D.E.V. and J.M.F.), 1142173 (to R.B.), 1204172, 1142041, 1043518 (to E.J.B.), 0839066 (to J.C.-D.), 0087345, 0944191 (to H.C. and E.D.W.), 0539232, 0537661 (to K.M.C.), 1142069, 1142115 (to N.W.D.), 0841135 (to IDDO), 0839093, 1142166 (to J.R.M.), 0440819, 1142164 (to K.C.M.), 1142178 (to P.B.P.), 0538657 (to J.P.S.), 1043500, 0944584 (to T.A.S.), 1043313 (to M.K.S), 0537930, 1043092 (to E.J.S.), 0230149, 0230396, 0440817, 0440819, 0944191, 0944348 (to K.C.T.), 0944266 (to M.S.T.), 0839137 (to K.C.W. and K.N.), 0537593 and 1043167 (to J.W.C.W.); the USGS Climate and Land Use Change Program (to G.D.C. and J.J.F.); the NOAA Climate and Global Change Felowship Program, administered by the University Corporation for Atmospheric Research (to C.B.); the Vilum Foundation (to M.W.); the Joint Institute for the Study of the Atmosphere and Ocean (to J.B.P., JISAO contribution no. 2343); and the Korea Polar Research Institute, grant PE15010 (to J.A.). The National Science Foundation Ofice of Polar Programs also funded the WAIS Divide Science Coordination Ofice at the Desert Research Institute of Nevada and University of New Hampshire for the colection and distribution of the WAIS Divide ice core and related tasks; the Ice Driling Program Ofice and Ice Driling Design and Operations group for coring activities; the NICL for curation of the core; Raytheon Polar Services for logistics support in Antarctica; and the 109th New York Air National Guard for airlift in Antarctica. Author ContributionsData analysis andDage modeling were performed by C.B.; annual-layer counting (dating) of upper 2,850 m by M.S., T.J.F., M.W., K.C.T. and K.C.M.; CH4synchronization (dating) of lower 555 m by C.B., K.M.C., J.P.S. and T.J.F.; age scalevalidation by N.W.D., N.I., K.C.W., K.N. and T.E.W.; discrete water isotope analysis by E.J.S., A.J.Sc. and S.W.S.; continuous water isotope analysis by J.W.C.W., T.R.J., B.H.V. and V.G.; discrete CH4analysis by T.A.S., L.E.M., J.E.L., J.S.E., J.L.R. and E.J.B.; continuousCH4analysis by R.H.R., E.J.B. and J.R.M.; CO2analysis by S.A.M., M.L.K., T.K.B., J.A. andE.J.B.;d15NofN2analysis by D.B., C.B., A.J.O. and J.P.S.; continuous-flow chemicalanalysis by M.S., O.J.M., N.J.C., D.R.P. and J.R.M.; discrete chemical analysis by J.C.-D., D.G.F., B.G.K., K.K. and G.J.W.; ice core physical properties by R.B.A., J.M.F., D.E.V., M.K.S. and J.J.F.; borehole logging by R.C.B. and G.D.C.; biological studies by J.C.P. and P.B.P.; temperature reconstructions by K.M.C. and G.D.C.; tephrochronology by N.W.D. and N.I.; firn studies by M.A., T.A.S. and S.G.;10Be analysis by K.C.W. and T.E.W.; field science oversight, D.E.V. and B.H.V.; site selection by H.C., E.D.W. and E.C.P.; science management and sample distribution by M.S.T. and J.M.S.; logistics support, planning and management by M.J.K.; driling management by A.J.Sh., C.R.B., D.A.L., and A.W.W.; deep dril design by A.J.Sh., J.A.J., N.B.M. and C.J.G.; driling field management by J.A.J., K.R.S. and N.B.M.; sample colection and dril operations by C.J.G., J.J.G., T.W.K. and P.J.S. The field sample handling leaders were A.J.O., B.G.K., P.D.N. and G.J.W.; sample curation, processing and distribution was performed by G.M.H., B.A., R.M.N., E.C. and B.B.B.; the overal WAIS Divide Project design and management, Chief Scientist and field leader was K.C.T. The manuscript was writen by C.B., E.J.S. and J.B.P. with assistance from J.P.S., B.R.M., E.J.B. and K.C.T; al authors discussed the results and contributed to improving the final manuscript. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the US Government. Author InformationReprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Readers are welcome to comment on the online version of the paper. Correspondence and requests for materials should be addressed to C.B. (buizertc@science. oregonstate.edu). WAIS Divide Project Members Christo Buizert1, Betty Adrian2, Jinho Ahn3, Mary Albert4, Richard B. Aley5, Daniel Baggenstos6, Thomas K. Bauska1, Ryan C. Bay7, Brian B. Bencivengo2,CharlesR. Bentley8, Edward J. Brook1, Nathan J. Chelman9,GaryD.Clow10, Jihong Cole-Dai11, Howard Conway12, Eric Cravens13, Kurt M. Cufey14, Nelia W. Dunbar15,JonS. Edwards1,JohnM.Fegyveresi5,DaveG.Ferris11, Joan J. Fitzpatrick16,T.J.Fudge12, Chris J. Gibson8, Vasileios Gkinis17,18, Joshua J. Goetz8, Stephanie Gregory4, Geofrey M. Hargreaves2,NelsIverson15, Jay A. Johnson8, Tyler R. Jones17,Michael L. Kalk1, Matthew J. Kippenhan19, Bess G. Kofman20,KarlKreutz21,TannerW.Kuhl8, Donald A. Lebar8, James E. Lee1, Shaun A. Marcott1,22, Bradley R. Markle12, Olivia J. Maseli9, Joseph R. McConnel9, Kenneth C. McGwire9, Logan E. Mitchel1, Nicolai B. Mortensen8, Peter D. Nef23, Kunihiko Nishizumi24, Richard M. Nunn2,AnaisJ. Orsi6,25, Daniel R. Pasteris9,JoelB.Pedro18,26, Erin C. Pettit27,P.BufordPrice7,John C. Priscu28, Rachael H. Rhodes1, Julia L. Rosen1, Andrew J. Schauer12,SpruceW. Schoenemann12, Paul J. Sendelbach8, Jefrey P. Severinghaus6, Alexander J. Shturmakov8, Michael Sigl9, Kristina R. Slawny8, Joseph M. Souney29,ToddA. Sowers5, Mathew K. Spencer30, Eric J. Steig12, Kendrick C. Taylor9,MarkS.Twickler29, Bruce H. Vaughn17, Donald E. Voigt5, Edwin D. Waddington12, Kees C. Welten24, Anthony W. Wendricks8, James W. C. White17, Mai Winstrup12,18, Giford J. Wong31& Thomas E. Woodruf32 1Colege of Earth, Oceanic and Atmospheric Sciences, Oregon State University, Corvalis, Oregon 97331, USA.2US Geological Survey National Ice Core Laboratory, Denver, Colorado 80225, USA.3School of Earth and Environmental Science, Seoul National University, Seoul 151-742, Korea.4Thayer School of Engineering, Dartmouth Colege, Hanover, New Hampshire 03755, USA.5Department of Geosciences, Pennsylvania State University, University Park, Pennsylvania 16802, USA.6Scripps Institution of Oceanography, University of California at San Diego, La Jola, California 92093, USA. 7Department of Physics, University of California at Berkeley, Berkeley, California 94720, USA.8Ice Driling Design and Operations, University of Wisconsin-Madison, Madison, Wisconsin 53706, USA.9Desert Research Institute, Nevada System of Higher Education, Reno, Nevada 89512, USA.10US Geological Survey, Boulder, Colorado 80309, USA. 11Department of Chemistry and Biochemistry, South Dakota State University, Brookings, South Dakota 57007, USA.12Department of Earth and Space Sciences, University of Washington, Seatle, Washington 98195-1310, USA.13ADC Management Services, Lakewood, Colorado 80226, USA.14Department of Geography, University of California at Berkeley, Berkeley, California 94709, USA.15Earth and Environmental Science Department, New Mexico Tech, Socorro, New Mexico 87801, USA.16US Geological Survey, Denver, Colorado 80225, USA.17Institute of Arctic and Alpine Research, University of Colorado, Boulder, Colorado 80309-0450, USA.18Centre for ice and climate, University of Copenhagen, DK-2100 Copenhagen Ø, Denmark.19Antarctic Support Contract, Lockheed Martin US Antarctic Program, Centennial, Colorado 80112, USA. 20Lamont-Doherty Earth Observatory, Columbia University, Palisades, New York 10964, USA.21Climate Change Institute and School of Earth and Climate Sciences, University of Maine, Orono, Maine 04469, USA.22University of Wisconsin-Madison, Madison, Wisconsin, Wisconsin 53706 USA. 23Antarctic Research Centre, Victoria University of Welington, Welington 6012, New Zealand.24Space Sciences Laboratory, University of California at Berkeley, Berkeley, California 94720, USA.25Laboratoire des Sciences du Climat et de l’Environnement, Institut Pierre Simon Laplace, 91191 Gif-Sur-Yvete, France.26Joint Institute for the Study of the Atmosphere and Ocean, University of Washington, Seatle, Washington 98195, USA.27Department of Geosciences, University of Alaska Fairbanks, Fairbanks, Alaska 99775, USA. 28Department of Land Resources and Environmental Sciences, Montana State University, Bozeman, Montana 59717, USA.29Institute for the Study of Earth, Oceans and Space, University of New Hampshire, Durham, New Hampshire 03824, USA. 30School of Physical Sciences, Lake Superior State University, Sault Sainte Marie, Michigan 49783, USA.31Department of Earth Sciences, Dartmouth Colege, Hanover, New Hampshire 03755, USA.32PRIME Laboratory, Purdue University, West Lafayete, Indiana 47907, USA. METHODS Core recovery and processing.The location of the WAIS Divide ice core (79.48uS, 112.11uW) was selected to provide the highest possible time-resolution record of Antarctic climate during the past 50 kyr or more, and to ensure that the record from the last deglaciation would not be in ice from the lower-quality britle- ice zone17. The site has a present-day ice accumulation rate of 22 cm yr21, divide flow and an ice thickness of 3,460 m (J. Paden, personal communication 2011). No stratigraphic irregularities were detected in the core using visual observation, multitrack electrical conductivity measurements or gas stratigraphy. Coring was stopped 50 m above the bed to prevent contamination of the unfrozen basal environment. Driling of the main core (WDC06A) started during the 2006/ 2007 field season and finished during the 2011/2012 field season, at the final depth of 3,405 m. Driling was done with the US Deep Ice Sheet Coring (DISC) dril, a new dril designed and operated by the Ice Driling Design and Operations group (engineering team at University of Wisconsin, Madison). Core handling and distribution of ice samples were performed by the National Ice Core Laboratory (US Geological Survey) and the Science Coordination Ofice (University of New Hampshire and Desert Research Institute). The dril and core handling methods had many innovations that improved core quality, especialy in the britle ice. These included improved mechanical characteristics of the dril, beter control of the driling operation, core handling procedures that minimized thermal and mechanical shock to core, neting to contain the britle ice, and temperature tracking of al ice samples from the dril to the laboratories. Data description.Water18O/16O composition (d18O) was measured at IsoLab, University of Washington, Seatle, Washington. Measurement procedures for WDC have been described elsewhere17,31. In short, measurements were made at 0.5 m depth averaged resolution, using laser spectroscopy (Picarro L2120-iwater isotope analyser).d18O data are reported relative to the VSMOW (Vienna Standard Mean Ocean water) standard, and normalized to the VSMOW-SLAP (Standard Light Antarctic Precipitation) scale. Two separate data sets of atmospheric methane exist for WDC. The first is based on discrete sampling of the core; it was measured jointly at Penn State University (0–68 kyr, 0.5–2 m resolution) and at Oregon State University (11.4–24.8 kyr, 1–2 m resolution). For this record the air was extracted using a melt–refreeze technique from,50 g ice samples, and analysed on a gas chromatograph with a flame ionization detector32. Corrections for gas solubility, blank and gravita- tional enrichment for the Oregon State University data were performed in accordance with ref. 33; corrections to the Penn State University data are described in ref. 17. A second data set (R. H. Rhodes, E. J. Brook, J. C. H. Chiang, T. Blunier, O. J. Maseli, J. R. McConnel, D. Romanini and J. P. Severinghaus, unpublished observations) is based on continuous flow analysis (CFA) in combination with laser spectroscopy34,35. This data set, which is of higher temporal resolution, was used to determine the timing of the abrupt DO transitions (Extended Data Fig. 2). Al CH4data are reported on the NOAA04 scale36. WDC Na concentrations were measured at the Ultra Trace Chemistry Laboratory at the Desert Research Institute, Reno, Nevada, by means of CFA. Sticks of ice were melted continuously on a heated plate. Potentialy contaminated meltwater from the outer part of the ice stick was discarded. Two inductively coupled plasma mass spectrometers were used to analyse the stream of melt- water37. The CFA chemistry and CH4data were measured simultaneously on the same ice samples, and al measurements are co-registered in depth. Age scale andDage.For WDC we use the WD2014 chronology, which is based on annual-layer counting down to 2,850 m (31.2 kyrBP), and on stratigraphic match- ing by means of globaly wel-mixed methane for the deepest part of the core (2,850–3,405 m)23. The deep section of the WD2014 chronology has been syn- chronized to a linearly scaled version of the layer-counted Greenland ice core chronology (GICC05). Multiplying the GICC05 chronology by 1.0063 brings the ages of abrupt DO events as observed in Greenland ice core records in agree- ment with the ages of these events observed in the U/Th-dated Hulu Cave spe- leothem record (which provides beter absolute age control)23. The NGRIP data in Fig. 1a have similarly been placed on the 1.00633GICC05 age scale; NGRIP data in Fig. 3a are ploted on the original GICC05 chronology. WDC permits precise interhemispheric synchronization for two reasons. First, as a result of the high accumulation rates,Dage remains relatively smal. WDC Dage is about an order of magnitude smaler thanDage in cores from the East Antarctic Plateau, and about one-third ofDage at other coastal sites that cover the last glacial period (Extended Data Fig. 1a)23. As the uncertainty inDage is pro- portional toDage itself, this leads to smaler dating uncertainties at WDC. Second, recent technological developments in coupling laser spectrometers with ice core continuous-flow analysis34have resulted in a WDC CH4record with the highest temporal resolution of al Antarctic ice cores so far. Extended Data Fig. 2c com- pares the WDC CH4record (grey,,2-year sampling resolution) with the EDML CH4record (orange,,100-year sampling resolution). The estimated uncertainty in WDCDage (ref. 23) is shown in Extended Data Fig. 1b–e. Note that the transition between the layer-counted and CH4-synchronized sections of the WD2014 chronology at 31.2 kyr does not influence our analysis of the bipolar seesaw phasing. In the stacking procedure we align each of the individual Antarctic events at the midpoint of the associated DO CH4transition. In doing so, we efectively synchronize al the WDC events to NGRIPd18O, not just the events between 68 and 31.2 kyr. In performing the cross-correlation between NGRIP and WDC (Fig. 2a), we only analyse the CH4-synchronized part of the record. A dynamical version of the firn densification model23,38is used to calculate the Dage of the WD2014 chronology23. We furthermore calculated WDCDage using an alternative firn densification model23,39,40(Extended Data Fig. 1b, blue curve). Using this alternativeDage scenario gives a NH lead time of 242 and 247 years for abrupt DO warming and DO cooling, respectively, in good agreement (within the 1suncertainty bound) with results presented in the text. Identification of abrupt DO transitions.We use the midpoint of the abrupt transitions as time markers to align individual events. The method for finding the midpoint is identical to that used in ref. 23, where further details can be found. The method is shown in Extended Data Fig. 2a, b for DO 17 and 16 (we use a recently proposed DO nomenclature41). For each transition we define a pre-transition level and a post-transition level (horizontal orange markers), and use linear interpola- tion to find the depth at which the tracer of interest (d18OorCH4) has completed 25%, 50% and 75% of the transition from the pre-transitional to post-transitional value. The 50% depth (red dots) is used to define the timing of the abrupt event; the 25% and 75% depths (blue dots) provide an uncertainty range on the midpoint23. It has been suggested that interhemispheric CH4synchronization is compli- cated by the fact that gas bubbles in ice cores have a gas age distribution, rather than a single age42. To investigate this efect at WDC, we construct a hypothetical gas age distribution (Extended Data Fig. 2d), using a (truncated) log-normal distribution as suggested elsewhere42,43, in which we conservatively sets51 andm5ln(50). The spectral width (or second moment) of this age distribution44 is 24.3 years; for comparison, the gas-age spectral width for present-day WDC firn is only 4.8 years45. Applying this filter to an atmospheric ramp in CH4(Extended Data Fig. 2e) shifts the midpoint of the transition (as recorded in the ice core) forwards by 5 years. Because we have conservatively chosen the age distribution to be very wide, this 5-year shift should be considered an upper bound. The influence of the gas age distribution on the transition midpoint identification is very smal in comparison with other sources of uncertainty, and is neglected in the remainder of the manuscript. Stacking procedure.We define a time vectortthat runs from timet521,200 to t51,200 in 1-year increments. Next, for each of the DO/AIM events, we set the midpoint of the NGRIPd18O transition tot50, and use linear interpolation to sample the diferent records onto time vectort. For each individual DO/AIM event this results in a record that spans from21,200 to11,200 years relative to the abrupt transition, at annual time steps; we shal refer to these as the contrib- utory records. To align the WDC records we set the midpoint in the WDC CH4 transition tot556 years, using the observation that the CH4transition slightly lags the Greenlandd18O signal23,24,46,47(the uncertainty in this time lag is evaluated in the sensitivity study). With al the individual events synchronized and resampled to identical time spacing, we can average them to obtain the stacked record. The spacing between DO events varies widely over time, ranging from several thousands of years to only hundreds of years. This means that in several cases the 2,400-year window we use contains more than a single event. We crop the con- tributory records whenever an adjacent event occurs within the sampling window (Extended Data Fig. 3). For widely spaced (that is, long-duration) events such as DO/AIM 12, no adjacent events fal within the time window, and no cropping is required (Extended Data Fig. 3a). For closely spaced events, such as DO/AIM 17.1, we need to crop neighbouring events (Extended Data Fig. 3b). The cropped parts of the contributory records are replaced with constant values that equal the bound- ary values (50-year averages) of the uncropped part of the record. The number of contributory records available as a function of time is shown in the botom half of Extended Data Fig. 4. The shortest contributory record is the AIM 17.1 record, which is 700 years long. Determining the breakpoint in the WDCd18O stack.We use an automated routine that is similar to the BREAKFIT algorithm48, the main diference being the use of a second-order polynomial (rather than a linear) fit to the data, to account for the fact that the WDCd18O stack is curved rather than linear. The algorithm finds the breakpoint that minimizes the root mean square deviation (r.m.s.d.) between the stack (that is, the data) and the polynomial fit in the2600 to1700 years time interval (Extended Data Fig. 4); the interval was chosen asymmetricaly around zero to account for the fact that the breakpoint is always found at positive values. The routine operates as folows. First, the user selects the range in which the algorithm is to look for the break- point, as wel as a time step. For example, if the user selects 0–400 years at a 100- year time step, the algorithm wil test the possibilitiest50, 100, 200, 300 and 400 years. We shal refer to this one-dimensional array of values as the input vector. Second, at each of the values in the input vector, thed18O stack is split into two pieces, and a second-order polynomial curve is fited to data on each side sepa- rately. The two fiting curves are merged at their point of intersection. Note that the point of intersection can difer from the point where the data series were split in two. Third, for each fit thus obtained, we calculate the r.m.s.d. to the stacked data over the2600 to1700 years time interval, and we select the best fit. The point of intersection of the two curves, rather than the point where the data series were split into two (that is, the values in the input vector), is considered the breakpoint in the d18O stack. For most of the results in this work we let the algorithm search in the 0–400-year range, at a 2-year time step. The exception is the stacks of randomly selected DO events (Extended Data Fig. 6), for which we used a wider range (2100 to 1500 years) to account for the larger spread in solutions. Note that the efective range of breakpoint detection is wider than the input range, because the point of intersection of the fiting curves is permited to lie outside the range specified by the input vector. The fiting curves generated by the algorithm are ploted in dark orange lines in Fig. 2b, c and Extended Data Figs 4a, b, 6a–d and 7b, c; the break- point is indicated with an orange dot. The MATLAB code implementing the fiting procedure is provided in Supplementary Information. The reason for using a quadratic rather than a linear fiting procedure is that Antarctic temperatures do not change in a linear fashion either before or after the abrupt transitions. As a test, we apply our fiting procedure to the WDCd18O stack (Fig. 2b), using either linear or quadratic fiting functions, in which we change the size of the data window from 300 years to 1,600 years. The window is applied symmetricaly, meaning that equal numbers of years are used before and after the detected breakpoint (800 years on each side for a 1,600-year window). Data faling outside this window are ignored in the fiting procedure. We find that the performance of the quadratic fiting procedure is superior in two important aspects (Extended Data Fig. 5). First, by using a quadratic fit, the breakpoint we detect is independent of the window size, which is not true for the linear fit. Second, the quadratic method provides a beter approximation to the data, as expressed by the r.m.s.d. between the data and the fiting curve. At smal window sizes of,700 years, both methods agree very wel, and the fit to the data is comparable in terms of the r.m.s.d. This is due to the fact that at smal window sizes the curvature of the time series is not very pronounced and is wel approximated by a linear trend. We have also tested the BREAKFIT routine48, which uses a linear fit on both sides of the breakpoint. We find a very good agreement between the breakpoints identified using the BREAKFIT routine and our linear fiting routine, with a mean diference of 0.2 years between them. Using a moving block bootstrap analysis48, the BREAKFIT algorithm suggests a 2suncertainty of 50.8 years in breakpoint identification. Monte Carlo sensitivity study.To assess how the diferent uncertainties in our method afect our conclusions, we perform a Monte Carlo sensitivity study. The largest source of uncertainty is the diference between gas age and ice age, orDage (ref. 49).Dage determines the relative timing of the WDC CH4andd18O records, and thereby the phasing we infer for WDC AIM events relative to Greenland climate2,3,50. IncreasingDage shifts WDC ice ages towards older ages, thereby decreasing the NH lead time that we infer; vice versa, decreasingDage wil shift the ice ages younger, increasing the observed NH lead time. TheDage reconstruc- tion we use here is obtained from a dynamical firn-densification model40,47,51–54 constrained byd15N-N2data, a proxy for past firn-column thickness55. The firn- densification model results are in good agreement23with the alternativeDdepth method developed recently56. For our current purpose we use an ensemble of 1,000 Dage histories (Extended Data Fig. 1b–e) generated by varying a suite of input parameters of the firn-densification model23. We furthermore use four diferent depth–age interpolation strategies23,57, providing a total of 1,0003454,000 alternative WDC chronologies. For each of these 4,000 potential WDC chronologies we further perform a Monte Carlo simulation (100 realizations) in which we perturb the alignment of the contributory records. The rationale behind the Monte Carlo study is that there is an uncertainty in the synchronization of the contributory records, and we want to assess how this uncertainty impacts our conclusions. In the Monte Carlo study, a random timing error is assigned to each of the contributory records, and thereby they are efectively shifted relative to one another in time. We distinguish between systematic and non-systematic errors, in which the former are correlated between the diferent AIM events, and the later are not. Sensitivity tests with the firn densification model show thatDage is a systematic error; if we have overestimated glacialDage by, for example, 50 years, this wil probably be true for al of the AIM events. We consider five independent sources of uncertainty in the synchroniza- tion procedure. First, the midpoint detection in the abrupt NGRIPd18O transitions has an uncertainty; we shift the record by a value drawn from a Gaussian distribution of zero mean and a width that equals the midpoint uncertainty as defined above. This error is non-systematic. Second, as in the first point above, but for the midpoint in the abrupt WDC CH4 transition. Third, although al four interpolation strategies synchronize the abrupt NH warming events between NGRIP and WDC, only two of them synchronize the NH cooling events23. For stacks of NH cooling we added an additional uncertainty drawn from a Gaussian distribution of zero mean and a width corresponding to the age diference between the cooling event as detected in NGRIPd18O and WDC CH4. This error is non-systematic. Fourth, a recent extensive analysis suggests a lag of CH4behind Greenlandd18O of 56638 years (2s)24. There are two types of uncertainty to consider: first, the CH4lag can difer between individual DO events (non-systematic error), and second, the central estimate of 56 years has an inherent uncertainty (systematic error). We here assume that both are equaly large, and in the stacking procedure we assign both a 27-year systematic error and a 27-year non-systematic error to the CH4lag time. Both are drawn independently from a Gaussian distribution, and their combined uncertainty is therefore (2721272)1/2538 years (that is, the 2s uncertainty bound24). Fifth, the breakpoint fiting procedure has an intrinsic uncertainty. We estimate this by using the moving-block bootstrap analysis of the BREAKFIT algorithm48, which suggests a 1suncertainty of 25.4 years (Extended Data Fig. 5a). We there- fore apply a 50.8-year (2s) systematic uncertainty (systematic because it applies to the fuld18O stack rather than to individual events). Each contributory record is then individualy shifted by an amount that corre- sponds to the sum of the randomly assigned non-systematic errors. After realign- ing the contributory records in this manner, the entire stack is shifted in time by an amount that corresponds to the sum of randomly assigned systematic errors. We detect the breakpoint by using the fiting algorithm. For each of the 4,000 chro- nologies we test 100 realizations of the Monte Carlo realignment, thus obtaining a total of 43105(1,000343100) estimates of the WDCd18O breakpoint. No statistical methods were used to predetermine the sample size used in the Monte Carlo analysis. The outcome of the sensitivity study is presented in the histograms of Fig. 2e, f and Extended Data Fig. 6e–h. The non-systematic errors, because they are independent of each other, accu- mulate slowly when added in quadrature. The systematic errors (such as theDage uncertainty or the breakpoint fiting procedure) therefore dominate the final uncertainty estimate. Running the Monte Carlo experiment withDage uncertain- ties withheld gives a 2suncertainty bound on the interpolar phasing of 60 and 68 years for the DO onset and DO termination stacks, respectively; using only Dage uncertainties gives a 2suncertainty bound of 69 years. TheDage uncertainty is therefore as important as al other sources of uncertainty combined. Code availability.The MATLAB code used to stack the records and perform the Monte Carlo analysis are provided in Supplementary Information. Alternative stacks.To exclude the possibility that timing of the WDCd18O stack is dominated by a few prominent AIM events we perform the folowing tests. First, we distinguish between major AIM events that are very prominent in the record, and minor events. The major AIM events are AIM 17, 14, 12, 8 and 4; al of these are preceded by NH Heinrich events. These large events have historicaly been caled A-events, and are clearly visible even in low-resolution Antarctic cores (AIM 8, 12, 14 and 17 have been referred to in earlier work as A1, 2, 3 and 4, respectively)2,3. The minor events are the remaining ones (AIM 3, 5.1, 5.2, 6, 7, 9, 10, 11, 13, 15, 16 and 18). When we compare a stack consisting of only the major AIM events with a stack of only the minor AIM events, we find no statisticaly significant diference in their timing (Extended Data Fig. 6a, b). The NH lead time for the major and minor AIM stacks is 230693 and 211694 years, respectively (Extended Data Fig. 6e, f); the timing of the ful stack (218692 years) is inter- mediate to these values, as expected (Fig. 2b). In a second experiment we randomly selected eight DO/AIM events in the stacking procedure for each of the 43105realizations of the Monte Carlo experi- ment; one such realization is shown in Extended Data Fig. 6c, d for abrupt NH warming and cooling, respectively. The NH lead time we find in this experi- ment is 2106113 years and 2086108 years for abrupt NH warming and cooling, respectively (Extended Data Fig. 6g, h). These mean values difer by only a few years from those found for the ful WDCd18O stack (Fig. 2). Note that the distribution width of NH lead times in this experiment is wider than the distribution width found for the fuld18O stacks. The reason for this may be threefold: first, a stack with fewer contributory records has a smaler signal-to- noise ratio, making the breakpoint detection less certain; second, if the NH lead time difers slightly between the AIM events, the average timing of a randomly selected subset may difer from the average timing of the ful set; and third, in the Monte Carlo procedure the contribution of non-systematic uncertainties scales roughly in proportion to!N,again suggesting a narrower distribution for higherN (the number of events in the stack). The wider distribution found on decreasing the number of events in the stack therefore does not represent a more realistic uncertainty bound on the phasing of the bipolar seesaw. Stacking the WDC Na record.Last, we investigate the phasing of the WDC sea- salt Na (ssNa) record relative to Greenland and Antarctic climate. The sea-ice surface is a major source of sea-salt aerosols to the atmosphere around Antarctica, and consequently ssNa has been interpreted as a proxy for the production and extent of regional sea ice25,58,59. Sea-salt Na is strongly anti-correlated with Antarctic climate on both orbital and milennial timescales58,60. Measured ssNa concentrations can be converted to a ssNa flux by multiplying with the accumulation rates at the site. On the East Antarctic plateau, the Na flux is expected to provide a beter metric for atmospheric concentrations, because dry deposition dominates as a result of the very low accumulation rates58. Accumulation rates are much higher at WAIS Divide, however, and Greenland ice cores provide a beter analogue than those on the East Antarctic plateau. It has been shown at Greenland Summit61that sea-salt concentrations in the ice reflect the atmospheric variations within 5%, whereas the sea-salt flux underestimates the atmospheric variations by about 50%. Thus, the observed variations in WDC ssNa concentrations should to first order reflect variations in the atmospheric loading. We investigate the relative timing of the ssNa record with the same tools as those that we applied to the WDCd18O record. An evaluation of the lagged correlation between NGRIPd18O and the time derivative of the various WDC records is shown in Extended Data Fig. 7a. For ssNa we obtain a centennial-scale NH lead time, comparable to that observed for WDCd18O. Stacking of the WDC impurity records yields the same picture (Extended Data Fig. 7b, c); the breakpoint in the ssNa stacks coincides with the breakpoint in the WDCd18O stack, and lags the abrupt NH transition by roughly two centuries. Note that becaused18O and ssNa are both measured in the ice phase, the uncertainty that we estimate for the timing of thed18O breakpoint (2sof about 90 years) also holds for the impurity records. The synchronicity of WDC ssNa andd18O variations suggests that Antarctic climate and sea-ice extent are closely linked on centennial, or perhaps even sub- centennial timescales. This may reflect a common forcing and/or a feedback between Southern Ocean surface temperatures and sea-ice extent. The sea-ice extent may furthermore be important in forcing Antarctic temperatures (and d18O)6,62, particularly at the WAIS Divide site, where the marine influence is stronger than on the East Antarctic plateau17. 31. Steig, E. J.et al.Recent climate and ice-sheet changes in West Antarctica compared with the past 2,000 years.Nature Geosci.6,372–375 (2013). 32. Sowers, T.et al.An interlaboratory comparison of techniques for extracting and analyzing trapped gases in ice cores.J. Geophys. Res.102,26527–26538 (1997). 33. Mitchel, L. E., Brook, E. J., Sowers, T., McConnel, J. R. & Taylor, K. Multidecadal variability of atmospheric methane, 1000–1800 C.E.J. Geophys. Res.116,G02007 (2011). 34. Stowasser, C.et al.Continuous measurements of methane mixing ratios from ice cores.Atmos. Meas. Tech.5,999–1013 (2012). 35. Rhodes, R. H.et al.Continuous methane measurements from a late Holocene Greenland ice core: atmospheric and in-situ signals.Earth Planet. Sci. Let.368, 9–19 (2013). 36. Dlugokencky, E.et al.Conversion of NOAA atmospheric dry air CH4mole fractionsto a gravimetricaly prepared standard scale.J. Geophys. Res.110,D18306 (2005). 37. Sigl, M.et al.A new bipolar ice core record of volcanism from WAIS Divide and NEEM and implications for climate forcing of the last 2000 years.J. Geophys. Res. Atmos.118,1151–1169 (2013). 38. Herron, M. M. & Langway, C. C. Firn densification: an empirical model.J. Glaciol.25, 373–385 (1980). 39. Arnaud, L., Barnola, J. M. & Duval, P. inPhysics of Ice Core Records(ed. Hondoh, T.) 285–305 (Hokkaido Univ. Press, 2000). 40. Goujon, C., Barnola, J. M. & Ritz, C. Modeling the densification of polar firn including heat difusion: application to close-of characteristics and gas isotopic fractionation for Antarctica and Greenland sites.J. Geophys. Res.108(D24), 4792 (2003). 41. Rasmussen, S. O.et al.A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: refining and extending the INTIMATE event stratigraphy.Quat. Sci. Rev. 106,14–28 (2014). 42. Ko¨hler, P. Rapid changes in ice core gas records. Part 1. On the accuracy of methane synchronisation of ice cores.Clim. Past Discuss.6,1453–1471 (2010). 43. Ko¨hler, P., Fischer, H. & Schmit, J. Atmosphericd13CO2and its relation to pCO2and deep oceand13C during the late Pleistocene.Paleoceanography25,PA1213 (2010). 44. Trudinger, C. M.et al.Reconstructing atmospheric histories from measurements of air composition in firn.J. Geophys. Res.107 (D24) 4780 (2002). 45. Batle, M. O.et al.Controls on the movement and composition of firn air at the West Antarctic Ice Sheet Divide.Atmos. Chem. Phys.11,11007–11021 (2011). 46. Huber, C.et al.Isotope calibrated Greenland temperature record over Marine Isotope Stage 3 and its relation to CH4.Earth Planet. Sci. Let.243,504–519 (2006).47. Rasmussen, S. O.et al.A first chronology for the North Greenland Eemian Ice Driling (NEEM) ice core.Clim. Past9,2713–2730 (2013). 48. Mudelsee, M. Break function regression.Eur. Phys. J. Spec. Top.174,49–63 (2009). 49. Schwander, J. & Staufer, B. Age diference between polar ice and the air trapped in its bubbles.Nature311,45–47 (1984). 50. Capron, E.et al.Synchronising EDML and NorthGRIP ice cores usingd18Oof atmospheric oxygen (d18Oatm) and CH4measurements over MIS5 (80–123 kyr).Quat. Sci. Rev.29,222–234 (2010). 51. Barnola, J. M., Pimienta, P., Raynaud, D. & Korotkevich, Y. S. CO2–climaterelationship as deduced from the Vostok ice core: a re-examination based on new measurements and on a re-evaluation of the air dating.Telus43,83–90 (1991). 52. Schwander, J.et al.Age scale of the air in the summit ice: implication for glacial– interglacial temperature change.J. Geophys. Res.102,19483–19493 (1997). 53. Seierstad, I.et al.Consistently dated records from the Greenland GRIP, GISP2 and NGRIP ice cores for the past 104 ka reveal regional milennial-scale isotope gradients with possible Heinrich Event imprint.Quat. Sci. Rev.106,29–46 (2014). 54. Buizert, C.et al.Greenland temperature response to climate forcing during the last deglaciation.Science345,1177–1180 (2014). 55. Sowers, T., Bender, M., Raynaud, D. & Korotkevich, Y. S.d15NofN2in air trapped inpolar ice: a tracer of gas transport in the firn and a possible constraint on ice age– gas age diferences.J. Geophys. Res.97,15683–15697 (1992). 56. Parrenin, F.et al.On the gas–ice depth diference (Ddepth) along the EPICA Dome C ice core.Clim. Past8,1239–1255 (2012). 57. Fudge, T. J., Waddington, E. D., Conway, H., Lundin, J. M. D. & Taylor, K. Interpolation methods for Antarctic ice-core timescales: application to Byrd, Siple Dome and Law Dome ice cores.Clim. Past Discuss.10,65–104 (2014). 58. Fischer, H., Siggaard-Andersen, M.-L., Ruth, U., Ro¨thlisberger, R. & Wolf, E. Glacial/ interglacial changes in mineral dust and sea-salt records in polar ice cores: sources, transport, and deposition.Rev. Geophys.45,RG1002 (2007). 59. de Vernal, A., Gersonde, R., Goosse, H., Seidenkrantz, M.-S. & Wolf, E. W. Sea ice in the paleoclimate system: the chalenge of reconstructing sea ice from proxies—an introduction.Quat. Sci. Rev.79,1–8 (2013). 60. Wolf, E. W.et al.Southern Ocean sea-ice extent, productivity and iron flux over the past eight glacial cycles.Nature440,491–496 (2006). 61. Aley, R.et al.Changes in continental and sea-salt atmospheric loadings in central Greenland during the most recent deglaciation: model-based estimates.J. Glaciol. 41,503–514 (1995). 62. Noone, D. & Simmonds, I. Sea ice control of water isotope transport to Antarctica and implications for ice core interpretation.J. Geophys. Res.109,D07105 (2004). 63. Veres, D.et al.The Antarctic ice core chronology (AICC2012): an optimized multi- parameter and multi-site dating approach for the last 120 thousand years.Clim. Past9,1733–1748 (2013). 64. Bazin, L.et al.An optimized multi-proxy, multi-site Antarctic ice and gas orbital chronology (AICC2012): 120–800 ka.Clim. Past9,1715–1731 (2013). 65. Kawamura, K.et al.Northern Hemisphere forcing of climatic cycles in Antarctica over the past 360,000 years.Nature448,912–916 (2007). 66. Schilt, A.et al.Atmospheric nitrous oxide during the last 140,000 years.Earth Planet. Sci. Let.300,33–43 (2010). 67. Ko¨hler, P., Knorr, G., Buiron, D., Lourantou, A. & Chappelaz, J. Abrupt rise in atmospheric CO2at the onset of the Bøling/Alerød: in-situ ice core data versustrue atmospheric signals.Clim. Past7,473–486 (2011). 68. Rosen, J. L.et al.An ice core record of near-synchronous global climate changes at the Bøling transition.Nature Geosci.7,459–463 (2014). 69. Ko¨hler, P., Knorr, G. & Bard, E. Permafrost thawing as a possible source of abrupt carbon release at the onset of the Bøling/Alerød.Nature Comm5,5520 (2014). 010203040506070 200 300 400 500 600 700800 1000 2000 3000 4000 5000 6000 7000 Vostok EDC Dome Fuji EDML TALDICE WDC Δa ge (y ear s) 0102030405060700 200 400 600 800 WD2014 Age (ka BP) Δa ge (y ear s) Al ± 2σ ± 1σ 200 250 300 3500 50 100 150 200 60 ka BP: 262 ± 50 year Fre qu en cy Δage (year) 250 300 350 400 450 50 100 150 200 40 ka BP: 351 ± 73 year Fre qu en cy Δage (year) 400 500 600 700 50 100 150 200 20 ka BP: 521 ± 120 year Fre qu en cy Δage (year) a b c d e Extended Data Figure 1|Diference between gas age and ice age (Dage) at WAIS Divide. a, Comparison of WDCDage with other Antarctic cores. Ice core abbreviations: EDC, EPICA Dome Concordia; EDML, EPICA Dronning Maud Land; TALDICE, Talos Dome; WDC, WAIS Divide.Dage values are taken from refs 23, 63–65. The vertical axis is on a logarithmic scale.b,Dage uncertainty bounds obtained from an ensemble of 1,000 alternativeDage scenarios; details are given elsewhere23.ADage scenario obtained with an alternative densification model (ref. 39 instead of ref. 38) is shown in blue. c–e, Histograms of the 1,000Dage scenarios at 20 kyrBP(c), 40 kyrBP (d) and 60 kyrBP(e); stated values give the distribution mean6the 2s standard deviation. 5858.55959.560 −44 −42 −40 −38 NG RIP δ 18 O ( ‰) 450 500 550 600 WD C C H 4 (p pb) 5858.55959.560 450 500 550 600 ED ML CH 4 ( pp b) Age (ka BP) −40−20 0 20 40 60 80 100120450 500 550 600 650 700 750 Time (years) CH 4 ( pp b) midpoint shifted by −5 years Atmosphere Ice core (filtered) −50 0 50 100 1500 0.01 0.02 Gas Age De nsit y ( ye ars −1 ) Gas age distributiona b c d e Extended Data Figure 2|Determining the timing of the abrupt DO transitions. a,b, DO 17.2, 17.1, 16.2 and 16.1 (from oldest to youngest41)as recorded in NGRIPd18O(a) and WDC CH4(b). Horizontal orange bars denotepre-transition and post-transition levels; the transition midpoint (50% of signal amplitude) is indicated by a red dot; the 25% and 75% signal amplitude markers are indicated with blue dots.c, Comparison of WDC CH4(grey) with EDMLCH4(orange)3,50,66.d, Hypothetical gas-age distribution for WDC due to firndensification and gradual bubble closure, using a truncated log-normal distribution67.e, Shift in transition midpoint induced by filtering of the atmospheric record in the firn column. −2000 −1000 0 1000 2000 −44 −42 −40 −38 −36 Time (year) NG RIP δ 18 O ( ‰) −41 −40 −39 −38 DO/AIM 12 WD C δ18 O ( ‰) −2000 −1000 0 1000 2000 −44 −42 −40 −38 Time (year) NG RIP δ 18 O ( ‰) −39 −38 −37 DO/AIM 17.1 WD C δ18 O ( ‰) a b Extended Data Figure 3|Cropping of individual records in the stack to prevent overlap of events. a, DO/AIM 12, where no cropping is needed. b, DO/AIM 17.1, where the most cropping is needed. Ful time series with five- point running average are ploted in grey, and the contributory records are ploted in blue and orange for NGRIP and WDC, respectively. The yelow vertical shading bar in background shows the NH lead time (200 years); the purple rectangle gives the21,200 to11,200 time window. −1000 −500 0 500 1000 −40 −39.5 −39 δ18 O ( ‰) −1000 −500 0 500 1000 0 5 10 15 20 Time (years) Nr. AI M r ec ord s −1000 −500 0 500 1000 −40 −39.5 −39 δ18 O ( ‰) −1000 −500 0 500 1000 0 5 10 15 20 Time (years) Nr. AI M r ec ord s a b Extended Data Figure 4|Number of records and fiting procedure. Number of contributory records to the WDCd18O stacks for abrupt NH warming (interstadial onset) (a) and for abrupt NH cooling (interstadial termination) (b). Blue rectangles indicate the time window over which the fiting procedure evaluates the fit to the data (2600 to1700 years); shaded vertical yelow bars show NH lead time. 400 600 800 1000 1200 1400 1600 80 100 120 140 160 180 200 220 NH le ad ti me (y ear s) Linear breakfit Quadratic breakfit Mudelsee (2009) breakfit 400 600 800 1000 1200 1400 1600 0.06 0.07 0.08 RM SD (‰ ) Window size (years) a b Extended Data Figure 5|Evaluating the performance of the breakpoint detection algorithm. a, Breakpoint detection as a function of data window size using both linear (blue) and quadratic (orange) functions, compared with the BREAKFIT algorithm48(grey dots with 1serror bars). The data window is applied symmetricaly, meaning that equal numbers of years (half the window size) are used before and after the detected breakpoint. Data faling outside this window are ignored in the fiting procedure.b, Root mean square deviation between the WDCd18O stack and the fiting curve. −500−250 0 250500750 0 0.5 1 1.5 2 Major AIM events Nor mali ze d si gn al Time (years) −250 0 250500750 Minor AIM events Time (years) −250 0 250500750 Random AIM events Time (years) −250 0 250500750 Random AIM events Time (years) 0 20k 40k 60k Fre qu en cy 0 100 200 300 230 ± 93 NH lead (years) 0 100 200 300 211 ± 94 NH lead (years) 0 100 200 300 210 ± 113 NH lead (years) 0 100 200 300 400 208 ± 108 NH lead (years) a b c d e f g h Extended Data Figure 6|Alternative stacking of AIM events. a, Stack of NGRIPd18O (blue), WDC CH4(green) and WDCd18O (orange with fit) forjust the major AIM events (4, 8, 12, 14 and 17), aligned at the abrupt NH warming.b,Asina, but for only the minor AIM events (3, 5.1, 5.2, 6, 7, 9, 10, 11, 13, 15, 16 and 18).c,Asina, but for eight randomly selected DO/AIM events. d,Asinc, but aligned at the abrupt NH cooling. Events are averaged with their original amplitudes and normalized after stacking for convenience of visualization.e–h, Histograms of NH lead time associated with a–d, respectively, generated by binning the 43105solutions from the sensitivity study. The distribution mean and 2suncertainty bounds are listed in the panels. Shaded vertical yelow bars (upper panels) show NH lead time. −750−500−250 0 250 500 750 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 SH leads NH leads ssNa δ18O Lag (years) Cor rel ati on −500−250 0 250 500 750 −1 0 1 2 Nor mali ze d si gn al Time (years) −250 0 250 500 750 Time (years) a b c Extended Data Figure 7|Timing of sea-salt sodium. a, Lagged correlation between NGRIPd18O and the time derivative of WDCd18O (orange), and between NGRIPd18O and the time derivative of WDC ssNa (grey). The dots indicate the maximum (anti-)correlation at 167-year and 229-year NH lead for WDCd18O and ssNa, respectively. A fourth-order Buterworth bandpass filter with a 500–10,000-year window is applied to the time series to isolate the milennial-scale variability.b, DO3–18 stack of NGRIPd18O (blue), WDC CH4 (green), WDCd18O (orange) and WDC ssNa (grey), aligned at the midpoint of the DO warming signal. The estimated breakpoint in the stacks (dots) occurs at t5218 and 195 years for WDCd18O and ssNa, respectively.c,Asinb, but for the abrupt NH cooling events, with the estimated breakpoint att5208 and 199 years for WDCd18O and ssNa, respectively. Shaded vertical yelow bars show NH lead times. Extended Data Table 1|Phasing of the bipolar seesaw during the last deglaciation For the deglaciation we use a Greenland-CH4phasing that is the average of the 56638-year CH4lag observed for the glacial period24,46, and the near-synchronous phasing (4.5624-year CH4lag) observed for the Bøling transition68. Al stated errors represent 2suncertainty bounds; uncertainty in the WDCd18O breakpoint is determined using the BREAKFIT algorithm48; uncertainty in the CH4transition is defined as the 25–75% range of the CH4transition (Extended Data Fig. 2);Dage uncertainty is found by using a firn densification model sensitivity study23(Extended Data Fig. 1); uncertainty in the Greenland-CH4phasing is taken to be the root mean square of the uncertainties in the cited studies24,68; the uncertainty in the NH climate is the root sum square of the uncertainties in the preceding columns. The potential role of CO2in delaying Antarctic cooling at the Bøling onset is discussed in ref. 69. {This estimate uses the onset, rather than the midpoint, of the Bøling–Alerød to Younger Dryas transition in WDC CH4in determining the phasing. {This estimate may be unreliable because of a double peak in WDCd18O at about this time that may reflect local climate, as wel as a unique and pronounced accumulation anomaly starting at about 11.9 kyrBP.