Fast and slow components of interstadial warming in the North Atlantic during the last glacial

The abrupt nature of warming events recorded in Greenland ice-cores during the last glacial has generated much debate over their underlying mechanisms. Here, we present joint marine and terrestrial analyses from the Portuguese Margin, showing a succession of cold stadials and warm interstadials over the interval 35–57 ka. Heinrich stadials 4 and 5 contain considerable structure, with a short transitional phase leading to an interval of maximum cooling and aridity, followed by slowly increasing sea-surface temperatures and moisture availability. A climate model experiment reproduces the changes in western Iberia during the final part of Heinrich stadial 4 as a result of the gradual recovery of the Atlantic meridional overturning circulation. What emerges is that Greenland ice-core records do not provide a unique template for warming events, which involved the operation of both fast and slow components of the coupled atmosphere–ocean–sea-ice system, producing adjustments over a range of timescales. Interstadial North Atlantic warming during the last glacial period involved the operation of both fast and slow components of the coupled atmosphere–ocean–sea-ice system, according to analyses from the Portuguese Margin and climate model simulations.

D uring Marine Isotopic Stage 3 (MIS 3; 59.5-24 thousand years ago [ka]), multiple high-amplitude air-temperature fluctuations (Dansgaard-Oeschger [DO] stadials and interstadials) have been recognized in Greenland ice-cores and linked to sea-surface temperature (SST) variations and intermittent iceberg discharges in the North Atlantic 1,2 . The structure of this variability in Greenland is characterized by a rapid decadal-scale warming at the onset of an interstadial, a gradual cooling within the interstadial and then an abrupt centennialscale cooling leading to stadial conditions 3 . In North Atlantic deep-sea cores, each stadial section contains ice-rafted detritus (IRD) from circum-North Atlantic ice-sheets, with the most prominent peaks (Heinrich layers) associated with iceberg discharges from the Laurentide ice-sheet through the Hudson Strait [4][5][6] . A distinction has been made between a Heinrich event (HE; the iceberg discharge through Hudson Strait), a Heinrich layer (HL; the physical manifestation of a HE in a marine sediment core) and a Heinrich stadial (HS; the cold period that contains the HE), underlining the fact that HEs are shorter in duration than the HS in which they occur 7,8 . Ocean circulation tracers indicate that the Atlantic meridional overturning circulation (AMOC) weakened during each stadial, with the most extreme AMOC and SST reductions associated with the largest iceberg discharges 9 .
Millennial-scale changes, bearing distinct similarities to DO variability, have been documented in speleothem and pollen archives downstream from the North Atlantic [10][11][12][13] , with intervals corresponding to HS showing the most extreme changes in local climate and vegetation. Atmospheric methane (CH 4 ) concentrations in ice-cores also show rapid variations at 50% of the glacialinterglacial amplitude 14 . The close correspondence with proxy records of low-latitude precipitation has pointed to a significant tropical wetland contribution to changes in atmospheric CH 4 concentration, linked to shifts in the mean latitudinal position of the Intertropical Convergence Zone (ITCZ) and associated tropical rain belts [15][16][17] . High-resolution analyses have also revealed abrupt and sustained CH 4 increases during the second part of HS 1, 2, 4 and 5, attributed to excess CH 4 production from Southern Hemisphere wetlands during extreme southerly displacements of the ITCZ 17 .
Various mechanisms have been proposed to account for the observed millennial-scale changes. Given the evidence for iceberg discharges and ocean changes, a prevalent view has been that freshwater input increased stratification in the North Atlantic and led to a weakening of AMOC with attendant cooling 18 . Meltwater experiments performed with climate models have suggested that the abrupt climate swings could be explained by switches between different AMOC equilibria [19][20][21][22][23] , potentially amplified by a decrease in Laurentide ice-sheet topography 24 . However, one argument for eschewing freshwater forcing has been the simulated slow (centennial-scale) recovery of AMOC and weak temperature response in Greenland, compared to the observed abrupt (decadal-scale) Greenland interstadial transitions 25 . Alternative hypotheses have therefore emphasized rapid changes in the coupled atmosphere-ice-ocean system, with simulated retreat of sea-ice in the Nordic Seas able to cause abrupt increases in Greenland air temperatures similar to proxy evidence 26 . Sea-ice retreat has been induced in coupled high-resolution climate models by increasing the height of the Laurentide ice-sheet 27 , changes in atmospheric CO 2 concentration 28 , or a saltoscillator 29 .
Establishing the rapidity of changes during interstadial warming transitions in different components of the climate system, therefore, could potentially inform on underlying mechanisms. The Portuguese Margin has been a prime location for investigating millennialscale variability, where analyses of palaeoceanographic and terrestrial proxies from the same samples in marine sequences allow an in situ assessment of the relative timing of changes [30][31][32] . This is a direct consequence of the geographic setting of the area, where the combined effects of major river systems and a narrow continental shelf lead to the rapid delivery of terrestrial material, including pollen, to the deep-sea environment. Here, we investigate changes in ocean conditions and hydroclimate in the MIS 3 section of the interval 35-57 ka in core MD01-2444 (refs. [32][33][34][35]. We focus on the sequence of changes during HS4 and HS5 and the transitions into the ensuing interstadials and compare with results from climate model experiments performed with the LOVECLIM Earth System model. The MD01-2444 records show that following peaks in IRD originating from the Hudson Strait during HS4 and HS5, SST and precipitation increased over several centuries, representing a gradual transition into interstadial conditions. The climate model experiment of HS4 reproduces the gradual evolution in SST and precipitation in western Iberia as a result of a slow AMOC recovery, but underestimates the rapidity and magnitude of warming in Greenland, which appears closely related to sea-ice retreat in the Nordic Seas. The results reveal diverse patterns and rates of warming, depending on location and underlying processes.    (Fig. 3). The X-ray fluorescence (XRF) Zr/Sr ratio (Fig. 1a), which is the inverse of Ca/Ti (Fig. 2a), has high values during stadials. The ratios reflect variations in the relative proportion of detrital (Zr, Ti) and biogenic (Sr, Ca) sediment supply, which is controlled by both carbonate production and dilution by terrigenous sediment 35 . During stadials, carbonate productivity decreased relative to the terrigenous supply, with contributions of material from the Tagus river drainage basin 37 , and also from lateral transport, as a result of the deepening and intensification of Mediterranean Outflow Water 38,39 . Input of IRD contributed secondarily to terrigenous sediment during HS4, but very little during HS5. Coincident with the IRD peaks, the abundance of the polar foraminifer Neogloboquadrina pachyderma reaches a maximum in HS5 and especially HS4 (Fig. 1d). The pollen record also shows contractions (expansions) of temperate tree populations (Fig. 1b) during stadials (interstadials), reflecting changes in moisture availability and temperature; semi-desert taxa (Artemisia, Chenopodiaceae/Amaranthaceae) (Fig. 1c) show the largest increases in HS5 and HS4. Following the initial expansion of tree populations, a minor reversal is observed during the early part of interstadials 12 and DO 8 ( Fig. 1b), that does not have any counterpart in the palaeoceanographic proxies.

MIS
Western Iberian climate variability in the context of wider MIS 3 changes. As reported from the Portuguese Margin 30 , MIS 3 changes in planktonic δ 18 O closely map onto those of the Greenland NGRIP δ 18 O ice record 40 (Fig. 4b). This has previously formed the basis of the MD01-2444 age model, whereby abrupt transitions in planktonic δ 18 O were aligned to those of the Greenland δ 18 O ice record 33 . Here, we retain the same principle for our age model, but use the modified WD2014 timescale 41 ("Methods" and Supplementary Fig. 1). The benthic δ 18 O record resembles the West Antarctic Ice Sheet Divide ice-core δ 18 O ice record 42 (Fig. 4a), both in its shape and phasing relative to Greenland and North Atlantic surface temperature records, reflecting changes in local deep-water δ 18 O, a significant portion of which may be due to changes in deep-water sourcing and/or source signature 43 . Accordingly, the asynchronous phasing between the planktonic and benthic δ 18 O has been interpreted as a fingerprint of ocean circulation changes associated with interhemispheric heat transport and the bipolar-seesaw 32 .
Comparison of the temperate tree pollen record with the Hulu speleothem record 44 from SE China (Fig. 4e) reveals a close similarity, suggesting coherent changes in the hydrological cycle. This is in line with a body of evidence from palaeorecords and climate modelling studies showing changes in hydroclimate across southern Eurasia associated with variations in AMOC strength, SST and latitudinal shifts in the ITCZ 12,16,21,22 . It is interesting to note that, as in the MD01-2444 pollen record, a minor reversal within the early part of interstadials 12 and DO 8 is observed in the Hulu δ 18 O record 44 as well as the δ 18 O and δ 13 C speleothem records from Sofular Cave 12 in Turkey. It is unclear whether these represent isolated local events or complex trends of precipitation increase across low-to-mid latitudes.
Variations in SST and shifts in the ITCZ with attendant changes in mid-and low-latitude hydroclimate and tropical wetland extent can also account for the good correspondence between changes in western Iberian vegetation and atmospheric CH 4 concentrations 17 (Fig. 4d). Interstadial expansions of tree populations are coherent with increases in CH 4 concentrations, although the amplitude of the a XRF Ca/Ti record 35 . b Weight percent CaCO 3 (ref. 35 ). c IRD abundance: IRD grains per gram (red) and IRD percent [(IRD grains/(IRD grains + planktonic foraminifera) × 100)] (ref. 33 ). d Bulk carbonate δ 18 O (ref. 35 ). e Planktonic (G. bulloides) foraminiferal δ 18 O (ref. 33 ). Vertical bars denote HS4, HS5 and HS5a. Triangles denote depths from which detrital carbonate grains were taken for stable isotope measurements (See Fig. 3). Fig. 3 Isotopic composition of carbonate grains from Heinrich layers HL4 and HL5. a δ 18 O and δ 13 C measured on single detrital carbonate grains from depths 10.06 m (HL4) in core MD01-2444 (circles) (this study) compared to the δ 18 O and δ 13 C composition of detrital carbonate grains 36 from HL4 in cores HU90-13-028 (triangles), U1303 (squares) and U1308 (diamonds) in the North Atlantic. b δ 18 O and δ 13 C measured on single detrital carbonate grains from depth 11.98 m (HL5) in core MD01-2444 (circles) (this study) compared to the δ 18 O and δ 13 C composition of detrital carbonate grains 36 from HL5 in cores HU90-13-028 (triangles), U1303 (squares) and U1308 (diamonds) in the North Atlantic.
changes is not always similar. Decreases in CH 4 concentration at the start of stadials are paralleled by contractions of temperate tree populations and increases in semi-desert taxa. The highest pollen values of semi-desert taxa occur during the first half of HS5 and HS4 and appear coeval with peaks in dust content (considered to originate from east Asian deserts) in the NGRIP ice-core 45 , pointing to increased aridity across Eurasia (Fig. 4f). The shifts towards higher CH 4 concentrations partway through HS5 and HS4 are temporally close to peaks in IRD and minimum SST in MD01-2444. Given the great distance from Hudson Strait to the Portuguese Margin, it is remarkable that icebergs could travel that far without melting during HS4 and HS5. Sea-ice in the North Atlantic must have been at its maximum extent at these times to permit trans-Atlantic transport 46 . This in turn could have led to the most southerly displacements of the ITCZ, which may account for the excess CH 4 production from Southern Hemisphere wetlands, as suggested in ref. 17 .
Anatomy of HS5 and HS4. While non-Heinrich stadials appear as relatively brief intervals of uniform cooling, HS5 and HS4 contain considerable structure. We therefore focus on the sequence of changes in different proxies during these intervals within MD01-2444 (Fig. 5). While the absolute ages depend on the alignment to the Greenland ice-core record, the relative phasing of our proxies within the same sediment sequence is independent of the age model. We demarcate the onset and end of HS5 and HS4 by using the mid-point of the most abrupt transitions in planktonic δ 18 O, which may reflect rapid temperature changes in the optimum habitat of G. bulloides. Similar boundaries emerge when we employ the statistical tools Rampfit/Breakfit to obtain an objective b NGRIP ice-core δ 18 O ice (ref. 40 ) and MD01-2444 δ 18 O planktonic . c MD01-2444 percent IRD and N. pachyderma (inverse Y-scales). d Atmospheric methane concentration in WDC (ref. 17 ). e Hulu speleothem δ 18 O calcite (ref. 44 ) and MD01-2444 temperate tree pollen percentages. f NGRIP dust content 45 and MD01-2444 semi-desert (Artemisia and Chenopodiaceae/Amaranthaceae) pollen percentages. All records are plotted on WD2014 timescale 41 , except Hulu, which is on its own U/Th chronology 44 . Vertical bars denote HS4 and HS5.
determination of the abrupt transition mid-points in the planktonic δ 18 O and the XRF Zr/Sr data ("Methods" and Supplementary Table 1). Overall, the proxies show a relatively fast cooling phase (200-300 years [Supplementary Table 2]) at the start of the stadials, leading to an interval (~750-900 years) characterized by maximum IRD and N. pachyderma values, minimum alkenone SST and maximum semi-desert pollen percentages, representing the most extreme impact of the iceberg discharge. This is followed by a 500-800-year interval (Supplementary Table 2) of slowly evolving conditions with superimposed oscillations towards more biogenic sediments, reduced IRD and N. pachyderma abundances, higher alkenone SST, increasing temperate tree populations. The presence of these gradual changes alongside the more abrupt shifts in planktonic δ 18 O and high sediment accumulation rates (20-24 cm kyr -1 ) suggest that these features are not artefacts of bioturbation. In addition, N. pachyderma abundance mirrors the gradual changes in alkenone SST, despite belonging to different size fractions, which suggests that particle-size induced differential bioturbation 39 is not significant.
A complex HS4 has also been observed in indicators of lowlatitude climate in ice-cores and Brazilian speleothems, with more extreme conditions occurring near the middle of the stadial and considered to represent the southernmost migration of the ITCZ 47,48 . Of particular interest is the gradual evolution towards interstadial conditions in western Iberian SST and precipitation following peak iceberg discharge, in line with Greenland ice-core 17 O-excess and Brazilian speleothem evidence indicating a slow recovery of Northern Hemisphere low-latitude precipitation associated with a gradual northward shift in the ITCZ 47,48 . This differs from atmospheric CH 4 concentrations (Fig. 4d), which resemble a square-wave, with a rapid increase in the middle of HS5 and HS4, followed by quasi-stable values until the next abrupt shift at the onset of the interstadial 17 . Given that atmospheric CH 4 is a globally integrated signal, it is possible that the quasi-stable values during a period that the ITCZ is inferred to have shifted gradually northward represent a balance between the deactivation and activation of different sources in the Southern and Northern Hemispheres 48 . Small warming trends (2.5°C from 47.8 to 47.3 ka; 3°C from 39 to 38.6 ka) are also observed in reconstructed temperatures from the NGRIP ice-core 49 , but are then followed by the abrupt DO warming events (~10-12°C) (Fig. 5c).
Climate model experiment. To explore the rate and phasing of changes occurring during the transition towards interstadial conditions as AMOC begins to recover, we perform a meltwater experiment with the Earth System model of intermediate complexity, LOVECLIM, forced with appropriate MIS 3 boundary conditions ("Methods").
From 40 ka, a triangular meltwater pulse with a maximum amplitude of 0.2 Sv reached after 400 years and decreasing thereafter for 400 years is added to the North Atlantic between 50°N and 60°N (Fig. 6a). This weakens the AMOC, leading to seaice advance in the Nordic Seas and the northern North Atlantic, and a 5.4°C cooling over Greenland (Fig. 6b-e). The reduced meridional heat transport to the North Atlantic induces a 3°C SST decrease at the Portuguese Margin and a 20% reduction in precipitation over western Iberia (Fig. 6f-g), while the ITCZ shifts to the south (Supplementary Fig. 2). To allow the AMOC to recover on its own, the Bering Strait is opened at model year 600 (39.40 ka). The opening of the Bering Strait leads to a reduction of the low salinity anomaly in the North Atlantic due to increased mixing with North Pacific waters, thus slowly reinstating deepocean convection south of the Norwegian Sea between model years 600 and 1080 (39.40 and 38.92 ka) and strengthening the transport of salty and warm tropical waters to the Northeast Atlantic. No significant changes in Greenland temperature are simulated during that time, while SST off the Portuguese Margin gradually increase by~1°C (Supplementary Table 3; Figs. 6c,f and 7a). Once deep convection reaches the Norwegian Sea, the AMOC increases from 10 to 27 Sv over 220 years (model years , the annual mean sea-ice area in the eastern side of the Nordic Seas gradually decreases by 17%, and Greenland temperature increases by 4.6°C (Fig. 6b-d). This is also associated with a 2°C increase in SST on the Portuguese Margin. As temperature increases over the North Atlantic, the ITCZ shifts northward, back to its original position at 40 ka ( Supplementary Fig. 2). Higher North Atlantic SST and a northward ITCZ shift lead to an increase in precipitation in western Iberia and the northern tropics, respectively (Figs. 6g and 7d), in agreement with paleo-records 16 . Precipitation over western Iberia displays strong decadal-scale variability (Fig. 6g). An AMOC overshoot occurs within 30 years (model years 1350-1380; 38.65-38.62 ka), whereby the AMOC increases by~6 Sv, associated with a major loss of sea-ice along the eastern shore of Greenland as well as in the Nordic Seas between~68°N and 78°N and a 2.5°C increase in Greenland temperature (Fig. 6c). This Greenland warming is mostly due to the sea-ice retreat, which coupled with enhanced deep-ocean convection leads to an increase of the ocean-to-atmosphere heat flux of 60 W m −2 . This last warming phase is particularly visible at high latitudes, but still induces a further 0.5°C SST increase on the Portuguese Margin, and enhanced precipitation over the Iberian peninsula and the northern tropics above their 40 ka level (Fig. 7e, f).
The experiment, therefore, simulates an initial slow warming in Greenland in HS4 and then a more abrupt warming, with temperature increases linearly scaling with the extent of sea-ice cover ( Supplementary Fig. 3c), in line with previous studies suggesting that changes in sea-ice and associated ocean-to-atmosphere heat flux directly influence Greenland temperature 25 . The modelled SST rise on the Portuguese Margin also linearly scales with the AMOC changes ( Supplementary  Fig. 3b) in agreement with the gradual SST increase observed in the MD01-2444 alkenone record. However, the simulation suggests an overall temperature increase of~7.1°C in Greenland in 300 years, whereas ice-core records suggest that the abrupt DO 8 warming was at least 10°C and took less than a hundred years 3,49 . If the AMOC recovery is induced by a negative freshwater input instead (Supplementary Methods and Notes), then the main temperature increase over Greenland (7.9°C) occurs in 150 years instead of 300 years ( Supplementary Fig. S4). Thus, both experiments show a centennial-scale gradual AMOC recovery, with the negative freshwater forcing inducing a relatively faster AMOC reinvigoration and a more abrupt warming in Greenland in the final stages. However, both experiments still underestimate the full amplitude and rate of warming in Greenland. The simulated Greenland temperature changes are highly dependent on the rate and amplitude of AMOC change, the location of deep-ocean convection, and the northern North Atlantic/Nordic Seas sea-ice cover. Given the close link between changes in sea-ice cover around Greenland and Greenland temperature 21,50 , this suggests that the present experiments most likely underestimate the stadial-interstadial changes in sea-ice cover. In addition, proxy records from the Nordic Seas highlight a significant sub-surface warming during stadials [51][52][53] , which may be underestimated in our simulations.

Discussion
Our Portuguese Margin records link millennial-scale variability in the high-latitude North Atlantic with lower-latitude hydroclimate changes and reveal that Greenland ice-core temperature records do not provide a unique template for the patterns and rates of warming observed at different locations. The similarity in the trajectory of changes following the maximum HS cooling between our data and modelling results suggests that a multi-centennial AMOC recovery and associated increased oceanic heat transport to the North Atlantic led to a gradual SST increase and a shift to wetter conditions over western Iberia. The AMOC recovery could be due to enhanced advection of salty low-latitude surface waters (a result of reduced precipitation as the ITCZ shifted southwards) to the North Atlantic 54 , potentially amplified by increased atmospheric CO 2 concentration 28 , or by an increase in the height of the Laurentide ice-sheet 27 . Rapid sea-ice retreat, on the other hand, highlights the important role of tipping points leading to large changes observed in Greenland air temperatures within a few decades. Instead of one dominant mechanism, what emerges is a combination of interacting processes that include both fast and slow components of the coupled atmosphere-ocean-sea-ice  10°W-5°W, 36°N-43°N). A 21year smoothing has been applied to the timeseries in (c-g). Bands refer to the phases discussed in the text and in Supplementary Table 3. system and accommodate a diversity of changes on a range of timescales.

Methods
Core location and study interval. Stable isotope analyses. Measurements of foraminiferal stable isotopes were carried out in the Godwin Laboratory for Palaeoclimate Research (University of Cambridge), using VG PRISM and VG SIRA mass spectrometers. On both instruments, samples were reacted sequentially using ISOCARB systems. Results are reported with reference to the international standard VPDB and the precision is better than ±0.08‰ for δ 18 O. For the planktonic isotope record, analyses were made on samples consisting of 30 specimens of G. bulloides in the 300-355 mm size fraction 33 . The benthic isotope record was generated by the late N.J. Bulk oxygen isotopes of carbonate were measured using a ThermoFisher GasBench II equipped with a CTC Combi-Pal autosampler and interfaced via continuous flow with a ThermoFisher MAT 253 isotope ratio mass spectrometer (IRMS) 35 . About 200-300 mg of bulk sample was used. Analytical precision is estimated to be ±0.08% for δ 18 O of an internal standard (Carrara marble). Results are reported in Supplementary Data 1.
Stable isotopes of detrital carbonate were measured on single grains taken from the peaks of HL4 (10.06 m) and HL5 (11.98 m) in MD01-2444. Detrital carbonate grains were picked from the >150 μm fraction and leached in a weak acid. Each detrital carbonate grain was ground and between 20 and 60 μg of powdered sample was loaded into Exetainer vials and sealed with silicone rubber septa using a screw cap. The samples were flushed with CP grade helium, then acidified with 104% orthophosphoric acid, left to react for 1 h at 70 degrees Celsius and then analysed using a Thermo Gasbench preparation system attached to a Thermo Delta V mass spectrometer in continuous flow mode. Each run of samples was accompanied by 10 reference carbonates (Carrara Z) and 2 control samples (Fletton Clay). Carrara Z has been calibrated to VPDB using the international standard NBS19. The results are reported with reference to the international standard VPDB and the precision is better than ±0.08 per mil for δ 13 C and ±0.10 per mil for δ 18 O. Results are reported in Supplementary Data 1.
Planktonic foraminifera fauna. Identification of N. pachyderma in the greater than 150 μm size fraction was carried out on split subsamples containing at least 300 whole specimens 33 . Results are reported in Supplementary Data 1.
Ice-rafted detritus. Lithic counts were undertaken on the same split fraction used for faunal identification 33 . Results are reported in Supplementary Data 1.
Alkenone SST. Sediment samples (2.5-2.57 g) were freeze-dried and extracted by sonication using dichloromethane. The extracts were dried under nitrogen and hydrolysed with 7% potassium hydroxide in methanol. Extraction with hexane yielded a fraction enriched in neutral compounds, which was cleaned with ultrapure water. The purified extracts were derivatised with BSTFA. Samples were injected diluted in toluene. They were then analysed with a Varian gas chromatograph Model 3400 equipped with a septum programmable injector and a flame ionisation detector. Absolute concentration errors were below 10%. Reproducibility tests showed that the uncertainty in the U k' 37 determinations was lower than 0.015 (ca. ±0.5°C). The SST reconstruction (U k' 37 -SST) was obtained from the relative composition of C37 unsaturated alkenones [U k' 37 = C37:2/(C37:2 + C37:3)], using the calibration equation U k' 37 = 0.033*SST + 0.044 (ref. 55 ). Detailed analytical procedures are described in in ref. 34 . Results are reported in Supplementary Data 1.
Pollen analysis. Sediment samples of 1.5-4 cm 3 were prepared for pollen analysis using the standard hot acid digestion technique. Fine sieving, through a mesh of 10 μm or less, was not used as it could result in a loss of pollen. Residues were mounted in silicone oil for microscopic analysis at magnifications of 400, 630 and 1000 times. Nomenclature follows ref. 56 . Abundances are expressed as percentages of the main sum, which includes all pollen except Pinus, Pteridophyte spores and aquatics. Pinus was excluded from the main sum, as it is over-represented in marine sediments 57 . A minimum of 100 pollen grains, excluding Pinus, spores and aquatics, was counted in each sample. Results are reported in Supplementary Data 1.
Pollen studies from continental shelf sequences suggest that transport to these areas is controlled primarily by fluvial and secondarily by aeolian processes 58 . Studies on modern pollen deposition in fluvial systems indicate the rapid incorporation of pollen to marine sediments 58 . On the Portuguese Margin Margin, aeolian pollen transport is limited by the direction of the prevailing offshore winds and pollen is mainly transported to the abyssal site by the sediments carried by the Tagus River 59 . Comparison of modern marine and terrestrial samples along western Iberia has shown that the marine pollen assemblages provide an integrated picture of the regional vegetation of the adjacent continent 59 .
XRF analysis. Archive halves of sections from Core MD01-2444 were analysed using an Avaatech XRF core scanner at the University of Cambridge 35 . The core surface was carefully scraped cleaned and covered with a 4-mm thin SPEXCerti-Prep Ultralene foil to avoid contamination and minimize desiccation. Each section was measured using a current of 0.2 mA at three different voltages: 10, 30 kV using a thin lead filter, and 50 kV using a copper filter. XRF data were collected every 2.5 mm. The length and width of the irradiated surface was 2.5 and 12 mm, respectively, with a count time of 40 s. Results are reported in Supplementary Data 1.
Carbonate content. Bulk sediment was acidified with 10% phosphoric acid using an AutoMateFX carbonate preparation system and evolved CO 2 was measured using a UIC (Coulometrics) 5011 CO 2 coulometer. Analytical precision is estimated to be ±1% by repeated measurement of a carbonate standard 35 . Results are reported in Supplementary Data 1.
Age model. Following refs. 30  Statistical analyses. To infer objectively the timing of the changes across HS4 and HS5 observed in the δ 18 O planktonic and XRF Zr/Sr records, we use two statistical tools, Rampfit 60 and Breakfit 61 . Rampfit is based on the simple assumption that shifts observed in the tracers can be characterized as a ramp, i.e., a linear change from one stable state to the other, or in other words, a continuous function, segmented in three parts, linearly connected between the transition end and start points, T1 and T2. Hence, only time intervals where the timeseries shaped like a ramp are suitable for the use of this tool, e.g., the timing of the end of HS4 and HS5 in the δ 18 O planktonic record and the timing of both the onset and end of HS4 and HS5 in the Zr/Sr record. To define the onset of HS4 and HS5 in δ 18 O planktonic record we use the Breakfit tool. Breakfit also uses a continuous function, but it consists of two linear parts connected at a specific point, T2 which represents the point where the slope of the data series changes. Both methods are based on (1) weighted least-squares regression to determine within a prescribed search time interval the ramp location with Rampfit and the break location with Breakfit and (2) a bootstrap simulation to estimate the uncertainty attached to the results.
Numerical experiment. A numerical experiment of HS4 was performed with the Earth system model of intermediate complexity, LOVECLIM. LOVECLIM includes an ocean general circulation model (3°× 3°, and 20 vertical levels), a dynamicthermodynamic sea-ice model, coupled to a quasi-geostrophic T21 atmospheric model, and a dynamic vegetation model 62 . The model was forced with appropriate MIS 3 boundary conditions: orbital parameters 63 , Northern Hemispheric ice-sheet geometry and albedo 64 , and a constant atmospheric CO 2 concentration of 195 ppm 65 . The global salinity is 1 psu larger than during pre-industrial conditions. The initial conditions at 40 ka are taken from ref. 22 . From 40 ka, a triangular meltwater pulse with a maximum amplitude of 0.2 Sv reached after 400 years and decreasing thereafter for 400 years is added into the North Atlantic between 50°N and 60°N. Given the global sea-level at the time of MIS 3 (ref. 66 ) and the possible sea-level increase during HS4 (ref. 67 ), the Bering strait is opened at year 39.4 ka to allow for the AMOC to recover.

Data availability
All data generated for this study are available in Supplementary materials and also from PANGAEA data publisher (https://doi.pangaea.de/10.1594/PANGAEA.919846). They are also available in the British Oceanographic Data Centre repository, where they are archived with the accession code: UCL200122 and can be requested by contacting BODC enquiries, see: https://www.bodc.ac.uk/about/contact_us/contact_details/. Results of the modelling experiment are available at: https://doi.org/10.26190/5ee170029cf6e. Received: 14 March 2020; Accepted: 18 June 2020;