Articles | Volume 16, issue 4
Special issue:
Model description paper
22 Feb 2023
Model description paper |  | 22 Feb 2023

Simulating marine neodymium isotope distributions using Nd v1.0 coupled to the ocean component of the FAMOUS–MOSES1 climate model: sensitivities to reversible scavenging efficiency and benthic source distributions

Suzanne Robinson, Ruza F. Ivanovic, Lauren J. Gregoire, Julia Tindall, Tina van de Flierdt, Yves Plancherel, Frerk Pöppelmeier, Kazuyo Tachikawa, and Paul J. Valdes

The neodymium (Nd) isotopic composition of seawater is a widely used ocean circulation tracer. However, uncertainty in quantifying the global ocean Nd budget, particularly constraining elusive non-conservative processes, remains a major challenge. A substantial increase in modern seawater Nd measurements from the GEOTRACES programme, coupled with recent hypotheses that a seafloor-wide benthic Nd flux to the ocean may govern global Nd isotope distributions (εNd), presents an opportunity to develop a new scheme specifically designed to test these paradigms. Here, we present the implementation of Nd isotopes (143Nd and 144Nd) into the ocean component of the FAMOUS coupled atmosphere–ocean general circulation model (Nd v1.0), a tool which can be widely used for simulating complex feedbacks between different Earth system processes on decadal to multi-millennial timescales.

Using an equilibrium pre-industrial simulation tuned to represent the large-scale Atlantic Ocean circulation, we perform a series of sensitivity tests evaluating the new Nd isotope scheme. We investigate how Nd source and sink and cycling parameters govern global marine εNd distributions and provide an updated compilation of 6048 Nd concentrations and 3278 εNd measurements to assess model performance. Our findings support the notions that reversible scavenging is a key process for enhancing the Atlantic–Pacific basinal εNd gradient and is capable of driving the observed increase in Nd concentration along the global circulation pathway. A benthic flux represents a major source of Nd to the deep ocean. However, model–data disparities in the North Pacific highlight that under a uniform benthic flux, the source of εNd from seafloor sediments is too non-radiogenic in our model to be able to accurately represent seawater measurements. Additionally, model–data mismatch in the northern North Atlantic alludes to the possibility of preferential contributions from “reactive” non-radiogenic detrital sediments.

The new Nd isotope scheme forms an excellent tool for exploring global marine Nd cycling and the interplay between climatic and oceanographic conditions under both modern and palaeoceanographic contexts.

1 Introduction

The neodymium (Nd) isotope composition of seawater shows a clear provinciality between different ocean basins and is often used as a water mass provenance tracer (e.g. Frank, 2002; Goldstein and Hemming, 2003). The measured 143Nd/144Nd ratio is denoted relative to the bulk Earth standard:

(1) ε Nd = 143 Nd / 144 Nd sample 143 Nd / 144 Nd CHUR - 1 × 10 4 ,

where (143Nd/144Nd)CHUR relates to the Chondritic Uniform Reservoir (CHUR; 0.512638: Jacobsen and Wasserburg, 1980). Distinct variations in the Nd isotope signal of water masses originate from different continental regions and their isotopic fingerprints and subsequent influence by ocean circulation, water mass mixing, and particle cycling, as well as interaction with sediments (e.g. Tachikawa et al., 2017; van de Flierdt et al., 2016 for recent reviews). Neodymium in the deep ocean has an estimated mean ocean residence time between 360 and 800 years (Arsouze et al., 2009; Rempfer et al., 2011; Gu et al., 2019; Pöppelmeier et al., 2020a, 2021b), which is shorter than the global overturning of the deep ocean. Unlike other tracers of ocean circulation (e.g. δ13C, Δ14C), Nd is not fractionated by marine biological cycling, giving rise to its promise as a carbon-cycle-independent ocean circulation tracer (Blaser et al., 2019a). Yet, a fundamental caveat in the application of εNd as a reliable oceanographic tracer is that knowledge of the Nd fluxes in and out of the ocean remains incomplete, and while internal cycling must occur, the efficiency of it with different particle types is poorly constrained (Abbott et al., 2015a; Haley et al., 2017; van de Flierdt et al., 2016).

Numerical models are useful tools for investigating Nd cycling since they can specify the processes that govern the spatiotemporal variability in Nd isotope distributions in the ocean. Neodymium isotopes have been simulated in a range of different modelling studies testing specific hypotheses relating to Nd fluxes and thermohaline redistribution (Pöppelmeier et al., 2020a; Arsouze et al., 2009; Rempfer et al., 2011; Siddall et al., 2008; Jones et al., 2008; Roberts et al., 2010; Tachikawa et al., 2003; Ayache et al., 2016; Gu et al., 2019; Oka et al., 2021, 2009; Pöppelmeier et al., 2022; Ayache et al., 2023; Pasquier et al., 2022; Du et al., 2020). However, recent work suggests that a seafloor-wide benthic flux, resulting from early diagenetic reactions, may dominate the marine Nd cycle (Haley et al., 2017; Abbott, 2019; Abbott et al., 2015a, b, 2019; Du et al., 2016). These observations, alongside an ever-growing body of high-quality and highly resolved measurements of dissolved seawater Nd concentrations ([Nd]) and εNd from the GEOTRACES programme (GEOTRACES Intermediate Data Product Group, 2021), present an opportunity to re-evaluate, revise, and explore constraints on the marine Nd cycle.

Initially, the predominant lithogenic fluxes of Nd to the ocean were believed to be mostly at the surface (aeolian dust and riverine fluxes; Goldstein et al., 1984; Goldstein and Jacobsen, 1987). Early modelling studies applying surface fluxes reproduced reasonable εNd in the North Atlantic (Bertram and Elderfield, 1993; Tachikawa et al., 1999). However, considering only dust and river fluxes alone led to an unrealistic calculated residence time of Nd in seawater on the order of 5000 years (Bertram and Elderfield, 1993; Jeandel et al., 1995). It was then found that considering dust and river surface inputs alone failed to balance both [Nd] and εNd, thus indicating a “missing source” of Nd to the ocean that accounted for ≈90 % of the Nd flux to the ocean (Tachikawa et al., 2003). This led to a new hypothesis relating to other Nd sources to seawater that could account for this missing source, including submarine groundwater discharge (SGD) (Johannesson and Burdige, 2007) and exchange of Nd between seawater and sediment deposited on the continental margins (Lacan and Jeandel, 2005). The term “boundary exchange” was then coined to describe significant modification of Nd isotopic composition by the co-occurrence of Nd release from sediment and boundary scavenging without substantially changing [Nd] (Lacan and Jeandel, 2005). Arsouze et al. (2007) simulated realistic global εNd distributions using sedimentary fluxes constrained to the continental margins as the only source–sink term, and since then particle–seawater exchange along the margins has represented the major flux of Nd to seawater in recent global Nd isotope models (Arsouze et al., 2009; Rempfer et al., 2011; Gu et al., 2019).

Nonetheless, seafloor–sediment interactions along the continental margins alone cannot fully reconcile the global marine Nd cycle. Specifically, it cannot explain the observed vertical profiles of [Nd], which are decoupled from εNd (i.e. the “Nd paradox”; Goldstein and Hemming, 2003), with low concentrations near the surface increasing with depth. This is a common characteristic of isotopes and elements (e.g. thorium) that are reversibly scavenged (i.e. where the element is scavenged onto sinking particles at the surface and is subsequently remineralised in the deep ocean) (Bertram and Elderfield, 1993; Bacon and Anderson, 1982). Siddall et al. (2008) first addressed a numerical hypothesis that the Nd paradox can be explained by a combination of lateral advection and reversible scavenging by applying the reversible scavenging model pioneered by Bacon and Anderson (1982) to Nd cycling. In their study, both [Nd] and εNd were modelled simultaneously and explicitly to explore internal cycling of Nd in the deep ocean. Their findings demonstrated that scavenging and remineralisation processes are important active components in the marine cycling of Nd, driving the increase of [Nd] with depth but still allowing εNd to act as an effective water mass tracer.

Although inclusion of reversible scavenging can explain aspects of marine Nd cycling, the use of oversimplified fixed surface [Nd] and εNd boundary conditions in the model by Siddall et al. (2008) limited what could be determined about the full marine cycling of Nd and hence the Nd paradox. The most comprehensive Nd-isotope-enabled ocean models to date now explicitly represent and quantify a wider range of distinct Nd fluxes that are both external and internal to the marine realm. For example, Arsouze et al. (2009) used a fully prognostic coupled dynamic and biogeochemical model to simulate [Nd] and εNd, considering dust fluxes, dissolved riverine sources, sediment sources along the continental margins, and reversible scavenging. In their study, a sediment source from the continental margins represented the major source of Nd to seawater (≈95 % of the total source). Rempfer et al. (2011) continued this work, undertaking a more detailed and comprehensive investigation of Nd sources and particle scavenging using a coarse-resolution intermediate-complexity model (Bern3D ocean model) and extensive sensitivity experiments. This later scheme was then closely followed by Gu et al. (2019) for the implementation of Nd isotopes in the ocean component of a more comprehensive Earth system model (ESM, specifically the Community Earth System Model, CESM1.3) to explore the changes to end-member εNd signatures in response to ocean circulation and climate changes in detail. Overall, these comprehensive models, capable of quantifying the major sources implicated in marine Nd cycling, indicated that dust and river fluxes were important for representing [Nd] and εNd distributions in the surface but that the main flux of Nd to seawater is via a sediment source operating along the continental margins.

Recent pore fluid concentration profiles measured on the Oregon margin in the Pacific Ocean indicate that there may be a benthic flux of Nd from sedimentary pore fluids, presenting a potential major seafloor-wide (i.e. no longer limited to the continental margins) source of Nd to seawater (Abbott et al., 2015b, a); however, there are observations of low [Nd] from southeastern equatorial Pacific pore waters, implying possibly heterogeneous benthic flux in relation to deep-sea sediment composition (Paul et al., 2019). Additional measurements from the Tasman Sea inferred the likely presence of a benthic source of similar magnitude to that estimated for the North Pacific, which may indicate that regions with dominantly calcareous sediment also contribute a significant benthic source of Nd (Abbott, 2019). Evidence of this previously overlooked benthic sedimentary source of Nd has led to a shifting paradigm that contends that the dominant addition of Nd to the ocean is from a diffuse sedimentary source at depth, rather than surface point sources from rivers and dust and the shallow to intermediate continental margins (Haley et al., 2017).

Simple box models have been employed to investigate, to a first order, the non-conservative effects from a benthic flux (Du et al., 2016, 2018, 2020; Haley et al., 2017; Pöppelmeier et al., 2020b), suggesting overprinting of deep-water mass εNd is linked to benthic flux exposure time and the difference between the Nd isotope composition of the benthic flux and the bottom water (Abbott et al., 2015a; Du et al., 2018). However, these models lack comprehensive descriptions of both the marine Nd cycle and of physical ocean circulation and climate interactions, limiting a clear interpretation of precisely how (and under what physical or environmental conditions) the benthic flux may determine global marine Nd distributions. Applying an intermediate complexity model, Pöppelmeier et al. (2021a) investigated the benthic flux hypothesis in more detail by updating the Nd-isotope-enabled Bern3D model (Rempfer et al., 2011) to represent recent observations that indicate a Nd flux from bottom waters could occur across the entire seafloor. This was done by removing the depth limitation of the sediment source (previously 3 km) and invoking a constant benthic flux that escapes from all sediment–water interfaces. The scheme was further extended by revising key source–sink parameterisations for subsequent investigation of non-conservative Nd isotope behaviour under different ocean circulation states (Pöppelmeier et al., 2022). The authors demonstrate substantial non-conservative effects occur even under strong circulation regimes with low benthic flux exposure times and are not strictly limited to the deep ocean. This work highlights the importance of downward vertical fluxes via reversible scavenging alongside the benthic flux to describe non-conservative marine Nd isotope behaviour. Nonetheless, the low horizontal resolution of the intermediate-complexity model limits full resolution of key circulation patterns such as distinct deep-water formation in the Labrador Sea and in the Nordic Seas, inhibiting the capability of the scheme to fully capture and explore water mass end-member εNd distributions.

Thus, despite substantial progress to explicitly describe seawater Nd budgets, it is apparent that outstanding questions remain, alongside divergent lines of argument amongst subject specialists, limiting a full, quantitative description of marine Nd cycling. The decoupled nature of marine [Nd] and εNd, which is yet to be fully understood, and new emerging observations (e.g. Stichel et al., 2020; Lagarde et al., 2020; Paffrath et al., 2021; Deng et al., 2022) emphasise the critical need to progress our understanding of modern marine Nd cycling in the light of continued use of εNd as a valuable tracer of ocean circulation. In this context, Nd-isotope-enabled ocean models remain an effective way to understand seawater measurements and to progress with testing and constraining processes governing marine Nd cycling, which can then feedback key information to the wider GEOTRACES and ocean tracer modelling community.

With this is mind, there is a current gap in the toolkit for modelling marine Nd cycling between the class of high-complexity, state-of-the-art ESMs (e.g. Gu et al., 2019) and the more efficient intermediate-complexity models (e.g. Rempfer et al., 2011; Pöppelmeier et al., 2020a). To bridge this gap, there is a need for a model with the full complexity of an atmosphere–ocean general circulation model (AOGCM), allowing the exploration of how Nd isotopes vary under changing climate conditions (including extensive palaeoceanographic applications), that is also capable of running very quickly to facilitate efficient parameter space exploration, performance optimisation and long integrations. Owing to a winning combination of speed and complexity, the FAMOUS general circulation model (GCM) fills this gap and is well suited to exploring Earth system interactions (Smith et al., 2008; Jones et al., 2005; Smith, 2012; Jones et al., 2008).

Here we present the new global marine Nd isotope scheme (Nd v1.0) implemented in FAMOUS. We utilise the new sediment εNd maps from Robinson et al. (2021) as boundary conditions for a mobile global sediment Nd flux, with the end goal of further constraining the major sources, sinks, and cycling of Nd isotopes and exploring instances of non-conservative behaviour related to changes in Nd source distributions and scavenging processes. We develop sensitivity experiments (Sect. 3) to isolate the physical effects of varying two key parameters that detail major Nd fluxes and cycling to the global ocean, verifying foremost that the new isotope scheme is responding to Nd source and sink and biogeochemical processes as expected and contextualising how and where reversible scavenging processes and the benthic flux can govern marine [Nd] and εNd distributions. Finally (Sects. 3 and 4), we evaluate overall model performance through comparison with modern measurements and assess the model's ability to simulate observed spatial and vertical gradients between ocean basins, encompassing underling structural uncertainty and model bias.

2 Methods

2.1 Model description

We use the FAMOUS GCM (Smith et al., 2008; Jones et al., 2005; Smith, 2012), a fast coupled AOGCM derived from the Met Office's Hadley Centre Coupled Model V3 (HadCM3) AOGCM (Gordon et al., 2000). The atmospheric component of FAMOUS is based on quasi-hydrostatic primitive equations and has a horizontal resolution of 5 latitude by 7.5 longitude, 11 vertical levels on a hybrid sigma-pressure coordinate system and a 1 h time step. The ocean component is a rigid-lid model, with a horizontal resolution of 2.5 latitude by 3.75 longitude and 20 vertical levels, spaced unequally in thickness from 10 m at the near-surface ocean to over 600 m at depth, and a 12 h time step. The ocean and atmosphere are coupled once per day.

FAMOUS's parameterisations of physical and dynamical processes are almost identical to those of HadCM3, its parent model, but it has approximately half the spatial resolution and a longer time step, allowing it to run 10 times faster. Thus, FAMOUS achieves a current speed of up to 650 model years per wall clock day on 16 processors, making it ideal for running large ensembles (Gregoire et al., 2011, 2016), more bespoke sensitivity studies (Smith and Gregory, 2009; Gregoire et al., 2015), and multi-millennial simulations, e.g. to examine ice sheet behaviours (Gregoire et al., 2012, 2016, 2015), ocean drifts (Dentith et al., 2019), and ocean biogeochemistry (Dentith et al., 2020). FAMOUS is calibrated to the performance of HadCM3, taking the philosophy that this is the most appropriate evaluation target and that it is unrealistic to expect the lower-resolution, lower-complexity model to outperform its parent model (Valdes et al., 2017; Smith and Gregory, 2009).

We added Nd isotopes (143Nd and 144Nd) as optional passive tracers into the ocean component of FAMOUS, using a version of the model with the Met Office Surface Exchange Scheme (MOSES) version 1 (Cox et al., 1999: FAMOUS–MOSES1). It was a pragmatic choice to avoid the more recent FAMOUS–MOSES2.2 configuration (Essery and Clark, 2003; Valdes et al., 2017; Williams et al., 2013; Essery et al., 2001) because evaluation of our new Nd scheme would be hindered by the collapsed Atlantic Ocean convection and strong deep Pacific meridional overturning circulation (MOC) produced by FAMOUS–MOSES2.2 under pre-industrial boundary conditions (see Dentith et al., 2019). Nonetheless, the Nd isotope code implementation presented here is directly transferable to other versions of the UK Met Office Unified Model (UM) version 4.5, including HadCM3/L or FAMOUS–MOSES2.2.

2.2 A new reference simulation for this study

The use of Nd isotopes as a water provenance tracer comes from measurements of distinct εNd signatures in different water masses. This is perhaps best demonstrated in the south Atlantic “zig-zag” depth profiles, where εNd displays large heterogeneity and distinguishes southward-flowing North Atlantic Deep Water (NADW) from northward-flowing Antarctic Intermediate Water (AAIW) and Antarctic Bottom Water (AABW) (Goldstein and Hemming, 2003; Jeandel, 1993). As such, it is desirable to implement the Nd isotope scheme in a version of FAMOUS that best positions these basinal water masses in the correct locations.

The standard pre-industrial FAMOUS setup (XFHCC; Smith, 2012) has certain known limitations, including over-ventilated abyssal Atlantic waters characterised by strong, over-deep NADW formation, with “North Atlantic” convection occurring only in the Norwegian and Greenland seas (there is no deep water formation in the Labrador Sea), and insufficient Atlantic sector AABW formation (Dentith et al., 2019; Smith, 2012). This known physical bias would dominate simulated Nd distributions, thus hampering validation of the new Nd isotope scheme against modern measurements. To mitigate this, we chose to employ a new reference simulation with improved basin-scale physical ocean circulation. This simulation was obtained from a perturbed parameter ensemble varying 13 physical tuning parameters (see Sect. S1 in the Supplement and Table S1 for brief description and for a list of perturbed parameters from this multi-sweep ensemble of FAMOUS MOSES1: Smith, 1990, 1993; Crossley and Roberts, 1995). We screened the resulting 549 simulations based on a set of pre-defined targets, with a particular focus on Atlantic Ocean structure and water mass composition. Specifically, we sought a simulation with appropriate modern Atlantic meridional overturning circulation (AMOC) strength (14–19 Sverdrup, Sv, where 1 Sv = 106 m3 s−1) (Frajka-Williams et al., 2019), AMOC structure (Talley et al., 2011), regions of AMOC convection as indicated by mixed-layer depth (specifically in both the Labrador Sea and the Nordic Seas as these represent key regions for NADW formation and the resultant end-member Nd isotope compositions; Lambelet et al., 2016), depth of maximum overturning (≈1000 m), and presence of AABW in the abyssal Atlantic (Frajka-Williams et al., 2019; Talley et al., 2011; Kuhlbrodt et al., 2007; Ferreira et al., 2018). Furthermore, we assessed the capabilities of the simulation to represent appropriate water mass structure in the Pacific (Talley et al., 2011). Through this approach, we identified four possible candidate simulations from the large ensemble as the basis for the new Nd isotope scheme, XPDAA, XPDAB, XPDAC, and XPDEA, here denoted by their unique five-letter Met Office UM identifier. See Table S2 in the Supplement for initial boundary conditions and Figs. S1–S4 for the Atlantic meridional stream function and mixed-layer depth for the four experiments.

These simulations were integrated for a further 5000 years to ensure adequate spin-up of the physical circulation. They were then assessed for their ability to simulate global modern observations of salinity and temperature (compared to the NOAA World Ocean Atlas salinity and temperature database; Locarnini et al., 2019; Zweng et al., 2019) as an objective and quantitative basis for selecting the best-performing parameter configuration to be our control (Fig. 1).

Figure 1Taylor diagrams summarising the performance of the four control candidate simulations (XPDAA, XPDAB, XPDAC, XPDEA) in terms of their correlation, centred root-mean-square error, and ratio of variance to the NOAA World Ocean Atlas (a) salinity and (b) temperature databases (Zweng et al., 2019). The simulations were selected from a large FAMOUS pre-industrial perturbed parameter ensemble following the circulation performance screening described in Sect. 2.2. XPCQX represents the initial experiment, upon which the four control candidates aim to improve.

From this analysis, XPDAA returned the lowest root-mean-square error (RMSE) (and hence best performance using this metric; Fig. 1) for simulating both salinity (RMSE of 0.93 PSU) and temperature (RMSE of 2.21 C) and was thus used as our control simulation (Fig. 2). For comparison, HadCM3, the parent model of FAMOUS (Gordon et al., 2000), and the Model for Interdisciplinary Research on Climate 4m (MIROC4m) AOGCM (Sherriff-Tadano et al., 2021) have a RMSE of 1.06 PSU and 1.47 C and 0.59 PSU and 1.42 C, respectively.

Figure 2(a) Salinity and (b) temperature profiles for the control simulation (centennial mean from final 100 years of a 5000-year simulation) along a transect crossing the Pacific–South Atlantic oceans.


The simulated steady-state AMOC strength in FAMOUS under fixed pre-industrial boundary conditions is approximately 17 Sv (Fig. 3a), which we consider to be in excellent agreement with direct modern AMOC observations from the RAPID AMOC array at 26.5 N of 17.2 Sv from April 2004 to October 2012, and the depth of maximum meridional stream function at 26.5 N is around 800 m (Fig. 3b), slightly shallower than RAPID observations of 900–1100 m (McCarthy et al., 2015; Sinha et al., 2018). In terms of the Atlantic circulation structure, the overturning cell of NADW descends to depths of 3000 m as it bridges into the South Atlantic, and AABW fills the bottom of the Atlantic basin with southern-sourced waters up to 20 N. For the Nd isotope implementation, all sensitivity studies and model tuning described subsequently are based upon this simulation (XPDAA).

Figure 3Atlantic meridional streamfunction for a 5000-year simulation using the control simulation (XPDAA), showing (a) the maximum streamfunction in time series (dashed red line indicates RAPID-AMOC 2004–2012 averaged AMOC strength at 26.5 N of 17.2 Sv; McCarthy et al., 2015) and (b) the zonal integration calculated from the last 100 years.

2.3 Neodymium isotope implementation in FAMOUS

In our simulated Nd isotope scheme (Nd v1.0), we represent three different global sources of Nd into seawater: aeolian dust flux, dissolved riverine input, and seafloor sediment. Neodymium (Nd) here refers to the sum of 143Nd and 144Nd, and the isotopic ratio (IR) relates to the ratio of 143Nd to 144Nd, as shown in the following equations:


By rearranging Eq. (2) and using the isotopic ratio (IR: Eq. 3), individual fluxes of each isotope can be calculated (Eqs. 4 and 5) using information about Nd fluxes from each specific source and their associated IR distributions:


Neodymium isotopes are thus simulated and transported individually and independently as two separate tracers in our scheme, explicitly resolving the concentration and distribution of each Nd isotope, allowing for [Nd] and εNd to be calculated offline from the model output.

The implementation of each source and sink term is described in detail in the following sections. To summarise these different components, Nd sources from aeolian dust fluxes and dissolved riverine input enter the ocean only via the uppermost near-surface ocean layer of the model. Seafloor sedimentary fluxes, an umbrella term that refers to a multitude of processes encompassing input from sediment deposited on the continental margins (Lacan and Jeandel, 2005), submarine groundwater discharge (Johannesson and Burdige, 2007), and an abyssal (i.e. below 3 km) benthic flux released from pore waters (Abbott et al., 2015a), are simulated via a combination of a sedimentary source applied across sediment–water interfaces together with a separate sink occurring via particle scavenging. Removal of Nd from the ocean model occurs when Nd scavenged or adsorbed onto sinking biogenic and lithogenic (dust) particles reaches the seafloor via vertical fluxes and undergoes sedimentation, acting as an active sink by removing the particle-associated Nd from the ocean.

In the model numerics, fluxes of each Nd isotope into the ocean (kg m−3 s−1) are multiplied by a factor of 1018. This technique is a common heuristic to assist human comprehension. It minimises the mathematical error associated with carrying small numbers (e.g. when manually checking the code inputs and outputs). Accordingly, concentrations of each Nd isotope in model units are in 1018 kg m−3. A full list of all variables described in the text and their abbreviations is given in Table 1.

Table 1Nd scheme model parameters, abbreviations, fixed model parameter values, and units.

Download Print Version | Download XLSX

2.3.1 Dust source

Surface dust deposition (Fdust) is prescribed in the model from the annual mean dust deposition for the pre-industrial era as simulated by the atmosphere component of the Hadley Centre Global Environmental Model version 2 (HadGEM2-A) GCM (Collins et al., 2011) (Fig. 4a). The dust deposition scheme (described by Woodward, 2011) has been shown to be in generally good agreement with observations, with concentrations in the Atlantic simulated well across the whole of the Saharan dust plume; however, some discrepancies occur, including an overestimation at some Pacific sites during spring (Collins et al., 2011). The simulation of pre-industrial climate conditions within HadGEM2 are described by Hopcroft and Valdes (2015), and the dust results specifically are described in full by Hopcroft et al. (2015). Based on these simulated dust fluxes, we apply an Nd source per volume (Sdust: g m−3 yr−1) in the uppermost layer of the ocean model, assuming a global mean concentration of Nd in dust Cdust (Cdust=20µg g−1) (Goldstein et al., 1984; Grousset et al., 1988, 1998), from which only a certain fraction βdust (βdust=0.02: Greaves et al., 1994) dissolves in seawater.

(6) S dust i , k = F dust i , 1 × C dust × β dust d z i , 1

Here i,k represents the horizontal and vertical indexing of model grid cell and dz is the grid cell's thickness (10 m in the uppermost surface layer where k=1). The total global flux of Nd from surface aeolian dust deposition to seawater (fdust) is 0.33 Gg yr−1. Although we use an updated dust deposition field compared to previous studies, prescribed Cdust and βdust are broadly consistent with earlier Nd isotope schemes, providing a comparable total Nd dust source (0.1–0.5 Gg yr−1; Arsouze et al., 2009; Gu et al., 2019; Pöppelmeier et al., 2020b; Rempfer et al., 2011; Tachikawa et al., 2003).

Figure 4(a) Annual dust deposition taken from the pre-industrial annual mean dust deposition simulated by HadGEM2-A GCM (Hopcroft et al., 2015). (b) The εNd signal from dust deposition following Tachikawa et al. (2003) and updated with information from Mahowald et al. (2006) and Blanchet (2019).

For the Nd isotope compositions of the dust flux, we started with the first-order estimate by Tachikawa et al. (2003), as follows: North Atlantic >50 N: εNd=-15; Atlantic <50 N: εNd=-12; North Pacific >44 N: εNd=-5; Indo-Pacific <44 N: εNd=-7; and remainder:
εNd=-8. This was revised with additional constraints, accounting for the dust plume expansions as reported by Mahowald et al. (2006), in combination with detrital Nd isotope data by Blanchet (2019) (Fig. 4b). For this, the mean detrital isotopic signature of a source region (e.g. the Sahara) was calculated and then applied to the region where this dust is deposited over the ocean based on the plume expansion by Mahowald et al. (2006) (for the example of the Sahara this is the North Atlantic). Regions in between major dust plumes were linearly interpolated. For the polar regions no constraints are available, and a default value has been assigned instead, but the lack of dust deposition in these regions ensures that this has no impact on the Nd cycle.

2.3.2 Dissolved riverine source

To represent the Nd source from dissolved river fluxes, we used the river outflow to the ocean simulated by FAMOUS (RIVER) as our river water discharge (g m−2 yr−1) and combined this outflow with both the Nd concentration (Criver; µg g−1) and isotopic concentration (used to calculate the flux from each Nd isotope using Eqs. 4 and 5) of river water dissolved material, as estimated by Goldstein and Jacobsen (1987; see Table 3). All river source Nd fields are shown in Fig. 5.

Figure 5(a) Simulated river outflow (RIVER) in FAMOUS, (b) major river εNd, (c) major river [Nd], and (d) the resulting riverine Nd source. The coastal grids in (b) and (c) are prescribed following average [Nd] and εNd estimates of dissolved river runoff to each of the oceans by Goldstein and Jacobsen (1987; see Table 3).

River outflow in FAMOUS is based on a routing scheme that instantaneously delivers terrestrial runoff (precipitation  [evaporation + soil moisture]) from the location that precipitation falls to the designated coastal grid cell that the runoff would reach due to river routing (i.e. the river mouth). Relating the riverine Nd flux to the model's prognostic river discharge allows the Nd river source to respond to different climatic conditions, making it a more dynamic and predictive tool for examining the impact of changes in the global hydrological balance, such as wetting and drying events or shifts in the monsoons. Essentially, river outflow plumes are “tagged” with the estimated εNd (which is provided as an input map along the coasts; Fig. 5b) according to the model's projected water discharge at that location and Criver (also provided as an input map along the coasts; Fig. 5c). For palaeoclimate or future climate applications where river routing is significantly different from today, the input maps controlling Criver tagging at the coast would need to be updated to reflect Nd dissolution from the modified fluvial pathway over land (the model should predict the rest).

Estuaries are important biogeochemical reactors of rare-earth elements (REEs), with sea salt driving flocculation of river-dissolved organic matter, which results in estuarine REE removal (Elderfield et al., 1990). This removal is known to be important in balancing the marine REE budgets; however, the complex processes involved are still not fully constrained (Elderfield et al., 1990). Rousseau et al. (2015) summarised published observations of Nd removal (%) in estuaries from observations of estuarine dissolved Nd dynamics (Rousseau et al., 2015; their Table 1). Published values (n=17) range from 40 % in Tamar neaps (Elderfield et al., 1990) to 97 % in the Amazon (Sholkovitz, 1993), with a mean Nd removal of 71 ± 16 % SD. Based upon this, and parallel to previous model schemes (e.g. Rempfer et al., 2011; Gu et al., 2019; Arsouze et al., 2009), we assume that 70 % of dissolved Nd from river systems are removed in estuaries (i.e. γriver=0.7).

The dissolved riverine source per unit volume of Nd (Sriver: g m−3 yr−1) in the uppermost layer of the ocean is therefore calculated as follows:

(7) S river i , k = RIVER i , 1 × C river × 1.0 - γ river d z i , 1 .

The total global flux of river-sourced dissolved Nd to seawater in the model (friver) is 0.44 Gg yr−1. Previous Nd isotope schemes have applied either fixed annual mean continental river discharge estimates from Goldstein and Jacobsen (1987) and Dai and Trenberth (2002), as applied in the Bern3D Nd isotope scheme (Rempfer et al., 2011; Pöppelmeier et al., 2021a), or, similar to this study, using the model's own river routing and discharge schemes, as with NEMO-OPA (Arsouze et al., 2009) and CESM1 (Gu et al., 2019). Our estimated global total dissolved riverine Nd source to the oceans sits within previous model estimates (0.26–1.7 Gg yr−1; Arsouze et al., 2009; Gu et al., 2019; Pöppelmeier et al., 2020b; Rempfer et al., 2011; Tachikawa et al., 2003).

There is a larger range in estimated riverine Nd flux to the ocean relative to the estimated dust flux ranges. It should be noted that the largest simulated Nd river source amongst these studies (1.7 Gg yr−1; Pöppelmeier et al., 2020b) in the updated Bern3D Nd isotope scheme applied a river scaling factor, used as a tuning parameter and based on findings by Rousseau et al. (2015), who suggest a globally significant release of Nd to seawater by dissolution of river-sourced lithogenic suspended sediments grounded upon observations in the Amazon estuary. Our model does not attempt to fully resolve all complex estuarine processes, and in this study we chose to represent the dissolved riverine flux as a single source to seawater. Furthermore, our sediment Nd source to the ocean (described Sect. 2.3.3), which occurs across sediment–water interfaces, utilises the continental margin and seafloor εNd distribution maps by Robinson et al. (2021), thus using the most recent compilation of published global observations of Nd isotope compositions of river sediment samples deposited on the continental shelf and slopes (alongside geological outcrops and marine sediment samples). It is therefore likely that this margin source encompasses (at least in part) the Nd isotope fingerprint from a river particle flux.

2.3.3 Continental margin and seafloor sediment source

The sediment source describes the flux of Nd into seawater entering each model grid cell adjacent to sediment. This source is not restricted to the uppermost surface layers and can be implemented across all vertical and horizonal sediment–water interfaces, as is the case for all experiments described in this study. As such, the sediment source represents (1) input from sediment deposited on the continental margins (Lacan and Jeandel, 2005), (2) submarine groundwater discharge, as suggested by Johannesson and Burdige (2007), which releases Nd to seawater via discharge of fresh groundwater to coastal seas and is mainly limited to the upper 200 m, and (3) a benthic flux, which specifically refers to a transfer of Nd from sediment pore water to seawater resulting from early diagenetic reactions and is not depth limited (Abbott et al., 2015b, a; Haley et al., 2017).

The total sediment Nd source per unit volume (Ssed: g m−3 yr−1) into any given ocean grid cell is dependent on the fraction of the surface area of that cell that is in contact with sediment. Similar to previous schemes (e.g. Gu et al., 2019; Pöppelmeier et al., 2020b; Rempfer et al., 2011), the total globally integrated sediment-associated Nd source to the ocean (fsed: g Nd yr−1) is used as a tuning parameter.

(8) S sed i , k = f sed × A i , k A total × 1 V i , k ,

where A(i,k) is the total area of the sediment surface in contact with seawater per ocean grid cell (m2), Atotal is the total global area of the sediment surface where a sediment source occurs (m2), and V(i,k) is the volume of water per ocean grid cell (m3).

There is still no true consensus on whether (and how) to apply spatial variability to sediment Nd sources, as reflected in the way previous schemes have adopted different approaches. Arsouze et al. (2009) limited the sediment flux to the upper 3000 m of the ocean, imposing a depth-scaling factor and considering estimated [Nd] distributions across the continental margins. Both Rempfer et al. (2011) and Gu et al. (2019) simplified this method by applying spatially uniform sediment Nd fluxes, also limited to the upper 3000 m. In more recent work, Pöppelmeier et al. (2020b) removed the depth limitation (as we have done), and incorporated a geographically varying scaling factor to try to capture more localised features through increased benthic fluxes (Abbott et al., 2015b; Blaser et al., 2020; Grenier et al., 2013; Lacan and Jeandel, 2005; Rahlf et al., 2020). Pasquier et al. (2022) presented the first inverse model of global marine cycling of Nd, and the strength of the sediment Nd flux to the ocean was imposed with an exponential depth function, resembling eddy kinetic energy and particulate organic matter fluxes, which are characteristically larger near the surface and in coastal regions.

For our scheme, we do not assume spatial variations in the sediment source of Nd (fsed), essentially avoiding making explicit inferences on the nature of the sedimentary Nd source. It has been proposed that preferential mobilisation of certain components of the sediment drive spatial variations in sediment fluxes (e.g. Abbott et al., 2015a; Wilson et al., 2013; Du et al., 2016; Abbott, 2019; Abbott et al., 2019) and that both detrital and authigenic phases likely exchange Nd within pore water during early diagenesis (Du et al., 2016; Blaser et al., 2019b). However, the elusiveness of marine Nd cycling alongside our limited knowledge of the specific mechanisms controlling sediment–water Nd interactions mean that determining generalisable rules for where and under what conditions (e.g. redox environments or fresh labile detrital material) preferential mobilisation may occur is unknown and challenging to resolve. Therefore, in accordance with Pöppelmeier et al. (2020a), Du et al. (2020), Gu et al. (2019), and Rempfer et al. (2011), we adopt a constant detrital sediment flux as a first-order approximation. In fact, we contend that applying this simpler method, as opposed to constructing a more complex source term that is arguably just as arbitrary (given the uncertainty in Nd cycling), allows for a more explicit quantification of differences between observed and simulated Nd distributions. As such, without overfitting our model, we allow for the clearest indication of those parts of the system that are well understood (and represented) and those which prove deficient. Under this framework, we may separate out and test the effect of many of the major Nd sources and sinks with dedicated sensitivity simulations, including the possibility in future work of incrementally modifying the sediment source distributions to increase the complexity of the scheme and assess the impact of our various assumptions.

The isotopic ratio of the sediment Nd flux to seawater is prescribed using the recently updated global gridded map of bulk detrital εNd at the continental margins and seafloor (Robinson et al., 2021, Fig. 6a). Using Eqs. (4) and (5), fluxes of each Nd isotope are calculated from this condition and bi-linearly regridded to the model's native resolution (Fig. 6b).

Figure 6(a) Map of the global εNd distributions at the sediment–ocean interface from Robinson et al. (2021) and (b) as used as a model input in this study (bi-linearly regridded from (a) onto the coarser FAMOUS ocean grid).

Pasquier et al. (2022) first applied the sedimentary εNd map from Robinson et al. (2021) in their recent global marine Nd isotope scheme, imposing positive (i.e. radiogenic) modifications to the Pacific sedimentary εNd and using a reactivity scaling factor (linked to sediment lability) that favours more extreme εNd signals. Here, we again adopt a simpler approach, imposing the unmodified sediment εNd distributions from Robinson at al. (2021), allowing for a presentation of our new scheme based on what we confidently know about Nd cycling without the complication of over-conditioning our model inputs. Future research may then use this work as a foundation to robustly explore different choices for model inputs (i.e. boundary conditions) to revisit these fundamental questions about Nd sediment source.

2.3.4 Internal cycling via reversible scavenging

Vertical cycling and removal of Nd from the water column via sinking particles (“reversible scavenging”) is parameterised using the same approach as previous Nd isotope implementations (Siddall et al., 2008; Arsouze et al., 2009; Rempfer et al., 2011; Gu et al., 2019; Pöppelmeier et al., 2020a; Oka et al., 2021). Based on the original scheme of Bacon and Anderson (1982), the method captures the physical process of absorption or incorporation and desorption or dissolution of Nd onto particle surfaces in seawater. The scheme assumes that particle-associated Nd is in dynamic equilibrium with falling particles throughout the water column, with continuous exchange between the particle and dissolved pools. This process redistributes Nd within the water column, acting as a net sink of dissolved Nd at shallower depths as they adsorb onto particle surfaces and a net source at greater depths where dissolution of particles releases dissolved Nd back to seawater. The scavenging is the only sink of Nd in our model; Nd associated with particles (particulate organic carbon, POC; calcium carbonate, CaCO3; opal; dust) reaching the seafloor is removed from the system, operating under a steady-state assumption that acts to balance the three external input sources (marine sediments, dust, and rivers). Previous studies (e.g. Siddall et al., 2008; Arsouze et al., 2009; Rempfer et al., 2011; Gu et al., 2019; Wang et al., 2021) have demonstrated that reversible scavenging is an active and important component in global marine Nd cycling and is necessary for successfully simulating both [Nd]d and εNd distributions.

Updating the approach employed by Siddall et al. (2008), we prescribed individual biogenic particle export fields based on satellite-derived primary production. FAMOUS does contain an optional ocean biogeochemistry module (Hadley Centre Ocean Carbon Cycle; HadOCC), which includes simplified nutrient, phytoplankton, zooplankton, and detritus (NPZD) classes (Palmer and Totterdell, 2001) and could instead be used as the basis for predicting vertical particle fluxes in the ocean (which was the approach adopted by Arsouze et al., 2009; Gu et al., 2019; and Rempfer et al., 2011). We favoured satellite-derived estimates in order to improve the accuracy of particle-associated cycling of Nd and reduce biases inherent to the intermediate-complexity biogeochemistry model (Dentith et al., 2020; Palmer and Totterdell, 2001). This approach also optimises the computational efficiency of our scheme.

In our scheme, biogenic particle fields (POC, CaCO3, and opal) are prescribed using gridded, global satellite-derived particle export productivity from Dunne et al. (2012, 2007). The euphotic zone (zeu) is set to a globally uniform depth in FAMOUS of 81 m, which is the closest bottom grid box depth to match that defined in the OCMIP II protocol (75 m) (Najjar and Orr, 1998). Below zeu, appropriate depth-dependent dissolution profiles, derived from assumptions of particle degradability and sinking speed (Martin et al., 1987; Laws et al., 2000; Behrenfeld and Falkowski, 1997) and widely used to model flux attenuation in ocean models (e.g. all models used in the Ocean Carbon-Cycle Model Intercomparison Project; Doney et al., 2004; Sarmiento and Le Quéré, 1996), were applied to the biogenic export fluxes.

Downward fluxes of POC (FPOC) follow the power law profile of Martin et al. (1987):

(9) F POC z = F POC z eu × z z eu - for z > z eu ,

where z is depth (m) and represents a dimensionless dissolution constant for POC set to 0.9 (Najjar and Orr, 1998). Although a widely used parameterisation of dissolution, it should be noted this so-called “Martin curve” is known to underestimate the flux to the sediment in the off-equatorial tropics and subtropics and overestimates the flux in subpolar regions, indicating particles penetrate deeper than the Martin curve in the tropics and shallower in sub-polar regions (Dunne et al., 2007).

Downward fluxes of opal (Fopal) and CaCO3 (FCaCO3) follow exponential dissolution profiles with particle-specific length scales, i.e. lopal (10 000 m) and lcalcite (3500 m) (Maier-Reimer, 1993; Henderson et al., 1999; Najjar and Orr, 1998), respectively, calculated as follows:


The sinking of dust is prescribed according to the pre-industrial annual mean dust deposition simulated by Hopcroft et al. (2015) (see Sect. 2.3.1). We assume that dust does not dissolve significantly with depth, and thus dust export fluxes are constant throughout the water column. In line with previous schemes (e.g. Arsouze et al., 2009; Pöppelmeier et al., 2020b; Rempfer et al., 2011; Siddall et al., 2008), a uniform settling velocity (ω) of 1000 m yr−1 is applied to all particle fields. Although we acknowledge this is a simplification, the uniform settling rate applied captures the mean particle flux to the seafloor.

Annual averaged export of biogenic fields is shown in Fig. 7, and the total annual export production of POC (9.6 Gt(C) yr−1), CaCO3 (0.45 Gt(C) yr−1), and opal (90 Tmol(Si) yr−1) are comparable with previous estimates, although export of CaCO3 and opal are at the lower end, as highlighted in Table 3 by Dunne et al. (2007). Note that the annual export of CaCO3 reflects the new optimised surface calcite parameterisation as described in Dunne et al. (2012).

Figure 7Particle export fields from the ocean surface (g m−2 yr−1). Biogenic particle fields are prescribed using satellite-derived export productivity fields for (a) POC, (b) CaCO3, and (c) opal (Dunne et al., 2007, 2012). The dust input fields (d) are annual mean dust deposition simulated for the pre-industrial era by the HadGEM2-A GCM (Hopcroft et al., 2015). Note the different scale used for panel (a).

Reversible scavenging (Ndrs) considers total Nd for each isotope (j) as the dissolved Nd (Ndd) and particulate Nd (Ndp) associated with the different particle fields (χ, where χ refers to POC, CaCO3, opal, and dust).

(12) Nd t j = Nd d j + Nd p j = Nd d j + χ Nd p , χ j

Ndp sinks in the water-column with the particles due to gravitational force. Dissolution of biogenic particles with increasing depth below the euphotic zone releases Nd incorporated or adsorbed onto particles back to seawater (i.e. the Nd is dissolved or desorbed). Thus, particles act as internal sinks for marine Ndd in shallower depths and as sources at greater depths. This combined process for reversible scavenging (Srs) in the model can be described by

(13) S rs i , k = - ω Nd p i , k z i , k ,

where [Nd]p can be determined within the model from total Nd for each isotope as follows:

(14) Nd p j i , k = Nd t j i , k × 1 - 1 1 + χ R χ i , k × K χ j .

Rχ(i,k) describes the dimensionless ratio between particle concentration for each particle type (Cχ(i,k)) and the average density of seawater (ρ: 1024.5 kg m−3), input as fixed boundary conditions in our scheme. Cχ is calculated from the prescribed particle fluxes (Fχ: Fig. 7) by assuming a globally uniform settling velocity (ω=1000 m yr−1; Arsouze et al., 2009; Dutay et al., 2009; Gu et al., 2019; Kriest and Oschlies, 2008; Rempfer et al., 2011), i.e. Cχ=Fχ/ω. The equilibrium between dissolved concentration ([Nd]d) and concentration associated with particle type ([Nd]p) can be described by the equilibrium partition coefficient (Kχj):

(15) K χ j = Nd p , χ j Nd d j × 1 R χ ,

where [Nd]p/[Nd]d represents the scavenging efficiency in the model. It is independent of particle type and is used as a tuning parameter. Rχ, however, is dependent on particle type (where Rχ=Cχ/ρ), and thus Kχ is different for different particles. Our global mean particle concentrations (Cχ; Table 2) are similar to those reported in previous schemes.

Table 2Global mean particle concentrations for each particle type used to calculate equilibrium scavenging coefficients following Eq. (15). In summary, export fluxes of POC, CaCO3, and opal in this study are from Dunne et al. (2007, 2012), while dust fluxes are from Hopcroft et al. (2015). Global mean particle concentrations reported in previous Nd isotope schemes are shown for comparison.

Download Print Version | Download XLSX

Isotopic fractionation during absorption or incorporation and desorption or dissolution is neglected due to similar masses of 143Nd and 144Nd, avoiding undue complexity arising from any assumption about preferential scavenging (Siddall et al., 2008). Adsorption occurs everywhere that particles are present and we do not allow for preferential scavenging onto different particle types, consistent with previous models (e.g. Rempfer et al., 2011).

Therefore, by including advection and diffusion processes (transport), the total conservation equation for each Nd isotope in the model scheme can be written as follows:

(16) δ Nd t j δ t = S dust j + S river j + S sed j + S rs j + transport .

2.4 Evaluation methods and data sets

To validate the new Nd isotope scheme (Nd v1.0) and to assess the model's performance, we compare the simulated [Nd]d and εNd to modern seawater measurements, with a focus on describing the ability of the model to represent key spatial and vertical distributions across ocean basins.

As part of this assessment, a basic indication of model skill is returned by the mean absolute error (MAE):

(17) MAE = 1 N k = 1 N obs k - sim k ,

where obsk and simk are measured and simulated [Nd]d or εNd, respectively, N is the number of observations, and k is an index over all observational data. For each measurement – based on its longitude, latitude, and depth – the value predicted by the model is extracted and the mean deviation of simulated and observed [Nd]d and εNd is presented (in pmol kg−1 and epsilon units, respectively). Here we chose specifically not to apply a grid box volume weighting to the MAE, which would act to emphasise abyssal Pacific results in our assessment of model skill, where there are few observations and relatively low variability in Nd distributions. The advantage of using an unweighted MAE is that the assessment metric better scrutinises regions with larger (spatial) gradients in both [Nd]d or εNd, i.e. at the surface and high latitudes. We also report root-mean-square error (RMSE) for each simulation as an additional indicator of model skill and to allow for additional comparability with other studies:

(18) RMSE = 1 N k = 1 N ( obs k - sim k ) 2 N .

The observational data used in this assessment are from the seawater REE compilation used by Osborne et al. (2017, 2015), augmented with more recent measurements including data in the GEOTRACES Intermediate Data Product 2021 (GEOTRACES Intermediate Data Product Group, 2021) from GEOTRACES cruises (GA02, GA08, GP12, GN02, GN03 and GIPY05). Combined, our observational database represents a total of 6048 [Nd]d and 3278 εNd measurements, making it the largest compilation of seawater Nd data used to date to validate the performance of an Nd-isotope-enabled model. Notably, we omit measurements of [Nd]d>100 pmol kg−1 from the model data comparison because they represent very localised signals that we do not attempt to resolve. These include extreme surface concentrations present in restricted basins, such as fjords, the Baltic Sea, and the Gulf of Alaska (Chen et al., 2013; Haley et al., 2014), or input from hydrothermal activity (Chavagnac et al., 2018), which is not believed to govern global marine Nd distributions due to the immediate removal of hydrothermal Nd at the vent site (Goldstein and O'Nions, 1981). The location and spatial distribution of all observational records used in this study are shown in Fig. 8, and full details of the seawater compilation including a full list of all the references for the data sources are provided in Table S3 in the Supplement.

Figure 8Location of marine observational records used in this study: (a) filled orange triangles show the location of dissolved Nd concentration records, while (b) filled sky-blue circles show the location of dissolved εNd records.

Neither the horizontal nor vertical distribution of global seawater [Nd]d and εNd observational data is even. Most [Nd]dmeasurements are in the Pacific Ocean. This is in contrast to εNd measurements, which are most frequent in the Atlantic Ocean. Both [Nd]d and εNd observational data are biased towards shallower depths. It is therefore important to highlight that (as with previous studies) a skew in the distribution of seawater Nd measurements will act to bias our assessment of model performance.

In some instances, the location of the compared measured and modelled Nd may not match due to the coarseness of the model grid near land grid cells. In such cases, we employ a nearest neighbour algorithm to extract the modelled value from the closest ocean model grid cell. Furthermore, if multiple measurements occur within one model grid cell, the arithmetic mean of the values is used for our comparison to model results, and as such n is 3471 and 2136 for the calculation of MAE[Nd] and MAEεNd, respectively.

To ensure that our evaluation is not overly reliant on the cost function analysis alone, and hence evaluated the spatial biases in our assessment from the geographically uneven spread of measured [Nd]d and εNd, we also assess the capability of the model to reproduce appropriate global Nd inventories and to simulate large-scale horizontal and vertical gradients. We compare our results briefly with findings from previous modelling studies, but we highlight that the purpose of this study is to understand the behaviour of our model and not to undertake a comprehensive calibration of its performance. Optimisation (or “tuning”) of the Nd scheme will follow this work, and this needs to be considered when comparing the cost function performance to previous schemes.

2.5 Sensitivity experiment design

We designed a number of sensitivity experiments (Table 3) to systematically vary individual model parameters describing the reversible scavenging efficiency ([Nd]p/[Nd]d) and the main Nd source (fsed). These two parameters ([Nd]p/[Nd]d and fsed) were chosen primarily as they represent important and largely unconstrained non-conservative processes that are understood to govern simulated global distributions of both seawater [Nd]d and εNd (Rempfer et al., 2011; Pöppelmeier et al., 2020a; Gu et al., 2019; Arsouze et al., 2009; Siddall et al., 2008). By isolating individual effects, the primary aim was to attain a detailed understanding of the model's sensitivity to different forcings, to identify which parameters are important for [Nd]d and/or εNd patterns, and to identify assumptions within the explored parameters that require further constraining (through further field campaigns, laboratory analyses, or model experimentation).

Table 3Suite of FAMOUS simulations designed to assess the sensitivity of simulated [Nd]d and εNd distributions to two systematically varied parameters: reversible scavenging efficiency ([Nd]p/[Nd]d) and the global rate of direct Nd transfer from sediment to ocean water (fsed). Simulation name refers to the title given to each sensitivity simulation in this paper, and the simulation identifier refers to the unique five-letter Met Office identifier (which, for example, can be used to call down full experiment details from the NERC PUMA facility: last access: 15 February 2023). Global mean absolute error (MAE) and root-mean-square error (RMSE) for [Nd]d and εNd are shown as an indicator of model skill.

Download Print Version | Download XLSX

Due to time constraints, all sensitivity experiments presented here were run in parallel. First, [Nd]p/[Nd]d is systematically varied in six sensitivity simulations ([Nd]p/[Nd]d ranging 0.001–0.006), these values are based upon results from similar modelling schemes (Rempfer et al., 2011; Arsouze et al., 2009; Gu et al., 2019) and consider the few direct observations of [Nd]p/[Nd]d (Jeandel et al., 1995; Stichel et al., 2020; Zhang et al., 2008; Paffrath et al., 2021; Bertram and Elderfield, 1993; Tachikawa et al., 1999). Here, based upon simulations undertaken when validating the scheme and estimates from previous optimised Nd isotope schemes (Arsouze et al., 2009; Rempfer et al., 2011; Gu et al., 2019; Pöppelmeier et al., 2020a), fsed is fixed at 4.5 Gg yr−1 throughout.

Second, fsed is varied in four sensitivity simulations (fsed ranging 1.5–6.0 Gg yr−1), using values based upon previous recent estimates of a global sediment flux to seawater (3.3–5.5 Gg yr−1; Gu et al., 2019; Pöppelmeier et al., 2020b; Rempfer et al., 2011) and encompassing a larger parameter space in order to explore the sensitivity of Nd distributions. Notably, our fsed sensitivity studies alter the percentage contribution of the sediment Nd flux to the total Nd flux to seawater from 66 % (where fsed=1.5 Gg yr−1) to 89 % (where fsed=6.0 Gg yr−1). In all simulations [Nd]p/[Nd]d is fixed at 0.003, which is in the middle of the range (0.001–0.006) discussed above. Following the [Nd]p/[Nd]d sensitivity experiment (Sect. 3.1), it is apparent that 0.004 would have been a better choice for [Nd]p/[Nd]d, and we propose this as a starting point for further performance refinement in subsequent work.

Dissolved seawater Nd in all simulations is initialised from zero and integrated for at least 9000 years under constant pre-industrial boundary conditions to allow the deep-ocean circulation and marine Nd cycle to reach steady state, which we define as being when the Nd inventory becomes (near) constant with time (<0.02 % change per 100 years). All the presented results refer to or show the centennial mean from the end of the 9000-year simulations.

3 Results and discussion

3.1 Model sensitivity to reversible scavenging efficiency ([Nd]p/[Nd]d)

The first set of sensitivity experiments tests the response of simulated [Nd]d and εNd to a systematic variation in the reversible scavenging tuning parameter ([Nd]p/[Nd]d). Despite total Nd flux to seawater being fixed, varying the scavenging efficiency ([Nd]p/[Nd]d) leads to different Nd inventories (see Fig. 9 and Fig. S7 in the Supplement for the percentage rate of change per 100 years) and residence times (Table 3), consistent with previous studies (Siddall et al., 2008; Arsouze et al., 2009; Rempfer et al., 2011; Gu et al., 2019). A higher [Nd]p/[Nd]d increases both the Nd scavenging efficiency and removal via sedimentation by enabling a larger fraction of seawater Nd to adsorb onto particles, in turn leading to a lower Nd inventory and a lower residence time (where; residence time is equal to Nd inventory and total Nd flux).

Figure 9Global Nd inventory (g) simulated with different values for the reversible scavenging tuning parameter, [Nd]p/[Nd]d, as indicated. Dashed line represents the estimated global marine Nd inventory of 4.2 × 1012 g by Tachikawa et al. (2003) used as a first-order target for our simulations.

By year 6000, experiments EXPT_RS3–EXPT_RS5 reach steady state. We therefore deem these simulations to have reached an acceptable equilibrium state. We explain why EXPT_RS1, EXPT_RS2, and EXPT_RS6 do not reach equilibrium below, but this is also because of their unrealistic condition after 9000 years. For example, EXPT_RS1 and EXPT_RS2 reach an Nd inventory far above the target inventory of 4.2 Tg (Tachikawa et al., 2003), while the high sink in EXPT_RS6 causes imbalances in surface [Nd]d (which tends towards zero). As a result, these simulations are largely omitted from further discussion and analysis.

Neodymium in EXPT_RS1, which has the lowest [Nd]p/[Nd]d, has a residence time of 3036 years. This is much larger than the global ocean overturning time of 1500 years, resulting in an Nd inventory of 16 Tg after 9000 years. The [Nd]p/[Nd]d, and consequently the sink in EXPT_RS1, is too small to balance the input of Nd from the sources, causing such high Nd to accumulate in the simulation, which as a result does not reach steady state (rate of change 0.16 % per 100 years at the end of the simulation; Fig. 9). Thus, EXPT_RS1 returns the largest MAE[Nd] and MAEεNd from these simulations, indicating the worst model–data fit, especially for representing [Nd]d measurements. EXPT_RS6, the simulation with the largest [Nd]p/[Nd]d, returns the lowest MAEεNd. However, the sink term is too strong, resulting in a Nd inventory of 2.95 Tg, and the simulation fails to reach steady state (rate of change 0.07 % per 100 years by the end of the run). On balance, we surmise that EXPT_RS4 has the most reasonable combination of prognostic skill in terms of simulated Nd inventory (4.51 Tg) and residence time (856 years). This simulation, with [Nd]p/[Nd]d of 0.004, shows the best balance between Nd accumulation and removal, and hence returns the lowest MAE[Nd] alongside a moderately well-performing MAEεNd that falls in the middle of the range of results.

In contrast to our model, the schemes by Rempfer et al. (2011) and Gu et al. (2019) required a lower [Nd]p/[Nd]d of 0.001 and 0.0009, respectively, for their optimised experiment with Nd inventories of ≈4.2 Tg. Despite having similar scavenging schemes, a direct comparison of the parameter values used in the different Nd isotope modelling studies is difficult to make. This is because the divergence in sensitivity to reversible scavenging efficiency can be attributed to a combination of the differing magnitude and spatial distributions of model biogeochemical particle fields and Nd inputs, which are also partly controlled by the different architecture and horizontal resolution of the physical models. We thus propose that a future modelling protocol for intercomparing different global Nd isotope schemes would be well suited to exploring these differing sensitivities comprehensively.

The different values of MAE[Nd] and MAEεNd across our sensitivity experiments (Table 3) illustrate the distinctive and uncoupled behaviour of [Nd]d and εNd within the ocean, as broadly described by the Nd paradox, indicating that different processes govern the global distributions of each.

Generally, simulated [Nd]d distributions in EXPT_RS4 match observational data well (Fig. 10). The lowest concentrations occur in the surface layers, and deep-water [Nd]d increases along the global circulation pathway and with increasing age of water masses (Bertram and Elderfield, 1993; van de Flierdt et al., 2016; Tachikawa et al., 2017). Thus, we may infer that the model scheme in FAMOUS does have the broad capability of representing the physical and biogeochemical processes governing global marine [Nd]d distributions. However, the scheme does tend to overestimate the global vertical [Nd]d gradient with depth (see Fig. 10 and Fig. S9 in the Supplement for major ocean basin averaged depth profiles), a feature reported in previous similar model schemes (e.g. Arsouze et al., 2009; Gu et al., 2019), indicating that the representation of processes governing vertical [Nd]d does not yet fully capture all processes occurring in the ocean. This overestimated vertical [Nd]d gradient with depth could be a result of the simplified uniform particle settling rate applied, which may be too high in the upper layers (due to lower settling velocities in the surface mixed layer resulting from enhanced turbulence and complexation with organic ligands; Chamecki et al., 2019; Noh et al., 2006) and too low at depth (due to higher settling velocities from particle aggregation).

Figure 10Maps of [Nd]d (left column) and εNd (right column) in simulation EXPT_RS4 split into four different depth bins: (a–b) shallow (0–200 m), (c–d) intermediate (200–1000 m), (e–f) deep (1000–3000 m), and (g–h) deep abyssal ocean (>3000 m). Water column measurements from within each depth bin (Osborne et al., 2017, 2015; GEOTRACES Intermediate Data Product Group, 2021) are superimposed as filled circles using the same colour scale.

Too high simulated [Nd]d at depth may also be caused by the benthic sediment source being too strong, coupled with the simplified assumption of a globally uniform sediment flux acting across the entire seafloor. This suggests that perhaps too much emphasis has been put on a widespread abyssal benthic flux to resolve the Nd paradox, and it may be that shallower and marginal regions pose a larger in magnitude and more distinct Nd source to seawater as opposed to an open-ocean benthic source. Recent pore water measurements by Deng et al. (2022) on the East China continental shelf revealed larger benthic fluxes across the dynamic hydraulic environment of the continental shelf compared to deep pore water, which acted as a sink, likely due to authigenic precipitation. As such, spatial variations in the benthic flux are likely affected by multiple complex depth-related processes, although deep-ocean regions may still dominate the total area-integrated benthic flux. A future evolution of the scheme could explore depth-dependent benthic flux rates.

Moreover, the overestimation of [Nd]d at depth could be a result of biases within the simulated biogenic particle fields. Comprehensive global observations of sinking particles fluxes – the central driver of ocean biogeochemical cycling – remain fundamentally difficult to obtain (Dunne et al., 2007), and our particle fluxes may be inaccurate as a consequence (see Sect. 2.3.4 for assumptions and respective limitations). It also seems likely that the reversible scavenging parameterisations, which are simple by design due to incomplete understanding (also Sect. 2.3.4), restrict the model's ability to precisely capture all aspects of the measured Nd distributions. For example, we know that the assumption of a globally uniform [Nd]p/[Nd]d is oversimplified, as particle measurements suggest [Nd]p/[Nd]d is not uniform in seawater (Stichel et al., 2020). High remineralisation rates, inducing increased exchange, have been observed in the Labrador Sea (Lagarde et al., 2020), whilst vast expanses of low biological productivity in the Pacific gyres mean [Nd]p levels are generally much lower here compared to the North Atlantic. Thus, the global [Nd]p/[Nd]d in our scheme may be skewed towards lower values. Nonetheless, it is necessary to adopt a pragmatic approach for the model development, and hence (similar to other model schemes) the globally uniform [Nd]p/[Nd]d remains in our implementation. Future work could explore regionally varied [Nd]p/[Nd]d, which would benefit from additional measurement constraints to sufficiently justify a more complex distribution. Additionally, other particles such as Fe and Mn oxides and hydroxides that are not considered here may play an important role in scavenging of Nd (e.g. Bayon et al., 2004; Lagarde et al., 2020; Du et al., 2022). Further observational evidence of the processes involved and their importance for Nd cycling by combined particulate and dissolved measurements and laboratory experiments (e.g. Stichel et al., 2020; Pearce et al., 2013; Rousseau et al., 2015; Wang et al., 2021; Paffrath et al., 2021; Lagarde et al., 2020), plus experimentation within a modelling framework, may help to improve this limitation.

However, the largest model–data disparities for [Nd]d occur in the shallow ocean (above 200 m) at specific locations close to continental margins where Nd is input to the ocean through major point sources that are not well resolved by the model. This includes continental margins in the Labrador Sea (where simulated [Nd]d is 3 pmol kg−1 compared to measured [Nd]d of 70 pmol kg−1) and the Sea of Japan (where simulated [Nd]d is 2 pmol kg−1 compared to measured [Nd]d of 50 pmol kg−1). Such low [Nd]d in the surface layers may be exacerbated by operational constraints in the scheme, such as the extensive and immediate dilution of point-sourced Nd across the whole of its containing grid cell combined with the instantaneous nature of simulated reversible scavenging, which may be much faster in the model than would occur normally. Furthermore, the scheme may be underestimating the magnitude of surface sources of Nd, for example the Nd sourced from dust may be too low in certain regions, possibly resulting from the globally fixed (2 %) solubility of Nd from dust sources applied. In the preliminary optimisation of the Nd isotope scheme in the GNOM v1.0 model, Pasquier et al. (2022) found that dust solubility needed to be increased to unrealistically high levels to best fit measurements; for example, the optimised parameter for North American dust solubility was 83 %. A future exploration of the regional solubility and amount of Nd coming from dust sources would be needed to explore this in detail. On the other hand, recent evidence suggests that there may be a major particulate riverine flux of Nd to the ocean (Rousseau et al., 2015; Rahlf et al., 2021), which implies that the magnitude of the sediment Nd sourced to seawater near river mouths could be too weak. Specific to our scheme, the vertical convection of water from the physical ocean model in FAMOUS may also be too intense in the North Atlantic, leading to surface underestimations of [Nd]d and overestimations at depth; a re-tuning of the physical model would help to ascertain the extent to which this limitation contributes towards the simulated Nd bias.

Nonetheless, consistent with compilations of water column measurements (Tachikawa et al., 2017; van de Flierdt et al., 2016), overall global distributions of simulated εNd are broadly most non-radiogenic in the North Atlantic and more radiogenic in the North Pacific (Fig. 10), with intermediate values in the Southern and Indian oceans. The most non-radiogenic εNd occurs in the surface layers of the Hudson Bay and Labrador Sea regions, and they closely match measured data (εNd=-18). However, the most radiogenic εNd, simulated in the surface layers of marginal regions in the North Pacific and equatorial western Pacific (εNd=-3), is significantly lower than measured (εNd=+3). In the central and North Pacific, particularly above 1000 m, simulated εNd is −7, but measurements are closer to −1. At the basin scale, the magnitude of the εNd gradient from Pacific to Atlantic is underestimated by the model and presents a familiar bias when compared to previous Nd isotope schemes (Arsouze et al., 2009; Rempfer et al., 2011; Pöppelmeier et al., 2020a; Jones et al., 2008; Gu et al., 2019). This is mainly due to the simulated Pacific being too non-radiogenic (basinal mean εNd of −7.5) compared to measured water samples (εNd=-4), while simulated and measured basinal mean Atlantic εNd values are in much better agreement (εNd=-11 and −12.5, respectively; Supplement Fig. S9). Possible explanations for why our simulated Pacific εNd is so much lower than what has been recorded in modern seawater include that the model boundary conditions, specifically the marine sediment source of Nd (taken directly from Robinson et al., 2021), may not be sufficiently radiogenic or that the sediment flux from the seafloor is not uniform; points we will return to later (Sect. 3.2).

Simulated [Nd]d with depth at distinct locations in all the reversible scavenging sensitivity experiments (Fig. 11a–f and h–m) generally (though not always) exhibit similar depth profiles to the observational data. The best model–data fit is seen within depth profiles in the Pacific and Southern oceans and notably more so under higher [Nd]p/[Nd]d. The largest model–data offsets across all sensitivity experiments in terms of both the magnitude and the depth gradients in [Nd]d occur in the North Atlantic Ocean. In particular, the large near-surface concentrations are not resolved, with the largest disparities occurring under higher scavenging parameterisations.

Figure 11The central panel (g) displays [Nd]d at the seafloor (i.e. lowermost box of the ocean model) in simulation EXPT_RS4 (100-year mean from the end of the run), with superimposed water column measurements (Osborne et al., 2017, 2015; GEOTRACES Intermediate Data Product Group, 2021) from ≥3000 m shown by filled coloured circles on the same colour scale. The surrounding panels (a)(f) and (h)(m) display depth profiles of simulated (coloured lines, one per sensitivity simulation with varied [Nd]p/[Nd]d) and measured (filled circles) [Nd]d, respectively. Larger shifts in the [Nd]d between simulations highlight the regions most sensitive to the efficiency of reversible scavenging.

For εNd, the response to varying scavenging efficiency has varied effects across depths and between ocean regions, indicating a more complex relationship between how reversible scavenging can delineate global εNd distributions in comparison to [Nd]d. Increased scavenging efficiency traps regional εNd provenance signals, resulting in larger inter-basin εNd gradients, the strength of which are favoured by the combination of strong scavenging efficiency and a short seawater residence time. Although the simulated εNd profiles with depth do not closely follow measurements, it is valuable to discuss the sensitivity of simulated εNd to variations in scavenging efficiency. The North Atlantic is the most sensitive basin to changes in reversible scavenging (registered by the greater εNd profile shifts between different experiments, Fig. 12a–d), particularly at depths below 2000 m, indicating that reversible scavenging, alongside convection at sites of deep-water formation, are important processes for governing the simulated εNd signal of deep-water masses. Thus, a stronger reversible scavenging efficiency is needed to set the non-radiogenic deep ocean signature in the North Atlantic and to trap these regional εNd provenance signals locally.

Figure 12The central panel (g) displays εNd at the seafloor (i.e. lowermost box of the ocean model) in simulation EXPT_RS4 (100-year mean from the end of the run), with superimposed water column measurements (Osborne et al., 2017, 2015; GEOTRACES Intermediate Data Product Group, 2021) from ≥3000 m shown by filled coloured circles on the same colour scale. The surrounding panels (a)(f) and (h)(m) display depth profiles of simulated (coloured lines, one per sensitivity simulation with varied [Nd]p/[Nd]d) and measured (filled circles) εNd, respectively. Larger shifts in the εNd between simulations highlight the regions most sensitive to the efficiency of reversible scavenging.

Typical depth profiles from measurements in the subtropical North Atlantic show contrasts in εNd that co-vary with the presence of major water masses (coloured circles in Fig. 12a–b). Across all sensitivity experiments, there is relative consistency with depth for the profiles in the subtropical North Atlantic (e.g. Fig. 12b), ranging from −10 to −12 under increasing scavenging efficiency, which acts to drive more localised non-radiogenic signals. This simulated uniformity may be due to the indiscriminate seafloor sedimentary flux here being too high and thereby dominating the deep seawater εNd signal. An additional explanation may be due to the lack of abyssal (i.e. below 3 km) AABW, which does not extend past  20 N (Fig. 3). This insufficient AABW production and penetration into the Atlantic are known limitations of FAMOUS (Dentith et al., 2019; Smith, 2012), although these biases are reduced in our control compared to the previous studies (Sect. 2.2). Thus, the seawater below the surface mixed layer comprises only North Atlantic water, and in so doing the model will not resolve the AABW signal inferred from measurements at depth in higher latitudes.

Higher up the water column, the more non-radiogenic NADW seen in the seawater measurements is conspicuously absent (note the εNd minima in subtropical North Atlantic measurements  1000–2000 m deep, Fig. 12a–b), even though we know NADW does reside there in the model (see Fig. 3 for verification). We may relate this to the high-latitude North Atlantic εNd being too radiogenic to tag NADW with its characteristically non-radiogenic εNd signal, particularly around the mouth of the Labrador Sea (Fig. 10a–b). This could be exacerbated by a dampening of the non-radiogenic NADW εNd from a relatively radiogenic seafloor benthic flux along the water flow path (Fig. 6). Consequently, despite a close correlation to seawater εNd in the subtropics and lower NADW (where measurements of εNd are −12.4 and simulated is −12), the upper NADW end-member εNd is not sufficiently non-radiogenic, even under the highest reversible scavenging efficiency ([Nd]p/[Nd]d=0.006), where simulated εNd is −12, which is in contrast to seawater records of −13.2 εNd (Lambelet et al., 2016).

In contrast, εNd in the Pacific Ocean is least sensitive to changes in [Nd]p/[Nd]d, particularly below 500 m (Fig. 12e–f, h–i). This is expected due to the absence of major ocean convection and ventilation, meaning the Pacific contains an older and more homogenised pool of water in comparison to the Atlantic, and thus water masses are far less distinct and the localised εNd signal does not get dispersed via convection as rapidly. The deep seafloor-wide sediment source governs simulated εNd distributions in the intermediate to deep Pacific due to its larger flux relative to that transported vertically via particle scavenging and dissolution or horizontally by the sluggish convection of the Pacific. One further aspect to note is that in some Pacific regions [Nd]d in the surface layers tends towards zero under a high sink ([Nd]p/[Nd]d=0.006). This causes numerical instabilities in modelled Nd ratios, and thus εNd is meaningless where [Nd]d<0.2 pmol kg−1 and is therefore masked out of the results shown by Fig. 12 to avoid misinterpretation.

The simulated depth profile in all experiments in the Indian Ocean matches the observed intermediate εNd signal of −8 between 500 and 2000 m (Fig. 12m). However, the more radiogenic εNd signal of −6 in the surface layers is not captured and neither is the shift in εNd below 2000 m to more non-radiogenic values reaching a minimum of −10.5 at 3000 m, with all experiments simulating a relatively uniform εNd with depth. Under higher reversible scavenging parameters for this depth profile (EXPT_RS4–EXPT_RS6), the surface concentrations are too low (simulated [Nd]d=3 pmol kg−1; observed [Nd]d=8 pmol kg−1), indicating that the model is not fully resolving either the concentration or the εNd from a surface flux here and that scavenging may be too intense at the surface. In the deeper ocean (≈3000 m), these simulations show the model represents the [Nd]d profiles better but is missing an non-radiogenic exchange of εNd at depth. This could point to an non-radiogenic sediment source or exchange that is misrepresented by the bulk sediment boundary conditions and sediment source assumptions (e.g. our application of a global seafloor sediment source of Nd irrespective of sedimentary characteristics).

In the Southern Ocean, simulations with [Nd]p/[Nd]d≥0.003 broadly match the general measured εNd at depths above 1000 m (decreasing with depth from −7.8 at the surface to −8.2 at 1000 m). Below this (and down to 3000 m), observed diversions in εNd are not captured in any of the sensitivity experiments. This discrepancy can be attributed to the simulation of quite homogenous AABW throughout the water column in the region, which represents a physical bias of FAMOUS (Dentith et al., 2019; Smith, 2012). Specifically, in the Pacific sector of the Southern Ocean (Fig. 12j), the model cannot resolve the measured non-radiogenic spike at 1500 m, which captures the distinct presence of lower Circumpolar Deep Water (CDW, εNd=-8.4± 1.6: Lambelet et al., 2018) formed from mixing of Atlantic-, Pacific-, and Indian-sourced waters.

Thus, we conclude that scavenging and removal via sedimentation is necessary to balance the simulated input sources and enables the scheme to reach equilibrium around reasonable Nd inventories. We find that reversible scavenging is an important process that enhances the εNd gradient between oceans by maintaining localised basinal εNd signals throughout the water column. The strength of this process is particularly important for maintaining the simulated non-radiogenic εNd in the well-ventilated North Atlantic Ocean (alongside deep-water formation) but less important for the more stagnant modern Pacific Ocean.

Nonetheless, parameterising reversible scavenging efficiency alone cannot account for the correct trends and magnitude within εNd gradients observed between basins and in depth profiles. Currently, the scheme assumes that all particles reaching the seafloor via reversible scavenging are buried in the sediment and as such are decoupled from the seafloor sediment source, whereas there is no Nd release and the particles are removed from the model. A future evolution of the scheme could explore the dissolution of authigenic sedimentary phases on the seafloor during diagenesis to investigate how this may influence the εNd distributions of the benthic flux. Moreover, our results demonstrate the importance of further constraining other aspects of marine Nd cycling under a holistic framework. Furthermore, we note that the scheme described here may not be fully resolving the end-member εNd of different water masses due to an imperfect representation of the sources of Nd to seawater (i.e. the model boundary conditions and strength of the source fluxes). Thus, in some instances, the scheme carries inappropriate or dampened εNd signals, coupled alongside particular structural model biases in the physical circulation (e.g. limited AABW intrusion in the North Atlantic). Although in the case of the latter point, we note that this is not a limitation of the presented scheme and that the implementation can be useful for identifying such physical biases.

3.2 Model sensitivity to Nd flux from the sediment (fsed)

The second set of sensitivity simulations tests the response of [Nd]d and εNd to systematically varying the total Nd flux from the sediment (fsed). In this experiment, Nd accumulates rapidly from the start of the simulations (Fig. 13 and Supplement: Fig. S8 for the percentage rate of change per 100 years) and tapers off thereafter to varying degrees depending on the rate of accumulation. By year 6000 all fsed sensitivity experiments have reached steady state (<0.02 % change per 100 years).

Figure 13Global Nd inventory (Tg) simulated with different values for the total sediment flux tuning parameter (fsed) as indicated. The dashed line represents the estimated global marine Nd inventory of 4.2 Tg from Tachikawa et al. (2003) used as an approximate target for the simulations. Note the logarithmic scale on the y axis.

Modifying the sediment flux while keeping the sink term [Nd]p/[Nd]d constant returns a varied equilibrium Nd inventory across the suite of simulations (Fig. 13). However, without changing the scavenging efficiency (i.e. Nd sink), there is a more moderate range of 40 years difference in the fsed simulations' Nd residence time (Table 3). In general, the relatively low scavenging efficiency of all of the fsed sensitivity simulations ([Nd]p/[Nd]d=0.003; see Sect. 3.1 for context) yields long residence times greater than 1000 years.

EXPT_SED1, which has the lowest fsed, corresponding to a sediment flux (per area) of 1.17 pmol cm2 yr−1, and hence the lowest total Nd flux to the ocean, returns the smallest total Nd inventory of 2.5 Tg (Table 3). Conversely, EXPT_SED4, which has the largest fsed, corresponding to a sediment flux of 4.69 pmol cm2 yr−1, results in the greatest Nd inventory (7.8 Tg), producing both the worst MAE[Nd] and worst MAEεNd.

We note that compared to varying the reversible scavenging efficiency, varying fsed drives relatively small changes in Nd distributions, as demonstrated by the minor differences in MAEεNd between sensitivity experiments. This makes sense since varying fsed only acts to change the fractional contribution from each specific Nd source, e.g. an enhanced fsed reduces the fraction of total Nd flux coming from dust and rivers (which are inputs constrained to the surface and point sources close to the continents), concentrating the flux across the global seafloor (and vice versa). However, in all simulations and consistent with previous studies, the sediment Nd source to seawater remains the major source, and fsed would need to be reduced much more to greatly influence MAEεNd.

Overall, from the simulations in this study, EXPT_SED2 demonstrates the best skill for [Nd]d compared to measurements, returning the lowest MAE[Nd]of 7.96 pmol kg−1 (Table 3) and achieving the most balanced Nd source and sink terms. Although the simulation does not represent the lowest MAEεNd, the range across the fsed experiment is small (2.71 to 2.93; Table 3), and it does simulate the highest percentage of simulated εNd within 3 εNd units of measurements (62 %). The sediment flux in EXPT_SED2 of 2.35 pmol cm2 yr−1 is considered similar to the few existing measurements of a benthic flux across different sedimentary environments, which are mostly on the order of 10 pmol cm2 yr−1 (Abbott, 2019; Abbott et al., 2015a; Haley et al., 2004; Du et al., 2018). In comparison with model estimates of the benthic flux, the sediment flux in EXPT_SED2 is most comparable to the global benthic flux estimated by Pöppelmeier et al. (2020a) of 5.4 pmol cm2 yr−1 but an order of magnitude lower than the value of 20–30 pmol cm2 yr−1 estimated by Du et al. (2020). Alternatively, the depth scaling of the benthic flux in the modelling study by Pasquier et al. (2022) resulted in surface fluxes of 83.7 pmol cm2 yr−1 and deep-ocean fluxes of 1.11 pmol cm2 yr−1.

The target global Nd inventory of 4.2 Tg, derived from the canonical work by Tachikawa et al. (2003), is based on a sparser array of seawater measurements, particularly in the Southern Ocean, compared to what is available today. Using our newer compilation of [Nd]d to identify the best-performing simulation for [Nd]d (EXPT_SED2), we can report an updated estimate of the global Nd inventory from the model of 3.89 Tg (Table 4), noting that this is likely an underestimate of the total since the model does not represent marginal seas. Furthermore, we can examine the basinal distribution of [Nd]d, (also Table 4), whereby the largest amount is found in the stagnant North Pacific, where the greatest area-integrated sediment flux occurs. Smaller Nd inventories are found in the more convective North Atlantic and Southern Ocean regions, which have less expansive seafloors and larger biogenic particle fluxes, resulting in more Nd removal via reversible scavenging and sedimentation. We propose that future work should take advantage of the greater abundance of seawater Nd measurements (e.g. as compiled in this study) to produce a detailed update to the estimate of oceanic inventories of Nd. These more accurate Nd budgets across diverse ocean regions and depths would help in improving constraints on the marine Nd cycle and would provide excellent model evaluation tools.

Table 4Nd inventory (Tg) for EXPT_SED2 within different ocean regions.

Download Print Version | Download XLSX

Altogether, simulated [Nd]d in EXPT_SED2 matches the general observational data trend of Nd concentration. However, in the deep layers of the North Pacific (below 3000 m) simulated [Nd]d is underestimated (36 pmol kg−1 compared to 50 pmol kg−1 from seawater measurements). This underestimation of [Nd]d at depth can be explained by a combination of having a seafloor sediment source at the lower end of our range (fsed is 3.0 Tg yr−1) and slow release from reversible scavenging ([Nd]p/[Nd]d of 0.003; our experiments suggest 0.004 may be a more suitable efficiency to use, Sect. 3.1), highlighting the importance of both non-conservative processes in governing deep [Nd]d distributions.

Figure 14Maps of [Nd]d (left) and εNd (right) in simulation EXPT_SED2 split into four different depth bins: (a–b) shallow (0–200 m), (c–d) intermediate (200–1000 m), (e–f) deep (1000–3000 m), and (g–h) deep abyssal ocean (>3000 m). Water column measurements from within each depth bin (Osborne et al., 2017, 2015; GEOTRACES Intermediate Data Product Group, 2021) are superimposed as filled circles using the same colour scale.

Moreover, the longer lifetime of simulated Nd in the EXPT_SED2 ocean (1032 years) compared to that estimated by previous work (360–800 years; Gu et al., 2019; Rempfer et al., 2011; Siddall et al., 2008; Tachikawa et al., 2003) means that Nd becomes well mixed in the deep ocean (below 1000 m), homogenising the εNd signal. This causes the observed inter-basin gradients (a critical feature in the use of Nd as an ocean circulation tracer) to become severely damped, particularly away from direct input of fresh reactive phases with distinctive εNd (Robinson et al., 2021; Abbott et al., 2019). Consistent with earlier studies (e.g. Arsouze et al., 2009; Jones et al., 2008), the largest model–data disparities occur in the North Pacific and equatorial Pacific, where simulated εNd is far too non-radiogenic compared to the observational data, pointing to a number of processes that may be better optimised in our Nd scheme. For example, a larger (and also more radiogenic) sediment source (i.e. greater fsed) may be needed, particularly around shallow to intermediate marginal settings, which is a suggestion also supported by the too low [Nd]d.

Total Nd concentrations and isotopic distributions show different responses to varying fsed. We find that [Nd]d is sensitive to fsed across all ocean basins (i.e. a wide divergence in the depth profiles shown in Fig. 15), mostly at depths below 500 m where there is no direct influence from river and dust inputs and the relatively large area of the deep abyssal seafloor as an Nd interface becomes important (particularly with low reversible scavenging).

Figure 15The central panel (g) displays [Nd]d at the seafloor (i.e. lowermost box of the ocean model) in simulation EXPT_SED2 (100-year mean from the end of the run), with superimposed water column measurements (Osborne et al., 2017, 2015; GEOTRACES Intermediate Data Product Group, 2021) from ≥3000 m shown by filled coloured circles on the same colour scale. The surrounding panels (a)(f) and (h)(m) display depth profiles of simulated (coloured lines, one per sensitivity simulation with varied fsed) and measured (filled circles) [Nd]d, respectively. Larger shifts in the [Nd]d between simulations highlight the regions most sensitive to the magnitude of the seafloor sediment source.

Typically, the best model–data fit for [Nd]d depth profiles is achieved under lower fsed (1.5 to 3.0 Gg yr−1) and in the Pacific, Indian, and Southern oceans in particular (Fig. 15). Correspondingly, under high fsed and in the deeper interior of the ocean in particular (below 3000 m), simulated [Nd]d diverges to higher concentrations than measured due to too much Nd being sourced to the deep ocean from the sediment. Additionally, in these high fsed scenarios, the strength of the deep sediment source obscures horizontal seafloor [Nd]d gradients across basins (see the near-uniform seafloor [Nd]d in Fig. 15g), masking the influence of reversible scavenging, which is controlled by the location and dissolution of particle fields and is important for governing [Nd]d patterns (Sect. 3.1).

Similar to the reversible scavenging sensitivity experiment (Sect. 3.1), the largest simulated offsets between simulated and measured [Nd]d under all fsed experiments occur in the North Atlantic and subtropical Atlantic at depths above 1000 m, where, even under the largest sediment fluxes, simulated [Nd]d is too low. The depth profile south of Greenland (Fig. 15c) shows the greatest sensitivity to varying fsed in the upper 1000 m, with EXPT_SED4's simulated surface concentrations of 10 pmol kg−1 (fsed=6.0 Gg yr−1) being closest to the measured concentrations of 22 pmol kg−1. By implication, the accurate representation of sediment fluxes is required to reproduce upper-ocean [Nd]d in this region, and an enhanced sediment flux alone cannot account fully for the observed high surface concentrations. Either the combination of the surface and near-surface fluxes resolved here is too diluted in the model (possibly due to grid box resolution) or the model is missing a significant surface or near-surface Nd source. Previous schemes have likewise simulated too low surface [Nd]d in the North Atlantic, likely also due to difficulties in representing highly localised and variable surface features in global models (Gu et al., 2019; Rempfer et al., 2011); more direct observations of surface fluxes would benefit future constraints in this region.

Whereas [Nd]d is sensitive to fsed globally, the height of εNd sensitivity is more regional (Fig. 16). The Atlantic, Indian, and Southern oceans have εNd values that are significantly more sensitive to changes in the bulk seafloor sediment Nd flux than the Pacific Ocean. Overall, differences in fsed tend to drive whole depth profile shifts of low magnitude in εNd. This contrasts the findings of Rempfer et al. (2011), who varied a margin-constrained fsed and reported that deep-water εNd (below 1000 m) values were affected less than in our results. If this conflicting difference in deep-ocean sensitivity to fsed is because of the applied spatial distributions of a sediment Nd source to seawater (<3000 m for Rempfer et al., 2011, and all depths for this study), then further constraints on sediment Nd fluxes across space and time are crucial when interpreting ocean circulation from εNd.

Figure 16The central panel (g) displays εNd at the seafloor (i.e. lowermost box of the ocean model) in simulation EXPT_SED2 (100-year mean from the end of the run), with superimposed water column measurements (Osborne et al., 2017, 2015; GEOTRACES Intermediate Data Product Group, 2021) from ≥3000 m shown by filled coloured circles on the same colour scale. The surrounding panels (a)(f) and (h)(m) display depth profiles of simulated (coloured lines, one per sensitivity simulation with varied fsed) and measured (filled circles) εNd, respectively. Larger shifts in the εNd between simulations highlight the regions most sensitive to the magnitude of the seafloor sediment source.

Once again, the most notable εNd model–data mismatch occurs within the depth profiles of the North Atlantic. Due to major rivers delivering a high [Nd] load to the Atlantic (Fig. 5c) in their vicinity, simulations with low fsed are conditioned towards riverine εNd, providing the typical non-radiogenic εNd signature of NADW (Fig. 16a–d). Conversely, enhancing fsed increases the fraction of Nd supplied to the ocean from the bulk seafloor sediment, which has more uniform, intermediate εNd values in the central North Atlantic (−12.5) in contrast to the more non-radiogenic εNd signals of the continental margins and riverine source in the Labrador (−28) and West Atlantic basins (−15.6), with localised near-surface sediment extreme minimums of −34 in the northern Labrador Sea (Robinson et al., 2021). Greater fsed thus acts to overprint and hence mix away the more distinct surface εNd signal of NADW gained at their sites of deep-water formation in favour of a more general intermediate εNd signal as it becomes exposed to Nd fluxes along its southward seafloor flow path. Model–measurement disparity at the mouth of the Labrador Sea and south of Greenland suggests that a specific fraction of the bulk sediment with more non-radiogenic εNd than is captured by Robinson et al. (2021) may be driving particle–seawater interactions (Fig. 17). Mechanical glacial erosion, which exposes large amounts of highly non-radiogenic and labile fine-grained crystalline detritus in the region surrounding the Labrador Sea likely drives enhanced localised extreme non-radiogenic benthic fluxes (Von Blanckenburg and Nägler, 2001; Anderson, 2005). Moreover, marine particle εNd measured in the North Atlantic has been shown to reflect the strong heterogeneity of the surrounding source rock, for example the provenance of extremely radiogenic particles (+6 to +9) at the Irminger Basin has been associated with mafic rocks found in eastern Greenland and Iceland (Stichel et al., 2020). It therefore follows that better spatial constraints on sediment fluxes are required for providing accurate εNd tagging of Atlantic seawater.

Figure 17Volume-weighted distributions of εNd in simulation EXPT_SED2 split into two different depth bins, (a) shallow (0–200 m) and (b) intermediate (200–1000 m), within the North Atlantic and Labrador Sea basins. Water column measurements from within each depth bin (Osborne et al., 2017, 2015; GEOTRACES Intermediate Data Product Group, 2021) are superimposed as filled circles on the same colour scale. Note that Iceland is not resolved in FAMOUS (Jones, 2003).

In contrast to the Atlantic, the Pacific Ocean εNd is the least sensitive to varying fsed (Fig. 16e–f, h–i). Characteristically, the Pacific encompasses vast open-ocean oligotrophic expanses with low biogenic particle export, in tandem with smaller relative dust and river sources, which means the seafloor sediment source already dominates the simulated Nd fluxes and distributions, even under lower sediment fluxes. In the North Pacific, the εNd signal from the bulk seafloor sediment cannot explain the radiogenic εNd observed here (Fig. 16), and thus the model is likely not fully capturing the correct source (or the prescribed boundary condition is not correctly representing the distribution and fraction of the sediment phase contributing to the sediment flux) of Nd to seawater. Authigenic precipitation of clay minerals during reverse weathering is prevalent in the pelagic sediment throughout the North Pacific. These top centimetres of sediment are typically assumed to be unreactive and act as an active sink for REE (Kato et al., 2011). Moreover, precipitation at hydrothermal vents also acts as an effective localised REE sink (German et al., 1990). Therefore, the red clays of the deep Pacific may act as an effective Nd sink and not a source. Alternatively, a recent reactive transport model by Du et al. (2022) suggests reverse weathering maintains under-saturation of primary silicates in deep Pacific pore water, driving enhanced silicate weathering and favouring preferential dissolution of highly radiogenic volcanic sediments (e.g. basaltic glass and clinopyroxene). Here, we highlight the need for more precise constraints on complex benthic processes in the Pacific, including exploring the environmental conditions that would drive spatially elevated and “reactivity-weighted” sediment fluxes.

To conclude, we surmise that a sediment–seawater flux represents a key major source of [Nd]d that is particularly fundamental to the intermediate and deep-ocean Nd budgets and as such plays an important role in governing marine Nd cycling. Indeed, a strong seafloor-wide uniform sediment source pushes the ocean towards a more globally uniform and too high [Nd]d than measurements from the deep ocean suggest. This could be interpreted as evidence against a globally widespread benthic-flux-driven model of the marine Nd cycle with a spatially constant flux across diverse sedimentary environments in favour of the more distinct [Nd]d distributions that may be achieved under a model of marine Nd cycling with larger and more heterogenous surface and near-surface Nd sources and a greater dominance of reversible scavenging. Future work should explore employing a more horizontally nuanced benthic flux tied to local environmental and sedimentary conditions, in line with recent modelling studies (Pöppelmeier et al., 2020a; Pasquier et al., 2022), to introduce more spatial patterning in simulated [Nd]d.

Certainly, these sensitivity simulations demonstrate that (at least under a spatially uniform benthic flux) the bulk εNd of seafloor detrital sediment (Robinson et al., 2021) cannot be considered fully representative of the εNd composition of the sediment that is interacting with seawater in all instances. They highlight the need for observational and experimental quantification of the broad mobile Nd phases globally and their εNd signal, as well as constraints on the spatial distribution of such a benthic flux (e.g. identifying where and under what environmental conditions a benthic flux occurs and at what strength). The model's response to fsed sets the stage for further testing of global sedimentary εNd on marine Nd cycling, providing the foundation for resolving the inherent complex multitude of processes.

4 Summary and conclusions

In this study, we describe the implementation of Nd isotopes (143Nd and 144Nd) into the ocean component of the FAMOUS GCM (Nd v1.0), providing a powerful tool designed for comprehensively exploring global marine Nd cycling, especially through representing explicit non-conservative processes to explore the extent of their influence on global seawater Nd distributions. Our Nd isotope scheme starts from previous Nd isotope implementations (Rempfer et al., 2011; Gu et al., 2019; Pöppelmeier et al., 2020a; Arsouze et al., 2009; Siddall et al., 2008) but revisits and updates Nd sources, sinks, and tracer transformation in line with increased observations and recent findings relating to global marine Nd cycling.

Model sensitivity to reversible scavenging efficiency demonstrates its importance for determining the increase in Nd concentration along the circulation pathway and acts to enhance regional basinal gradients in simulated εNd by maintaining the localised provenance signal. On the other hand, a widespread seafloor sedimentary flux presents a major deep-ocean source of Nd, the magnitude of which governs horizontal seafloor Nd concentrations across ocean basins. In the deep North Pacific, simulated εNd is too non-radiogenic compared to measurements, highlighting that the εNd of the sediment flux, as captured by the bulk εNd, is not a true representation of the highly radiogenic labile sediment phase interacting with seawater. Moreover, authigenic precipitation may pose a missing effective Nd sink. Model–data mismatch at the mouth of the Labrador Sea suggests that detrital contributions from highly non-radiogenic sediments disproportionally govern particle–seawater interactions. Overall, our results demonstrate that more precise constraints are needed to resolve the complex benthic processes that would drive spatially elevated and reactivity-weighted sediment fluxes.

Exploring the behaviour of simulated [Nd]d and εNd distributions in detail also highlighted some of the structural limitations of the model (e.g. difficulties representing highly localised and surface features) and influential biases in the physical ocean circulation (e.g. limited northward intrusion of AABW in the North Atlantic).

This study provides the groundwork for a future comprehensive optimisation of the marine Nd isotope scheme in FAMOUS. In the first instance, we suggest calibration of the key tuning parameters ([Nd]p/[Nd]d and fsed) to best represent modern seawater measurements. Additionally, it would be beneficial to obtain additional observational constraints on the broad labile sediment εNd interacting with seawater across different seafloor regions, including constraining the exchange between authigenic and detrital sediment phases during early diagenesis. Future sensitivity studies could also focus on the influence of river particulate and continental marginal sources on marine Nd to provide further insight into (and possibly constrain) the relative importance of these inputs, as opposed to a predominantly benthic seafloor-wide source.

This new model scheme can aid in the delivery of more robust applications of εNd as a modern and palaeo-tracer. It provides a platform for dynamic modelling under different modes of marine Nd cycling or varied climatic and oceanographic conditions, enabling current hypotheses to be tested rigorously with the aim of constraining Nd cycling under a complex and partially understood marine geochemical system.

Code availability

The code detailing the advances for the marine Nd isotope scheme (Nd v1.0) described in this paper is available via the Research Data Leeds Repository (Robinson et al., 2022: under a Creative Commons Attribution 4.0 International (CC-BY 4.0) license. These files are known as code modification (i.e. “mod”) files and should be applied to the original FAMOUS model code, which is protected under UK Crown Copyright and can be obtained from the National Centre for Atmospheric Science (NCAS) Computational Modelling Services (CMS;; last access: 12 February 2013​​​​​​​), with specific FAMOUS documentation available at (NCAS, 2023). All files and corresponding information that needs to be applied to configure each individual simulation presented here are available from the same DOI (; note that to complete the setup of these simulations, line 4909 in each simulation's tracer.f file needs updating with the corresponding RS_TUNE value as listed in Table 3. A complete version of the modified code for the EXPT_RS4 simulation using FAMOUS–MOSES1, including Nd isotope implementation, is archived at the Research Data Leeds Repository and linked from the DOI above.

Control candidate simulations for new FAMOUS reference are as follows:

  • XPDAA control simulation (0–5000 years),

  • XPDAB control candidate simulation (0–5000 years),

  • XPDAC control candidate simulation (0–5000 years),

  • XPDEA control candidate simulation (0–5000 years).

Reversible scavenging efficiency ([Nd]p/[Nd]d) sensitivity simulations are as follows:

  • XPDAI [Nd]p/[Nd]d=0.001 (0–9000 years),

  • XPDAD [Nd]p/[Nd]d=0.002 (0–9000 years),

  • XPDAH [Nd]p/[Nd]d=0.003 (0–9000 years),

  • XPDAE [Nd]p/[Nd]d=0.004 (0–9000 years),

  • XPDAF [Nd]p/[Nd]d=0.005 (0–9000 years),

  • XPDAG [Nd]p/[Nd]d=0.006 (0–9000 years).

Total Nd source from sediment (fsed) sensitivity simulations are as follows:

  • XPDAL fsed=1.5 Gg yr−1 (0–9000 years),

  • XPDAM fsed=3.0 Gg yr−1 (0–9000 years),

  • XPDAH fsed=4.5 Gg yr−1 (0–9000 years),

  • XPDAN fsed=6.0 Gg yr−1 (0–9000 years).

Data availability

The data are available via the Research Data Leeds Repository (Robinson et al., 2022:


The Supplement related to this article is to be associated with a DOI and is currently attached as a .zip, containing the following files. Supplementary Information. Documented supplementary text, figures, and tables (Sect. S1, Tables S1–S2, Figs. S1–S11). Table S3. A spreadsheet of seawater Nd concentration and isotope measurements and references used to validate model scheme. simulation_files/; a folder with all of the simulation files needed to run each simulation. standard_famous_mods/. A folder containing all standard FAMOUS GCM modification files. The supplement related to this article is available online at:

Author contributions

RFI designed the project with SR and supervised its completion. SR wrote and implemented the code with technical input from JT, RFI, and LJG. SR completed the experiment design, ran the simulations, analysed results, and prepared the manuscript with input from all co-authors. All co-authors contributed scientific expertise throughout. PJV designed and ran the pre-industrial FAMOUS GCM perturbed physics ensemble. FP updated the dust source εNd boundary conditions. YP provided the seawater REE compilation, which was updated by SR.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Suzanne Robinson was funded by the Natural Environment Research Council (NERC) SPHERES Doctoral Training Partnership (grant no. NE/ L002574/1). Lauren J. Gregoire was funded by UKRI Future Leaders Fellowship “SMBGen” (grant no. MR/S016961/1). Julia Tindall was supported through the Centre for Environmental Modelling and Computation (CEMAC) at the University of Leeds. Frerk Pöppelmeier was supported by the European Union's Horizon 2020 programme (grant no. 101023443). Tina van de Flierdt acknowledges funding from NERC (grant no. NE/P019080/1). Our heartfelt thanks go to Jamie D. Wilson for valuable discussions when implementing the reversible scavenging scheme and to Andy Ridgwell for helpful discussions during the design of code developments. We would like to thank all reviewers for taking the time and effort necessary to review the manuscript. We sincerely appreciate all valuable comments and suggestions, which helped us to improve the quality of the manuscript.

Financial support

This research has been supported by the Natural Environment Research Council (grant no. NE/L002574/1).

Review statement

This paper was edited by Andrew Yool and reviewed by Benoit Pasquier, Ed Hathorne, Catherine Jeandel, and Torben Stichel.


Abbott, A. N.: A benthic flux from calcareous sediments results in non-conservative neodymium behavior during lateral transport: A study from the Tasman Sea, Geology, 47, 363–366,, 2019. 

Abbott, A. N., Haley, B. A., and McManus, J.: Bottoms up: Sedimentary control of the deep North Pacific Ocean's εNd signature, Geology, 43, 1035–1038,, 2015a. 

Abbott, A. N., Haley, B. A., McManus, J., and Reimers, C. E.: The sedimentary flux of dissolved rare earth elements to the ocean, Geochim. Cosmochim. Ac., 154, 186–200,, 2015b. 

Abbott, A. N., Löhr, S., and Trethewy, M.: Are clay minerals the primary control on the oceanic rare earth element budget?, Front. Mar. Sci., 6, 2296–7745,, 2019. 

Anderson, S. P.: Glaciers show direct linkage between erosion rate and chemical weathering fluxes, Geomorphology, 67, 147–157,, 2005. 

Arsouze, T., Dutay, J.-C. C., Lacan, F., and Jeandel, C.: Modeling the neodymium isotopic composition with a global ocean circulation model, Chem. Geol., 239, 165–177,, 2007. 

Arsouze, T., Dutay, J.-C., Lacan, F., and Jeandel, C.: Reconstructing the Nd oceanic cycle using a coupled dynamical – biogeochemical model, Biogeosciences, 6, 2829–2846,, 2009. 

Ayache, M., Dutay, J.-C., Arsouze, T., Révillon, S., Beuvier, J., and Jeandel, C.: High-resolution neodymium characterization along the Mediterranean margins and modelling of εNd distribution in the Mediterranean basins, Biogeosciences, 13, 5259–5276,, 2016. 

Ayache, M., Dutay, J.-C., Tachikawa, K., Arsouze, T., and Jeandel, C.: Neodymium budget in the Mediterranean Sea: evaluating the role of atmospheric dusts using a high-resolution dynamical-biogeochemical model, Biogeosciences, 20, 205–227,, 2023. 

Bacon, M. P. and Anderson, R. F.: Distribution of Thorium Isotopes Between Dissolved and Particulate Forms in the Deep Sea, J. Geophys. Res., 87, 2045–2056,, 1982. 

Bayon, G., German, C. R., Burton, K. W., Nesbitt, R. W., and Rogers, N.: Sedimentary Fe-Mn oxyhydroxides as paleoceanographic archives and the role of aeolian flux in regulating oceanic dissolved REE, Earth Planet. Sc. Lett., 224, 477–492,, 2004. 

Behrenfeld, M. J. and Falkowski, P. G.: Photosynthetic rates derived from satellite-based chlorophyll concentration, Limnol. Oceanogr., 1–20,, 1997. 

Bertram, C. J. and Elderfield, H.: The geochemical balance of the rare earth elements and neodymium isotopes in the oceans, Geochim. Cosmochim. Ac., 57, 1957–1986,, 1993. 

Blanchet, C. L.: A database of marine and terrestrial radiogenic Nd and Sr isotopes for tracing earth-surface processes, Earth Syst. Sci. Data, 11, 741–759,, 2019. 

Blaser, P., Frank, M., and van de Flierdt, T.: Revealing past ocean circulation with neodymium isotopes, Past Glob. Chang. Mag., 27, 54–55,, 2019a. 

Blaser, P., Pöppelmeier, F., Schulz, H., Gutjahr, M., Frank, M., Lippold, J., Heinrich, H., Link, J. M., Hoffmann, J., Szidat, S., and Frank, N.: The resilience and sensitivity of Northeast Atlantic deep water εNd to overprinting by detrital fluxes over the past 30,000 years, Geochim. Cosmochim. Ac., 245, 79–97,, 2019b. 

Blaser, P., Gutjahr, M., Pöppelmeier, F., Frank, M., Kaboth-Bahr, S., and Lippold, J.: Labrador Sea bottom water provenance and REE exchange during the past 35,000 years, Earth Planet. Sc. Lett., 542, 116299,, 2020. 

Chamecki, M., Chor, T., Yang, D., and Meneveau, C.: Material Transport in the Ocean Mixed Layer: Recent Developments Enabled by Large Eddy Simulations, Rev. Geophys., 57, 1338–1371,, 2019. 

Chavagnac, V., Saleban Ali, H., Jeandel, C., Leleu, T., Destrigneville, C., Castillo, A., Cotte, L., Waeles, M., Cathalot, C., Laes-Huon, A., Pelleter, E., Nonnotte, P., Sarradin, P. M., and Cannat, M.: Sulfate minerals control dissolved rare earth element flux and Nd isotope signature of buoyant hydrothermal plume (EMSO-Azores, 37 N Mid-Atlantic Ridge), Chem. Geol., 499, 111–125,, 2018. 

Chen, T. Y., Stumpf, R., Frank, M., Bełdowski, J., and Staubwasser, M.: Contrasting geochemical cycling of hafnium and neodymium in the central Baltic Sea, Geochim. Cosmochim. Ac., 123, 166–180,, 2013. 

Collins, W. J., Bellouin, N., Doutriaux-Boucher, M., Gedney, N., Halloran, P., Hinton, T., Hughes, J., Jones, C. D., Joshi, M., Liddicoat, S., Martin, G., O'Connor, F., Rae, J., Senior, C., Sitch, S., Totterdell, I., Wiltshire, A., and Woodward, S.: Development and evaluation of an Earth-System model – HadGEM2, Geosci. Model Dev., 4, 1051–1075,, 2011. 

Cox, P. M., Betts, R. A., Bunton, C. B., Essery, R. L. H., Rowntree, P. R., and Smith, J.: The impact of new land surface physics on the GCM simulation of climate and climate sensitivity, Clim. Dynam., 15, 183–203,, 1999. 

Crossley, J. F. and Roberts, D. L.: The Thermodynamic/Dynamic Sea Ice Model, Unified Model Doc. Pap. No 45, 1995. 

Dai, A. and Trenberth, K. E.: Estimates of Freshwater Discharge from Continents: Latitudinal and Seasonal Variations, J. Hydrometeorol., 3, 660–687, 2002. 

Deng, K., Yang, S., Du, J., Lian, E., and Vance, D.: Dominance of benthic flux of REEs on continental shelves: implications for oceanic budgets, Geochemical Perspect. Lett., 22, 26–30,, 2022. 

Dentith, J. E., Ivanovic, R. F., Gregoire, L. J., Tindall, J. C., and Smith, R. S.: Ocean circulation drifts in multi-millennial climate simulations: the role of salinity corrections and climate feedbacks, Clim. Dynam., 52, 1761–1781,, 2019. 

Dentith, J. E., Ivanovic, R. F., Gregoire, L. J., Tindall, J. C., and Robinson, L. F.: Simulating stable carbon isotopes in the ocean component of the FAMOUS general circulation model with MOSES1 (XOAVI), Geosci. Model Dev., 13, 3529–3552,, 2020. 

Doney, S. C., Lindsay, K., Caldeira, K., Campin, J. M., Drange, H., Dutay, J. C., Follows, M., Gao, Y., Gnanadesikan, A., Gruber, N., Ishida, A., Joos, F., Madec, G., Maier-Reimer, E., Marshall, J. C., Matear, R. J., Monfray, P., Mouchet, A., Najjar, R., Orr, J. C., Plattner, G. K., Sarmiento, J., Schlitzer, R., Slater, R., Totterdell, I. J., Weirig, M. F., Yamanaka, Y., and Yool, A.: Evaluating global ocean carbon models: The importance of realistic physics, Global Biogeochem. Cy., 18, 3017,, 2004. 

Du, J., Haley, B. A., and Mix, A. C.: Neodymium isotopes in authigenic phases, bottom waters and detrital sediments in the Gulf of Alaska and their implications for paleo-circulation reconstruction, Geochim. Cosmochim. Ac., 193, 14–35,, 2016. 

Du, J., Haley, B. A., Mix, A. C., Walczak, M. H., and Praetorius, S. K.: Flushing of the deep Pacific Ocean and the deglacial rise of atmospheric CO2 concentrations, Nat. Geosci., 11, 749–755,, 2018. 

Du, J., Haley, B. A., and Mix, A. C.: Evolution of the Global Overturning Circulation since the Last Glacial Maximum based on marine authigenic neodymium isotopes, Quaternary Sci. Rev., 241, 106396,, 2020. 

Du, J., Haley, B. A., Mix, A. C., Abbott, A. N., McManus, J., and Vance, D.: Reactive-transport modeling of neodymium and its radiogenic isotope in deep-sea sediments: The roles of authigenesis, marine silicate weathering and reverse weathering, Earth Planet. Sc. Lett., 596, 117792,, 2022. 

Dunne, J. P., Sarmiento, J. L., and Gnanadesikan, A.: A synthesis of global particle export from the surface ocean and cycling through the ocean interior and on the seafloor, Global Biogeochem. Cy., 21, GB4006,, 2007. 

Dunne, J. P., Hales, B., and Toggweiler, J. R.: Global calcite cycling constrained by sediment preservation controls, Global Biogeochem. Cy., 26, GB3023,, 2012. 

Dutay, J. C., Lacan, F., Roy-Barman, M., and Bopp, L.: Influence of particle size and type on 231Pa and 230Th simulation with a global coupled biogeochemical-ocean general circulation model: A first approach, Geochem. Geophy. Geosy., 10, GB3023,, 2009. 

Elderfield, H., Upstill-Goddard, R., and Sholkovitz, E. R.: The rare earth elements in rivers, estuaries, and coastal seas and their significance to the composition of ocean waters, Geochim. Cosmochim. Ac., 54, 971–991,, 1990. 

Essery, R. and Clark, D. B.: Developments in the MOSES 2 land-surface model for PILPS 2e, Glob. Planet. Change, 38, 161–164,, 2003. 

Essery, R. L. H., Best, M. J., and Cox, P. M.: MOSES 2.2 technical documentation, Hadley Centre Technical Note 30, (last access: 15 February 2023), 2001. 

Ferreira, D., Cessi, P., Coxall, H. K., De Boer, A., Dijkstra, H. A., Drijfhout, S. S., Eldevik, T., Harnik, N., McManus, J. F., Marshall, D. P., Nilsson, J., Roquet, F., Schneider, T., and Wills, R. C.: Atlantic-Pacific Asymmetry in Deep Water Formation, Annu. Rev. Earth Planet. Sc., 46, 327–352,, 2018. 

Frajka-Williams, E., Ansorge, I. J., Baehr, J., Bryden, H. L., Chidichimo, M. P., Cunningham, S. A., Danabasoglu, G., Dong, S., Donohue, K. A., Elipot, S., Heimbach, P., Holliday, N. P., Hummels, R., Jackson, L. C., Karstensen, J., Lankhorst, M., Le Bras, I. A., Susan Lozier, M., McDonagh, E. L., Meinen, C. S., Mercier, H., Moat, B. I., Perez, R. C., Piecuch, C. G., Rhein, M., Srokosz, M. A., Trenberth, K. E., Bacon, S., Forget, G., Goni, G., Kieke, D., Koelling, J., Lamont, T., McCarthy, G. D., Mertens, C., Send, U., Smeed, D. A., Speich, S., van den Berg, M., Volkov, D., and Wilson, C.: Atlantic meridional overturning circulation: Observed transport and variability, Front. Mar. Sci., 6, 260,, 2019. 

Frank, M.: Radiogenic isotopes: Tracers of past ocean circulation and erosional input, Rev. Geophys., 40, 1001,, 2002. 

GEOTRACES Intermediate Data Product Group: The GEOTRACES Intermediate Data Product 2021 (IDP2021), NERC EDS British Oceanographic Data Centre NOC [data set],, 2021. 

German, C. R., Klinkhammer, G. P., Edmond, J. M., Mitrat, A., and H.Elderfield: Hydrothermal scavenging of rare-earth elements in the ocean, Nature, 345, 516–518, 1990. 

Goldstein, S. J. and Jacobsen, S. B.: The Nd and Sr isotopic systematics of river-water dissolved material: Implications for the sources of Nd and Sr in seawater, Chemical Geology: Isotope Geoscience Section, 66, 245–272,, 1987. 

Goldstein, S. L. and Hemming, S. R.: Long-lived Isotopic Tracers in Oceanography, Paleoceanography, and Ice-Sheet Dynamics, Treatise of Geochemistry, 6, 453–483,, 2003. 

Goldstein, S. L. and O'Nions, R. K.: Nd and Sr isotopic relationships in pelagic clays and gerromanganese deposits, Nature, 292, 324–327, 1981. 

Goldstein, S. L., O'Nions, R. K., and Hamilton, P. J.: A Sm-Nd isotopic study of atmospheric dusts and particulates from major river systems, Earth Planet. Sc. Lett., 70, 221–236,, 1984. 

Gordon, C., Cooper, C., Senior, C. A., Banks, H., Gregory, J. M., Johns, T. C., Mitchell, J. F. B., and Wood, R. A.: The simulation of SST, sea ice extents and ocean heat transports in a version of the Hadley Centre coupled model without flux adjustments, Clim. Dynam., 16, 147–168,, 2000. 

Greaves, M. J., Statham, P. J., and Elderfield, H.: Rare earth element mobilization from marine atmospheric dust into seawater, Mar. Chem., 46, 255–260,, 1994. 

Gregoire, L. J., Valdes, P. J., Payne, A. J., and Kahana, R.: Optimal tuning of a GCM using modern and glacial constraints, Clim. Dynam., 37, 705–719,, 2011. 

Gregoire, L. J., Payne, A. J., and Valdes, P. J.: Deglacial rapid sea level rises caused by ice-sheet saddle collapses, Nature, 487, 219–222,, 2012. 

Gregoire, L. J., Valdes, P. J., and Payne, A. J.: The relative contribution of orbital forcing and greenhouse gases to the North American deglaciation, Geophys. Res. Lett., 42, 9970–9979,, 2015. 

Gregoire, L. J., Otto-Bliesner, B., Valdes, P. J., and Ivanovic, R.: Abrupt Bølling warming and ice saddle collapse contributions to the Meltwater Pulse 1a rapid sea level rise, Geophys. Res. Lett., 43, 9130–9137,, 2016. 

Grenier, M., Jeandel, C., Lacan, F., Vance, D., Venchiarutti, C., Cros, A., and Cravatte, S.: From the subtropics to the central equatorial Pacific Ocean: Neodymium isotopic composition and rare earth element concentration variations, J. Geophys. Res.-Oceans, 118, 592–618,, 2013. 

Grousset, F. E., Biscaye, P. E., Zindler, A., Prospero, J., and Chester, R.: Neodymium isotopes as tracers in marine sediments and aerosols: North Atlantic, Earth Planet. Sc. Lett., 87, 367–378,, 1988. 

Grousset, F. E., Parra, M., Bory, A., Martinez, P., Bertrand, P., Shimmield-M, G., and Ellamn, R. M.: Saharan Wind Regimes Traced by the Sr-Nd Isotopic Composition of Subtropical Atlantic Sediments: Last Glacial Maximum vs Today, Quaternary Sci. Rev., 17, 395–409,, 1998. 

Gu, S., Liu, Z., Jahn, A., Rempfer, J., Zhang, J., and Joos, F.: Modeling Neodymium Isotopes in the Ocean Component of the Community Earth System Model (CESM1), J. Adv. Model. Earth Syst., 11, 624–640,, 2019. 

Haley, B. A., Klinkhammer, G. P., and McManus, J.: Rare earth elements in pore waters of marine sediments, Geochim. Cosmochim. Ac., 68, 1265–1279,, 2004. 

Haley, B. A., Frank, M., Hathorne, E., and Pisias, N.: Biogeochemical implications from dissolved rare earth element and Nd isotope distributions in the Gulf of Alaska, Geochim. Cosmochim. Ac., 126, 455–474,, 2014. 

Haley, B. A., Du, J., Abbott, A. N., and McManus, J.: The Impact of Benthic Processes on Rare Earth Element and Neodymium Isotope Distributions in the Oceans, Front. Mar. Sci., 4, 426,, 2017. 

Henderson, G. M., Heinze, C., Anderson, R. F., and Winguth, A. M. E.: Global distribution of the 230Th flux to ocean sediments constrained by GCM modelling, Deep-Res. Pt. I, 46, 1861–1893,, 1999. 

Hopcroft, P. O. and Valdes, P. J.: Last glacial maximum constraints on the Earth System model HadGEM2-ES, Clim. Dynam., 45, 1657–1672,, 2015. 

Hopcroft, P. O., Valdes, P. J., Woodward, S., and Joshi, M. M.: Last glacial maximum radiative forcing from mineral dust aerosols in an Earth system model, J. Geophys. Res., 120, 8186–8205,, 2015. 

Jacobsen, S. B. and Wasserburg, G. J.: Sm-Nd evolution of chondrites, Earth Planet. Sc. Lett., 50, 139–155,, 1980. 

Jeandel, C.: Concentration and isotopic composition of Nd in the South Atlantic Ocean, Earth Planet. Sc. Lett., 117, 581–591,, 1993. 

Jeandel, C., Bishop, J. K., and Zindler, A.: Exchange of neodymium and its isotopes between seawater and small and large particles in the Sargasso Sea, Geochim. Cosmochim. Ac., 59, 535–547,, 1995. 

Johannesson, K. H. and Burdige, D. J.: Balancing the global oceanic neodymium budget: Evaluating the role of groundwater, Earth Planet. Sc. Lett., 253, 129–142,, 2007. 

Jones, C.: A fast ocean GCM without flux adjustments, J. Atmos. Ocean. Technol., 20, 1857–1868,<1857:AFOGWF>2.0.CO;2, 2003. 

Jones, C., Gregory, J., Thorpe, R., Cox, P., Murphy, J., Sexton, D., and Valdes, P.: Systematic optimisation and climate simulation of FAMOUS, a fast version of HadCM3, Clim. Dynam., 25, 189–204,, 2005. 

Jones, K. M., Khatiwala, S. P., Goldstein, S. L., Hemming, S. R., and van de Flierdt, T.: Modeling the distribution of Nd isotopes in the oceans using an ocean general circulation model, Earth Planet. Sc. Lett., 272, 610–619,, 2008. 

Kato, Y., Fujinaga, K., Nakamura, K., Takaya, Y., Kitamura, K., Ohta, J., Toda, R., Nakashima, T., and Iwamori, H.: Deep-sea mud in the Pacific Ocean as a potential resource for rare-earth elements, Nat. Geosci., 4, 535–539,, 2011. 

Kriest, I. and Oschlies, A.: On the treatment of particulate organic matter sinking in large-scale models of marine biogeochemical cycles, Biogeosciences, 5, 55–72,, 2008. 

Kuhlbrodt, T., Griesel, A., Montoya, M., Levermann, A., Hofmann, M., and Rahmstorf, S.: On the driving processes of the Atlantic meridional overturning circulation, Rev. Geophys., 45, 2001,, 2007. 

Lacan, F. and Jeandel, C.: Neodymium isotopes as a new tool for quantifying exchange fluxes at the continent-ocean interface, Earth Planet. Sc. Lett., 232, 245–257,, 2005. 

Lagarde, M., Lemaitre, N., Planquette, H., Grenier, M., Belhadj, M., Lherminier, P., and Jeandel, C.: Particulate rare earth element behavior in the North Atlantic (GEOVIDE cruise), Biogeosciences, 17, 5539–5561,, 2020. 

Lambelet, M., van de Flierdt, T., Crocket, K., Rehkämper, M., Kreissig, K., Coles, B., Rijkenberg, M. J. A., Gerringa, L. J. A., de Baar, H. J. W., and Steinfeldt, R.: Neodymium isotopic composition and concentration in the western North Atlantic Ocean: Results from the GEOTRACES GA02 section, Geochim. Cosmochim. Ac., 177, 1–29,, 2016. 

Lambelet, M., van de Flierdt, T., Butler, E. C. V., Bowie, A. R., Rintoul, S. R., Watson, R. J., Remenyi, T., Lannuzel, D., Warner, M., Robinson, L. F., Bostock, H. C., and Bradtmiller, L. I.: The Neodymium Isotope Fingerprint of Adélie Coast Bottom Water, Geophys. Res. Lett., 45, 11247–11256,, 2018. 

Laws, E. A., Falkowski, P. G., Smith, W. O., Ducklow, H., and McCarthy, J. J.: Temperature effects on export production in the open ocean, Global Biogeochem. Cy., 14, 1231–1246,, 2000. 

Locarnini, R. A., Mishonov, A. V., Baranova, O. K., Boyer, T. P., Zweng, M. M., Garcia, H. E., Reagan, J. R., Seidov, D., Weathers, K., Paver, C. R., and Smolyar, I.: World Ocean Atlas 2018, Volume 1: Temperature, edited by: Mishonov, A., NOAA Atlas NESDIS 81, 52 pp., 2019. 

Mahowald, N. M., Muhs, D. R., Levis, S., Rasch, P. J., Yoshioka, M., Zender, C. S., and Luo, C.: Change in atmospheric mineral aerosols in response to climate: Last glacial period, preindustrial, modern, and doubled carbon dioxide climates, J. Geophys. Res.-Atmos., 111, D10202,, 2006. 

Maier-Reimer, E.: Geochemical cycles in an ocean general circulation model. Preindustrial tracer distributions, Global Biogeochem. Cy., 7, 645–677,, 1993. 

Martin, J. H., Knauer, G. A., Karl, D. M., and Broenkow, W. W.: VERTEX: carbon cycling in the northeast Pacific, Deep-Sea Res. Pt. I, 34, 267–285,, 1987. 

McCarthy, G. D., Smeed, D. A., Johns, W. E., Frajka-Williams, E., Moat, B. I., Rayner, D., Baringer, M. O., Meinen, C. S., Collins, J., and Bryden, H. L.: Measuring the Atlantic Meridional Overturning Circulation at 26 N, Prog. Oceanogr., 130, 91–111,, 2015. 

Najjar, R. and Orr, J. C.: Design of OCMIP-2 simulations of chlorofluorocarbons, the solubility pump and common biogeochemistry, Internal OCMIP Report, LSCE/CEA Saclay, Gif-sur-Yvette, France 1998. 

National Centre for Atmospheric Science (NCAS): Standalone UM system and FAMOUS [code],, last access: 15t February 2023. 

Noh, Y., Kang, I. S., Herold, M., and Raasch, S.: Large eddy simulation of particle settling in the ocean mixed layer, Phys. Fluids, 18, 85109,, 2006. 

Oka, A., Hasumi, H., Obata, H., Gamo, T., and Yamanaka, Y.: Study on vertical profiles of rare earth elements by using an ocean general circulation model, Global Biogeochem. Cy., 23, GB4025,, 2009. 

Oka, A., Tazoe, H., and Obata, H.: Simulation of global distribution of rare earth elements in the ocean using an ocean general circulation model, J. Oceanogr., 1, 413–430,, 2021. 

Osborne, A. H., Haley, B. A., Hathorne, E. C., Plancherel, Y., and Frank, M.: Rare earth element distribution in Caribbean seawater: Continental inputs versus lateral transport of distinct REE compositions in subsurface water masses, Mar. Chem., 177, 172–183,, 2015. 

Osborne, A. H., Hathorne, E. C., Schijf., J., Plancherel, Y., Boning, P., and Frank, M.: The potential of sedimentary foraminiferal rare earth element patterns to trace water masses in the past, Geochem. Geophy. Geosy., 18, 1550-1568,, 2017. 

Paffrath, R., Pahnke, K., Böning, P., Rutgers van der Loeff, M., Valk, O., Gdaniec, S., and Planquette, H.: Seawater-Particle Interactions of Rare Earth Elements and Neodymium Isotopes in the Deep Central Arctic Ocean, J. Geophys. Res.-Oceans, 126, e2021JC017423,, 2021. 

Palmer, J. R. and Totterdell, I. J.: Production and export in a global ocean ecosystem model, Deep. Res.-Pt. I, 48, 1169–1198,, 2001. 

Pasquier, B., Hines, S. K. V., Liang, H., Wu, Y., Goldstein, S. L., and John, S. G.: GNOM v1.0: an optimized steady-state model of the modern marine neodymium cycle, Geosci. Model Dev., 15, 4625–4656,, 2022. 

Paul, S. A. L., Haeckel, M., Bau, M., Bajracharya, R., and Koschinsky, A.: Small-scale heterogeneity of trace metals including rare earth elements and yttrium in deep-sea sediments and porewaters of the Peru Basin, southeastern equatorial Pacific, Biogeosciences, 16, 4829–4849,, 2019. 

Pearce, C. R., Jones, M. T., Oelkers, E. H., Pradoux, C., and Jeandel, C.: The effect of particulate dissolution on the neodymium (Nd) isotope and Rare Earth Element (REE) composition of seawater, Earth Planet. Sc. Lett., 369–370, 138–147,, 2013. 

Pöppelmeier, F., Scheen, J., Blaser, P., Lippold, J., Gutjahr, M., and Stocker, T. F.: Influence of Elevated Nd Fluxes on the Northern Nd Isotope End Member of the Atlantic During the Early Holocene, Paleoceanogr. Paleoclimatology, 35, e2020PA003973,, 2020a. 

Pöppelmeier, F., Blaser, P., Gutjahr, M., Jaccard, S. L., Frank, M., Max, L., and Lippold, J.: Northern Sourced Water dominated the Atlantic Ocean during the Last Glacial Maximum, Geology, 48, 826–829,, 2020b. 

Pöppelmeier, F., Janssen, D. J., Jaccard, S. L., and Stocker, T. F.: Modeling the marine chromium cycle: new constraints on global-scale processes, Biogeosciences, 18, 5447–5463,, 2021a. 

Pöppelmeier, F., Scheen, J., Jeltsch-Thömmes, A., and Stocker, T. F.: Simulated stability of the Atlantic Meridional Overturning Circulation during the Last Glacial Maximum, Clim. Past, 17, 615–632,, 2021b. 

Pöppelmeier, F., Lippold, J., Blaser, P., Gutjahr, M., Frank, M., and Stocker, T. F.: Neodymium isotopes as a paleo-water mass tracer: A model-data reassessment, Quaternary Sci. Rev., 279, 107404,, 2022. 

Rahlf, P., Hathorne, E., Laukert, G., Gutjahr, M., Weldeab, S., and Frank, M.: Tracing water mass mixing and continental inputs in the southeastern Atlantic Ocean with dissolved neodymium isotopes, Earth Planet. Sc. Lett., 530, 115944,, 2020. 

Rahlf, P., Laukert, G., Hathorne, E. C., Vieira, L. H., and Frank, M.: Dissolved neodymium and hafnium isotopes and rare earth elements in the Congo River Plume: Tracing and quantifying continental inputs into the southeast Atlantic, Geochim. Cosmochim. Ac., 294, 192–214,, 2021. 

Rempfer, J., Stocker, T. F., Joos, F., Dutay, J.-C., and Siddall, M.: Modelling Nd-isotopes with a coarse resolution ocean circulation model: Sensitivities to model parameters and source/sink distributions, Geochim. Cosmochim. Ac., 75, 5927–5950,, 2011. 

Roberts, N. L., Piotrowski, A. M., McManus, J. F., and Keigwin, L. D.: Synchronous deglacial overturning and water mass source changes, Science, 327, 75–78,, 2010. 

Robinson, S., Ivanovic, R., van de Flierdt, T., Blanchet, C. L., Tachikawa, K., Martin, E. E., Cook, C. P., Williams, T., Gregoire, L., Plancherel, Y., Jeandel, C., and Arsouze, T.: Global continental and marine detrital εNd: An updated compilation for use in understanding marine Nd cycling, Chem. Geol., 567, 120119,, 2021. 

Robinson, S., Ivanovic, R., Gregoire, L., Tindall, J., van de Flierdt, T., Plancherel, Y., Pöppelmeier, F., Tachikawa, K., and Valdes, P.: Simulating marine neodymium isotope distributions using Nd v1.0 coupled to the ocean component of the FAMOUS-MOSES1 climate model: sensitivities to reversible scavenging efficiency and benthic source distributions, University of Leeds [data set],, 2022. 

Rousseau, T. C. C., Sonke, J. E., Chmeleff, J., Van Beek, P., Souhaut, M., Boaventura, G., Seyler, P., and Jeandel, C.: Rapid neodymium release to marine waters from lithogenic sediments in the Amazon estuary, Nat. Commun., 6, 7592,, 2015. 

Sarmiento, J. L. and Le Quéré, C.: Oceanic carbon dioxide uptake in a model of century-scale global warming, Science, 274, 1346–1350,, 1996. 

Sherriff-Tadano, S., Abe-Ouchi, A., and Oka, A.: Impact of mid-glacial ice sheets on deep ocean circulation and global climate, Clim. Past, 17, 95–110,, 2021. 

Sholkovitz, E. R.: The geochemistry of rare earth elements in the Amazon River estuary, Geochim. Cosmochim. Ac., 57, 2181–2190,, 1993. 

Siddall, M., Khatiwala, S., van de Flierdt, T., Jones, K., Goldstein, S. L., Hemming, S., and Anderson, R. F.: Towards explaining the Nd paradox using reversible scavenging in an ocean general circulation model, Earth Planet. Sc. Lett., 274, 448–461,, 2008. 

Sinha, B., Smeed, D. A., McCarthy, G., Moat, B. I., Josey, S. A., Hirschi, J. J. M., Frajka-Williams, E., Blaker, A. T., Rayner, D., and Madec, G.: The accuracy of estimates of the overturning circulation from basin-wide mooring arrays, Prog. Oceanogr., 160, 101–123,, 2018. 

Smith, R. N. B.: A scheme for predicting layer clouds and their water content in a general circulation model, Q. J. Roy. Meteor. Soc., 116, 435–460,, 1990. 

Smith, R. N. B.: Subsurface, surface and boundary layer processes. Unified Model Documentation Paper 24, The Met. Office, Bracknell, Berkshire UK, 64 pp., 1993. 

Smith, R. S.: The FAMOUS climate model (versions XFXWB and XFHCC): description update to version XDBUA, Geosci. Model Dev., 5, 269–276,, 2012. 

Smith, R. S. and Gregory, J. M.: A study of the sensitivity of ocean overturning circulation and climate to freshwater input in different regions of the North Atlantic, Geophys. Res. Lett., 36, L15701,, 2009. 

Smith, R. S., Gregory, J. M., and Osprey, A.: A description of the FAMOUS (version XDBUA) climate model and control run, Geosci. Model Dev., 1, 53–68,, 2008. 

Stichel, T., Kretschmer, S., Geibert, W., Lambelet, M., Plancherel, Y., Rutgers van der Loeff, M., and van de Flierdt, T.: Particle–Seawater Interaction of Neodymium in the North Atlantic, ACS Earth Sp. Chem., 4, 1700–1717,, 2020. 

Tachikawa, K., Jeandel, C., Vangriesheim, A., and Dupré, B.: Distribution of rare earth elements and neodymium isotopes in suspended particles of the tropical Atlantic Ocean (EUMELI site), Deep-Res. Pt. I, 46, 733–755,, 1999. 

Tachikawa, K., Athias, V., and Jeandel, C.: Neodymium budget in the modern ocean and paleo-oceanographic implications, J. Geophys. Res., 108, 3254,, 2003. 

Tachikawa, K., Arsouze, T., Bayon, G., Bory, A., Colin, C., Dutay, J. C., Frank, N., Giraud, X., Gourlan, A. T., Jeandel, C., Lacan, F., Meynadier, L., Montagna, P., Piotrowski, A. M., Plancherel, Y., Pucéat, E., Roy-Barman, M., and Waelbroeck, C.: The large-scale evolution of neodymium isotopic composition in the global modern and Holocene ocean revealed from seawater and archive data, Chem. Geol., 457, 131–148,, 2017. 

Talley, L. D., Pickard, G. L., Emery, W. J., and Swift, J. H.: Descriptive physical oceanography: An introduction: Sixth edition, Academic Press, ISBN 9780750645522,, 2011. 

Valdes, P. J., Armstrong, E., Badger, M. P. S., Bradshaw, C. D., Bragg, F., Crucifix, M., Davies-Barnard, T., Day, J. J., Farnsworth, A., Gordon, C., Hopcroft, P. O., Kennedy, A. T., Lord, N. S., Lunt, D. J., Marzocchi, A., Parry, L. M., Pope, V., Roberts, W. H. G., Stone, E. J., Tourte, G. J. L., and Williams, J. H. T.: The BRIDGE HadCM3 family of climate models: HadCM3@Bristol v1.0, Geosci. Model Dev., 10, 3715–3743,, 2017. 

van de Flierdt, T., Grifths, A. M., Lambelet, M., Little, S. H., Stichel, T., and Wilson, D. J.: Neodymium in the oceans: A global database, a regional comparison and implications for palaeoceanographic research, 374, 20150293,, 2016. 

Von Blanckenburg, F. and Nägler, T. F.: Weathering versus circulation-controlled changes in radiogenic isotope tracer composition of the Labrador Sea and North Atlantic Deep Water, Paleoceanography, 16, 424–434,, 2001. 

Wang, R., Clegg, J. A., Scott, P. M., Larkin, C. S., Deng, F., Thomas, A. L., Zheng, X. Y., and Piotrowski, A. M.: Reversible scavenging and advection – Resolving the neodymium paradox in the South Atlantic, Geochim. Cosmochim. Ac., 314, 121–139,, 2021. 

Williams, J. H. T., Smith, R. S., Valdes, P. J., Booth, B. B. B., and Osprey, A.: Optimising the FAMOUS climate model: inclusion of global carbon cycling, Geosci. Model Dev., 6, 141–160,, 2013. 

Wilson, D. J., Piotrowski, A. M., Galy, A., and Clegg, J. A.: Reactivity of neodymium carriers in deep sea sediments: Implications for boundary exchange and paleoceanography, Geochim. Cosmochim. Ac., 109, 197–221,, 2013. 

Woodward, S.: Mineral Dust in HadGEM2, Hadley Centre Technical Note 87 Mineral Dust in HadGEM 2 March 2011, (last access: 15 February 2023), 2011. 

Zhang, Y., Lacan, F., and Jeandel, C.: Dissolved rare earth elements tracing lithogenic inputs over the Kerguelen Plateau (Southern Ocean), Deep-Res. Pt. II, 55, 638–652,, 2008.  

Zweng, M. M., Reagan, J. R., Seidov, D., Boyer, T. P., Locarnini, R. A., Garcia, H. E., Mishonov, A. V., Baranova, O. K., Weathers, K., Paver, C. R., and Smolyar, I.: World Ocean Atlas 2018, Volume 2: Salinity, edited by: Mishonov, A., NOAA Atlas NESDIS 82, 50 pp., (last access: 15 February 2023​​​​​​​) 2019. 

Short summary
We present the implementation of neodymium (Nd) isotopes into the ocean model of FAMOUS (Nd v1.0). Nd fluxes from seafloor sediment and incorporation of Nd onto sinking particles represent the major global sources and sinks, respectively. However, model–data mismatch in the North Pacific and northern North Atlantic suggest that certain reactive components of the sediment interact the most with seawater. Our results are important for interpreting Nd isotopes in terms of ocean circulation.
Special issue