Rewrite Prevalent North Atlantic Deep Water during the Last Glacial Maximum and Heinrich Stadial 1 as a headline for a science magazine post, using no more than 8 words

Carbon and oxygen stable isotope data from the subpolar North Atlantic and Arctic Mediterranean

Sediment records with carbon and oxygen stable isotope data that span the LGM and HS1 were compiled from around the Arctic Mediterranean and the subpolar North Atlantic for Fig. 1 (Supplementary Table 10). Whereas most of these records are based on epifaunal benthic foraminifera of the Cibicidoides genus, we included three records that are based on infaunal foraminifera species, as indicated in Fig. 1 (sites U1302, MD95-2010 and HH15-1252PC). For infaunal foraminifera, the carbon and oxygen isotopic signatures do not directly reflect the dissolved inorganic carbon (DIC) of local bottom water and we thus applied constant interspecies corrections to account for species-specific fractionation in the oxygen and carbon isotopes. We assumed that isotopes in the shells of Cibicidoides wuellerstorfi are precipitated without biological fractionation and applied constant fractionation factors of −0.64‰ for δ¹⁸O_b (ref. ⁶⁵) and +0.9‰ for the δ¹³C_b values⁶⁶ of Uvigerina peregrina (site U1302)⁶⁷. Cassidulina neoteretis data from site HH15-1252PC (ref. ⁵³) and Cassidulina teretis data from site MD95-2010 (ref. ⁵¹) were assumed to be equally offset from equilibrium and were corrected by −0.64‰ for δ¹⁸O_b values⁵¹ and +1.5‰ for δ¹³C_b. The correction for δ¹³C_b was obtained by alignment of data from site MD95-2010 and Cibicidoides-based data from site PS2644 across the available data overlap between 11.4 to 12.1 ka (20 data points of MD95-2010). This last correction certainly bears the largest potential for bias, but it is noteworthy that the overall interpretations related to Fig. 1 and the definition of the different components of NADW do not exclusively rely on data from non-Cibicidoides species and would not change significantly if only data from Cibicidoides records were considered.

Another potential problem is a low abundance of target foraminifera species such as Cibicidoides, which can promote biases through bioturbation moving these shells vertically in the sediment column. This has been suggested for the glacial section of site PS1243 (ref. ⁴⁷). However, the data agree well with those of site HM52-43 and shallower site PS2644 (which has more abundant glacial Cibicidoides), potentially suggesting a well-preserved palaeoenvironmental signal.

New carbon and oxygen stable isotope time-slice data

Cibicidoides δ¹⁸O_b and δ¹³C_b values from 19 sites across the Atlantic are derived from new records that were produced during the ACCLIMATE project and consistently dated by integrating radiocarbon ages and stratigraphic tie points using the Undatable software² (age models can be downloaded from https://www.seanoe.org/data/00484/59554/). Epifaunal benthic foraminifers of the Cibicidoides genus were hand-picked in the >150 µm size fraction. The C. wuellerstorfi samples were picked when possible.

Cibicidoides oxygen and carbon isotope ratios for these samples were measured at the Laboratory for Climate and Environmental Sciences using a MicroMass IsoPrime100 mass spectrometer on samples of 1–5 specimens using the NBS-19 standard relative to the Vienna PeeDee Belemnite (VPDB). The mean external reproducibility of carbonate standards was 0.05‰ for δ¹⁸O and 0.03‰ (one standard deviation) for δ¹³C; the measured NBS-18 δ¹⁸O was −23.27 ± 0.10‰ VPDB, and δ¹³C was −5.01 ± 0.03‰ VPDB.

Glacial Atlantic water mass proxy database

Proxy data for Atlantic deep water stable isotope signatures (δ¹³C_b and δ¹⁸O_b), [CO₃²⁻] values inferred from B/Ca ratios and radiocarbon ventilation ages (¹⁴C_b–atm)—all measured on benthic foraminifera calcite—as well as radiogenic Nd isotope signatures (εNd), extracted from authigenic sediment phases via dissolution of foraminifera or acid-reductive bulk sediment leaching, were compiled from this study and several original publications and compilations^{14,15,16,18,21,22,31,68,69,70,71,72,73,74} (Extended Data Figs. 4 and 6 and Supplementary Table 13). For carbon and oxygen isotopes and B/Ca ratios, only foraminifera from the Cibicidoides genus, preferably C. wuellerstorfi, were used for this compilation. The Nd isotope data of five sites from the eastern subpolar North Atlantic were omitted because it has been suggested that they are compromised by localized non-conservative effects^75,76,77. The data were averaged for each site for the LGM (23–19 ka bp), HS1 (17.5–14.6 ka bp) and LH (5–0 ka bp) on the basis of the existing age models from the same literature or updated age models within the ACCLIMATE project as described above. For completeness, we briefly describe each proxy in Supplementary Text 1. Carbon and oxygen stable isotope data from non-Cibicidoides foraminifera were not used for this compilation, but only for Fig. 1 as indicated in Fig. 1b.

Ice-volume correction and core-top-offset correction of δ¹⁸O data

All δ¹⁸O_b values reported are corrected for global ice volume (δ¹⁸O_b,ivc) by converting the contemporary sea level relative to modern to a related change in global marine δ¹⁸O. To this end, we used the sea-level curve from ref. ⁷⁸ and assumed a sensitivity of 1.05‰ per 134.3 m sea-level change to enable a comparison across time periods (Extended Data Fig. 1 and Supplementary Figs. 1–3). The accuracy and consistency of this correction depends on the quality of the age models. This is particularly true for time intervals during which global sea level underwent drastic changes, such as the transition between Marine Isotope Stages 2 and 1. In the time intervals around the LGM and HS1, inaccuracies in the age models of, for example, 0.5 and 1.0 ka can lead to biases in the δ¹⁸O_b,ivc signals of up to 0.13 and 0.21‰, respectively (Supplementary Fig. 1). Whereas this source of uncertainty is comparable to measurement and offset correction uncertainties, it is still minor compared with the overall range of the δ¹⁸O_b,ivc that we estimate for the different glacial source waters (2.2–4.0‰).

It has been shown that analyses of oxygen isotope data from foraminifera may suffer from systematic biases due to gas mixing in the mass spectrometer source or other non-ideal instrument performance^31,79. Assuming that such biases can be considered to be constant along a given δ¹⁸O_b record, we corrected the measured signals via a constant site-specific offset that minimizes the difference between the LH (here younger than 4 ka bp if available, otherwise <6 ka or <8 ka) data and the equilibrium Cibicidoides δ¹⁸O_b values computed from local seawater δ¹⁸O and temperature according to equation (9) from ref. ⁵⁴. Local seawater δ¹⁸O, in turn, was inferred from local seawater salinity and basin- and depth-specific linear regressions of seawater δ¹⁸O versus salinity (Supplementary Table 1 and Supplementary Fig. 2). The regressions were generated from the seawater δ¹⁸O dataset of the Goddard Institute for Space Studies^80,81,82. Local seawater temperature and salinity were interpolated from the World Ocean Atlas 2013 gridded global dataset^83,84.

The offsets between LH foraminifera δ¹⁸O_b and equilibrium Cibicidoides δ¹⁸O_b average to a slightly positive value of 0.19 ± 0.56‰ (two standard deviations, n = 104). By adding these constant site-specific offsets, glacial δ¹⁸O_b,ivc values are therefore in average shifted towards slightly higher values and low outliers are reduced (Supplementary Fig. 3). The data thus appear more consistent. In particular, the data spread in water depths between 2 and 4 km is reduced. For example, uncorrected LGM δ¹⁸O_b,ivc data across the Atlantic below a 2 km water depth average to 3.27 ± 0.54‰, and offset-corrected data (δ¹⁸O_b,ivoc) to 3.60 ± 0.49‰. Importantly, the reduced data spread tends to decrease the contribution of low-δ¹⁸O_b (mode 2) source waters in the mixing results.

Multi-proxy mixing calculations

We estimated the relative contributions of different source waters in the deep Atlantic from the multi-proxy dataset using the ‘simmr’ package in the R programming language^43,85 (Extended Data Fig. 7). The simmr package is a Bayesian stable isotope mixing model that uses Gibbs sampling and Markov chain Monte Carlo simulations, and it was originally developed for isotopic mixing calculations in ecological feeding studies but can be directly applied to other mixing scenarios. Starting from an a priori source probability distribution, simmr repeatedly samples the proxy space semi-randomly and tries to find mixing proportions of defined sources that agree with the observation(s). Proxy uncertainties of the sources are included, but not those of individual observations. Prior distributions can be used in the form of suggested source water probability distributions to improve the calculations with additional knowledge of the mixing system. Fixed proxy concentrations can be included and are here used in the form of DIC for δ¹³C and the ¹⁴C ventilation age and Nd concentrations for εNd in the different source waters. The method can cope with an arbitrary number of sources and proxies, but the larger the number of sources compared with proxies the more uncertain the results will be. A posteriori combination of sources can be used to reduce the uncertainty again, which we use, for example, for the estimation of bulk NADW and SSW. The simmr results are given as probability distributions from which we calculated the summary statistics. The choice of sources is critical (see below) and systematically affects the resulting mixing proportions (Supplementary Figs. 11, 14 and 15).

To estimate Atlantic-wide source water volumes we subdivided the Atlantic into three roughly latitudinal sections (Extended Data Fig. 4) and generated depth-smoothed profiles for each time slice and each proxy that was sufficiently available (Extended Data Fig. 6). We solve the simmr model for each latitudinal section and water depth independently, and the results between boxes are only linked via the intrinsic connection in the proxy data. Smoothing the proxy data emphasizes general trends in proxy and source water distributions. However, there is a limited amount of subjectivity related to the choice of smoothing parameters. Thus, we compare the results with an alternative and more objective approach, which we call box-pooling (Extended Data Fig. 7). For this approach, we subdivide the deep Atlantic into seven depth layers (from a water depth of 2 to 5 km in 0.5 km steps, and one box of >5 km) and increase the resolution by splitting each latitudinal section into east and west (Extended Data Fig. 4), resulting in 32 boxes containing observations during the LGM and HS1. This box-pooling approach is more objective because it does not include any other data treatment (such as smoothing or averaging across cores). Both approaches agree well overall (Extended Data Fig. 9), with the most obvious difference being that the box-pooling method leads to more mode-2 NADW in the LGM. We suggest that this is due to the proxy data (especially δ¹⁸O_b) being more variable and including more outliers.

The model ensemble contains 3,000 different parametrizations for each region and water depth, picked semi-randomly from the three time slices and different combinations of modifications to incorporate variations of the model systematics, non-conservative proxy behaviour and the exclusion of individual proxies or source waters (Supplementary Table 6). We regard the final model ensemble as representative of a large range of potential past source water distributions that generally encompasses the limited knowledge about past non-conservative effects, sampling biases, source water properties and transient changes within each time period. From the whole ensemble we select the best 20% of model realizations in each time slice for the source water mixing estimates given in the main text and Figs. 3 and 4. See Extended Data Fig. 9 for a synthesis of the different mixing model results and Supplementary Figs. 4 and 5 for model quality assessment via Taylor diagrams.

Validity of the mixing model

The principal validity of the multi-proxy mixing model was assessed via a direct comparison of its performance in estimating NADW abundance from oceanographic parameters (corrected for remineralization using the Redfield ratios given in ref. ³) with estimates from an optimized multi-parameter analysis (OMPA; Extended Data Fig. 5 and ref. ³). The direct comparison shows the very good agreement of both methods, but it also shows that extreme values (close to 0 or 100% NADW) agree less well, which presumably is rooted in the methodological differences between our Bayesian approach and the OMPA approach⁸⁶. The mean absolute error between both methods for bulk NADW mixing fraction was 5%.

Another assessment of the validity of our approach in using a multi-proxy mixing model comes from the comparison between LH estimates of NADW abundance of 57% (95% range: 52–78%) compared to that from OMPA (75%; Extended Data Fig. 9). This comparison could indicate a bias towards low values in our proxy-based estimates. However, whether this bias is consistent between the Holocene and the two glacial time slices cannot easily be evaluated. Importantly, our LH estimates suffer from the fact that modern source water proxy signatures span only small ranges in proxy spaces, which decreases the quality of the source water mixing estimates (Extended Data Fig. 8). We thus focus our discussion on the comparison between the two glacial time slices, and note that the comparison with the Holocene and the absolute quantification of NADW abundance is less certain. Similarly, estimates of individual NSW abundances are less precise than for bulk NADW due to their similarity in proxy signatures, especially during the LH (Extended Data Fig. 8), whereas the quantification of bulk NADW (that is, the a posteriori combination of its tributary source waters) is considerably more precise.

Fit quality of the mixing model

To evaluate the model ‘goodness of fit’, we rely on Taylor diagrams (Supplementary Figs. 4 and 5). To simplify the choice of best fitting models, we chose the 20% of model runs for each time slice that were closest to the observations in each Taylor diagram. The model errors are generally similar to those of the proxy observations themselves (Supplementary Figs. 6–8; note that in these figures some of the biases to high values in δ¹³C_b, [CO₃²⁻] and εNd are deliberately caused by the model modifications, as described in Supplementary Table 6).

Source waters

The choice of relevant source waters and their characteristics, such as proxy signatures and concentrations, is decisive for the outcome of the proxy mixing model. Here we defined six source waters for the LGM and HS1 Atlantic, in addition to three source waters for the LH and modern (Figs. 1 and 2, Extended Data Fig. 8 and Supplementary Tables 4 and 5). Of the six glacial source waters, three have already been described in considerable detail in the literature (AABW, PDW and u-NADW), although not all relevant proxy signatures have been ascribed (see the main text and, for example, refs. ^{16,18,21,22,34,36,40,68,72,87}). Source water u-NADW-2 has essentially been described in the same literature, although it was not generally considered to be an actual source water or was taken as glacial Antarctic Intermediate Water. Its characteristic deviations in proxy signatures from u-NADW-1 have been mostly explained by some combination of increased carbon remineralization, or meteoric water admixture, and warming, with the potential of additional SSW admixture^34,50,52,87.

We estimated the remaining NSW carbon and oxygen isotope signatures as described in the main text, and the remaining proxy signatures from additional literature^14,35,40,88 and glacial proxy data distributions (Extended Data Figs. 2 and 3). Radiocarbon ventilation ages for u-NADW and AMW are taken from the ‘mid-depth’ and ‘deep’ Atlantic data compilations from ref. ⁸⁸. There is some potential for bias for the definition of some source water proxy signatures, in particular for the radiocarbon ages, [CO₃²⁻] and εNd of the two AMW modes. We translate this into larger uncertainties of our proxy signature estimates, which are fully considered in the mixing calculations^14,76,89 (Fig. 2). AMW εNd signatures are taken from a smoothed estimate of ref. ¹⁴. For the mixing calculations, we generally allow all source waters to contribute, that is, AABW, l-NADW and u-NADW for the modern and LH, and AABW, PDW, l-NADW-1, l-NADW-2, u-NADW-1 and u-NADW-2 for the LGM and HS1, albeit some source waters are removed by some model modifications (Supplementary Table 6 and Supplementary Figs. 14 and 15). Note that AABW and PDW were assigned slightly different proxy signatures for the LGM and HS1, as they are reflected by the defining proxy records in the abyssal South Atlantic and in the deep Pacific, respectively (Extended Data Fig. 2 and Supplementary Table 4). Radiocarbon ventilation ages for AABW and PDW are taken from the ‘deep’ Southern and Pacific Ocean data compilations in ref. ⁸⁸.

Apart from the actual proxy signatures of the source waters, the concentrations with which the proxies are transported in the respective source waters affect the mixing results. This is particularly relevant for Nd, whose concentration in intermediate to deep waters today varies by roughly a factor of two and can presumably be affected by climatically induced changes in continental weathering^45,90. To a much lesser degree it also affects DIC (the relevant concentration for δ¹³C_b and radiocarbon), which varies by roughly 10% in today’s open oceans. In addition, none of these concentrations can currently be directly reconstructed, and they are thus essentially unknown for the considered source waters, although several studies have estimated past DIC (for example, refs. ^22,91,92). We therefore initially assume modern-like concentrations for all source waters. For DIC, we adopted the suggested concentrations for PDW, AABW and u-NADW from ref. ²². Furthermore, we suggest that the concentration of Nd in AMW is the least constrained, as the production of this source water and the weathering regime in its source region were presumably the most different from today⁴⁵. Hence, we incorporated a series of modifications, varying the Nd concentration of AMW or alternatively equalizing the Nd concentration for all source waters (Supplementary Table 6).

The following nomenclature for the volumetric contributions of the different source waters is used in this study:

NADW = u-NADW + l-NADW = NADW-1 + NADW-2

u-NADW = u-NADW-1 + u-NADW-2

l-NADW = l-NADW-1 + l-NADW-2

NADW-1 = u-NADW-1 + l-NADW-1

NADW-2 = u-NADW-2 + l-NADW-2

AMW = AMW-1 + AMW-2

SSW = AABW + PDW

NADW + SSW = 1

See Supplementary Text 2 for a discussion of the connection between glacial AMW and l-NADW⁸.