Articles | Volume 14, issue 9
Research article
18 Sep 2018
Research article |  | 18 Sep 2018

Relative timing of precipitation and ocean circulation changes in the western equatorial Atlantic over the last 45 kyr

Claire Waelbroeck, Sylvain Pichat, Evelyn Böhm, Bryan C. Lougheed, Davide Faranda, Mathieu Vrac, Lise Missiaen, Natalia Vazquez Riveiros, Pierre Burckel, Jörg Lippold, Helge W. Arz, Trond Dokken, François Thil, and Arnaud Dapoigny

Thanks to its optimal location on the northern Brazilian margin, core MD09-3257 records both ocean circulation and atmospheric changes. The latter occur locally in the form of increased rainfall on the adjacent continent during the cold intervals recorded in Greenland ice and northern North Atlantic sediment cores (i.e., Greenland stadials). These rainfall events are recorded in MD09-3257 as peaks in ln(Ti  Ca). New sedimentary Pa  Th data indicate that mid-depth western equatorial water mass transport decreased during all of the Greenland stadials of the last 40 kyr. Using cross-wavelet transforms and spectrogram analysis, we assess the relative phase between the MD09-3257 sedimentary Pa  Th and ln(Ti  Ca) signals. We show that decreased water mass transport between a depth of ∼1300 and 2300 m in the western equatorial Atlantic preceded increased rainfall over the adjacent continent by 120 to 400 yr at Dansgaard–Oeschger (D–O) frequencies, and by 280 to 980 yr at Heinrich-like frequencies.

We suggest that the large lead of ocean circulation changes with respect to changes in tropical South American precipitation at Heinrich-like frequencies is related to the effect of a positive feedback involving iceberg discharges in the North Atlantic. In contrast, the absence of widespread ice rafted detrital layers in North Atlantic cores during D–O stadials supports the hypothesis that a feedback such as this was not triggered in the case of D–O stadials, with circulation slowdowns and subsequent changes remaining more limited during D–O stadials than Heinrich stadials.

1 Introduction

Rapid changes in ocean circulation and climate have been observed in marine sediments and polar ice cores over the last glacial and deglacial period (e.g., Johnsen et al., 1992; Vidal et al., 1997). These observations demonstrate that the ocean's current mode of circulation is not unique and can rapidly switch between dramatically different states, in conjunction with climate changes. Furthermore, these observations highlight the non-linear character of the climate system.

Documenting the precise timing and sequence of events in proxy records is a prerequisite for understanding the processes responsible for rapid climate changes and improving climate models' predictive skills. However, the task is complicated by the difficulty of deriving precise age models for marine sediment cores. When marine cores are radiocarbon dated, uncertainties can arise from bioturbation biases (e.g., Lougheed et al., 2018) and changes in past surface reservoir ages (Waelbroeck et al., 2001; Thornalley et al., 2011). In the best cases, when changes in past surface reservoir ages and bioturbation biases remain limited, dating uncertainties mainly derive from the calibration of radiocarbon ages into calendar ages.

In these cases, errors are less than 150 yr for the time interval spanning 0–11 calendar kyr BP (noted ka), about 400 yr for the 11–30 ka interval, and 600 to 1100 yr for the 30–40 ka interval (Reimer et al., 2013). Thus, minimum relative dating errors between records from different marine sediment cores, or between marine and ice cores records, reach 500 yr at the end of the last deglaciation and increase from 500 to 1500 yr, for increasing ages between 11 and 40 ka. Therefore, it is not possible to quantify leads or lags of less than 500 yr between records from different marine cores, or between marine and ice core records.

Figure 1(a) Map showing the position of the main Brazilian rivers and surface currents that could influence the terrigenous input at the study site; the study site is indicated by a white square. The orange area represents the catchment area of the local rivers directly delivering sediments to the study site, and the green area represents the catchment area of the Sao Franzisco River (Milliman et al., 1975). NBC stands for North Brazil Current, SSEC stands for Southern South Equatorial Current and BC stands for Brazil Current. (b) Salinity section showing the core site and the main water masses in the modern Atlantic Ocean. NADW stands for North Atlantic Deep Water, AAIW stands for Antarctic Intermediate Water and AABW stands for Antarctic Bottom Water.


Here we take advantage of the fact that the northern Brazilian margin core MD09-3257 records both ocean circulation and atmospheric changes. On the one hand, we reconstruct ocean circulation changes based on new sedimentary Pa  Th data and on epifaunal benthic isotopic ratios. On the other hand, sediment Ti  Ca measured by X-ray fluorescence reflects past changes in rainfall on the adjacent continent (Arz et al., 1998; Jaeschke et al., 2007). Because Pa  Th and Ti  Ca are recorded in the same core, their relative phasing can be examined with virtually no relative dating uncertainty.

We first present the new sedimentary Pa  Th data and their relation to changes in mid-depth water transport in the western equatorial Atlantic over the last 45 kyr. We then precisely assess the relative phasing between the changes in rainfall and ocean circulation recorded in core MD09-3257.

2 Material and methods

2.1 Core locations

Core MD09-3257 (0414.7 S, 3621.2 W, 2344 m) was recovered in 2009 from the northern Brazilian margin during the R/V Marion Dufresne cruise MD173/RETRO3 at approximately the same position as core GeoB3910-2 (0414.7 S, 3620.7 W, 2362 m) (Arz et al., 2001; Jaeschke et al., 2007). The improved recovery of deep-sea sediments with little or no deformation of sediment layers was achieved thanks to the systematic use of the CINEMA software (Bourillet et al., 2007; Woerther and Bourillet, 2005). This software computes the amplitude and duration of the aramid cable elastic recoil, as well as the piston displacement throughout the coring phase, accounting for the length of the cable (water depth) and total weight of the coring system.

At present, the northern Brazilian margin is bathed by southward flowing upper North Atlantic Deep Water (NADW) at these depths (Lux et al., 2001; Schott et al., 2003; Rhein et al., 2015) (Fig. 1). The southward advection of dense waters formed at higher northern latitudes is channeled through the western boundary current (Rhein et al., 2015), meaning that our sediment cores are ideally located to detect changes in the transport of northern-sourced waters above a depth of 2500 m.

2.2 X-ray fluorescence spectrometry

Elemental composition was measured by employing nondestructive, profiling X-ray fluorescence (XRF) spectrometry. The measurements were made using an AVAATECH XRF core scanner at the Bjerkness Centre for Climate Research, Bergen (Norway) at intervals of 0.5 mm on core MD09-3257, and using a CORTEX XRF scanner at the Bremen Integrated Ocean Drilling Program core repository at intervals of 0.4 cm on core GeoB3910-2 (Jaeschke et al., 2007). This automated scanning method allows for a rapid qualitative determination of the geochemical composition of the sediment at very high resolution (Croudace and Rothwell, 2015).

2.3 Chronology

Ti  Ca records from core MD09-3257 and core GeoB3910-2 exhibit marked peaks corresponding to increased terrigenous input due to enhanced precipitation and runoff from the continent (Arz et al., 1998; Jaeschke et al., 2007) (Fig. 2). These precipitation events are also recorded in South American speleothems, and have been shown to correspond to North Atlantic cold stadial periods (Cheng et al., 2013). The core GeoB3910-2 radiocarbon (14C) age model shows that increases in sedimentary Ti  Ca are indeed synchronous with decreases in South American speleothem δ18O (Burckel et al., 2015). Based on the observed synchronicity, composite age models of core MD09-3257 and GeoB3910-2 have been developed using 14C dating for the past 35 kyr, combined with the alignment of sediment Ti  Ca increases with decreases in speleothem δ18O for the older portion of the cores (Vazquez Riveiros et al., 2018); thus, speleothem ages beyond 35 ka could be transferred to the marine cores. The chronology of core GeoB3910-2 is based on 17 monospecific radiocarbon dates between 0 and 31 ka (Burckel et al., 2015; Jaeschke et al., 2007). The Ti  Ca record of core GeoB3910-2 was aligned to that of core MD09-3257 in order to transfer the radiocarbon dates of GeoB3910-2 over the interval from 12 to 36 ka to this nearby core. In addition, five monospecific radiocarbon dates over 1–21 ka were directly obtained from core MD09-3257. Speleothem tie points were used to derive the chronology of this core over the period from 38 to 48 ka (Tables S1 and S2 in the Supplement) (Vazquez Riveiros et al., 2018). All radiocarbon dates were converted to calendar dates using the OxCal 4.2 software, the IntCal13 calibration curve (Reimer et al., 2013), and a surface water reservoir age of 550±50 yr over 0–18 ka (Key et al., 2004), and of 750±250 yr over 18–31 ka (Freeman et al., 2016). The final age models of cores GeoB3910-2 and MD09-3257 were obtained using a P_Sequence depositional model (Bronk Ramsey, 2008), i.e., a Bayesian algorithm producing posterior probability distributions for each core depth (Tables S1 and S2) (Vazquez Riveiros et al., 2018).

In the present study, the GeoB3910-2 age scale for the 32–50 ka interval was further adjusted by precise alignment of GeoB3910-2 XRF to the MD09-3257 XRF signal (Fig. S1 in the Supplement), which produced a composite record from these two nearby cores. Given that both XRF signals are virtually identical and measured at very high resolution (sampling step ≤0.5 cm), the mean relative dating uncertainty between the two cores is extremely small (less than 105 yr; Fig. S1).

Here, we use XRF ln(Ti  Ca) rather than Ti  Ca because log-ratios provide a unique measure of sediment composition, in contrast to simple ratios, which are asymmetric (i.e., conclusions based on evaluation of A/B cannot be directly translated into equivalent statements about B/A) and suffer from statistical intractability (Weltje and Tjallingii, 2008). We adopt the same terminology as Burckel et al. (2015) and define the larger ln(Ti  Ca) peaks as precipitation events PE0 to PE5, with PE0 occurring during the Younger Dryas, and PE1 to PE5 occurring during Heinrich stadials 1 to 5 (Fig. 2). What we refer to as Heinrich stadials are strictly those stadials characterized by the occurrence of iceberg discharges in the mid- to high-latitude North Atlantic. We refer to smaller ln(Ti  Ca) peaks corresponding to D–O stadials using the Greenland stadial numbering system, as defined by Rasmussen et al. (2014).

Figure 2MD09-3257 Pa  Th, ln(Ti  Ca) and composite C. wuellerstorfi δ13C (δ13CCw) records versus MD09-3257 age scale, independent from the North Greenland Ice Core Project (NGRIP) age scale. (a) NGRIP air temperature versus the Greenland Ice Core Chronology 2005 (GICC05) age scale (Kindler et al., 2014), transposed from kyr  b2k (before 2000) to ka. Greenland interstadials are numbered according to Rasmussen et al. (2014). To avoid overcrowding the figure, Greenland stadials (GS) and interstadials (GI) are only explicitly named in the case of GS-8 and GI-8. (b) The MD09-3257 core Pa  Th record. Empty symbols denote data points that may be affected by terrigenous fluxes and should be interpreted with caution; crosses denote replicate measurements; the red line connects average values (filled symbols). Pa  Th could not be measured over the first half of PE2 because of the occurrence of two small sand layers (Burckel et al., 2015). (c) The MD09-3257 core and GeoB3910-2 core composite δ13CCw record (Vazquez Riveiros et al., 2018); crosses denote replicate measurements. (d) MD09-3257 ln(Ti  Ca). Diamonds above the x axis indicate calibrated radiocarbon dates in MD09-3257 (filled symbols) and GeoB3910-2 (empty symbols). Triangles indicate alignment tie points to South American speleothem δ18O (Vazquez Riveiros et al., 2018). Grey bands delineate precipitation events recorded in MD09-3257 ln(Ti  Ca).


2.4 Benthic isotopes

Epifaunal benthic foraminifers of the Cibicides wuellerstorfi species were handpicked in the > 150 µm size fraction (Vazquez Riveiros et al., 2018). Core MD09-3257 C. wuellerstorfi 13C 12C (δ13C, expressed in ‰ versus Vienna Pee Dee Belemnite, VPDB) was measured at the LSCE on Finnigan DELTAplus and Elementar Isoprime mass spectrometers on samples of one to three specimens. VPDB is defined with respect to the NBS-19 calcite standard (δ18O =2.20 and δ13C =+1.95 ‰). The mean external reproducibility (1σ) of carbonate standards is ±0.03 ‰ for δ13C; measured NBS-18 δ18O is -23.27±0.10 and δ13C is -5.01±0.03 ‰ VPDB. Core GeoB3910-2 C. wuellerstorfi gδ13C was measured at the University of Bremen, Germany, on a Finnigan MAT 252 mass spectrometer on samples of one to five specimens (Heil, 2006), with a mean external reproducibility (1σ) for carbonate standards of ±0.05 ‰ for δ13C. A composite high-resolution benthic isotopic record was generated by combining isotopic data from the upper 294 cm of core MD09-3257 (covering the last 32 kyr) with isotopic data from the interval from 246 to 451 cm in core GeoB3910-2 for the older part of the record (Vazquez Riveiros et al., 2018).

The 13C/12C isotopic ratio (δ13C, expressed in per mil –  ‰ – versus VPDB) of the epifaunal benthic foraminifer C. wuellerstorfi has been shown to record the δ13C of bottom-water dissolved inorganic carbon (DIC) with minor isotopic fractionation (Duplessy et al., 1984; Zahn et al., 1986; Schmittner et al., 2017). A water mass' initial DIC isotopic concentration is governed by surface productivity in its formation region (i.e., the preferential consumption of 12C by primary productivity, which increases dissolved δ13C), as well as temperature dependent air–sea exchanges (Lynch-Stieglitz et al., 1995). DIC δ13C subsequently decreases as deep water ages, due to the progressive remineralization at depth of relatively 13C-depleted biogenic material. As a result, DIC δ13C largely follows water mass structure and circulation in the modern ocean, and C. wuellerstorfi δ13C (δ13CCw hereafter) has been used to trace water masses and as a proxy of bottom-water ventilation (Duplessy et al., 1988 and numerous subsequent studies). A recent study further highlighted that DIC δ13C more faithfully follows water oxygen content than phosphate content (Eide et al., 2017), lending strong support to the use of δ13CCw as a proxy for bottom-water ventilation – the term ventilation here refers to the transmission of oxygen-rich, atmosphere-equilibrated water to the ocean interior.

2.5 New sedimentary Pa  Th data

New sedimentary (231Paxs,0/230Thxs,0) measurements (excess activity ratio at the time of deposition, Pa  Th hereafter) were produced in core MD09-3257 in order to extend the Pa  Th record of Burckel et al. (2015) and to cover the entire time interval from 10 to 43 ka (Table S3). The excess activity corresponds to the fraction of each radioisotope produced in the water column by uranium (U) decay and is transferred to the sediment by adsorption onto particles sinking in the water column. 230Th and 231Pa excess activities are calculated from bulk sediment measurement by correcting for the contribution of the detrital and authigenic fractions (François et al., 2004; Henderson and Anderson, 2003) using a detrital (238U 232Th) value of 0.5±0.1 (2σ) (Missiaen et al., 2018). These excess activities are then further corrected for radioactive decay since the time of sediment deposition. Bulk sediment measurements were performed by isotopic dilution mass spectrometry on the LSCE MC-ICP-MS (NeptunePlus, Thermo Fisher), following a method derived from Guihou et al. (2010). Error bars (2 standard deviations) on Pa  Th measurements were computed by Monte Carlo runs (Missiaen et al., 2018), accounting for the uncertainties in Pa, Th and U measurements, as well as those of the detrital (238U 232Th) value, spike calibrations and dating.

Sedimentary Pa  Th can be used to reconstruct changes in the renewal rate of water masses overlying the core site. This tracer has been successfully used to reconstruct past changes in deep Atlantic circulation intensity (Burckel et al., 2015 and references therein). 231Pa and 230Th are produced at a constant Pa  Th activity ratio of 0.093 by dissolved uranium, which is homogeneously distributed in oceans. 230Th is much more particle reactive than 231Pa, as reflected by their respective residence time in the ocean (30–40 yr for 230Th and 100–200 yr for 231Pa, François, 2007). Therefore, 230Th is rapidly removed from the water column to the underlying sediment, while 231Pa can be advected by oceanic currents. Thus, when averaged over an entire ocean basin, high (low) flow rates result in high (low) 231Pa export and a low (high) sedimentary Pa  Th ratio. In contrast to δ13CCw, which records the DIC δ13C of bottom waters at the core site, sedimentary Pa  Th does not reflect the flow rate at the seabed but that of a water layer of a few hundreds to more than 1000 m above the seafloor (Thomas et al., 2006).

Several potential caveats of the proxy were tested. In particular, 231Pa has a higher affinity for opal than for the other types of particles (Chase et al., 2002), meaning that high opal fluxes can result in high sedimentary Pa  Th values even in the presence of lateral advection. Similarly, areas of very high vertical particle flux, such as the Atlantic off western Africa, are characterized by high Pa  Th values (Yu et al., 1996; François, 2007; Lippold et al., 2012). Recent studies have shown that the caveats that may apply to this proxy in some areas do not apply to the western tropical Atlantic region. More specifically, a study including core top material from the western tropical Atlantic margin and using a 2-D model (Luo et al., 2010) showed that the measured Pa  Th vertical profile is consistent with a dominant role of the overturning circulation, rather than particle scavenging; this demonstrates that Pa  Th can be used to record changes in water mass overturning rates in that region (Lippold et al., 2011). However, because there are large increases in terrigenous material deposition on the northeastern Brazilian margin during the last glacial, we carefully evaluated/assessed if increased terrigenous deposition may have impacted Pa  Th values.

The 230Th-normalized 232Th flux, hereafter simply referred to as the 232Th flux, is indicative of the vertical terrigenous flux to the core site. As 232Th is a trace element that is mostly contained in the continental crust (Taylor and McLennan, 1985), it is commonly used as a geochemical tracer for material of detrital origin (e.g., Anderson et al., 2006). Besides the main precipitation events PE0 to PE4, there is no significant correlation between the Pa  Th ratio and the 232Th flux (r=0.21, p=0.07) (Figs. S2 and S3). In contrast, because the correlation between Pa  Th and the 232Th flux becomes significant (r=0.57, p≪0.001) when including the main precipitation events (Fig. S3), the high Pa  Th values observed during PE0 to PE4 could be partly caused by increased terrigenous flux and should be interpreted with caution (empty symbols in Fig. 2). Note that a possible terrigenous influence during the main precipitation events does not preclude that the high Pa  Th values during these periods reflect an almost halted oceanic circulation above the core site. Indeed, Pa scavenging by boundary scavenging can be intensified in times of reduced overturning circulation due to boundary scavenging becoming the main control on sedimentary Pa  Th.

Another source of possible biases in Pa  Th results from variations in opal flux (Chase et al., 2002). However, the northern Brazilian margin is known for its low siliceous primary production (Arz et al., 1998). This is confirmed by 230Th-normalized opal flux measurements in MD09-3257, which are below 0.06 g cm−2 kyr−1 (Fig. S3). Moreover, outside of precipitation events PE0 to PE4, there is no correlation between Pa  Th and opal flux (Fig. S3). In conclusion, we may consider that outside of the main precipitation events, our Pa  Th record can be interpreted in terms of changes in the strength of overturning circulation above the MD09-3257 coring site.

2.6 Cross-correlation and wavelet analysis

Assuming that there is a constant phase shift between two time series over their entire length, one can perform a simple cross-correlation analysis and compute how the correlation coefficient between the two time series varies as a function of the time lag between the two series (e.g., Davis, 1986).

We normalized (i.e., subtracted the mean and divided by the standard deviation) and resampled the time series Pa  Th, δ13CCw and ln(Ti  Ca) to a common age scale using scenarios with constant time steps varying between 50 and 500 yr. We then used the R function cor.test (R package stats version 3.2.2) for correlation between paired series (R script in Supplement) to compute the Spearman correlation coefficient between all pairs of the three time series, after having shifted one with respect to the other by increments of the time step.

Another approach consists of classical spectral analysis methods that examine the coherence and phase between two time series in frequency space, such as Fourier transforms. Fourier transforms involve decomposing a signal into infinite-length oscillatory functions (such as sine waves). As such, these methods also rely on the assumption that the decomposition of each signal into characteristic frequencies is valid over its entire length, i.e., that the underlying processes are stationary in time.

In contrast, wavelet analysis can be used to decompose a time series into “time–frequency” space, rather than frequency space, that is, to determine both the dominant modes of variability and how these modes vary in time (Torrence and Compo, 1998). To do so, the wavelet transform decomposes the signal into a sum of small wave functions of finite length that are highly localized in time. Thus, wavelet transform can describe changes in frequencies along the studied time series and are particularly relevant for dealing with climatic signals, since they are in essence not stationary in time, but in constant evolution in response to external forcing (i.e., insolation changes), and as a result of internal climate variability.

Given two times series X and Y, with wavelet transforms WX and WY, the cross-wavelet spectrum is defined as WXY=WXWY*, where WY* is the complex conjugate of WY (Torrence and Compo, 1998). Similarly to Fourier coherency, which is used to identify frequency bands in which two time series are related, the wavelet coherency was developed to identify both frequency bands and time intervals over which the two time series are related. The wavelet coherence between two time series is defined as the square of the smoothed cross-wavelet spectrum normalized by the smoothed individual wavelet power spectra (Torrence and Webster, 1999). This definition resembles that of a traditional correlation coefficient, i.e., wavelet coherence ranges between 0 and 1, and may be viewed as a localized correlation coefficient in time–frequency space (Grinsted et al., 2004).

Analogous to Fourier cross-spectral analysis, the phase difference between two time series can also be computed using a cross-wavelet spectrum. The complex argument arg(WXY) can be interpreted as the local relative phase between X and Y in time–frequency space (Grinsted et al., 2004).

In the present study we use the software developed by Grinsted et al. (2004) to compute the cross-wavelet spectrum, coherence and relative phase between our time series that were normalized and resampled as previously described. To test for the persistence of regions of high cross-wavelet coherence, we ran all cross-wavelet analyses 1000 times for each dataset pair (i.e., a Monte Carlo approach). For each of the 1000 runs, each time data point was randomly sampled, whereby a Gaussian distribution of each data point's value (based on the measurement uncertainty) was used to weight the random sampling. Mean and standard deviation values for the coherence and phase direction were calculated using the 1000 runs.

3 Results

3.1 Ocean circulation proxy records

The Pa  Th record of core MD09-3257 now covers the entire 10–43 ka time interval, encompassing the Younger Dryas (YD) and the last four Heinrich stadials (Fig. 2). We have increased its temporal resolution over the time interval from 31 to 38 ka comprising Dansgaard–Oeschger (D–O) events 5 to 8, with respect to the rest of the study period, in order to examine Atlantic circulation dynamics during D–O events.

Pa  Th data exhibit systematic increases in conjunction with stadials, even if Pa  Th data points that are potentially biased towards elevated values by increased terrigenous input (empty symbols) are discarded. Thus, Pa  Th data indicate that the renewal rate of the water mass overlying the site decreased during stadials. More specifically, transport of the overlying water mass decreased not only during the YD and Heinrich stadials, but also during practically all D–O stadials. Among D–O stadials, the Pa  Th increase is well marked for GS-7, GS-8 and GS-11, but the signal is too noisy to provide a clear picture for GS-6. This noisy Pa  Th signal is very likely due to sediment reworking, given that the δ13CCw record is also noisy over this section of the core. Also, it is noteworthy that no precipitation event is recorded in MD09-3257 or GeoB3910-2 during GS-10 (Fig. 2). There is no clear decrease in the well-dated El Condor (Cheng et al., 2013) speleothem δ18O records associated with GS-10 either, which is in contrast with the other Greenland stadials (Burckel et al., 2015). It would seem that there was no apparent increase in precipitation during GS-10 over tropical South America, in contrast with all other GS of the past 40 kyr. Overall, longer stadials seem to be associated with larger increases in Pa  Th than shorter stadials.

The δ13CCw composite record varies in concert with Pa  Th, with high values indicating the presence of well-ventilated waters during the Holocene and interstadials, and low values indicating a marked reduction in water ventilation during stadials at ∼2350 m in the western equatorial Atlantic (Vazquez Riveiros et al., 2018).

3.2 Relative timing of Pa  Th, δ13CCw and Ti  Ca

Pa  Th, δ13CCw and Ti  Ca are recorded in the same core or in two cores from the same location, which could be precisely aligned through high resolution XRF signals. This situation provides ideal conditions to examine the relative phasing of one proxy with respect to another. Pa  Th and Ti  Ca are recorded in the same core, so their relative phasing can consequently be examined with the smallest possible relative dating uncertainty, whereby the only remaining source of uncertainty is bioturbation. The situation is practically the same when examining δ13CCw versus Ti  Ca or δ13CCw versus Pa  Th. Apart from the unavoidable uncertainty introduced by bioturbation, the relative dating uncertainty between the δ13CCw composite record and any MD09-3257 record is null over 0–32 ka, and amounts to 102 yr on average over the 32–50 ka time interval (Fig. S1).

In what follows, we assess the relative phasing between Pa  Th, δ13CCw and Ti  Ca, using all Pa  Th data points (including Pa  Th values susceptible to being partially impacted by large particle fluxes or boundary scavenging resulting from slower overturning circulation) in order to have sufficient data to examine periodicities ranging from 1000 to 6000 yr. In doing so, we assume that changes in particle fluxes may affect the amplitude of the Pa  Th changes, rather than the timing of these changes. In the following text, we show that excluding the Pa  Th values susceptible to being partially impacted by large particle fluxes does not change our conclusions concerning D–O periodicities (i.e., 1000 to 3000 yr).

Figure 3Spearman correlation coefficient of δ13CCw versus Pa  Th (blue curves), of δ13CCw versus ln(Ti  Ca) (red curves), and of Pa  Th versus ln(Ti  Ca) (green curves), as a function of the time lag. A positive time lag means that series 1 lags series 2 (e.g., δ13CCw lags Pa  Th); a negative time lag means that series 1 leads series 2 (e.g., Pa  Th leads ln(Ti  Ca)). Bold lines correspond to the calculation over the entire time interval from 10.6 to 42.6 ka, thin lines to the calculation over the time interval from 10.6 to 26.6 ka, and thin dashed lines to the calculation over the time interval from 26.6 to 42.6 ka. Vertical dashed lines indicate the time lags corresponding to the maximum correlation coefficients for the three pairs of series over the entire time interval.


3.2.1 Average relative phases

We first apply the simple stationary cross-correlation approach to examine how the correlation coefficients of Pa  Th versus ln(Ti  Ca), of δ13CCw versus ln(Ti  Ca), and of δ13CCw versus Pa  Th, vary as a function of the lag between the different time series (Fig. 3). Prior to computing the correlation coefficients, the three time series were resampled with a time step of 100 yr and normalized.

Taken at face value, these results indicate that Pa  Th leads ln(Ti  Ca) (or Ti  Ca) by 200±100 yr, that there is no significant phase shift between δ13CCw and Ti  Ca, and that δ13CCw lags Pa  Th by 200±100 yr (Table S4). The uncertainty of ±100 yr directly results from the adopted sampling step of 100 yr. In addition, in order to assess the robustness of these results, we applied the same approach to the upper half and lower half of the records. In all cases, we obtained δ13CCw lags over Pa  Th of 200 yr, and Pa  Th leads over ln(Ti  Ca) of 200 or 300 yr, while the phase shift between δ13CCw and Ti  Ca remained between 100 and +100 yr.

Figure 4Cross-wavelet transform of MD09-3257 ln(Ti  Ca) versus Pa  Th. (a, b) Wavelet coherence and phase direction computed using Grinsted et al. (2004) software. The thick contour line corresponds to the 95 % confidence level against red noise. Phase direction is only computed for coherences higher than 0.5. (c, d) Mean coherence and phase direction computed from 1000 Monte Carlo simulations. (e, f) Standard deviation around the mean coherence and phase direction computed from these 1000 Monte Carlo simulations.


Although this simple method has been applied to climatic time series in previous studies (Langehaug et al., 2016; Henry et al., 2016), such results must be interpreted with caution, as the method has been designed for signals that are stationary in time and is therefore not suitable for climatic signals.

3.2.2 Wavelet transforms

The non-stationary character of climatic signals over the last 40–45 kyr is particularly pronounced. Different typical pseudo-periodicities can be identified for Heinrich and D–O stadials. In the case of Heinrich stadials (corresponding to our main precipitation events), the interval from 11.7 to 49 ka comprises five pseudo-cycles that are ∼6 to 9 kyr long (Fig. 2), such that Heinrich stadials over the studied interval are characterized by an average pseudo-periodicity of about 7 kyr. Concerning D–O events, the interval located between HS3 and HS4 (32.5–38.1 ka) comprises three pseudo-cycles that are ∼1.2, 1.5 and 3 kyr long (Fig. 2), yielding an average pseudo-periodicity of about 1.8 kyr.

We computed the cross-wavelet spectrum, coherence and phase between ln(Ti  Ca) and Pa  Th (Fig. 4), between δ13CCw and Pa  Th (Fig. 5), and between δ13CCw and ln(Ti  Ca) (Fig. 6), using the software from Grinsted et al. (2004). The 95 % confidence level against red noise is shown as a thick contour line. Relative phases are only plotted for coherences higher than 0.5 (< 0.5 is masked as dark blue). Note that the shaded areas in Figs. 4–6 correspond to the region of the wavelet transform graphs where the edge effects due to the finite length of the time series limit the ability to carry out cross-wavelet analysis. These regions are not considered in our interpretations.

To assess the robustness of our results, we repeated the cross-wavelet transform for different interpolation resolutions ranging from 50 to 500 yr; therefore, we could verify that the features corresponding to the 95 % confidence level against red noise for a time step of 100 yr are still present at roughly the same time and frequency for other time steps (e.g., see Fig. S4 for results obtained for a time step of 400 yr).

Table 1Relative phases over regions of the cross-wavelet graphs corresponding to coherences > 0.5.

* Within these time intervals, only results from the unshaded region of the wavelet graphs are taken into account.

Download Print Version | Download XLSX

Moreover, we ran a spectrogram analysis in order to confirm our wavelet results and avoid any overinterpretation (see Fig. S5 and explicative caption). Unlike the wavelet, the spectrogram analysis is based on a finite time Fourier transform that spans different periods. Therefore, it provides an alternative base to check wavelet-based results. These tests confirmed the wavelet results for periods between 1 and 6 kyr. Beyond 6 kyr, wavelet results could not be confirmed by spectrograms due to the short duration of the analyzed records. Thus, we do not discuss periodicities longer than 6 kyr in what follows.

With this in mind, the following regions of significant mean coherence and well-defined mean relative phases can be identified in the cross-wavelet graphs between Pa  Th and ln(Ti  Ca) produced by 1000 Monte Carlo runs (Fig. 4, middle panels): a coherence higher than 0.5 is found for periodicities around 2000 yr (ranging from ∼1000 to 3000 yr) over the time interval from ∼28 to 40 ka, and for periodicities around 5000 yr (∼4000 to 6000 yr) over ∼25–40 ka. Computing the average phases over each of these two regions, we find that Pa  Th leads ln(Ti  Ca) by 259±140 yr (1σ) for periodicities of 1000 to 3000 yr over 28–40 ka, and by 631±345 yr (1σ) for periodicities of 4000 to 6000 yr over 15–40 ka (Table 1).

The cross-wavelet graph between δ13CCw and Pa  Th displays slightly different regions of high mean coherence (Fig. 5, middle panels). Examining the same frequency bands as for Pa  Th versus ln(Ti  Ca), we find mean coherences higher than 0.5 for periodicities around 2000 yr over ∼28–40 ka, and for periodicities around 5000 yr over ∼15–40 ka. Average phases for these regions indicate that δ13CCw lags Pa  Th by 279±244 yr (1σ) for periodicities of 1000 to 3000 yr over 28–40 ka, but that the lag of δ13CCw with respect to Pa  Th for periodicities of 4000 to 6000 yr over ∼15–40 ka is not significant (Table 1).

Finally, the regions characterized by mean coherences higher than 0.5 between δ13CCw and ln(Ti  Ca) are similar to those observed in the graph for Pa  Th and ln(Ti  Ca) (Fig. 6, middle panels). However, the average phases between δ13CCw and ln(Ti  Ca) over these regions are not significantly different from zero (Fig. 6d and Table 1), indicating that decreases in δ13CCw are in phase with increases in ln(Ti  Ca) within uncertainties.

Figure 5Cross-wavelet transform of δ13CCw composite record versus MD09-3257 Pa  Th. δ13CCw values have been multiplied by −1 to allow a straightforward reading of the relative phase between a decrease in δ13C and an increase in Pa  Th. (a)(f) as in Fig. 4.


The uncertainties of the leads and lags (Table 1) are computed assuming Gaussian error propagation of the two following independent uncertainties: (i) the standard deviation of the mean relative phases over the given time–frequency region (Figs. 4–6d), and (ii) the median value of the standard deviation computed by 1000 Monte Carlo runs over the same time–frequency region (Figs. 4–6f). In the case of relative phases between Pa  Th or ln(Ti  Ca) and δ13CCw, we also accounted for the additional error due to the combining of the MD09-3257 and GeoB3910-2 δ13CCw records.

Finally, we applied the aforementioned cross-wavelet method to the subset of Pa  Th data points not affected by large particle fluxes. For the periodicities between ∼1000 and 4000 yr, the results obtained using this subset (Fig. S6) are virtually unchanged with respect to the results obtained using the entire dataset. For longer periodicities, coherence decreases as expected because the suppressed data points are all located in the main precipitation events (i.e., the YD and Heinrich stadials).

Figure 6Cross-wavelet transform of δ13CCw composite record versus MD09-3257 ln(Ti  Ca). δ13CCw values have been multiplied by −1 to allow a straightforward reading of the relative phase between a decrease in δ13C and an increase in ln(Ti  Ca). (a–f) as in Figs. 4 and 5.


4 Discussion

4.1 Reconstructed ocean circulation changes over the last 45 kyr

Oceanographic studies have shown that the southward transport of northern-sourced waters in the equatorial Atlantic mainly takes place between a depth of ∼1300 and 4000 m in a ∼100 km wide Deep Western Boundary Current (DWBC) (Lux et al., 2001; Rhein et al., 2015). Using hydrographic, geochemical and direct velocity measurements acquired in 1993 to inverse an ocean circulation model, Lux et al. (2001) estimated that the volumetric flow of upper NADW occupying water depths between ∼1300 and 2300 m at 4.5 S within the DWBC is 11.2 Sv (1 Sv = 106 m3 s−1). This estimate is in good agreement with the 10.9 Sv estimated by Schott et al. (2003) based on data from 13 shipboard current-profiling sections taken during the World Ocean Circulation Experiment period (1990–2002).

Our data show that outside of the main precipitation events, the total vertical particle flux did not vary much (remaining within 25.1±3.6 g m−2 yr−1, 1σ) (Fig. S2). The Pa  Th values of these interstadials are similar or slightly higher than those of the late Holocene (Fig. 2), suggesting that the transport of the water mass overlying the MD09-3257 core site was also ∼10 Sv during these interstadials.

It is more difficult to translate the observed increases in Pa  Th during stadials into quantified decreases in water mass transport. However, our new data bring additional observational constraints on the Atlantic circulation changes associated with last glacial millennial climate changes. Two recent studies have indicated that decreases in northern-sourced deep water flow took place during each stadial. On the one hand, increases in Pa  Th during each stadial of the last glacial have been observed at a very deep western North Atlantic site located at ∼42 N and a depth of 4500 m (Henry et al., 2016). On the other hand, reconstructions of water corrosiveness in a South Atlantic core located at ∼44 S and a depth of 3800 m indicate the absence of northern-sourced deep water at that site during stadials, whereas nearly all interstadials of the last 60 kyr are characterized by incursions of northern-sourced deep water into the deep South Atlantic (Gottschalk et al., 2015). Together with these independent results, our results indicate that decreases in both the flow rate and extension of northern-sourced deep waters during stadials were not limited to very dense waters circulating at 3800 m or deeper, but also affected water mass transport above 2350 m in the western equatorial Atlantic.

4.2 Relative timing of Pa  Th, δ13CCw and Ti  Ca

4.2.1 Stationary cross-correlation versus cross-wavelet results

At the MD09-3257 site, cross-wavelet graphs (Figs. 4–6) show that significant coherence and well-defined relative phases between δ13CCw, Pa  Th and Ti  Ca can only be identified in some regions of the time–frequency space. For instance, when examining the relative phase between δ13CCw and Pa  Th over the interval from 10 to 43 ka, a meaningful relative phase can only be identified over ∼28–40 ka at D–O frequencies (i.e., periodicities of 1000 to 3000 yr) (Fig. 5, Table 1). Furthermore, cross-wavelet results indicate that δ13CCw lags Pa  Th by 279±244 yr at D–O frequencies, and that decreases in δ13CCw are in phase with increases in Pa  Th for periodicities of 4000 to 6000 yr (i.e., closer to Heinrich periodicities). This is in contrast with the constant 200±100 lag of δ13CCw with respect to Pa  Th obtained by cross-correlation between the two same time series (Fig. 3), and confirms that the latter method yields imprecise and unreliable results when applied to non-stationary climatic signals.

Nevertheless, cross-correlation has recently been applied to climatic signals (Langehaug et al., 2016), including Pa  Th and δ13CCw records from the last glacial (Henry et al., 2016). In the latter study, cross-correlation between marine records from two deep Bermuda Rise cores and the NGRIP ice oxygen isotopic record was used to infer that deep Bermuda Rise δ13CCw led NGRIP by approximately two centuries, and that Pa  Th was approximately in phase with NGRIP over the interval from 25 to 60 ka (Henry et al., 2016). The authors further inferred that Pa  Th lags δ13CCw by two centuries at their deep Bermuda Rise site. However, as shown here, cross-correlation is not a suitable method to analyze non-stationary climatic signals such as those of the last glacial. Moreover, the inferred relative phases between the marine and NGRIP records are much smaller than the dating error for each individual time series; therefore they are also much smaller than the relative dating error of one time series with respect to the other. In summary, the application of stationary cross-correlation techniques and incomplete consideration of geochronological uncertainty casts doubt on the conclusions of the aforementioned studies.

4.2.2 Lead of Pa  Th with respect to ln(Ti  Ca)

Our cross-wavelet results show that MD09-3257 Pa  Th leads ln(Ti  Ca) by 259±140 yr (1σ) for periods of 1000 to 3000 yr during the ∼28–40 ka time interval, and by 631±345 yr (1σ) for periods of 4000 to 6000 yr during ∼15–40 ka (Table 1). Periods of 1000 to 3000 yr correspond to pseudo-periodicities typical of D–O stadials, while periods of 4000 to 6000 yr are close to those of Heinrich stadials. It can be noted that the cross-wavelet results for D–O periodicities are only significant for the ∼28–40 ka time interval, which indeed corresponds to the interval of our records for which D–O events are best recorded.

It is important to examine if the observed relative phases could be an artifact due to bioturbation. It has been shown that smaller particles are more likely to be transported by bioturbation than larger particles (Wheatcroft, 1992; McCave, 1988; DeMaster and Cochran, 1982), and that this results in fine particles having apparent younger ages than coarse particles from the same depth in a core (Brown et al., 2001; Sepulcre et al., 2017).

Sedimentary Pa  Th is measured on bulk sediment samples, with dissolved Pa and Th being more readily adsorbed on small particles because of their higher surface to volume ratio (Chase et al., 2002). It has been shown that 50 %–90 % of 230Th excess inventory is found in particles smaller than 10 µm (Kretschmer et al., 2010; Scholten et al., 1994; Thomson et al., 1993). Therefore, it is reasonable to assume that the Pa  Th signal is mostly carried by small particles (< 100 µm).

Assessing the size fraction corresponding to the Ti  Ca signal is more complicated. XRF measurements show that the marked changes in ln(Ti  Ca) recorded in MD09-3257 result from sharp changes in both Ca and Ti concentration in the sediment. Ca is a component of marine calcite and aragonite, and is thus mainly carried by large size fractions of the sediment (> 60 µm). However, previous studies have shown that changes in marine carbonate production and dissolution between 2000 and 3000 m in the western tropical Atlantic were relatively small over the last glacial (Rühlemann et al., 1996; Gerhardt et al., 2000). Therefore, the sharp decreases in Ca concentration during stadials result from the dilution of marine carbonates by the increased input of terrigenous material. Therefore, the ln(Ti  Ca) is driven by changes in terrigenous input, rather than changes in marine carbonate production or dissolution. It is difficult to assess in which particle size fraction Ti is mostly concentrated. Knowing that Rb and K are typical constituents of clays, and thus characteristic of small grain sizes, we verified if a phase shift could be detected between the XRF Ti signal and the XRF Rb and K signals. We found no relative offset between Ti and Rb and almost no relative offset between Ti and K, with the inflexion point in the K signal taking place 0.05 cm deeper than in the Ti signal. Given that core MD09-3257 sedimentation rates range from 6 to 14 cm kyr−1, 0.05 cm corresponds to 5 to 10 yr; thus, it is completely negligible with respect to the observed phase shifts between ln(Ti  Ca) and Pa  Th. Therefore, we may consider that ln(Ti  Ca) and Pa  Th are both carried by small particles and that the observed phase shifts between these two signals are not the result of bioturbation.

Finally, if, against all likelihood, bioturbation were responsible for a lead of Pa  Th with respect to ln(Ti  Ca), such a lead would be independent of the examined periodicity. Therefore, we may reasonably assume that the observed lead of Pa  Th with respect to ln (Ti  Ca) is not an artifact resulting from bioturbation.

We compute a 631±345 yr (1σ) lead for Pa  Th over ln(Ti  Ca) by cross-wavelet analysis for frequencies close to those characterizing Heinrich stadials. This lead is comparable to the relative phase previously estimated between MD09-3257 Pa  Th and Ti  Ca at the onset of HS4 (690±180 yr) and HS2 (1420±250 yr), respectively, based on the identification of the transition in the Pa  Th and Ti  Ca signals at the beginning of these two stadials (Burckel et al., 2015). The large lead of Pa  Th with respect to ln(Ti  Ca) is clearly visible for the YD and all Heinrich stadials, except HS1 (Fig. 2). The apparent synchronicity of the Pa  Th and ln(Ti  Ca) signals at the onset of HS1 in core MD09-3257, as also recently observed in another core from the northern Brazilian margin (Mulitza et al., 2017), suggests that the sequence of events was different at the beginning of HS1 from those at the beginning of the YD and other Heinrich stadials. Such a different sequence of events seems to indicate that the increase in rainfall over tropical South America during HS1 was not a response to a decrease in Atlantic overturning circulation. Instead, a southward shift of the low-latitude atmospheric convection zone (Intertropical Convergence Zone, ITCZ), along with its associated maximum in precipitation, could have occurred in response to extended Northern Hemisphere ice sheets and sea ice cover without any change in ocean circulation (Chiang et al., 2003). This atmospheric mechanism would have prevailed at the beginning of HS1 because ice sheets reached their maximum extent around that time.

Our results also indicate that a significant lead of Pa  Th with respect to ln(Ti  Ca) is present at D–O frequencies. Moreover, the lead of Pa  Th with respect to ln(Ti  Ca) is markedly shorter at D–O frequencies (259±140 yr) than at Heinrich frequencies (631±345 yr).

Climate models simulate a southward shift of the ITCZ in response to a slowdown of the Atlantic meridional overturning circulation (AMOC), but after just a few years (Dong and Sutton, 2002). In contrast, our results indicate that rainfall increases in the region adjacent to MD09-3257 occurred several hundred years after the increase in sedimentary Pa  Th at our core site. Furthermore, this lead of sedimentary Pa  Th over ln(Ti  Ca) should be taken as a minimum lead of AMOC over ln(Ti  Ca) because a change in AMOC does not instantaneously translate into a change in sedimentary Pa  Th. A delay between a change in AMOC and the resulting change in sedimentary Pa  Th is indeed expected, which depends on the propagation time of the circulation change to the core site and on the response time of dissolved Th and Pa in the water column overlying the core site (i.e., 30–40 for 230Th and 100–200 yr for 231Pa, François, 2007). However, increases or decreases in sedimentary Pa  Th should be measurable before the dissolved Th and Pa have fully adjusted to the new circulation regime, especially at sites with high sedimentation rates such as our study site. Thus, we expect this additional delay to be less than 100 yr and much smaller than the computed lead of MD09-3257 sedimentary Pa  Th over ln(Ti  Ca).

A mechanism has been proposed by Burckel et al. (2015) to explain the large lead of AMOC slowdowns during Heinrich and D–O stadials with respect to precipitation events over tropical South America. In this scenario, AMOC slowdowns are progressively amplified through a positive feedback linking the decrease in deep water formation to subsurface warming at high northern latitudes (Mignot et al., 2007; Alvarez-Solas et al., 2013), leading in the case of Heinrich stadials, to erosion of ice shelves and iceberg discharges, which in turn reinforce the initial AMOC slowdown. In contrast, AMOC slowdowns associated with D–O stadials would not trigger such a positive feedback loop and would consequently remain limited.

Alternatively, or in addition to an actual lead of the changes in AMOC with respect to precipitation events over tropical South America, another factor could induce a lead of Pa  Th with respect to ln(Ti  Ca) in core MD09-3257. It has been shown that the North Brazil Current (NBC) is able to transport terrigenous material laterally (Allison et al., 2000). Also, different studies have shown that a weakening of the AMOC is associated with a decrease of NBC transport, taking place not only on decadal timescales (Zhang et al., 2011), but also during the YD and HS1 (Arz et al., 1999; Wilson et al., 2011). Based on this evidence, a recent study suggested that a reduced NBC during HS1 allowed the enhanced input of terrigenous material to settle on the continental margin offshore of northeastern Brazil, instead of being transported northward (Zhang et al., 2015). Thus, it seems possible that terrestrial input would be deviated northward as long as the NBC was vigorous and reached the core site only once the NBC and AMOC were sufficiently reduced, thereby yielding a time-delayed peak in ln(Ti  Ca). If this were the case, the lag of the terrestrial input signal with respect to the Pa  Th signal would be partially or totally caused by the impact of the NBC on terrigenous material deposition (Zhang et al., 2015). Therefore, the exceptional synchronicity of the onset of terrigenous influx and AMOC slowdown at the beginning of HS1 could be due to the exceptionally large fluxes of large grain size material eroded from the proximal exposed shelf during low eustatic sea level, which would have rained down through the water column, even before full reduction of the NBC.

However, in the absence of direct measurements of the NBC velocity and vertical particle flux on the northeastern Brazilian margin, the actual delay of terrestrial input with respect to NBC slowdown remains speculative.

4.2.3 Lag of δ13CCw with respect to Pa  Th

Our cross-wavelet results show that δ13CCw lags Pa  Th by 279±244 yr (1σ) at D–O frequencies over ∼28–40 ka at the MD09-3257 site (Table 1).

δ13CCw is measured using > 150 µm foraminifera; thus, it is carried by much larger particles than the Pa  Th signal. Therefore, differential bioturbation mixing processes would lead to Pa  Th being carried by sediment material younger than the epibenthic foraminifera sampled within the same depth interval. Thus, bioturbation may induce an artificial lead of Pa  Th with respect to δ13CCw. Knowing that the sedimentation rates of MD09-3257 and GeoB3910-2 vary between 6 and 14 cm over the interval from 28 to 40 ka, a 280 yr lead translates to a downward shift of 2 to 4 cm in the sediment column, which seems plausible for the effect of differential bioturbation.

In conclusion, the lag of δ13CCw with respect to Pa  Th at D–O frequencies during ∼28–40 ka is likely an artifact resulting from the differential bioturbation of fine and coarse particles. The same differential bioturbation processes likely also affect the relative phase between δ13CCw and ln (Ti  Ca). Thus, we will not discuss the results of the cross-wavelet analyses involving δ13CCw any further.

5 Conclusions

New sedimentary Pa  Th data from core MD09-3257 located on the northern Brazilian margin (∼4 S, 36 W) at a depth of ∼2350 m indicate decreases in water mass transport above the core site during all Greenland stadials of the last 45 kyr. Together with two other recent studies (Gottschalk et al., 2015; Henry et al., 2016), these results demonstrate that all stadials of the last 45 kyr were not only characterized by decreases in flow rate and extension of northern-sourced waters below a depth of 3800 m, but also by decreases in mid-depth water mass transport in the western equatorial Atlantic.

Due to its exceptional location, core MD09-3257 records both ocean circulation and atmospheric changes. Ocean circulation changes induce changes in sedimentary Pa  Th and δ13CCw, whereas changes in precipitation over the adjacent continent induce changes in marine sediments Ti  Ca.

Using cross-wavelet transforms and spectrogram analysis, we were able to precisely and robustly assess the relative phase between MD09-3257 sedimentary Pa  Th and ln(Ti  Ca) signals over the interval from 10 to 43 ka with minimal uncertainty. This is owing to the fact that both signals are recorded in the same sediment core. We show that Pa  Th leads ln(Ti  Ca) by 259±140 yr (1σ) at D–O frequencies over 28–40 ka, and by 631±345 yr (1σ) for periodicities close to Heinrich periodicities (4000 to 6000 yr) over 15–40 ka.

In other words, our cross-wavelet transforms and spectrogram analysis results show that changes in water mass transport between a depth of ∼1300 and 2300 m in the western equatorial Atlantic (i.e., within a ∼1000 m water layer above MD09-3257 core site) preceded changes in precipitation over the adjacent continent by 110 to 400 yr at D–O frequencies, and by 280 to 980 yr at Heinrich-like frequencies.

We suggest that the large lead of ocean circulation changes with respect to tropical South American precipitation changes at Heinrich-like and D–O frequencies is likely related to the action of a positive feedback in the case of Heinrich stadials, in agreement with Burckel et al. (2015). In that case, an AMOC slowdown would lead to subsurface warming at high northern latitudes, inducing ice-sheet calving and iceberg discharges that would in turn reinforce the initial AMOC slowdown. In contrast, the absence of marked ice rafted detritus layers in North Atlantic sediments during D–O stadials suggests that in the case of D–O stadials, AMOC slowdowns did not trigger such a positive feedback and, consequently, remained limited (Burckel et al., 2015).

Finally, the relative lead of Pa  Th over ln(Ti  Ca) is visible for the YD and for all Heinrich stadials, except HS1. In the case of HS1, the southward shift of the ITCZ may have been an atmospheric response to the maximum extent in northern high-latitude ice sheets and sea ice cover (Chiang et al., 2003) around that time, rather than a progressive response to a slowdown of the AMOC, as is the case for the other stadials. These different atmospheric and oceanic scenarios remain to be tested by numerical experiments performed over several thousands of years in glacial conditions, whereby climate models compute water and calcite δ18O, DIC δ13C and sedimentary Pa  Th.

Data availability

Data related to this article are available as a Supplement file and on Pangaea.


The supplement related to this article is available online at:

Author contributions

CW and SP designed the research. EB, PB, JL, FT and AD performed the sedimentary Pa  Th measurements. BCL performed the wavelet analyses. DF and LV contributed expert advice on statistical results and performed the spectrogram analyses. LM produced the sedimentary Pa  Th values and error bars from MC-ICP-MS output. NVR improved the age models of the two cores. CW, NVR and TD participated in the 2009 RETRO coring cruise. HWA contributed expert knowledge on the Brazilian margin. CW and TD obtained funding. CW and BL wrote the manuscript.

Competing interests

The authors declare that they have no conflict of interest.


This is a contribution to the ACCLIMATE ERC project; the research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 339108. Core MD09-3257 was collected on board R/V Marion Dufresne during the 2009 RETRO coring cruise, supported by IPEV, ANR project ANR-09-BLAN-0347 and the ESF EUROMARC project “RETRO”. We thank the IPEV team, crew members of R/V Marion Dufresne and all scientists who participated in the 2009 RETRO cruise. We also thank Matthieu Roy-Barman for advice regarding Pa  Th measurements from the LSCE MC-ICP-MS. We acknowledge Vincent Scao and Jørund Strømsøe for XRF measurements, Christophe Moreau, Jean-Pascal Dumoulin, and the UMS ARTEMIS for AMS 14C dates, as well as Gülay Isguder, Lucile Mauclair and Fabien Dewilde for invaluable technical assistance. We are grateful to Roger François and one anonymous reviewer for their helpful comments on an earlier version of this article. This paper is LSCE contribution 6408.

Edited by: Luc Beaufort
Reviewed by: Roger Francois and one anonymous referee


Allison, M. A., Lee, M. T., Ogston, A. S., and Aller, R. C.: Origin of Amazon mudbanks along the northeastern coast of South America, Mar. Geol., 163, 241–256, 2000. 

Alvarez-Solas, J., Robinson, A., Montoya, M., and Ritz, C.: Iceberg discharges of the last glacial period driven by oceanic circulation changes, P. Natl. Acad. Sci. USA, 110, 16350–16354,, 2013. 

Anderson, R., Fleisher, M., and Lao, Y.: Glacial–interglacial variability in the delivery of dust to the central equatorial Pacific Ocean, Earth Planet. Sc. Lett., 242, 406–414, 2006. 

Arz, H. W., Pätzold, J., and Wefer, G.: Correlated Millennial-Scale Changes in Surface Hydrography and Terrigenous Sediment Yield Inferred from Last-Glacial Marine Deposits off Northeastern Brazil, Quaternary Res., 50, 157–166, 1998. 

Arz, H. W., Pätzold, J., and Wefer, G.: The deglacial history of the western tropical Atlantic as inferred from high resolution stable isotope records off northeastern Brazil, Earth Planet. Sc. Lett., 167, 105–117, 1999. 

Arz, H. W., Gerhardt, S., Pätzold, J., and Röhl, U.: Millennial-scale changes of surface- and deep-water flow in the western tropical Atlantic linked to Northern Hemisphere high-latitude climate during the Holocene, Geology, 29, 239–242, 2001. 

Bourillet, J.-F., Damy, G., Dussud, L., Sultan, N., Woerther, P., and Migeon, S.: Behaviour Of A Piston Corer From Accelerometers And New Insights On Quality Of The Recovery., Proc. 6th Int. Off shore Site Investig, Geotech. Conf., Confronting New Challenges Shar, Knowledge, 11–13 September 2007, London, UK, available at: (last access: 14 September 2018), 2007. 

Bronk Ramsey, C.: Deposition models for chronological records, Quaternary Sci. Rev., 27, 42–60, 2008. 

Brown, L., Cook, G. T., MacKenzie, A. B., and Thomson, J.: Radiocarbon age profiles and size dependency of mixing in northeast Atlantic sediments, Radiocarbon, 43, 929–937, 2001. 

Burckel, P., Waelbroeck, C., Gherardi, J.-M., Pichat, S., Arz, H., Lippold, J., Dokken, T., and Thil, F.: Atlantic Ocean circulation changes preceded millennial tropical South America rainfall events during the last glacial, Geophys. Res. Lett., 42, 411–418, 2015. 

Chase, Z., Anderson, R. F., Fleisher, M. Q., and Kubik, P. W.: The influence of particle composition and particle flux on scavenging of Th, Pa and Be in the ocean, Earth Planet. Sc. Lett., 204, 215–229, 2002. 

Cheng, H., Sinha, A., Cruz, F. W., Wang, X., Edwards, R. L., d'Horta, F. M., Ribas, C. C., Vuille, M., Stott, L. D., and Auler, A. S.: Climate change patterns in Amazonia and biodiversity, Nat. Commun., 4, 1411–1417,, 2013. 

Chiang, J. C. H., Biasutti, M., and Battisti, D. S.: Sensitivity of the Atlantic Intertropical Convergence Zone to Last Glacial Maximum boundary conditions, Paleoceanography, 18, 1–18,, 2003. 

Croudace, I. W. and Rothwell, R. G. (Eds.): Micro-XRF Studies of Sediment Cores: Applications of a non-destructive tool for the environmental sciences, Dev. in Paleoenviron. Res., Springer, 17, 656 pp., 2015. 

Davis, J. C.: Statistics and data analysis in geology, J. Wiley, New York, 646 pp., 1986. 

DeMaster, D. J. and Cochran, J. K.: Particle mixing rates in deep-sea sediments determined from excess 210Pb and 32Si profiles, Earth Planet. Sc. Lett., 61, 257–271, 1982. 

Dong, B. W. and Sutton, R.: Adjustment of the coupled ocean–atmosphere system to a sudden change in the thermohaline circulation, Geophys. Res. Lett., 29, 18-1–18-4, 2002. 

Duplessy, J.-C., Shackleton, N. J., Matthews, R. K., Prell, W., Ruddiman, W. F., Caralp, M., and Hendy, C. H.: 13C record of benthic foraminifera in the last interglacial ocean: implications for the carbon cycle and the global deep water circulation, Quaternary Res., 21, 225–243, 1984. 

Duplessy, J.-C., Shackleton, N. J., Fairbanks, R. G., Labeyrie, L., Oppo, D., and Kallel, N.: Deepwater source variations during the last climatic cycle and their impact on the global deepwater circulation, Paleoceanography, 3, 343–360, 1988. 

Eide, M., Olsen, A., Ninnemann, U. S., and Johannessen, T.: A global ocean climatology of preindustrial and modern ocean δ13C, Global Biogeochem. Cy., 31, 515–534, 2017. 

François, R.: Chapter Sixteen Paleoflux and Paleocirculation from Sediment 230Th and 231Pa 230Th, Dev. Mar. Geol., 1, 681–716, 2007. s 

François, R., Frank, M., Rutgers van der Loeff, M. M., and Bacon, M. P.: 230Th normalization: An essential tool for interpreting sedimentary fluxes during the late Quaternary, Paleoceanography, 19, 1–16, 2004. 

Freeman, E., Skinner, L. C., Waelbroeck, C., and Hodell, D.: Radiocarbon evidence for enhanced respired carbon storage in the Atlantic at the Last Glacial Maximum, Nat. Commun., 7, 1–8,, 2016. 

Gerhardt, S., Groth, H., Rühlemann, C., and Henrich, R.: Aragonite preservation in late Quaternary sediment cores on the Brazilian Continental Slope: implications for intermediate water circulation, Int. J. Earth Sci., 88, 607–618, 2000. 

Gottschalk, J., Skinner, L. C., Misra, S., Waelbroeck, C., Menviel, L., and Timmermann, A.: Abrupt changes in the southern extent of North Atlantic Deep Water during Dansgaard–Oeschger events, Nat. Geosci., 8, 950–954,, 2015. 

Grinsted, A., Moore, J. C., and Jevrejeva, S.: Application of the cross wavelet transform and wavelet coherence to geophysical time series, Non. Proc. Geophys., 11, 561–566, 2004. 

Guihou, A., Pichat, S., Nave, S., Govin, A., Labeyrie, L., Michel, E., and Waelbroeck, C.: Late slowdown of the Atlantic Meridional Overturning Circulation during the Last Glacial Inception: new constraints from sedimentary (231Pa 2303Th), Earth Planet. Sc. Lett., 289, 520–529, 2010. 

Heil, G.: Abrupt climate shifts in the western tropical to subtropical Atlanic region during the last glacial, PhD Thesis, Bremen University, 121 pp., 2006. 

Henderson, G. M. and Anderson, R. F.: The U-series toolbox for paleoceanography, Rev. Mineral. Geochem., 52, 493–531, 2003. 

Henry, L., McManus, J. F., Curry, W. B., Roberts, N. L., Piotrowski, A. M., and Keigwin, L. D.: North Atlantic ocean circulation and abrupt climate change during the last glaciation, Science, 353, 470–474, 2016. 

Jaeschke, A., Rühlemann, C., Arz, H., Heil, G., and Lohmann, G.: Coupling of millennial-scale changes in sea surface temperature and precipitation off northeastern Brazil with high-latitude climate shifts during the last glacial period, Paleoceanography, 22, PA4206,, 2007. 

Johnsen, S., Clausen, H. B., Dansgaard, W., Fuhrer, K., Gundestrup, N., Hammer, C. U., Iversen, P., Jouzel, J., and Stauffer, B.: Irregular glacial interstadials recorded in a new Greenland ice core, Nature, 359, 311–313, 1992. 

Key, R. M., Kozyr, A., Sabine, C. L., Lee, K., Wanninkhof, R., Bullister, J. L., Feely, R. A., Millero, F. J., Mordy, C., and Peng, T. H.: A global ocean carbon climatology: Results from Global Data Analysis Project (GLODAP), Global Biogeochem. Cy., 18, 1–23, 2004. 

Kindler, P., Guillevic, M., Baumgartner, M., Schwander, J., Landais, A., and Leuenberger, M.: Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core, Clim. Past, 10, 887–902,, 2014. 

Kretschmer, S., Geibert, W., van der Loeff, M. M. R., and Mollenhauer, G.: Grain size effects on 230Th xs inventories in opal-rich and carbonate-rich marine sediments, Earth Planet. Sc. Lett., 294, 131–142, 2010. 

Langehaug, H., Mjell, T., Otterå, O., Eldevik, T., Ninnemann, U., and Kleiven, H.: On the reconstruction of ocean circulation and climate based on the “Gardar Drift”, Paleoceanography, 31, 399–415, 2016. 

Lippold, J., Gherardi, J.-M., and Luo, Y.: Testing the 231PA 230Th paleocirculation proxy: A data versus 2D model comparison, Geophys. Res. Lett., 38, 1–7,, 2011. 

Lippold, J., Mulitza, S., Mollenhauer, G., Weyer, S., Heslop, D., and Christl, M.: Boundary scavenging at the East Atlantic margin does not negate use of 231PA 230Th to trace Atlantic overturning, Earth Planet. Sc. Lett., 333, 317–331, 2012. 

Lougheed, B. C., Metcalfe, B., Ninnemann, U. S., and Wacker, L.: Moving beyond the age–depth model paradigm in deep-sea palaeoclimate archives: dual radiocarbon and stable isotope analysis on single foraminifera, Clim. Past, 14, 515–526,, 2018. 

Luo, Y., Francois, R., and Allen, S. E.: Sediment 231PA 230Th as a recorder of the rate of the Atlantic meridional overturning circulation: insights from a 2-D model, Ocean Sci., 6, 381–400,, 2010. 

Lux, M., Mercier, H., and Arhan, M.: Interhemispheric exchanges of mass and heat in the Atlantic Ocean in January–March 1993, Deep-Sea Res. Pt. I, 48, 605–638, 2001. 

Lynch-Stieglitz, J., Stocker, T., Broecker, W. S., and Fairbanks, R. G.: The influence of air-sea exchange on the isotopic composition of oceanic carbon: Observations and modeling, Global Biogeochem. Cy., 9, 653–665, 1995. 

McCave, I.: Biological pumping upwards of the coarse fraction of deep-sea sediments, J. Sediment. Res., 58, 148–158, 1988. 

Mignot, J., Ganopolski, A., and Levermann, A.: Atlantic subsurface temperatures: Response to a shutdown of the overturning circulation and consequences for its recovery, J. Clim., 20, 4884–4898, 2007. 

Milliman, J. D., Summerhayes, C. P., and Barretto, H. T.: Quaternary sedimentation on the Amazon continental margin: a model, Geol. Soc. Am. Bull., 86, 610–614, 1975. 

Missiaen, L., Pichat, S., Waelbroeck, C., Douville, E., Bordier, L., Dapoigny, A., Thil, F., Foliot, L., and Wacker, L.: Downcore variations of sedimentary detrital (238U 232Th) ratio: implications on the use of 230Thxs and 231Paxs to reconstruct sediment flux and ocean circulation, Geochem. Geophy. Geosy., 19, 1–14,, 2018 

Mulitza, S., Chiessi, C. M., Schefuß, E., Lippold, J., Wichmann, D., Antz, B., Mackensen, A., Paul, A., Prange, M., and Rehfeld, K.: Synchronous and proportional deglacial changes in Atlantic Meridional Overturning and northeast Brazilian precipitation, Paleoceanography, 32, 622–633, 2017. 

Rasmussen, S. O., Bigler, M., Blockley, S. P., Blunier, T., Buchardt, S. L., Clausen, H. B., Cvijanovic, I., Dahl-Jensen, D., Johnsen, S. J., and Fischer, H.: A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: refining and extending the INTIMATE event stratigraphy, Quaternary Sci. Rev., 106, 14–28, 2014. 

Reimer, P., Bard, E., Bayliss, A., Beck, J. W., Blackwell, P. G., Bronk Ramsey, C., Buck, C. E., Cheng, H., Edwards, R. L., Friedrich, M., Grootes, P. M., Guilderson, T. P., Haflidason, H., Hajdas, I., Hatté, C., Heaton T. J., Hoffmann, D. L., Hogg, A. G., Hughen, K. A., Kaiser, K. F., Kromer, B., Manning, S. W., Niu, M., Reimer, R. W., Richards, D. A., Scott, E. M., Southon, J. R., Staff, R. A., Turney, C. S. M., and van der Plicht, J.: IntCal13 and Marine13 radiocarbon age calibration curves 0–50,000 years cal BP, Radiocarbon, 55, 1869–1887, 2013. 

Rhein, M., Kieke, D., and Steinfeldt, R.: Advection of North Atlantic Deep Water from the Labrador Sea to the southern hemisphere, J. Geophys. Res.-Ocean., 120, 2471–2487, 2015. 

Rühlemann, C., Frank, M., Hale, W., Mangini, A., Mulitza, S., Müller, P., and Wefer, G.: Late Quaternary productivity changes in the western equatorial Atlantic: Evidence from 230Th-normalized carbonate and organic carbon accumulation rates, Mar. Geol., 135, 127–152, 1996. 

Schmittner, A., Bostock, H. C., Cartapanis, O., Curry, W. B., Filipsson, H. L., Galbraith, E. D., Gottschalk, J., Herguera, J. C., Hoogakker, B., Jaccard, S. L., Lisiecki, L. E., Lund, D. C., Martínez-Méndez, G., Lynch-Stieglitz, J., Mackensen, A., Michel, E., Mix, A. C., Oppo, D. W., Peterson, C. D., Repschläger, J., Sikes, E. L., Spero, H. J., and Waelbroeck, C.: Calibration of the carbon isotope composition (δ13C) of benthic foraminifera, Paleoceanography, 32, 512–530,, 2017. 

Scholten, J., Botz, R., Paetsch, H., and Stoffers, P.: 230Thex flux into Norwegian-Greenland Sea sediments: Evidence for lateral sediment transport during the past 300,000 years, Earth Planet. Sc. Lett., 121, 111–124, 1994. 

Schott, F. A., Dengler, M., Brandt, P., Affler, K., Fischer, J., Bourles, B., Gouriou, Y., Molinari, R. L., and Rhein, M.: The zonal currents and transports at 35 W in the tropical Atlantic, Geophys. Res. Lett., 30, 1–4, 2003. 

Sepulcre, S., Durand, N., and Bard, E.: Large 14C age offsets between the fine fraction and coexisting planktonic foraminifera in shallow Caribbean sediments, Quaternary Geochronol., 38, 61–74, 2017. 

Taylor, S. and McLennan, S.: The continental crust: its composition and evolution, Oxford, Blackwell Press, 312 pp., 1985. 

Thomas, A. L., Henderson, G. M., and Robinson, L. F.: Interpretation of the 231Pa 230Th paleocirculation proxy: New water-column measurements from the southwest Indian Ocean, Earth Planet. Sc. Lett., 241, 493–504, 2006. 

Thomson, J., Colley, S., Anderson, R., Cook, G., MacKenzie, A., and Harkness, D.: Holocene sediment fluxes in the northeast Atlantic from 230Thexcess and radiocarbon measurements, Paleoceanography, 8, 631–650, 1993. 

Thornalley, D. J. R., Barker, S., Broecker, W. S., Elderfield, H., and McCave, I. N.: The Deglacial Evolution of North Atlantic Deep Convection, Science, 331, 202–205, 2011. 

Torrence, C. and Compo, G. P.: A practical guide to wavelet analysis, B. Am. Meteorol. Soc., 79, 61–78, 1998. 

Torrence, C. and Webster, P. J.: Interdecadal changes in the ENSO–monsoon system, J. Clim., 12, 2679–2690, 1999. 

Vazquez Riveiros, N., Waelbroeck, C., Roche, D. M., Moreira, S., Burckel, P., Dewilde, F., Skinner, L., Böhm, E., Arz, H. W., and Dokken, T.: Northern origin of western tropical Atlantic deep waters during Heinrich Stadials, submitted, 2018. 

Vidal, L., Labeyrie, L., Cortijo, E., Arnold, M., Duplessy, J. C., Michel, E., Becqué, S., and van Weering, T. C. E.: Evidence for changes in the North Atlantic Deep Water linked to meltwater surges during the Heinrich events, Earth Planet. Sc. Lett., 146, 13–26, 1997. 

Waelbroeck, C., Duplessy, J.-C., Michel, E., Labeyrie, L., Paillard, D., and Duprat, J.: The timing of the last deglaciation in North Atlantic climate records, Nature, 412, 724–727, 2001. 

Weltje, G. J. and Tjallingii, R.: Calibration of XRF core scanners for quantitative geochemical logging of sediment cores: theory and application, Earth Planet. Sc. Lett., 274, 423–438, 2008.  

Wheatcroft, R. A.: Experimental tests for particule size-dependent bioturbation in the deep ocean, Limnol. Oceanogr., 37, 90–104, 1992. 

Wilson, K. E., Maslin, M. A., and Burns, S. J.: Evidence for a prolonged retroflection of the North Brazil Current during glacial stages, Palaeogeogr. Palaeocl., 301, 86–96, 2011. 

Woerther, P. and Bourillet, J. F.: Exploitation des mesures faites avec les accéléromètres sur le carottier CAPYPSO-Mission SEDICAR4-ALIENOR, Ifremer, Brest, 47 pp., 2005. 

Yu, E.-F., François, R., and Bacon, M.: Similar rates of modern and last-glacial ocean thermohaline circulation inferred from radiochemical data, Nature, 379, 689–694, 1996. 

Zahn, R., Winn, K., and Sarnthein, M.: Benthic foraminiferal δ13C and accumulation rates of organic carbon: Uvigerina peregrina group and Cibicidoides wuellerstorfi, Paleoceanography, 1, 27–42, 1986. 

Zhang, D., Msadek, R., McPhaden, M. J., and Delworth, T.: Multidecadal variability of the North Brazil Current and its connection to the Atlantic meridional overturning circulation, J. Geophys. Res.-Ocean., 116, 1–9, 2011. 

Zhang, Y., Chiessi, C. M., Mulitza, S., Zabel, M., Trindade, R. I., Hollanda, M. H. B., Dantas, E. L., Govin, A., Tiedemann, R., and Wefer, G.: Origin of increased terrigenous supply to the NE South American continental margin during Heinrich Stadial 1 and the Younger Dryas, Earth Planet. Sc. Lett., 432, 493–500, 2015. 

Short summary
Recording the precise timing and sequence of events is essential for understanding rapid climate changes and improving climate model predictive skills. Here, we precisely assess the relative timing between ocean and atmospheric changes, both recorded in the same deep-sea core over the last 45 kyr. We show that decreased mid-depth water mass transport in the western equatorial Atlantic preceded increased rainfall over the adjacent continent by 120 to 980 yr, depending on the type of climate event.