Reconstructing the Nd oceanic cycle using a coupled dynamical – biogeochemical model

Reconstructing the Nd oceanic cycle using a coupled dynamical – biogeochemical model T. Arsouze, J.-C. Dutay, F. Lacan, and C. Jeandel Laboratoire des Sciences du Climat et de l’Environnement (LSCE), IPSL, CEA/UVSQ/CNRS, Orme des Merisiers, Gif-Sur-Yvette, Bat 712, 91191 Gif sur Yvette cedex, France Laboratoire d’Etudes en Géophysique et Océanographie Spatiale (LEGOS), UPS/CNES/CNRS/ IRD, Observatoire Midi-Pyrénées, 14 av. E. Belin, 31400 Toulouse, France now at: Lamont-Doherty Earth Observatory (LDEO), P.O. Box 1000 61 Route 9W, Palisades, NY 10964-1000, USA Received: 11 May 2009 – Accepted: 19 May 2009 – Published: 10 June 2009 Correspondence to: T. Arsouze (arsouze@ldeo.columbia.edu) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
Variations in neodymium isotopic composition (hereafter referred as Nd IC 1 ) observed within the ocean reflect influences from both lithogenic inputs of the element (whose IC varies as a function of age and geological composition of T. Arsouze et al.: Reconstructing the Nd oceanic cycle the continent, Jeandel et al., 2007) and the subsequent redistribution by oceanic circulation.ε Nd data show that these variations are closely linked to water mass distribution at depth, and that far from any continental sources, Nd IC behaves quasi-conservatively (Piepgras and Wasserburg, 1982;Jeandel, 1993;von Blanckenburg, 1999;Goldstein and Hemming, 2003).Hence, the main water masses display characteristic ε Nd values (e.g.NADW: ε Nd ≈ −13.5, AAIW and AABW: ε Nd ≈ −8).This water mass tracer property has been recently explored in the modern ocean, by measuring dissolved Nd IC and concentrations (Jeandel, 1993;Piepgras and Wasserburg, 1980;Piepgras andWasserburg, 1982, 1987;Shimizu et al., 1994;Lacan 2004), but has also been used for reconstructing past ocean circulation from measurements in the authigenic fraction of sediments (Rutberg et al., 2000;Piotrowski et al., 2004Piotrowski et al., , 2005;;Gutjahr et al., 2008).The water mass tracer property of ε Nd is commonly accepted, however our understanding of the complete Nd oceanic cycle is far from sufficient to allow a reliable use of this proxy as a paleocirculation tracer (Arsouze et al., 2008).In particular, the nature and the relative importance of the different Nd sources and sinks, and the dissolved/particulate interactions within the water column remain unconstrained.
Previous studies show a decoupling behaviour between Nd concentration and its IC in the open ocean (Bertram and Elderfield, 1993;Jeandel et al., 1995;Tachikawa et al., 1999Tachikawa et al., , 2003;;Lacan and Jeandel, 2001).This feature has been named the "Nd paradox" (Tachikawa et al., 2003;Goldstein and Hemming, 2003;Lacan and Jeandel, 2005), and results in different properties characterizing the Nd concentration and its IC distribution in the dissolved phase.Vertical profiles of Nd concentration are similar to that of all Rare Earth Elements (excluding Ce), showing low values in surface waters which increase with depth, suggesting the influence of vertical cycling (i.e. the element is scavenged at the surface, sinks with the particles, and is subsequently remineralized at depth).Furthermore, Nd concentrations increase along the thermohaline circulation, which is a typical distribution for nutrients like silicate (Elderfield, 1988), whose characteristic residence time is on the order of ∼2×10 4 years (Broecker and Peng, 1982).In contrast, pronounced ε Nd variations between each oceanic basin indicate that the Nd residence time is shorter than the global oceanic mixing time estimated to ∼10 3 years (Broecker and Peng, 1982).The conservativity of ε Nd within the major oceanic water masses suggests at first sight that ε Nd may vary in the open ocean only by water-mass mixing, excluding vertical cycling.
It has been demonstrated that Nd oceanic budgets that consider only dissolved river and atmospheric dust inputs fail to balance both concentration and Nd IC (Bertram and Elderfield, 1993;Tachikawa et al., 2003;van de Flierdt et al., 2004).Other sources have been suggested in order to reconcile the budget of both quantities, like submarine groundwaters (Johannesson and Burdige, 2007) or input from continental margin inputs (Tachikawa et al., 2003) subsequently leading to the notion of "Boundary Exchange" (BE, strong interactions between continental margin and water masses by co-occurrence of sediment dissolution and Boundary Scavenging, Lacan and Jeandel, 2005).This last process has the advantage of including both a source (Boundary Source) and a sink (Boundary Scavenging) for the element, affecting the IC without changing the observed concentration when both fluxes are similar.
Resolving the "Nd Paradox" hence resides in 1) finding the vertical processes responsible for both Nd concentration increase with depth and the Nd IC conservative property in the water column, and 2) constraining the missing source which may explain the important observed Nd IC gradient observed along the thermohaline circulation together with the relatively moderated increase of Nd concentrations.
Recent modelling of oceanic ε Nd , schematically presented in Fig. 1, has helped to improve our understanding of the Nd oceanic cycle.Arsouze et al. (2007, Fig. 1a), taking into account only BE as a source/sink term, simulated a realistic global ε Nd distribution, suggesting that this process plays a major role in the oceanic cycle of the element.However, this first approach was based on the simulation of the Nd IC only (the Nd concentration being considered constant), using a simple relaxing term parameterization for the BE process.Given that no source flux is explicitly considered in this method, the authors could not address the "Nd paradox" and quantify the processes acting in the oceanic Nd cycle.Jones et al. (2008, Fig. 1b) considered no external sources, but prescribed surface Nd IC with observations, to conclude that ε Nd behaves conservatively in the ocean (changing only by water mass mixing).However they needed to invoke an input to the deep North Pacific, which could represent an input directly to the deep ocean (e.g. by boundary exchange), or vertical cycling.We underline here that prescribing ε Nd at the surface must be considered as an implicit source of Nd, in contradiction to the conclusion of these authors upon the role of mixing.Lastly, Siddall et al. (2008, Fig. 1c) have modelled explicitly both Nd concentration and IC using a reversible scavenging model to test the influence of vertical cycling in the oceanic distribution of this element.These authors suggest that both scavenging and remineralisation processes are important for explaining Nd concentration and IC profiles, consistent with the conclusions of Bertram and Elderfield (1993) and Tachikawa et al. (1999Tachikawa et al. ( , 2003)).If the results of Siddall et al. (2008) provide important implications for solving the "Nd paradox", the authors can not in return estimate the different fluxes involved in the Nd cycle.Indeed, this study by Siddall et al. (2008) is prescribing Nd concentration and ε Nd to observations in surface waters, and leaves the detailed simulation of Nd oceanic sources and sinks for further work.This parameterization, although appropriate to initiate studies using coupled dynamical/biogeochemical models, is insufficient for comprehending the full Nd oceanic cycle and investigating all parts of the "Nd paradox".Also, prescribing Nd IC surface distribution (Jones et al., 2008;Siddall et  Prescribed by data Process studied Relaxing term Arsouze et al. 2007Siddall et al., 2008 This study Jones et al., 2008 Continental margin Fig. 1.Summary of the main published Nd oceanic modelling efforts using an OGCM.In each case, the different processes explicitly simulated are represented.Arrows and legends in blue underline the processes specifically studied, while red components refer to parameters prescribed by the data.(A) Arsouze et al. (2007), modelling only Nd IC, focused on the role of Boundary Exchange using a relaxing term and defined on the first 3000 m of the continental margin: this work highlighted the importance of this process in the element's global cycle.
(B) Jones et al. (2008), prescribing a surface ε Nd value estimated by the data suggested that Nd IC distribution can simply be explained by water mass mixing, except in the North Pacific region, where some external radiogenic inputs are required.(C) Using prescribed surface ε Nd and Nd concentration, as well as particle fields determined by satellite observations and then interpolated within the water column, Siddall et al. (2008) showed the major role of vertical cycling (reversible scavenging process) to reproduce both ε Nd and Nd concentration oceanic distribution.(D) In this study, besides confirming, as in Siddall et al. (2008), the role of vertical cycling in Nd oceanic distribution, we aimed to determine and quantify the different sources involved in the Nd oceanic cycle.The Boundary Exchange process is implicitly taken into account via sediment redissolution along the continental margins (source term) as well as by particle sedimentation (sink term).
al., 2008) and Nd concentration, as well as particle distribution (imposed by satellite data and then extrapolated into the water column, Siddall et al., 2008) limits any potential paleo-application.
This study continues the modelling work initiated by Tachikawa et al. (2003), Arsouze et al. (2007) and Siddall et al. (2008) using a coupled dynamical (Ocean Global Circulation Model, that generates dynamical fields)/biogeochemical (dedicated to carbon cycle and ecosystems studies, that simulates particle fields in the ocean) model.The vertical cycling is simulated using a reversible scavenging model developed for the simulation of trace elements with this coupled model (Dutay et al., 2009).Although these authors identified some limitations of this coupled model for geochemical tracer modelling, as we will see hereafter, our approach is dictated by the need for collaboration with carbon cycling modellers to initiate further studies aimed at developing the biogeochemical model.Such developments may lead to mutual advance for both carbon cycle and geochemical (Pa, Th, δ 13 C, Nd, etc. . . ) modelling (Dutay et al., 2009), and also to a significant improvement on our understanding of the element oceanic cycles in the long term, mostly in the context of the increasing database that will result from the interna-tional GEOTRACES effort (GEOTRACES, 2005).Also, this model is fully prognostic so that it can be easily used for paleo applications or even future climate scenarios.
We explicitly represent and quantify the different sources implied in the oceanic cycle of the element (Fig. 1d).This allows investigation from an independent approach 1) if the Boundary Source could globally represent the primordial flux of the tracer into the ocean, as suggested by Arsouze et al. (2007), and 2) if the Boundary Source and Boundary Scanvenging fluxes of the BE process can be associated with the "missing flux" proposed by Tachikawa et al. (2003) in order to balance the Nd oceanic budget.
Besides being a potential component of primary importance for the Nd oceanic cycle, the determination of trace element sources is essential to better constrain the general lithogenic flux to the ocean.It is important to understand the effects of erosion fluxes on the ocean chemistry, and more particularly on primary production and sedimentation.Studying the BE process may provide substantial information on the general geochemical cycle of the elements, such as carbon and iron, which have a direct impact on climate change.
2 Modelling the Nd oceanic cycle

The dynamical model NEMO-OPA
The dynamical model used is the NEMO-OPA model (IPSL/LOCEAN, Madec, 2006).It includes the sea-ice model LIM (Louvain-La Neuve, Fichefet and Maqueda, 1997), in its low-resolution configuration ORCA2.The definition of the mesh is based on a 2 • ×2 • ×cos (latitude) MER-CATOR grid, with poles defined on the continents so as to get rid of singularities near the North Pole.Meridional resolution increases to 0.5 • near the equator, in order to account for the specific local dynamics.Vertical resolution varies with depth, from 10 m at the surface (12 levels included in the first 125 m) to 500 m at the bottom (31 levels in total).The model is forced at the surface by heat and freshwater fluxes obtained from bulk formulae and ERS satellite data from the tropics, and NCEP/NCAR data from polar regions.Surface salinity is readjusted every 40 days to monthly WA01 data to prevent model drift (Timmermann et al., 2005).The Turbulent Kinetic Energy closure is applied to the mixing layer (Blanke and Delecluse, 1993), and subscale physics is parameterized using the Gent and McWilliams scheme (1990).The same model, in similar configurations, has previously been used for other geochemical tracer simulations (Dutay et al., 2002(Dutay et al., , 2004(Dutay et al., , 2009;;Doney et al., 2004;Arsouze et al., 2007).Despite classical shortcomings of low resolution models (boundary currents too weak, crude representation of sinking of dense water during deep water masses formation), this model satisfyingly simulates the main structures of the global thermohaline circulation.

The biogeochemical model PISCES
The biogeochemical model PISCES developed for carbon cycle modelling was coupled to NEMO.It is a prognostic ecosystem and oceanic carbon cycle model, based on the Hamburg Model of Carbon Cycle version 5 (HAMOCC5, Aumont et al., 2003), that represents the biogeochemical cycles of carbon, oxygen, and five nutrients of primary production (phosphates, nitrates, silicates, ammonium and iron).Redfield ratios are set constant, and are based on phytoplankton growth limited by nutrient availability (Monod, 1942).
The model features two classes of phytoplankton (nanophytoplankton and diatoms), and two classes of zooplankton (microzooplankton and mesozooplankton), as well as three non-living components, including Dissolved Organic Carbon (DOC), small and large particles.Small particles (particles between 2 and 100 µm in size) have a sinking velocity of 3 m/day (determined using the relationship between particle size and the sinking velocity established by Kriest, 2002), and consist of Particulate Organic Carbon (POCs).Large particles include Particulate Organic Carbon with a diameter larger than 100 µm (POCb), biogenic silica (BSi), calcite (CaCO 3 ) and lithogenic particles (atmospheric dust), sinking with a velocity varying from 50 m/day at the surface to 300 m/day at depth.The two classes of POC interact via the processes of aggregation/disaggregation.Also a remineralization process depending on temperature (with a Q 10 of about 1.9) is represented between POCb and POCs, and between POCs and DOC (Aumont and Bopp, 2006).The reference of 100 µm used to differentiate small and large carbonate particles (POCs and POCb, respectively) does not correspond to the definition used by experimentalists, who usually refer to large particles as particles collected in traps or by large volume filtration and has having a diameter larger than 50 µm.This being so, the size range of particles trapped remains unknown.In addition, the particles observed by the Underwater Video Profilers (UVP, Gorsky, 2000) have a diameter larger than 100 µm meaning the same limit as prescribed in the model.
A more detailed description of the model, as well as the equations used, are available as supplementary material of Aumont and Bopp (2006).
An evaluation of the particle fields generated by the PISCES model with available data collected by particle traps, satellite observations and estimations, was performed by Gehlen et al. (2006) and Dutay et al. (2009).The small particle pool represents the main particle stock at the surface (at least one order of magnitude higher than the large particle concentration).The simulated small particle concentrations in surface waters are in agreement with observations, mostly in regions of high productivity (coastal upwellings, Equatorial Pacific, Austral Ocean).In contrast, PISCES generates exaggerated vertical variations of the small particle distributions at depth in the water column leading to small particles concentrations that are too low, due to remineralization and aggregation processes.This POCs vertical gradient largely overestimates the available observations (only a factor of 50 between surface and depth, compared to a factor of 1000 for the model).In contrast, the few CaCO 3 and BSi vertical profiles observed in the water column are successfully replicated by the model (Dutay et al., 2009).

The reversible scavenging model
The increase in Nd concentration with depth suggests a reversible exchange between the dissolved and the particulate phases (Nozaki and Alibo, 2003).Nd is adsorbed onto sinking particles in the surface layers, and is redissolved at depth (we refer to this combined process as vertical cycling).In order to represent this process, we use the approach developed by Nozaki et al. (1981) and Bacon and Anderson (1982), and reformulated by Henderson et al. (1999), Siddall et al. (2005) and Dutay et al. (2009) for 231 Pa and 230 Th modelling, and Siddall et al. (2008) for Nd isotope modelling.This approach considers dissolved and particulate phases as in equilibrium within the ocean, and their relative contribution is set using an equilibrium partition coefficient K.
Considering Nd dissolved concentrations (Nd d ) and Nd particulate concentrations (concentration of Nd by unity of particle mass, Nd p ), the equilibrium partition coefficient K is defined as: where C p is the mass of particles per mass of water.This coefficient K is defined for each type of particles represented in the model: small (POCs) and big (POCb) Particulate Organic Carbon, Biogenic Silica (BSi), calcite (CaCO 3 ) and lithogenic atmospheric dust (litho).
We model the two 143 Nd and 144 Nd isotopes independently and calculate ε Nd and Nd concentration afterward.Observations do not suggest any fractionation between 143 Nd and 144 Nd dissolved, colloidal and particulate phases (Dahlqvist, 2005), as they are two isotopes of the same element, and their masses are quite similar.Partition coefficients (K) are thus assumed as being identical for the two isotopes for each particle type.
In the model, we transport for both 143 Nd and 144 Nd tracers the total concentration (Nd T ), defined as the sum of dissolved concentration (Nd d ), small (POCs: Nd ps ) and big (POCb, BSi, CaCO, litho: Nd pg ) particulate concentration: Applying Eq. ( 1) to the particulate pools to express total concentration as a function of dissolved Nd concentration, we obtain: that leads to: This allows Nd concentrations in small and big particles to be defined a posteriori as a function of partition coefficients and total Nd concentration.The main advantage of this method is that we simulate explicitly the total concentration of the two isotopes (two tracers 143 Nd T and 144 Nd T ) rather than concentration in every phase (dissolved, small particles and all big particles, meaning 12 tracers), which implies a substantial gain of computational cost.
The evolution of the total concentration of the tracer is equal to the sum of all sources (A), influence of vertical cycling (B) and physical transport by advection and diffusion (C).The conservation equation of the tracer can therefore be written: where S(Nd T ) is the Source term of the tracer (cf.Sect.2.4).
The vertical cycling represents the scavenging of Nd by the particles (w s and w b are the sinking velocities of small and big particles, respectively).The simulations are performed "off-line" using pre-calculated dynamical: velocity (U ) and mixing coefficient (k), and particle (POCs, POCb, BSi, and CaCO 3 ) distributions in order to reduce computational costs, which allows us to perform some sensitivity tests.

Description of Nd sources and sink
One of our main objectives is to study the relative influence of the different sources of Nd to the ocean.We therefore explicitly represent the different source of Nd in the ocean in our simulations.
The source of the BE process (Boundary Source) is assumed to be the dissolution of a small percentage of the sediments deposited along the continental margin.We specify this source in the model by imposing an input flux on each continental margin grid point of the model between 0 and 3000 m (Boillot and Coulon, 1998).This Nd Boundary Source is written as: where F sed is the source flux of sedimentary Nd to the ocean (in g(Nd)/m 2 /yr).Sediment flux is assumed here as geographically constant.This hypothesis may be not verified in the ocean and disparities could play an important role on a global scale.However, considering the lack of knowledge of the processes acting on this source, it seems reasonable to set it as constant as a first approximation.F sed is then determined for both 143 Nd and 144 Nd isotopes by multiplying this sediment flux to the concentration along the margin.
Nd IC of margins (in ε Nd ) Nd IC of rivers (in ε Nd ) Nd IC of atmospheric dusts (in ε Nd ) Nd concentration of margins (in pmol/kg) Nd concentration of rivers (in μg/g) River runo (in 10 15 g/yr) Fig. 2. Input maps applied to the model.(A) ε Nd map along the continental margin determined by Jeandel et al. (2007).This input is applied to the sedimentary remobilization source F sed .(B) Nd concentration map along continental margin (in µg/g) determined by Jeandel et al. (2007).This input is also applied to the sedimentary remobilization source F sed .(C) ε Nd map of river runoff given by Goldstein and Jacobsen (1987).(D) Nd concentration of river runoff (in ng/g) Goldstein and Jacobsen (1987); scale is non linear.(E) Interpolated ε Nd map of atmospheric dust (Grousset et al., 1988;Grousset et al., 1998;Jeandel et al., 2007).Dust particle fields are provided by Tegen and Fung (1995), and Nd concentration is set constant to 20 µg/g.(F) Runoff map prescribed by NEMO OGCM (in 10 15 g/an).Scale is non linear.
layer depth.Also, we multiply F sed by a vertical function equal to 1 in the first 1000 m, then exponentially decreasing at depth (f (z)).This vertical variation has been applied in accordance with iron modelling (Aumont and Bopp, 2006) and previous ε Nd modelling studies (Arsouze et al., 2007), because we suppose that the sediment dissolution is more important in surface waters than at depth (at the surface, dynamic activity is more vigorous, and biological processes are more important).The only available global estimation of the Nd Boundary Source flux is the "missing flux" calculated by Tachikawa et al. (2003).We therefore use their value of 8.0×10 9 g(Nd)/yr as a reference for our simulations.
In addition to Boundary Sources, dissolved river discharge (defined on continental margin points only) and atmospheric dust deposit are taken into account as Nd inputs in surface waters (first vertical level).

S(Nd
where F surf is the Nd flux of these two sources to the ocean (in g(Nd)/m 2 /yr).
For dissolved river discharge we use the climatological runoff applied to the dynamical model (Fig. 2f).Nd IC (Fig. 2c) and concentrations (Fig. 2d) in river inputs are estimated using the compilation of data provided by Goldstein and Jacobsen (1987).Using the runoff estimation provided by the NEMO model, and a subtraction of 70% of material in the estuaries (Sholkovitz, 1993;Elderfield et al., 1990;Nozaki and Zhang, 1995), we obtain a Nd dissolved flux from rivers of 2.6×10 8 g(Nd)/yr (Table 1).This value is lower, but the same order of magnitude, than the flux used by Tachikawa et al. (2003), equal to 5×10 8 g(Nd)/yr.Atmospheric dust flux is determined using monthly maps provided by Tegen and Fung (1995), in agreement with the method used for the PISCES model for nutrients (Aumont and Bopp, 2006).Nd IC of this source has been established using available data (Grousset et al., 1988(Grousset et al., , 1998)).However, in the areas where no data were available, the ε Nd of the dust is determined according to the ε Nd value for the region of origin of the dusts (Fig. 2e).These data are available from Jeandel et al. (2007).Nd concentration in atmospheric dusts was set constant at 20 µg/g (Goldstein et al., 1984).We considered a 2% Nd dissolution rate (Greaves et al., 1994), as used by Tachikawa et al. (2003), i.e. 98% of the Nd brought by atmospheric dusts to the ocean sinks with the lithogenic particles to the sediment without interacting with the scavenging cycle.It is important to note that some uncertainty remains concerning dissolution rates of atmospheric dusts, as published data vary from 2 to 50% (Tachikawa et al., 1999;Greaves et al., 1994).However, Table 1.Main characteristics of equilibrium partition coefficients and source fluxes for each experiment.Residence time of Nd in the ocean is calculated using the sum of flux (source or sink) and the total quantity of Nd simulated in the ocean: τ =Q Nd /(S(Nd T ) sed +S(Nd T ) surf ).Also, for comparison purpose, values of equilibrium partition coefficient and residence time from Siddall et al. (2008)  The re-dissolution of solid material deposited by rivers at estuary mouths is potentially a source of great importance (Sholkovitz, 1993;von Blanckenburg et al., 1996).However, this source is implicitly taken into account via the Boundary Source (F sed ), as we consider that all material deposited along the continental margin is made available for sediment re dissolution.
To balance the budget of Nd in the ocean, the only Nd sink considered is the burying of particles in the sediment at the bottom of the water column, i.e. scavenging.The particle-associated neodymium sinking at the deepest level thus leaves the model.Simulations are run until 143 Nd and 144 Nd oceanic concentrations reach an equilibrium state, which is when the sink balances global sources.

Simulation description
Currently the available observations of Nd concentrations in particulate phases, which allow constraint of the equilibrium partition coefficients (K) are scarce.In addition, the pool of particles generated by the PISCES model is biogenic, whereas data suggest that particulate Nd, is mainly adsorbed on the hydroxide phases coating the biogenic particles (Sholkovitz et al., 1994).Nd concentrations in iron or manganese hydroxides are at least an order of magnitude higher than those observed for biogenic particles (Bayon et al., 2004).The fact that PISCES does not take into account these particles might be critical to simulate realistic Nd dissolved-particulate interactions.Nevertheless, observations in the Atlantic and Austral basins, give an estimation of Nd p /Nd d values between 0.05 and 0.1 (Jeandel et al., 1995;Tachikawa et al., 1997;Zhang et al., 2008).
Using the same model for their 231 Pa and 230 Th simulations, Dutay et al. (2009) have shown the necessity to impose a variation of the partition coefficient as a function of particle stock, in order to simulate realistic vertical profiles for both the dissolved and particulate phases.In addition, the scavenging of these tracers had to be controlled by the small particle pool, and thus the value of their partition coefficients with this pool had to be greater than those applied to the big particles.We adopt this approach for our Nd simulation and therefore set equilibrium partition coefficients for the big particle pools lower than those for the small particle pool.
Lastly, contrary to 231 Pa and 230 Th isotopes, for which scavenging coefficients can vary from several orders of magnitude depending on the particulate pool (Dutay et al., 2009), there is no current evidence from Nd data that preferential scavenging occurs when biogenic silica or carbonate dominate the particle pool.Thus we only differentiated small and big particles.The same equilibrium partition coefficient value is taken for each pool of big particles, considering their global averaged concentrations (POCb, BSi, CaCO 3 and litho).Characteristics of sensitivity tests performed on K values performed are summarized in Table 1 (EXP3, EXP4 and EXP5).
The first experiment (EXP1) takes into account only dissolved river discharge and atmospheric dust (surface sources), which are the classical sources generally considered in trace element budget studies.In addition only small and lithogenic particles are taken into account to simulate the vertical cycle.We take K POMs =1.4×10 7 and K litho = 2.3×10 6 that correspond to Nd P /Nd d =0.001 and 1×10 −4 respectively, and K=0 for remaining particle pools (POMb, CaCO 3 , BSi, cf.Table 1).Although these values are much lower than the available data (Tachikawa et al., 1997), they are taken close to the same order of magnitude as in Siddall et al. (2008), and normalized to unit Nd P /Nd d value.
In a second experiment (EXP2), we consider the same configuration as the first simulation, but we added the sediment remobilization source (Boundary Source).
EXP3 and EXP4 test the sensitivity of the vertical cycling, considering the same sources of Nd as those considered in EXP2.For the third experiment (EXP3), we enhanced the small particle role relative to EXP2 by increasing the value of the equilibrium partition coefficient on small particles to K POMs =7.7×10 7 (which is equivalent to Nd P /Nd d =0.005).In the fourth experiment (EXP4) we have inserted the big particles (POMb, CaCO 3 , BSi) that were neglected in the first three experiments, the other equilibrium partition coefficients being identical to those of EXP2.The equilibrium partition coefficient value for big particles in EXP4 is set two orders of magnitude lower than for small particles, considering their respective concentrations in the PISCES model, i.e.K POCb =2.6×10 4 , K BSi =1.8×10 4 and K CaCO 3 =7.8×10 4 that correspond to Nd p /Nd d =1×10 −5 (Table 1).
Finally, in the last experiment (EXP5) we optimized the characteristics of both our dynamical and biogeochemical models according to the information gained through the previous tests, in order to produce the most realistic simulation.BE flux is adjusted to F(Nd T ) sed =1.1×10 10 g(Nd)/yr and equilibrium partition coefficients are slightly reduced compared to EXP4 for POMb, CaCO 3 , BSi and litho, and unchanged for the small particles (Table 1).Depending on the estimations, only between 3 and 5% (Milliman and Syvitski, 1992) dissolution of Nd material is required to explain our total Boundary Source flux S(Nd T ) sed =1.1×10 10 g(Nd)/yr.

Results
For each experiment, both Nd IC and concentration are shown along meridional sections in the Pacific Ocean (averaged zonally between 120 • W and 160 • W) and along the western boundary of the Atlantic Ocean (Figs. 3 and 5), where the water mass ε Nd have been characterized.The first experiment (EXP1) produces a satisfactory ε Nd distribution in surface waters (Fig. 7), but is much too homogeneous in the deep ocean (Figs. 3,5 and 8), particularly in the Atlantic basin, where the signature of the main water masses (NADW, AAIW, AABW) are not correctly sim-ulated.Moreover, the Nd concentrations are at least an order of magnitude too low compared to the observations (average global value of [Nd] (model)=2.3pmol(Nd)/kg, [Nd] (data)=21.5pmol(Nd)/kg and only 12% of modelled Nd concentrations fit into the ±10 pmol/kg envelop when compared to the data, Fig. 5 and Table 1).
In the second experiment (EXP2, BE added), simulated ε Nd values in surface waters are still in agreement with the observed data, and its distribution in the deep Atlantic Ocean is closer to the observed data.Indeed, the water masses in this basin are now simulated with very realistic ε Nd values (Fig. 3 and 6).It should also be noted that this simulation produces a realistic increase of concentrations along the thermohaline circulation, in agreement with the observations.
Increasing the value of the equilibrium partition coefficient on small particles in the third experiment (EXP3) leads to a strengthening of the vertical cycling effect on simulated ε Nd and Nd concentrations.It results in an undesirable vanishing of NADW isotopic signature along its pathway in the Atlantic Ocean (Fig. 3).However it improves the ε Nd distribution in the deep Pacific Ocean where values becomes more radiogenic (Fig. 8), but also leads to unrealistic high values at intermediate depths.Increasing vertical cycling enhances the sedimentation process (sink of Nd), and excessively reduces the simulated Nd concentrations ([Nd] (model) = 6.8 pmol(Nd)/kg and only 21% fit into the [Nd] (model) = [Nd] (data) ±10 pmol/kg envelop, Fig. 5 and Table 1).The decrease of total Nd concentration with a constant source flux leads to a reduction of the residence time τ (defined as τ =(F(Nd) surf +F(Nd) sed )/Nd T ) from ∼700 years in EXP1 and 2 to 150 years in EXP3 (Table 1).
Taking into account the influence of fast sinking biogenic big particles (POMb, CaCO 3 , bSi; EXP4) yields results close to those observed when increasing the equilibrium partition coefficient for small particles (EXP3): the oceanic residence time of the tracer decreases (125 years), and both Nd concentration and Nd IC in the water column are more homogeneous.Simulated Nd IC distribution in the deep Pacific Ocean is in better agreement with the observations than in EXP2, but Nd concentration becomes too low.Most of all,  (Piepgras et Wasserburg, 1982;Piepgras et Wasserburg, 1983;Stordal et Wasserburg, 1986;Piepgras et Wasserburg, 1987;Piepgras et Jacobsen, 1989;Jeandel, 1993;Jeandel et al., 1998;Lacan et Jeandel, 2001, Amakawa et al., 2004).Colour scales are the same for all simulations and data, but different for each section.For each basin, characteristic profiles are numbered and located on the associated map, and at the bottom of each transect, for each experiment.Main characteristics for each experiment are summarized in Table 1.
this simulation succeeds in keeping a realistic isotopic signature of AAIW and NADW in the Atlantic basin, although generating an AABW which is slightly too non-radiogenic (Fig. 3).
Finally in the last experiment (EXP5) in which intensity of the BE flux and values of partition coefficients of large particles have been adjusted, we simulate ε Nd distribution in good agreement with the data (71% fit into the ε Nd (model)=ε Nd (data) ±3 ε Nd envelop, Figs. 3,4,7 and 8).In particular, inter-basin ε Nd gradients (Fig. 8) and isotopic signature of the main water masses (Fig. 3) are correctly reproduced, as well as surface isotopic signatures (Fig. 7).The increase of the Nd sink has been compensated by an increase of the sediment remobilization process (F sed ) to produce more coherent concentrations (Figs. 5,6 and 10,average   in the [Nd] (data) ±10 pmol/kg envelop).As observed in EXP4, the scavenging onto big particles leads to more coherent vertical gradients, even if concentrations in the surface layer remain systematically too low (Figs. 5 and 9) while concentrations generated at high northern latitudes are too high (Figs. 5,6 and 9).Residence time is 360 years.

Sensitivity of Nd oceanic cycle to sources
The above results provide valuable information about the role (or impact) of the different Nd sources on the oceanic Nd IC and concentration distributions.This information allows a better understanding of one important aspects of the Nd paradox: how can we explain the observed Nd IC gradient along the global thermohaline circulation despite a relatively small increase in Nd concentration?
In our first experiment, we first applied dissolved river discharge and atmospheric dust (F surf -EXP1).This lead to simulated ε Nd values that are close to the existing data at surface depths (0-200 m, Fig. 7), suggesting that these sources could control the isotopic Nd distribution in the upper ocean.However, it produces an ε Nd distribution that is too homogeneous in the deep ocean, too negative in the deep Pacific, and provides a very poor simulation of water masses signature in the Atlantic Ocean.A possible way to improve the results in the Pacific Ocean would be to enhance vertical cycling (either by increasing the equilibrium partition coefficient for the small particles or by inserting the big particles in the scavenging model), in order to export a radiogenic surface waterlike signature at depth.However, this action yielded an increase in Nd sedimentation and subsequently a decrease in global Nd concentration.Since the simulated Nd concentrations in EXP1 are already an order of magnitude lower than the observed concentrations, increasing the vertical cycling does not help to reconcile both Nd IC and concentration distributions, as much as the sources are restricted to the surface only.
Another option for improving the simulated deep Pacific and Atlantic ε Nd distributions was to consider an additional source of Nd to the oceanic reservoir.Sediment redissolution effects along the continental margin (Boundary Source) have been tested in the second experiment (EXP2).Sediment flux (S(Nd T ) sed ) intensity was estimated using first the "missing flux" proposed by Tachikawa et al. (2003) and then adjusted in EXP5.Taking into account this source largely improved our simulation of the Nd oceanic cycle, generating Nd concentrations closer to the observed data, though still an order of magnitude lower when compared to the data in surface and intermediate waters.It also greatly improved the simulation of the isotopic composition of the main water masses in the Atlantic Ocean.The Nd oceanic budget simulated here suggests a large predominance (about twenty times higher than cumulated other sources) of the role of the sediment re-dissolution source (S(Nd T ) sed ), meanwhile surface sources (S(Nd T ) surf ) appeared predominant for constraining surface water isotopic signatures.Other experiments, not reported here, showed that considering BE only, excluding dissolved river and dust inputs, have led to less realistic simulations of the surface ε Nd distribution.As a source that influences Nd concentration and Nd IC at depth, we can preclude submarine groundwater as acting in the Boundary Sources, because their influence is mainly limited to the upper 200 m (Johannesson and Burdige, 2007).
Changing the solubility of atmospheric dust entering seawater or reducing the subtraction of dissolved material in the estuaries (which would contradict the field or experimental results) may change the S(Nd T ) surf value.However, this latter value still remains small compared to S(Nd T ) sed value, and the change does not significantly modify the Nd distribution (sensitivity tests not presented here).Boundary Source can therefore be considered as a major source of Nd in the ocean that is fundamental to simulate a realistic Nd oceanic cycle.Rickli et al. (2009).Colour scale is non linear.Main characteristics for each experiment are summarized in Table 1.[Nd]

Sensitivity of Nd oceanic cycle to vertical cycling: small vs. big particles
The second experiment successfully simulated reasonable distributions of both Nd IC and concentrations.However, some discrepancies still remain, particularly in the deep Pacific where simulated deep waters had underestimated Nd IC and overestimated concentrations.Enhancing vertical cy-cling in order to homogenize the vertical column in the Pacific is therefore pertinent, though it should maintain the outcome in the Atlantic Ocean.As the lateral isopycnal transport is more intense in the Atlantic than in the Pacific Ocean, vertical processes in this last basin might be more sensitive.
Basically, the equilibrium partition coefficient establishes the relative influence between dynamical lateral transport and vertical scavenging.Enhancing the influence of vertical [Nd] (in pmol/kg) cycling, either by increasing the partition coefficient for small particles (EXP3) or by inserting big particles in the reversible scavenging model (EXP4), increases the removal of Nd from the water column and reduces our global Nd concentration to more unrealistic values (cf.Table 1 ,Figs. 5,6,9 and 10).Both EXP3 and EXP4 successfully generated the required reduction of Nd concentration, and the formation of more radiogenic waters in the deep Pacific Ocean (Figs. 3,4 and 8).However they differed in their performance in the Atlantic Ocean.EXP3 generated a significant degradation of simulated ε Nd characteristics in the Atlantic Ocean, while EXP4 preserved and even improved agreement with the data in the Atlantic basin.
These results are consistent with those of Siddall et al. (2008), which demonstrated the importance of vertical cycling (reversible scavenging process) in reconciling Nd IC and concentration distributions, and therefore resolving the "Nd paradox".Also, this highlights the importance of different pools of particles size in this model for setting the characteristics of the Nd isotopic distribution, and more precisely the potential role of fast sinking particles.

The "Nd paradox"
Based on information gained from previous sensitivity tests on sources and vertical cycling (EXP1 to EXP4), the configuration of our last experiment (EXP5) was adjusted to obtain the best agreement between simulations and data.
We have considered the surface sources of EXP1 (S(Nd T ) surf =2.6×10 8 g(Nd)/yr) whereas the sediment remobilization has been slightly increased (S(Nd T ) sed =1.1×10 10 g(Nd)/yr), in order to compensate for the decrease in Nd concentration caused by the insertion of big particles in the reversible scavenging model.Associated with this later source, the sink of Nd along the continental margin (corresponding to the Boundary Scavenging process), even if less important as the source, represents as much as 64% of the global Nd sink (Fig. 11), confirming that BE acts as both a major source and sink for oceanic Nd, as deduced from the observations (Lacan and Jeandel, 2005).
Due to the limited number of sensitivity experiments that could be conducted due to time constraints, this version is not necessarily optimal, but it shows some improvements compared to the four other experiments.The resulting simulated Nd IC is in excellent agreement with the observed data (Figs. 3,4,7 and 8), though concentrations in the surface ocean are low (Figs. 5,6 and 9).We attribute this bias to the particulate fields generated by the biogeochemical model PISCES.In this work, as in the 231 Pa and 230 Th simulations proposed by Dutay et al. (2009), and consistent with either experimental or field observations, the small particle pool drives the vertical cycling in the water column.It is characterized here with a higher equilibrium coefficient for small particles than for the big particle pool.However, small particle fields generated by PISCES are underestimated at depth (Dutay et al., 2009) leading to an overestimation of the Nd concentrations (as Nd is less scavenged).Improving the small particle concentration at depth in the PISCES model would likely help to smooth the vertical modelled gradient, recover deep Pacific radiogenic waters in our simulation (Figs. 3 and 8), and better simulate element cycles in the ocean, which is an important goal.Nevertheless, in our model, big particles also play an important role in reconciling deep ocean Nd concentration while keeping the Nd IC as a water-mass tracer property.
The residence time of Nd in the ocean in this last experiment is 360 years, which is in the lower range of previous estimations of ∼300 to ∼600 years (Piepgras and Wasserburg, 1983;Jeandel et al., 1995;Tachikawa et al., 1999Tachikawa et al., , 2003;;Siddall et al., 2008).This value is still consistent with the conservativity of ε Nd within the major oceanic water masses.All together, these factors suggest that i) the Boundary Source is a major primordial source for Nd oceanic cycle (∼95% of the global sources), ii) Boundary Scavenging, resulting from vertical cycling in margin areas, is a major sink (64% of the global sink), and iii) vertical cycling in the open ocean is essential for the redistribution of Nd and its IC within the ocean interior.While the two first points provide an explanation for the BE observed in field data, the three points provide a means to reconcile the "Nd paradox".The first estimation of the order of magnitude of the "missing flux" proposed by Tachikawa et al. (2003) is confirmed by our model results.The reversible scavenging is also essential to explain the Nd nutrient-like profile in the water column and the water mass tracer property of Nd IC, in agreement with recent REE modelling studies (Nozali and Alibo, 2003;Siddall et al., 2008;Oka et al., 2008).
However, if Boundary Sources are found to explain the large Nd IC variation observed along the global thermohaline circulation, without a large Nd concentration change (due to the associated Boundary Scavenging, Fig. 11), the question of why there is an apparent contradiction between conservative Nd IC and non conservative Nd concentrations remains to be solved.This critical point should be addressed in another dedicated study, so as to fully resolve the "Nd paradox".

Conclusions
The objective of this study was to use an ocean model as a tool to better understand the Nd oceanic cycle and to work toward the resolution of the "Nd paradox", (i.e.estimating the sources and sink of the element, and qualifying the processes acting on the element's distribution within the reservoir).We simultaneously simulated both Nd isotopic composition and concentration using a prognostic coupled dynamical/biogeochemical model and used a reversible scavenging model to simulate the Nd transformation and sink into the ocean.We have also implemented for the first time a realistic calculation of the Nd sources, with an explicit representation of sedimentary remobilization along continental margins (source of the Boundary Exchange process) as well as dissolved river discharge and atmospheric dust sources (surface sources).
In accordance with previous results of Siddall et al. (2008) andOka et al. (2008), vertical cycling needs to be invoked (surface removal combined with remineralisation at depth) to correctly simulate Nd isotopic composition and concentration distributions.As yet, we are unable to verify the validity of our partition coefficients because very little data are available.A strong recommendation deduced from this work is to improve our knowledge of the dissolved/particle distribution of the geochemical tracers, which the GEOTRACES program should make possible.Our performance in simulating Nd concentration is also still limited by our models.In particular, the low concentration in the small particles field simulated at depth by the PISCES model is very likely the cause for overestimated Nd concentration surface to depth gradients.This work also suggests that regarding the parameterization used and the particle distribution generated by the biogeochemical model, it is important to consider different kinds of particles, and especially their sinking velocity, in setting the characteristics of the Nd isotopic distribution.This study can therefore be used as a basis for improving the representation of the complexity of observed particle field distribution in the model, such as considering a whole continuous spectrum of particle size (Gehlen et al., 2006).
The simulations performed strongly suggest that sediment dissolution is the main Nd source to the oceanic reservoir, with a flux of 1.1×10 10 g(Nd)/yr (95% of the total Nd source), whereas the associated Boundary Scavenging process represents up to 64% of the total Nd sink.Combination of these dominant source and sink at the margins are consistent with the concept of Boundary Exchange deduced from field observations.We suggest that the Boundary source is likely the "missing source term" mentioned by Tachikawa et al. (2003) that could reconcile both Nd isotopic composition and concentration oceanic budgets.The fact that this potential source has only been recently considered, despite the large amount of Nd available from solid material from river discharge deposited on the continental margin, may be explained by the very small amount of The high computational cost for one simulation was a limiting factor for sensitivity tests on sources and equilibrium partition coefficients.Optimal values of parameters of the model may not be those used for experiment EXP5, but we are confident that our main results are robust.One easily envisaged improvement could be the use of the Transport Matrix Method (Khatiwala et al., 2005) for this model, to better optimize these parameters.
This study highlights the predominant role of Boundary Exchange on Nd oceanic cycle, and suggests Nd isotopes as a powerful tool to quantify margin source to the ocean, which may have significance for the general ocean chemistry.However, we have no indication of chemical, biological or physical processes that act on this sediment re-dissolution.Particularly, one of the main shortcomings resides in the assumption that the input from the margin is set geographically constant, whereas some factors can act on the source (oxygen concentration, currents velocity, turbulence, temperature, organic matter deposition, etc...), leading to potentially large variations.It then appears to be quite important for the geochemical community to pay attention at continental margins/open ocean interfaces to determine if the "Boundary Scavenging" (observed for Be and Pa/Th, Anderson et al., 1990) is compensated by lithogenic inputs for other chemical tracers, and how this input is materialized (e.g.cold hydrothermalism on margins, remobilization and dissolution of sediments after resuspension, early diagenesis or both, water penetration in sediments, etc. . .).
This implies the need for a multi-tracer approach, in particular with highly reactive elements such as 231 Pa and 230 Th.Also, it seems inevitable to improve and adapt biogeochemical models to the simulation of trace elements, including iron manganese hydroxides and oxide crusts that are hypothesized to be main carriers of these tracers.
512638 is the present averaged earth value(Jacobsen and Wasserburg, 1980).
Nd and [Nd] Figures 4 and 6 are scatter plots showing model and observations and Figs. 7, 8, 9 and 10 display horizontal maps of ε Nd in surface (0-200 m) and at depth (2500-4000 m), and Nd concentration averaged over the same depth ranges, respectively.

Fig. 3 .
Fig. 3. Vertical ε Nd sections for all simulations, in the Pacific basin averaged between 120 • W and 160 • W (upper panel), and along the section of the Atlantic basin (lower panel).Superimposed circles the data from the compilation by F. Lacan and T. van de Flierdt; http://www.legos.obs-mip.fr/fr/equipes/geomar/results/databasemay06.xls(Piepgras et Wasserburg, 1982;Piepgras et Wasserburg, 1983;Stordal et Wasserburg, 1986;Piepgras et Wasserburg, 1987; Piepgras et Jacobsen, 1989;Jeandel, 1993;Jeandel et al., 1998;Lacan et Jeandel, 2001, Amakawa et al., 2004).Colour scales are the same for all simulations and data, but different for each section.For each basin, characteristic profiles are numbered and located on the associated map, and at the bottom of each transect, for each experiment.Main characteristics for each experiment are summarized in Table1.

Fig. 4 .
Fig. 4. ε Nd model/data cross plot for each experiment as a function of depth (colour code).Circles represent the data located in the Atlantic basin; squares are data in the Indian basin and triangles are data in the Pacific basin.Red line is the linear interpolation line.Black lines represent the lines ε Nd (model)=ε Nd (data) and ε Nd (model)=ε Nd (data)±3 ε Nd , which represents a characteristic value used to define a good agreement between the model and the data.Main characteristics for each experiment are summarized inTable 1.Data are from the compilation by F. Lacan and T. van de Flierdt; http://www.legos.obs-mip.fr/fr/equipes/geomar/results/databasemay06.xls(see references therein).

Fig. 7 .
Fig. 7.Horizontal ε Nd maps averaged between 0 and 200 m, for all the experiments.Superimposed circles represent the data at the same depth, from the compilation by F. Lacan and T. van de Flierdt; http://www.legos.obs-mip.fr/fr/equipes/geomar/results/databasemay06.xls(see references therein) andRickli et al. (2009).Colour scale is non linear.Main characteristics for each experiment are summarized in Table1.
have been added.
Duce et al. (1991)2003)ll, this value is lower thanTachikawa et al. (2003)who considered the earlierDuce et al. (1991)atmospheric dust deposit estimation, leading to a flux of 5×10 8 g(Nd)/yr.Additional sensitivity tests on both dissolution rate ratios for atmospheric dusts and Nd dissolved river discharge (not shown here) do not significantly change the results and conclusions presented here.