# Poleward shift in the Southern Hemisphere westerly winds synchronous with the deglacial rise in CO2

William Gray, Casimir de Lavergne, Robert Wills, Laurie Menviel, Paul

Spence, Mark Holzer, Masa Kageyama, Elisabeth Michel

► To cite this version:

William Gray, Casimir de Lavergne, Robert Wills, Laurie Menviel, Paul Spence, et al. Poleward shift in the Southern Hemisphere westerly winds synchronous with the deglacial rise in CO2. 2023. [hal-03861919](https://hal.science/hal-03861919)

HAL Id: hal-03861919

https://hal.science/hal-03861919

Preprint submitted on 7 Jun 2023

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
This is a non-peer reviewed pre-print submitted to EarthArXiv. This manuscript is currently
under review at the journal *Nature*.

Subsequent versions of this manuscript may have slightly different content.

We welcome feedback; please contact William Gray (william.Gray@lsce.ipsl.fr).
# Poleward shift in the Southern Hemisphere westerly winds synchronous with the deglacial rise in CO₂

William R. Gray¹*, Casimir deLavergne², Robert C. J. Wills³, Laurie Menviel⁴, Paul Spence⁵, Mark Holzer⁶, Masa Kageyama¹, Elisabeth Michel¹

¹Laboratoire des Sciences du Climat et de l'Environnement (LSCE/IPSL), Université Paris-Saclay, Gif-sur-Yvette, France

²LOCEAN Laboratory, Sorbonne Université-CNRS-IRD-MNHN, Paris, France

³Department of Atmospheric Sciences, University of Washington, Seattle, WA 98195, USA

⁴Climate Change Research Centre, University of New South Wales, NSW 2052 Sydney, Australia

⁵School of Geosciences, University of Sydney, Sydney, Australia

⁶School of Mathematics and Statistics, University of New South Wales, NSW 2052 Sydney, Australia

*To whom correspondence should be addressed: william.Gray@lsce.ipsl.fr

The Southern Hemisphere westerly winds strongly influence deep ocean circulation and carbon storage¹. While the westerlies are hypothesised to play a key role in regulating atmospheric CO₂ over glacial-interglacial cycles²⁻⁴, past changes in their position and strength remain poorly constrained⁵⁻⁷. Here, we use a compilation of planktic foraminiferal δ¹⁸O from across the Southern Ocean and constraints from an ensemble of climate models to reconstruct changes in the westerlies over the last deglaciation. We find a 4.7° (2.9-6.9°, 95% confidence interval) equatorward shift and about a 25% weakening of the westerlies during the Last Glacial Maximum (about 20,000 years ago) relative to the mid-Holocene (about 6,000 years ago). Our reconstruction shows that the poleward shift in the westerlies over deglaciation closely mirrors the rise in atmospheric CO₂. Experiments with a 0.25° resolution ocean-sea-ice-carbon model demonstrate that shifting the westerlies equatorward substantially reduces the overturning rate of the abyssal ocean, leading to a suppression of CO₂ outgassing from the Southern Ocean. Our results establish a central role for the westerly winds in driving the deglacial CO₂ rise, and suggest natural CO₂ outgassing from the Southern Ocean is likely to increase as the westerlies shift poleward due to anthropogenic warming⁸⁻¹⁰.

The Southern Hemisphere westerly winds play a key role in returning deep ocean waters to the surface and thus largely govern the rate at which the deep oceanic reservoirs of heat and carbon communicate with the surface ocean and atmosphere¹,¹¹. South of ~47°S the modern westerly winds drive divergent northward Ekman transports in the near-surface ocean that contribute to lift deepwaters and tilt density surfaces¹¹. Although mesoscale eddies work to flatten out the steep isopycnals, the counteraction of the wind-driven circulation by the
eddies is incomplete, resulting in a residual circulation which brings macro-nutrient and carbon rich deepwaters to the surface¹,¹². Due to iron¹³ and light limitation¹⁴ the upwelled nutrients are not completely utilised before buoyancy loss close to the Antarctic continent causes some of the upwelled waters to sink as Antarctic Bottom Water, filling the deep ocean with 'preformed' nutrients¹⁵. This 'leak' in the biological pump, largely caused by the over-supply of nutrients to the surface ocean by the wind-driven upwelling, leads to the hypothesis that changes in the Southern Hemisphere westerly winds could regulate atmospheric CO₂ over glacial-interglacial cycles²⁻⁴. Modelling studies¹⁵⁻¹⁷ indicate a tight coupling between the oceans' preformed nutrient inventory and atmospheric CO₂, while proxy data⁴,¹⁸ indicate large changes in nutrient supply to the surface of the Southern Ocean over glacial-interglacial cycles.

Past changes in the position and strength of the Southern Hemisphere westerly winds are poorly constrained, making it difficult to assess their role in driving past carbon cycle changes. A recent compilation of diverse proxy data⁵ suggested signals of an equatorward shift in the westerlies during the Last Glacial Maximum (LGM, 19-23 ka) relative to the Holocene. However, relating changes in the measured proxies (i.e. terrestrial moisture, marine productivity) back to the position of the westerlies is challenging, both quantitatively and qualitatively, such that even the direction of change during the LGM (i.e. poleward versus equatorward) is debated⁶,⁷. Furthermore, while climate models show a relatively clear and consistent signal of an equatorward shift in the Northern Hemisphere near-surface westerlies under glacial forcings¹⁹⁻²¹ in good agreement with proxy data¹⁹, they show little consistency in the magnitude or sign of change in the Southern Hemisphere⁶,²⁰,²². Ice core data suggest abrupt shifts in the westerlies during the millennial scale atmospheric CO₂ variability of the last glacial period²³, but there is currently very little constraint on how or when the westerlies shifted over the last deglaciation (20 - 10 ka), as atmospheric CO₂ rose by ~90 ppmv²⁴.

To reconstruct changes in the position of the Southern Hemisphere westerly winds over the last deglaciation we exploit the coupling between the latitude of the maximum zonal wind stress (hereafter referred to as *wind latitude*) and the latitude of the maximum meridional sea surface temperature (SST) gradient (hereafter referred to as *SST front latitude*) in the Southern Ocean (Fig. 1). This two-way coupling arises through the winds' dependence on the atmospheric temperature gradient and through the winds' control on ocean circulation, which together shape the meridional SST profile²⁵⁻²⁷. The position of the peak in westerly wind-speed coincides with the position of the SST front around much of the modern Southern Ocean, except where bathymetric steering of ocean currents causes local decoupling²⁵,²⁸ (Fig 1). An ensemble of models from the Paleoclimate Modelling Intercomparison Project (PMIP3 and PMIP4) and Coupled Model Intercomparison Project (CMIP5 and CMIP6) show a tight
77

relationship between the wind latitude and the SST front latitude across pre-industrial, LGM, and 4xCO₂ simulations (Methods; Fig. 1c). This relationship, stemming from large-scale atmosphere-ocean coupling, indicates we can use changes in the SST front latitude to reconstruct past shifts in the wind latitude.

78

79

80

81

To track changes in the latitude of the Southern Ocean SST front over the deglaciation we use a basin-wide compilation of δ¹⁸O in planktic foraminiferal calcite (δ¹⁸Ocalcite). Although δ¹⁸Ocalcite is a function of both temperature and the δ¹⁸O of seawater (δ¹⁸Owater), the effect of temperature is around six times greater than the effect of δ¹⁸Owater at the basin scale (Extended Data Fig. 1). As no physical mechanism exists to drive such large changes in δ¹⁸Owater at the basin scale, the meridional pattern of δ¹⁸Ocalcite will always be dominated by temperature (Methods). Meridional profiles of δ¹⁸Ocalcite thus allow us to identify the latitude of the SST front and ultimately that of the westerly winds.

82

83

84

85

86

87

88

89

90

We compiled existing records of planktic foraminiferal δ¹⁸O from near-surface dwelling species from core sites across the Southern Ocean and generated new data from the Kerguelen plateau and southeast Pacific, resulting in a dataset of 64 planktic foraminiferal δ¹⁸O records spanning the last deglaciation (Fig. 1a; Extended Data Figs. 2 and 3; Methods). Meridional shifts in the SST front (ΔLatSST) are calculated by finding the latitudinal shift that minimises the difference between the Holocene δ¹⁸O profile and the profile at each time step within a 10° latitudinal window that includes the steepest part of the meridional δ¹⁸O profile¹⁹ (Methods; Fig. 1b and Extended Data Figs. 2 and 3). We account for whole ocean changes in δ¹⁸Owater and the global-mean SST change (δ¹⁸Oivc-gtc; Fig. 1a; Extended Data Figs. 2 and 3), quantifying uncertainties via bootstrapping and Monte-Carlo simulation (Methods). We compute ΔLatSST across the entire Southern Ocean, as well as separately in the Indian-Pacific and Atlantic sectors (Fig. 2; Extended Data Fig. 4).

91

92

93

94

95

96

97

98

99

100

101

102

103

The data reveal an equatorward shift in the SST front during the LGM (20ka) relative to 10 ka, indicative of an equatorward shift in the westerly winds (Fig. 1b; Extended Data Figs. 2 and 3). Mapping the LGM δ¹⁸Oivc-gtc anomalies shows a large mid-latitude cooling during the LGM across the Indian and Pacific sectors (Fig. 1a; equivalent to a cooling of 4-5°C beyond the global-mean SST change). Our Indian-Pacific ΔLatSST reconstruction shows a 4.8° (3.6-6.1° 95% confidence interval [CI]) equatorward shift in the SST front during the LGM relative to 10 ka (Fig. 2c; Extended Data Fig. 4). We perform a 'leave-one-out' analysis of the Indian-Pacific dataset which shows that no single record is contributing more than 5% of the total

104

105

106

107

108

109

110

111
reconstruction and indicates that Indian-Pacific ΔLat$_{SST}$ reconstruction primarily reflects a mid-latitude signal (Methods; Extended Data Fig. 4).

By contrast, we find a slight warming anomaly (relative to the global-mean LGM SST change) at all latitudes in the western Atlantic (Fig. 1), possibly a signal of a weakened Atlantic Meridional Overturning Circulation (AMOC)$^{29,30}$, and no significant change in the SST front latitude within the Atlantic sector over deglaciation (Extended Data Fig. 4). This result is consistent with strong bathymetric steering of currents in the Atlantic sector$^{28}$, and the lack of a modelled relationship between the SST front latitude and wind latitude in this sector, particularly in the eastern Atlantic where the vast majority of the mid-latitude cores are located (Fig. 1; Extended Data Fig. 5bc). Although we attribute the lack of change in the SST front latitude within the Atlantic to the modelled decoupling of the westerlies and SST front there (Extended Data Fig. 5bc), we cannot rule out that the westerlies did not shift substantially over deglaciation within the Atlantic.

We apply the multi-model relationship between the wind latitude and the SST front latitude in the Indian-Pacific sector (Fig. 2; Extended Data Fig. 5a) to our Indian-Pacific ΔLat$_{SST}$ reconstruction and quantify shifts in the zonal-mean wind latitude over the deglaciation (ΔLat$_{wind}$). A regional subset test demonstrates that our method is able to track shifts in the zonal-mean wind latitude using changes in the latitude of regional SST fronts (Fig. 2; Methods); although different regions yield differing magnitudes of deglacial change in the SST front latitude (Fig. 2b), they yield almost identical time series of the zonal-mean wind latitude (Fig. 2d), given the region-specific multi-model relationship between the two variables (Fig. 2c). Measured shifts in the SST front latitude thus provide a robust constraint on the wind latitude. We find a 6.2° (4.4-8.4° 95% CI) equatorward shift in the wind latitude during the LGM (20 ka) relative to 10 ka (Figs. 2d and 3a). The evolution of the wind latitude over deglaciation closely mirrors, and is highly correlated with, the evolution of atmospheric CO$_2$ (ref$^{24}$) (Fig. 3). A lagged correlation suggests changes in the wind latitude may have led the rise in atmospheric CO$_2$ and global temperature$^{31}$ by a few hundred years during certain intervals over the deglaciation, although difficulties in fully quantifying age uncertainties make this result tentative (Methods; Fig. 3c).

Extending our analysis to span a longer time interval (22 to 6.5 ka) indicates an early Holocene extremum in the poleward position of the westerlies, followed by a ~1.5° equatorward shift in the winds over 10-6.5 ka (Extended Data Fig. 6). Despite the larger uncertainties in the early Holocene reconstruction, it agrees well with the analysis of the full dataset (i.e., 10-20 ka) in the overlapping sections (Extended Data Fig. 6). Our results thus
indicate a 4.7° (2.9-6.9° 95% CI) equatorward shift of the westerlies during the LGM (20 ka) relative to the mid Holocene (6.5 ka). Assuming no substantial shift in the westerlies from the mid- to late-Holocene (Methods), the magnitude of the reconstructed wind shift is substantially greater than that predicted by any of the models within the PMIP3/4 ensemble between preindustrial and LGM states (Extended Data Fig. 7). The models show a relationship between the wind latitude and the maximum magnitude of the zonal-mean zonal wind stress (wind strength; Extended Data Fig. 5e), such that the reconstructed shift in wind latitude also implies a reduction in wind strength of ~25% during the LGM relative to the mid Holocene (Methods; Extended Data Fig. 8).

The similarity of the changes in wind latitude and atmospheric CO₂ over the deglaciation (Fig. 3) reinforces the hypothesis of their coupling through Southern Ocean circulation and carbon cycling. While modelling studies typically show a consistent increase in oceanic carbon storage following a weakening of the westerlies, the impact of shifts in the latitude of the westerlies is more uncertain¹⁷,³². To better understand how changes in the latitude of the westerlies may affect the oceanic overturning circulation and carbon cycle, we performed two experiments with a global ocean-sea-ice-carbon model with 0.25° horizontal resolution³³,³⁴: A Control experiment is forced by climatological atmospheric forcing representative of the recent instrumental period. A Perturbed experiment uses the same forcing except for a uniform 4° equatorward shift of the Southern Hemisphere westerlies, with no change in their magnitude (Methods). Both experiments are run for 125 years, starting from the near-equilibrium control state. As the applied wind stress forcing does not include the implied 25% reduction in wind strength (and is smaller than our reconstructed LGM shift) this simulation is not designed to represent our 'best estimate' of LGM Southern Hemisphere wind stress, but rather represents a conservative estimate which enables us to assess the broad impacts of an equatorward shift in the westerlies alone. The 125-year transient response of the system forced by an equatorward shift in the westerlies does not allow quantification of the equilibrium response of the deep ocean nutrient and carbon cycles¹⁷. It nevertheless reveals clear trends in circulation and biogeochemistry which provide an indication of how the rapidly responding Ekman-driven transport reorganizes the residual circulation, allowing us to assess how an equatorward shift in the westerlies may qualitatively impact the carbon cycle on longer timescales.

We find a complete suppression of CO₂ outgassing south of 60°S in the Perturbed experiment (Fig. 4), with only a partial compensation further north. As such, there is an anomalous uptake of 27 GtC by the Southern Ocean south of 35°S over the course of the Perturbed experiment, equivalent to a CO₂ decrease of 13 ppm (Fig. 4). Deepwater upwelling
and surface nutrient and carbon concentrations are substantially reduced south of 60°S (Fig. 4), indicating that reduced exposure of nutrient and carbon rich deepwaters in the polar Southern Ocean underpins the simulated carbon cycle response to equatorward-shifted westerlies (Methods). As the winds shift equatorward relative to their modern position, northward Ekman transports become more divergent north of about 60°S, but less divergent south of 60°S (Fig. 4). In our simulation, the effect is a substantial decrease in upwelling within the polar Southern Ocean south of 60°S (Fig. 4) and a pronounced slowdown in the residual circulation of the abyssal ocean globally (Fig. 5), with the largest reduction within the Pacific (Extended Data Fig. 9). The decrease in abyssal overturning results in increased storage of carbon and regenerated nutrients throughout the deep ocean below ~1.5 km, concurrent with a decrease in dissolved oxygen (Fig. 5c; Extended Data Fig. 10). Conversely, we see an increase in upwelling north of 60°S and increased overturning at intermediate depths, such that carbon concentrations decrease in the upper ~1.5 km (Fig. 5c; Extended Data Fig. 10). Hence, although shifting the winds equatorward increases the overall Ekman divergence across the Southern Ocean (Fig. 4b), it focuses the wind's energy away from isopycnals outcropping carbon-rich abyssal waters, toward isopycnals outcropping intermediate waters containing relatively less carbon (Fig. 5b); these changes are synergistic such that the concentration of carbon and nutrients within the intermediate depths decreases as the overturning at these depths spins up. The net effect in the model is an increase in oceanic carbon storage (Fig 5c). In contrast, some previous studies using coarser-resolution models simulated a decrease in oceanic carbon content in response to an equatorward shift in the westerlies¹⁷,³². We attribute this difference to the response of the residual circulation in the Southern Ocean, the representation of which should improve at higher resolution³⁵,³⁶.

Crucially, the reduction in the upwelling of nutrients to the surface of the polar Southern Ocean lowers the preformed nutrient concentration of Antarctic Bottom Water (Fig. 4; Extended Data Fig. 10). Preformed nutrient concentrations also decrease in the upper cell (Methods), such that preformed nutrient concentrations begin to decrease throughout the deep ocean (Extended Data Fig. 10). The simulated circulation changes with equatorward-shifted winds thus imply a substantial long-term increase in the efficiency of the biological pump and decrease in atmospheric CO₂ (Methods; refs¹⁵,¹⁶). In addition to the overall increase in regenerated nutrients and carbon, the model indicates a redistribution of the regenerated nutrient and carbon pools toward the abyss (Fig. 5c; Extended Data Fig. 10), such that the vertical carbon gradient increases. This implies a further reduction atmospheric CO₂ via carbonate compensation¹⁶,³⁷ (Methods).
Despite the 4° equatorward shift in the Southern Hemisphere westerlies being the only perturbation applied, and the short duration of the Perturbed experiment, the sign of the simulated circulation and carbon cycle changes concurs with proxy observations from the LGM of a more sluggish abyssal circulation³⁸, an increase in regenerated nutrients and carbon within the deep ocean³⁹⁻⁴¹ and a redistribution of regenerated nutrients and carbon toward the abyss³⁹,⁴⁰,⁴², as well as a shoaling of the AMOC²⁹ (Fig. 5; Extended Data Fig. 9). The simulated decrease in nutrient upwelling and export production within the polar Southern Ocean, and the increase further north (Fig. 4; Extended Data Fig. 10), are also in good agreement with LGM proxy data³,⁵,⁴³.

The sign of the simulated trends is unlikely to be contingent on the absence of a simulated change in atmospheric temperature or on model approximations such as parameterized eddy effects (Methods). Indeed, the trends directly ensue from changes in Ekman transports—themselves a function only of wind stress and geometry—and their position relative to the underlying bathymetry. The northward shift of the Ekman divergence (Fig. 4b,c) explains the northward shift in upwelling and CO₂ outgassing (Fig. 4d,e; Methods). This shift does not merely draw the overturning northward but shoals the upper cell and weakens the abyssal cell (Fig. 5b; Extended Data Fig. 9). This response stems from the influence of bathymetry on the vertical extent of wind-driven upwelling¹¹,⁴⁴. Specifically, the presence of a zonally continuous channel at Drake Passage latitudes (56-60°S) above 2 km depth favors deeper waters as the mass replacement for the surface divergence to its south¹¹,⁴⁵. Displacement of the Ekman divergence from south of 60°S to lower latitudes thus suppresses the privileged upwelling pathway of abyssal waters (> 2 km depth) and instead favors overturning in the upper ocean. Although the representation of surface buoyancy forcing and mesoscale eddies can affect the response of the overturning to changing winds¹²,⁴⁴, the simulated slowdown of the abyssal overturning is a robust consequence of its bathymetry-driven sensitivity to Ekman flow in the polar Southern Ocean.

Our results thus suggest that as the winds shifted poleward through the deglaciation, their ability to lift abyssal waters to the surface increased in approximate proportion to the northward Ekman transport at 60°S. Using the relationship between wind latitude and the zonal-mean wind stress at 60°S within the PMIP/CMIP ensemble (Extended Data Fig. 5d; Extended Data Fig. 8) we estimate that over deglaciation this transport increased from an LGM minimum of around $4 \times 10^6 \text{ m}^3 \text{ s}^{-1}$ to a maximum of $18 \times 10^6 \text{ m}^3 \text{ s}^{-1}$ at 11 ka, before decreasing back to $14 \times 10^6 \text{ m}^3 \text{ s}^{-1}$ by the mid-Holocene (Fig. 6a; Methods). This invigoration of the abyssal overturning over deglaciation would have driven carbon out of the abyssal ocean and back to the intermediate depths, surface, and atmosphere. This concurs with
records of deep Pacific oxygen⁴¹ and deep Southern Ocean pH³⁹ which suggest of a loss of regenerated nutrients and carbon from the abyssal ocean over deglaciation (Fig. 6c). Furthermore, a decrease in the global-mean deep-intermediate carbon isotope gradient⁴², and other carbon cycle tracers³⁹,⁴⁶, suggest the vertical gradient of regenerated nutrients and carbon weakened over deglaciation (Fig. 6b). Finally, records of nutrient utilisation¹⁸ and export production⁴ support an increase in nutrient upwelling to the Southern Ocean’s surface over deglaciation, with boron isotope records⁴⁷ demonstrating a concurrent increase in CO₂ outgassing from the Southern Ocean. The early Holocene maximum in poleward position of the winds may have led to an overshoot of the oceanic carbon cycle, with enhanced upwelling and CO₂ outgassing persisting until ~6 ka⁴,⁴⁷. However, atmospheric CO₂ decreased by ~10 ppm from 10 to 6 ka suggesting uptake by another sink, likely the terrestrial biosphere⁴⁸.

Our findings indicate that shifts in the Southern Hemisphere westerlies are likely to have played a key role in driving the deglacial rise in atmospheric CO₂ (ref³³), and thus may be an important mechanism underlying glacial-interglacial CO₂ variations²,³. Given that atmospheric CO₂ and global temperature may also influence the latitude of the westerly winds⁸⁻¹⁰, a deglacial feedback mechanism has been proposed². The apparent temporal lead of shifts in the westerlies over atmospheric CO₂ and global temperature during certain intervals of the deglaciation indicated by our reconstruction (Fig 3c) suggests that some initial change in the winds, perhaps driven by obliquity¹⁸,⁴⁹, could have initiated a cascade of increasing CO₂, global warming, and poleward shifting winds. The tight coupling demonstrated here between the latitude of the westerlies and atmospheric CO₂ over the last deglaciation, together with the sensitivity of the abyssal circulation and carbon cycle to the latitude of the westerlies in our ocean-sea-ice-carbon model experiments, suggest that future poleward shifts in the position of the westerly winds⁸⁻¹⁰ are likely to drive a positive feedback on anthropogenic warming—via a decrease in the efficiency of the biological pump and an increase in natural CO₂ outgassing from the Southern Ocean⁵⁰.

## Acknowledgements

Some of the collaborations in this study were initiated through the PAGES QUIGS workshop 2020. EM received financial support from French-Swedish project on SO VR-349-2012-6278 and CNRS-INSU project INDIEN-SUD. RCJW was funded by NSF Grant AGS-1929775. LM acknowledges funding from Australian Research Council grant FT180100606. PS was supported by the Australian Research Council grant FT190100413 and the Australian National Computing Infrastructure. MH acknowledges funding from Australian Research Council grant DP210101650. We acknowledge the World Climate Research Programme's Working Group
295

on Coupled Modelling for the coordination of CMIP and thank the climate modelling groups
for producing and making available their model output.

296

297

298

**Author contributions**

299

WRG, EM, LM and CdL designed the study. EM generated the new δ¹⁸O data. WRG compiled the δ¹⁸O data and performed the statistical analysis. RCJW and WRG analysed the coupled climate model output with input from MK. LM and PS conducted the ocean-carbon model simulations and analysed the output with CdL. MH provided an analysis of deep-to-surface transport in the modern ocean to aid interpretation. WRG, CdL and LM developed the interpretation with input from all authors. WRG led the writing of the manuscript. All authors contributed to the preparation of the manuscript.

300

301

302

303

304

305

306

307

**Competing interests**

308

The authors declare no competing interests.

309

310

**Correspondence and requests for material** should be addressed to WRG

311

312
<figure><img src="image_2.png" /><figcaption>**Fig. 1:** δ<sup>18</sup>O data and relationship between wind latitude and SST front latitude in the model ensemble. (a) LGM-Holocene δ<sup>18</sup>O<sub>ivc-gtc</sub> (δ<sup>18</sup>O<sub>calcite</sub> corrected for ice volume and global-mean SST changes; Methods) at the core sites with climatological meridional δ<sup>18</sup>O gradient represented by background shading (Extended Data Fig. 1; darkest shade represents 0.25 ‰/°Lat, equivalent to ~1 °C/°Lat). Modern position of maximum westerlies<sup>51</sup> (solid red line) and the estimated LGM position of the westerlies based on a uniform shift (dashed red line). (b) LGM and Holocene δ<sup>18</sup>O<sub>ivc-gtc</sub> meridional profiles. The grey box shows latitudinal window in which ΔLat<sub>SST</sub> is calculated (Methods). Symbols on (a) and (b) distinguish foraminiferal species. (c) Relationship between positions of the wind latitude (latitude of maximum τ<sub>u</sub>) and SST front latitude (latitude of maximum ∂SST/∂Lat) within the PMIP3/CMIP5 and PMIP4/CMIP6 ensemble (Methods).</figcaption></figure>

313

314

315

316

317

318

319

320

321

322

323

324

325
<figure><img src="image_3.png" /></figure>

**Fig. 2:** SST front and westerly wind changes based on regional subsets (a) Map showing regional subsets of data. Light grey dashed line corresponds to all core sites. Black includes all Indian-Pacific sites, blue has eastern Pacific sites removed, while pink has western Indian sites removed. Note that given the paucity of data from south of 65 °S we include Antarctic marginal sites from all sectors in all regional subsets. (b) Relationship between wind latitude (latitude of maximum zonal-mean $τ_u$) and regional SST front latitude (latitude of maximum $∂SST/∂Lat$) within the model ensemble. (c) Reconstructed change in SST front latitude ($ΔLat_{SST}$) within the regional subsets. (d) Reconstructed zonal-mean wind latitude ($ΔLat_{wind}$) calculated from the regional $ΔLat_{SST}$ reconstruction (c) and the relationships between regional SST front and zonal-mean wind latitude given in (b).
<figure><img src="image_4.png" /><figcaption>**Fig. 3: Deglacial shifts in the westerlies and atmospheric CO₂.** (a) Deglacial change in the position of the wind latitude (ΔLat<sub>wind</sub>, lines show the 5th, 32nd, 50th, 68th, and 95th percentiles) and atmospheric CO₂ (ref<sup>24</sup>) over time. HS1, ACR, and YD are Heinrich Stadial 1, Antarctic Cold Reversal, and Younger-Dryas, respectively. (b) Correlation between ΔLat<sub>wind</sub> and atmospheric CO₂ over deglaciation. (c) Lead-lag between changes in ΔLat<sub>wind</sub> and changes in atmospheric CO₂ and global temperature<sup>31</sup> over deglaciation (Methods).</figcaption></figure>

337

338

339

340

341

342

343

344
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4:** Modelled impact of shifted westerlies on Southern Ocean CO₂ outgassing. Zonally averaged (a) zonal wind stress, (b) northward Ekman transport, (c) Ekman divergence and, (d) ocean-atmosphere CO₂ flux (positive flux for ocean outgassing). Black curves correspond to the Control equilibrium state (CTR) and red curves to the average over years 116-126 of the Perturbed experiment (EQ). The grey box on (a-d) indicates the latitude of the Drake Passage (56-60°S). 1 Sv is 10⁶ m³ s⁻¹. Map of (e) ocean-atmosphere CO₂ flux (mol/m²/yr, positive flux for ocean outgassing) and (f) nitrate (mmol/m³, average over upper 149 m depth) anomalies, obtained as difference between Perturbed and Control simulations.</figcaption></figure>

345

346

347

348

349

350

351

352

353
<figure><img src="image_7.png" /><figcaption>**Fig. 5: Modelled impact of shifted westerlies on deep ocean circulation and carbon storage. (a)** Northward Ekman transport in the Control (Black) and Perturbed (red) simulations. 1 Sv is $10^6 \text{ m}^3 \text{ s}^{-1}$. (b) Global MOC anomaly (Perturbed-Control) after 125 years, normalised for sign such that blue is a decrease in overturning rate and red is an increase. Contours show Control (solid) and Perturbed (dashed) isopycnals (c) Dissolved Inorganic Carbon (DIC) and Oxygen anomalies (Perturbed-Control) horizontally integrated globally. Zonally averaged concentration anomalies are shown in Extended Data Fig. 10). Grey box on (a) and (b) indicates the Drake Passage.</figcaption></figure>

354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390
<figure><img src="image_8.png" /><figcaption>**Fig. 6: Deglacial changes in northward Ekman transport at 60°S and deep ocean carbon cycling.**

**(a)** Reconstructed change in northward Ekman transport at 60°S over deglaciation (lines show the 5th, 32nd, 50th, 68th, and 95th percentiles). **(b)** Global volume-weighted deep-intermediate $\delta^{13}$C gradient$^{42}$, broadly indicative of the vertical gradient in regenerated nutrients/carbon. **(c)** $O_2$ in the deep Pacific$^{41}$ (fit with a LOESS smooth), indicative of regenerated carbon storage.</figcaption></figure>
## Materials and methods

### Planktic foraminiferal $\delta^{18}$O records from the Kerguelen plateau and southeast Pacific

We generated new planktic foraminiferal $\delta^{18}$O data from sediment cores spanning the last deglaciation located in the mid-latitudes of the Indian and Pacific sectors of the Southern Ocean. We generated new records from two sediment cores retrieved from the Kerguelen Plateau during cruise OSCAR INDIEN-SUD (MD12-3396CQ, -47.73 °N, 87.69 °E; MD12-3401CQ, -44.68°S, 80.39°E). Furthermore, we extended/increased the resolution of two previously published records from cores located on Kerguelen (MD02-2488, -51.07 °N, 67.73 °E) and in the southeast Pacific (MD07-3119, -46.08 °N, -76.1 °E). We analysed $\delta^{18}$O on either *G. bulloides* or *N. pachyderma* using a GV Isoprime 100 and an OPTIMA, and a Finnigan MAT251 and a Δ+ at CFR/LSCE laboratory. The measurements are reported versus Vienna Pee Dee Belemnite standard (VPDB) defined with respect to the NBS19 standard. The mean external reproducibility ($1\sigma$) of carbonate standards is ±0.06% for $\delta^{18}$O; the different mass spectrometers are regularly inter-calibrated and the data are corrected, depending on the devices, for nonlinearity and the common acid bath. Within this internal calibration, NBS18 is -23.2±0.2% VPDB for $\delta^{18}$O and -5.0±0.1% VPDB for $\delta^{13}$C. Age models for all the cores are based on radiocarbon dating, and further details of the age models can be found in ref⁵² for core MD07-3119 and ref⁵³ for core MD12-3396CQ. Reservoir age changes for the Kerguelen area followed recent results by ref⁵³ to establish the age model of core MD12-3401CQ. The new data are provided in Table S1 and are available on Pangaea (DOI pending).

### Southern Ocean planktic foraminiferal $\delta^{18}$O compilation

We compiled all available existing $\delta^{18}$O records for near-surface dwelling planktic foraminifera species (*G. ruber*, *G. bulloides*, *N. pachyderma*) spanning the last deglaciation (10-20 ka) from across the Southern Ocean (refs⁵⁴⁻⁹⁷). All records are kept on the original age model of publication. The compilation contains 64 deglacial records of planktic foraminiferal $\delta^{18}$O. All $\delta^{18}$O data are given in table Table S1 and are available on Pangaea (DOI pending).

### Using planktic foraminiferal $\delta^{18}$O to track the latitude the SST front

Although $\delta^{18}$O$_{calcite}$ is a function of both temperature and $\delta^{18}$O$_{water}$, at the basin scale the effect of temperature dominates over $\delta^{18}$O$_{water}$ (Extended Data Fig. 1). Using the Southern Ocean salinity-$\delta^{18}$O$_{water}$ relationship of ref⁹⁸ a meridional salinity difference of greater than 25 PSU would be required to equal the meridional temperature changes across the basin. As no physical mechanism exists to drive such changes, the meridional pattern of $\delta^{18}$O$_{calcite}$ will always be dominated by temperature, enabling us to use meridional profiles of $\delta^{18}$O$_{calcite}$ to identify the latitude of the SST front¹⁹.

We quantify changes in the position of the SST front through time ($\Delta Lat_{SST}$) from following the method of ref¹⁹. Briefly, we first interpolate the $\delta^{18}$O data to 250-year time steps from 20 ka to 10 ka using a Generalised Additive Model (GAM)⁹⁹, with the smoothing term determined by restricted maximum likelihood (REML)¹⁰⁰. The reader is referred to ref¹⁰¹ for an overview of GAMs. Only foraminiferal $\delta^{18}$O records that span the entire time period of the reconstruction are utilised such that our analysis is always comparing relative changes in the same cores through time. The mean resolution of the individual records over deglaciation is about 1 point
per 250 yrs and we only include δ¹⁸Oₗcalcite records with a minimum of 1 point per 2 ka over the deglaciation.

We model the δ¹⁸O data at each time step (first correcting for whole ocean effects, see below) as a function of latitude using a GAM⁹⁹, with the smoothing term determined by restricted maximum likelihood (REML)¹⁰⁰ (Extended Data Figs. 2 and 3). We compute the shift in latitude which minimises the Euclidean distance (L²) between the GAM fit at each time step relative to 10 ka, within a 10° latitude band centred around the steepest part of the Holocene meridional SST/δ¹⁸Oₗcalcite profile (-50 to -40 °N; grey box in Fig. 1b and Extended Data Figs. 2 and 3). The 10 ka reference time is chosen to maximise the number of records spanning the deglaciation. Our analysis thus tracks changes in the position of the steepest part of meridional SST profile across a 10° latitudinal window, rather than the position of any individual front.

To minimise changes in the meridional δ¹⁸Oₗcalcite profile between different time steps that arise from whole-ocean changes rather than local dynamics we correct the δ¹⁸Oₗcalcite data for the whole ocean change in δ¹⁸Oₗwater (arising from ice sheet growth/retreat) and the global-mean SST change (δ¹⁸Oᵢₗ-gtct). For the whole-ocean change in δ¹⁸Oₗwater we scale the LGM-Holocene change of 1±0.1 ‰ (2σ) (ref¹⁰²) to the sea level curve of ref¹⁰³. For the global-mean SST change we scale the -1.7±1.0 °C (2σ) area-weighted global-mean LGM-preindustrial change in SST from the PMIP3/4 ensemble (see below) to the global temperature record of ref³¹, using the water-calcite temperature fractionation (δ¹⁸Oₗcalcite-water) of ref¹⁰⁴, which agrees well with the pooled-foraminiferal species sensitivity of ref¹⁰⁵. We propagate all of the uncertainty from these whole-ocean corrections through to our final results (see below). We note these whole-ocean corrections make the reconstructed changes in the position of the winds more conservative; removing the global-mean SST change correction entirely results in an Indian-Pacific ΔLatₙSST of ~7° between 10-20 ka (c.f. ~5° including the correction). Furthermore, our leave-one-out analysis shows that our reconstruction is primarily driven by mid-latitude sites, indicating these whole ocean corrections are unlikely to be biasing our results.

We derive uncertainties via bootstrapping¹⁰⁶ (10,000 iterations), accounting for the age and analytical uncertainties on individual records as well as the uncertainties in the whole-ocean change in δ¹⁸Oₗwater and the global-mean SST change with Monte-Carlo simulation. We ascribe conservative age uncertainties of ±1000 years (2σ), and analytical uncertainties of ±0.08 ‰ (2σ) to all planktic foraminiferal δ¹⁸O records. The R code and data used to perform the analysis is provided on Github (https://github.com/willyrgray/SOd18O).

We repeat the analysis over a longer time period (22-6.5 ka) which necessitates using a smaller subsection of cores and thus results in larger uncertainties (Extended Data Fig. 6). Still, the results show excellent agreement with the complete set of cores (i.e. spanning 10-20 ka) (Extended Data Fig. 6), with a ratio of 0.95 during the overlapping interval. We apply this correction to the longer reconstruction, though the correction makes a negligible difference and does not impact our conclusions (Extended Data Fig. 6). Performing the analysis further into the Holocene (22-2.5 ka) results in substantially higher uncertainties due to the lack of cores spanning this interval, but the results suggest the westerlies are relatively stable through the late Holocene (Extended Data Fig. 6), in agreement with qualitative indicators of the winds²³,¹⁰⁷.
Leave-one-out analysis

We perform a jackknife resampling of the dataset to determine the contribution of each planktic foraminiferal $\delta^{18}$O record to the Indian-Pacific $\Delta\text{Lat}_{\text{SST}}$ reconstruction. We sequentially remove each record from the dataset, recompute $\Delta\text{Lat}_{\text{SST}}$, and determine the contribution of that record as the time-integrated absolute difference from $\Delta\text{Lat}_{\text{SST}}$ computed using the entire dataset (expressed as a percentage of the absolute cumulative change in $\Delta\text{Lat}_{\text{SST}}$ over deglaciation; Extended Data Fig. 4). This analysis shows that no single record is contributing more than 5% of the total signal, and that the cores with the highest weighting are all located in the mid-latitudes. Thus, we conclude that our reconstruction is primarily tracking a mid-latitude signal.

Seasonality of planktic foraminifera

Our approach assumes that any changes in the seasonal bias of foraminifera relating to their habitat preference are small relative to the change in temperature due to the movement of the SST front. The validity of this approach is supported by $\delta^{18}$O records measured on species of foraminifera with different temperature/seasonal habitats at the same (or nearby) mid-latitude sites within the compilation. Foraminiferal species with different habitat temperature (and thus seasonal) preferences (*G. bulloides* and *N. pachyderma*, ref $^{108}$) show very similar Holocene-LGM changes (compare circles and triangles on Fig. 1a). Furthermore, the leave-one-out analysis (Extended Data Fig. 4) shows that records of *G. bulloides* and *N. pachyderma* $\delta^{18}$O (which should have different seasonal biases) are both contributing highly to the $\Delta\text{Lat}_{\text{SST}}$ reconstruction. This suggests the impacts of seasonality on the reconstruction are likely to be small. Crucially, the relatively warm temperature preference of *G. bulloides* (the dominant species in the mid-latitudes) means that global SST cooling during the LGM would shift the annual SST distribution further away from its preferred habitat temperature, shifting its seasonal bias further towards the summer and minimising the degree of glacial cooling it records. This would thus make the anomalous mid-latitude cooling we see in the LGM more conservative relative to the annual mean cooling, thus making our reconstruction of changes in the wind latitude more conservative.

PMIP3/4 and CMIP5/6 ensemble

We use an ensemble of models from the Paleoclimate Modelling Intercomparison Project (PMIP3$^{109}$ and PMIP4$^{20,110}$) and Coupled Model Intercomparison Project (CMIP5$^{111}$ and CMIP6$^{112}$); data available at https://esgf-node.llnl.gov/projects/esgf-llnl/. We use all CMIP5/PMIP3 and CMIP6/PMIP4 models for which both SST and zonal surface wind fields are available for the LGM and preindustrial simulations, and also include the 4xCO$_2$ simulations from these models where available. Annual-mean climatologies are calculated from the final 100 years of each simulation and interpolated to a common 2-degree analysis grid.

To convert the reconstructed changes in the SST front latitude ($\Delta\text{Lat}_{\text{SST}}$) to changes in the latitude of the westerly winds we use the relationship between the SST front latitude and the wind latitude within the ensemble. We calculate the SST front latitude as the latitude of maximum meridional SST gradient ($\partial\text{SST}/\partial\text{Lat}$) in each simulation within a 10° latitudinal averaging window (to emulate the 10° latitudinal window used to calculate $\Delta\text{Lat}_{\text{SST}}$ from the $\delta^{18}\text{O}_{\text{ivc-gtc}}$ data downcore. We do this using zonal-mean SST across the Southern Ocean (Fig. 1c) and across regional subsets (Fig. 2; Extended Data Fig. 5). We calculate the wind latitude as the latitude of maximum zonal-mean zonal wind stress ($\tau_u$). We regress the SST front
latitude against wind latitude (Fig. 2; Extended Data Fig. 5) and then apply this relationship to our $\Delta\text{Lat}_{\text{SST}}$ reconstruction in order to determine changes in the position of the wind latitude ($\Delta\text{Lat}_{\text{wind}}$), propagating the uncertainty in the relationship through to our final estimates of $\Delta\text{Lat}_{\text{wind}}$ using Monte-Carlo simulation. We exclude the MPI model from the regressions for the Indian-Pacific sector, as this model sits as an outlier from the remainder of the ensemble (however we note the slope between wind latitude and SST front latitude within the MPI model is consistent with the rest of the ensemble). The peculiarity of the MPI model is possibly due to SST biases relating to meridional heat transport in the Indian sector within the model$^{113}$, where we find the largest difference relative to the other models. Including the MPI model in the Indian-Pacific sector regression has a negligible effect on our results, slightly increasing the reconstructed change in the wind latitude ($\Delta\text{Lat}_{\text{wind}}$) between 10-20 ka from 6.2° to 6.7° and increasing the uncertainty by 0.5° at the 95% CI. The $\Delta\text{Lat}_{\text{wind}}$ reconstruction is provided in Table S2.

### Regional subset test

To ascertain that our method is able to track zonal-mean shifts in the wind latitude from changes in the SST front latitude (which is also affected by bathymetry and ocean eddies$^{25,28,114}$), we perform a regional subset test. We sequentially remove regional subsets of cores from the compilation, re-compute $\Delta\text{Lat}_{\text{SST}}$, and recalculate the relationship between the zonal-mean wind latitude and the regional SST front latitude (Fig. 2). The results demonstrate that the difference in zonal-mean $\Delta\text{Lat}_{\text{wind}}$ reconstructed from the regional subsets is markedly smaller than the difference in $\Delta\text{Lat}_{\text{SST}}$ between the regional subsets. This indicates that (i) the model ensemble and (ii) our method to track the SST front downcore, capture the regional SST front response to zonal-mean shifts in the wind latitude. We can thus confidently apply the relationship observed in the model ensemble to estimate shifts in the zonal-mean latitude of the westerlies using the Indian-Pacific $\Delta\text{Lat}_{\text{SST}}$ reconstruction over deglaciation.

### Reconstructing peak wind stress

Using our $\Delta\text{Lat}_{\text{wind}}$ reconstruction and the relationship between wind latitude and the magnitude of peak wind stress in the model ensemble (wind strength, Extended Data Fig. 5), we estimate changes in the magnitude of peak wind stress over deglaciation (Extended Data Fig. 8). The tendency for the magnitude of peak wind stress to decrease as the winds move equatorward (and vice-versa) has been previously noted$^{115,116}$. Our reconstructed equatorial shift in the wind latitude implies a weakening of the peak westerlies by 0.034 N m$^{-2}$ (about 25%) during the LGM relative to the mid Holocene, resulting in a LGM wind strength of 0.106 (0.085-0.12, 95% CI) N m$^{-2}$, assuming mid-Holocene wind stress is equal to the modern climatology (0.14 N m$^{-2}$; Fig. 4a). This assumption is supported by qualitative tracers of the westerlies$^{23,107}$ which indicate little change between ~6.5 ka and the present day. Furthermore, running our analysis further into the Holocene suggests little change in the position of the winds (Extended Data Fig. 6), although uncertainties are very large given the few core sites that span this interval. Although the shape of the zonal-mean westerly wind profile may have changed together with the latitude and magnitude of its peak, the model ensemble suggests that the meridional wind profile changes little between PI and LGM states; there is a -0.25±7% and 0.5±3.5% ensemble mean change in peak width at 50% and 15% peak height, respectively, between LGM and preindustrial. Furthermore, we observe no significant relationship between the wind latitude and the peak width (at either 50% or 15% peak height) within the ensemble, such that we do not expect substantial changes in the
shape of the wind profile as the winds shift. The wind strength (max τ_u) reconstruction is
provided in Table S2.

Reconstructing northward Ekman transport at 60°S

To calculate changes in zonal wind stress at 60°S (60S τ_u) we use the relationship between the wind latitude and the zonal-mean zonal wind stress at 60°S in the model ensemble (Extended Data Fig. 5d), and apply this relationship to our reconstruction of ΔLat<sub>wind</sub>, propagating through the uncertainty in the relationship (Extended Data Fig. 8). We assume mid-Holocene wind stress equal to the modern climatology (0.09 N m<sup>-2</sup>; Fig. 4a). We calculate northward Ekman transport (Fig. 6a) as the zonal integral of τ_u/(ρ<sub>0</sub>*f), where ρ<sub>0</sub> is the density of seawater (1027 kg/m<sup>3</sup>) and f is the Coriolis parameter. The reconstruction of 60S τ_u and northward Ekman transport at 60°S is provided in Table S2. The ensemble shows that the latitude of the peak in zonal wind stress is a better predictor of wind stress at 60°S (R<sup>2</sup>= 0.9; Extended Data Fig. 5d) than the magnitude of the peak in wind stress (R<sup>2</sup>= 0.7).

Lead-lag analysis

We calculate the lead-lag between the Indian-Pacific ΔLat<sub>SST</sub> reconstruction and the change in atmospheric CO<sub>2</sub> (or in global temperature) over deglaciation as the time offset that maximises the correlation between the two time series, broadly following the approach of ref<sup>31</sup>. We account for age uncertainties in the CO<sub>2</sub> record using the typical gas age uncertainty in the WAIS divide ice core over deglaciation<sup>117</sup> (±75 yrs, 1σ). For the age uncertainty in the global mean temperature stack we take the uncertainty in the lag of global temperature over CO<sub>2</sub> (±340 yrs, 1σ) in ref<sup>31</sup>. To account for age uncertainties in the Indian-Pacific ΔLat<sub>SST</sub> reconstruction we repeat the lead-lag analysis using each boot strap/Monte-Carlo iteration of the Indian-Pacific ΔLat<sub>SST</sub> reconstruction. This approach assumes age uncertainties within the δ<sup>18</sup>O compilation are uncorrelated, which is unlikely to strictly hold given e.g. reservoir age changes (although methodological differences between studies add a source of random 'human behavioural' error), but provides a first-order assessment of leads and lags between the different time series. As a sensitivity test, we repeat the lead-lag analysis with CO<sub>2</sub>, including varying degrees of 'structural' uncertainty in the Indian-Pacific ΔLat<sub>SST</sub> reconstruction. This suggests that the lead in the change in the winds over CO<sub>2</sub> is significant at the 95% level until more than ~30% of the age uncertainty is correlated across the compilation. Assuming 100% of the age uncertainty is correlated (i.e. perfect covariance) results in uncertainties of ±860 yrs (95% CI). To test whether the lead of ΔLat<sub>SST</sub> relative to CO<sub>2</sub> holds in the early deglaciation, we repeat this analysis for 20-14 ka (cf. 20-10 ka) and find a lead of 160 yrs (-10 to 330 yrs 95% CI), cf. 330 yrs (100 to 560 yrs 95% CI).

Experiments with a 0.25° ocean-sea-ice-carbon model

We use the ocean-sea-ice-carbon model MOM5-SIS-Wombat with a Mercator horizontal resolution of 0.25° (~11 km grid spacing at 65°S), and 50 vertical levels<sup>33,34</sup>. The model is initialised with modern-day temperature and salinity distributions, and biophysical fields derived from an observation-based climatology (GLODAP v2, 2016)<sup>118</sup>. The model is then spun-up for 700 years with version 2 of the Coordinated Ocean-ice Reference Experiments (CORE) Normal Year Forcing (NYF) reanalysis data<sup>119</sup>, representative of a 'normal year' during the recent instrumental period. To study the impact of equatorward shifted southern hemispheric westerlies, a 4° equatorward shift (with no change in magnitude) is applied to the near surface wind speeds between 25°S and 70°S. The perturbation simulation is run for 125
years. Given the carbon and nutrient cycles have around an order of magnitude longer time scales in the deep Pacific¹⁷, the perturbation simulation is far from equilibrium, but it does provide a strong indication of the impact on circulation and the carbon cycle. Anomalies are calculated as the difference, averaged over the last 10 years, between the Perturbed experiment and the concomitantly extended Control experiment to remove the influence of linear model drift.

The model includes a parameterization for the effects of mesoscale eddies, with isopycnal tracer¹²⁰ and thickness¹²¹ diffusivities set to 600 m² s⁻¹. The meridional overturning circulation (Fig. 5b and Extended Data Fig. 9) includes both resolved and parameterized advection. It is calculated in neutral density¹²² coordinate and reprojected onto the depth coordinate¹²³ in order to eliminate adiabatic recirculations and avoid spurious effects due to vertical inversions in potential density fields.

As a result of the equatorward shifted winds, the Ekman pumping is suppressed south of 60°S, and the ACC weakens by 6% (122 Sv instead of 130 Sv). South of 60°S, the surface DIC concentration decreases by 13 mmol/m³ (Extended Data Fig. 10), thus leading to a surface pCO₂DIC decrease of 25 ppm, which is partly compensated by a decrease of surface alkalinity of 8.6 mmol/m³ (Delta pCO₂alk = +14 ppm). A small decrease in surface salinity is compensated by a small increase in SST, so that the solubility contribution to pCO₂sol is negligible south of 60°S. As a result, the CO₂ outgassing over the polar Southern Ocean (south of 60°S) is completely suppressed (Fig. 4d). On the other hand, the equatorward shift in the westerlies and associated increased Ekman pumping north of 60°S leads to a greater CO₂ outgassing in that region. Nevertheless, if integrated over the Southern Ocean (south of 35°S), there is an anomalous uptake of 27 GtC by the Southern Ocean over the course of the Perturbed experiment, equivalent to a CO₂ decrease of 13 ppm. While the model also displays an increase in Southern Ocean sea ice extent in the Perturbed experiment (not shown), likely an impact of the reduction in upwelling of relatively warm deepwaters, the decrease in CO₂ flux out of the polar ocean is associated with decreased surface ocean pCO₂ rather than increased air-sea disequilibria. Hence, changes in wind-driven carbon supply, rather than sea ice driven disequilibria, cause the reduction in CO₂ outgassing from the polar Southern Ocean. Furthermore, the CO₂ flux anomalies are largely decoupled from the sea ice anomalies around the basin. Both pCO₂ and CO₂ flux anomalies are, however, clearly linked to changes in Ekman divergence.

As a result of shifted winds, the strength of the lower overturning cell (maximum stream flow below 2500 m within 20-50°S) weakens by 3 Sv, from 16 to 13 Sv and the mixed layer shoals in the polar Southern Ocean (Fig. 5b; Extended Data Fig. 9). The DIC concentration increases in the Southern Ocean below the sub-surface south of 60°S by about 10 mmol/m³, and the slowdown of the lower cell causes DIC to increase throughout the deep ocean below 1.5km (Fig. 5c; Extended Data Fig. 10), even if the magnitude of the increase is small due to the short duration of the experiment. The increase in upwelling between 60°S and 50°S leads to a DIC increase in the upper 2000m depth in that region, but causes the overturning of the upper cell to increase by 6.4 SV (maximum overturning between 500-1000m and 20°S-40°S; Fig. 5b; Extended Data Fig. 9) which results in a DIC decrease in the intermediate depths throughout the ocean north of 35°S (Fig. 5c; Extended Data Fig. 10). These DIC changes are accompanied by similar changes in remineralised nitrogen, and inverse changes in dissolved O₂, highlighting the role of changes in oceanic circulation.
The reduction in upwelling within the polar Southern Ocean causes the preformed nitrate concentration of Antarctic Bottom Water to decrease. This may be enhanced by a shoaling of the polar Southern Ocean mixed layer depth within the simulation (Extended Data Fig. 9), which would reduce light limitation, providing a potential mechanism to further increase the utilisation of the upwelled nutrients¹³. Driven by a reduction in water column mixing due to a shoaling of the mixed layer in deepwater formation regions of the North Atlantic (not shown), the preformed nitrate concentration of the northern end member also decreases. As such, the preformed nitrate concentration begins to decrease throughout the deep ocean, with a 1.3% increase in global mean N* (N* = regenerated NO₃/total NO₃) after 125 years (Extended Data Fig. 10). The short duration of the Perturbed experiment inhibits a full quantification of the CO₂ response. However, by extrapolating the initial changes in endmember preformed NO₃ (Nₚre, given in mmol/m³) based on the relative volume of the ocean they represent (V) we can broadly estimate the magnitude of CO₂ change implicated by the initial changes in endmember preformed nitrate.

We take the Southern Ocean endmember as an average through the water column between 60-80°S, and the North Atlantic endmember as an average through the water column between 60-70°N in the North Atlantic. In the control run Nₚre_SO_ctr = 22.12 and Nₚre_NA_ctr = 8.5. Global mean Nₚre in the control run (Nₚre_ctr) is 15.45. Using the global mean and endmember values we calculate the volumetric contribution of the Southern Ocean endmember (Vₙₐ) as 0.51 and the North Atlantic endmember (Vₙₐ) as 0.49. Total nitrate (Nₜt) is 33.84 mmol/m³ such that globally averaged N* in the control run is 54%.

After 125 years of the Perturbed experiment Nₚre in both the endmembers decreases: Nₚre_SO_125 = 20.31 and Nₚre_NA_125 = 7.3. Assuming the same volumetric contributions as the control experiment, we calculate the expected change in global mean Nₚre in the Perturbed run (Nₚre EQSH) as,
Nₚre EQSH = Vₙₐ * Nₚre SO_125 + Vₙₐ * Nₚre_NA_125 = 13.9, equivalent to a globally averaged N* value of 59%.

Based on the initial changes in Nₚre within the endmember regions we would thus expect a global mean N* increase of ~5% in the Perturbed experiment relative to the Control, once these anomalies had propagated through the deep ocean. Applying the sensitivity of atmospheric CO₂ to global preformed nutrients of ref¹⁵, this increase in N* within the Perturbed experiment equates to an atmospheric CO₂ decrease of ~15 ppm. However, timeseries of the endmember Nₚre values indicate they are not yet equilibrated and are still decreasing after 125 years such that this likely represents a conservative estimate of the increase N* and associated lowering of CO₂ we would expect if the Perturbed experiment was run to equilibrium.

Further to this effect, the model indicates a redistribution of the regenerated nutrient and carbon pools from the intermediate depths toward the abyss. The deepening of the regenerated nutrient and carbon pools would drive a further CO₂ decrease via carbonate compensation¹⁶,³⁷. The higher DIC content at depth would lower the carbonate ion concentration, increasing CaCO₃ dissolution within the ocean. This is turn would result in a transient increase in alkalinity, lowering atmospheric CO₂. While the scaling between preformed nutrients and atmospheric CO₂ used above (ref³⁷,¹²⁴) accounts for a linear approximation of CaCO₃ dissolution following the total increase in regenerated carbon within
728

the ocean, the effect of nutrient deepening is not accounted for by this scaling. Based on the
initial changes we see in deep ocean DIC in the model (Extended Data Fig. 10), and the
scaling between intermediate-deep DIC and atmospheric CO₂ given in given refs³⁷,¹²⁴ we
would thus expect a further substantial decrease in CO₂.

729

730

731

732

Data availability

733

734

The new and compiled δ¹⁸O data are given in table Table S1 and are available on Pangaea (DOI pending). The westerly wind reconstructions generated in this study are given in table Table S2. The PMIP and CMIP data are available from https://esgf-node.llnl.gov/projects/esgf-llnl/.

735

736

737

738

Code availability

739

740

The R code and data used to perform the ΔLat<sub>SST</sub> analysis is provided on Github
(https://github.com/willyrgray/SOd18O).

741
<figure><img src="image_9.png" /><figcaption>Lat (°N)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Lat (°N)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Lat (°N)</figcaption></figure>

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

742

743

744

745

746

747

748

**Extended Data 1. Climatological δ¹⁸O.** Climatological δ¹⁸O<sub>water</sub> (ref⁹⁸), δ¹⁸O<sub>calcite-water</sub> (calculated using refs¹⁰⁴,¹²⁵), and δ¹⁸O<sub>calcite</sub> (note, colour scale is the same for panels). Symbols show location of core sites and species of planktic foraminifera (circles = *N. pachyderma*, triangles = *G. bulloides*, squares = *G. ruber*). Black line shows the position of westerlies as determined by maximum zonal wind speed (ref⁵¹) and grey line shows zero sea surface height anomaly ref(¹²⁶). Zonal-mean climatological δ¹⁸O<sub>water</sub>, δ¹⁸O<sub>calcite-water</sub>, and δ¹⁸O<sub>calcite</sub> (as shown on Fig. 1), zonal wind speed, and sea surface height.
<figure><img src="image_13.png" /><figcaption>**Extended Data 2. Indian-Pacific meridional $\delta^{18}$O profiles.** Indian-Pacific meridional $\delta^{18}$O$_{ivc-gtc}$ data with GAM fits at 500-year time steps. The GAM fit at 10 ka is shown in black. The grey box is the window in which $\Delta\text{Lat}_{\text{SST}}$ is calculated and the thick lines show the portion of the curves falling within this window.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Extended Data 3. Atlantic meridional $\delta^{18}\text{O}$ profiles. As above, but for the Atlantic.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Extended Data 4. ΔLat$_{SST}$ within different sectors and core weightings in the Indian-Pacific sectors.** Change in the SST front latitude (ΔLat$_{SST}$) from 20-10 ka using (a) all data from across the Southern Ocean, and in the (b) Indian-Pacific and (c) Atlantic sectors separately. The 5th, 32nd, 50th, 68th, and 95th percentiles are indicated. Contribution of each record to the Indian-Pacific ΔLat$_{SST}$ reconstruction based on leave-one-out analysis (Methods). (d-e) 10-20 ka reconstruction (f-g) 6.5-22 ka reconstruction. Note Antarctic marginal sites from the Atlantic sector are also included given the paucity of data from south of 65°S. Symbols distinguish species of planktic foraminifera (circles = *N. pachyderma*, triangles = *G. bulloides*, squares = *G. ruber*).</figcaption></figure>

757

758

759

760

761

762

763

764

765

766

767

768

769

770

771

772

773

774

775

776

777

778
<figure><img src="image_16.png" /><figcaption>**Extended Data 5. PMIP3/4 and CMIP5/6 ensemble.** (a) Relationship between the wind latitude (latitude of maximum zonal-mean $τ_u$) and SST front latitude (latitude of maximum $∂SST/∂Lat$, averaged over a 10° latitudinal window) in the model ensemble across the Indian and Pacific sectors. The MPI model sits as an outlier from this the rest of the ensemble, possibly due to known SST temperature biases in the Indian Ocean¹¹³, and is excluded from the regression. Including the MPI model has a negligible impact on our results (Methods). (b) As (a) but for the Atlantic sector. (c) As (a) but for the eastern Atlantic sector, where the vast majority of the mid-latitude sites in the Atlantic are located. (d) Relationship between the wind latitude and zonal-mean wind stress at 60°S (60S $τ_u$) in the model ensemble (e) Relationship between the wind latitude and the wind strength (maximum zonal-mean $τ_u$) in the model ensemble.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Extended Data 6. 22-6.5 ka analysis. (top)** Changes in the wind latitude (ΔLat<sub>wind</sub>) from 22-6.5 ka. The higher uncertainties compared to the 20-10 ka reconstruction is due to the lower number of sites than span this interval. The 5<sup>th</sup>, 32<sup>nd</sup>, 50<sup>th</sup>, 68<sup>th</sup>, and 95<sup>th</sup> percentiles are indicated. The grey dashed line shows the same analysis extended to 2.5 ka; while the uncertainties are very large due to the limited number of cores than span this interval, the results indicate little change in the position of the wind latitude during the Holocene (bottom). Comparison of reconstructed Indian-Pacific ΔLat<sub>SST</sub> using the 10-20 ka subset and 6.5-22 ka subset of cores. Open circles and the dashed grey curve correspond to the 6.5-22 ka reconstruction uncorrected; filled circles and solid grey line correspond to the 6.5-22 ka reconstruction with a 1/0.95 correction applied (see Methods). The pink curve shows the 10-20 ka reconstruction.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Extended Data 7. LGM-PI change in the PMIP3/4 ensemble.** LGM-PI change in position of maximum zonal-mean zonal wind stress in the PIMP3 (circles) and PMIP4 (triangles) ensemble, compared to the reconstructed change in $\Delta\text{Lat}_{\text{wind}}$ between 20-6.5 ka.</figcaption></figure>

802

803

804

805

806
<figure><img src="image_19.png" /></figure>

Age (ka)

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

Age (ka)

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

Age (ka)

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

Age (ka)

<caption>Extended Data 8. Deglacial changes in wind latitude, wind strength, wind stress at 60°S and northward Ekman transport at 60°S. Reconstructed changes in (a) the wind latitude (ΔLat<sub>wind</sub>) (b)</caption>

the wind strength (max τ<sub>u</sub>) (c) zonal-mean wind stress at 60°S (60S τ<sub>u</sub>) (d) northward Ekman transport at 60°S. (b) and (c) are calculated using the reconstructed changes in wind latitude (a) and the relationships between wind latitude and these parameters in the model ensemble (Extended Data Fig. 5). Northward Ekman transport at 60°S is calculated using zonal wind stress at 60°S (Methods).
<figure><img src="image_23.png" /><figcaption>**Extended Data 9. MOM5-SIS-Wombat circulation results. (a-f) Meridional overturning streamfunction in the Control (a,c,e) and Perturbed experiments (b,d,f), for the global ocean (a,b), Atlantic (c,d) and Indo-Pacific (e,f) oceans. (g) Zonally averaged mixed layer depth (MLD) in the Control (black) and Perturbed (red) simulations.**</figcaption></figure>

815

816

817

818

819

820
821

<figure><img src="image_24.png" /><figcaption>Regenerated Nitrate anomalies</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Preformed Nitrate anomalies</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Dissolved O2 anomalies</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>DIC anomalies</figcaption></figure>

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

822

823

824

825

826

827

828

**Extended Data 10. MOM5-SIS-Wombat biogeochemical results.** Zonally averaged anomalies (Perturbed-Control) in (a) regenerated NO₃ (mmol/m³), (b) preformed NO₃ (mmol/m³), (c) O₂ (mmol/m³), and (d) DIC(mmol/m³). (e) Vertically integrated (over the upper 149 m) gross phytoplankton production (molC/m²/yr) anomalies (Perturbed-Control).
829

References

830

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

831

832

2. Toggweiler, J. R., Russell, J. L. & Carson, S. R. Midlatitude westerlies, atmospheric CO2, and climate change during the ice ages. *Paleoceanography* **21**, (2006).

833

834

3. Sigman, D. M. & Boyle, E. A. Glacial/interglacial variations in atmospheric carbon dioxide. *Nature* **407**, 859–869 (2000).

835

836

4. Anderson, R. F. et al. Wind-Driven Upwelling in the Southern Ocean and the Deglacial Rise in Atmospheric CO2. *Science* **323**, 1443–1448 (2009).

837

838

5. Kohfeld, K. E. et al. Southern Hemisphere westerly wind changes during the Last Glacial Maximum: paleo-data synthesis. *Quat. Sci. Rev.* **68**, 76–95 (2013).

839

840

6. Sime, L. C. et al. Southern Hemisphere westerly wind changes during the Last Glacial Maximum: model-data comparison. *Quat. Sci. Rev.* **64**, 104–120 (2013).

841

842

7. Sime, L. C. et al. Sea ice led to poleward-shifted winds at the Last Glacial Maximum: the influence of state dependency on CMIP5 and PMIP3 models. *Clim. Past* **12**, 2241–2253 (2016).

843

844

8. Yin, J. H. A consistent poleward shift of the storm tracks in simulations of 21st century climate. *Geophys. Res. Lett.* **32**, (2005).

845

846

9. Chen, G., Lu, J. & Frierson, D. M. W. Phase Speed Spectra and the Latitude of Surface Westerlies: Interannual Variability and Global Warming Trend. *J. Clim.* **21**, 5942–5959 (2008).

847

848

10. Goyal, R., Gupta, A. S., Jucker, M. & England, M. H. Historical and Projected Changes in the Southern Hemisphere Surface Westerlies. *Geophys. Res. Lett.* **48**, e2020GL090849 (2021).

849

850

11. Toggweiler, J. R. & Samuels, B. Effect of drake passage on the global thermohaline circulation. *Deep Sea Res. Part Oceanogr. Res. Pap.* **42**, 477–500 (1995).

851

852

12. Abernathey, R., Marshall, J. & Ferreira, D. The Dependence of Southern Ocean Meridional Overturning on Wind Stress. *J. Phys. Oceanogr.* **41**, 2261–2278 (2011).

853

854

13. Martin, J. H., Gordon, R. M. & Fitzwater, S. E. Iron in Antarctic waters. *Nature* **345**, 156–158 (1990).

855

856

14. Mitchell, B. G., Brody, E. A., Holm-Hansen, O., McClain, C. & Bishop, J. Light limitation of phytoplankton biomass and macronutrient utilization in the Southern Ocean. *Limnol. Oceanogr.* **36**, 1662–1677 (1991).

857

858
859

15. Ito, T. & Follows, M. J. Preformed phosphate, soft tissue pump and atmospheric CO2. *J. Mar. Res.* **63**, 813–839 (2005).

860

861

16. Toggweiler, J. R. Variation of atmospheric CO2 by ventilation of the ocean's deepest water. *Paleoceanography* **14**, 571–588 (1999).

862

863

17. Lauderdale, J. M., Williams, R. G., Munday, D. R. & Marshall, D. P. The impact of Southern Ocean residual upwelling on atmospheric CO2 on centennial and millennial timescales. *Clim. Dyn.* **48**, 1611–1631 (2017).

864

865

866

18. Ai, X. E. et al. Southern Ocean upwelling, Earth's obliquity, and glacial-interglacial atmospheric CO2 change. *Science* **370**, 1348–1352 (2020).

867

868

19. Gray, W. R. et al. Wind-Driven Evolution of the North Pacific Subpolar Gyre Over the Last Deglaciation. *Geophys. Res. Lett.* **47**, e2019GL086328 (2020).

869

870

20. Kageyama, M. et al. The PMIP4-CMIP6 Last Glacial Maximum experiments: preliminary results and comparison with the PMIP3-CMIP5 simulations. *Clim. Past Discuss.* 1–37 (2020) doi:https://doi.org/10.5194/cp-2019-169.

871

872

873

21. Li, C. & Battisti, D. S. Reduced Atlantic Storminess during Last Glacial Maximum: Evidence from a Coupled Climate Model. *J. Clim.* **21**, 3561–3579 (2008).

874

875

22. Chavaillaz, Y., Codron, F. & Kageyama, M. Southern westerlies in LGM and future (RCP4.5) climates. *Clim. Past* **9**, 517–524 (2013).

876

877

23. Buizert, C. et al. Abrupt ice-age shifts in southern westerly winds and Antarctic climate forced from the north. *Nature* **563**, 681–685 (2018).

878

879

24. Bereiter, B. et al. Revision of the EPICA Dome C CO2 record from 800 to 600 kyr before present. *Geophys. Res. Lett.* **42**, 542–549 (2015).

880

881

25. Dong, S., Sprintall, J. & Gille, S. T. Location of the Antarctic Polar Front from AMSR-E Satellite Sea Surface Temperature Measurements. *J. Phys. Oceanogr.* **36**, 2075–2089 (2006).

882

883

26. Yang, H. et al. Tropical Expansion Driven by Poleward Advancing Midlatitude Meridional Temperature Gradients. *J. Geophys. Res. Atmospheres* **125**, e2020JD033158 (2020).

884

885

27. Gnanadesikan, A. & Hallberg, R. W. On the Relationship of the Circumpolar Current to Southern Hemisphere Winds in Coarse-Resolution Ocean Models. *J. Phys. Oceanogr.* **30**, 2013–2034 (2000).

886

887
888

28. Gordon, A. L., Molinelli, E. & Baker, T. Large-scale relative dynamic topography of the Southern Ocean. *J. Geophys. Res. Oceans* **83**, 3023–3032 (1978).

889

890

29. Gherardi, J.-M. et al. Glacial-interglacial circulation changes inferred from 231Pa/230Th sedimentary record in the North Atlantic region. *Paleoceanography* **24**, (2009).

891

892

30. Stocker, T. F. & Johnsen, S. J. A minimum thermodynamic model for the bipolar seesaw. *Paleoceanography* **18**, (2003).

893

894

31. Shakun, J. D. et al. Global warming preceded by increasing carbon dioxide concentrations during the last deglaciation. *Nature* **484**, 49–54 (2012).

895

896

32. Gottschalk, J. et al. Mechanisms of millennial-scale atmospheric CO2 change in numerical model simulations. *Quat. Sci. Rev.* **220**, 30–74 (2019).

897

898

33. Menviel, L. et al. Southern Hemisphere westerlies as a driver of the early deglacial atmospheric CO2 rise. *Nat. Commun.* **9**, 2503 (2018).

899

900

34. Hogg, A. M., Spence, P., Saenko, O. A. & Downes, S. M. The Energetics of Southern Ocean Upwelling. *J. Phys. Oceanogr.* **47**, 135–153 (2017).

901

902

35. Spence, P., Saenko, O. A., Eby, M. & Weaver, A. J. The Southern Ocean Overturning: Parameterized versus Permitted Eddies. *J. Phys. Oceanogr.* **39**, 1634–1651 (2009).

903

904

36. Hallberg, R. & Gnanadesikan, A. The Role of Eddies in Determining the Structure and Response of the Wind-Driven Southern Hemisphere Overturning: Results from the Modeling Eddies in the Southern Ocean (MESO) Project. *J. Phys. Oceanogr.* **36**, 2232–2252 (2006).

905

906

907

37. Boyle, E. A. Vertical oceanic nutrient fractionation and glacial/interglacial CO2 cycles. *Nature* **331**, 55–56 (1988).

908

909

38. Du, J., Haley, B. A. & Mix, A. C. Evolution of the Global Overturning Circulation since the Last Glacial Maximum based on marine authigenic neodymium isotopes. *Quat. Sci. Rev.* **241**, 106396 (2020).

910

911

912

39. Rae, J. W. B. et al. CO2 storage and release in the deep Southern Ocean on millennial to centennial timescales. *Nature* **562**, 569–573 (2018).

913

914

40. Jaccard, S. L. & Galbraith, E. D. Large climate-driven changes of oceanic oxygen concentrations during the last deglaciation. *Nat. Geosci.* **5**, 151–156 (2012).

915

916

41. Hoogakker, B. A. A. et al. Glacial expansion of oxygen-depleted seawater in the eastern tropical Pacific. *Nature* **562**, 410–413 (2018).

917
42. Peterson, C. D. & Lisiecki, L. E. Deglacial carbon cycle changes observed in a compilation of 127 benthic δ¹³C time series (20–6&thinsp;ka). *Clim. Past* **14**, 1229–1252 (2018).

43. Jaccard, S. L. et al. Two Modes of Change in Southern Ocean Productivity Over the Past Million Years. *Science* **339**, 1419–1423 (2013).

44. Bishop, S. P. et al. Southern Ocean Overturning Compensation in an Eddy-Resolving Climate Simulation. *J. Phys. Oceanogr.* **46**, 1575–1592 (2016).

45. Dufour, C. O. et al. Role of Mesoscale Eddies in Cross-Frontal Transport of Heat and Biogeochemical Tracers in the Southern Ocean. *J. Phys. Oceanogr.* **45**, 3057–3081 (2015).

46. Burke, A. & Robinson, L. F. The Southern Ocean's Role in Carbon Exchange During the Last Deglaciation. *Science* **335**, 557–561 (2012).

47. Martínez-Botí, M. A. et al. Boron isotope evidence for oceanic carbon dioxide leakage during the last deglaciation. *Nature* **518**, 219–222 (2015).

48. Menviel, L. & Joos, F. Toward explaining the Holocene carbon dioxide and carbon isotope records: Results from transient ocean carbon cycle-climate simulations. *Paleoceanography* **27**, (2012).

49. Fogwill, C. J., Turney, C. S. M., Hutchinson, D. K., Taschetto, A. S. & England, M. H. Obliquity Control On Southern Hemisphere Climate During The Last Glacial. *Sci. Rep.* **5**, 11673 (2015).

50. Gruber, N., Landschützer, P. & Lovenduski, N. S. The Variable Southern Ocean Carbon Sink. *Annu. Rev. Mar. Sci.* **11**, 159–186 (2019).

51. Kalnay, E. et al. The NCEP/NCAR 40-Year Reanalysis Project. *Bull. Am. Meteorol. Soc.* **77**, 437–472 (1996).

52. Haddam, N. A., Michel, E., Siani, G., Licari, L. & Dewilde, F. Ventilation and Expansion of Intermediate and Deep Waters in the Southeast Pacific During the Last Termination. *Paleoceanogr. Paleoclimatology* **35**, e2019PA003743 (2020).

53. Gottschalk, J. et al. Glacial heterogeneity in Southern Ocean carbon storage abated by fast South Indian deglacial carbon release. *Nat. Commun.* **11**, 6192 (2020).

54. Sarnthein, M. et al. Changes in East Atlantic Deepwater Circulation over the last 30,000 years: Eight time slice reconstructions. *Paleoceanography* **9**, 209–267 (1994).
946

55. Wang, Y. V. et al. Northern and southern hemisphere controls on seasonal sea surface temperatures in the Indian Ocean during the last deglaciation. *Paleoceanography* **28**, 619–632 (2013).

947

948

56. Stuut, J.-B. W. et al. A 300-kyr record of aridity and wind strength in southwestern Africa: inferences from grain-size distributions of sediments on Walvis Ridge, SE Atlantic. *Mar. Geol.* **180**, 221–233 (2002).

949

950

951

57. Levi, C. et al. Low-latitude hydrological cycle and rapid climate changes during the last deglaciation. *Geochem. Geophys. Geosystems* **8**, (2007).

952

953

58. Schneider, R. R., Müller, P. J. & Ruhland, G. Late Quaternary surface circulation in the east equatorial South Atlantic: Evidence from Alkenone sea surface temperatures. *Paleoceanography* **10**, 197–219 (1995).

954

955

956

59. Stuut, J.-B. W. et al. A 5.3-Million-Year History of Monsoonal Precipitation in Northwestern Australia. *Geophys. Res. Lett.* **46**, 6946–6954 (2019).

957

958

60. Bostock, H. C., Opdyke, B. N., Gagan, M. K. & Fifield, L. K. Carbon isotope evidence for changes in Antarctic Intermediate Water circulation and ocean ventilation in the southwest Pacific during the last deglaciation. *Paleoceanography* **19**, (2004).

959

960

961

61. Santos, T. P. et al. Prolonged warming of the Brazil Current precedes deglaciations. *Earth Planet. Sci. Lett.* **463**, 1–12 (2017).

962

963

62. Caley, T. et al. High-latitude obliquity as a dominant forcing in the Agulhas current system. *Clim. Past* **7**, 1285–1296 (2011).

964

965

63. Portilho-Ramos, R. C. et al. Methane release from the southern Brazilian margin during the last glacial. *Sci. Rep.* **8**, 5948 (2018).

966

967

64. Pereira, L. S., Arz, H. W., Pätzold, J. & Portilho-Ramos, R. C. Productivity Evolution in the South Brazilian Bight During the Last 40,000 Years. *Paleoceanogr. Paleoclimatology* **33**, 1339–1356 (2018).

968

969

970

65. Carlson, A. E. et al. Subtropical Atlantic salinity variability and Atlantic meridional circulation during the last deglaciation. *Geology* **36**, 991–994 (2008).

971

972

66. Scussolini, P. & Peeters, F. J. C. A record of the last 460 thousand years of upper ocean stratification from the central Walvis Ridge, South Atlantic. *Paleoceanography* **28**, 426–439 (2013).

973

974

975
976

67. Chiessi, C. M. et al. Variability of the Brazil Current during the late Holocene. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **415**, 28–36 (2014).

977

68. Winn, K. Carbon and oxygen isotope ratios on planktonic foraminifera in subtropical Southeast Pacific core GIK17747-2. 333 data points (2013) doi:10.1594/PANGAEA.815877.

979

69. Mohtadi, M. et al. Deglacial pattern of circulation and marine productivity in the upwelling region off central-south Chile. *Earth Planet. Sci. Lett.* **272**, 221–230 (2008).

980

981

70. Dyez, K. A., Zahn, R. & Hall, I. R. Multicentennial Agulhas leakage variability and links to North Atlantic climate during the past 80,000 years. *Paleoceanography* **29**, 1238–1248 (2014).

982

983

71. Schiraldi, B., Sikes, E. L., Elmore, A. C., Cook, M. S. & Rose, K. A. Southwest Pacific subtropics responded to last deglacial warming with changes in shallow water sources. *Paleoceanography* **29**, 595–611 (2014).

984

985

986

72. Martínez-Méndez, G. et al. Contrasting multiproxy reconstructions of surface ocean hydrography in the Agulhas Corridor and implications for the Agulhas Leakage during the last 345,000 years. *Paleoceanography* **25**, (2010).

987

988

989

73. Calvo, E., Pelejero, C., Deckker, P. D. & Logan, G. A. Antarctic deglacial pattern in a 30 kyr record of sea surface temperature offshore South Australia. *Geophys. Res. Lett.* **34**, (2007).

990

991

74. Nelson, C. S., Hendy, I. L., Neil, H. L., Hendy, C. H. & Weaver, P. P. E. Last glacial jetting of cold waters through the Subtropical Convergence zone in the Southwest Pacific off eastern New Zealand, and some geological implications. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **156**, 103–121 (2000).

992

993

994

995

75. Gottschalk, J. et al. Past Carbonate Preservation Events in the Deep Southeast Atlantic Ocean (Cape Basin) and Their Implications for Atlantic Overturning Dynamics and Marine Carbon Cycling. *Paleoceanogr. Paleoclimatology* **33**, 643–663 (2018).

996

997

998

76. Lamy, F. et al. Antarctic Timing of Surface Water Changes off Chile and Patagonian Ice Sheet Response. *Science* **304**, 1959–1962 (2004).

999

1000

77. Haddam, N. A. et al. Changes in latitudinal sea surface temperature gradients along the Southern Chilean margin since the last glacial. *Quat. Sci. Rev.* **194**, 62–76 (2018).

1001

1002

78. Hodell, D. A., Venz, K. A., Charles, C. D. & Ninnemann, U. S. Pleistocene vertical carbon isotope and carbonate gradients in the South Atlantic sector of the Southern Ocean. *Geochem. Geophys. Geosystems* **4**, 1–19 (2003).

1003

1004

1005
1006 79. Pichon, J.-J. et al. Surface water temperature changes in the high latitudes of the southern
hemisphere over the Last Glacial-Interglacial Cycle. *Paleoceanography* **7**, 289–318 (1992).

1007

1008 80. Rickaby, R. E. M. & Elderfield, H. Planktonic foraminiferal Cd/Ca: Paleonutrients or
paleotemperature? *Paleoceanography* **14**, 293–303 (1999).

1009

1010 81. Gottschalk, J., Skinner, L. C. & Waelbroeck, C. Contribution of seasonal sub-Antarctic surface
water variability to millennial-scale changes in atmospheric CO2 over the last deglaciation and
Marine Isotope Stage 3. *Earth Planet. Sci. Lett.* **411**, 87–99 (2015).

1011

1012

1013 82. Sikes, E. L. et al. Southern Ocean seasonal temperature and Subtropical Front movement on the
South Tasman Rise in the late Quaternary. *Paleoceanography* **24**, (2009).

1014

1015 83. Sicre, M. A. et al. Mid-latitude Southern Indian Ocean response to Northern Hemisphere Heinrich
events. *Earth Planet. Sci. Lett.* **240**, 724–731 (2005).

1016

1017 84. Labeyrie, L. et al. Hydrographic changes of the Southern Ocean (southeast Indian Sector) Over
the last 230 kyr. *Paleoceanography* **11**, 57–76 (1996).

1018

1019 85. Govin, A. et al. Evidence for northward expansion of Antarctic Bottom Water mass in the
Southern Ocean during the last glacial inception. *Paleoceanography* **24**, (2009).

1020

1021 86. Bostock, H. C., Hayward, B. W., Neil, H. L., Sabaa, A. T. & Scott, G. H. Changes in the position of
the Subtropical Front south of New Zealand since the last glacial period. *Paleoceanography* **30**,
824–844 (2015).

1022

1023

1024 87. Hodell, D. A., Charles, C. D. & Ninnemann, U. S. Comparison of interglacial stages in the South
Atlantic sector of the southern ocean for the past 450 kyr: implifications for Marine Isotope Stage
(MIS) 11. *Glob. Planet. Change* **24**, 7–26 (2000).

1025

1026

1027 88. Charles, C. D., Froelich, P. N., Zibello, M. A., Mortlock, R. A. & Morley, J. J. Biogenic opal in
Southern Ocean sediments over the last 450,000 years: Implications for surface water chemistry
and circulation. *Paleoceanography* **6**, 697–728 (1991).

1028

1029

1030 89. *Use of Proxies in Paleoceanography: Examples from the South Atlantic.* (Springer-Verlag, 1999).
doi:10.1007/978-3-642-58646-0.

1031

1032 90. Caniupán, M. et al. Millennial-scale sea surface temperature and Patagonian Ice Sheet changes
off southernmost Chile (53°S) over the past ~60 kyr. *Paleoceanography* **26**, (2011).

1033

1034 91. Hasenfratz, A. P. et al. The residence time of Southern Ocean surface waters and the 100,000-
year ice age cycle. *Science* **363**, 1080–1084 (2019).
1035

1036
1036

92. *Proceedings of the Ocean Drilling Program, 177 Scientific Results*. vol. 177 (Ocean Drilling Program, 2003).

1037

1038

93. Mashiotta, T. A., Lea, D. W. & Spero, H. J. Glacial-interglacial changes in Subantarctic sea surface temperature and δ18O-water using foraminiferal Mg. *Earth Planet. Sci. Lett.* **170**, 417–432 (1999).

1039

1040

1041

94. Crosta, X., Sturm, A., Armand, L. & Pichon, J.-J. Late Quaternary sea ice history in the Indian sector of the Southern Ocean as recorded by diatom assemblages. *Mar. Micropaleontol.* **50**, 209–223 (2004).

1042

1043

1044

95. *Carbon Cycling in the Glacial Ocean: Constraints on the Ocean's Role in Global Change: Quantitative Approaches in Paleoceanography*. (Springer-Verlag, 1994). doi:10.1007/978-3-642-78737-9.

1045

1046

1047

96. Grobe, H. & Mackensen, A. Late Quaternary climatic cycles as recorded in sediments from the Antarctic continental margin. in *Antarctic Research Series* (eds. Kennett, J. P. & Warkne, D. A.) vol. 56 349–376 (American Geophysical Union, 1992).

1048

1049

1050

97. Lu, Z. et al. Oxygen depletion recorded in upper waters of the glacial Southern Ocean. *Nat. Commun.* **7**, 11146 (2016).

1051

1052

98. LeGrande, A. N. & Schmidt, G. A. Global gridded data set of the oxygen isotopic composition in seawater. *Geophys. Res. Lett.* **33**, L12604 (2006).

1053

1054

99. Wood, S. N. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. *J. R. Stat. Soc. Ser. B Stat. Methodol.* **73**, 3–36 (2011).

1055

1056

1057

100. Wood, S. N., Pya, N. & Säfken, B. Smoothing Parameter and Model Selection for General Smooth Models. *J. Am. Stat. Assoc.* **111**, 1548–1563 (2016).

1058

1059

101. Simpson, G. L. Modelling Palaeoecological Time Series Using Generalised Additive Models. *Front. Ecol. Evol.* **6**, (2018).

1060

1061

102. Schrag, D. P. et al. The oxygen isotopic composition of seawater during the Last Glacial Maximum. *Quat. Sci. Rev.* **21**, 331–342 (2002).

1062

1063

103. Lambeck, K., Rouby, H., Purcell, A., Sun, Y. & Sambridge, M. Sea level and global ice volumes from the Last Glacial Maximum to the Holocene. *Proc. Natl. Acad. Sci.* **111**, 15296–15303 (2014).

1064

1065
1066 104. Kim, S.-T. & O'Neil, J. R. Equilibrium and nonequilibrium oxygen isotope effects in synthetic carbonates. *Geochim. Cosmochim. Acta* **61**, 3461–3475 (1997).

1068 105. Malevich, S. B., Vetter, L. & Tierney, J. E. Global Core Top Calibration of δ18O in Planktic Foraminifera to Sea Surface Temperature. *Paleoceanogr. Paleoclimatology* **34**, 1292–1315 (2019).

1071 106. Efron, B. Bootstrap Methods: Another Look at the Jackknife. *Ann. Stat.* **7**, 1–26 (1979).

1072 107. Saunders, K. M. et al. Holocene dynamics of the Southern Hemisphere westerly winds and possible links to CO 2 outgassing. *Nat. Geosci.* **11**, 650–655 (2018).

1074 108. Lombard, F. et al. Modelling planktic foraminifer growth and distribution using an ecophysiological multi-species approach. *Biogeosciences* **8**, 853–873 (2011).

1076 109. Braconnot, P. et al. Evaluation of climate models using palaeoclimatic data. *Nat. Clim. Change* **2**, 417–424 (2012).

1078 110. Kageyama, M. et al. The PMIP4 contribution to CMIP6 – Part 1: Overview and over-arching analysis plan. *Geosci. Model Dev.* **11**, 1033–1057 (2018).

1080 111. Taylor, K. E., Stouffer, R. J. & Meehl, G. A. An Overview of CMIP5 and the Experiment Design. *Bull. Am. Meteorol. Soc.* **93**, 485–498 (2012).

1082 112. Eyring, V. et al. Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. *Geosci. Model Dev.* **9**, 1937–1958 (2016).

1084 113. Fathrio, I. et al. Assessment of western Indian Ocean SST bias of CMIP5 models. *J. Geophys. Res. Oceans* **122**, 3123–3140 (2017).

1086 114. Chapman, C. C., Lea, M.-A., Meyer, A., Sallée, J.-B. & Hindell, M. Defining Southern Ocean fronts and their influence on biological and physical processes in a changing climate. *Nat. Clim. Change* **10**, 209–219 (2020).

1089 115. Barnes, E. A. & Polvani, L. Response of the Midlatitude Jets, and of Their Variability, to Increased Greenhouse Gases in the CMIP5 Models. *J. Clim.* **26**, 7117–7135 (2013).

1091 116. McGraw, M. C. & Barnes, E. A. Seasonal Sensitivity of the Eddy-Driven Jet to Tropospheric Heating in an Idealized AGCM. *J. Clim.* **29**, 5223–5240 (2016).

1093 117. 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 (2016).
1095 118. Olsen, A. et al. The Global Ocean Data Analysis Project version 2 (GLODAPv2) – an internally consistent data product for the world ocean. *Earth Syst. Sci. Data* **8**, 297–323 (2016).

1097 119. Griffies, S. M. et al. Coordinated Ocean-ice Reference Experiments (COREs). *Ocean Model*. **26**, 1–46 (2009).

1099 120. Redi, M. H. Oceanic Isopycnal Mixing by Coordinate Rotation. *J. Phys. Oceanogr.* **12**, 1154–1154 (1982).

1101 121. Gent, P. R. & Mcwilliams, J. C. Isopycnal Mixing in Ocean Circulation Models. *J. Phys. Oceanogr.* **20**, 150–155 (1990).

1103 122. Jackett, D. & McDougall, T. A Neutral Density Variable for the World’s Oceans. (1997) doi:10.1175/1520-0485(1997)027<0237:ANDVFT>2.0.CO;2.

1105 123. Zika, J. D., England, M. H. & Sijp, W. P. The Ocean Circulation in Thermohaline Coordinates. *J. Phys. Oceanogr.* **42**, 708–724 (2012).

1107 124. Boyle, E. A. The role of vertical chemical fractionation in controlling late Quaternary atmospheric carbon dioxide. *J. Geophys. Res. Oceans* **93**, 15701–15714 (1988).

1109 125. Boyer, T. P. et al. World ocean database 2013. (2013) doi:10.7289/V5NZ85MT.

1110 126. AVISO: Satellite derived Sea Surface Height above Geoid | NCAR - Climate Data Guide.
https://climatedataguide.ucar.edu/climate-data/aviso-satellite-derived-sea-surface-height-above-

1112 geoid.