Southern Hemisphere westerlies as a driver of the early deglacial atmospheric CO2 rise

The early part of the last deglaciation is characterised by a ~40 ppm atmospheric CO2 rise occurring in two abrupt phases. The underlying mechanisms driving these increases remain a subject of intense debate. Here, we successfully reproduce changes in CO2, δ13C and Δ14C as recorded by paleo-records during Heinrich stadial 1 (HS1). We show that HS1 CO2 increase can be explained by enhanced Southern Ocean upwelling of carbon-rich Pacific deep and intermediate waters, resulting from intensified Southern Ocean convection and Southern Hemisphere (SH) westerlies. While enhanced Antarctic Bottom Water formation leads to a millennial CO2 outgassing, intensified SH westerlies induce a multi-decadal atmospheric CO2 rise. A strengthening of SH westerlies in a global eddy-permitting ocean model further supports a multi-decadal CO2 outgassing from the Southern Ocean. Our results highlight the crucial role of SH westerlies in the global climate and carbon cycle system with important implications for future climate projections.

T he natural climate variability of the last 800,000 years is dominated by glacial-interglacial cycles, with atmospheric CO 2 variations providing a major positive feedback 1 . However, the sequence of events leading to deglacial CO 2 rise remains poorly constrained and a combination of mechanisms has been invoked to explain the full~90 ppm amplitude. These include reduced CO 2 solubility, global ocean alkalinity decrease 2 , reduced iron fertilisation [3][4][5] , increased Southern Ocean ventilation 6,7 and poleward shift of the Southern Hemisphere (SH) westerlies 8,9 . Modelling studies trying to tackle the problem of glacial-interglacial CO 2 changes mostly involve idealised sensitivity studies, often performed under constant preindustrial 8,[10][11][12] or LGM boundary conditions 3,4 . A detailed study of the deglaciation would provide a more direct link between changes in the climate and carbon cycle and would allow a direct model-data comparison, thus further constraining the processes responsible for atmospheric CO 2 changes. HS1 (~17. 6-14.7 ka), at the beginning of the last deglaciation, is an important period to understand as it represents a major phase of atmospheric CO 2 rise and the transition out of the glacial period. Paleoproxy records suggest that North Atlantic Deep Water (NADW) formation weakened significantly during HS1 13 , effectively reducing the meridional heat transport to the North Atlantic and leading to cold and dry conditions over Greenland 14 and the North Atlantic 15 . In contrast, paleoproxies indicate that Antarctic surface air temperature and Southern Ocean surface waters experienced a warming of~5°C and~3°C 16,17 , respectively. This warming is partly due to increased heat content in the South Atlantic, subsequent advection of warm waters through the Antarctic Circumpolar Current and the concurrent 40 ppm atmospheric CO 2 increase. However, the mechanisms leading to HS1 CO 2 rise are still poorly constrained and an overarching mechanism linking this CO 2 rise to the North Atlantic cooling and high southern latitude warming is still missing.
Recent high-resolution Antarctic ice core records 7,18,19 show that atmospheric CO 2 rose in two major phases at~17.2 and 16.2 ka, each associated with a~0.2‰ decrease in the atmospheric carbon isotopic composition (δ 13 CO 2 ). In addition, atmospheric radiocarbon content (Δ 14 C) declined by~112‰ between 17.6 and 15 ka 20 . Explaining the atmospheric CO 2 increase thus also requires attributing the concurrent δ 13 CO 2 and Δ 14 C declines. δ 13 CO 2 integrates changes in terrestrial carbon, marine export production, oceanic circulation, and air-sea gas exchange 21 , while Δ 14 C is controlled by atmospheric 14 C production and carbon exchange between the atmosphere and abyssal ocean carbon or old terrestrial carbon.
A number of processes have been put forward to explain the early deglacial atmospheric CO 2 increase. Co-variations of iron flux and nutrient utilisation in the sub-Antarctic 5 suggest that iron fertilisation could exert a significant control on atmospheric CO 2 during HS1 through its modulation of the Southern Ocean biological pump efficiency. Modelling studies show that reduced iron fertilisation could lead to a millennial atmospheric CO 2 increase of~10 ppm coupled with a 0.1 ‰ δ 13 CO 2 decrease 3,4,19,22 , during the early deglaciation. However, iron fertilisation does not affect atmospheric Δ 14 C or ocean ventilation ages, and variations in atmospheric iron deposition to the Southern Ocean are ultimately controlled by the exposure of the continental shelves, SH hydrology and winds. Indeed, a deglacial decline in South Atlantic ventilation ages has been observed 6 , indicating a possible role of Southern Ocean ventilation in driving the deglacial CO 2 increase. This change in Southern Ocean ventilation could be modulated by the strength and position of the SH westerlies [7][8][9] . Idealised modelling studies performed under constant pre-industrial conditions have shown that stronger or poleward shifted SH westerlies could enhance deep ocean ventilation thereby leading to an atmospheric CO 2 rise and δ 13 CO 2 decline 10,11,21 . While all numerical experiments performed show that stronger SH westerlies lead to an atmospheric CO 2 increase 10,11,21,23 , the impact of changes in their latitudinal position is more ambiguous and could depend on their initial latitudinal position 23,24 . However, the latitudinal position of the SH westerlies at the LGM remains unclear 25,26 . In addition, a recent modelling study 12 , also performed under constant preindustrial boundary conditions, concluded that the SH westerlies did not lead to significant changes in atmospheric CO 2 during HS1. Instead, and even though no abrupt atmospheric CO 2 rise was simulated, the conclusion was that the HS1 atmospheric CO 2 rise was solely due to a reduced efficiency of the biological pump, resulting from a weaker Atlantic Meridional Overturning Circulation.
As the magnitude and rate associated with an oceanic carbon release to the atmosphere during the deglaciation have been questioned, a deglacial transfer of carbon from the terrestrial to atmospheric reservoir has been put forward, either as thawing of Northern Hemisphere permafrost 27 during HS1, or as a Northern Hemisphere terrestrial carbon release due to a southward shift of the Intertropical Convergence Zone (ITCZ) at 16.2 ka 19,28 . However, the global terrestrial carbon reservoir increased over the deglaciation 29,30 and the timing and magnitude of the permafrost carbon contribution remain poorly constrained.
So far, no three-dimensional transient simulation was able to reproduce the changes in atmospheric CO 2 , its isotopic composition, as well as oceanic δ 13 C and ventilation ages across HS1. Here, we explore the processes leading to the two-stage CO 2 increase during HS1, and their links to NADW weakening and Antarctic warming by performing a suite of transient experiments of HS1 with the carbon-isotopes enabled Earth System Model LOVECLIM 21 . This suite of simulations assesses the impact of Southern Ocean ventilation changes, including the potential role of buoyancy and dynamic forcing, such as meltwater and SH westerlies, in driving the rapid atmospheric CO 2 increase during HS1. The ocean carbon response to changes in SH westerlies is further assessed in a global eddy-permitting ocean model. We show that enhanced Southern Ocean convection and upwelling of Circumpolar Deep Water, driven by intensified SH westerlies, lead to an atmospheric CO 2 rise, δ 13 CO 2 and Δ 14 C decrease in agreement with paleo-records 7,19,20 .

Results
Simulating atmospheric CO 2 and δ 13 CO 2 during HS1. The transient simulation is initialised from a Last Glacial Maximum (LGM) state constrained by oceanic δ 13 C and ventilation age distributions 31 . The LGM ocean circulation is characterised by shallow NADW, relatively weak North Pacific Intermediate Water (NPIW) and very weak Antarctic Bottom Water (AABW), obtained by adding a meltwater flux into the Southern Ocean and weakening the SH windstress by 20% 31 . By forcing the model with changes in orbital parameters, Northern Hemisphere icesheet extent and albedo as well as freshwater input in the North Atlantic (Methods), the deep ocean convection in the Norwegian Sea is suppressed, thus resulting in very weak NADW formation during HS1 (Fig. 1a). Reduced moisture transport from the Atlantic to the Pacific and a deepening of the Aleutian low due to NADW cessation leads to stronger NPIW formation in agreement with paleoproxy records 32,33 (Fig. 1j).
The NADW weakening induces a southward shift of the ITCZ (Fig. 2 13 , thus indicating very weak NADW transport. However, another phase of NADW weakening 36 , probably associated with the disintegration of the Laurentide ice-sheet and Heinrich event 1, occurred at~16.2 ka 37 . A further intensification of the SH westerlies and reduced Southern Ocean freshwater flux at 16.2 ka enhances Southern Ocean convection in our simulation, resulting in stronger Antarctic Intermediate Water (AAIW) and AABW formation during HS1 (Fig. 3d, e), consistent with reduced ventilation ages in the South Atlantic 6 and the Pacific 38 (Fig. 1h, i) as well as a peak in the Southern Ocean opal flux 9 .
Enhanced formation of AAIW and AABW and the associated upwelling of Circumpolar Deep Water decrease the oceanic carbon content below~2000 m depth and particularly in the deep South Pacific, leading to CO 2 outgassing in the Southern Ocean (Fig. 2). As a result, the simulated atmospheric CO 2 increases in close agreement with high-resolution Antarctic ice core records ( Fig. 1d) 7 , with a simulated 19 ppm pCO 2 increase between 17.2 and 16.2 ka and an abrupt 9 ppm rise at 16.2 ka. Over the course of the experiment the global ocean carbon content decreases bỹ 100 GtC, while the terrestrial carbon content increases by~50 GtC ( Fig. 3g red line, Supplementary Table 1). This terrestrial carbon increase is mostly due to enhanced carbon storage in vegetation, roots and soils of the southern tropics due to a southward shift of the ITCZ, and the fertilisation effect linked to the atmospheric CO 2 increase ( Supplementary Fig. 1).
In line with ice core records 19 , each phase of atmospheric CO 2 rise is associated with~0.25‰ δ 13 CO 2 decrease (Fig. 1e). Each of these drops in δ 13 CO 2 is however~0.05‰ higher than the ones recorded in ice-core records, thus leading to a significant overshoot at~16 ka. These δ 13 CO 2 decreases are primarily due to enhanced ventilation of deep and intermediate waters with low δ 13 C signatures and upwelling in the Southern Ocean (Fig. 4a, b) 21 . This reduces the respired oceanic carbon content (Supplementary Fig. 3) and results in positive δ 13 C anomalies along the AAIW and AABW pathways, as recorded in benthic δ 13 C ( Fig. 4a, b, Supplementary Table 2). During HS1, ventilation of "old" deep waters decreases atmospheric Δ 14 C by~150‰, consistent with atmospheric Δ 14 C 20 and marine Δ 14 C reconstructions 6,32,33,38 (Fig. 1f,h-j and 4c, d, Supplementary  Fig. 4).
Finally, due to the concurrent atmospheric CO 2 rise and enhanced meridional heat transport towards high southern   latitudes, Southern Ocean sea surface temperature (SST) and Antarctic air temperature increase by up to 3°C and 4°C, respectively, between 17.6 ka and 15 ka, and the Southern Ocean sea-ice cover is significantly reduced (Figs. 1c, 2 and 3f red line). This is in agreement with idealised experiments performed with a global eddy-permitting coupled ocean sea-ice model, which shows that stronger AABW induces a Southern Ocean SST increase through enhanced poleward heat transport 39 . Enhanced AABW leads to a positive feedback between SST increase and atmospheric CO 2 rise through the solubility effect ( Fig. 5 green and red), thus explaining a major part of the early deglacial warming at mid-and high southern latitudes.
To investigate the mechanisms driving atmospheric CO 2 increase and δ 13 CO 2 decrease during HS1, three additional transient simulations (LH1, LH1-SO, LH1-SHW) are performed with LOVECLIM ( Fig. 3, Methods). Similar to the simulation presented in Fig. 1 (LH1-SO-SHW, red), the experiments start from a LGM state featuring meltwater input in the Southern Ocean and weaker SH westerlies, and include a meltwater input in the North Atlantic during HS1 (Fig. 3a-c). For LH1, Southern Ocean meltwater input and weaker SH westerlies are kept at the LGM level thus maintaining weak AAIW and AABW formation (Fig. 3, magenta). Enhanced NPIW decreases the carbon content at intermediate depths in the North Pacific (~700-2500 m depth), but NADW cessation leads to carbon accumulation in the Atlantic Ocean, mostly as respired carbon ( Supplementary Fig. 3) at intermediate depths in the North Atlantic (Fig. 2, Supplementary Fig. 5). As a result, there is little change in atmospheric CO 2 , δ 13 CO 2 or Δ 14 C (Fig. 3h,  In LH1-SO, the meltwater input in the Southern Ocean is stopped, resulting in stronger AABW, but SH westerlies and AAIW stay weak (Fig. 3d, e, green). Enhanced AABW steepens Southern Ocean isopycnals ( Supplementary Fig. 5) and strengthens deep ocean ventilation, thus decreasing the ocean carbon content below 2000 m depth, particularly in the Pacific Ocean (Supplementary Table 1) 40 . This leads to a gradual 15 ppm atmospheric CO 2 increase through Southern Ocean CO 2 outgassing ( Supplementary Fig. 7). In addition, the ventilation of 13 C and 14 C depleted deep ocean causes slow atmospheric δ 13 C and Δ 14 C decreases of 0.016 and 130‰, respectively 21 ( Fig. 3h and Supplementary Fig. 6). However, the simulated rates of change of CO 2 and δ 13 CO 2 are smaller than recorded in Antarctic ice cores and constant throughout HS1, despite a step-wise forcing. As ventilating the deep ocean is a multi-centennial process, even relatively rapid AABW changes cannot reproduce the fast atmospheric CO 2 increase at 16.2 ka.
Finally, when the SH westerly windstress is intensified but the Southern Ocean meltwater input is kept at LGM level (LH1-SHW, Fig. 3b-e, blue), both AAIW and AABW strengthen and the isopycnal slopes steepen in the Southern Ocean with deep and intermediate waters outcropping south of 50°S (Supplementary Fig. 5) 41 . This results in very abrupt CO 2 outgassing in the Southern Ocean ( Supplementary Fig. 7), and leads to a 8 ppm atmospheric CO 2 rise and 0.2 ‰ δ 13 CO 2 decline within~100 years at 16.2 ka.  Table 3) A further weakening of the oceanic circulation from the LGM into HS1 as simulated in experiment LH1 only leads to small changes in the carbon reservoirs (Supplementary Table 1), with a slight global increase in the remineralized PO 4 content, and a 2% decrease in the preformed PO 4 content ( Supplementary Fig. 3). In contrast, enhanced Southern Ocean ventilation leads to a carbon loss in the Pacific Ocean associated with a global decrease in remineralized PO 4 and a~6% increase in preformed PO 4 , thus implying a reduced efficiency of the biological pump in all experiments with enhanced AABW formation (LH1-SO, LH1-SHW and LH1-SO-SHW).
However, the processes leading to the atmospheric CO 2 increase also depend on the nature of the forcing: buoyancy (LH1-SO) or dynamic (LH1-SHW) forcing. Enhanced AABW via reduced surface freshwater flux leads to a global SST increase coupled to an increase in surface salinity, which induces a 8 ppm CO 2 rise through the solubility effect (Fig. 5, green triangle). In contrast, no notable CO 2 increase occurs through the solubility effect in the wind only forcing case (Fig. 5, blue triangle). While stronger deep ocean convection triggered by changes in the buoyancy forcing increases the southward baroclinic flow and thus the oceanic meridional heat transport to high southern latitudes 39 , the SST response to changes in windstress is more complex. Stronger windstress over the Southern Ocean increases the Ekman transport and enhances the oceanic heat loss to the atmosphere, thus triggering an initial cooling in some areas of the Southern Ocean. However, after~20 years, the eddy-driven poleward flow and enhanced upwelling of relatively warm Circumpolar Deep Water reverse the cooling trend into a warming trend 42 . In addition, an intensification of SH westerlies leads to an increase in the ventilation of low salinity AAIW, resulting in a global sea surface salinity (SSS) decrease.
The CO 2 rise in the wind forced simulation is due to an imbalance between Dissolved Inorganic Carbon (DIC) and alkalinity (ALK) increase (Fig. 5, blue diamond). As seen in Fig. 2, stronger AABW and enhanced upwelling in the Southern Ocean lead to a DIC transfer from deep and intermediate depths of the Pacific and Southern Oceans up to the surface of the Southern Ocean. While an increase in surface DIC leads to an atmospheric CO 2 rise, an increase in surface ALK has an opposite effect. DIC and ALK thus tend to compensate each other as their oceanic distributions are similar (Supplementary Fig. 9) and their fractional impacts on CO 2 are of similar magnitude but opposite (Methods). When AABW is enhanced through reduced buoyancy forcing, this imbalance leads to a 5 ppm CO 2 increase, while it causes the 12 ppm CO 2 increase in the wind forcing case (Fig. 5, blue and green diamonds). This difference is mainly due to enhanced formation of AAIW in the simulation with stronger SH westerlies (Fig. 3e). Stronger AAIW formation entrains intermediate depth waters, characterised by low alkalinity, to the surface. In addition, this process acts on a faster timescale than deep ocean ventilation through enhanced AABW.
Both δ 13 CO 2 and Δ 14 C decrease are abrupt and of higher amplitude in the wind forcing experiment because of ventilation of intermediate waters and enhanced air-sea gas exchange, which is particularly important for carbon isotopes. When multimillennial deep ocean ventilation and solubility decrease are combined with decadal scale changes in intermediate waters (Fig. 5 red, LH1-SO-SHW), the atmospheric CO 2 increase is in close agreement with ice core records (Fig. 1d).
The atmospheric Δ 14 C changes shown in Fig. 1f and S6 were obtained by keeping the atmospheric 14 C production rate constant at LGM levels (Methods). Our results thus suggest that most, if not all, of the atmospheric Δ 14 C changes can be attributed to changes in ocean circulation and air-sea gas exchange with a smaller contribution from a varying atmospheric 14 C production rate (Supplementary Fig. 4). The relatively large changes in atmospheric Δ 14 C simulated here are in line with a deglacial experiment performed with CLIMBER-2 43 showing the dominant role of reduced Southern Ocean stratification in decreasing atmospheric Δ 14 C across HS1. However, the simulated atmospheric Δ 14 C decrease is larger than the one simulated by global carbon cycle box models 44,45 , probably because of the large oceanic Δ 14 C gradient present in the initial LGM state, itself resulting from the weak LGM oceanic circulation. While the simulated changes in atmospheric δ 13 CO 2 and Δ 14 C are in very good agreement with paleo-records, the pCO 2 increase is slightly underestimated over the last 1000 years of HS1. This could be due to an overestimation of the surface alkalinity increase in these simulations (Fig. 5 squares). It is important to note that the global ocean alkalinity inventory was kept constant during these transient simulations. However, a rising sea-level and an increase in deep ocean carbonate ion saturation during HS1 29 could have enhanced carbonate sedimentation, thus leading to a reduced global alkalinity content. Assuming a total glacial-interglacial alkalinity change of 96 μmol L −1 (Methods), an ocean alkalinity decrease of 6 μmol L −1 per 1000 years would induce a steady atmospheric CO 2 rise of~5 ppm per 1000 years (orange line in Fig. 1d).
Recent studies 19,28 raised the possibility that the abrupt CO 2 increase occurring at 16.2 ka could be due to a Northern Hemisphere terrestrial carbon release resulting from a southward shift of the ITCZ. Our model-data comparison shows instead that an increase in Southern Ocean ventilation during HS1 cause an abrupt atmospheric CO 2 rise and δ 13 CO 2 decline in agreement with ice core and marine sediment records. In addition, our simulated terrestrial carbon content increases during this period due to a higher CO 2 content and warmer conditions (Fig. 3g). However, given the short duration of the 16.2 ka event and the relative low resolution of marine sediment records, a terrestrial carbon release cannot be ruled out (Supplementary Fig. 8).
Our simulations suggest that enhanced ventilation of AABW and AAIW played a crucial role in driving the atmospheric CO 2 increase during HS1. Stronger AABW transfers carbon from the deep ocean to the atmosphere, and concurrently leads to a warming south of 30°S on a millennial timescale. In contrast, increased AAIW formation caused by SH westerlies leads to a multi-decadal CO 2 rise.
Rapid ocean carbon release in a global eddy-permitting model. The strength of SH westerly winds exerts a strong control on the slope of Southern Ocean isopycnals and on the strength of the oceanic circulation. Mesoscale eddies are ubiquitous in the Southern Ocean and have a tendency to compensate for many of the wind-driven circulation changes that appear in non-eddying models 46 . To test whether the Southern Ocean CO 2 outgassing response to a strengthening of SH westerlies shown above is a robust feature, we perform an experiment under fixed modernday forcing with a global high-resolution ocean sea-ice carbon cycle model that resolves most of the ocean mesoscale energy (Methods). The model is perturbed with a poleward intensifying Southern Ocean wind scenario based on observed and projected 21st century wind trends (Supplementary Fig. 10). The model's response is dominated by a large polynya in the Weddell Sea, that is similar in scale to an observed 1970's Weddell Sea polynya 47 and lasts for~6 years (Supplementary Fig. 10). The polynya intensifies deep convection in the Weddell Sea, thus leading to a strengthening of AABW from 19 to 30 Sv and a CO 2 outgassing in the Southern Ocean (Fig. 6a, b). Open-ocean convection was likely the dominant source of AABW formation during glacial periods because grounded ice covered most of the Weddell and Ross Seas 48 , where AABW is primarily formed today. Enhanced ventilation of intermediate and bottom waters, and associated increased upwelling of Circumpolar Deep Water lead to a total ocean carbon loss of 42 GtC over 50 years, corresponding to an upper estimate of~20 ppm atmospheric CO 2 increase. Ventilation of Southern Ocean is rapid and leads to large (~−60 μmol L −1 ) DIC anomalies at intermediate depths (Fig. 6c). At 4300 m depth, negative DIC anomalies spread from the Southern Ocean towards the western Pacific and Atlantic basins, reaching about 20°N after 50 years (Fig. 6d). Enhanced chlorofluorocarbon content in bottom and intermediate waters further confirms enhanced AABW and AAIW ventilation (Supplementary Fig. 11).
Similar to the coarse resolution experiment, most of the CO 2 flux is caused by a surface-water DIC increase, due to enhanced Ekman pumping of DIC-rich waters, with a strong compensation effect coming from the associated alkalinity increase (Supplementary Fig. 12a, b). However, warming of surface waters south of Australia and in the South Atlantic ( Supplementary Fig. 12c), due to a poleward shift of the subtropical front, also significantly contributes to the CO 2 flux. In agreement with the simulations performed with LOVECLIM, global export production increases due to enhanced nutrient upwelling. Even though this simulation was performed under fixed modern-day forcing and the carbon cycle response to a latitudinal shift of the SH westerlies could depend on their initial position 24 , this high-resolution simulation supports the significant role played by intensified SH westerlies in driving abrupt Southern Ocean CO 2 outgassing.

Discussion
Reduced NADW formation during HS1 13 weakens the meridional heat transport to the North Atlantic and induces a Southern Ocean warming, but our simulations suggest that the magnitude of this direct heat redistribution is small. In addition, reduced North Atlantic ventilation increases the carbon content in the Atlantic Ocean (Supplementary Table 1). Through oceanic and atmospheric teleconnections, NADW cessation enhances the formation of NPIW, thus ventilating the intermediate North Pacific. However, the decrease in North Pacific oceanic carbon content is compensated by the carbon increase in the North Atlantic, with negligible net effect on atmospheric CO 2 (Fig. 3i,  Supplementary Fig. 5). Other mechanisms impacting SH climate and/or carbon cycle must have therefore been at play.
Modelling studies 34,35 have shown that the North Atlantic cooling can strengthen and shift the SH westerlies poleward via a southward shift of the ITCZ (Fig. 2). During the first phase of weak NADW transport at~17.2 ka 13 , an intensification of the SH westerlies and associated enhanced Southern Ocean deep convection could have led to a rise in atmospheric CO 2 concentration and initiated the deglaciation at high southern latitudes. Our simulations show that the resulting enhanced transport of AABW leads to a transfer of carbon from the deep Pacific to the surface of the Southern Ocean, thus inducing a millenial-scale atmospheric CO 2 increase. In addition, enhanced AAIW formation and the associated increased upwelling of Circumpolar Deep Water, lead to a multi-decadal pCO 2 increase. Both stronger SH westerlies and buoyancy loss near the Antarctic continent act to steepen Southern Ocean isopycnals, thus increasing the southward baroclinic flow. This results in a stronger southward meridional heat transport 39 , decreasing Southern Ocean sea-ice concentration and further warming the mid-to high southern latitudes. Our simulations thus suggest that through its impact on both atmospheric CO 2 and temperature, enhanced Southern Ocean ventilation played a significant role during the last deglaciation.
Changes in SH westerlies, Southern Ocean circulation and the resulting impact on the hydrological cycle and terrestrial biosphere could have reduced the aeolian iron input into the Southern Ocean, thus potentially contributing to the atmospheric CO 2 rise during HS1 5 . However, changes in iron fertilisation alone cannot explain the observed variations in atmospheric Δ 14 C, ocean ventilation and high southern latitude warming. We therefore suggest that changes in iron fertilisation mostly respond to changes in sea-level and Southern Ocean ventilation, possibly providing a positive feedback. In addition, while marine export production decreased north of the polar front 5 , marine sediment cores south of the polar front display an increase in opal flux 9 , thus indicating a limited geographic extent of iron fertilisation changes.
Sediment records from the Iberian margin have shown that 16.2 ka probably corresponds to the beginning of Heinrich event 1 stricto sensu 37 , and to another phase of NADW weakening 36 . North Pacific records suggest that after a period of strong NPIW formation during the first phase of HS1, NPIW weakened significantly at~16.2 ka 33,49 . In addition, Chinese speleothems record a shift to significantly drier conditions at~16.1 ka 50 ( Supplementary Fig. 2). Both NADW and NPIW weakening at 16.2 ka would cool the Northern Hemisphere, thus shifting the ITCZ southward, including in the Pacific sector 51 (Supplementary Fig. 1). This would further strengthen/shift the SH westerlies, enhance Southern Ocean deep convection and lead to the abrupt atmospheric CO 2 increase and δ 13 CO 2 decrease at~16.2 ka. A rapid CO 2 outgassing in the Southern Ocean due to intensified SH westerlies is confirmed by a simulation performed with a global ocean eddy-permitting model (Fig. 6).
Taking smoothing and dating uncertainties (~1 kyr) into account, the above sequence of events is in agreement with radiocarbon records, which suggest increased ventilation of the deep South Atlantic 6 and Pacific 52 and enhanced mixing between deep and intermediate waters at high southern latitudes 38,53 during HS1 (Figs. 1h, i and 4c, d). In addition, enhanced ventilation of deep and intermediate waters is consistent with positive benthic δ 13 C anomalies across HS1 south of 30°S and below 1000 m depth 38,54,55 (Figs. 1g and 4a, b) as well as high opal production in the Antarctic zone 9 .
Our study highlights the crucial role of SH westerlies in driving abrupt atmospheric CO 2 rise and associated global climate changes. Given the projected poleward intensification of SH westerlies over the 21st Century, and the fact that the Southern Ocean has absorbed~10% of anthropogenic CO 2 emissions 56,57 , our results suggest a future reduction in CO 2 sequestration in the Southern Ocean, with significant impacts on future atmospheric CO 2 and climate change.

Methods
Carbon isotopes-enabled earth system model and LGM state. Transient experiments were performed with the carbon isotope-enabled ( 13 C and 14 C) Earth system model LOVECLIM 58 . This coupled model consists of a free-surface primitive equation ocean model (3°× 3°, 20 vertical levels), a  dynamic-thermodynamic sea ice model, an atmospheric model based on quasigeostrophic equations of motion (T21, three vertical levels), a land surface scheme,  a dynamic global vegetation model 59 and a marine carbon cycle model 21,60 . The initial LGM state was selected amongst 28 LGM experiments based on its representation of oceanic δ 13 C and ventilation age distributions 31 . It was obtained by equilibrating LOVECLIM under 35 ka B.P. boundary conditions, namely appropriate orbital parameters 61 , Northern Hemisphere ice-sheet extent, topography and albedo 62 , an atmospheric CO 2 content of 190 ppm, a δ 13 CO 2 of −6.46‰ and Δ 14 C of 393‰. After a 10,000 years long equilibration phase, the model was run transiently until 20 ka with prognostic atmospheric CO 2 , δ 13 CO 2 and Δ 14 C. During the equilibration phase, an atmospheric 14 C production rate of 2.05 atoms cm −2 s −1 was diagnosed and then subsequently applied for all transient simulations, except for the deglacial simulations presented in Supplementary Fig. 4, in which the 14 C production rate is varied according to Hain et al. 44 . This rate is higher than Holocene and present-day 14 C production rate estimates 63 61 and Northern Hemispheric ice-sheet extent, topography and albedo 62 . In addition, meltwater is added to the northern North Atlantic (grey line, Fig. 3a) to obtain a nearly collapsed NADW. To explore the impact of changes in AABW on the global carbon cycle and climate, additional experiments are performed whereby AABW is enhanced by decreasing the buoyancy forcing (LH1-SO, LH1-SO-SHW) in the Southern Ocean and/or enhancing the simulated southern hemispheric westerly windstress (LH1-SHW, LH1-SO-SHW; Fig. 3b, c). The magnitudes of these changes in windstress are within estimates of probable windstress changes. Indeed some PMIP3 models display 10-20% weaker SH westerly winds at the LGM compared to pre-industrial times 65 and SH westerlies have increased by~8% from 1990 to 2010 (equivalent to a~15% increase in windstress) 66 and are forecast to keep on strengthening by~15% over the coming century 67 .
In all transient experiments presented here, the atmospheric 14 C production rate is kept constant and so is the total ocean alkalinity content except for the simulation shown by the orange line in Fig. 1d, where the total ocean alkalinity was decreased at a rate of 6 μmol L −1 per 1000 years starting at 16.2 ka. The magnitude of that decrease is equivalent to a linear global alkalinity decrease of 96 μmol L −1 over a time interval of 16,000 years. This corresponds to a first order approximation based on total alkalinity conservation and on changes in ocean volume (V, linked to sea-level changes): The timeseries of AABW and AAIW transport refer to the maximum stream function of the zonally integrated meridional transport south of 60°S for AABW and at 1225 m between 30 and 60°S for AAIW.
Decomposition of pCO 2 changes. Changes in surface water pCO 2 , which exert a dominant control on atmospheric CO 2 , arise due to changes in DIC, ALK and solubility (SSS and SST) 68 . pCO 2 changes can thus be decomposed as follows (Fig. 5): For DIC, ALK and SSS, ΔpCO 2 can be expressed as where pCO 2Ref is the pCO 2 value at 19 ka; X represents the mean surface DIC, ALK or salinity value and γ DIC , γ ALK , γ SSS are Revelle factors equal to 10, −9.4 and 1, respectively 68 . The temperature contribution is derived from where γ SST is equal to 0.0423 ( 68 ).
Eddy-permitting global ocean sea-ice carbon cycle model. Experiments are conducted with the eddy-permitting global ocean, sea-ice model MOM5, which is based on the Geophysical Fluid Dynamics Laboratory CM2.4 and CM2.5 coupled climate models 69,70 . The model has a 1/4°Mercator horizontal resolution with~11 km grid spacing at 65°S. The MOM5 ocean model has 50 vertical levels and is coupled to the Sea Ice Simulator dynamic/thermodynamic sea-ice model. The atmospheric forcing is derived from version 2 of the Coordinated Ocean-ice Reference Experiments Normal Year Forcing (CORE-NYF) reanalysis data 71,72 . CORE-NYF provides a modern-day climatological mean atmospheric state at 6-h intervals for 1 year and includes synoptic variability. The model is coupled to the Whole Ocean Model with Biogeochemistry and Trophic-dynamics (WOMBAT) model, a Nutrient-Phytoplankton-Zooplankton-Detritus (NPZD) model 73,74 . WOMBAT includes DIC, alkalinity, oxygen, phosphate and iron, which are linked to the phosphate uptake and remineralisation through a constant Redfield ratio. Phytoplankton growth is limited by light, phosphate and iron, with the minimum of these three terms limiting growth. The biogeochemical parameters are slightly modified from the values used in the ACCESS-ESM simulations 75,76 . Two important changes needed for the highresolution simulations were to increase detritus sinking rate to 20 m d −1 and the background iron concentration was set to 0.3 μmol Fe m −3 to reduce nutrient trapping and improve export production in the tropical East Pacific. The formation of calcium carbonate is a constant fraction of organic carbon production. The air-sea exchange of carbon dioxide is a function of wind speed 77 and sea ice concentration. Initial conditions for the biophysical fields are derived from an observation-based climatology 78 .
The model was initialised with modern-day temperature and salinity distributions and equilibrated for 180 years of Normal Year Forcing (CORE-NYF). A control run and a wind perturbation experiment were then run for 50 years. The wind perturbation experiment includes a poleward intensifying wind forcing, namely a 4°southward shift and 15% increase in 10 m wind speeds between 30°S and 65°S ( Supplementary Fig. 10). This wind forcing is based on projected SH wind changes in CMIP5 business as usual scenarios 67 . The model setup and experimental design are similar to the one employed in a previous study 79 , except that previous simulations did not use any neutral physics ocean parameterisations. In the simulations presented here, neutral physics parameterisations are used, based on options from the ACCESS Ocean Model 80 , with Redi diffusivity (600 m 2 s −1 ) and Gent McWilliams skew diffusion (600 m 2 s −1 ). These neutral physics options improved the simulated AABW transport and the distribution of ocean biogeochemical tracers relative to observations, for instance, the oxygen in the Southern Ocean and alkalinity of bottom waters penetrating into ocean basins ( Supplementary Fig. 9).