Neodymium isotopes in the ocean model of the Community Earth System 1 Model ( CESM 1 . 3 ) 2

Abstract. Neodymium (Nd) isotope ratio (eNd) is a quasi-conservative water mass tracer and has been used increasingly as paleoclimate proxy to indicate the past evolution of ocean circulation. However, there are many uncertainties in interpreting eNd reconstructions. For the purposes of direct comparison between climate models and proxy reconstructions, we implement Nd isotopes (143Nd and 144Nd) in the ocean model of the Community Earth System Model (CESM). Two versions of Nd tracers are implemented: one is the "abiotic" Nd in which the particle fields are prescribed as the particle climatology generated by the marine ecosystem module of the CESM under present day forcing; the other is the "biotic" Nd that is coupled with the marine ecosystem module. Under present day climate forcing, our model is able to simulate both Nd concentrations and eNd in good agreement with available observations. Also, Nd concentration and eNd in our model show similar sensitivities to the total boundary source and the ratio between particle related Nd and dissolved Nd as in previous modeling study (Rempfer et al., 2011). Therefore, our Nd-enabled ocean model provides a promising tool to study past changes in ocean and climate.


Introduction
Radiogenic 143 Nd is produced by the radioactive decay of 147 Sm with decay halflife of 106 billion years (Lugmair, 1974).During magma formation, Nd is more likely to enter magma than Sm, therefore, continents have lower Sm/Nd or 143 Nd/ 144 Nd compared to mantle (melt residue) and the bulk of earth.The difference of 143 Nd/ 144 Nd between continents and the bulk of earth increases with the age of the continent as 143 Nd/ 144 Nd in younger continents is more similar to the mantle.Therefore, younger (older) continents have higher (lower) 143 Nd/ 144 Nd, which is more radiogenic (unradiogenic) (Goldstein and Hemming, 2003).Nd isotopic ratio ( 143 Nd/ 144 Nd) relative to the "bulk earth" value is reported as εNd: , where ( 143 Nd/ 144 Nd)bulkearth is 0.512638 (Jacobsen and Wasserburg, 1980).Due to the different ages of continental crust, εNd in continental crust varies geographically (Albarède and Goldstein, 1992).The general feature consists of the two extremes, with the most unradiogenic values (minimum) in the North Atlantic (-10 to-14), the most radiogenic values (maximum) in the Pacific (-3 to -4), and intermediate values in the Indian and Southern Ocean (-7 to -10).Seawater derives its εNd value mainly through weathering and erosion of continental crust (Piepgras et al., 1979).
Therefore, different water masses form from different locations have different εNd values.For example, εNd of North Atlantic Deep Water (NADW) is around -13.5, whereas εNd of Antarctic Intermediate Water (AAIW) and Antarctic Bottom Water (AABW) is around -8.In the Atlantic, εNd covaries with salinity (von Blanckenburg, 1999) and behaves as quasi-conservative water mass mixing tracer (Goldstein and Hemming, 2003;Piepgras and Wasserburg, 1982).
Unlike quasi-conservative εNd, Nd concentration shows a nutrient-like behavior as it increases with depth and also along the circulation pathway (Bertram and Elderfield, 1993).The decoupling of εNd and Nd concentration, or the so-called "Nd paradox", can be explained by reversible scavenging (Bacon and Anderson, 1982;Siddall et al., 2005) in internal Nd cycling (Siddall et al., 2008).
" Nd = [( ( 143 Nd/ 144 Nd) sample ( 143 Nd/ 144 Nd) bulkearth ) 1] ⇥ 10 4 εNd has been increasingly used in paleoceanographic studies (e.g.Piotrowski et al. 2004;Gutjahr et al. 2008;Roberts et al. 2010;Piotrowski et al. 2012) because of its ability to trace different water masses.Also, biological fractionation of Nd isotopes are negligible (Goldstein and Hemming, 2003).However, our knowledge about Nd is limited for a reliable interpretation of εNd for past ocean changes.For example, interpretation of the Atlantic εNd reconstructions is based on the assumption of the stable north (NADW) and south (AAIW and AABW) εNd endmembers.NADW is a mixture of low εNd water from the Labrador Sea (<-20) and high εNd water from the Norwegian and Greenland Sea (-7 to -10).Therefore, small changes in deep water formation during the last deglaciation, which is highly uncertain (Crocket et al., 2011;Dokken and Jansen, 1999;Labeyrie et al., 1992), will result in large changes in εNd of NADW (van de Flierdt et al., 2016).In addition, the magnitude and isotopic composition of Nd in sources, which have been suggested to be changing in the past (e.g.: Grousset et al. 1998;Harris and Mix 1999;Amakawa et al. 2000;Lézine et al. 2005;Wolff et al. 2006;Rickli et al. 2010), may also influence εNd in seawater (Tachikawa et al., 2003).Therefore, incorporating Nd isotopes into climate models can help to improve our understanding of Nd cycling.Previous modeling efforts of Nd have made much progress (Arsouze et al., 2009;Rempfer et al., 2011;Siddall et al., 2008).Modeling studies also suggest that εNd end-member changes are relatively small compared with εNd changes resulted from water mass distribution (Rempfer et al., 2012a) and glacial-deglacial εNd variations are hard to be obtained by changes in Nd sources alone (Rempfer et al., 2012b).
Currently, many uncertainties and controversies in our understanding of past ocean evolution involve the interpretation of Nd reconstructions (e.g.Huang et al., 2014;Pahnke et al., 2008;Xie et al., 2012).Therefore, it is crucial to incorporate Nd isotopes into climate models such that model simulation and proxy record can be compared directly.This direct model-data comparison will help us to better interpret the proxy records and, furthermore, understand past ocean circulation changes.This paper is the documentation of the implementation of neodymium isotopes, 143 Nd and 144 Nd, in the ocean model of the Community Earth System Model (CESM) (Hurrell et al., 2013).2006), with a successful simulation of both Nd concentration and εNd in good agreement with observations.The parameters tuned in Rempfer et al., (2011) are based on the compilation of observations up to September 2011 by Lacan et al., (2012).Nd isotopes are included in the GEOTRACES program (Mawji et al., 2014).
van de Flierdt et al., (2016) complied available Nd data up to January 2016, which includes data collected by the GEOTRACES program and additionally approximately 1,000 published data points collected outside GEOTRACES.This compilation is more than double the amount of the previous data compilation by Lacan et al., (2012).Our study uses the new database to tune model parameters.Also, our Nd module is coupled with the marine ecosystem model of the CESM (eco_Nd).In addition, we also implement an abiotic Nd (abio_Nd, similar to Rempfer et al., (2011)).Using a prescribed particle flux field, the abio_Nd can be run without the marine ecosystem module and thus has a much-reduced computation cost.Most importantly, the abio_Nd can be compared with the eco_Nd to separate the effect of circulation change and biological change on Nd.These two Nd implementations will be added to the code trunk of the current ocean model of the CESM, which will make them available to other scientist and will allow them to be maintained as CESM evolves.
This paper serves as a reference for future studies using Nd isotopes in the CESM.We will describe the model and the details of the implementation of the Nd isotopes in Section 2. The experimental design of the test simulations is described in Section 3. The results of the parameter tuning process, the comparison between the simulated Nd concentrations and εNd with observations, and the model sensitivities to two parameters are discussed in Section 4.

Physical Ocean model
The implementation of Nd is based on the code of CESM, version 1.3.CESM is a state-of-art coupled model and many of the papers describing model component  (Danabasoglu et al., 2012).The experiments in this study are carried out using the fully active and isotope-enabled POP2 coupled to the data atmosphere, land, ice and river runoff under the normal year forcing from CORE-II data (Large and Yeager, 2008).The model has a nominal horizontal resolution of 3° and 60 vertical layers, with a 10-m resolution in the upper 200m, increasing to 250m below 3000m.The ocean-alone model at 3° resolution is used due to its low computational cost, allowing us to carry out extensive parameter test simulations.Future applications of the Nd isotopes should use the scientifically validated 1° resolution of the CESM.

Biogeochemical component
The biogeochemical variables used in the Nd isotopes implementation (particle fluxes: CaCO3, opal, POC, and dust fluxes) are generated by the marine ecosystem model in the CESM (Moore et al., 2013) through the ecosystem driver (Jahn et al., 2015).Simulated annual mean particle (CaCO3, opal, and POC) fluxes leaving the euphotic zone at 105m (Fig. 1, a~c) show patterns and magnitudes similar to those in satellite observations (Fig. 7.2.5 and 9.2.2 in Sarmiento and Gruber 2006).Surface dust deposition is taken from the ecosystem module, which is prescribed monthly surface dust flux from Luo et al., (2003) (Fig. 1d).The remineralization scheme of particle is based on the ballast model of Armstrong et al., (2002).Detailed parameterizations for particle remineralization are documented in Moore et al., (2004) with temperature dependent remineralization length scales for POC and opal.

Nd isotopes implementation
The Nd isotopes ( 143 Nd and 144 Nd) are added as optional tracers, which can be turned on at case build time as some other passive tracers (e.g., ideal age, carbon isotopes (Jahn et al., 2015) and water isotopes (Zhang, 2016)).We implement both abio_Nd and eco_Nd, the latter of which is coupled with the marine ecosystem model and therefore requires the ecosystem model to be turned on at the same time.The only difference between abio_Nd and eco_Nd is that abio_Nd uses a set of prescribed annually averaged dust, opal, POC, and CaCO3 fields that are generated from the ecosystem module offline (Fig. 1), while eco_Nd uses these fields simultaneously computed from the ecosystem module.
The Nd module is implemented following Rempfer et al., (2011)

Nd sources
Dust deposition over the ocean surface is one of the Nd sources to the ocean.
The surface dust source, Sdust(i,j) (g m -3 s -1 ), is applied to the surface layer of ocean grid, and can be calculated as: (3) Here, surface dust flux, Fdust(i,j) (g cm -2 s -1 ), is obtained from the ecosystem module of the CESM (Fig. 1d); global mean Nd concentration in dust is 20 μg/g (Cdust) (Goldstein et al., 1984;Grousset et al., 1988Grousset et al., , 1998)); 2% (βdust) of the total Nd in the dust is released into ocean (Greaves et al., 1994); dz1 (m) is the thickness of the surface layer of the ocean grid.The annual total Nd from dust, fdust, is 2.1×10 8 g River runoff also provides Nd to the ocean.The river source of Nd is applied at the surface layer of the ocean.River source, Sriver(i,j) (g m -3 s -1 ), can be obtained from: (4) Here, ROFF(i,j) (kg m -2 s -1 ) is the river runoff from the coupler of the CESM.The simulated global annual river discharge is 41,584 km 3 /yr, similar to the observational estimate of 42,439 km 3 /yr in (Goldstein and Jacobsen, 1987).Nd concentration in river runoff, Criver (g kg -1 ) is extrapolated from the river Nd concentration data (Goldstein and Jacobsen, 1987); we assume 70% (γriver) of Nd in rivers is removed in estuaries, following Rempfer et al., (2011); dz1 (m) is the thickness of the surface layer of ocean grid.The annual total Nd source from river runoff, friver, is 1.3×10 9 g/yr, which is larger than the reported values of 5×10 8 g/yr (Goldstein and Jacobsen, 1987).The difference may be caused by the extrapolation of Nd concentration from the original data from Goldstein and Jacobsen, 1987.Individual river sources for 143 Nd and 144 Nd can be calculated from the prescribed IR field following Jeandel et al. ( 2007) (Fig. 2b).
Weathering of continental crust is another source of Nd.This boundary source of Nd is applied to all continental margin grids above 3,000 m (Fig. 2b).We assume a globally uniform boundary source per unit area (fboundary/Atot), where fboundary (g/yr) is the total boundary source and Atot (m 2 ) is the total area of continental margin.fboundary is a tuning parameter, as in Rempfer et al., (2011).
Boundary source used in Arsouze et al., (2009) is assumed to be exponentially decreasing with depth but observations from GEOTRACES data suggest no obvious depth dependence (van de Flierdt et al., 2016).Boundary source, Sboundary(x,y,z) can be calculated by Eq. 5, where dzk is the thickness of the ocean grid layers.Individual ( 5)

Reversible scavenging and Nd sink
Reversible scavenging is the process of Nd adsorption onto sinking particles (POC, opal, CaCO3, and dust) and desorption during particle dissolution at depth, which transports Nd downwards (Siddall et al., 2008).Total Nd can be separated into dissolved Nd phase ([Nd]d) and particle associated Nd phase ([Nd]p) (Eq.( 6)).
Particle associated Nd can be further separated into Nd associated with different particle types (POC, CaCO3, opal, and dust)(Eq.( 7)).At the bottom grid, Nd associated with undissolved particles is removed from the ocean through sedimentation, which is the sink for Nd budget.
(6) (7) The ratio of between dissolved [Nd]d and [Nd]p is given by the "equilibrium scavenging coefficient", K: where i refers to different particle types (i = POC, CaCO3, opal, and dust), j refers to the different Nd isotopes (j = 143 Nd and 144 Nd), and here is another tuning parameter. ! is the ratio between the global average particle concentration ( !, Table1) and the average density of seawater (1024.5 kg m -3 ).We assume the dissolved Nd and the particle associated Nd are in equilibrium as in other studies (Arsouze et al., 2009;Rempfer et al., 2011;Siddall et al., 2008).Therefore, in each grid cell, the ratio between [Nd]p and [Nd]d can be obtained from Eq. ( 9).
where Ri(x,y,z) is the ratio between the particle concentration, Ci(x,y,z), and the density of seawater.Ci(x,y,z) can be calculated from particle fluxes Fi(x,y,z), which are provided by the ecosystem module, by applying a settling velocity (w) (Ci = Fi/w).We assume a uniform settling velocity of 1000 m yr -1 for all four kinds of particles.This is the velocity of the small particles, which drives the vertical cycling of isotopes (Arsouze et al., 2009;Dutay et al., 2009;Kriest, 2002).Isotopic fractionation between 143 Nd and 144 Nd during the reversible scavenging process is neglected as in Rempfer et al., (2011) because of similar molecule mass of 143 Nd and 144 Nd.We, therefore, apply the same Ki to 143 Nd and 144 Nd.
The reversible scavenging process acts as the internal cycling of Nd, which transports Nd from shallow layers to deep layers.This process can be quantified as a source term in the Nd equation where w is the settling velocity of particles (1000m yr -1 ) and [Nd]p is Nd associated with particles, which can be calculated at every time step using Eq. ( 6), ( 7) and ( 9) (combined as Eq. ( 11)). (11) Therefore, the conservation equation for Nd can be written as such that the Nd concentration change is determined by three source terms in Eq.

Experiments
Geosci Following Rempfer et al., (2011), our Nd model is tuned with two parameters: fboundary and , in the abio_Nd implementation under present-day climate forcing.The tuning in the abio_Nd implementation gives us a great computational efficiency because the ecosystem module can be turned off.Yet, the parameters tuned for abio_Nd should also apply to the eco_Nd since the particle fields used in the reversible scavenging process for abio_Nd are the climatology taken from the equilibrium ecosystem module under the same climate forcing.
The Nd concentrations ( 143 Nd and 144 Nd) are initialized from zero and each experiment is integrated for 4,000 model years (experiments with fboundary = 1×10 9 and 2×10 9 are initiated from 3,000 model years of the experiment with fboundary = 3×10 9 and then integrated for another 1,300 model years each).Nd inventory has reached equilibrium in most of the experiments at the end of the simulation.Those that do not reach equilibrium show unreasonable Nd concentrations or εNd and drift further and further from observation as the model integrates (e.g.cases with ), and therefore are terminated at some point.

Results
Geosci

Cost function minimum cases
We first study the dependence of the cost function J[Nd]d on fboundary and (Fig. 3, solid line contours).Given a fixed fboundary, for an increase in reaches the minimum at fboundary = 2×10 9 g yr -1 and -4 , which is marked by the pink square in Fig. 3.This pair of parameters of optimal J[Nd]d produces a [Nd]d pattern in reasonable agreement with the observation, as shown in a transect across the Atlantic and the Pacific in Fig. 4a (the transect is indicated in Fig. 2a).
With this set of parameters, 79% of the [Nd]d observations are modeled within ±10 pmol kg -1 of the observational value.The global Nd inventory (Fig. 3, dash contours) is 4.6×10 12 g, which is close to the current estimate of 4.2×10 12 g (Arsouze et al., 2009;Tachikawa et al., 2003).In this set-up, model εNd is also similar to the observations (Fig. 4b) and 76% of the points are within ±3 εNd unit of the observational value.However, the interbasin gradient of εNd is slightly smaller than observation (Fig. 4b).The Nd residence time (τNd) in this setup is 1150 yr.
We then examine the dependence of JεNd on the two parameters in our model (Fig. 3, color shading).JεNd is more dominated by the change of than fboundary,  et al., 2009;Tachikawa et al., 2003).
As pointed out by Rempfer et al., (2011), it is not possible to have a setup of fboundary and (EXP1) or decrease by 1×10 -4 (EXP2) from CTRL, we get similar results as CTRL (Table 2).The detailed results of CTRL are shown in the next section.
Our show good agreement with observations (Fig. 5a, with the track indicated by the black line in Fig. 2a), the same as in Rempfer et al., (2011).In addition to the good agreement with observations of [Nd]d at the seafloor (Fig. 9), the simulated [Nd]d agrees with observations at other depths as well (Fig. 5a).
[Nd]d increases with depth at all latitudes, as in the observations, due to the reversible scavenging by particles transport Nd isotopes downward.The vertical gradient is small in the high latitude North Atlantic because of the deep convection there.In addition, the than observation at shallow depth, as pointed out previously.This phenomenon is especially obvious in the region downstream of the Sahara desert plume, probably due to the overestimation of dust scavenging ability in the model (Fig. 5a).All of these features are consistent with Rempfer et al., (2011).
Cross-sections of εNd also agree well with observations (Fig. 5b).CTRL can simulate the εNd signatures of different water masses: NADW is -12 to -13; AAIW and AABW are -7 to -9 and the North Pacific is larger than -4.The "zig-zag" pattern of εNd in the Atlantic is also well simulated (Goldstein and Hemming, 2003) (Fig. 10 profile 5 and 7), showing different water masses dominate in different depth ranges, alternating from AAIW to NADW to AABW from the surface to the bottom.In the Atlantic, the simulated εNd covaries with salinity (Fig. 11), indicating northward and southward flowing water masses.However, for the same water mass, the εNd tongue seems to be shifted to deeper depth compared with salinity, which is attributed to the effect of scavenging by sinking particles (van de Flierdt et al., 2016;Rempfer et al., 2011).In general, our simulated εNd captures the main features of εNd in observation.However, compared with observation, the surface Pacific εNd in CTRL is lower (Fig. 5b and Fig. 7c), which has also been observed in Rempfer et al., (2011).
Furthermore, the vertical profile of εNd in the North Atlantic (Fig. 10 profile 4) in CTRL is quite different from observations.The simulation shows a more uniform εNd profile but observation shows large vertical gradients, which indicates that these regional model-data discrepancies may be caused by simplification of the Nd sources as well as the coarse resolution of the model (Rempfer et al., 2011).
Overall, our CTRL shows a simulation of both Nd concentration and εNd consistent with observations.The fbounday and in CTRL are near the values used in Rempfer et al., (2011).If we use the exact fbounday and 2), which suggests that the parameters tuned in Rempfer et al., (2011) can also be used in other models.Also, CTRL here is tuned based on observations from van de Flierdt et al., (2016).We have also tried using the observational data from Lacan et al., (2012) to tune the parameters as in Rempfer et al., (2011) and the results of the optimal parameter combination (CTRL_old) are shown in Table 2, which show similar performance as in Rempfer et al. (2011) (CTRL_R*), as well as CTRL, suggest the robustness of the overall range of the two parameters (fbounday and ).
We also apply the tuned parameters of fbounday and in abio_Nd to eco_Nd and the fields of both [Nd]d and εNd in eco_Nd are similar to the abio_Nd.The magnitude of the difference of [Nd]d and εNd between eco_Nd and abio_Nd is quite small compare the magnitude of [Nd]d and εNd in CTRL (Fig. 12).Therefore, parameters tuned in abio_Nd also apply to eco_Nd.

Sensitivity to fbounday and
[] [] In this section, we show model sensitivity to fbounday and by studying cases with double or half fbounday and from CTRL.BS05 is the half fbounday case and BS20 is the double fbounday case.PD05 is the half case (note: we don't have the exact half case) and PD20 is the double 2).Here, BS05, BS20, and PD20 have reached their equilibrium while PD05 has not.Although PD05 has not reached equilibrium, it does not matter for the purpose of showing model sensitivity, which will be discussed below.
A change of fbounday mainly effects Nd concentration, with only modest impact on εNd, especially in the deep ocean.In comparison with CTRL, the total number of points within ±10 pmol kg -1 of observational [Nd]d (J1) is much reduced in BS05 and BS20 (Table 2).The [Nd]d increase with depth and circulation pathway is still (BS20)(Table 2).If the boundary source is the only Nd source, double (half) fbounday will lead to a double (half) of [Nd]d everywhere.However, the difference of εNd between BS05 or BS20 and CTRL is small (Fig. 14).The number of points within ±3 εNd unit of observational εNd (J2) is comparable to CTRL in both BS05 and BS20.The increases (decreases) with the increase (decrease) of , which is consistent with the shorter (longer) residence time in PD20 (PD05) (Table 2).In PD20, the North Atlantic becomes more unradiogenic while the North Pacific becomes more radiogenic compared to CTRL (opposite in PD05), consistent with Rempfer et al., (2011).With a larger , more efficient scavenging by sinking particles transports more Nd from the water column to the ocean floor locally, where it is buried as sediment.This leaves less Nd that carries the εNd signature from its source region to be transported by the circulation, resulting in a larger interbasin εNd gradient.In addition, the slightly deeper εNd tongue than salinity for the same water mass mentioned previously appears more (less) prominent in PD20 (PD05) as AAIW and NADW are all shifted to slightly deeper (shallower) depth than CTRL in PD20 (PD05), another indication that this mismatch is due to the reversible scavenging (Rempfer et al., 2011).It should be noted that PD05 is not completely equilibrated yet as stated previously and if run until complete equilibrium, the difference described above will be even larger.
. Model Dev.Discuss., doi:10.5194/gmd-2017-40,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 13 March 2017 c Author(s) 2017.CC-BY 3.0 License.To evaluate the performance of each parameter combination, we compare the cost function of [Nd]d and εNd as in Rempfer et al., (2011each observational point, N is the total number of observational points, obsk is the observational [Nd]d or εNd and modelk is the model [Nd]d or εNd at the observational location.The cost function J measures the average deviation of the simulated [Nd]d or εNd from observation.In addition to the cost function, we also examine the distributions of [Nd]d and εNd compared with observations.
in fboundary, J[Nd]d also first decreases and then increases.As such, J[Nd]d Dev. Discuss., doi:10.5194/gmd-2017-40,2017   Manuscript under review for journal Geosci.Model Dev. Discussion started: 13 March 2017 c Author(s) 2017.CC-BY 3.0 License.again similar to inRempfer et al., (2011).JεNd reaches the minimum at fboundary = 1×10 9 g yr -1 and[!"] ![!"] != 16×10-4 , which is indicated by the red triangle in Fig.3.In this setup of optimal JεNd, εNd in the model shows a good agreement with the observation (Fig.4d, Fig.6e~h), with 90% of the points simulated within ±3 εNd unit of the observational value.The residence time τNd is now 351 yr, much shorter than that in the optimal J[Nd]d.Therefore, the interbasin gradient of εNd is more prominent than that in the optimal J[Nd]d (Fig.4 b and d).However,[Nd]d in this setup is much smaller than the observations (Fig.4c) and the model only simulates 19% of the points within ±10 pmol kg -1 of the observational Nd concentration.The Nd inventory is only 0.9×10 12 g, much smaller than the estimate of 4.2×10 12 g (Arsouze both J[Nd]d and JεNd simultaneously.In our model, when J[Nd]d reaches the minimum, the overall performance of εNd is not good enough and when JεNd reaches the minimum, [Nd]d is too far away from the observations.Therefore, we further examined the distributions of [Nd]d and εNd as well as the cost functions in other parameter settings as we did for the two cases with minimum J[Nd]d or JεNd.We first pick the combinations that can produce both reasonable [Nd]d and εNd distributions.Then we pick the one simulation that have relative better performance in simulating both [Nd]d and εNd among those combinations, as our control experiment (CTRL).The simulation is not very sensitive to the parameter combination around CTRL, for example, if we increase fboundary by 1×10 9 g yr -1 Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-40,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 13 March 2017 c Author(s) 2017.CC-BY 3.0 License.CTRL can simulate 72% of the [Nd]d within ±10 pmol kg -1 of the observational value and 82% of the εNd within ±3 εNd unit of the observational value.The performance of simulated [Nd]d in CTRL is comparable with the optimal J[Nd]d case but with a much CTRL can simulate the distributions of [Nd]d and εNd in good agreement with the observations.The global and each basin distributions of [Nd]d and εNd in CTRL are overall similar to the observations (Fig. 5-7).The correlation between model [Nd]d and observation is 0.50 and the correlation between model εNd and observation is 0.83, both are significant at 0.01 significance level.However, [Nd]d in CTRL is smaller than observations at shallow depth and the largest and smallest εNd found in observations are not simulated by CTRL (Fig. 6 and 8), similar to Rempfer et al. (2011).Seafloor maps of [Nd]d and εNd in CTRL are consistent with observations (Fig.9 and 10).[Nd]d shows a general increase along the circulation pathway from the North Atlantic to the North Pacific, in agreement with observations, except for some localized points in the North Atlantic (Fig. 9).[Nd]d at the seafloor is near 20 pmol/kg in the North Atlantic, around 30 pmol/kg in the Southern Ocean and almost 40 pmol/kg in the North Pacific.The seafloor map of εNd in CTRL shows an interbasin gradient as in the observation (Fig. 10).The lowest εNd values occur in the North Atlantic (<-14); the highest εNd values occur in the North Pacific (~0), with the intermediate εNd values in the Indian and Southern Ocean (-7 to -10), similar to observations.Cross-sections of [Nd]d from the North Atlantic to the North Pacific in CTRL Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-40,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 13 March 2017 c Author(s) 2017.CC-BY 3.0 License.increase along the circulation pathway is also clear in the cross-section of [Nd]d.[Nd]d increases from the surface North Atlantic to the deep Southern Ocean to the deep Pacific.A comparison of the simulation and observations in selected vertical profiles of [Nd]d also show good agreement (Fig. 9), although model [Nd]d is lower Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-40,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 13 March 2017 c Author(s) 2017.CC-BY 3.0 License.reproduced in both BS05 and BS20, but the vertical sections of [Nd]d along the track from the North Atlantic to the North Pacific show much decreased (increased) [Nd]d in BS05 (BS20) (Fig. 13), consistent with smaller (larger) Nd inventory in BS05 pattern and magnitude of the difference between BS05 (BS20) and CTRL are similar to that inRempfer et al. (2011).The change of εNd in BS05 and BS20 is small because the εNd in Nd sources are fixed at prescribed values and there is no differentiation between 143 Nd and 144 Nd during the reversible scavenging.Changing fbounday only changes the relative contributions of the boundary source and the surface sources (dust and river), whose influences are limited to the upper 1km(Rempfer et al., 2011).If the boundary source is the only Nd source, changing fbounday will not change the εNd distribution at all, since εNd measures the ratio between 143 Nd and 144 Nd.Nd]d and εNd significantly (Fig.15 and 16).J[Nd]d in PD05 and PD20 are much larger than that in CTRL and the percentages of [Nd]d within ±10 pmol kg -1 of observational [Nd]d are greatly reduced (Table 2).Increasing [!"] ![!"] !decreases [Nd]d as it increases the compensation of [Nd]p and, in turn, sedimentation, and vice versa (Fig. 15).Increased [!"] ![!"] ! will also reduce the vertical gradient of [Nd]d and the gradient along the circulation pathway such that [Nd]d becomes more homogeneous.This is because a larger [!"] ![!"] !favors the vertical transport of Nd by particles and therefore Nd is less subject to the transport by circulation.Hence, the Nd transport by circulation is overwhelmed by vertical scavenging.In addition, in the region under the Saharan dust plume, the low [Nd]d signature in CTRL is weakened in PD05 but intensified in PD20, and the low [Nd]d signature reaches almost the ocean bottom in PD20.The εNd also shows a large difference between PD05 or PD20 and CTRL (Fig. 16).The interbasin gradient of εNd Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-40,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 13 March 2017 c Author(s) 2017.CC-BY 3.0 License.
both a fully biotic (eco_Nd) and an approximate abiotic (abio_Nd) version.Extensive sensitivity experiments are performed to test the model simulation with respect to fboundary and [!"] ![!"] !, which are the tuning parameters in the Nd model.With the parameters we found under present climate forcing for our CTRL, our model is able to simulate the major features of the present-day global distribution of Nd and εNd, with reasonable agreement with the observation.However, the model [Nd]d is smaller than observations at shallow depth and the model does not simulate the extremes of εNd that are seen in the observations.These biases are similar to the biases in the Nd implementation of Rempfer et al. (2011) in the Bern3D model.Despite these shortcomings, the simulated εNd is a useful tracer for paleoclimate studies as simulated εNd in CTRL captures the different εNd signatures for different Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-40,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 13 March 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 3 .
Figure 3. Contours of cost function of [Nd]d (J[Nd]d: solid lines, unit: pmol kg -1 ) and εNd (JεNd: color shading, unit: εNd units) and total Nd inventory (dashed lines, unit: 10 12 g) for different combinations of fboundary and [!"] ![!"] ! at the end of each experiment.J[Nd]d reaching minimum is indicated by the pink square; JεNd reaching minimum is indicated by the red triangle; CTRL is indicated by the yellow star.

Figure 4 .
Figure 4. Vertical sections of [Nd]d (left) and εNd (right) along a track from the North Atlantic to the North Pacific when J[Nd]d reaches minimum (a and b) and JεNd reaches minimum (c and d).Color contours are model results and observations are attached as filled cycles using the same color map.

Figure 5 .
Figure 5. Vertical sections of [Nd]d (a) and εNd (b) along the track (indicated in Fig.2 (a)) from the North Atlantic to the North Pacific in CTRL.Color contours are model results and observations are attached as filled cycles using the same color map.

820 Figure 6 .Figure 7 .
Figure 6.Histograms of observational (white) [Nd]d and model values at 821 observation location in CTRL (grey).First column shows the distribution of the 822 global ocean, second column for the Atlantic, third column for the Pacific and the 823 forth column for the Indian and Southern Ocean.Each ocean basin has been 824 separated into different depth ranges: first row for 0-126m, second row for 126m-825 1039m, third row for 1039m-2755m and the forth row for ocean deeper than 826 2755m.827 828 829

Figure 12 .
Figure 12.Difference of [Nd]d and εNd between eco_Nd and abio_Nd.(a) [Nd]d 863 difference at the seafloor.(b) [Nd]d difference along the track (indicated in Fig.2 864 (a)).(c) εNd difference at the seafloor.(d) εNd difference along the track.865 866 Figure 15.Vertical sections of [Nd]d along the track from the North Atlantic to the North Pacific (indicated in Fig.2 (a)) for the sensitivity on [!"] ![!"] !: (a) PD05 (b) difference between PD05 and CTRL.(c) PD20 (d) difference between PD20 and CTRL.Observations are attached as filled cycles using the same color map in (a) and (c).
Figure 16.Vertical sections of εNd along the track from the North Atlantic to the North Pacific (indicated in Fig.2 (a)) for the sensitivity on [!"] ![!"] !: (a) PD05 (b) difference between PD05 and CTRL.(c) PD20 (d) difference between PD20 and CTRL.Observations are attached as filled cycles using the same color map in (a) and (c).
. Nd has three sources: river source, dust source, and boundary source.Sedimentation of Nd is the only sink for Nd budget.143Ndand 144re modeled as two separate tracers.In addition to 143 Nd and 144 Nd, Nd also has other stable isotopes and the sum of143Nd and 144 Nd accounts for 36% of total Nd(Magill et al., 2006).Since we use 36% of the total Nd fluxes as fluxes for 143 Nd and 144 Nd, we need to scale the simulated

Table 2 .
Rempfer et al. 2011 uses 12 g/mol for CaCO3 (g-C) Parameters and general performance for different experiments.CTRL is the parameters tuned for our model.J[Nd]d min and JεNd min are cost function of [Nd]d Rempfer et al., (2011)combination if we use the observational data fromLacan et    al., (2012)as inRempfer et al., (2011).J1 is the percentage of the observation, which model [Nd]d is within ±10 pmol kg -1 of observation.J2 is the percentage of the observation, which model εNd is within ±3 εNd unit of observation.