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

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 (04^∘14.7^′ S, 36^∘21.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 (04^∘14.7^′ S, 36^∘20.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 (¹⁴C) age model shows that increases in sedimentary Ti ∕ Ca are indeed synchronous with decreases in South American speleothem δ¹⁸O (Burckel et al., 2015). Based on the observed synchronicity, composite age models of core MD09-3257 and GeoB3910-2 have been developed using ¹⁴C dating for the past 35 kyr, combined with the alignment of sediment Ti ∕ Ca increases with decreases in speleothem δ¹⁸O 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 2 MD09-3257 Pa ∕ Th, ln(Ti ∕ Ca) and composite C. wuellerstorfi δ¹³C (δ¹³C_Cw) 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 δ¹³C_Cw 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 δ¹⁸O (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 ¹³C ∕ ¹²C (δ¹³C, 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 (δ¹⁸O = −2.20 and δ¹³C $=+\mathrm{1.95}$ ‰). The mean external reproducibility (1σ) of carbonate standards is ±0.03 ‰ for δ¹³C; measured NBS-18 δ¹⁸O is $-\mathrm{23.27}\pm \mathrm{0.10}$ and δ¹³C is $-\mathrm{5.01}\pm \mathrm{0.03}$ ‰ VPDB. Core GeoB3910-2 C. wuellerstorfi gδ¹³C 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 δ¹³C. 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 ¹³C/¹²C isotopic ratio (δ¹³C, expressed in per mil –  ‰ – versus VPDB) of the epifaunal benthic foraminifer C. wuellerstorfi has been shown to record the δ¹³C 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 ¹²C by primary productivity, which increases dissolved δ¹³C), as well as temperature dependent air–sea exchanges (Lynch-Stieglitz et al., 1995). DIC δ¹³C subsequently decreases as deep water ages, due to the progressive remineralization at depth of relatively ¹³C-depleted biogenic material. As a result, DIC δ¹³C largely follows water mass structure and circulation in the modern ocean, and C. wuellerstorfi δ¹³C (δ¹³C_Cw 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 δ¹³C more faithfully follows water oxygen content than phosphate content (Eide et al., 2017), lending strong support to the use of δ¹³C_Cw 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 (²³¹Pa${}_{xs,\mathrm{0}}{/}^{\mathrm{230}}$Th_xs,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. ²³⁰Th and ²³¹Pa 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 (²³⁸U ∕ ²³²Th) 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 (Neptune^Plus, 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 (²³⁸U ∕ ²³²Th) 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). ²³¹Pa and ²³⁰Th are produced at a constant Pa ∕ Th activity ratio of 0.093 by dissolved uranium, which is homogeneously distributed in oceans. ²³⁰Th is much more particle reactive than ²³¹Pa, as reflected by their respective residence time in the ocean (30–40 yr for ²³⁰Th and 100–200 yr for ²³¹Pa, François, 2007). Therefore, ²³⁰Th is rapidly removed from the water column to the underlying sediment, while ²³¹Pa can be advected by oceanic currents. Thus, when averaged over an entire ocean basin, high (low) flow rates result in high (low) ²³¹Pa export and a low (high) sedimentary Pa ∕ Th ratio. In contrast to δ¹³C_Cw, which records the DIC δ¹³C 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, ²³¹Pa 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 ²³⁰Th-normalized ²³²Th flux, hereafter simply referred to as the ²³²Th flux, is indicative of the vertical terrigenous flux to the core site. As ²³²Th 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 ²³²Th flux (r=0.21, p=0.07) (Figs. S2 and S3). In contrast, because the correlation between Pa ∕ Th and the ²³²Th 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 ²³⁰Th-normalized opal flux measurements in MD09-3257, which are below 0.06 g cm⁻² kyr⁻¹ (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, δ¹³C_Cw 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 W^X and W^Y, the cross-wavelet spectrum is defined as $W^{XY}=W^X{W}^{Y*}$, where W^Y* is the complex conjugate of W^Y (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(W^XY) 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 δ¹³C_Cw 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 δ¹⁸O 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 δ¹³C_Cw 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, δ¹³C_Cw and Ti ∕ Ca

Pa ∕ Th, δ¹³C_Cw 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 δ¹³C_Cw versus Ti ∕ Ca or δ¹³C_Cw versus Pa ∕ Th. Apart from the unavoidable uncertainty introduced by bioturbation, the relative dating uncertainty between the δ¹³C_Cw 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, δ¹³C_Cw 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 3 Spearman correlation coefficient of δ¹³C_Cw versus Pa ∕ Th (blue curves), of δ¹³C_Cw 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., δ¹³C_Cw 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 δ¹³C_Cw versus ln(Ti ∕ Ca), and of δ¹³C_Cw 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 δ¹³C_Cw and Ti ∕ Ca, and that δ¹³C_Cw 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 δ¹³C_Cw 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 δ¹³C_Cw and Ti ∕ Ca remained between −100 and +100 yr.

Figure 4 Cross-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 δ¹³C_Cw and Pa ∕ Th (Fig. 5), and between δ¹³C_Cw 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 1 Relative 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.

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 δ¹³C_Cw 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 δ¹³C_Cw lags Pa ∕ Th by 279±244 yr (1σ) for periodicities of 1000 to 3000 yr over 28–40 ka, but that the lag of δ¹³C_Cw 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 δ¹³C_Cw 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 δ¹³C_Cw and ln(Ti ∕ Ca) over these regions are not significantly different from zero (Fig. 6d and Table 1), indicating that decreases in δ¹³C_Cw are in phase with increases in ln(Ti ∕ Ca) within uncertainties.

Figure 5 Cross-wavelet transform of δ¹³C_Cw composite record versus MD09-3257 Pa ∕ Th. δ¹³C_Cw values have been multiplied by −1 to allow a straightforward reading of the relative phase between a decrease in δ¹³C 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 δ¹³C_Cw, we also accounted for the additional error due to the combining of the MD09-3257 and GeoB3910-2 δ¹³C_Cw 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 6 Cross-wavelet transform of δ¹³C_Cw composite record versus MD09-3257 ln(Ti ∕ Ca). δ¹³C_Cw values have been multiplied by −1 to allow a straightforward reading of the relative phase between a decrease in δ¹³C 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 = 10⁶ m³ s⁻¹). 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⁻² yr⁻¹, 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, δ¹³C_Cw 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 δ¹³C_Cw, 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 δ¹³C_Cw 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 δ¹³C_Cw lags Pa ∕ Th by 279±244 yr at D–O frequencies, and that decreases in δ¹³C_Cw 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 δ¹³C_Cw 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 δ¹³C_Cw 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 δ¹³C_Cw 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 δ¹³C_Cw 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 ²³⁰Th 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⁻¹, 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 ²³⁰Th and 100–200 yr for ²³¹Pa, 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 δ¹³C_Cw with respect to Pa ∕ Th

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

δ¹³C_Cw 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 δ¹³C_Cw. 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 δ¹³C_Cw 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 δ¹³C_Cw and ln (Ti ∕ Ca). Thus, we will not discuss the results of the cross-wavelet analyses involving δ¹³C_Cw 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 δ¹³C_Cw, 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 δ¹⁸O, DIC δ¹³C and sedimentary Pa ∕ Th.