# 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, Tachikawa Kazuyo, Paul J. Valdes

► To cite this version:

Suzanne Robinson, Ruza F. Ivanovic, Lauren J. Gregoire, Julia Tindall, Tina van de Flierdt, et al. 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. Geoscientific Model Development, 2023, 16, pp.1231-1264. 10.5194/gmd-16-1231-2023. insu-04131495

HAL Id: insu-04131495

https://insu.hal.science/insu-04131495

Submitted on 16 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# 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<fnref n="1" />, Ruza F. Ivanovic<fnref n="1" />, Lauren J. Gregoire<fnref n="1" />, Julia Tindall<fnref n="1" />, Tina van de Flierdt<fnref n="2" />, Yves Plancherel<fnref n="2" />, Frerk Pöppelmeier<fnref n="3" />, Kazuyo Tachikawa<fnref n="4" />, and Paul J. Valdes<fnref n="5" />

**Correspondence:** Suzanne Robinson (ee14s2r@leeds.ac.uk)

Received: 6 July 2022 – Discussion started: 21 September 2022
Revised: 10 January 2023 – Accepted: 24 January 2023 – Published: 22 February 2023

**Abstract.** 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 ($\varepsilon_{Nd}$), presents an opportunity to develop a new scheme specifically designed to test these paradigms. Here, we present the implementation of Nd isotopes ($^{143}$Nd and $^{144}$Nd) 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 $\varepsilon_{Nd}$ distributions and provide an updated compilation of 6048 Nd concentrations and 3278 $\varepsilon_{Nd}$ measurements to assess model performance. Our findings support the notions that reversible scavenging is

a key process for enhancing the Atlantic–Pacific basinal $\varepsilon_{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 $\varepsilon_{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 $^{143}$Nd/$^{144}$Nd

<footnote n="1">¹School of Earth and Environment, University of Leeds, Leeds, LS2 9JT, UK</footnote>

<footnote n="2">²Department of Earth Science and Engineering, Imperial College London, London, SW7 2AZ, UK</footnote>

<footnote n="3">³Climate and Environmental Physics, Physics Institute,
University of Bern, 3012 Bern, Switzerland</footnote>

<footnote n="4">⁴Aix Marseille Univ., CNRS, IRD, INRAE, Coll France, CEREGE, Aix-en-Provence, 13545, France</footnote>

<footnote n="5">⁵School of Geographical Sciences, University of Bristol, Bristol, BS8 1RL, UK</footnote>
ratio is denoted relative to the bulk Earth standard:

$$ \varepsilon_{Nd} = \left( \frac{\left( \frac{^{143}\text{Nd}}{^{144}\text{Nd}} \right)_{\text{sample}}}{\left( \frac{^{143}\text{Nd}}{^{144}\text{Nd}} \right)_{\text{CHUR}}} - 1 \right) \times 10^4, \quad (1) $$

where $(^{143}\text{Nd}/^{144}\text{Nd})_{\text{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. $\delta^{13}\text{C}$, $\Delta^{14}\text{C}$), 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 $\varepsilon_{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 $\varepsilon_{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 $\varepsilon_{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 $\varepsilon_{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 $\varepsilon_{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 $\varepsilon_{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 $\varepsilon_{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 $\varepsilon_{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 $\varepsilon_{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 $\varepsilon_{Nd}$, considering dust fluxes,
dissolved riverine sources, sediment sources along the con-
tinental 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, under-
taking 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 comprehen-
sive Earth system model (ESM, specifically the Community
Earth System Model, CESM1.3) to explore the changes to
end-member ε<sub>Nd</sub> 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 ε<sub>Nd</sub> 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, present-
ing 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 mea-
surements 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 previ-
ously 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
ε<sub>Nd</sub> is linked to benthic flux exposure time and the differ-
ence 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 circula-
tion and climate interactions, limiting a clear interpretation
of precisely how (and under what physical or environmen-
tal conditions) the benthic flux may determine global ma-
rine 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 wa-
ters could occur across the entire seafloor. This was done by
removing the depth limitation of the sediment source (pre-
viously 3 km) and invoking a constant benthic flux that es-
capes from all sediment-water interfaces. The scheme was
further extended by revising key source-sink parameterisa-
tions for subsequent investigation of non-conservative Nd
isotope behaviour under different ocean circulation states
(Pöppelmeier et al., 2022). The authors demonstrate sub-
stantial non-conservative effects occur even under strong cir-
culation regimes with low benthic flux exposure times and
are not strictly limited to the deep ocean. This work high-
lights the importance of downward vertical fluxes via re-
versible 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 ε<sub>Nd</sub> distri-
butions.

Thus, despite substantial progress to explicitly describe
seawater Nd budgets, it is apparent that outstanding ques-
tions 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
ε<sub>Nd</sub>, 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 crit-
ical need to progress our understanding of modern marine
Nd cycling in the light of continued use of ε<sub>Nd</sub> as a valu-
able 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 con-
straining 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 com-
plexity 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 run-
ning very quickly to facilitate efficient parameter space ex-
ploration, performance optimisation and long integrations.
Owing to a winning combination of speed and complex-
ity, 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 bio-geochemical 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 ($^{143}$Nd and $^{144}$Nd) 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 = 10⁶ m³ s⁻¹) (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).

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.

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).

## 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 ¹⁴³Nd and ¹⁴⁴Nd, and the isotopic ratio (IR) relates to the ratio of ¹⁴³Nd to ¹⁴⁴Nd, as shown in the following equations:

$$Nd = ¹⁴³Nd + ¹⁴⁴Nd, \qquad (2)$$

$$IR = \frac{¹⁴³Nd}{¹⁴⁴Nd}. \qquad (3)$$

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:

$$¹⁴³Nd = \frac{Nd}{(1 + (1/IR))}, \qquad (4)$$

$$¹⁴⁴Nd = \frac{Nd}{(IR + 1)}. \qquad (5)$$

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⁻³ s⁻¹) are multiplied by a factor of 10¹⁸. 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 10¹⁸ kg m⁻³. A full list of all variables described in the text and their abbreviations is given in Table 1.
<table><caption>**Table 1.** Nd scheme model parameters, abbreviations, fixed model parameter values, and units.</caption><thead><tr><th>Variable</th><th>Abbreviation</th><th>Fixed parameter value</th><th>Unit</th></tr></thead><tbody><tr><td>Total Nd concentration</td><td>[Nd]<sub>t</sub></td><td>–</td><td>pmol kg<sup>-1</sup></td></tr><tr><td>Dissolved Nd</td><td>[Nd]<sub>d</sub></td><td>–</td><td>pmol kg<sup>-1</sup></td></tr><tr><td>Particle-associated Nd</td><td>[Nd]<sub>p</sub></td><td>–</td><td>pmol kg<sup>-1</sup></td></tr><tr><td>Nd source, total</td><td>f<sub>total</sub></td><td>–</td><td>g Nd yr<sup>-1</sup></td></tr><tr><td>Nd source, density</td><td>S<sub>total</sub></td><td>–</td><td>g Nd m<sup>-3</sup> yr<sup>-1</sup></td></tr><tr><td>Dust source, total</td><td>f<sub>dust</sub></td><td>0.33</td><td>Gg Nd yr<sup>-1</sup></td></tr><tr><td>Dust source, density</td><td>S<sub>dust</sub></td><td>–</td><td>g Nd m<sup>-3</sup> yr<sup>-1</sup></td></tr><tr><td>Flux of dust</td><td>F<sub>dust</sub></td><td>2D horizontal global field from Hopcroft et al. (2015)</td><td>g m<sup>-2</sup> yr<sup>-1</sup></td></tr><tr><td>Nd concentration dust</td><td>C&#773;<sub>dust</sub></td><td>20</td><td>µg g<sup>-1</sup></td></tr><tr><td>Nd dust dissolution</td><td>β<sub>dust</sub></td><td>0.02</td><td></td></tr><tr><td>Riverine source, total</td><td>f<sub>river</sub></td><td>0.44</td><td>Gg Nd yr<sup>-1</sup></td></tr><tr><td>Riverine source, density</td><td>S<sub>river</sub></td><td>–</td><td>g Nd m<sup>-3</sup> yr<sup>-1</sup></td></tr><tr><td>River discharge</td><td>RIVER</td><td>–</td><td>g m<sup>-2</sup> yr<sup>-1</sup></td></tr><tr><td>Riverine scaling factor</td><td>α<sub>river</sub></td><td>1</td><td></td></tr><tr><td>Nd concentration river</td><td>C<sub>river</sub></td><td>–</td><td>µg g<sup>-1</sup></td></tr><tr><td>Nd removal, estuaries</td><td>y&#773;<sub>river</sub></td><td>0.7</td><td></td></tr><tr><td>Sediment source, total</td><td>f<sub>sed</sub></td><td>–</td><td>g Nd yr<sup>-1</sup></td></tr><tr><td>Sediment source, density</td><td>S<sub>sed</sub></td><td>–</td><td>g Nd m<sup>-3</sup> yr<sup>-1</sup></td></tr><tr><td>Total sediment surface</td><td>A<sub>total</sub></td><td>–</td><td>m<sup>2</sup></td></tr><tr><td>Grid box sediment surface</td><td>A(i, k)</td><td>–</td><td>m<sup>2</sup></td></tr><tr><td>Grid box volume</td><td>V(i, k)</td><td>–</td><td>m<sup>3</sup></td></tr><tr><td>Thickness of euphotic layer</td><td>z<sub>eu</sub></td><td>81</td><td>m</td></tr><tr><td>Penetration depth of opal</td><td>l<sub>opal</sub></td><td>10 000</td><td>m</td></tr><tr><td>Penetration depth of CaCO<sub>3</sub></td><td>l<sub>calcite</sub></td><td>3500</td><td>m</td></tr><tr><td>Particle settling velocity</td><td>ω</td><td>1000</td><td>m yr<sup>-1</sup></td></tr><tr><td>Ratio [Nd]<sub>p</sub> to [Nd]<sub>d</sub></td><td>[Nd]<sub>p</sub>/[Nd]<sub>d</sub></td><td>–</td><td></td></tr><tr><td>Global average density of seawater</td><td>ρ</td><td>1024.5</td><td>kg m<sup>-3</sup></td></tr><tr><td>Reversible scavenging, density</td><td>S<sub>rs</sub></td><td>–</td><td>g Nd m<sup>-3</sup> yr<sup>-1</sup></td></tr><tr><td>Ratio between average POC concentration and density of seawater</td><td>R<sub>POC</sub></td><td>2.93 × 10<sup>-9</sup></td><td></td></tr><tr><td>Ratio between average CaCO<sub>3</sub> concentration and density of seawater</td><td>R<sub>CaCO<sub>3</sub></sub></td><td>6.27 × 10<sup>-9</sup></td><td></td></tr><tr><td>Ratio between average opal concentration and density of seawater</td><td>R<sub>opal</sub></td><td>5.21 × 10<sup>-9</sup></td><td></td></tr><tr><td>Ratio between average dust concentration and density of seawater</td><td>R<sub>dust</sub></td><td>1.73 × 10<sup>-9</sup></td><td></td></tr><tr><td>Total Nd inventory after equilibrium</td><td>Nd(I)</td><td>–</td><td>Tg</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Taylor 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.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

### 2.3.1 Dust source

Surface dust deposition ($F_{dust}$) 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-
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Atlantic 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.</figcaption></figure>

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 ($S_{dust}$: g m⁻³ yr⁻¹) in the uppermost layer of the ocean model, assuming a global mean concentration of Nd in dust $\bar{C}_{dust}$ ($\bar{C}_{dust} = 20 \mu g g^{-1}$) (Goldstein et al., 1984; Grousset et al., 1988, 1998), from which only a certain fraction $\beta_{dust}$ ($\beta_{dust} = 0.02$: Greaves et al., 1994) dissolves in seawater.

$$S_{dust}(i, k) = \frac{F_{dust}(i, 1) \times \bar{C}_{dust} \times \beta_{dust}}{dz(i, 1)} \quad (6)$$

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 ($f_{dust}$) is 0.33 Gg yr⁻¹. Although we use an updated dust deposition field compared to previous studies, prescribed $\bar{C}_{dust}$ and $\beta_{dust}$ are broadly consistent with earlier Nd isotope schemes, providing a comparable total Nd dust source (0.1–0.5 Gg yr⁻¹; Arsouze et al., 2009; Gu et al., 2019; Pöppelmeier et al., 2020b; Rempfer et al., 2011; Tachikawa et al., 2003).

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: $\varepsilon_{Nd} = -15$; Atlantic < 50°N: $\varepsilon_{Nd} = -12$; North Pacific > 44°N: $\varepsilon_{Nd} = -5$; Indo-Pacific < 44°N: $\varepsilon_{Nd} = -7$; and remainder:

<figure><img src="image_6.png" /><figcaption>**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 $\varepsilon_{Nd}$ signal from dust deposition following Tachikawa et al. (2003) and updated with information from Mahowald et al. (2006) and Blanchet (2019).</figcaption></figure>
$ε_{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⁻² yr⁻¹) and combined this outflow with both the Nd concentration ($C_{river}$; µg g⁻¹) 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.

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 $C_{river}$ (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 $C_{river}$ 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

<figure><img src="image_7.png" /><figcaption>**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).</figcaption></figure>

($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. $\bar{\gamma}_{river} = 0.7$).

The dissolved riverine source per unit volume of Nd ($S_{river}$: g m⁻³ yr⁻¹) in the uppermost layer of the ocean is therefore
calculated as follows:

$$S_{\text{river}}(i, k) = \frac{\text{RIVER}(i, 1) \times C_{\text{river}} \times (1.0 - \bar{y}_{\text{river}})}{\text{dz}(i, 1)} \quad (7)$$

The total global flux of river-sourced dissolved Nd to seawater in the model ($f_{river}$) 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 $\epsilon_{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 horizontal 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 ($S_{sed}$: 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 ($f_{sed}$: g Nd yr$^{-1}$) is used as a tuning parameter.

$$S_{\text{sed}}(i, k) = f_{\text{sed}} \times \frac{A(i, k)}{A_{\text{total}}} \times \frac{1}{V(i, k)}, \quad (8)$$

where $A(i, k)$ is the total area of the sediment surface in contact with seawater per ocean grid cell (m$^2$), $A_{total}$ is the total global area of the sediment surface where a sediment source occurs (m$^2$), and $V(i, k)$ is the volume of water per ocean grid cell (m$^3$).

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 ($f_{sed}$), 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
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Map of the global ε<sub>Nd</sub> 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).</figcaption></figure>

is unknown and challenging to resolve. Therefore, in accor-
dance 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, in-
cluding the possibility in future work of incrementally mod-
ifying the sediment source distributions to increase the com-
plexity 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 ε<sub>Nd</sub> 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).

Pasquier et al. (2022) first applied the sedimentary ε<sub>Nd</sub>
map from Robinson et al. (2021) in their recent global marine
Nd isotope scheme, imposing positive (i.e. radiogenic) mod-
ifications to the Pacific sedimentary ε<sub>Nd</sub> and using a reac-
tivity scaling factor (linked to sediment lability) that favours

more extreme ε<sub>Nd</sub> signals. Here, we again adopt a simpler ap-
proach, imposing the unmodified sediment ε<sub>Nd</sub> 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, CaCO₃; 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, CaCO₃, and opal) are prescribed using gridded, global satellite-derived particle export productivity from Dunne et al. (2012, 2007). The euphotic zone ($z_{eu}$) 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 $z_{eu}$, 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 ($F_{POC}$) follow the power law profile of Martin et al. (1987):

$$F_{POC}(z) = F_{POC}(z_{eu}) \times \left( \frac{z}{z_{eu}} \right)^{-\alpha} \quad (\text{for } z > z_{eu}), \quad (9)$$

where $z$ is depth (m) and $\alpha$ 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 ($F_{opal}$) and CaCO₃ ($F_{CaCO₃}$) follow exponential dissolution profiles with particle-specific length scales, i.e. $l_{opal}$ (10 000 m) and $l_{calcite}$ (3500 m) (Maier-Reimer, 1993; Henderson et al., 1999; Najjar and Orr, 1998), respectively, calculated as follows:

$$F_{opal}(z) = F_{opal}(z_{eu}) \times \exp\left(\frac{z_{eu} - z}{l_{opal}}\right), \quad (10)$$

$$F_{CaCO_3}(z) = F_{CaCO_3}(z_{eu}) \times \exp\left(\frac{z_{eu} - z}{l_{CaCO_3}}\right). \quad (11)$$

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öpkelmeier et al., 2020b; Rempfer et al., 2011; Siddall et al., 2008), a uniform settling velocity ($\omega$) of 1000 m yr⁻¹ 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⁻¹), CaCO₃ (0.45 Gt(C) yr⁻¹), and opal (90 Tmol(Si) yr⁻¹) are comparable with previous estimates, although export of CaCO₃ and opal are at the lower end, as highlighted in Table 3 by Dunne et al. (2007). Note that the annual export of CaCO₃ reflects the new optimised surface calcite parameterisation as described in Dunne et al. (2012).

Reversible scavenging ($Nd_{rs}$) considers total Nd for each isotope (j) as the dissolved Nd ($Nd_d$) and particulate Nd ($Nd_p$) associated with the different particle fields ($\chi$, where $\chi$ refers to POC, CaCO₃, opal, and dust).

$$[Nd]_t^j = [Nd]_d^j + [Nd]_p^j = [Nd]_d^j + \sum_{\chi} [Nd]_{p,\chi}^j \quad (12)$$

$Nd_p$ 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 $Nd_d$ in shallower depths and as sources at greater depths. This combined process for reversible scavenging ($S_{rs}$) in the model can be described by

$$S_{rs}(i, k) = \frac{-\partial(\omega[\mathrm{Nd}]_p(i, k))}{\partial z(i, k)}, \quad (13)$$

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

$$[Nd]_p^j(i, k) = [Nd]_t^j(i, k) \times \left( 1 - \frac{1}{1 + \sum_{\chi} R_{\chi}(i, k) \times K_{\chi}^j} \right). \quad (14)$$

$R_{\chi}(i, k)$ describes the dimensionless ratio between particle concentration for each particle type ($C_{\chi}(i, k)$) and the average density of seawater ($\rho$: 1024.5 kg m⁻³), input as fixed boundary conditions in our scheme. $C_{\chi}$ is calculated from the prescribed particle fluxes ($F_{\chi}$: Fig. 7) by assuming a globally uniform settling velocity ($\omega$ = 1000 m yr⁻¹; Arsouze et al., 2009; Dutay et al., 2009; Gu et al., 2019; Kriest and Oschlies, 2008; Rempfer et al., 2011), i.e. $C_{\chi} = F_{\chi}/\omega$. The equilibrium between dissolved concentration ($[Nd]_d$) and concentration associated with particle type ($[Nd]_p$) can be described by the equilibrium partition coefficient ($K_{\chi}^j$):

$$K_{\chi}^j = \frac{[Nd]_{p,\chi}^j}{[Nd]_d^j} \times \frac{1}{R_{\chi}}, \quad (15)$$

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_{\chi}$, however, is dependent on particle type (where $R_{\chi} = C_{\chi}/\rho$), and thus $K_{\chi}$ is different for different particles. Our global mean particle concentrations ($\bar{C}_{\chi}$; Table 2) are similar to those reported in previous schemes.

Isotopic fractionation during absorption or incorporation and desorption or dissolution is neglected due to similar
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Particle export fields from the ocean surface (g m⁻² yr⁻¹). Biogenic particle fields are prescribed using satellite-derived export productivity fields for (a) POC, (b) CaCO₃, 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).</figcaption></figure>

<table><caption>**Table 2.** Global mean particle concentrations for each particle type used to calculate equilibrium scavenging coefficients following Eq. (15). In summary, export fluxes of POC, CaCO₃, 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.</caption><thead><tr><th rowspan="2">Particle type</th><th rowspan="2">Acronym</th><th colspan="3">Global mean particle concentration (kg m<sup>-3</sup>)</th></tr><tr><th>This study</th><th>Rempfer et al. (2011)</th><th>Gu et al. (2019)</th></tr></thead><tbody><tr><td>POC</td><td>C&#773;<sub>POC</sub></td><td>3.00 &times; 10<sup>-6</sup></td><td>3.30 &times; 10<sup>-6</sup></td><td>2.66 &times; 10<sup>-6</sup></td></tr><tr><td>CaCO<sub>3</sub></td><td>C&#773;<sub>CaCO<sub>3</sub></sub></td><td>6.43 &times; 10<sup>-6</sup></td><td>1.60 &times; 10<sup>-6</sup></td><td>9.53 &times; 10<sup>-6</sup></td></tr><tr><td>Opal</td><td>C&#773;<sub>opal</sub></td><td>5.33 &times; 10<sup>-6</sup></td><td>5.90 &times; 10<sup>-6</sup></td><td>8.30 &times; 10<sup>-6</sup></td></tr><tr><td>Dust</td><td>C&#773;<sub>dust</sub></td><td>1.78 &times; 10<sup>-6</sup></td><td>1.30 &times; 10<sup>-6</sup></td><td>1.23 &times; 10<sup>-6</sup></td></tr></tbody></table>

masses of ¹⁴³Nd and ¹⁴⁴Nd, 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:

$$ \frac{\delta[\mathrm{Nd}]_t^j}{\delta t} = S_{\mathrm{dust}}^j + S_{\mathrm{river}}^j + S_{\mathrm{sed}}^j + S_{\mathrm{rs}}^j + \mathrm{transport}. \quad (16) $$

## 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):

$$ \mathrm{MAE} = \frac{1}{N} \sum_{k=1}^{N} |\mathrm{obs}_k - \mathrm{sim}_k|, \quad (17) $$

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⁻¹ 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:

$$RMSE = \sqrt{\frac{1}{N} \sum_{k=1}^{N} \frac{(\text{obs}_k - \text{sim}_k)^2}{N}}. \quad (18)$$

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⁻¹ 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.

Neither the horizontal nor vertical distribution of global seawater [Nd]d and εNd observational data is even. Most [Nd]d measurements 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

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Location 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.</figcaption></figure>

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öpkelmeier 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 re-
quire further constraining (through further field campaigns, laboratory analyses, or model experimentation).

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; Ar-souze 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 (Ar-souze et al., 2009; Rempfer et al., 2011; Gu et al., 2019; Pöp-pelmeier et al., 2020a), *f*<sub>sed</sub> is fixed at 4.5 Gg yr<sup>-1</sup> through-out.

Second, *f*<sub>sed</sub> is varied in four sensitivity simulations (*f*<sub>sed</sub> ranging 1.5–6.0 Gg yr<sup>-1</sup>), using values based upon previous recent estimates of a global sediment flux to seawater (3.3–5.5 Gg yr<sup>-1</sup>; 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 *f*<sub>sed</sub> sensitivity studies alter the percentage contribution of the sediment Nd flux to the total Nd flux to seawater from 66% (where *f*<sub>sed</sub> = 1.5 Gg yr<sup>-1</sup>) to 89% (where *f*<sub>sed</sub> = 6.0 Gg yr<sup>-1</sup>). 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.,

<table><caption>**Table 3.** Suite 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 (*f*<sub>sed</sub>). 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: https://cms.ncas.ac.uk/puma/ 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.</caption><thead><tr><th>Simulation name</th><th>Simulation identifier</th><th>f<sub>sed</sub><br>(Gg yr<sup>-1</sup>)</th><th>[Nd]<sub>p</sub>/[Nd]<sub>d</sub></th><th>Total Nd flux<br>(Gg yr<sup>-1</sup>)</th><th>Nd inventory<br>(Tg)</th><th>Residence time<br>(years)</th><th>MAE<sub>[Nd]</sub></th><th>RMSE<sub>[Nd]</sub></th><th>% within<br>10 pmol kg<sup>-1</sup></th><th>MAE<sub>εNd</sub></th><th>RMSE<sub>εNd</sub></th><th>% within<br>3 εNd-units</th></tr></thead><tbody><tr><td colspan="13" style="text-align: left; font-weight: bold;">Varying [Nd]<sub>p</sub>/[Nd]<sub>d</sub></td></tr><tr><td>EXPT_RS1</td><td>XPDAI</td><td>4.5</td><td>0.001</td><td>5.27</td><td>16</td><td>3036</td><td>44.46</td><td>48.7</td><td>6</td><td>3.34</td><td>4.12</td><td>48</td></tr><tr><td>EXPT_RS2</td><td>XPDAI</td><td>4.5</td><td>0.002</td><td>5.27</td><td>8.91</td><td>1691</td><td>15.96</td><td>19.78</td><td>41</td><td>3.11</td><td>3.89</td><td>53</td></tr><tr><td>EXPT_RS3</td><td>XPDAH</td><td>4.5</td><td>0.003</td><td>5.27</td><td>6.03</td><td>1145</td><td>10.13</td><td>13.42</td><td>58</td><td>2.88</td><td>3.67</td><td>57</td></tr><tr><td>EXPT_RS4</td><td>XPDAE</td><td>4.5</td><td>0.004</td><td>5.27</td><td>4.51</td><td>856</td><td>9.66</td><td>13.01</td><td>62</td><td>1.67</td><td>3.45</td><td>60</td></tr><tr><td>EXPT_RS5</td><td>XPDAF</td><td>4.5</td><td>0.005</td><td>5.27</td><td>3.58</td><td>679</td><td>10.52</td><td>13.77</td><td>56</td><td>2.46</td><td>3.19</td><td>64</td></tr><tr><td>EXPT_RS6</td><td>XPDAG</td><td>4.5</td><td>0.006</td><td>5.27</td><td>2.95</td><td>559</td><td>11.7</td><td>14.63</td><td>50</td><td>2.27</td><td>3.00</td><td>70</td></tr><tr><td colspan="13" style="text-align: left; font-weight: bold;">Varying f<sub>sed</sub></td></tr><tr><td>EXPT_SED1</td><td>XPDAL</td><td>1.5</td><td>0.003</td><td>2.27</td><td>2.5</td><td>1110</td><td>10.71</td><td>13.82</td><td>53</td><td>2.71</td><td>3.55</td><td>60</td></tr><tr><td>EXPT_SED2</td><td>XPDAM</td><td>3.0</td><td>0.003</td><td>3.77</td><td>3.89</td><td>1032</td><td>7.96</td><td>11.40</td><td>73</td><td>2.82</td><td>3.68</td><td>62</td></tr><tr><td>EXPT_SED3</td><td>XPDAH</td><td>4.5</td><td>0.003</td><td>5.72</td><td>6.03</td><td>1145</td><td>10.13</td><td>13.42</td><td>58</td><td>2.88</td><td>3.67</td><td>57</td></tr><tr><td>EXPT_SED4</td><td>XPDAN</td><td>6.0</td><td>0.003</td><td>6.77</td><td>7.79</td><td>1150</td><td>13.23</td><td>17.19</td><td>47</td><td>2.93</td><td>3.71</td><td>55</td></tr></tbody></table>
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).

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<sub>εNd</sub> 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<sub>εNd</sub>. 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<sub>εNd</sub> 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<sub>εNd</sub> 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).

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 parti-
cles fluxes – the central driver of ocean biogeochemical cy-
cling – remain fundamentally difficult to obtain (Dunne et al.,
2007), and our particle fluxes may be inaccurate as a conse-
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Global 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 \times 10^{12}$ g by Tachikawa et al. (2003) used as a first-order target for our simulations.</figcaption></figure>

quence (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 \text{ pmol kg}^{-1}$ compared to measured $[Nd]_d$ of $70 \text{ pmol kg}^{-1}$) and the Sea of Japan (where simulated $[Nd]_d$ is $2 \text{ pmol kg}^{-1}$ compared to measured $[Nd]_d$ of $50 \text{ 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
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Maps 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.</figcaption></figure>

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 ex-
planations for why our simulated Pacific εNd is so much
lower than what has been recorded in modern seawater in-
clude that the model boundary conditions, specifically the
marine sediment source of Nd (taken directly from Robin-
son 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 param-
eterisations.

For εNd, the response to varying scavenging efficiency has
varied effects across depths and between ocean regions, indi-
cating a more complex relationship between how reversible
scavenging can delineate global εNd distributions in compar-
ison 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 combina-
tion of strong scavenging efficiency and a short seawater res-
idence 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 re-
versible 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.

Typical depth profiles from measurements in the subtropi-
cal 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 con-
sistency 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 seawa-
ter εNd signal. An additional explanation may be due to the
lack of abyssal (i.e. below 3 km) AABW, which does not ex-
tend past ~20° N (Fig. 3). This insufficient AABW produc-
tion 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 previ-
ous 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 Pa-
cific 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⁻¹ and is
therefore masked out of the results shown by Fig. 12 to avoid
misinterpretation.

The simulated depth profile in all experiments in the In-
dian 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 uni-
form ε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⁻¹; observed [Nd]d = 8 pmol kg⁻¹), indicating that
the model is not fully resolving either the concentration or
<figure><img src="image_13.png" /></figure>

**Figure 11.** The 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.

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 pro-
files 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 sed-
iment 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 sur-
face to −8.2 at 1000 m). Below this (and down to 3000 m),
observed diversions in εNd are not captured in any of the sen-
sitivity experiments. This discrepancy can be attributed to the
simulation of quite homogenous AABW throughout the wa-
ter 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 Cir-
cumpolar 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 sedi-
mentation is necessary to balance the simulated input sources
and enables the scheme to reach equilibrium around reason-
able Nd inventories. We find that reversible scavenging is
an important process that enhances the εNd gradient between
oceans by maintaining localised basinal εNd signals through-
out the water column. The strength of this process is particu-
larly 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 stag-
nant modern Pacific Ocean.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** The central panel (g) displays $\varepsilon_{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 $\ge$ 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) $\varepsilon_{Nd}$, respectively. Larger shifts in the $\varepsilon_{Nd}$ between simulations highlight the regions most sensitive to the efficiency of reversible scavenging.</figcaption></figure>

Nonetheless, parameterising reversible scavenging efficiency alone cannot account for the correct trends and magnitude within $\varepsilon_{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 $\varepsilon_{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 $\varepsilon_{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 $\varepsilon_{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 ($f_{sed}$)

The second set of sensitivity simulations tests the response of [Nd]$_d$ and $\varepsilon_{Nd}$ to systematically varying the total Nd flux from the sediment ($f_{sed}$). 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 $f_{sed}$ sensitivity ex-
periments have reached steady state (< 0.02 % change per 100 years).

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 cm² yr⁻¹, 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 cm² yr⁻¹, 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⁻¹ (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 cm² yr⁻¹ 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 cm² yr⁻¹ (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 cm² yr⁻¹ but an order of magnitude lower than the value of 20–30 pmol cm² yr⁻¹ 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 cm² yr⁻¹ and deep-ocean fluxes of 1.11 pmol cm² yr⁻¹.

<table><caption>**Table 4.** Nd inventory (Tg) for EXPT_SED2 within different ocean regions.</caption><thead><tr><td>Ocean region</td><td>Nd inventory (Tg)</td></tr></thead><tbody><tr><td>Global</td><td>3.89</td></tr><tr><td>Arctic Ocean</td><td>0.05</td></tr><tr><td>North Atlantic</td><td>0.33</td></tr><tr><td>South Atlantic</td><td>0.45</td></tr><tr><td>North Pacific</td><td>0.96</td></tr><tr><td>South Pacific</td><td>0.76</td></tr><tr><td>Indian Ocean</td><td>0.63</td></tr><tr><td>Southern Ocean</td><td>0.28</td></tr></tbody></table>

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.

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⁻¹ compared to 50 pmol kg⁻¹ 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⁻¹) 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.

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)
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Global Nd inventory (Tg) simulated with different values for the total sediment flux tuning parameter ($f_{sed}$) 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.</figcaption></figure>

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 $f_{sed}$) 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 $f_{sed}$. We find that [Nd]d is sensitive to $f_{sed}$ 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).

Typically, the best model-data fit for [Nd]d depth profiles is achieved under lower $f_{sed}$ (1.5 to 3.0 Gg yr$^{-1}$) and in the Pacific, Indian, and Southern oceans in particular (Fig. 15). Correspondingly, under high $f_{sed}$ 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 $f_{sed}$ 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 $f_{sed}$ 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 $f_{sed}$ in the upper 1000 m, with EXPT_SED4's simulated surface concentrations of 10 pmol kg$^{-1}$ ($f_{sed}$ = 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 $f_{sed}$ 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
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Maps 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.</figcaption></figure>

$f_{sed}$ 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 $f_{sed}$ 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 $f_{sed}$ 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.

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 $f_{sed}$ are conditioned towards riverine εNd, providing the typical non-radiogenic εNd signature of NADW (Fig. 16a–d). Conversely,

enhancing $f_{sed}$ 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 $f_{sed}$ 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
<figure><img src="image_17.png" /><figcaption>**Figure 15.** The 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 $f_{sed}$) 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.</figcaption></figure>

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 $\varepsilon_{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 $\varepsilon_{Nd}$ tagging of Atlantic seawater.

In contrast to the Atlantic, the Pacific Ocean $\varepsilon_{Nd}$ is the least sensitive to varying $f_{sed}$ (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 $\varepsilon_{Nd}$ signal from the bulk seafloor sediment cannot explain the radiogenic $\varepsilon_{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-
<figure><img src="image_18.png" /></figure>

**Figure 16.** The central panel (g) displays $\varepsilon_{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 $\ge$ 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 $f_{sed}$) and measured (filled circles) $\varepsilon_{Nd}$, respectively. Larger shifts in the $\varepsilon_{Nd}$ between simulations highlight the regions most sensitive to the magnitude of the seafloor sediment source.

<figure><img src="image_19.png" /></figure>

**Figure 17.** Volume-weighted distributions of $\varepsilon_{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).
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 f_sed 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 ($^{143}$Nd and $^{144}$Nd) into the ocean component of the FA-MOUS 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 disproportionately 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 FA-MOUS. In the first instance, we suggest calibration of the key tuning parameters ([Nd]p/[Nd]d and f_sed) 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: https://doi.org/10.5518/1136) 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; https://cms.ncas.ac.uk/; last access: 12 February 2013), with specific FAMOUS documentation available at https://cms.ncas.ac.uk/miscellaneous/um-famous/ (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 (https://doi.org/10.5518/1136); 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 ($f_{sed}$) sensitivity simulations are as follows:

- XPDAL $f_{sed} = 1.5 \text{ Gg yr}^{-1}$ (0–9000 years),

- XPDAM $f_{sed} = 3.0 \text{ Gg yr}^{-1}$ (0–9000 years),

- XPDAH $f_{sed} = 4.5 \text{ Gg yr}^{-1}$ (0–9000 years),

- XPDAN $f_{sed} = 6.0 \text{ Gg yr}^{-1}$ (0–9000 years).

**Data availability.** The data are available via the Research Data Leeds Repository (Robinson et al., 2022: https://doi.org/10.5518/1136.)

**Supplement.** 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: https://doi.org/10.5194/gmd-16-1231-2023-supplement.

**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 $\epsilon_{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.

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

**Acknowledgements.** 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.
## References

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, https://doi.org/10.1130/g45904.1, 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, https://doi.org/10.1130/G37114.1, 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, https://doi.org/10.1016/j.gca.2015.01.010, 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, https://doi.org/10.3389/fmars.2019.00504, 2019.

Anderson, S. P.: Glaciers show direct linkage between erosion rate and chemical weathering fluxes, Geomorphology, 67, 147–157, https://doi.org/10.1016/j geomorph.2004.07.010, 2005.

Arsouze, T., Dutay, J.-C., Lacan, F., and Jeandel, C.: Modeling the neodymium isotopic composition with a global ocean circulation model, Chem. Geol., 239, 165–177, https://doi.org/10.1016/j.chemgeo.2006.12.006, 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, https://doi.org/10.5194/bg-6-2829-2009, 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, https://doi.org/10.5194/bg-13-5259-2016, 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, https://doi.org/10.5194/bg-20-205-2023, 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, https://doi.org/10.1029/JC087iC03p02045, 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, https://doi.org/10.1016/j.epsl.2004.05.033, 2004.

Behrenfeld, M. J. and Falkowski, P. G.: Photosynthetic rates derived from satellite-based chlorophyll concentration, Limnol. Oceanogr., 1–20, https://doi.org/10.4319/lo.1997.42.1.0001, 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, https://doi.org/10.1016/0016-7037(93)90087-D, 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, https://doi.org/10.5194/essd-11-741-2019, 2019.

Blaser, P., Frank, M., and van de Flierdt, T.: Revealing past ocean circulation with neodymium isotopes, Past Glob. Chang. Mag., 27, 54–55, https://doi.org/10.22498/pages.27.2.54, 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, https://doi.org/10.1016/j.gca.2018.10.018, 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, https://doi.org/10.1016/j.epsl.2020.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, https://doi.org/10.1029/2019RG000655, 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, https://doi.org/10.1016/j.chemgeo.2018.09.021, 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, https://doi.org/10.1016/j.gca.2013.09.011, 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, https://doi.org/10.5194/gmd-4-1051-2011, 2011.

Cox, P. M., Betts, R. A., Bunton, C. B., Essery, R. L. H., Rown-tree, 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, https://doi.org/10.1007/s003820050276, 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, https://doi.org/10.7185/geochemlet.2223, 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, https://doi.org/10.1007/s00382-018-4243-y, 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, https://doi.org/10.5194/gmd-13-3529-2020, 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, https://doi.org/10.1029/2003GB002150, 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, https://doi.org/10.1016/j.gca.2016.08.005, 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 CO₂ concentrations, Nat. Geosci., 11, 749–755, https://doi.org/10.1038/s41561-018-0205-6, 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, https://doi.org/10.1016/j.quascirev.2020.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, https://doi.org/10.1016/j.epsl.2022.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, https://doi.org/10.1029/2006GB002907, 2007.

Dunne, J. P., Hales, B., and Toggweiler, J. R.: Global calcite cycling constrained by sediment preservation controls, Global Biogeochem. Cy., 26, GB3023, https://doi.org/10.1029/2010GB003935, 2012.

Dutay, J. C., Lacan, F., Roy-Barman, M., and Bopp, L.: Influence of particle size and type on ²³¹Pa and ²³⁰Th simulation with a global coupled biogeochemical-ocean general circulation model: A first approach, Geochem. Geophy. Geosy., 10, GB3023, https://doi.org/10.1029/2008GC002291, 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, https://doi.org/10.1016/0016-7037(90)90432-K, 1990.

Essery, R. and Clark, D. B.: Developments in the MOSES 2 land-surface model for PILPS 2e, Glob. Planet. Change, 38, 161–164, https://doi.org/10.1016/S0921-8181(03)00026-2, 2003.

Essery, R. L. H., Best, M. J., and Cox, P. M.: MOSES 2.2 technical documentation, Hadley Centre Technical Note 30, https://jules.jchmr.org/sites/default/files/HCTN_30.pdf (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, https://doi.org/10.1146/annurev-earth-082517-010045, 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, https://doi.org/10.3389/fmars.2019.00260, 2019.

Frank, M.: Radiogenic isotopes: Tracers of past ocean circulation and erosional input, Rev. Geophys., 40, 1001, https://doi.org/10.1029/2000RG000094, 2002.

GEOTRACES Intermediate Data Product Group: The GEO-TRACES Intermediate Data Product 2021 (IDP2021), NERC EDS British Oceanographic Data Centre NOC [data set], https://doi.org/10.5285/cf2d9ba9-d51d-3b7c-e053-8486abc0f5fd, 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 systems of river-water dissolved material: Implications for the sources of Nd and Sr in seawater, Chemical Geology: Isotope Geoscience Section, 66, 245–272, https://doi.org/10.1016/0168-9622(87)90045-5, 1987.

Goldstein, S. L. and Hemming, S. R.: Long-lived Isotopic Tracers in Oceanography, Paleooceanography, and Ice-Sheet Dynamics, Treatise of Geochemistry, 6, 453–483, https://doi.org/10.1016/B0-08-043751-6/06179-X, 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, https://doi.org/10.1016/0012-821X(84)90007-4, 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, https://doi.org/10.1007/s003820050010, 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, https://doi.org/10.1016/0304-4203(94)90081-7, 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, https://doi.org/10.1007/s00382-010-0934-8, 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, https://doi.org/10.1038/nature11257, 2012.
Gregoire, L. J., Valdes, P. J., and Payne, A. J.: The relative con-
tribution of orbital forcing and greenhouse gases to the North
American deglaciation, Geophys. Res. Lett., 42, 9970–9979,
https://doi.org/10.1002/2015GL066005, 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, https://doi.org/10.1002/2016GL070356, 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, https://doi.org/10.1029/2012JC008239, 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, https://doi.org/10.1016/0012-821X(88)90001-5, 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, https://doi.org/10.1016/S0277-3791(97)00048-6, 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, https://doi.org/10.1029/2018MS001538, 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, https://doi.org/10.1016/j.gca.2003.09.012, 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, https://doi.org/10.1016/j.gca.2013.11.012, 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, https://doi.org/10.3389/fmars.2017.00426, 2017.

Henderson, G. M., Heinze, C., Anderson, R. F., and Winguth, A. M. E.: Global distribution of the ²³⁰Th flux to ocean sediments constrained by GCM modelling, Deep-Res. Pt. I, 46, 1861–1893, https://doi.org/10.1016/S0967-0637(99)00030-8, 1999.

Hopcroft, P. O. and Valdes, P. J.: Last glacial maximum constraints on the Earth System model HadGEM2-ES, Clim. Dynam., 45, 1657–1672, https://doi.org/10.1007/s00382-014-2421-0, 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, https://doi.org/10.1002/2015JD023742, 2015.

Jacobsen, S. B. and Wasserburg, G. J.: Sm-Nd evolution of chondrites, Earth Planet. Sc. Lett., 50, 139–155, https://doi.org/10.1016/0012-821X(80)90125-9, 1980.

Jeandel, C.: Concentration and isotopic composition of Nd in the South Atlantic Ocean, Earth Planet. Sc. Lett., 117, 581–591, https://doi.org/10.1016/0012-821x(93)90104-h, 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, https://doi.org/10.1016/0016-7037(94)00367-U, 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, https://doi.org/10.1016/j.epsl.2006.10.021, 2007.

Jones, C.: A fast ocean GCM without flux adjustments, J. Atmos. Ocean. Technol., 20, 1857–1868, https://doi.org/10.1175/1520-0426(2003)020<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, https://doi.org/10.1007/s00382-005-0027-2, 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, https://doi.org/10.1016/j.epsl.2008.05.027, 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, https://doi.org/10.1038/ngeo1185, 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, https://doi.org/10.5194/bg-5-55-2008, 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, https://doi.org/10.1029/2004RG000166, 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, https://doi.org/10.1016/j.epsl.2005.01.004, 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, https://doi.org/10.5194/bg-17-5539-2020, 2020.

Lambelet, M., van de Flierdt, T., Crocket, K., Rehkämper, M., Kreissig, K., Coles, B., Rijkenberg, M. J. A., Geringa, 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 GEO-TRACES GA02 section, Geochim. Cosmochim. Ac., 177, 1–29, https://doi.org/10.1016/j.gca.2015.12.019, 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, https://doi.org/10.1029/2018GL080074, 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, https://doi.org/10.1029/1999GB001229, 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, https://doi.org/10.1029/2005JD006653, 2006.

Maier-Reimer, E.: Geochemical cycles in an ocean general circulation model. Preindustrial tracer distributions, Global Biogeochem. Cy., 7, 645–677, https://doi.org/10.1029/93GB01355, 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, https://doi.org/10.1016/0198-0149(87)90086-0, 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, https://doi.org/10.1016/j.pocean.2014.10.006, 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], https://cms.ncas.ac.uk/miscellaneous/um-famous/, 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, https://doi.org/10.1063/1.2337098, 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, https://doi.org/10.1029/2008GB003353, 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, https://doi.org/10.1007/s10872-021-00600-x, 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, https://doi.org/10.1016/j.marchem.2015.03.013, 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, https://doi.org/10.1002/2016GC006782, 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, https://doi.org/10.1029/2021JC017423, 2021.

Palmer, J. R. and Totterdell, I. J.: Production and export in a global ocean ecosystem model, Deep. Res.-Pt. I, 48, 1169–1198, https://doi.org/10.1016/S0967-0637(00)00080-7, 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, https://doi.org/10.5194/gmd-15-4625-2022, 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, Biogeoosciences, 16, 4829–4849, https://doi.org/10.5194/bg-16-4829-2019, 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, https://doi.org/10.1016/j.epsl.2013.03.023, 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, https://doi.org/10.1029/2020PA003973, 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, https://doi.org/10.1130/G47628.1, 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, https://doi.org/10.5194/bg-18-5447-2021, 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, https://doi.org/10.5194/cp-17-615-2021, 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, https://doi.org/10.1016/JQUASCIREV.2022.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, https://doi.org/10.1016/j.epsl.2019.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, https://doi.org/10.1016/j.gca.2020.11.017, 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, https://doi.org/10.1016/j.gca.2011.07.044, 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, https://doi.org/10.1126/science.1178068, 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 $\epsilon_{Nd}$: An updated compilation for use in understanding marine Nd cycling, Chem. Geol., 567, 120119, https://doi.org/10.1016/j.chemgeo.2021.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], https://doi.org/10.5518/1136, 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, https://doi.org/10.1038/ncomms8592, 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, https://doi.org/10.1126/science.274.5291.1346, 1996.

Sherriff-Tadano, S., Abe-Ouchi, A., and Oka, A.: Impact of midglacial ice sheets on deep ocean circulation and global climate, Clim. Past, 17, 95–110, https://doi.org/10.5194/cp-17-95-2021, 2021.

Sholkovitz, E. R.: The geochemistry of rare earth elements in the Amazon River estuary, Geochim. Cosmochim. Ac., 57, 2181–2190, https://doi.org/10.1016/0016-7037(93)90559-F, 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, https://doi.org/10.1016/j.epsl.2008.07.044, 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, https://doi.org/10.1016/J.POCEAN.2017.12.001, 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, https://doi.org/10.1002/qj.49711649210, 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, https://doi.org/10.5194/gmd-5-269-2012, 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, https://doi.org/10.1029/2009GL038607, 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, https://doi.org/10.5194/gmd-1-53-2008, 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, https://doi.org/10.1021/acsearthspacechem.0c00034, 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, https://doi.org/10.1016/S0967-0637(98)00089-2, 1999.

Tachikawa, K., Athias, V., and Jeandel, C.: Neodymium budget in the modern ocean and paleo-oceanographic implications, J. Geophys. Res., 108, 3254, https://doi.org/10.1029/1999JC000285, 2003.

Tachikawa, K., Arsouze, T., Bayon, G., Bory, A., Colin, C., Duyat, 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, https://doi.org/10.1016/j.chemgeo.2017.03.018, 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, https://doi.org/10.1016/C2009-0-24322-4, 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, https://doi.org/10.5194/gmd-10-3715-2017, 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, https://doi.org/10.1098/rsta.2015.0293, 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, Paleooceanography, 16, 424–434, https://doi.org/10.1029/2000PA000550, 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, https://doi.org/10.1016/j.gca.2021.09.015, 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, https://doi.org/10.5194/gmd-6-141-2013, 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 paleo-oceanography, Geochim. Cosmochim. Ac., 109, 197–221, https://doi.org/10.1016/j.gca.2013.01.042, 2013.
Woodward, S.: Mineral Dust in HadGEM2, Hadley Centre Technical Note 87 Mineral Dust in HadGEM 2 March 2011, https://sds-was.aemet.es/forecast-products/dust-forecasts/ Woodward_2011_HadGEM2.pdf (last access: 15 Febru- ary 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, https://doi.org/10.1016/j.dsr2.2007.12.029, 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., http://www.nodc.noaa.gov/ (last access: 15 February 2023) 2019.