Introduction

Ocean circulation constitutes a major component of global surface heat transport (Ganachaud and Wunsch 2000) and, as such, plays a significant role in determining the Earth’s climate. However, despite improvements in our understanding of the role of the oceans in climate, significant uncertainty remains regarding how the oceans behave as a climate feedback (Marshall et al. 2015), particularly in the context of future climate predictions in response to anthropogenic forcings (e.g. Cheng et al. 2013). To resolve this uncertainty, one approach is to unravel the role ocean circulation plays in abrupt climate variability of the past. Of particular interest is variability in the Atlantic meridional overturning circulation (AMOC), a major influencer of the climate of the North Atlantic and the rest of the globe (e.g. Jackson et al. 2015), as AMOC transport can be significantly reduced under high freshwater fluxes to regions of deep water formation (e.g. Zhang and Delworth 2005; Stouffer et al. 2006).

At present, the zonally integrated AMOC comprises a strong upper overturning cell ventilated from the north by deep water formation in the subpolar North Atlantic and marginal seas (North Atlantic Deep Water, NADW), overlaying a weaker overturning cell ventilated by dense bottom waters from the Southern Ocean (present in the eastern North Atlantic as Lower Deep Water, LDW, a highly modified Antarctic Bottom Water), with a boundary between the two zonally averaged cells in the North Atlantic at a depth of c. 3.5–4 km (Buckley and Marshall 2015), although it is somewhat shallower in the eastern North Atlantic at about 3 km depth (van Aken 2000). When considering deep and bottom water masses in the eastern North Atlantic, one can therefore distinguish between water masses originating from the north (northern-sourced water, NSW) in the form of NADW today and those originating from the south (southern-sourced water, SSW) in the form of LDW. It is important to note that there is considerable horizontal structure to the AMOC that is obfuscated by this zonally averaged view, with the main component of meridional transport in the North Atlantic occurring along the western boundary (Buckley and Marshall 2015).

Two time periods of particular interest for unravelling the role of ocean circulation in abrupt climate variability are the Last Glacial Maximum (LGM) and Heinrich Stadial 1 (HS1). Despite considerable ongoing discussion, there is a general consensus that Atlantic overturning at the LGM was shallow (e.g. Curry and Oppo 2005; Marchitto and Broecker 2006; Lippold et al. 2016) and strong (e.g. Praetorius et al. 2008; Böhm et al. 2015; Freeman et al. 2016), with an active, shoaled upper cell and an expanded, sluggish lower cell (e.g. Freeman et al. 2016; Lippold et al. 2016). Estimates for the extent of shoaling vary from a kilometre or more (Oppo and Lehman 1993; Marchitto and Broecker 2006) to unchanged from Holocene values (Gebbie 2014).

The extent of different water masses can be tracked through palaeotracers, principally in the form of nutrient tracers such as (foraminiferal) δ13C (Curry and Oppo 2005), radiocarbon (e.g. Thornalley et al. 2011) and Cd/Ca ratios (e.g. Rickaby and Elderfield 2005). Recently, however, an inorganic tracer in the form of neodymium isotopes (expressed using the εNd notation) has come into widespread use (e.g. Arsouze et al. 2009; Piotrowski et al. 2012; Rempfer et al. 2012; Roberts and Piotrowski 2015; Pöppelmeier et al. 2018). In contrast to nutrient tracers, εNd is assumed to behave quasi-conservatively in the oceans (Jones et al. 2008; van de Flierdt et al. 2016).

However, this assumption of quasi-conservative behaviour is equivocal. Neodymium enters the oceans in dissolved and particulate phases through continental runoff and atmospheric fallout, and, like many other trace elements, is scavenged in the water column through adsorption onto particulates. This results in a mean ocean residence time on the order of hundreds of years (Tachikawa et al. 2003; Arsouze et al. 2007), although a recent study puts the figure at closer to 100 years (Hayes et al. 2018). However, adsorbed neodymium can also undergo desorption, resulting in continuous exchange between adsorbed (Ndads) and dissolved (Ndd) neodymium reservoirs as particles settle through the water column, in a process called reversible scavenging (Bacon and Anderson 1982).

The potential significance of localised highly non-conservative behaviour of neodymium isotopes for palaeoceanography has been raised by recent geochemical studies on sediment cores from the eastern North Atlantic (Roberts and Piotrowski 2015; Blaser et al. 2019). The extensive sea-ice in the North Atlantic during the last glaciation resulted in the widespread deposition of ice-rafted debris (IRD). Layers of exceptionally high IRD content originating from the Laurentide Ice Sheet were termed Heinrich Layers, representing depositional events during the last glacial and deglaciation known as Heinrich Stadials (Heinrich 1988; Broecker et al. 1992), most recently Heinrich Stadial 1 and 2 at c. 17 ka and 24 ka BP respectively (Fig. 1). During these Heinrich Stadials (HS), the IRD flux in the northeast Atlantic exceeded the ambient flux during the last glacial by up to an order of magnitude (McManus et al. 1998; Haapaniemi et al. 2010; Obrochta et al. 2012). Furthermore, there is a strong consensus that overturning in the upper (northern) cell was considerably reduced during HS1 (McManus et al. 2004; Bradtmiller et al. 2014; Lynch-Stieglitz et al. 2014).

Fig. 1
figure 1

Lithic grain counts from DSDP 609 (eastern North Atlantic) across HS1 and HS2 compared to the intervening LGM (dashed) and a 2 ka moving average (solid) (adapted from Bond et al. 2012)

Full size image

Roberts and Piotrowski (2015) found very high (radiogenic) εNd throughout the water column in the eastern North Atlantic during deglaciation. Blaser et al. (2019) extended this study and suggested that the underlying SSW appeared to have developed a Nd isotopic signature similar to that of overlying NSW. Both authors concluded that these observations are challenging to interpret under assumptions of conservative behaviour (Fig. 2), with Roberts and Piotrowski (2015) suggesting that dissolution of radiogenic IRD could be the cause. Blaser et al. (2019) found that vertical εNd profiles during exceptionally high IRD flux regimes such as Heinrich Stadials (HSs) deviated even further from a hypothetical conservative profile and argued that IRD may instead act as a vertical shuttle for neodymium, leading to the overprinting of abyssal SSW neodymium isotopic signatures with signals from further up in the water column through an acceleration of reversible scavenging. This suggestion is bolstered by additional recent observations of apparently non-conservative behaviour of εNd in high sediment flux regions, e.g. in the deep waters in the modern Panama Basin (Grasse et al. 2017) and Bay of Bengal (Yu et al. 2017).

Fig. 2
figure 2

εNd reconstructions from the Dreizack Seamount at the LGM (dark blue) and HE1b/2b (light blue), adapted from Blaser et al. (2019) Fig. 8. The predicted locations and water mass boundary between northern and southern component waters (North Atlantic Deep Water [NADW] and Atlantic Bottom Water [ABW] respectively) and a predicted ‘conservative’ εNd profile are included

Full size image

Since most interpretations of neodymium isotopes as a palaeotracer operate under assumptions of dominantly conservative behaviour, it is crucial to check the viability of this assumption. The key hypothesis we investigate is therefore this—could geologically reasonable fluxes of IRD result in a dominant vertical transport of Nd and thus the overprinting of bottom water εNd signatures? In this paper, we present new palaeoenvironmental data from the Me68-89 sediment core (Blaser et al. 2019). We then present results from a simple proof-of-concept box model for neodymium reversible scavenging in order to constrain the hypothetical magnitude of scavenging required to alter the εNd signal accordingly.

Materials and methods

Sediment cores

As detailed in Blaser et al. (2019), three cores from the Dreizack seamount (Me68-91, Me68-89 and Po08-23) are combined with BOFS 5K, 8K, 10K and 11K, and IODP Site U1308 at a morphological high in the Maury Channel, spanning a depth range of 3547–4470 m between 47–53° N. HS1 and HS2 are prominent in all cores and this paper specifically focuses on HS1, which has been further subdivided into a Hudson Strait event (represented by Heinrich Layer 1a) and a ‘European’ precursor (HL1b) (Blaser et al. 2019), which is identifiable as a detrital layer in these cores. εNd leachate data from Me68-89 (Blaser et al. 2019) reveals the large radiogenic excursion centred towards the end of HL1a. As discussed by Blaser et al. (2019), the εNd recorded at Me68-89 during HS1 is comparable to the records from other cores from below 3.5 km depth, as part of an apparent homogenisation of the abyssal εNd signature.

Here, to complement εNd measurements from Me68-89, IRD and planktonic foraminifera counts were made from the washed > 63 μm fraction at 1 cm intervals spanning HS1. Sediment core Me68-89 (47° 26.10′ N, 19° 33.40′ W, 4260 m water depth) is a core of intermediate depth at the Dreizack Seamount and is therefore relatively unaffected by turbidite flows in comparison to core Me68-91 (Blaser et al. 2019). IRD counts were subdivided into the four key lithologies identified: unstained quartz clasts, haematite-stained quartz clasts, mafic volcanic clasts and carbonate clasts. The total number of planktonic foraminifera (pf) was counted, as well as Neogloboquadrina pachyderma sinistral (NPS).

Model design

Overview

The processes included in our model and key domain parameters and boundary conditions are summarised in Fig. 3 and Table 1. The model broadly follows the approach of Luo et al. (2010). The model represents an area of dimensions 600 × 600 km in the eastern North Atlantic. The neodymium concentration [Nd] (pmol kg−1) and neodymium isotopic composition εNd are calculated for a vertical array of boxes representing the water column. In each box, there are two reservoirs of neodymium, dissolved (Ndd) and adsorbed (Ndads), which exchange with adsorption and desorption rate constants of K and K− (s−1) respectively, and are also influenced by advection. The adsorbed reservoir sinks through the water column due to particulate settling, transferring the Ndads signal downwards. There is no vertical component of advection, but vertical exchange of the dissolved reservoirs is possible through vertical diffusion. Sediment settling out of the water column enters a sediment box, where it can then exchange with a pore water–dissolved neodymium reservoir through a particle concentration factor fc and particulate dissolution fraction fp. This pore water can interact with the water column through a benthic ‘piston velocity’ w∗ (m s−1), used by Abbott et al. (2015) to relate pore water fluxes into the water column to overlying [Nd]. Finally, sediment entering the water column from the surface can be ‘pre-loaded’ with neodymium. These processes are schematically represented in Fig. 4. The model integrates the equations below to steady state, with constants and variables defined in Table 1.

$$ \frac{\mathbf{\partial}}{\mathbf{\partial}\boldsymbol{t}}{\left[\mathbf{Nd}\right]}_{\boldsymbol{i}}^{\boldsymbol{d}}=\frac{v_i^{\ast }}{\hat{y}}\kern0.2em \left\{{\left[\mathrm{Nd}\right]}_b^d-{\left[\mathrm{Nd}\right]}_i^d\right\}+{K}^{-}{\left[\mathrm{Nd}\right]}_i^{ads}-K{\left[\mathrm{Nd}\right]}_i^d+\frac{1}{\hat{z}}\left({\hat{J}}_i-{\overset{\check{} }{J}}_i\right)+{\delta}_{z=0}\frac{w^{\ast }}{\hat{y}}{\left[\mathrm{Nd}\right]}_b^d $$
$$ \frac{\mathbf{\partial}}{\mathbf{\partial}\boldsymbol{t}}{\left[\mathbf{Nd}\right]}_{\boldsymbol{i}}^{\boldsymbol{ads}}=\frac{v_i^{\ast }}{\hat{y}}\ \left\{{\left[\mathrm{Nd}\right]}_b^{ads}-{\left[\mathrm{Nd}\right]}_i^{ads}\right\}+\frac{s_i^{\ast }}{\hat{z}}\ \left\{{\left[\mathrm{Nd}\right]}_{i-1}^{ads}-{\left[\mathrm{Nd}\right]}_i^{ads}\right\}+K{\left[\mathrm{Nd}\right]}_i^d-{K}^{-}{\left[\mathrm{Nd}\right]}_i^{ads}+{\delta}_{seabed}\kern0.55mm \frac{s^{\ast }}{\hat{z}}{\left[\mathrm{Nd}\right]}_b^{ads} $$
$$ {\left[\mathrm{Nd}\right]}_i^d\frac{\mathbf{\partial}}{\mathbf{\partial}\boldsymbol{t}}\boldsymbol{\varepsilon} \mathbf{N}{\mathbf{d}}_{\boldsymbol{i}}^{\boldsymbol{d}}=\frac{v_i^{\ast }}{\hat{y}}\ \left\{{\left[\mathrm{Nd}\right]}_b^d\varepsilon N{d}_b^d-{\left[\mathrm{Nd}\right]}_i^d\varepsilon {\mathrm{Nd}}_i^d\right\}+{K}^{-}{\left[\mathrm{Nd}\right]}_i^{ads}\varepsilon {\mathrm{Nd}}_i^{ads}-K{\left[\mathrm{Nd}\right]}_i^d\varepsilon {\mathrm{Nd}}_i^d+\frac{1}{\hat{z}}{\hat{J}}_i\varepsilon {\mathrm{Nd}}_{i-{1}^{+},{i}^{-}}^d-\frac{1}{\hat{z}}{\overset{\check{} }{J}}_i\varepsilon {\mathrm{Nd}}_{i^{+},i+{1}^{-}}^d+{\delta}_{z=0}\frac{w^{\ast }}{\hat{y}}{\left[\mathrm{Nd}\right]}_b^d\varepsilon {\mathrm{Nd}}_b^d+{\Delta}_i^d $$
$$ {\left[\mathrm{Nd}\right]}_i^{ads}\frac{\boldsymbol{\partial}}{\boldsymbol{\partial t}}\boldsymbol{\varepsilon} {\mathbf{Nd}}_{\boldsymbol{i}}^{\boldsymbol{ads}}=\frac{v_i^{\ast }}{\hat{y}}\ \left\{{\left[\mathrm{Nd}\right]}_b^{ads}\varepsilon {\mathrm{Nd}}_b^{ads}-{\left[\mathrm{Nd}\right]}_i^{ads}\varepsilon {\mathrm{Nd}}_i^{ads}\right\}+\frac{s^{\ast }}{\hat{z}}\left\{{\left[\mathrm{Nd}\right]}_{i-1}^{ads}\varepsilon {\mathrm{Nd}}_{i-1}^{ads}-{\left[\mathrm{Nd}\right]}_i^{ads}\varepsilon {\mathrm{Nd}}_i^{ads}\right\}+K{\left[\mathrm{Nd}\right]}_i^d\varepsilon {\mathrm{Nd}}_i^d-{K}^{-}{\left[\mathrm{Nd}\right]}_i^{ads}\varepsilon {\mathrm{Nd}}_i^{ads}++{\delta}_{seabed}\frac{s^{\ast }}{\hat{z}}{\left[\mathrm{Nd}\right]}_b^{ads}\varepsilon {\mathrm{Nd}}_b^{ads}+{\Delta}_i^{ads} $$

where \( {\Delta}_i^x=-\varepsilon {\mathrm{Nd}}_i^x\frac{\partial }{\partial t}{\left[\mathrm{Nd}\right]}_i^x\sim 0 \) at equilibrium.

Fig. 3
figure 3

a The key dynamical processes involved in the cycling of neodymium in the open glacial and deglaciating ocean, integrated into the model. b Boundary conditions relevant to the eastern North Atlantic during the glacial and deglaciating ocean, including model domain dimensions and a schematic velocity profile

Full size image
Table 1 Model variables (and constant values) for simulations
Full size table
Fig. 4
figure 4

The model consists of 50 vertical boxes of dimensions [5 × 103m × 6 × 106m × 6 × 106 m]. Within each box is a dissolved (Ndd) and adsorbed (Ndads) neodymium reservoir, for which concentration [Nd] (pmol kg−1water) and isotopic composition εNd are stored. See Table 1 for fluxes

Full size image

Velocity profile

A synthetic velocity (advection) profile is used, representing the key characteristics of mean advection profiles during the LGM. This profile is presented in Fig. 5. Constraints on the LGM water column velocity profile in the eastern North Atlantic are poor, so the ‘LGM-like’ profile reflects the hypothesis that overturning during the LGM was strong and shoaled versus the modern (e.g. Lippold et al. 2012). There is however debate regarding the extent of shoaling. Many studies, particularly those based on nutrient proxies, have argued for a significant shoaling of the NSW-SSW interface up to 2000–2500 m (e.g Curry and Oppo 2005; Marchitto and Broecker 2006; Lippold et al. 2016). However, it has been argued that nutrient proxies exaggerate the extent of shoaling (Gherardi et al. 2009; Lippold et al. 2012; Gebbie 2014; Howe et al. 2016a), and that there was strong stratification between NSW and SSW, contrary to simple mass-balance estimates from δ13C (Yu et al. 2008; Gebbie 2014). It is also likely that the extent of shoaling may have been spatially variable, for instance a shallower boundary at the Rockall Plateau than the central basin (e.g. van Aken 2000). We therefore take a more conservative view here with a glacial NSW-SSW boundary at 3500 m, albeit with the recognition that this figure is uncertain (as with the entire advection profile). For comparison, a synthetic ‘shoaled’ LGM profile has also been produced with a NSW-SSW boundary at 2500 m. These LGM profiles are normalised to represent the same overturning strength as the late Holocene profile, in line with the general consensus that net LGM overturning strength was similar to today (Yu et al. 1996; McManus et al. 2004).

Fig. 5
figure 5

Velocity profiles used in the model runs in this paper

Full size image

Constraints on circulation during Heinrich Stadials are similarly poor. Whilst early evidence was interpreted as representing an AMOC cessation (e.g. McManus et al. 2004), it has since become clear that overturning continued, albeit significantly weakened (Bradtmiller et al. 2014; Lynch-Stieglitz et al. 2014; Wilson et al. 2014; Chen et al. 2015; Wei et al. 2016). To simulate a HS, the LGM-like advection velocity profile is scaled by a factor of 0.2 to represent qualitatively similar but significantly weakened circulation (Fig. 5).

Scavenging rates

Exchange between the dissolved and adsorbed reservoirs is assumed to be proportional to adsorption and desorption rate constants, following the approach of Luo et al. (2010). At equilibrium, in a static water column:

$$ K{\left[\mathrm{Nd}\right]}_{\mathrm{d}}={K}^{-}{\left[\mathrm{Nd}\right]}_{\mathrm{ads}} $$
$$ {K}_{eqm}=\frac{K}{K^{-}}=\frac{{\left[\mathrm{Nd}\right]}_{\mathrm{ads}}}{{\left[\mathrm{Nd}\right]}_{\mathrm{d}}} $$

(concentrations are per unit mass of water)

The ratio Keqm therefore gives the distribution of Nd between the dissolved and adsorbed reservoirs at equilibrium. For example, if K = K−, at equilibrium Nd will be equally distributed between the dissolved and adsorbed reservoirs. Note that since all concentrations are per unit mass of water, the concentration of Nd in the adsorbed phase will be \( \frac{1}{\alpha }{\left[\mathrm{Nd}\right]}_{\mathrm{ads}} \) where α is the adsorbed particulate mass fraction. In the absence of empirical constraints, it is assumed that K = K−.

Since scavenging is assumed to be the only Nd sink in this model, at steady state and in the absence of advection, the scavenging rates correspond to the residence time of Nd in the water column, i.e.

$$ {\left[\mathrm{Nd}\right]}_{\mathrm{d}}(t)={\left[\mathrm{Nd}\right]}_{\mathrm{d}}^0\ {e}^{- Kt} $$

It is important to note that K is not a material property, it is a bulk parameterisation for all the factors that impact the rate of exchange of Nd between the dissolved and adsorbed reservoirs. Therefore, K will be a function of material adsorption and desorption rates for the various components that make up IRD, as well as the IRD concentration in the water column. K and K− are not necessarily equal and may vary as a function of time there are no constraints on this variability. Furthermore, the scavenging rate is assumed to be directly proportional to particle concentration, although this is debatable (Tachikawa et al. 1997). It is assumed that particle concentration is proportional to particle flux, which relies on constant settling rate and therefore constant particle size. Here, only one average type of particle is modelled. Given the mineralogy of IRD (e.g. Hodell et al. 2017), IRD particles are assumed not to dissolve over settling timescales, that is to say that Nd exchange between particles and the water column is restricted to scavenging interactions on particle surfaces as opposed to bulk particle dissolution.

Boundary conditions

Boundary conditions are based on Blaser et al. (2019). NSW and SSW are assumed to have constant εNd of − 5 and − 10 respectively, the former signature derived from Roberts and Piotrowski (2015) and Howe et al. (2016b). Concentration profiles are based on modern profiles from van de Flierdt et al. 2012. Specifically, both NSW and SSW [Nd]d increase linearly with depth from 15 pmol kg−1 at the surface to 25 pmol kg−1 at the seafloor. It is uncertain whether this would also be the case at the LGM but there is no good way of reconstructing palaeo-[Nd]d at present and the model εNd sensitivity to [Nd]d is quite low.

In the model runs below, it is assumed that all adsorbed Nd is sourced from adsorption within the model domain, i.e. no pre-loaded particle-reactive neodymium exists in particles entering the domain. The justifications for this assumption are (i) that the rate of scavenging within the IRD belt is significantly greater than outside the IRD belt and (ii) that scavenging only occurs at the surface of particles so there is no exchange between the dissolved Nd reservoir and detrital particulate Nd.

Results from Me68-89

A peak in the detrital carbonate fraction as well as a minimum in the pf:IRD ratio is centred around the previously established HL1a interval, validating these earlier results and representing a significant detrital carbonate IRD flux of Laurentide origin.

The planktonic foraminiferum Neogloboquadrina pachyderma sinistral (NPS) is an established proxy for polar sea-surface temperature (e.g. Ericson 1959; Kohfeld et al. 1996). The NPS (measured as a fraction of the total pf count) data from Me68-89 (Fig. 6) is noisy and does not appear to characterise HL1a as is the case for the pf:IRD ratio. If NPS faithfully represents cold SSTs and the pf:IRD ratio represents the sediment dilution by the accompanying high IRD flux, these data would be expected to be negatively correlated; however, there is no statistically significant correlation. This could indicate that SST perturbations at site Me68-89 were not significant enough during HS1 to significantly impact the NPS abundance. There is nevertheless a generally declining trend from the peak glacial conditions preceding HL1b to the deglacial conditions postdating HL1a, likely representing the ongoing trend of deglaciation in the eastern North Atlantic.

Fig. 6
figure 6

IRD mineralogy in core Me68-89, normalised to unstained quartz (qz), haematite-stained quartz (haem), volcanic (volc) and carbonate (carb) fractions; εNd (Blaser et al. 2019); N. pachyderma S proportion (%nps); IRD-normalised planktonic foram count (#pf/#IRD). Also plotted for comparison is the 231Pa/230Th record from the Iberian margin SU81-18 site (Gherardi et al. 2005) and IRD counts from the eastern North Atlantic DSDP 609 site (Bond et al. 2012), which have been scaled to match the HS1 intervals identified from the three cores. However, these data are provided for qualitative comparison only as there is no precise age-model for Me68-89

Full size image

A simple scavenging model for the glacial and deglacial eastern North Atlantic

Using our ‘deep’ (hereon default) LGM-like advection profile (Fig. 5) with the widely-assumed residence time of centuries (Tachikawa et al. 2003), a sharp boundary between northern and southern component water masses is reproduced, reflecting conservative behaviour (Figs. 2 and 7a), deviating significantly from the observations. It is possible to recreate the key features of the bottom water εNd profile as reported by Blaser et al. (2019), but only by setting the scavenging intensity to an equivalent residence time of τ = 8 years (Fig. 7b). This results in a gradual transition in εNdd from close to − 5 (the εNdd of the NSW, i.e. [glacial] North Atlantic Deep Water) at the NSW-SSW interface to − 10 at a depth of 5000 m. The εNd of the adsorbed reservoir (εNdads) retains a NSW-like signature for longer, only reaching − 7 at 5000 m. Model-data agreement is more challenging using the significantly shoaled LGM velocity profile, with the model struggling to achieve a steep enough εNdd gradient at abyssal depths (Fig. 7c) and requiring a greater scavenging rate to do so.

Fig. 7
figure 7

Modelled steady-state εNd profiles, a with an LGM-like velocity profile (Fig. 5) with a residence time (scavenging intensity) τ = 300 years; b as with a but with τ = 8 years, tuned to fit the profile compiled by Blaser et al. (2019) and Roberts and Piotrowski (2015) (data points); c with a shoaled LGM-like velocity profile (Fig. 5) and τ = 2 years; and d with a HS-like velocity profile (Fig. 5) and τ = 1 year

Full size image

To homogenise the bottom water εNdd and therefore recreate the reconstructed profiles for HS1/2 (Blaser et al. 2019), the scavenging intensity must be increased by approximately an order of magnitude, to an equivalent residence time of τ = 1 year (Fig. 7d). With such a low residence time, the dissolved and adsorbed reservoirs are close to equilibrium at all depths. Whilst these low residence times differ dramatically from assumed whole-ocean Nd residence times (Tachikawa et al. 2003), localised low residence times may be possible (see discussion, “Benthic Nd flux”).

The effects of modifying the two key physical parameters of scavenging intensity and advection can be further explored by plotting the difference in εNdd (ΔεNdd) across the SSW layer (3.5–5.0 km) as a function of these two parameters (Fig. 8). In the LGM- and HS-like advection profiles, NSW and SSW carry a εNdd of − 10 and − 5 respectively, so perfectly conservative behaviour would result in a trans-SSW ΔεNdd of − 5 (e.g. close to Fig. 7b). Similarly, completely overprinting the bottom water εNdd with NSW-like εNd would result in a trans-SSW ΔεNdd of 0 (e.g. close to Fig. 7d). Figure 8 shows that in order to begin deviating from conservative behaviour, irrespective of advection strength, the residence time τ must fall below a few decades. Sensitivity to τ is greatest in the range 100 − 101 years—within this range, relatively small changes in scavenging intensity can result in very significant changes in the water column εNdd profile. Within the range plotted in Fig. 8, εNdd has a limited sensitivity to advection strength apart from at very high scavenging intensity.

Fig. 8
figure 8

Parameter space showing the effects of residence time (scavenging intensity) and scaled water velocity on bottom water overprinting. Scaled water velocity is relative to the LGM-like profile, i.e. vLGM = 1 and vHS = 0.2. Contours show the εNd difference across the SSW, where − 5 represents the advected profile (i.e. no overprinting, entirely conservative behaviour) and 0 represents total overprinting of the ABW εNd. Note that neither end-member is reached in this parameter space. Results from model runs in Fig. 7b, d are plotted in green and blue respectively

Full size image

Data-model misfit can also be explored within this parameter space, as shown in Fig. 9. These diagrams again strongly implicate an approximately order of magnitude increase in scavenging rate to explain the differences between the reconstructed LGM and HS profiles. The advection strength corresponding to the lowest model-data misfit is not significantly different between the LGM and HS; however, at these points, the misfit is relatively insensitive to advection strength.

Fig. 9
figure 9

Data-model RMS misfit (εNd) plotted as a function of residence time (scavenging intensity) and scaled water velocity, as in Fig. 8, for a LGM data and b HS data (Fig. 7). The BOFS 10K site (2777 m) is excluded as this point is difficult to reconcile with all profiles

Full size image

A similar method is used to investigate the assumptions made on the NSW/SSW interface depth (see “Velocity profile”) by plotting data-model misfit as a function of interface depth and scavenging intensity (Fig. 10). The LGM data reveal a trade-off between interface depth and scavenging intensity. It is possible to get good data-model agreement across all reasonable interface depths; however, a shoaling of 1 km relative to 3.5 km requires an approximately order of magnitude increase in scavenging intensity, eroding any LGM-HS difference. The HS data also permit good data-model agreement across most reasonable interface depths and there is a weaker trade-off with scavenging intensity, unsurprising due to the simple structure implied by the HS data. These data are therefore not fundamentally incompatible with a shoaled NSW/SSW interface; however, the implication is that a shoaled interface would require significantly greater scavenging rates at the LGM, which may be difficult to justify.

Fig. 10
figure 10

Data-model RMS misfit (εNd) plotted as a function of residence time (scavenging intensity) and NSW-SSW interface depth, using the LGM-like and HS-like velocity profile scalings (justified due to the relative insensitivity to velocity, e.g. Fig. 9), for a LGM data and b HS data as in Fig. 9. The NSW/SSW interface depth was modified by simply linearly shifting the default interface up or down, using the top/bottom water velocity to fill in the blanks. Note that this simple approach may generate unrealistic velocities in the shallowest/deepest waters

Full size image

The bathymetry in the eastern North Atlantic is not flat and water column εNdd can only be recorded where a sediment-water interface exists. To test whether diagenetic remobilisation of adsorbed εNd within the sediment could alter the above results, a new diagenetic εNd profile is produced, εNddiagenetic. For this metric, IRD is assumed to be concentrated by a factor of 103 in the sediment versus the water column and a remobilisation of 1% of adsorbed Nd to produce an altered diagenetic εNd signature. Blaser et al. (2019) did not find evidence of vertical migration of remobilised Nd within sediments, so we do not consider the vertical mixing of pore fluids. The results are presented in Fig. 11a, b. The εNdads has a greater ‘inertia’ with depth than the εNdd profile (i.e. it retains signals from overlying water masses for longer), so the εNddiagenetic profile is an intermediate between the two. The best model-data fits for εNddiagenetic require weaker scavenging than εNdd (τ= 40 and 5 years for the LGM-like and HS-like profiles respectively) but the difference between these two residence times remains similar to the εNdd case, about an order of magnitude.

Fig. 11
figure 11

Modelled steady-state εNd profiles, a, b for the LGM and HS cases with a hypothetical diagenetic εNd profile (dashed) added, assuming 1% particle dissolution in the sediment and a sediment:water column concentration factor of 103; c, d for the LGM case with τ = 300 years and κD = 10−5 and 10−3 m2 s−1 respectively; e, f, g with benthic flux switched on with κD = 5×10−4 with a piston velocity of 3×10−7 m s−1 and [Nd]pf = 1000pM (q.v. Abbott et al. 2015)

Full size image

Explanations for reconstructed εNd profiles

The results of our model (Figs. 7, 8 and 9) indicate that, for IRD fluxes to drive significant non-conservative εNd behaviour, (1) the residence time of Nd in the glacial and deglaciating eastern North Atlantic would have to be on the order of decades at most and (2) an increase in scavenging intensity of approximately an order of magnitude is needed to transition from an LGM-like profile to a HS-like profile. In the following, we discuss to what extent IRD is capable of causing this behaviour and whether the model results are reasonable and applicable to the glacial eastern North Atlantic.

Conservative behaviour

The null hypothesis is that εNd is behaving conservatively, and that the changes in the εNd profile between the LGM and HSs are therefore due to the movement of water masses and/or are an inherited property obtained from upstream changes in the εNd profile. However, advection of water masses cannot explain the highly gradational LGM εNd profile (Fig. 2). The conservative interpretation of Fig. 2 would require a significant increase in the strength and depth of NSW in the eastern North Atlantic during Heinrich Stadials versus the LGM, the opposite of what other proxies suggest (e.g. McManus et al. 2004). Whilst there is evidence for significant NSW end-member εNd variability during the LGM and deglaciation (Roberts and Piotrowski 2015), this still cannot explain the partial and total absence of a SSW signature in bottom water εNd without either an expansion in NSW or SSW end-member εNd changing in lockstep. Therefore, it is difficult to find a reasonable justification for complete conservative interpretations for the observed glacial and deglacial εNd profiles at the Dreizack seamount of Blaser et al. (2019).

Diffusion

Whilst the LGM profile in Fig. 2 superficially resembles a diffusive boundary between two water masses, the (kilometre) scale of this trend cannot be reconciled with the strength of advection that acts to reinforce a sharp water mass boundary. Vertical diffusivity in the ocean is highly spatially variable, ranging from order 1 × 10−5 m2 s−1 in the open ocean to over 5 × 10−4 m2 s−1 above rough topography (Polzin et al. 1997; Toole et al. 1997). Since the εNd profile of Blaser et al. (2019) was collected from a seamount, it is feasible that diffusivity could reach the upper limits of this range. However, as demonstrated in Fig. 11c, d, even extreme values for diffusivity are unable to recreate the reconstructed LGM profile of Blaser et al. (2019). In any case, diffusion is unable to explain the homogenised water column reconstructed for HSs (Fig. 2). In the context of this problem, the only relevance of vertical diffusion is as a control on benthic Nd flux propagation into the water column (“Benthic Nd flux”).

Benthic Nd flux

It has been suggested that fluxes of pore fluid into the water column (benthic flux) may act as an additional control on bottom water εNd in the Pacific (Abbott et al. 2015, 2016). Benthic flux may also influence the here discussed core locations, but due to the stronger circulation in the Atlantic Ocean compared to the Pacific Ocean and the fact that benthic Nd fluxes would be expected to result in the most non-conservative behaviour at the sediment-water interface (in contrast to the findings of Blaser et al. (2019)), we expect it to be less relevant in the Atlantic. Additionally, whilst turbidites are prevalent in this region of the eastern North Atlantic and can result in the local overprinting of pore water εNd during early diagenesis, Blaser et al. (2019) found no evidence of vertical migration of Nd within sediments, nor any evidence of a turbidite-related benthic flux, indicating that this is only a concern for turbidite layers.

Nevertheless, the effects of benthic Nd fluxes on water column εNd have been tested to verify whether the very high scavenging rates implied by the above simulations are necessary. With a residence time of τ = 300 years (stronger scavenging intensities would overpower the effects of any benthic Nd flux) under the HS velocity profile, a εNdpore water of − 6.5 results in a εNdd profile that is within the range of the reconstructed HS profile of Blaser et al. (2019) (Fig. 11e). However, under the same parameters, switching to the LGM velocity profile results in poor model-data fit (Fig. 11f). Modifying the εNdpore water for the LGM (Fig. 11g) still fails to reproduce the reconstructed εNd gradient. In any case, we cannot suggest any justification for significantly modifying pore water properties between the LGM and HSs, at least not in a way that results in good model-data fit. The above simulations also assume extreme values for diffusivity, [Nd]pore fluid and piston velocity—under more realistic estimates, the effects of benthic Nd fluxes are significantly reduced.

Therefore, whilst our model results suggest that Nd fluxes can feasibly modify bottom water εNdd, it cannot explain the differences between the Dreizack LGM and HS εNd profiles (Fig. 2), suggesting that it is not the main cause of this non-conservative behaviour.

Reversible scavenging as a source of non-conservative εNd behaviour

In the modern ocean, the residence time of Nd is on the order of several centuries (Tachikawa et al. 2003). That Nd appears to behave quasi-conservatively in the present (van de Flierdt et al. 2016; Tachikawa et al. 2017) is therefore consistent with our model results (Fig. 7a). As discussed in “Scavenging rates”, the scavenging rates used by the model, K and K−, can be considered as a function of the material adsorption and desorption rates for IRD components, and the concentration of these components in the water column. In particular, these reactions are a function of surface area rather than volume. One can therefore simplistically assume that:

$$ K=\sum {A}_i{S}_i $$
$$ {K}^{-}=\sum {D}_i{S}_i $$

where Ai and Di are the specific adsorption and desorption rates of Nd with component i per m2 (m−2 s−1) and Si is the surface area of component i per m3 of water (m−1). For the sake of simplicity, the average geometry, volume and speciation of IRD particles is assumed to be temporally and spatially constant, which in turn results in a constant settling rate. Under these assumptions, K and K− are proportional to IRD accumulation rate.

The model prediction is that K(=K−) should increase by a factor of approximately 10 between the LGM and Heinrich Stadials. Using the above assumptions, one would therefore expect a factor of 10 increases in IRD mass flux during Heinrich Stadials versus the LGM, and this is indeed within the range observed in sediment cores (Haapaniemi et al. 2010; McManus et al. 1998).

Despite this agreement, the problem remains that the model still predicts residence times for Nd in the glacial eastern North Atlantic significantly below the residence time of Nd in the global ocean. To test whether this is realistic, it would be useful to know reasonable values for Ai and Di for specific IRD components (or even lithics in general) as described above, but to our knowledge these parameters do not yet exist. Obtaining these constraints is complicated by the fact that IRD entering the modern North Atlantic is highly heterogenous (Maslov et al. 2018) and lacks the detrital carbonate contribution from the Laurentide Ice Sheet (Broecker et al. 1992).

An additional complication is the significantly different εNd profiles that can be obtained if diagenetic contributions from adsorbed Nd are allowed (Fig. 11a, b), particularly with higher residence times when εNdd and εNdads are significantly offset (e.g. Fig. 11a). In this case, matching reconstructed εNd profiles with the model becomes a curve fitting exercise. The purpose of Fig. 11 is therefore not to suggest that the modelled diagenetic system (1000-fold particle concentration in the sediment versus the water column and 1% dissolution) is realistic, but to demonstrate that the return of Nd from the adsorbed phase to diagenetic fluids could significantly alter the archived εNdd signal.

Modern analogues

Recent studies from the Panama Basin (Grasse et al. 2017) and Bay of Bengal (Yu et al. 2017) find evidence of variability in εNd of deep waters on annual timescales that is not associated with changes in hydrography. If this is indeed the case, whilst these two cases are clearly in a very different setting to the glacial IRD belt, this provides two modern examples of significantly non-conservative behaviour in the εNd of deep waters in the modern ocean and the suggestion of short (annual to decadal) residence times in the glacial and deglaciating eastern North Atlantic becomes more believable. Both settings of the studies, the Panama Basin and Bay of Bengal, are characterised by high particle fluxes (biogenic and lithogenic). However, the mineralogy of these particles differs from that of typical Ruddiman-belt IRD and in both of these cases; the authors argue that particles are either acting as Nd sources or sinks rather than implicating reversible scavenging.

Conclusion

Bottom water εNd profiles for the LGM and Heinrich Stadials 1 and 2 in the eastern North Atlantic, obtained from sediment cores from the Dreizack Seamount (Roberts and Piotrowski 2015; Blaser et al. 2019) show patterns that are difficult to reconcile with assumptions of conservative εNd behaviour. Using a simple water column box model, we have demonstrated that reversible scavenging processes are capable of producing simulated εNd profiles that fit well with geochemical reconstructions for the LGM and Heinrich Stadials (Fig. 7), necessitating a (localised) residence time on the order of 1–10 years within the North Atlantic IRD belt, significantly lower than the centennial-scale residence times for the ocean as a whole (Tachikawa et al. 2003). Furthermore, the predicted relative increase in scavenging intensity between the LGM and Heinrich Stadials is in agreement with the relative increase in IRD intensity. Whilst it is unclear how compatible this theory is with a significantly shoaled glacial NSW/SSW interface in the eastern North Atlantic, mechanisms that do not include reversible scavenging still fail to explain the observed εNd gradients. Whilst we cannot fully rule out alternative explanations such as conservative behaviour or benthic Nd fluxes, both explanations would require physical or chemical changes in the glacial and deglaciating eastern North Atlantic that are difficult to justify. We emphasise that this is an idealised proof-of-concept study—the idealised model used in this study clearly does not incorporate the full complexity of the biogeochemical cycling of Nd in the north Atlantic and relies on many unconstrained parameters. This model nevertheless enables the broad quantification of the scavenging rates required to reproduce observations, and these scavenging rates may be physically reasonable. We therefore argue that the suggestion that IRD could cause non-conservative behaviour in εNd is worth taking seriously.

Given the importance of reconstructing past ocean circulation for understanding the role of the oceans in abrupt climate variability and the potential applications of εNd as a palaeotracer, it is essential to have a good understanding of the geochemical cycling of Nd in the oceans and whether assumptions of conservative behaviour, as appears to be the case in most of the modern ocean (albeit with the exceptions detailed in “Modern analogues”) holds true for the past as well. If it is correct that processes such as IRD deposition (and potentially others such as productivity, dust deposition and sediment plumes from rivers) can cause significant localised non-conservative behaviour in palaeotracers, then this has serious implications for the use of these tracers in palaeoceanographical reconstructions, particularly where the modification of end-member properties is possible (as is the case in the North Atlantic). We therefore emphasise the importance of multi-proxy approaches to palaeoceanography, as well as the careful consideration of the processes involved in palaeotracer systems rather than simply considering them as simple two-component systems with fixed end-members.

Key uncertainties remain regarding the use of Nd as a palaeotracer, principally a lack of fundamental chemical constraints on adsorption-desorption rates of Nd with IRD components and a poor understanding of the effects of diagenesis on the εNd records in different sedimentary phases. We nevertheless conclude that the assumption of conservative εNd behaviour in this oceanographical setting is sufficiently uncertain that it warrants further investigation.