# Abrupt ice-age shifts in southern westerly winds and Antarctic climate forced from the north

Christo Buizert, Michael Sigl, Mirko Severi, Bradley Markle, Justin Wettstein, Joseph A Mcconnell, Joel Pedro, Harald Sodemann, Kumiko Goto-Azuma, Kenji Kawamura, et al.

► To cite this version:

Christo Buizert, Michael Sigl, Mirko Severi, Bradley Markle, Justin Wettstein, et al. Abrupt ice-age shifts in southern westerly winds and Antarctic climate forced from the north. Nature, 2018, 563 (7733), pp.681-685. 10.1038/s41586-018-0727-5. hal-02395936

HAL Id: hal-02395936

https://hal.science/hal-02395936v1

Submitted on 1 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Abrupt Ice Age Shifts in Southern Westerlies and Antarctic Climate Forced from
the North

3 Christo Buizert¹, Michael Sigl², Mirko Severi³, Bradley R. Markle⁴, Justin J. Wettstein¹,⁵,

4 Joseph R. McConnell⁶, Joel B. Pedro⁷,⁸, Harald Sodemann⁵, Kumiko Goto-Azuma⁹, Kenji Kawamura⁹,

5 Shuji Fujita⁹, Hideaki Motoyama⁹, Motohiro Hirabayashi⁹, Ryu Uemura¹⁰, Barbara Stenni¹¹,

6 Frédéric Parrenin¹², Feng He¹,¹³, T.J. Fudge⁴ and Eric J. Steig⁴

7 ¹College of Earth, Ocean and Atmospheric Sciences, Oregon State University, Corvallis OR 97331, USA

8 ²Laboratory of Environmental Chemistry, Paul Scherrer Institute, Villigen, Switzerland

9 ³Department of Chemistry “Ugo Schiff”, University of Florence, Florence, Italy

10 ⁴Department of Earth and Space Science, University of Washington, Seattle WA 98195, USA

11 ⁵Geophysical Institute and Bjerknes Centre for Climate Research, University of Bergen, 5020 Bergen,

12 Norway

13 ⁶Desert Research Institute, Nevada System of Higher Education, Reno, NV 89512, USA

14 ⁷Centre for Ice and Climate, Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark

15 ⁸Antarctic Climate & Ecosystems Cooperative Research Centre, University of Tasmania, Private Bag 80,

16 Hobart, Tasmania 7001, Australia

17 ⁹National Institute for Polar Research, Tachikawa, Tokyo, Japan

18 ¹⁰Department of Chemistry, Biology and Marine Science, University of the Ryukyus, Okinawa, Japan

19 ¹¹Department of Environmental Sciences, Informatics and Statistics, Ca' Foscari University of Venice,

20 30172 Venezia, Italy

21 ¹²Université Grenoble Alpes, CNRS, IRD, IGE, 38000 Grenoble, France

22 ¹³Center for Climatic Research, Nelson Institute for Environmental Studies, University of Wisconsin-

23 Madison, Madison WI 53706, USA
The Southern Hemisphere (SH) mid-latitude westerly winds play a central role in the global climate system via Southern Ocean upwelling¹, carbon exchange with the deep ocean², Agulhas Leakage³, and possibly Antarctic ice sheet stability⁴. Meridional shifts of the SH westerlies have been hypothesized in response to abrupt North Atlantic Dansgaard-Oeschger (DO) climatic events of the last ice age⁵,⁶, in parallel with the well-documented shifts of the intertropical convergence zone (ITCZ)⁷. Shifting moisture pathways to West Antarctica⁸ are consistent with this view, but may represent a Pacific teleconnection pattern⁹. The full SH atmospheric-circulation response to the DO cycle and its impact on Antarctic temperature remain unclear¹⁰. Here we use five volcanically-synchronized ice cores to show that the Antarctic temperature response to the DO cycle can be understood as the superposition of two modes: a spatially homogeneous oceanic “bipolar seesaw” mode that lags Northern Hemisphere (NH) climate by about 200 years, and a spatially heterogeneous atmospheric mode that is synchronous with NH abrupt events. Temperature anomalies of the atmospheric mode are similar to those associated with present-day Southern Annular Mode (SAM) variability, rather than the Pacific South America (PSA) pattern. Moreover, deuterium excess records suggest a zonally coherent migration of the SH westerlies over all ocean basins in phase with NH climate. Our work provides a simple conceptual framework for understanding the circum-Antarctic temperature response to abrupt NH climate change. We provide observational evidence for abrupt shifts in the SH westerlies, which has previously-documented¹⁻³ ramifications for global ocean circulation and atmospheric CO₂. These coupled changes highlight the necessity of a global, rather than a purely North Atlantic, perspective on the DO cycle.
During the glacial DO cycle, abrupt variations in northward heat transport by the Atlantic Meridional Overturning Circulation (AMOC) affect Greenland and Antarctic temperature oppositely (Fig. 1), via an oceanic teleconnection called the bipolar seesaw⁶,¹¹. Antarctica warms during Greenland cold phases (stadials), and cools during Greenland warmth (interstadials), with the gradual nature of Antarctic climate change reflecting buffering by a large heat reservoir¹¹ – likely the global ocean interior⁶. The DO cycle affects atmospheric circulation also; the ITCZ shifts southwards during stadials, and northwards during interstadials⁷. General Circulation Model (GCM) simulations suggest parallel shifts of the SH westerlies⁵,⁶,¹², but the available observational evidence (a deuterium excess record from West Antarctica⁸) cannot distinguish between such shifts and Pacific-only teleconnections⁹. Furthermore, the impact of the atmospheric circulation changes on Antarctic climate remains unknown, and models are inconclusive on this question¹⁰,¹³.

We use water stable isotope ratios, a proxy for site temperature¹⁴, from five Antarctic ice cores: WAIS (West Antarctic Ice Sheet) Divide (WDC), EPICA (European Project for Ice Coring in Antarctica) Dronning Maud Land (EDML), EPICA Dome C (EDC), Dome Fuji (DF) and Talos Dome (TAL). WDC is synchronized to Greenland ice cores at high precision via atmospheric methane (Fig. 1a, b)¹⁵; here we synchronize WDC to the other cores in the 10-57 ka before present (BP) interval using volcanic markers (Methods; Extended Data Figure 1), greatly improving our ability to study the timing of regional Antarctic climate variations relative to Greenland. The Antarctic response to DO events is investigated using a stacking technique, in which 19 individual events are aligned at the midpoint of their abrupt methane transition in the WDC core, and averaged to obtain the shared climatic signal (Methods).
Antarctica cools in response to DO warming (Fig. 2a, b), consistent with the bipolar seesaw theory⁶,¹¹. In the Antarctic mean δ¹⁸O stack, the cooling onset occurs about two centuries after the abrupt NH event, providing validation of earlier results from West Antarctica¹⁵. There is a spatial pattern to the Antarctic response, however. A step-like divergence from the mean signal is seen around t ≈ 0 yr (i.e., synchronous with NH climate), with the interior East Antarctic Plateau sites (DF and EDC) warming, and EDML cooling (Fig. 2c). This instantaneous warming over the Plateau is particularly pronounced at DO events 1, 8, 12 and 14 (Fig. 1d, red curve).

Using principal component analysis (PCA, see Methods), we find that two modes of variability explain over 96% of signal variance in the five stacked records (Fig. 2d). The first principal component (PC1, 83% of variance explained) has the triangular shape of the Antarctic Isotope Maximum events – the classic thermal bipolar seesaw signal¹¹ – with a spatially homogeneous expression (Fig. 2f). The two-century lag behind Greenland warming identifies PC1 as an ocean-propagated response¹⁵.

The second principal component (PC2, 13% of variance explained) is a step-like function with a heterogeneous spatial pattern (Fig. 2g). This mode is very different from the bipolar seesaw. The PC2 response is synchronous with NH warming within precision (28 ± 40 year lag); this timing, and additional evidence presented below, suggest this mode represents an atmospheric teleconnection. The PCA does not necessarily separate physical processes. We assume two underlying teleconnections: oceanic (two-century lag) and atmospheric (synchronous). Some amount of each process is included in PC1, as evident by some immediate warming around t=0. We perform a rotation of the PCA vectors (Methods) to isolate the “purely” oceanic and atmospheric responses (Fig. 2e). The associated estimate of the atmospherically-forced
temperature anomaly (Fig. 2h) is cooling at EDML, warming at DF, EDC and TAL, and a negligible response at WDC; this pattern is robustly reproduced using different methods (Extended Data Fig. 6). The magnitude of the Antarctic atmospheric response is roughly proportional to the Greenland ice core $\delta^{18}$O perturbation (Extended Data Fig. 4).

The Antarctic response to DO cooling is qualitatively similar to the DO warming case. The ocean seesaw warming response is delayed by 226 ± 44 years and the EOF2 spatial pattern has the opposite sign – i.e. additional warming at EDML, and cooling on the interior East Antarctic Plateau (Extended Data Fig. 7). The atmospheric signal over Antarctica is much weaker for the DO cooling case, with PC2 explaining only 9% of variance. This difference is likely due to the fact that DO warmings are more abrupt and of larger magnitude than DO coolings.

To better understand the atmospheric mode, we turn to deuterium excess (d), a proxy for vapor source conditions¹⁶ commonly used to identify changes in atmospheric circulation and vapor transport pathways⁸,¹⁷,¹⁸. In isotope-enabled GCM simulations, Antarctic d is anti-correlated with the Southern Annular Mode (SAM) index⁸,¹⁹. This anti-correlation can be understood conceptually: when the SH westerlies are displaced equatorward (negative SAM phase), Antarctic moisture will originate further north where sea-surface temperature (SST) is higher and relative humidity lower (Extended Data Fig. 8b), both of which act to make d more positive¹⁶. We use the logarithmic definition of deuterium excess ($d_{ln}$), which better preserves isotopic moisture source information than the linear definition⁸,²⁰.

The Antarctic mean $d_{ln}$ response (Figs. 3a, b) lags NH climate by 8 ± 48 years for DO warming, and 9 ± 42 years for DO cooling, respectively, consistent with previous results for WDC⁸. The
observed $d_{\text{In}}$ response is consistent with a shift in the meridional position of the SH westerly winds and vapor origin, such that they move equatorward in response to NH warming, and poleward in response to NH cooling. The timing suggests propagation to the SH high-latitudes via an atmospheric teleconnection. The $d_{\text{In}}$ response is largest for the interior Plateau sites (DF, EDC), possibly because their vapor source areas are more distant from confounding local effects such as the sea-ice edge<sup>21</sup>. The response is weak or absent at EDML; possibly because SH westerlies' variability is relatively weak in the Atlantic sector (Extended Data Fig. 9), or because of regional effects such as wind-driven changes to the sea ice, gyre circulation or Weddell Sea deep convection<sup>22</sup>. Critically, the four cores that do show a clear $d_{\text{In}}$ response collectively sample water vapor from all ocean basins (Extended Data Fig. 8a), suggesting the changes to the SH atmospheric circulation are zonally coherent and involve all ocean basins (rather than just the Pacific basin as demonstrated previously with WDC).

Figure 4a compares the two independent signals we attribute to a change in atmospheric circulation: PC2 of the $\delta^{18}$O response, and the Antarctic mean $d_{\text{In}}$ response. Their time evolution is nearly identical, suggesting they are distinct but consistent manifestations of the atmospheric circulation change. The SAM and Pacific-South American (PSA) pattern are the leading modes of large-scale SH atmospheric variability with strong influence on Antarctic temperature<sup>9,23</sup>. We focus our analysis on East Antarctica, where we infer the largest response. The SAM (Fig. 4b) clearly impacts East Antarctic surface air temperature (SAT) strongly (correlation $|r|$ up to 0.65), and with the correct sign to explain the warming seen at EDC, DF and TAL. East Antarctic warming is seen for a more negative SAM index, driven primarily by anomalous atmospheric heat advection<sup>24</sup> (the observed cooling at EDML is discussed below). The PSA (Fig. 4c), on the
other hand, is not meaningfully correlated with SAT at the East Antarctic sites ($|r|$ at or below 0.15). We further create a synthetic index that is the projection of the atmospheric loadings (Fig. 2h) onto the reanalysis SAT anomaly at the core locations (Methods). The patterns in SAT and geopotential height associated with this index (Fig. 4d) closely resemble those of the SAM, with warming in East Antarctica, and roughly annular geopotential height anomalies.

These tests suggest that the SAM is the closest present-day analog to the temperature response we identify in the ice core record, corroborating our independent evidence from the $d_{ln}$ data. While the PSA pattern may have been active during the DO cycle, it does not dominate the Antarctic response.

Our data-based inferences on the timing and sign of changes to the SH westerlies/SAM are consistent with coupled atmosphere-ocean GCM simulations in which AMOC transitions are induced by North Atlantic fresh-water forcing$^{6,12,25}$. Such model simulations show a positive shift in SAM index in response to AMOC shutdown and vice versa (Fig. 3a, b); this shift is synchronous with the applied forcing within uncertainty (Extended Data Table 1). Our observed atmospheric response is more gradual than the model-simulated SAM shift, possibly because of (multi-decadal) data resolution in some cores and the fact that the $d_{ln}$ signal integrates over a large moisture source area extending to 20°S.

Next, we address differences between the ice-core data and the modern-day correlation pattern (Fig. 4b), most notably at EDML. The reanalysis correlation pattern captures the SAT response to monthly internal SAM variability, representing atmospheric heat advection anomalies$^{24}$. The ice cores, on the other hand, record the response to a persistent long-term
shift in SAM<sup>13,26</sup>, driving changes in SST, stratification and sea ice extent<sup>22,26</sup>. We speculate that
on longer timescales the oceanic influence of the Weddell Sea drives the cooling at EDML, due
to e.g. enhanced sea ice cover<sup>22</sup> and stratification, and a weakening of the wind-driven Weddell
Gyre. The negligible warming response at WDC is consistent with the relatively weak influence
of the SAM in West Antarctic seen in monthly reanalysis (Fig. 4b). Our observations may help
constrain the long-term response to a persistent SAM shift, on which GCMs disagree<sup>13</sup>.

Last, we want to highlight additional structure in the Antarctic δ<sup>18</sup>O stacks that is currently not
part of scientific discourse on interhemispheric climate coupling. Most notably, Antarctic
warming appears to slow down around t=−400 yr (Fig. 2b), forming a secondary change point
that precedes the abrupt DO warming events<sup>27</sup>. Likewise, the rate of Antarctic cooling appears
to increase 200 years prior to the abrupt DO cooling events (Extended Data Fig. 7b). These
secondary change points are subtler than the ones analyzed in this work, have no apparent
corresponding features in Antarctic *d*<sub>ln</sub> or Greenland climate, and remain unexplained.

In conclusion, our results show that Antarctica is influenced by NH abrupt climate change on
two distinct time scales, representing a slow oceanic and a fast atmospheric teleconnection. In
particular, we provide observational evidence for zonally-coherent meridional shifts in the SH
westerly winds in phase with Greenland DO events, and its impact on Antarctic temperature.
Such shifts have implications for global ocean circulation, Southern Ocean upwelling and
productivity, and atmospheric CO<sub>2</sub><sup>1-3</sup>. It is therefore paramount to consider the DO cycle from a
global, rather than a purely North Atlantic perspective.

174
175 References:

176 Marshall, J. & Speer, K. Closure of the meridional overturning circulation through Southern Ocean upwelling. *Nature Geosci* **5**, 171-180 (2012).

178 Toggweiler, J. R., Russell, J. L. & Carson, S. R. Midlatitude westerlies, atmospheric CO₂, and climate change during the ice ages. *Paleoceanography* **21**, PA2005, doi:10.1029/2005pa001154 (2006).

181 Biastoch, A., Boning, C. W., Schwarzkopf, F. U. & Lutjeharms, J. R. E. Increase in Agulhas leakage due to poleward shift of Southern Hemisphere westerlies. *Nature* **462**, 495-U188, doi:10.1038/nature08519 (2009).

184 Pritchard, H. D. et al. Antarctic ice-sheet loss driven by basal melting of ice shelves. *Nature* **484**, 502-505, doi:10.1038/nature10968 (2012).

186 Lee, S. Y., Chiang, J. C., Matsumoto, K. & Tokos, K. S. Southern Ocean wind response to North Atlantic cooling and the rise in atmospheric CO₂: Modeling perspective and paleoceanographic implications. *Paleoceanography* **26** (2011).

189 Pedro, J. B. et al. Beyond the bipolar seesaw: Toward a process understanding of interhemispheric coupling. *Quat. Sci. Rev.* **192**, 27-46, doi:10.1016/j.quascirev.2018.05.005 (2018).

192 Schneider, T., Bischoff, T. & Haug, G. H. Migrations and dynamics of the intertropical convergence zone. *Nature* **513**, 45-53 (2014).

194 Markle, B. R. et al. Global atmospheric teleconnections during Dansgaard-Oeschger events. *Nature Geosci* **10**, 36-40, doi:10.1038/ngeo2848 (2017).

196 Ding, Q., Steig, E. J., Battisti, D. S. & Kuttel, M. Winter warming in West Antarctica caused by central tropical Pacific warming. *Nature Geosci* **4**, 398-403, doi:10.1038/ngeo1129 (2011).

198 Buiron, D. et al. Regional imprints of millennial variability during the MIS 3 period around Antarctica. *Quat. Sci. Rev.* **48**, 99-112, doi:10.1016/j.quascirev.2012.05.023 (2012).

200 Stocker, T. F. & Johnsen, S. J. A minimum thermodynamic model for the bipolar seesaw. *Paleoceanography* **18**, 1087, doi:10.1029/2003pa000920 (2003).

202 Rind, D. et al. Effects of glacial meltwater in the GISS coupled atmosphere-ocean model: 2. A bipolar seesaw in Atlantic Deep Water production. *Journal of Geophysical Research: Atmospheres* (1984-2012) **106**, 27355-27365 (2001).

205 Kostov, Y. et al. Fast and slow responses of Southern Ocean sea surface temperature to SAM in coupled climate models. *Clim. Dyn.* **48**, 1595-1609, doi:10.1007/s00382-016-3162-z (2017).

207 Jouzel, J. et al. Validity of the temperature reconstruction from water isotopes in ice cores. *J. Geophys. Res.* **102**, 26471-26487, doi:10.1029/97jc01283 (1997).

209 WAIS-Divide-Project-Members. Precise interpolar phasing of abrupt climate change during the last ice age. *Nature* **520**, 661-665, doi:10.1038/nature14401 (2015).

211 Merlivat, L. & Jouzel, J. Global climatic interpretation of the deuterium-oxygen 18 relationship for precipitation. *J. Geophys. Res.* **84**, 5029-5033, doi:10.1029/JC084iC08p05029 (1979).

213 Masson-Delmotte, V. et al. GRIP Deuterium Excess Reveals Rapid and Orbital-Scale Changes in Greenland Moisture Origin. *Science* **309**, 118-121 (2005).

215 Masson-Delmotte, V. et al. Abrupt change of Antarctic moisture origin at the end of Termination II. *Proc. Natl. Acad. Sci. U. S. A.* **107**, 12091-12094, doi:10.1073/pnas.0914536107 (2010).

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

219 Uemura, R. et al. Ranges of moisture-source temperature estimated from Antarctic ice cores stable isotope records over glacial-interglacial cycles. *Clim. Past* **8**, 1109-1125, doi:10.5194/cp-8-1109-2012 (2012).
222 21 Sodemann, H. & Stohl, A. Asymmetries in the moisture origin of Antarctic precipitation.

223 *Geophys. Res. Lett.* **36**, L22803, doi:10.1029/2009gl040242 (2009).

224 22 Lefebvre, W., Goosse, H., Timmermann, R. & Fichefet, T. Influence of the Southern Annular

225 Mode on the sea ice–ocean system. *J. Geophys. Res.* **109**, C09005, doi:10.1029/2004JC002403

226 (2004).

227 23 Thompson, D. W. J. & Wallace, J. M. Annular Modes in the Extratropical Circulation. Part I:

228 Month-to-Month Variability*. *J. Clim.* **13**, 1000-1016, doi:10.1175/1520-

229 0442(2000)013<1000:amitec>2.0.co;2 (2000).

230 24 Sen Gupta, A. & England, M. H. Coupled ocean-atmosphere-ice response to variations in the

231 Southern Annular Mode. *J. Clim.* **19**, 4457-4486 (2006).

232 25 Liu, Z. *et al*. Transient Simulation of Last Deglaciation with a New Mechanism for Bolling-Allerod

233 Warming. *Science* **325**, 310-314, doi:10.1126/science.1171041 (2009).

234 26 Ferreira, D., Marshall, J., Bitz, C. M., Solomon, S. & Plumb, A. Antarctic Ocean and Sea Ice

235 Response to Ozone Depletion: A Two-Time-Scale Problem. *J. Clim.* **28**, 1206-1226,

236 doi:10.1175/JCLI-D-14-00313.1 (2014).

237 27 Raisbeck, G. M. *et al*. An improved north–south synchronization of ice core records around the

238 41 kyr 10Be peak. *Clim. Past* **13**, 217-229, doi:10.5194/cp-13-217-2017 (2017).

239 28 Grootes, P. M., Stuiver, M., White, J. W. C., Johnsen, S. & Jouzel, J. Comparison of oxygen

240 isotope records from the GISP2 and GRIP Greenland ice cores. *Nature* **366**, 552-554 (1993).

241 29 Rhodes, R. H. *et al*. Enhanced tropical methane production in response to iceberg discharge in

242 the North Atlantic. *Science* **348**, 1016-1019, doi:10.1126/science.1262005 (2015).

243 30 Dee, D. P. *et al*. The ERA-Interim reanalysis: configuration and performance of the data

244 assimilation system. *Quarterly Journal of the Royal Meteorological Society* **137**, 553-597,

245 doi:10.1002/qj.828 (2011).

246 **Acknowledgements.** This work is funded through the U.S. National Science Foundation, grants

248 ANT-1643394 (to C.B. and J.J.W.), ANT-1643355 (to T.J.F. and E.J.S), AGS-1502990 (to F.H.); the

249 Swiss National Science Foundation grant 200021_143436 (to H.S.); CNRS/INSU/LEFE projects

250 IceChrono and CO2Role (to F.P.); JSPS KAKENHI grants 15H01731 (to K.G.-A., H.M. and M.H.),

251 15KK0027 (to K.K.) and 26241011 (to K.K., S.F. and H.M.); MEXT KAKENHI grant 17H06320 (to

252 K.K., H.M. and R.U.); European Research Council under the European Community’s Seventh

253 Framework Programme (FP7/2007–2013)/ERC grant agreement 610055 (to J.B.P.) ; and the

254 NOAA Climate and Global Change Postdoctoral Fellowship program, administered by the

255 University Corporation for Atmospheric Research (F.H.). We acknowledge high-performance

256 computing support from Yellowstone (ark:/85065/d7wd3xhc) provided by NCAR's
Computational and Information Systems Laboratory, sponsored by the NSF. This research used
resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory,
which is supported by the Office of Science of the U.S. Department of Energy under Contract №
DE-AC05-00OR22725. This is TALDICE publication № 52.

**Author Contributions.** Data analysis by C.B., M.Se., M.Si. and J.J.W.; manuscript preparation by C.B.; volcanic ice core synchronization by M.Si, M.Se., C.B., J.R.M., F.P., S.F. and T.J.F.; GCM simulations and interpretation by F.H., J.B.P. and J.J.W.; Moisture tagging/tracing experiments by B.R.M. and H.S.; Ice core water isotope analysis by K.G-A, K.K., H.M., M.H., R.U., B.S. and E.J.S.; all authors discussed the results and contributed towards improving the final manuscript.

**Author Information.** Reprints and permissions information is available at
www.nature.com/reprints. The authors declare no competing financial interests.
Correspondence and requests for materials should be addressed to Christo Buizert (buizertc@oregonstate.edu).

**Figure 1 | Records of abrupt glacial climate variability.** a, Greenland Summit (average of GISP2 and GRIP²⁸) ice core δ¹⁸O. b, WDC methane²⁹. c, Antarctic 5-core average *d*ln anomaly. d, Antarctic δ¹⁸O anomaly at EDML (blue), the Antarctic Plateau (average of DF and EDC, red) and 5-core average (black); offset for clarity. All records are synchronized to the WAIS Divide WD2014 chronology; Antarctic data are shown as anomalies relative to present. DO interstadial periods marked in grey and numbered; Heinrich stadials marked in blue. Isotope ratios are on the VSMOW (Vienna Standard Mean Ocean Water) scale. Thin curves show records at original
resolution (ranging from ~5 to ~50 years), with the thick lines a moving average (300 and 150 year window for Antarctic and Greenland data, respectively).

**Figure 2** | The Antarctic climate response to DO warming. a, Stack of NGRIP $\delta^{18}$O. b, Stack of Antarctic $\delta^{18}$O at indicated locations, with "Plateau" the average of DF and EDC. c, As in b, but with 5-core mean subtracted. d, First two principal components of the Antarctic $\delta^{18}$O stacks (1500 yr window), with percentage of variance explained (offset for clarity). PC1 is strongly correlated ($r$ = 0.998) to the Antarctic mean. Linear fit to PC1 ($t$ = -400 to $t$=0 interval) is shown to highlight the response around $t$ = 0. e, Rotated PC1 and PC2 vectors representing proposed oceanic and atmospheric modes, with fits from change point analysis (Methods). The oceanic mode lags by 211 ± 42 years; the atmospheric mode by 28 ± 40 years (1$\sigma$ bounds; Extended Data Table 1). f-g, Empirical Orthogonal Functions EOF1 and EOF2 associated with PC1 and PC2 in d, scaled to show the magnitude in permil. h, Spatial pattern associated with the atmospheric mode as shown in e, scaled to permil. Isotope ratios are on the VSMOW scale.

**Figure 3** | Deuterium excess and the SH westerlies. a, DO warming: Greenland $\delta^{18}$O stack (turquoise); 5-core average Antarctic $d_{\text{In}}$ stack (orange with BREAKFIT result, see Extended Data Table 1); SAM index (here the leading principal component of sea level pressure variability south of 20°S) following a freshwater-forced AMOC perturbation in CCSM3 (Community Climate System Model version 3) model simulations (grey with fit from change-point analysis, see Extended Data Table 1). b, as a, but for DO cooling. c, Magnitude of Antarctic $d_{\text{In}}$ response to DO
298 warming in permil. The weak $d_{ln}$ trend before and after the abrupt jump likely reflects the SST

299 of SH vapor source waters following the thermal bipolar seesaw$^{8,11}$. **d**, as **c**, but for DO cooling.

300 Isotope ratios are on the VSMOW scale.

301

302 **Figure 4 | Attribution of the atmospheric mode of Antarctic temperature variability.** **a,**

303 Comparison of PC2 of the 5 Antarctic δ$^{18}$O stacks as in Fig. 2d (pink, left axis) and the Antarctic

304 mean $d_{ln}$ stack as in Fig. 3a (black, right axis). Isotope ratios are on the VSMOW scale. **b,**

305 Correlation between a standardized monthly SAM index and SAT (2-meter temperature) in ERA-

306 Interim$^{30}$ for 1979-2017 (shading, scale bar on right) with superimposed 850 hPa geopotential

307 height regressions (10 m contours) and the ice core atmospheric temperature mode from Fig.

308 2h (circles, scale bar from Fig. 2). Note that regions of anti-correlation are colored red (i.e.

309 warming in response to negative SAM shift). **c**, as panel **b**, but for a standardized PSA index. The

310 SAM and PSA are here taken to be PC1 and PC2 of the 850 hPa geopotential height field south

311 of 20°S, respectively. **d**, as panel **b**, but for a synthetic index of the atmospheric mode created

312 by regressing ERA-Interim SAT anomalies at the ice core sites onto the coefficients in Fig. 2h

313 (Methods).

314
315

# Methods

316 **Volcanic ice core synchronization:** Volcanic reference horizons provide the most precise way to 317 synchronize ice-core age scales³¹⁻³⁵. Within the last decade, progress has been made in 318 volcanically linking the EPICA (European Project for Ice Coring in Antarctica) Dome C (EDC) ice 319 core to the EPICA Dronning Maud Land (EDML) core³¹, the Talos Dome (TAL) core³², and the 320 Dome Fuji (DF) core³³. Here we provide new volcanic stratigraphic links between the WAIS 321 (West Antarctic Ice Sheet) Divide ice core (WDC) and the EDML, EDC and TAL cores, based on 322 pattern matching of volcanic peaks in high-resolution records of either sulfur (WDC) or sulfate 323 (EDML, EDC, TAL). Extended Data Fig. 1b summarizes the various synchronizations, with 324 previously published ones indicated with grey arrows, and the new ones presented here 325 indicated with colored arrows.

326 Volcanic synchronization via sulfur/sulfate and/or ECM (electrical conductivity measurements) 327 records is a commonly-used technique for Antarctic ice cores, and we rely on previously 328 established methods described in detail elsewhere³¹⁻³⁵. Matches are made by identifying 329 sequences of sulfur peaks that have the same relative spacing in both cores³¹⁻³³,³⁵. Additional 330 confidence in the match points comes from approximately proportional acid concentration 331 levels, smooth variations in the resulting annual layer thickness between stratigraphic tie- 332 points, and in some cases a distinctive shape of the common signals in the different ice cores. 333 We identify 773 volcanic ties between WDC and EDML, 396 between WDC and EDC, and 425 334 between WDC and TAL (Source Data).

335 Stratigraphic matching was performed independently by two authors (M.Si. and M.Se.), and 336 then compiled and compared by a third author (C.B.). Both analysts use an iterative approach, 337 in which they first identify the major, unambiguous events. After marking these events (or 338 clusters thereof), they re-align the data sets, and replot them with the marked events now 339 overlapping. At this point, usually several of the smaller events are nearly on top of each other. 340 These events are then marked, and the data is replotted with the newly marked events 341 overlapping, et cetera.

342 We distinguish three cases: “doubly”, “singly”, and “inconsistently” identified events. Doubly 343 identified matches are cases where both workers identify the same stratigraphic match 344 between two cores for a given volcanic event (within a margin of a few cm). Singly identified 345 matches are cases where only one of the two workers identified a stratigraphic match. 346 Inconsistently identified matches are cases where a single volcanic event in one core is linked to 347 two different volcanic events in the other core. Around 99% of the singly identified events were 348 found by M.Si., demonstrating a difference in event detection threshold. For example, in the 349 WDC-EDC synchronization the median volcanic peak sizes in WDC non-sea salt Sulfur (nssS) 350 were 67.7 and 31.4 ppb for the doubly and singly identified matches, respectively, 351 demonstrating that analyst M.Se. is more conservative in assigning match points. For 352 comparison, the background nssS level is 15±4 ppb.
Here, all the doubly identified events are all assumed correct, and retained. In the (relatively rare) case of an inconsistently matched event, the stratigraphic links suggested by both workers are discarded to avoid ambiguity. Sequences of singly identified events are retained only if they occur in between two doubly identified match points. In the case of an inconsistently identified event (which is discarded), all singly identified matches adjacent to it are discarded also, until another doubly identified match is encountered. In all three synchronizations, the vast bulk of the singly identified matches come from the same worker, suggesting a difference in style and event detection threshold between the two workers. We give some examples below of hypothetical sequences of tie points, and how they are dealt with. Let "d", "s" and "i" denote doubly, singly, and inconsistently identified tie points, respectively.

**d-s-s-s-d:** This is a hypothetical series of three s tie points in between two d tie points. Because both workers agree on the tie points on either end, there is no reason to assume the s tie points are incorrect; it simply reflects the fact that one worker (M.Se) is much more conservative in assigning tie points than the other worker (M.Si). Therefore all tie points are retained in the final synchronization (**d-s-s-s-d**). Singly identified ties are retained in such cases irrespective of how many are in the series (for example a series of 10 s ties would be retained if bracketed on either side by a d).

**d-i-d:** In this case the i tie is removed, but the d ties are retained resulting in the sequence d-d in the final synchronization.

**d-s-s-i-s-d:** The i tie is removed in all cases. However, in this example the s tie points occur adjacent to an i tie point, which casts doubt on their reliability. Therefore they are removed together with the inconsistent tie, and only the sequence d-d is retained in the final synchronization.

The matching is described below for the 10-61ka interval of interest (the WDC-EDC synchronization extends back only to 57 ka BP).

In synchronizing WDC and EDC, the two workers identified 473 matches, out of which 103 were doubly, 8 inconsistently, and 362 singly identified (all but two of which by M.Si.). The final selection has retained the 103 doubly identified matches, and 293 of the singly identified matches; 69 singly identified matches were discarded as they were bracketed on either side by an inconsistent match. The 8 inconsistently identified matches all differed by less than 70 cm in EDC (or about 65 years).

In synchronizing WDC and EDML, the two workers identified a total of 793 matches, out of which 247 were doubly identified, 3 matches were inconsistently identified (all adjacent, and not separated by doubly identified events), and 543 were singly identified (all by M.Si.). The final selection has retained the 247 doubly identified matches, and 526 of the singly identified matches; 17 singly identified matches were discarded as they were bracketed on either side by an inconsistent match. The 3 inconsistently identified matches all differed by less than 1.6 m in EDML (or about 65 years).
Talos Dome proved to be the most difficult of the cores to synchronize, presumably because the layers are more strongly compressed in this intermediate depth core. The first attempt at synchronization yielded 55 doubly and 5 inconsistently identified matches, with most of the errors in the 770 to 905 m depth range (up to 5 m offsets in TAL). In light of this inconsistency, the workers reviewed their volcanic ties throughout the core, with particular focus on the problematic interval. In the revised synchronization, the workers identified a total of 437 matches, out of which 253 were doubly identified, 4 matches were inconsistently identified (all adjacent, and not separated by doubly identified events), and 180 were singly identified (all but 9 by M.Si). The final selection has retained the 253 doubly identified matches, and 172 of the singly identified matches; 8 singly identified matches were discarded. The 4 inconsistently identified matches all differed by less than 75 cm in TAL (or about 75 years). So while the final synchronization shows good agreement between both workers, we feel obliged to also report the first, less successful attempt. The TAL depth range that is hardest to synchronize to WDC (between 770 and 905 m TAL depth, or 14.9 to 25.2 ka BP) also yielded no matches to EDC in a published study³² (see yellow triangles at bottom of Extended Data Fig. 1a). Note that this problematic 15-25 ka time interval does not influence the main results of this work – none of the abrupt DO events used in the stacking lie in this interval (DO 2 is excluded from the stacking because the absence of an abrupt CH4 signal precludes precise synchronization to Greenland¹⁵).

The number of doubly and inconsistently identified events provides one way to assess the reliability of the volcanic synchronization; the percentage of inconsistently identified events (out of the pool of doubly and inconsistently identified events) ranges from 1% to 7% for the three cores. Errors tend to occur in clusters of adjacent picks, due to the misidentification of a sequence of peaks; seen in this light the WDC-EDML and WDC-TAL synchronizations each only contain a single inconsistently identified sequence. The observed inconsistencies were always relatively small, and on the order of a few decades. A second method of assessing the reliability is via the redundancy offered by having multiple cores. Whenever three ice cores in Extended Data Fig. 1b are connected via three independent synchronizations (i.e., whenever the arrows form a triangle), this offers the possibility to test the internal consistency of the synchronization. Over the age interval of interest (the last 61 ka) 213 ties had previously been identified between EDML and EDC³⁶, as well as 102 ties between TAL and EDC. This introduces a degree of redundancy that allows testing the internal consistency of the synchronization. EDC, TAL and EDML are volcanically synchronized within the AICC2012 (Antarctic Ice Core Chronology)³⁶, while WDC uses the independent WD2014 time scale³⁷,³⁸. For each volcanic tie point, the difference in WD2014 age minus the AICC2012 age is shown in Extended Data Fig. 1a. EDC and EDML are volcanically synchronized over the last 60ka (blue triangles), and therefore the excellent agreement between the blue (WDC age minus EDML age) and red (WDC age minus EDC age) curves in Extended Data Fig. 1a shows the volcanic framework to be internally consistent. Likewise, for the period 25-42 ka BP and <13 ka BP, where TAL and EDC are volcanically synchronized (yellow triangles), the yellow (WDC age minus TAL age) and red (WDC age minus EDC age) curves agree well, suggesting internal consistency. Beyond 42 ka BP the
431 yellow curve deviates from the blue and red ones, suggesting the TAL core is imperfectly
432 synchronized within the AICC2012 chronology (due to an absence of volcanic ties at the time).

433 All ice cores used in this study were synchronized to the WAIS Divide WD2014 chronology<sup>37,38</sup>.
434 For the four non-WDC cores, we start from their original ice core chronologies; this is the
435 AICC2012 chronology<sup>36</sup> for the EDML, EDC and TAL ice cores, and the DF2006 chronology for
436 the DF core<sup>39</sup>. For each core we add a time-variable offset to the WD2014 chronology that is
437 obtained using linear interpolation between the volcanic tie points. For the EDML, EDC and TAL
438 cores we have direct synchronizations to WDC as described above. For the Dome Fuji core,
439 synchronization was indirect via the EDC core (Extended Data Fig. 1b). Previously, 297 tie points
440 have been identified<sup>33</sup> between EDC and DF in the interval of interest (mean spacing of 173
441 years). These volcanic ties are transferred to WDC using the WDC-EDC synchronization
442 described above.

443 Over the time interval of interest, the offset with the AICC2012 cores (EDML, EDC, TAL) ranges
444 roughly from -330 years to + 430 years, with an average offset of -10 years (with negative
445 values meaning WD2014 is younger than AICC2012). For the DF core the range is from -230 to
446 1884 years, with an average of +739 years.

447 **Uncertainty in volcanic synchronization:** There are two types of uncertainty to consider. First,
448 the volcanic ties themselves may be incorrect, and second, in between the ties the
449 interpolation strategy introduces an error. The first type is difficult to quantify. Either the ties
450 have been correctly identified, and the relative age uncertainty is zero, or the ties are false, and
451 the relative age uncertainty is infinite (i.e., we have learned nothing). Past studies have
452 sometimes assigned Gaussian errors to volcanic tie points; while this is a practical necessity for
453 applications that optimize the fit to multiple age constraints<sup>36,40-42</sup>, it does not reflect the true
454 uncertainty meaningfully and is not applied here.
455

455 We have high confidence in the correctness of the volcanic ties because of the internal
456 consistency of the new volcanic ties with previously published ones (Extended Data Fig. 1a), and
457 the fact that doubly identified ties greatly outnumber the inconsistently identified ones. We
458 have removed inconsistent matches from the synchronization, and we here assume the
459 remaining matches to be correct.

460 The second type of uncertainty is due to interpolation between volcanic ties. This introduces an
461 age uncertainty that depends on the age difference between adjacent ties, L. We estimate the
462 interpolation uncertainty using the layer-counted section of the WD2014 chronology, which
463 goes back to 31.2 ka BP. To estimate the interpolation uncertainty for two volcanic markers
464 that are, say, L=100 years apart, we can randomly pick thousands of 100-year intervals from the
465 WAIS Divide WD2014 chronology. Within each of these, the age evolution deviates from the
466 assumed linear interpolation; 1σ standard deviation of this age deviation is used as the
467 interpolation uncertainty. Typical results are shown in Extended Data Fig. 2a for several values
468 of L. We find that the maximum uncertainty scales as ∝L<sup>2</sup>. Compared to East Antarctica, West
Antarctica receives a larger contribution of its snowfall from storm systems and synoptic activity, making accumulation rates more variable<sup>43,44</sup>; the estimates given here should thus be considered conservative when used in interior Antarctica. The volcanic interpolation uncertainty for the four cores is plotted in Extended Data Fig. 2b. For DF, synchronization to WDC is done via EDC as an intermediary core, and therefore the two synchronization errors are added in quadrature.

In our synchronization we use both the singly and doubly identified tie points, and treat them equally. We acknowledge, however, that the doubly identified ties are more reliable. Therefore, we have repeated the analyses described in the main text using only the doubly identified tie points (as opposed to both singly and doubly identified tie points). We find that the conclusions of this work do not depend on this choice of tie points. Those tests are not shown here, but the alternative chronologies that use only the doubly identified tie points, and alternative versions of the manuscript figures showing those analyses were presented to the reviewers and are available from the corresponding author upon request.

**Water stable isotope data:** A combination of previously published and unpublished ice core water isotope data (δ<sup>18</sup>O and δD) are used in this study. Deuterium excess (d<sub>ln</sub>) is calculated from the δ<sup>18</sup>O and δD isotope ratios using the logarithmic definition of Uemura et al.<sup>20</sup>:

$$d_{\text{ln}} = \ln(1+\delta\text{D}) - 8.47 \ln(1+\delta^{18}\text{O}) + 0.0285 [\ln(1+\delta^{18}\text{O})]^2$$

For WDC we use previously published δ<sup>18</sup>O and δD data<sup>8,15,45</sup>, measured using laser spectroscopy. Data have a typical depth resolution of 1m for the 0 to 2.3 ka BP interval, of 0.5 m for 2.3 to 56 ka BP, and of 0.25 m for 56-68 ka BP; this corresponds to an average time resolution of 17 years for the study period (11 to 61 ka BP). Using the cm-scale CFA (continuous flow analysis) record of WDC δ<sup>18</sup>O instead gives identical results to those presented here<sup>46</sup>.

For EDML we use previously published δ<sup>18</sup>O and δD data<sup>47,48</sup>, measured using conventional IRMS (isotope ratio mass spectrometry). Data have a typical depth resolution of 0.5 m, corresponding to an average time resolution of 24 years for the study period.

For EDC we use previously published δ<sup>18</sup>O and δD data<sup>48,49</sup>, measured using conventional IRMS. Data have a typical depth resolution of 0.55 m, corresponding to an average time resolution of 44 years for the study period.

For DF we use new and published water isotope data<sup>39,50,51</sup>. Two data sets are used. The first is a data set of δ<sup>18</sup>O measured using IRMS in the 300 to 1151 m depth range (10 to 71 ka BP) at 0.5 m resolution<sup>50</sup>. The second is a data set of δ<sup>18</sup>O and δD measured using IRMS in the 550 to 849 m depth range (23 to 45 ka BP); this data set was measured at 0.1 m resolution, and averaged into 0.5 m bins. Note that d<sub>ln</sub> is only available from the second data set, which spans DO/AIM (Antarctic Isotopic Maximum) evens 2 through 11. In the depth range of overlap, δ<sup>18</sup>O data from both data sets are averaged (after correcting the second data set by + 0.213‰ to
account for a calibration offset) to produce the final time series. The combined $\delta^{18}$O record has an average time resolution of 35 years for the study period.

For TAL we use a combination of new and previously published$^{10,52,53}$ data. Bag average, 1 m resolution $\delta^{18}$O and $\delta$D data measured using IRMS are available for the entire core. High-resolution (0.1 m) $\delta^{18}$O data measured using IRMS are available in the 598 to 786 m (10 to 16 ka BP) and 1030 to 1282 m (37 to 65 ka BP) depth ranges. High-resolution $\delta^{18}$O are averaged into 0.5 m bins, and combined with bag-average data for remaining depth intervals. The combined $\delta^{18}$O record has an average time resolution of 50 years for the study period.

**Stacking procedure:** The stacking procedure used to investigate the Antarctic climate response to abrupt DO variability is described in detail elsewhere$^{15}$. In short, the individual Greenland events are aligned at the midpoint of their abrupt $\delta^{18}$O transitions (either DO warming or DO cooling). All Antarctic events (on their volcanically-synchronized WD2014 timescales) are aligned at the midpoints of the abrupt WDC $CH_4$ transitions. We then average over events to obtain the shared climatic signal; to derive the north-south phasing we use the independently established $CH_4$ delay of 56 ± 19 years (1$\sigma$) behind $\delta^{18}$O in Greenland$^{54}$.

For one of the abrupt events (DO-10 / AIM-10), improved inter-polar synchronization data are available from $^{10}$Be variations during the Laschamp event (41 ka BP) between the Greenland NGRIP core, and the Antarctic EDC and EDML cores$^{27}$; these timing constraints were incorporated into our stacking procedure (DO-10 only). The EDML, EDC, DF and TAL ice cores have much higher $\Delta$age uncertainty and lower resolution $CH_4$ records than WDC, therefore the north-south phasing precision cannot be improved by considering $CH_4$ synchronization for those cores also.

In this work we consider DO-0 (i.e. the YD-Holocene transition) through DO-16; DO-17 falls outside the volcanic synchronization for the EDC and DF cores. DO-2 is omitted due to the absence of a clear $CH_4$ response, precluding synchronization. All stacked records are shown in Extended Data Fig. 3.

Uncertainty in the timing of the stacked records comes from these sources: (1) The gas age-ice age difference ($\Delta$age) in the WDC core$^{37}$; (2) DO midpoint detection in the abrupt NGRIP $\delta^{18}$O record; (3) DO midpoint detection in the abrupt WDC $CH_4$ record; (4) Interpolation of the WDC chronology between tie points; (5) The climatic lag of atmospheric $CH_4$ behind Greenland $\delta^{18}$O of 56 ± 19 years (1$\sigma$)$^{54}$; (6) Trend analysis in the BREAKFIT$^{55}$ and RAMPFIT$^{56}$ fitting routines; and (7) volcanic synchronization onto the WD2014 chronology.

The combined uncertainty due to the first 5 items was assessed previously using a Monte-Carlo routine, suggesting a 1$\sigma$ timing uncertainty of 38 and 41 years in the WDC stacks for DO warming and DO cooling, respectively$^{15}$. The uncertainty in the trend analysis is given in Extended Data Table 1. The uncertainty in the volcanic synchronization is shown in Extended
Data Fig. 2; averaged over the stacked events the 1σ synchronization uncertainty is 0 years at WDC, 2 years at EDML, 6 years at EDC, 13 years at DF, and 2 years at TAL.

By stacking only the most prominent DO/AIM events (those following Heinrich events; i.e., DO-4, 8, 12, 14 and 17) or just the minor ones (the remainder), we find that the magnitude of the atmospherically-forced Antarctic response is larger for the former, suggesting proportionality with the climate perturbation (Extended Data Fig. 4a-4f). Proportionality of the atmospheric response is further seen for individual events (Extended Data Fig. 4g); see the figure caption for details.

**Principal Component Analysis:** Principal component analysis (PCA) allows different climatic modes to be identified in (paleoclimatic) time series from different locations⁵⁷. Here we perform PCA on the stacked δ¹⁸O records at the 5 individual sites using the MATLAB function "pca". The δ¹⁸O stacks in the main manuscript combine 17 individual events – all those that fall within the volcanic synchronization interval. To get a sense for the sensitivity to including or excluding individual events, we perform additional experiments in which we stack only *n* events (rather than all 17). The *n* events are randomly sampled without replacement (i.e. any given event cannot be picked twice for each stacking). We then perform PCA of these stacked records (the same events are stacked for each of the cores), and standardize the PC vectors by taking the z-score (or standard score). Extended Data Fig. 5 shows typical results for *n*=2 and *n*=8, where for each *n* we repeat the experiment 50,000 times to obtain reliable statistics. Extended Data Figs. 5a and 5b show PC1 and PC2, respectively, with the solid line showing the mean of 50,000 experiments, and the shaded envelope the associated +/- 1σ standard deviation. Extended Data Fig. 5c shows a histogram of the percentage of variance explained by each of the modes.

We find that even by stacking as few as just *n*=2 events the method can, on average, identify the oceanic and atmospheric components described in the main manuscript. Perhaps not surprisingly, we find that when including fewer events the signal-to-noise ratio decreases: with fewer events the estimated signal amplitude is smaller in both PC1 and PC2, and the uncertainty envelope is wider. As more events are included in the stacking, the percentage of variance explained by PC1 increases as the coherence between the various Antarctic cores increases due to an improved signal-to-noise ratio.

Analysis in the main manuscript uses a 1500-year window (centered around *t* = 100 yr), which is chosen because it corresponds to the recurrence time of the shortest DO cycles⁵⁸⁻⁶¹. The variance explained by the oceanic and atmospheric modes depends on the window length, as shown in Extended Data Fig. 6a. For window lengths exceeding ~750 years, the “oceanic” mode explains most of the variance (PC1), with the “atmospheric” mode explaining less. However, at short window lengths (<750 year) the atmospheric mode explains most of the variance, making it PC1. Comparing PC1 at 400-year window to PC2 at 2000-year window (Extended Data Fig. 6b) illustrates the ability of the PCA to identify the atmospheric mode at different window lengths.
The cross-over behavior (i.e. the atmospheric mode shifts from being PC2 to PC1 as a function of window length) is due to the fact that the signal variance of the step-like atmospheric mode occurs chiefly within the $t = 0$ to $t = 100$ interval; signal variance of the oceanic mode depends strongly on the window length, due to its gradual nature.

Lag time analysis of PC1 and PC2 using the RAMPFIT⁵⁶ and BREAKFIT⁵⁵ routines confirms the cross-over behavior. At short window-length (<700 years), PC1 is characterized by the instantaneous /fast response of the hypothesized atmospheric teleconnection, whereas at large window-length (>800 years) it shows the 200-year delayed hypothesized oceanic teleconnection. PC2 shows the opposite behavior, though note that for PC2 at 400-year window length no meaningful solution can be found with either routine. At window lengths >700 years the lag times remain stable and vary only within the uncertainty bound.

Unless specified otherwise, a 1500-year window is used in this work.

### Robustness of the atmospheric spatial pattern:

The spatial pattern we identify for the atmospheric mode is one of the main results of this work, and we here test its robustness. In Extended Data Figs. 6d—6f we compare three alternative ways of identifying the spatial pattern; the Pearson correlation coefficient ($r$) between the shown pattern and the pattern identified in the main text (EOF 2 at 1500 yr window) is given for each. Details on the three methods are given in the figure caption. Correlation coefficients ranging from $r = 0.92$ to $r = 0.999$ suggest that identification of the atmospheric pattern is both qualitatively and quantitatively robust.

**Rotated PCA vectors:** PCA aims to explain the largest amount of variance, whereas our goal is to distinguish between the oceanic and atmospheric modes. While PC1 is clearly dominated by the 200-year-delayed oceanic bipolar seesaw (Fig. 2d), it appears that PC1 also captures some abrupt warming around $t = 0$, apparently due to the fact that atmospherically-induced warming is more prevalent over Antarctica than cooling (NH DO warming case). We therefore construct (admittedly somewhat subjective) oceanic and atmospheric response functions (Fig. 2e), which are derived from the principal components in the following way. Let PC1 and PC2 be the first two principal components through time; these vectors are perpendicular. We let the atmospheric response ATM be identical to PC2. The oceanic response OCE is found by rotating the PC1 vector in the plane spanned by vectors PC1 and PC2 over an angle of -13 degrees:

$$OCE = \cos(-13^{\circ}) \times PC1 + \sin(-13^{\circ}) \times PC2$$

The rotation angle of -13 degrees is picked such that the $\delta^{18}$O shift at $t = 0$ is minimized. The spatial pattern associated with the ATM response (Fig. 2h) is found by multiple linear regression of the $\delta^{18}$O stacks at the individual sites to OCE and ATM using the MATLAB function "regress".

**SAM, PSA and synthetic “atmospheric” indices.** The SAM and PSA indices were calculated from climate model and reanalysis data as respectively the first and second mode of variability in
PCA/EOF analysis (MATLAB function "pca") of sea-level pressure (model) or 850-hPa
geopotential height (reanalysis) south of 20°S, after subtracting the long term mean and scaling
the anomalies by the square root of the cosine of latitude to account for the decreased surface
area closer to the poles. A synthetic index of the inferred atmospheric teleconnection is created
by projecting reanalysis SAT anomalies at the ice core locations onto the atmospheric loading
coefficients from Fig. 2h, the spatial signature of which is shown in Fig. 4d. Mathematically, this
projection is done by multiplying the SAT anomalies with the loading coefficients, and summing
over them at each monthly time step.

It is worth noting that all these inferences with respect to reanalysis data are reliant upon just
five atmospheric loading coefficients from the limited (from the perspective of the large-scale
circulation) spatial domain of Antarctica, so it is difficult to rigorously exclude non-SAM
atmospheric influences from the temperature pattern alone.

In Extended Data Fig. 9 we compare variability in the SH westerly winds in the glacial CCSM3
(Community Climate System Model version 3) climate simulations<sup>25,62</sup> to the ERA-Interim
reanalysis<sup>30</sup>. Both show greater variability in the Indian and Pacific sectors, relative to the
Atlantic sector. Compared to ERA-interim, the CCSM3 SAM is more zonal/annular in its
structure; the CCSM3 westerlies also appear to have smaller variability (some of this difference
could be due to comparing annual mean with decadal mean data). In CCSM3 the forced
response of the SH westerlies (at 19ka in response to increased North Atlantic freshwater, right
panel) is very similar to the internal variability of the SH westerlies (middle panel), suggesting
that the change in the SH atmospheric circulation induced by freshwater forcing in the North
Atlantic is analogous to the existing mode of internal variability.

To estimate the magnitude of SAM shifts of the DO cycle, we analyze the changes in central
East Antarctica where the signal is largest and most consistent with the present-day observed
relationship (Fig. 4). Using an isotope sensitivity of 0.8 ‰ K<sup>-1</sup>, the atmospherically-induced
temperature anomaly in central East Antarctica (DF, EDC) is in the range of 0.20°C to 0.45°C
during a DO-warming (lower and upper bound reflecting typical minor and major DO/AIM
events, respectively). Regression of ERA-interim 2m temperature at DF and EDC to the SAM
index shows a slope of around -1.2°C per unit of normalized SAM (the normalized SAM index
time series has a standard deviation of 1), implying a shift in SAM index of around 0.2 to 0.4
normalized (modern-day) SAM units (rounded to one significant figure). This estimate assumes
(1) a linear isotope-temperature response using the modern-day spatial slope, and (2) that the
monthly SAM-SAT regression from monthly internal SAM variability also applies to persistent
SAM shifts during the glacial; both assumptions are subject to uncertainties that we do not
address here. The CCSM3 model simulates a persisting SAM shift of the same magnitude as the
internal SAM variability in decadal averaged data (Fig.3); because internal SAM variability will
be larger in monthly than in decadally-averaged data, the model and data-based estimates may
be in agreement. The reanalysis time period is too short to derive robust estimates of decadally
averaged SAM variability.
**GCM simulations.** We used the TraCE-21k transient climate model simulations done with the Community Climate System Model version 3 (CCSM3)²⁵,⁶²⁻⁶⁵. AMOC collapse and resumption are triggered using freshwater forcing in the North Atlantic. The “DGL19k” run is used for the AMOC collapse, the “DGL-overshoot-C” run for the AMOC resumption case⁶⁴.

**Moisture origin analysis.** We use two separate experiments that trace moisture origin of the precipitation at the coring sites.

The first method uses a Lagrangian moisture source diagnostic⁶⁶ based on a previously published data set⁶⁷. Using the winds, temperature and humidity of the ERA-Interim reanalysis data set³⁰ covering the years 1980-2013, 5 million air parcel trajectories have been calculated covering the global atmosphere at 6-hourly resolution using the Lagrangian particle dispersion model FLEXPART⁶⁸. From the analysis of specific humidity changes over time along the air parcel trajectories, moisture sources have been identified whenever specific humidity in the air parcels increased more than a threshold value of 0.1 g kg⁻¹ over a 6-hour period. The fractional contribution of each moisture source to the final precipitation at the target location (an area of 300×300 km centered on each ice core site) is obtained from calculating the amount contributed by a moisture source to the humidity already present in an air parcel. During precipitation en route, the previous sources' contributions are proportionally discounted. This results in a quantitative estimate of the contribution of surface areas to the precipitation in the target region in units of evaporation (water mass per unit area per time), including their position in terms of latitude and longitude. The moisture source contributions for each site and precipitation event have been composited into mass-weighted annual mean values and scaled with respect to the total amount of water deposited at the target region.

The second method uses water tagging in a 50-year simulation with the Community Atmosphere Model (CAM), with prescribed seasonally varying SST and modern boundary conditions. The experiment is set up to evaluate the meridional moisture source distribution, with further details and figures given in Ref. ⁸. In short, evaporation is tagged in 11 bins. One bin is the Antarctic continent (re-evaporation) and ocean and ice shelfs south of 70°S; 10 are zonal bins of 5° latitude each, ranging from 20°S to 70°S. For each core, the moisture source distribution is found by evaluating the relative contributions from each of the bins for the last 30 years of the run. Two moisture source distributions were created, one for all year in which the mean annual SAM index was positive, and one for all years in which it was negative (Extended Data 8b).

**Change point detection.** We use two well-documented and widely-used change-point detection methods: BREAKFIT⁵⁵ and RAMPFIT⁵⁶, with results given in Extended Data Table 1. The choice of which one to use is based on the shape of the time series x(t). BREAKFIT finds a single change-point, and fits a linear slope on either side; these features make it suitable for the oceanic mode / PC1 in the main manuscript. RAMPFIT finds two change-points; it is assumed x(t) has constant value x₁ for t<t₁, then a ramp up or down until it reaches value x₂ at time t₂, after which it remains constant at value x₂ for t>t₂. These features make RAMPFIT suitable for the
atmospheric mode / PC2 in the main manuscript. To find the two change-points in the $d_{ln}$ stacks, we applied both the RAMPFIT algorithm, and the BREAKFIT algorithm twice (once for each change-point). The results are comparable (Extended Data Table 1), and in the main text we report the average of both methods.

## Methods References:

700 31 Severi, M. et al. Synchronisation of the EDML and EDC ice cores for the last 52 kyr by volcanic signature matching. *Clim. Past* 3, 367-374, doi:10.5194/cp-3-367-2007 (2007).

702 32 Severi, M., Udisti, R., Becagli, S., Stenni, B. & Traversi, R. Volcanic synchronisation of the EPICA-DC and TALDICE ice cores for the last 42 kyr BP. *Clim. Past* 8, 509-517, doi:10.5194/cp-8-509-2012 (2012).

705 33 Fujita, S., Parrenin, F., Severi, M., Motoyama, H. & Wolff, E. W. Volcanic synchronization of Dome Fuji and Dome C Antarctic deep ice cores over the past 216 kyr. *Clim. Past* 11, 1395-1416, doi:10.5194/cp-11-1395-2015 (2015).

708 34 Sigl, M. et al. Insights from Antarctica on volcanic forcing during the Common Era. *Nature Clim. Change* 4, 693-697, doi:10.1038/nclimate2293 (2014).

710 35 Parrenin, F. et al. Volcanic synchronisation between the EPICA Dome C and Vostok ice cores (Antarctica) 0-145 kyr BP. *Clim. Past.* 8, 1031-1045 (2012).

712 36 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. Past* 9, 1733-1748, doi:10.5194/cp-9-1733-2013 (2013).

715 37 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 difference. *Clim. Past.* 11, 153-173, doi:10.5194/cp-11-153-2015 (2015).

718 38 Sigl, M. et al. The WAIS Divide deep ice core WD2014 chronology – Part 2: Annual-layer counting (0-31 ka BP). *Clim. Past* 12, 769-786, doi:10.5194/cp-12-769-2016 (2016).

720 39 Kawamura, K. et al. Northern Hemisphere forcing of climatic cycles in Antarctica over the past 360,000 years. *Nature* 448, 912-916, doi:10.1038/nature06015 (2007).

722 40 Bazin, L. et al. An optimized multi-proxy, multi-site Antarctic ice and gas orbital chronology (AICC2012): 120&ndash;800 ka. *Clim. Past* 9, 1715-1731, doi:10.5194/cp-9-1715-2013 (2013).

724 41 Parrenin, F. et al. IceChrono1: a probabilistic model to compute a common and optimal chronology for several ice cores. *Geosci. Model Dev.* 8, 1473-1492, doi:10.5194/gmd-8-1473-2015 (2015).

727 42 Lemieux-Dudon, B. et al. Consistent dating for Antarctic and Greenland ice cores. *Quat. Sci. Rev.* 29, 8-20, doi:10.1016/j.quascirev.2009.11.010 (2010).

729 43 Fudge, T. J. et al. Variable relationship between accumulation and temperature in West Antarctica for the past 31,000 years. *Geophys. Res. Lett.* 43, 2016GL068356, doi:10.1002/2016gl068356 (2016).

732 44 Fortuin, J. & Oerlemans, J. Parameterization of the annual surface temperature and mass balance of Antarctica. *Ann. Glaciol* 14, 78-84 (1990).

734 45 WAIS Divide Project Members. Onset of deglacial warming in West Antarctica driven by local orbital forcing. *Nature* 500, 440-444, doi:10.1038/nature12376 (2013).

736 46 Jones, T. R. et al. Water isotope diffusion in the WAIS Divide ice core during the Holocene and last glacial. *Journal of Geophysical Research: Earth Surface*, 2016JF003938, doi:10.1002/2016JF003938 (2017).
739 47 EPICA Community Members. One-to-one coupling of glacial climate variability in Greenland and Antarctica. *Nature* **444**, 195-198, doi:10.1038/nature05301 (2006).

741 48 Stenni, B. *et al.* The deuterium excess records of EPICA Dome C and Dronning Maud Land ice cores (East Antarctica). *Quat. Sci. Rev.* **29**, 146-159, doi:10.1016/j.quascirev.2009.10.009 (2010).

742 49 Jouzel, J. *et al.* Orbital and millennial Antarctic climate variability over the past 800,000 years. *Science* **317**, 793-796, doi:10.1126/science.1141038 (2007).

743 50 Kawamura, K. *et al.* State dependence of climatic instability over the past 720,000 years from Antarctic ice cores and climate modeling. *Science Advances* **3**, doi:10.1126/sciadv.1600446 (2017).

748 51 Uemura, R. *et al.* Asynchrony between Antarctic temperature and CO2 associated with obliquity over the past 720,000 years. *Nature Communications* **9**, 961, doi:10.1038/s41467-018-03328-3 (2018).

751 52 Stenni, B. *et al.* Expression of the bipolar see-saw in Antarctic climate records during the last deglaciation. *Nat. Geosci.* **4**, 46-49, doi:10.1038/ngeo1026 (2011).

753 53 Landais, A. *et al.* A review of the bipolar see-saw from synchronized and high resolution ice core water stable isotope records from Greenland and East Antarctica. *Quat. Sci. Rev.* **114**, 18-32, doi:10.1016/j.quascirev.2015.01.031 (2015).

756 54 Baumgartner, M. *et al.* NGRIP CH4 concentration from 120 to 10 kyr before present and its relation to a δ15N temperature reconstruction from the same ice core. *Clim. Past* **10**, 903-920, doi:10.5194/cp-10-903-2014 (2014).

759 55 Mudelsee, M. Break function regression. *The European Physical Journal-Special Topics* **174**, 49-63 (2009).

761 56 Mudelsee, M. Ramp function regression: a tool for quantifying climate transitions. *Computers & Geosciences* **26**, 293-307 (2000).

763 57 Clark, P. U., Pisias, N. G., Stocker, T. F. & Weaver, A. J. The role of the thermohaline circulation in abrupt climate change. *Nature* **415**, 863-869, doi:10.1038/415863a (2002).

765 58 Grootes, P. M. & Stuiver, M. Oxygen 18/16 variability in Greenland snow and ice with 10-3- to 105-year time resolution. *J. Geophys. Res.* **102**, 26455-26470, doi:10.1029/97jc00880 (1997).

767 59 Alley, R. B., Anandakrishnan, S. & Jung, P. Stochastic resonance in the North Atlantic. *Paleoceanography* **16**, 190-198, doi:10.1029/2000pa000518 (2001).

769 60 Buizert, C. & Schmittner, A. Southern Ocean control of glacial AMOC stability and Dansgaard-Oeschger interstadial duration. *Paleoceanography* **30**, 2015PA002795, doi:10.1002/2015pa002795 (2015).

772 61 Ditlevsen, P. D., Kristensen, M. S. & Andersen, K. K. The recurrence time of Dansgaard-Oeschger events and limits on the possible periodic component. *J. Clim.* **18**, 2594-2603, doi:10.1175/jcli3437.1 (2005).

775 62 He, F. *et al.* Northern Hemisphere forcing of Southern Hemisphere climate during the last deglaciation. *Nature* **494**, 81-85, doi:10.1038/nature11822 (2013).

777 63 Liu, Z. *et al.* Younger Dryas cooling and the Greenland climate response to CO2. *Proc. Natl. Acad. Sci. U. S. A.* (2012).

779 64 He, F. *Simulating Transient Climate Evolution of the Last deglaciation with CCSM3*, University of Wisconsin-Madison 2010, (2010).

781 65 Pedro, J. B. *et al.* The spatial extent and dynamics of the Antarctic Cold Reversal. *Nature Geosci* **9**, 51-55, doi:10.1038/ngeo2580 (2016).

783 66 Sodemann, H., Schwierz, C. & Wernli, H. Interannual variability of Greenland winter precipitation sources: Lagrangian moisture diagnostic and North Atlantic Oscillation influence. *J. Geophys. Res.* **113**, D03107, doi:10.1029/2007JD008503 (2008).
786 67 Läderach, A. & Sodemann, H. A revised picture of the atmospheric moisture residence time.
787 Geophys. Res. Lett. 43, 924-933, doi:10.1002/2015GL067449 (2016).
788 68 Stohl, A., Forster, C., Frank, A., Seibert, P. & Wotawa, G. Technical note: The Lagrangian particle
dispersion model FLEXPART version 6.2. Atmos. Chem. Phys. 5, 2461-2474, doi:10.5194/acp-5-
789 790 2461-2005 (2005).
791 69 Locarnini, R. A. et al. in NOAA Atlas NESDIS 73 (eds S. Levitus & A. Mlshonov) (2013).
792 70 Kalnay, E. et al. The NCEP/NCAR 40-year reanalysis project. Bulletin of the American
793 meteorological Society 77, 437-471 (1996).

794
795 **Data Availability.** Source Data (WDC sulfur data, volcanic tie points and water isotope data on
796 synchronized chronologies) and derived products (stacks, PCA results, etc.) are available with
797 the online version of the paper, and via the NOAA paleoclimate data archive.

798
799 **Code Availability.** MATLAB code used for the stacking procedure can be found in the
800 supplementary information of ref.¹⁵, and is available from the corresponding author upon
801 request.

802
803 **Extended Data Figure 1 | Volcanic synchronization of Antarctic ice cores.** a, Age offset
804 between the WD2014³⁷,³⁸ (WDC) and AICC2012³⁶ (TAL, EDML, EDC) age scales, with each dot
805 representing a volcanic tie point. Yellow and blue triangles denote the timing of TAL-EDC and
806 EMDL-EDC volcanic ties³²,³⁶, respectively. b, Overview of synchronizations between the ice
807 cores used in this study. Previously published synchronizations are in grey; synchronization
808 presented here are in color. Synchronizations within Antarctica are based purely on volcanic
809 links; synchronization between WDC and NGRIP (Greenland) are based on atmospheric CH₄
810 (green arrow).

811
812 **Extended Data Figure 2 | Age uncertainty due to volcanic synchronization.** a, Interpolation
813 uncertainty (1σ) for 4 different values of *L* (the spacing between two adjacent volcanic tie
814 points), based on the layer-counted WD2014 age scale³⁸. b, Interpolation uncertainty in
815 synchronizing the 4 ice cores to WAIS Divide. Grey vertical lines give the timing of DO events.

816
817 **Extended Data Figure 3 | Site-specific stacks of δ¹⁸O and dln.** a, Stack of NGRIP δ¹⁸O (teal, left
818 axis) and WDC CH₄ (green, right axis) during DO warming. b, As in a, but for DO cooling. c, Stack
819 of Antarctic δ¹⁸O at indicated locations (see color legend) during DO warming. d, As in c, but for
820 DO cooling. **e**, Stack of Antarctic *d*<sub>ln</sub> at indicated locations (see color legend) during DO

821 warming. **f**, As in **e**, but for DO cooling. Isotope ratios are on the VSMOW scale.

822
823

**Extended Data Figure 4** | Proportionality of the atmospheric response. Panels a-f compare stacks of just the major DO/AIM events (those following Heinrich events, namely DO/AIM0, DO/AIM1, DO/AIM4, DO/AIM8, DO/AIM12 and DO/AIM14, left panels), and just the minor DO/AIM events (the remainder, right panels). a,b, stacks of NGRIP $\delta^{18}$O (left axis) and CH4 (right axis). c,d, Stacks of Antarctic $\delta^{18}$O at indicated locations. e,f, Same as panels c and d, but with Antarctic mean subtracted. g, Proportionality of the atmospheric response for individual events (numbered). NGRIP event size is found via regression of individual NGRIP events to the multi-event NGRIP $\delta^{18}$O stack normalized to unit variance (Fig. 2a). The Antarctic atmospheric response is found via multiple linear regression of single-site, individual events to the atmospheric and oceanic modes (Fig. 2e). Shown are average (dots) and standard deviation (grey error bars) of the response at EDC, DF and EDML (the cores with the strongest atmospheric response); the EDML response is multiplied by -1 as it has the opposite sign of the response at DF and EDC. Red and blue dots denote the major and minor DO/AIM events, respectively. Isotope ratios are on the VSMOW scale.

838

**Extended Data Figure 5** | Reducing the number of events in the $\delta^{18}$O stacks. a, PC1 when stacking 2 or 8 randomly selected events; thick line and shaded area represent the mean and +/- 1$\sigma$ standard deviation of 50,000 runs, respectively. Vertical yellow bars denote the 200-year period after the abrupt DO event at $t = 0$. b, as panel a, but for PC2. c, Histogram of percentage of signal variance explained by PC1 and PC2, when stacking 2 or 8 randomly selected events. Color coding as in panels a and b.

845

**Extended Data Figure 6** | Robustness of the atmospherically-forced warming pattern. a, Principal Component analysis as function of window length, with percent variance explained by PC1 and PC2. b, Comparison of PC1 at 400 year window length, to PC2 at 2000 year window length, to show the cross-over of the atmospheric response (i.e. from PC1 to PC2) as a function of window length. c, Lag time of the Antarctic PC1 and PC2 response as a function of window length, assessed using BREAKFIT$^{55}$ and RAMPFIT$^{56}$ routines (see Methods for explanation of which one is used). d, EOF1 at 400 year window length scaled to permil units. e, EOF2 at 2000 year window length scaled to permil units. f, Slope of linear fit to $\delta^{18}$O stacks in the $t = 0$ to $t = +200$ year interval, shown as the change in permil during the 200 years. Isotope ratios are on the VSMOW scale.

856

**Extended Data Figure 7** | The Antarctic climate response to DO cooling. a, Stack of NGRIP $\delta^{18}$O. b, Stack of Antarctic $\delta^{18}$O at indicated locations. c, As in b, but with Antarctic mean subtracted. d, First two principal components of the Antarctic $\delta^{18}$O stacks, with percentage of
variance explained (offset for clarity) Lines show the BREAKFIT (PC1) and RAMPFIT (PC2) fits. **e-f**, Empirical Orthogonal Functions EOF1 and EOF2 associated with PC1 and PC2 in **d**, scaled to show the magnitude in permil. Isotope ratios are on the VSMOW scale.

**Extended Data Figure 8 | Moisture sources of Antarctic ice core and the SAM.** **a**, Mass-weighted probability distribution functions (PDFs) of Antarctic moisture sources for the five ice cores of interest (5×10⁻⁵ deg⁻² contour lines; the area integrated PDF equals 1). Distributions are calculated from reanalysis data³⁰,⁶⁷ using a Lagrangian source diagnostic²¹,⁶⁶ (Methods). **b**, Parallels are plotted in 15° increments of latitude; meridians in 45° increments of longitude. **c**, Sea surface temperature (SST, black)⁶⁹ and relative humidity (RH, grey)⁷⁰ as a function of latitude. The colored curves give the latitudinal source distribution during negative SAM phase (solid curves, SAM index < 0) and positive SAM phase (dashed curves, SAM index > 0). The solid and open dots show the first moment of the source distribution during negative and positive SAM phase, respectively. Note that during a positive SAM phase moisture sources for all core locations are located closer to the Antarctic continent. Source distribution data were obtained using water tagging experiments⁸ in the Community Atmosphere Model (Methods).

**Extended Data Figure 9 | Southern Annular Mode-like variability in zonal near-surface winds.** **a**, ERA-Interim reanalysis (1979-2016 annual means) zonal wind speed at 10 m height regressed onto SAM index (here the first principal component of sea level pressure variability south of 20°S), expressed in units of m s⁻¹ per standard deviation in the index. **b**, Same as panel **a**, but for internal variability in the CCSM3 SynTraCE model simulation²⁵,⁶² during glacial climate prior to freshwater forcing of Heinrich stadial 1 (19.5ka to 19.01ka BP, decadal means). **c**, Same as panel **a**, but for the response to North Atlantic freshwater forcing in the CCSM3 SynTraCE model (19.1ka to 18.9ka BP, decadal means, with the freshwater forcing applied at 19ka).

**Extended Data Table 1 | Change-point analysis on Antarctic response.** Change-points are found using the BREAKFIT⁵⁵ or RAMPFIT⁵⁶ routine as indicated; the parameter *t*<sub>OCE</sub> is the single change-point of the oceanic mode; *t*<sub>ATM1</sub> and *t*<sub>ATM2</sub> are the two change-points of the atmospheric mode, representing the onset and ending of the shift, respectively (Methods). Stated uncertainties give the 1σ value in the fitting routine only, found using a Monte Carlo moving block bootstrap analysis with 1000 iterations⁵⁵,⁵⁶; the full uncertainty in the combined interpolar CH<sub>4</sub> synchronization and stacking procedure¹⁵ is estimated to be around 40 years (1σ), which can mostly be attributed to uncertainty in the WAIS Divide ice age-gas age difference Δ*age*³⁷. The δ<sup>18</sup>O PC2 and δ<sup>18</sup>O atmospheric modes are identical.
<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

NGRIP

<figure><img src="image_4.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_8.png" /></figure>

EOF1

<figure><img src="image_9.png" /></figure>

EOF2

<figure><img src="image_10.png" /></figure>

Atmospheric

$\delta^{18}\text{O} (\text{‰})$

Time (year)
<figure><img src="image_11.png" /><figcaption>Time (year)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>d</figcaption></figure>
Antarctic δ¹⁸O PC2

<figure><img src="image_14.png" /></figure>

Antarctic d₁ₙ (%)

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>