1 Introduction

Seawater contains various trace elements, some of which are known to play an important role in controlling the ocean biogeochemical cycles. However, because strict clean sampling method and highly sensitive analytical method are required to determine trace element concentrations, the observational coverage of such trace elements remained very poor until recently. The international GEOTRACES program is now correcting the poor data coverage, providing us with a valuable global dataset about trace elements (Anderson 2020).

Neodymium (Nd) isotope has been recognized as a useful chemical tracer for investigating the modern and past ocean circulation (Lacan and Jeandel 2004a; Piotrowski et al. 2004); it was therefore selected as one of key parameters in the GETRACES program (van de Flierdt et al. 2012). Nd belongs to the group of rare earth elements (REEs). REEs have similar chemical characteristics; however, they also exhibit systematic differences because their ionic radii decrease with an increase in atomic number (Hathorne et al. 2015). Consequently, a small but systematic variation in particulate affinity exists among REEs (Sholkovitz et al. 1994; Nozaki 2001; Hara et al. 2009). An investigation of the distribution of REEs in the ocean and their differences among REEs can provide an opportunity to gain deeper and more comprehensive understanding about what kinds of source and sink are important for controlling the cycles of trace elements.

In the periodic table reported by Nozaki (2001), the vertical profiles of dissolved REEs measured in the North Pacific Ocean are similar for different REEs in that their concentrations increase with depth for all REEs (except Ce and Pm). However, they also exhibit systematic differences among REEs: light REEs (LREEs) increase almost linearly with depth, whereas heavy REEs (HREEs) exhibit more nutrient-like profiles. A previous study (Oka et al. 2009) investigated the cause of the systematic difference in the vertical profiles among REEs found in the North Pacific Ocean. Using ocean general circulation model (OGCM) simulations under various types of sources and sinks in the ocean, Oka et al. (2009) demonstrated that including the reversible scavenging process can explain the observed vertical profiles of both LREEs and HREEs, whereas other processes (e.g., biological uptake) have difficulty in explaining the profiles of both LREEs and HREEs. The reversible scavenging has been recognized as an important process for the vertical transport of thorium (Nozaki et al. 1981, 1987; Bacon and Anderson 1982), where a simple vertical one-dimensional balance between the radioactive source and removal by the scavenging can successfully explain the linearly increasing vertical profile of thorium. For REEs, Oka et al. (2009) demonstrated that the balance between three-dimensional physical transport and the reversible scavenging controls the vertical profiles of REEs.

Knowledge of the external source of REEs is also important in terms of utilizing the Nd isotope as a proxy for ocean circulation because the amount of external Nd input to the ocean and its isotopic composition controls the Nd isotope distribution in the ocean. From previous studies on Nd and its isotope ratio, continental shelf areas were identified as the most important suppliers of Nd to the ocean (Tachikawa et al. 1997; Lacan and Jeandel 2005; Jeandel et al. 2007). By considering source from the coastal regions, several studies succeeded in reproducing the global distribution of Nd and its isotope ratio in the ocean (Arsouze et al. 2007, 2008; Jones et al. 2008; Rempfer et al. 2012). However, although global modeling of Nd and its isotope ratio has been performed in several studies as described above, there has been no global simulation of multiple REE abundances except that performed in Oka et al. (2009). Oka et al. (2009) demonstrated that the reversible scavenging is an important process for explaining the distribution of REEs and the systematic differences in their vertical profiles in the North Pacific Ocean; however, the reproducibility in the other areas and the importance of the external source were only briefly discussed therein and have not yet been thoroughly investigated. One of reason for this was due to insufficient coverage of observational data on REEs; however, as mentioned, the international GEOTRACES program is now addressing such poor data coverage. Therefore, it is timely to attempt a global simulation of multiple REEs and a comprehensive comparison of the global distributions of REEs between the model and observations.

In this study, toward more realistic simulation of REEs in the global ocean than that reported in Oka et al. (2009), we focus on the external source of REEs and discuss its role in controlling the global distribution of REEs in the ocean. As mentioned earlier, boundary input from coastal regions has been proposed to be the most important external source of Nd; therefore, in this study, we extend our modeling by explicitly considering this source. We will compare the results of our model with globally covered observations, including recently obtained dataset from the GEOTRACES program. Through REE simulations and additional sensitivity experiments, we discuss the role of the internal cycle (i.e., the reversible scavenging) and the external source (i.e., REE input from the coastal areas) in controlling the global distribution of REEs and the differences in their distributions in the ocean.

The remainder of this paper is organized as follows. The model and observational data used in this study are described in Sect. 2. The design of the model simulations is explained in Sect. 3. The results and discussion are presented in Sects. 4 and 5, respectively. Finally, the summary and concluding remarks are provided in Sect. 6.

2 Method

2.1 Observational data of REEs

The observational data of REEs referenced in this study are presented in Fig. 1. The REE observations taken from GETRACES Intermediate Data Product 2017 (Schlitzer et al. 2018) are plotted as red circles in Fig. 1. Of the GEOTRACES IDP2017 data, we utilized REE data from the mid-latitude sections of the North Atlantic (GA03; Stichel et al. 2015), the tropical South Atlantic (GAc01; Zheng et al. 2016), the western South Atlantic (GIPY04; Garcia-Solsona et al. 2014), the Atlantic sectors of the Southern Ocean (GIPY05; Stichel et al. 2012a, b), the central North Pacific (GPpr05; Fröllje et al. 2016), and the western North Pacific (GPpr04; Behrens et al. 2018a, b). For Nd, the data were also included from the eastern North Atlantic (GA02; Lambelet et al. 2016), the eastern tropical North Atlantic (GA11; Zieringer et al. 2019), the Pacific sectors of the Southern Ocean (GIPY06; Lambelet et al. 2018), and the polar South Pacific (GPc02; Basak et al. 2015). We also used pre-GEOTRACES observations plotted as black open circles in Fig. 1, which are compiled from data reported in pre-GEOTRACES literature (Elderfield and Greaves 1982; De Baar et al. 1983, 1985; Sholkovitz and Schneider 1991; Piepgras and Jacobsen 1992; Bertram and Elderfield 1993; Sholkovitz et al. 1994; Zhang et al. 1994, 2008; Shimizu et al. 1994; German et al. 1995; Sholkovitz 1996; Zhang and Nozaki 1996, 1998; Tachikawa et al. 1999b; Alibo and Nozaki 2000; Nozaki and Alibo 2003; Lacan and Jeandel 2004b; Hongo 2005).

Fig. 1
figure 1

Location of reported observations about REEs. GEOTRACES data are shown with red circles and open circles represent pre-GEOTRACES observations. Sites used in Fig. 4 are indicated by the characters “P”, “S”, and “A” (color figure online)

Full size image

2.2 Physical model

The ocean general circulation model (OGCM) used in this study is COCO (Hasumi 2006). The model has 120 and 128 grids in the horizontal directions and 40 vertical layers. The model setup is identical to TideNF used by Oka and Niwa (2013) except that the background vertical diffusivity (Kb) is set to be 0.01 [cm2/s] instead of 0.1 [cm2/s] (Komuro 2014). In this model, the tidal mixing parameterization is applied for vertical diffusivity. This parameterization was shown to improve the reproducibility of ∆14C distribution (Oka and Niwa 2013), which means that our physical model can realistically simulate water mass age in the ocean.

2.3 Tracer model

The concentrations of ocean tracers (i.e., REEs in this study) are calculated using the following tracer equation:

$$\frac{\partial C}{{\partial t}} = - \nu \cdot \nabla C + K_{{\text{H}}} \nabla_{{\text{H}}}^{2} C + K_{{\text{V}}} \frac{{\partial^{2} C}}{{\partial z^{2} }} + Sc,$$
(1)

where C is the concentration of the tracer under consideration, ν is the velocity, KH and KV are the horizontal and vertical diffusivities, respectively, and Sc represents source/sink term of the tracer due to external flux and/or internal cycles. Note that although isopycnal diffusion and layer thickness diffusion are not explicitly expressed in Eq. (1), they are actually calculated in the model. The calculation of Eq. (1) is performed by using an offline method (Oka et al. 2008), and the physical fields (i.e., ν, KH, KV, and isopycnal/layer-thickness diffusivities) are prescribed from the result of the above-mentioned OGCM simulation. The treatment of the source/sink term is described in Sect. 2.4.

2.4 Source/sink term

For the source/sink term (Sc) of REEs, we consider the reversible scavenging process (\(S{c}_{\mathrm{scav}}\)) and external input (\(S{c}_{\mathrm{in}}\)):

$$Sc=S{c}_{\mathrm{scav}}+ S{c}_{\mathrm{in}}.$$
(2)

The first term represents the internal cycle, which does not affect the globally averaged concentration of REEs. The second term represents the external process, which does affect the globally averaged concentration. For the internal cycle of REEs, we previously demonstrated that the vertical profile of REEs and their systematic differences among REEs in the North Pacific Ocean can be explained clearly by the reversible scavenging process (Oka et al. 2009). Therefore, we take this process into account as the first term in the present study. For the external input of REEs, although our previous approach (Oka et al. 2009) treated this process in a very simple manner, we use a more realistic approach in the present study. The detailed treatment is described below.

Following the previous studies (Siddall et al. 2005; Oka et al. 2009), the reversible scavenging process is parameterized as follows:

$$S{c}_{\mathrm{scav}} = {w}_{\mathrm{s}}\frac{\partial {A}_{\mathrm{p}}}{\partial z},$$
(3)

where \({A}_{\mathrm{p}}\) is the concentration of particulate-phase REEs and \({w}_{\mathrm{s}}\) is the constant sinking velocity (1000 m/yr) of particles. The particulate-phase REEs (\({A}_{\mathrm{p}}\)) are assumed to be diagnosed from the following partitioning relationship:

$${K}_{\mathrm{P}} = \frac{{A}_{\mathrm{p}}}{{A}_{\mathrm{d}}{C}_{\mathrm{P}}},$$
(4)

where \({K}_{\mathrm{P}}\) is a dimensionless partitioning coefficient, \({A}_{\mathrm{d}}\) is the concentration of dissolved-phase REEs, and \({C}_{\mathrm{P}}\) is a dimensionless ratio of the particle mass density to seawater density. From Eq. (4) and the fact that the concentration of REEs (\(C\)) in Eq. (1) is defined as the sum of Ap and Ad,

$$C={A}_{\mathrm{p}}+{A}_{\mathrm{d}},$$
(5)

where the concentration of particulate-phase REEs (\({A}_{\mathrm{p}}\)) can be diagnosed as:

$${A}_{\mathrm{p}}=\frac{C}{{1+K}_{\mathrm{P}}{C}_{\mathrm{P}}}.$$
(6)

The partitioning coefficient (\({K}_{\mathrm{P}}\)) represents the affinity of REEs to particles and its appropriate value depends on individual REEs, which we discuss in Sect. 3.2. The distribution of the particle concentration (\({C}_{\mathrm{P}}\)) is estimated from the sum of biogenic particles (organic carbon, calcite, and opal) and lithogenic particles. The former is diagnosed based on the satellite observations (Behrenfeld and Falkowski 1997; Dunne et al. 2005), while the latter is obtained from the results of a dust-transport model simulation (Takemura et al. 2005). The obtained distribution of \({C}_{\mathrm{P}}\) is similar to those reported in (Siddall et al. 2005) and (Oka et al. 2009), and the detailed procedure about estimation of \({C}_{\mathrm{P}}\) is provided in these previous studies. In summary, after specifying \({K}_{\mathrm{P}}\), Eqs. (3) and (6) are used to calculate the reversible scavenging process. For Eq. (3), the boundary conditions at the sea surface and bottom needs to be specified. In this study, the zero flux is assumed at the sea surface and all particulate REEs reaching the sea bottom are assumed to be removed from the ocean:

$$w_{s} A_{p} = \left\{ {\begin{array}{*{20}l} {0,} & {{\text{at}}\;{\text{the}}\;{\text{sea}}\;{\text{surface}}} \\ {w_{s} A_{{p\left( {{\text{bottom}}} \right)}} ,} & {{\text{at}}\;{\text{the}}\;{\text{sea}}\;{\text{bottom}}} \\ \end{array} } \right.,$$
(7)

where \({A}_{\mathrm{p}(\mathrm{bottom})}\) is the particulate REE concentration at the sea bottom.

For the external input (\(S{c}_{\mathrm{in}}\)), we consider the boundary flux of REEs around the coastal regions because the boundary flux around these regions has been recognized as the most important source area of REEs into the ocean (Jeandel et al. 2007). In this study, we specify this source in the form of surface REE flux as follows:

$$Sc_{{{\text{in}}}} = \left\{ {\begin{array}{*{20}c} {\frac{{F_{{{\text{REE}}}} }}{{\Delta z_{1} }},} & {k = 1} \\ {0, } & { k > 1} \\ \end{array} } \right.,$$
(8)

where \({F}_{\mathrm{REE}}\) is the two-dimensional REE input flux per unit area, \(\Delta {z}_{1}\) is the thickness of the model’s surface (top) layer, and k is the index of the model’s vertical layers; k = 1 corresponds to the model’s top layer, while k > 1 represents layers below the top layer. To obtain the distribution of \({F}_{\mathrm{REE}}\), we first estimate the Nd flux by following previous modeling approaches (Arsouze et al. 2007), which are described in Sect. 3.

The processes related to sources and sinks of REEs considered in this study are summarized in Fig. 2a. The term \(S{c}_{\mathrm{scav}}\) in Eq. (2) represents for reversible scavenging, as previously mentioned. This term also implicitly includes the particle dissolution process because \({A}_{\mathrm{p}}\) is not assumed to be constant but has vertical profiles, which means that particle dissolution occurs according to Eq. (3). Further, the sedimentation process is represented again by the term \(S{c}_{\mathrm{scav}}\) in the form of its boundary condition at the sea bottom (see Eq. (7)). Finally, the term \(S{c}_{\mathrm{in}}\) represents for the source from the continental shelf.

Fig. 2
figure 2

Processes related to REEs cycles in the ocean. a Processes considered in our CTL experiment. b Processes considered in our RCYC experiment

Full size image

3 Experimental design

3.1 Preliminary Nd simulation

Before REE simulations, we first performed a preliminary Nd simulation to obtain the boundary REE flux (i.e., the distribution of \({F}_{\mathrm{REE}}\) in Eq. (8)). In this preliminary simulation, the boundary \({F}_{\mathrm{REE}}\) of Nd (\({F}_{\mathrm{Nd}}\)) was determined from

$$F_{{{\text{Nd}}}} = \left\{ {\begin{array}{*{20}l} {r\left( {C_{{{\text{sfc}}}} - C_{{{\text{obs}}}} } \right),} & {{\text{where}}\; C_{{{\text{obs}}}} \;{\text{is}}\;{\text{available}}} \\ {0,} & {{\text{elsewhere}}} \\ \end{array} } \right.,$$
(9)

where \({C}_{\mathrm{sfc}}\) is the simulated Nd concentration at the model’s top layer, \({C}_{\mathrm{obs}}\) is the observed Nd concentration near the coastal regions, and r is the restoring time constant (set to 30 day−1). This approach is similar to that in a previous modeling study by Arsouze et al. (2007); however, our treatment is simpler in that we assume a constant value for \(r\) and the Nd flux to be located at the surface rather than from the sea bottom. The distribution of \({C}_{\mathrm{obs}}\) is shown in Fig. 3a, which includes compiled data provided by Jeandel et al. (2007). For conversion into the model grid, we first converted individual point data from the dataset (i.e., longitude/latitude/concentration database) into the model ocean grid closest to the data point. Next, to avoid a single-grid point source, the same concentration was also specified in neighboring model grids (i.e., in total, individual single-point observational data were converted into nine-point model grids). For a model grid containing multiple single-point data in Jeandel et al. (2007), their average concentration was used. From this procedure, \({C}_{\mathrm{obs}}\) defined in the model grid was obtained, and the source flux of Nd was calculated based on Eq. (9) (see Fig. 3b).

The value of \({K}_{\mathrm{P}}\) should be dependent on individual REEs; after performing simulations for several values of \({K}_{\mathrm{P}}\), we selected 8.0 × 105 as the most appropriate value of \({K}_{\mathrm{P}}\) for Nd simulation. The value of \({K}_{\mathrm{P}}\) controls the particle fraction of Nd (\({\sim A}_{\mathrm{p}}/{A}_{\mathrm{d}}\)), as expressed in Eq. (4). Our choice (i.e., \({K}_{\mathrm{P}}\) = 8.0 × 105) is equivalent to assuming that the particle fraction of Nd averaged above the thermocline (\(\sim 1000 \mathrm{m})\) becomes around 3% under our particulate (i.e., \({C}_{\mathrm{P}}\)) distributions. This is consistent with the observational estimate reported by Nozaki (2001): the globally averaged particle fraction of Nd was estimated to be 3.31% for Nd.

Starting from a constant zero concentration, the model reached a steady state after 3000-year integration, and the last 100 years were averaged and used for analysis. The distribution of \({F}_{\mathrm{Nd}}\) averaged over the last 100 years in the preliminary Nd simulation is presented in Fig. 3b; this flux was utilized for subsequent REE simulations.

Fig. 3
figure 3

a Nd concentration at the sea bottom compiled by Jeandel et al. (2007). Unit is pmol/kg. b Nd boundary flux obtained from our preliminary Nd simulation. Unit is pmol/m2/s (color figure online)

Full size image

3.2 REE simulation (CTL experiment)

In this study, we focused on three REEs as representative light REE (LREE), medium REE (MREE) and heavy REE (HERR), i.e., neodymium (Nd), dysprosium (Dy), and lutetium (Lu), respectively. After performing the preliminary Nd simulation, the obtained Nd flux (i.e. \({F}_{\mathrm{Nd}})\) was used for REE simulations. In the REE simulations, \({F}_{\mathrm{REE}}\) was assumed to be proportional to \({F}_{\mathrm{Nd}}\) depending on individual REEs:

$${F}_{\mathrm{REE}}={f}_{\mathrm{REE}}{F}_{\mathrm{Nd}},$$
(10)

where \({f}_{\mathrm{REE}}\) is a constant ratio of REE flux to Nd flux. The values of \({f}_{\mathrm{REE}}\) and \({K}_{\mathrm{P}}\) should be dependent on individual REEs, which is discussed in detail below.

To select the appropriate values of \({f}_{\mathrm{REE}}\) for each REE, the globally averaged concentration and mean residence time estimated by Nozaki (2001) were utilized. The total inflow of REE (i.e., globally integrated \({F}_{\mathrm{REE}}\)) should be balanced with the total outflow of REE. The total outflow is proportional to the residence time (\({t}_{\mathrm{REE}}\)) and mean concentration (\(\overline{C} _{{{\text{REE}}}}\): \({F}_{\mathrm{REE}}\propto {t}_{\mathrm{REE}}\overline{C} _{{{\text{REE}}}}\)). Therefore, \({f}_{\mathrm{REE}}\) can be roughly estimated from

$$f_{{{\text{REE}}}} = \frac{{F_{{{\text{REE}}}} }}{{F_{{{\text{Nd}}}} }} = \frac{{t_{{{\text{REE}}}} \overline{C} _{{{\text{REE}}}} }}{{t_{{{\text{Nd}}}} \overline{C} _{{{\text{Nd}}}} }},$$
(11)

where \({t}_{\mathrm{Nd}}\) is the residence time of Nd and \(\overline{C} _{{{\text{Nd}}}}\) is the globally averaged concentration of Nd. From the observational data, the mean residence time was estimated to be 400 years for Nd, 740 years for Dy, and 2890 years for Lu, whereas the globally averaged concentration is 20 pmol/kg for Nd, 6.5 pmol/kg for Dy, and 1.3 pmol/kg for Lu (Nozaki 2001). From these values and Eq. (11), the values of \({f}_{\mathrm{REE}}\) were estimated to be 0.176 and 0.009 for Dy and Lu, respectively (and, as a matter of course, the value of \({f}_{\mathrm{REE}}\) was 1.0 for Nd).

The value of \({K}_{\mathrm{P}}\) must also be specified depending on each REE. The appropriate value of \({K}_{\mathrm{P}}\) can be obtained from observational information of dissolved-phase and particulate-phase REE concentrations. From Eq. (4), the ratio of \({K}_{\mathrm{P}}\) between Nd and REE becomes

$$\frac{{K}_{\mathrm{P}}}{{K}_{\mathrm{P}}(Nd)} = \frac{{A}_{\mathrm{p}}/{A}_{\mathrm{d}}}{{A}_{\mathrm{P}}(Nd)/{A}_{\mathrm{d}}(Nd)},$$
(12)

where the values with Nd in parentheses indicate those for Nd. In Table 2 of Nozaki (2001), global averages of the particle fraction of REEs (\({\sim A}_{\mathrm{p}}/{A}_{\mathrm{d}}\)) were reported to be 3.31% for Nd, 1.80% for Dy, and 0.46% for Lu. From these values and the fact that the value of 8.0 × 105 was selected for \({K}_{\mathrm{P}}(\mathrm{Nd})\), the values of \({K}_{\mathrm{P}}\) were estimated to be 4.4 × 105 for Dy and 1.1 × 105 for Lu.

In summary, by specifying different values of \({f}_{\mathrm{REE}}\) and \({K}_{\mathrm{P}}\) depending on individual REEs (summarized in Table 1), we performed Nd, Dy, and Lu simulations. Starting from a constant zero concentration, the model reached a steady state after a 9000-year time integration, and the last 100 years were averaged and used for analysis. Hereafter, we refer to the above-mentioned REE simulations as CTL experiment (CTL is an abbreviated word for “control”) to distinguish from RCYC experiment (RCYC is an abbreviated word for “recycle”) which is described in Sect. 3.3.

Table 1 Model parameter values which are specified depending on individual REEs
Full size table

3.3 REE simulation (RCYC experiment)

To investigate the role of the external source and sink (i.e., processes that affect the globally averaged concentration of REEs) in controlling the distribution of REEs, we performed an additional RCYC experiment. In the RCYC experiment, the globally averaged concentration of REEs was assumed to be conservative. For this purpose, both the term \(S{c}_{\mathrm{in}}\) in Eq. (2) and the boundary REE flux from the coastal regions represented by Eq. (7) were set to be zero. Namely, in the RCYC experiment, the Eqs. (7) and (8) used in the CTL experiment were replaced by

$$w_{{\text{s}}} A_{{\text{p}}} \left( {{\text{RCYC}}} \right) = \left\{ {\begin{array}{*{20}c} {0,} & {{\text{at }}\;{\text{the }}\;{\text{sea}}\;{\text{surface}}} \\ {0,} & {{\text{at}}\;{\text{the}}\;{\text{sea }}\;{\text{bottom}}} \\ \end{array} } \right.,$$
(13)

and

$$S{c}_{\mathrm{in}}(\mathrm{RCYC})=0,$$
(14)

respectively. The processes considered in the RCYC experiment are summarized in Fig. 2b. With the exception that there was no external boundary REE flux, the internal cycles of REEs associated with reversible scavenging and particle dissolution were treated in the same way as in the CTL experiment. The model parameter values (e.g., \({K}_{\mathrm{P}}\)) were also the same in both CTRL and RCYC experiments.

As in the CTL experiment, three REEs (i.e., Nd, Dy, and Lu) were simulated in the RCYC experiment. Starting from the final state of the CTL experiment, the model reached a steady state after a 5000-year time integration, and the last 100 years were averaged and used for analysis.

4 Results

4.1 Vertical profiles at Site P, Site A, and Site S

Figure 4 presents the vertical profiles of Nd, Dy, and Lu at the North Pacific site (Site P), the Southern Ocean site (Site S), and the North Atlantic site (Site A). At Site P, the observational data (dots in Fig. 4) indicate that the vertical profiles of REEs systematically vary from LREE to HREE; LREE (Nd) almost linearly increases with depth, whereas HREE (Lu) exhibits a more nutrient-like profile (Nozaki 2001). Our CTL experiment (bold lines in Fig. 4) successfully reproduced the shape of the vertical profiles of Nd, Dy, and Lu, and their systematic differences. We previously demonstrated that this systematic difference can be explained by differences in the values of \({K}_{\mathrm{P}}\) between REEs (Oka et al. 2009); lighter REEs adsorb onto particles more strongly than heavier REEs (i.e., larger \({K}_{\mathrm{P}}\) for lighter REEs), which leads to differences in the vertical profiles. It is worth noting that realistic treatment of the external source of REEs (i.e., the term \(S{c}_{\mathrm{in}}\)) in this study led to the improved vertical profile of Nd; the linearly increasing profile was more clearly reproduced in our CTL simulation than in results reported in Oka et al. (2009). Although the model tended to underestimate the concentration of LREE to some degrees, there was good agreement between the observation and the model, especially for HREE at the North Pacific site.

Fig. 4
figure 4

Vertical profiles of REEs: a Nd, b Dy, and c Lu at Site P (140E, 40N); d Nd, e Dy, and f Lu at Site S (140E, 54S); g Nd, h Dy, and i Lu at Site A (64W, 31N). Locations of Sites P, S, and A are shown in Fig. 1. Results from CTL and RCYC experiments are shown by bold and thin lines, respectively. Thin broken lines show preformed REE concentration in CTL experiment. Dots represent observations from Nozaki (2001) for Site P, Hongo (2005) for Site S, and Sholkovitz et al. (1994) for Site A

Full size image

At the North Atlantic site (Site A), the observational data indicate that the concentrations increase with depth as in Site P but exhibit more uniform vertical profiles than those at the Site P. It is also observed that the observational data exhibit a local maximum around 1500 m for all REEs. Our CTL experiment successfully reproduced the vertical profiles at Site A, including the local maxima around 1500 m, although the surface concentrations of Nd and Dy were somewhat underestimated and the local maxima of Lu were slightly overestimated in the model.

At the Southern Ocean site (Site S), the observational data suggest that the vertical profiles of Nd, Dy, and Lu generally resemble each other; the concentrations simply increase with depth for all the REEs. Quantitatively speaking, this increase with depth is more significant for lighter REEs, as expected from the fact that the vertical transport associated with reversible scavenging becomes stronger for lighter REEs. Although the model somewhat overestimated this tendency (e.g., the model underestimated the surface concentrations of Nd and Dy while also underestimating the concentration of Lu in the deep ocean), our CTL experiment successfully reproduced the vertical profiles observed at the Site S.

4.2 Inter-basin distribution of Nd, Dy, and Lu

To evaluate the inter-basin distribution of REEs and compare the model and the observations, the vertical sections of Nd, Dy, and Lu along the line transecting the Atlantic, Southern, and Pacific Oceans are displayed in Fig. 5. A transect similar to Fig. 5a has been used for examining the inter-basin distribution of other ocean tracers (Sarmiento and Gruber 2006; Oka 2020), which we use to discuss the distribution of REEs.

Fig. 5
figure 5

a Map showing location of vertical section along 30W in the Atlantic, 60S in the Southern Ocean, and 150W in the Pacific Ocean. The distribution of b Nd, c Dy, and d Lu concentrations along the section shown by (a). Contours and colors show results of CTL. Colored circles represent observational data located near the section (indicated by blue circles in a). Unit is pmol/kg (color figure online)

Full size image

Figure 5d indicates that the simulated distribution of Lu is similar to that of nutrients such as nitrate, phosphate, and silicate; the concentration tends to increase from the Atlantic Ocean to the Pacific Ocean, and the maximum concentration is located at a depth of approximately 2000 m in the northern North Pacific Ocean. Available observations also support this nutrient-like distribution of Lu; the observed concentrations are less than 1 pmol/kg in the deep Atlantic Ocean, increase to approximately 1.4 pmol/kg in the Southern Ocean, and reach more than 2.0 pmol/kg in the deep North Pacific Ocean. In general, the model well reproduced the observed distributions of Lu, although it somewhat overestimated the concentrations in the North Atlantic Ocean, especially around 1000 m.

The simulated distribution of Nd (Fig. 5b) is similar to that of Lu; for example, the maximum concentration appears in the deep North Pacific Ocean as seen in Lu. However, notable differences exist between Nd and Lu. First, as illustrated in Fig. 4, the concentration of Nd increases almost linearly with depth in the Pacific Ocean, and the peak of the concentration is located at the sea bottom, which is much deeper than for Lu. Second, as also indicated in Fig. 4, the surface concentrations of Nd are more severely depleted than those of Lu; for example, in the Southern Ocean, the surface concentration is depleted to some degree for Nd, whereas a relatively high concentration is maintained for Lu. Finally, the inter-basin distribution of Nd is somewhat different from the nutrient-like distribution; the concentration in the North Atlantic Ocean is higher than that in the Southern Ocean. This implies that a non-conservative process is more important for LREEs than HREEs, which we will discuss in greater detail when we examine our RCYC experiment (see Sect. 5).

The simulated distribution of Dy (Fig. 5c) has features between those of the distributions of Nd and Lu; the depth of the maximum concentration in the Pacific Ocean is deeper than that of Lu and slightly shallower than that of Nd, and the depletion of the surface concentration in the Southern Ocean is greater than that of Lu but moderate compared with that of Nd. In addition, the inter-basin contrast between the Atlantic and Pacific Oceans is weaker than that for Lu and stronger than that for Nd. These features were also confirmed by the observations.

Although our model generally captured the global distribution of the observed REE concentrations, we found some systematic differences between the model and the observations, especially for surface concentrations. In Fig. 6, we present the surface concentrations of Nd, Dy, and Lu along the sections displayed in Fig. 5a. The model significantly underestimated the concentrations of Nd and Dy in the Atlantic Ocean. Similar underestimation occurred also for the Pacific Ocean, whereas the concentrations in the Sothern Ocean were generally well reproduced by the model. Underestimation of surface Nd and Dy concentrations was commonly observed in the global oceans except the Southern Ocean as illustrated Figs. 7a and 8a. On the other hand, the surface concentration of Lu was well reproduced by the model with the exception that the model somewhat overestimated the concentration in the Southern Ocean (Figs. 6c,  9a).

Fig. 6.
figure 6

The surface (averaged between 0 and 100 m) concentrations of b Nd, c Dy, and d Lu along the section shown with red line in Fig. 5a. Unit is pmol/kg. The solid curves represent results from CTL experiment and observational surface (located between 0 and 100m) data displayed in Fig. 1a are also shown with black dots

Full size image
Fig. 7
figure 7

a Nd concentration averaged between 0 and 200 m. b Nd concentration at 1000 m. c Nd concentration at 2000 m. Contours and colors show results of CTL. Colored circles represent observational data. Unit is pmol/kg (color figure online)

Full size image
Fig. 8
figure 8

a Dy concentration averaged between 0 and 200 m. b Dy concentration at 1000 m. c Dy concentration at 2000 m. Contours and colors show results of CTL. Colored circles represent observational data. Unit is pmol/kg (color figure online)

Full size image

4.3 Global distribution of Nd, Dy, and Lu: horizontal distribution

In this subsection, we further discuss the distribution of REEs in terms of their horizontal distributions. The horizontal distributions of Nd, Dy, and Lu are presented in Figs. 7, 8, and 9, respectively. As already discussed in Sect. 4.2, the overall features of the observed REEs distribution were well reproduced by the model; however, there were some discrepancies between the model and the observations, as discussed below.

Fig. 9
figure 9

a Lu concentration averaged between 0 and 200 m. b Lu concentration at 1000 m. c Lu concentration at 2000 m. Contours and colors show results of CTL. Colored circles represent observational data. Unit is pmol/kg (color figure online)

Full size image

For Nd, the model tended to underestimate the surface concentration (Fig. 7a), especially in the North Pacific and the Atlantic, as illustrated in Fig. 6a. At depths of 1000 m (Fig. 7b) and 2000 m (Fig. 7c), the simulated concentration was also underestimated in the North Pacific, whereas overestimation was observed near the coast of the North Atlantic.

In contrast to Nd, much better model–data agreement was obtained for Lu. For the surface concentration (Fig. 9a), the model reproduced the observations well except for slight underestimation in the North Pacific regions and slight overestimation around the Polar Front (around 55°S) in the Southern Ocean. The model also exhibited good agreement with the observations at a depth of 2000 m, although concentrations in the Indian Ocean appeared to be slightly overestimated (Fig. 9c). At a depth of 1000 m (Fig. 9b), the model tended to overestimate the observations both in the Pacific and Atlantic Oceans. This appeared to be related to the fact that the simulated depth of the maximum concentration in the North Pacific was somewhat shallower than the observation (Fig. 4c).

The simulated distribution of Dy (Fig. 8) exhibits an intermediate feature between those of Nd and Lu, as illustrated in Figs. 4 and 5. Compared with the observations, the surface concentrations (Fig. 8a) were somewhat underestimated, as for Nd. At a depth of 1000 m (Fig. 8b), significant underestimation was also observed in the North Pacific, similar to Nd. Although it was not clear from the observations, the model displayed systematic variation in the deep-water concentration among LREE, MREE, and HREE; a high concentration was located primarily in the Pacific for Lu (Fig. 9c), a high concentration was located in both the Atlantic and Pacific for Nd (Fig. 7c), and a moderate inter-basin contrast was simulated for Dy (Fig. 8c).

5 Discussion

5.1 Role of internal and external processes

In Sect. 4, we demonstrated that our CTL experiment successfully reproduced the overall pattern of the observed REE concentrations. To examine how the distribution of REEs is determined in the CTL experiment, we compared the RCYC and CTL experiments. The vertical profiles of REEs simulated in the RCYC experiment are represented by thin lines in Fig. 4 (bold lines indicate the CTL experiment). For Lu, Fig. 4c, f, and i demonstrate that differences between the CTL and RCYC experiments are generally small. This means that the distribution of Lu is controlled mainly by the internal cycle (Fig. 1b; the term \(S{c}_{scav}\) in Eq. (2)) similarly to nutrients, and the external source/sink (the term \(S{c}_{in}\) in Eq. (2)) has only limited effects.

On the other hand, for Nd, Fig. 4a, d, and g indicate that there are significant differences between the CTL and RCYC experiments. When compared with CTL, Fig. 4a demonstrates that at Site P, the concentrations are more severely depleted above the thermocline, while extreme enrichment occurs in the deep ocean in RCYC. This extreme enrichment is caused by very efficient vertical transport of Nd because of the large \({K}_{\mathrm{P}}\) value for Nd. Because the globally averaged concentration is assumed unchanged in RCYC, this extreme enrichment in the North Pacific leads to extremely depleted concentrations above 2000 m in the entire Atlantic, as illustrated in Fig. 4g. The difference between the CTL and RCYC experiments is caused by the external source from the continental shelf (the term \(S{c}_{in}\) in Eq. (2)) and sink associated with sedimentation (Eq. (7)). Therefore, unlike for HREEs, external processes have a significant role in controlling the global distribution of LREEs. The role of external processes is also significant for MREEs, as indicated in Fig. 4b, e, and h.

5.2 Distribution of preformed REEs

To further examine the distribution of REEs and the differences in the distribution of individual REEs, we separated REE concentration simulated in the CTL experiment into two components: “preformed” and “regenerated” concentrations. The concept of preformed concentration is widely used for discussing nutrient and carbon cycles (Ito and Follows 2005). When surface nutrients are not completely utilized by phytoplankton at the surface ocean, the unutilized nutrients can be advected into the deeper ocean and increase the nutrient concentration in the interior ocean. In addition, the nutrient concentration in the interior ocean is directly increased by the remineralization of organic matter. The distribution of the preformed nutrients is controlled by the former process, while the distribution of the regenerated nutrient is controlled by the latter process.

Here, we apply the similar analysis to REE concentrations by separating REE concentration into preformed and regenerated components. In this study, the distribution of the preformed REEs was obtained from an additional OGCM simulation, in which the surface concentration was strongly restored to that obtained in the CTL experiment, and no source/sink term was considered (i.e., the term \(Sc\) is treated to be zero) in the interior ocean. The difference between the preformed concentration and the actual concentration corresponds to the regenerated concentration, which is interpreted as the enrichment of the REE concentration in the interior ocean by the vertical transport of REEs associated with the scavenging term (i.e., \(S{c}_{\mathrm{scav}}\)) and the external source term (i.e., \(S{c}_{\mathrm{in}}\)).

The vertical profiles of the preformed REEs are represented by thin broken lines in Fig. 4. From the definition of preformed concentration, the surface concentration becomes equivalent to the actual concentration (bold lines). At all sites, the preformed concentration tends to be larger in the deeper ocean than that in the surface ocean although there is no vertical transport of REEs associated with the scavenging. This is because the relatively high surface concentrations over deep-water formation areas are transported to the interior ocean. Because the surface concentration is much higher around the high-latitude Southern Ocean than the northern North Atlantic Ocean (Figs. 7a, 8a, 9a), the deep-water enhancement becomes more significant in the North Pacific Ocean (Site P in Fig. 4) than in the North Atlantic Ocean (Site A in Fig. 4). However, when compared with the preformed concentration (thin broken lines), the actual concentration (i.e., concentration simulated in the CTL experiment; bold lines) is much higher in the deep ocean. The regenerated REE concentration causes this enhancement, which is controlled mainly by the vertical transport of REEs due to the scavenging term. The enhancement of the regenerated concentration in the deep ocean is more significant for LREE (Nd) because of the stronger scavenging and vertical transport. This is consistent with the results from the RCYC experiment; more significant enrichment of the concentration occurs for LREE (Nd) than HREE (Lu).

Previous observational studies have reported that vertical profiles of HREEs are similar to those of silicate. Because our CTL experiment successfully reproduced the observed Lu distribution in the global ocean, here, we discuss the similarity in the distributions of Lu and silicate together with their differences from Dy and Nd. For this purpose, the distributions of preformed REEs along the sections (red line in Fig. 4a) are displayed in Fig. 10. This figure indicates that a fraction of the preformed component (indicated by contour values in Fig. 10) is larger for HREE (Lu) than for LREE (Nd). For example, in most of the Atlantic Ocean, the fraction of the preformed concentration is more than 70% for Lu, while it is less than 50% for Nd. This is expected, as stronger scavenging occurs for LREE, reducing the surface concentration and enhancing the deep-water concentration more significantly than for HREE, which leads to a smaller portion of the preformed component of LREE.

Fig. 10
figure 10

The distribution of b preformed Nd, c preformed Dy, and d preformed Lu along the section shown by Fig. 5a. Color represents for concentration of preformed REE [pmol/kg] and contour shows fraction of preformed REE concentration to the actual REE concentration [%] (color figure online)

Full size image

In addition to the scavenging efficiency, the meridional distribution of the surface concentration around the Southern Ocean appears also important. Previous studies examining the differences between the phosphate and silicate cycles indicated that significant parts of global thermocline phosphate are supplied by its preformed component over the Polar Front Zone (around 50\(^\circ\)S) in the Southern Ocean, whereas the absence of preformed silicate there causes the depletion of the thermocline silicate compared with phosphate (Sarmiento et al. 2004). This is because the Polar Front Zone is located around the formation area of Subantarctic Mode Water (SAMW), and the significant portions of the thermocline water are ventilated by SAMW (Sloyan and Rintoul 2001). Similar differences were also observed between REEs simulated in our CTL experiment. As demonstrated in Fig. 6, there are notable differences in concentrations around 50\(^\circ\)S between REEs; the concentration is significantly depleted for Nd, but retains higher values for Lu. This leads to the notable differences in the distribution of preformed Nd and Lu; more severely depleted preformed Nd and relatively abundant preformed Lu in the thermocline are obtained.

With respect to the similarity to silicate, it has been demonstrated that the similarity between zinc and silicate is caused by their similar meridional profiles of surface concentration around the Southern Ocean (Vance et al. 2017; de Souza et al. 2018). An analogical discussion can also be applied to the similarity between silicate and Lu; the similarity to silicate can be partly explained by the similarity of their preformed distributions. In fact, the high surface concentration of Lu in the Southern Ocean begins to decrease around the Polar Front Zone and reaches a low concentration around 40\(^\circ\)S (Fig. 6c), which is similar to the distributions of silicate and zinc. In addition, the scavenging process also affects the distribution of Lu. In the case of zinc, a previous study proposed that this process is additionally required for explaining the vertical profile in the North Pacific (Weber et al. 2018). This appears to be true for Lu as well because the fraction of the preformed Lu is less than 60% in the North Pacific (Fig. 10c), and the vertical profile of the preformed component is different from that of the actual Lu concentration (Fig. 4c), which implies that additional enrichment of the concentration by scavenging is necessary to reproduce the observed Lu concentration in the North Pacific.

In our model, surface Nd and Dy concentrations around 50\(^\circ\)S (especially in the Atlantic) were significantly underestimated compared with the observations, which suggests the possibility that the above-mentioned differences in the preformed concentrations of REEs may have been somewhat exaggerated in our CTL experiment. It is also important to note that the separation of preformed and regenerated REEs was possible for the simulated REEs but this separation is difficult for the observational REEs because apparent oxygen utilization (AOU), which can be used for diagnosing the regenerated nutrients, is not an accurate measure for regenerated REEs. Therefore, the distribution of preformed REEs simulated in our model (Fig. 10) is difficult to be directly compared with the observations at the current stage. To discuss the importance of the preformed component more accurately, the model bias in the surface REEs needs to be improved and more careful comparison of REEs concentrations between model and observational data will be required for confirming whether the importance of the preformed component found in our model simulation is also applicable to the real ocean.

5.3 Discussion of possible causes of model-data discrepancies

Although our CTL simulation generally reproduced the observed global distribution of REE concentrations well, there was a discrepancy between the simulated and observed REE concentrations. Here we discuss the reasons for the mode-data discrepancies.

First, although we selected the value of \({K}_{\mathrm{P}}\) based on the estimate by Nozaki (2001), there is an uncertainty in the estimate and this choice might not be ideal for reproducing the observed REE concentrations. In fact, for Lu, when we selected a slightly larger \({K}_{\mathrm{P}}\) than in our CTL experiment (i.e., we changed the value of \({K}_{\mathrm{P}}\) from 1.1 × 105 to 1.5 × 105), we can obtain a better Lu simulation than in the CTL experiment. Overestimation of the surface concentration in the Southern and Pacific Oceans in the CTL experiment (Fig. 6c) was significantly corrected, and the depth of the maximum concentration in the North Pacific Ocean (Fig. 5d) was shifted slightly downward to more closely resemble the profile of silicate (not shown). The global average of the particle fraction of REEs (\({\sim A}_{\mathrm{p}}/{A}_{\mathrm{d}}\)) reported by Nozaki (2001) (e.g., 0.46% for Lu) was utilized for selecting \({K}_{\mathrm{P}}\) in our CTL experiment; however, there was a large uncertainty in this estimate (e.g., \(0.44\pm 0.21\) % in the Atlantic Ocean and in \(0.48\pm 0.34\) % in the Indian Ocean from Table 2 in Nozaki, 2001). The choice of 1.5 × 105 for the value of \({K}_{\mathrm{P}}\) means that the value of \({A}_{\mathrm{p}}/{A}_{\mathrm{d}}\) is assumed to be 0.63% which is still within the uncertainty of the estimate by Nozaki (2001). On the other hand, for Nd and Dy, changing the value of \({K}_{\mathrm{P}}\) from that in the CTL experiment did not necessarily lead to improvement. If we chose a smaller \({K}_{\mathrm{P}}\) than that in the CTL experiment, the underestimation of the surface concentration in the Atlantic and Pacific in the CTL experiment (Fig. 6a,b) was improved while causing significant overestimation in the Southern Ocean.

The second potential reason for the model-data discrepancies is that in our CTL experiment, the external source of REEs was assumed to originate solely from the continental shelf areas presented in Fig. 3; however, additional sources not captured in Fig. 3 may also play a role. For example, it was reported that a Nd isotope anomaly was supplied from remote islands in the Pacific (Lacan and Jeandel 2001; Fröllje et al. 2016), which suggests the possibility that there are external REE inputs from small islands not captured in Fig. 3. Because underestimation of the surface concentration was globally observed in the Nd and Dy simulations, additional external sources from small islands may improve such the underestimation. Concentrations in the North Pacific regions are particularly underestimated in our CTL experiment, parts of which may originate from insufficient external REE flux there. Another possibility is that REE inputs from dust deposition (Tachikawa et al. 1999a) may also contribute to the enhancement of surface REEs concentrations in the open ocean, especially in the Atlantic Ocean. More sophisticated treatment of the external REE input should be considered in future studies.

Finally, our formulation of REE cycles described in Sect. 2.4 may require improvement, and additional processes not included in our formulation may also be necessary to more realistically reproduce the observed distributions of REE concentrations. In particular, underestimation of the surface concentration of Nd and Dy suggests that our treatment of scavenging is too strong near the surface, which cannot be improved by simply modifying the value of \({K}_{\mathrm{P}}\) as discussed above. To resolve this underestimation with minimal modification of our formulation, one approach is to treat \({K}_{\mathrm{P}}\) as a function of depth, in which the value of \({K}_{\mathrm{P}}\) is specified to have a smaller value near the surface than in the deep ocean, although such treatment appears somewhat artificial. Another approach is to incorporate the so-called “particle effect” which controls the value of \({K}_{\mathrm{P}}\). In the case of Th, it was reported that the value of \({K}_{\mathrm{P}}\) depends on the concentration of particles; for a higher particle concentration, the value of \({K}_{\mathrm{P}}\) decreases (Honeyman et al. 1988). Although the mechanism of this “particle effect” is not yet fully understood (Honeyman et al. 1988; Henderson et al. 1999; Hayes et al. 2015) and it is unclear whether the same effect is applicable to REEs, this dependence on particle concentration generally leads to smaller \({K}_{\mathrm{P}}\) at the surface than in the deep ocean and is expected to reduce the underestimation of Nd and Dy in our present model. In addition, recent studies have reported that the bottom scavenging is important for controlling the distribution of Th and Pa (Roy-Barman 2009; Okubo et al. 2012; Rempfer et al. 2017; Okubo 2018). If the bottom scavenging process is also important for REEs, the amount of the external REEs specified in our model should be increased to balance the input and output of REEs. Therefore, the inclusion of the bottom scavenging leads to an increase in the external REE flux, which is expected to increase the surface concentration and reduce the underestimation of surface Nd and Dy. With respect to possible role of additional processes not included in our formulation, recent studies have proposed the possibility that REEs are incorporated and carried with diatom opal (Akagi 2013; Nishino and Akagi 2019). If this process actually occurs in the ocean, the model needs to additionally incorporate biological uptake of REEs by diatoms. In the case of zinc, Weber et al. (2018) reported that both biological uptake and scavenging play a role in controlling the vertical distribution of zinc in the Pacific Ocean; the similar discussion may be relevant for Lu in such a model.

6 Summary and concluding remarks

In this study, from OGCM simulations, we successfully demonstrated that the global distribution of REEs can be reproduced by considering the internal cycle associated with reversible scavenging and external REEs inputs around continental regions. In our simulations, the external REE inputs were assumed to be proportional to Nd flux such that the total REE input was consistent with the estimate from observational data (Nozaki 2001). The Nd flux was obtained from the observed Nd database around coastal areas (Jeandel et al. 2007) in the same way as previous modeling of Nd and its isotopes (Arsouze et al. 2007; Jones et al. 2008). Reversible scavenging was incorporated into an OGCM, as described in Oka et al., 2009. The parameter controlling the strength of the affinity to particles (\({K}_{\mathrm{P}}\) in Eq. (4)) depends on individual REEs, and we set the value of \({K}_{\mathrm{P}}\) to be consistent with the observed ratios of particulate REEs to dissolved REEs reported in Nozaki (2001). By specifying different values of \({K}_{\mathrm{P}}\) for different REEs (Table 1), the model can represent different affinities to particles depending the REEs: lighter REEs can have a stronger affinity to particles by setting the value of \({K}_{\mathrm{P}}\) to be larger than that for heavier REEs.

As reported in a previous study (Oka et al. 2009), the vertical profile of REEs in the North Pacific is strongly controlled by the reversible scavenging process, and the systematic differences in the vertical profiles of REEs can be reproduced by the model with an appropriate choice of the parameter \({K}_{\mathrm{P}}\). We newly demonstrated that the external REE input around the coastal regions also plays a key role in controlling the vertical profiles and their inter-basin differences. The role of the external inputs is especially important for LREEs (e.g., Nd); the linear increase in Nd concentration in the North Pacific cannot be reproduced only by the reversible scavenging (thin line in Fig. 4a), but rather by a combination of the reversible scavenging and the external inputs (bold line in Fig. 4a). By contrast, the distribution of HREEs (e.g., Lu) can be broadly reproduced without external inputs, suggesting that Lu behaves as a conservative nutrient-like tracer. In terms of the similarity between the distribution of Lu and silicate, our analysis of the preformed components confirmed that the similarity in their surface distributions around the Southern Ocean was one of the reasons. In addition, our analysis suggested that the reversible scavenging process is also required to explain this similarity especially in the Pacific Ocean.

Our CTL experiment successfully reproduced the observations in the REE database (compiled from both the recently obtained GEOTRACES dataset and pre-GEOTRACES observations; see Fig. 1) in terms of the inter-basin distribution (Fig. 5). The observations suggest that the vertical profiles of REEs are not the same among the basins; our model simulations demonstrated that the observed distributions can be clearly reproduced by considering the reversible scavenging and external REE inputs from the coastal regions. However, our simulations had limitations: the surface concentrations of Nd and Dy tended to be underestimated in the Atlantic and the lower-latitude Pacific Oceans, while enhanced concentrations in the North Pacific Ocean were not sufficiently reproduced. In addition, local high concentrations reported in the observations were not always well reproduced. To improve the model simulations, modification of our formulation for the reversible scavenging and more sophisticated treatment regarding the external REEs input should be addressed in future work.

Our REE simulations demonstrated that the global distributions of REEs in the ocean could be explained by a combination of the reversible scavenging and the REE input around the coastal regions. As a natural extension of our REE simulations, we can simulate the Nd isotope ratio by individually specifying the external inputs of 144Nd and 143Nd from the coastal regions in the same way as REEs. The Nd isotope ratio has been used as a paleo proxy for the ocean circulation (Böhm et al. 2015) and its simulation is important for understanding not only the modern ocean but also for the past ocean. For the modern ocean, the amount of observed data on the Nd isotopes has recently increased owing to the GEOTRACES program, and we plan to perform Nd isotope simulation and a comparison with the recent observations in a future study. As another extension of our study, because the other trace elements are supplied from the coastal regions in the same way as REEs, our obtained REEs fluxes (illustrated in Fig. 3b) may be useful for simulating other trace elements. For example, there has been large uncertainty about the source of iron in the global iron simulation (Tagliabue et al. 2016), and our obtained REEs fluxes could be utilized as a boundary source of iron, although information on the ratio of iron to REEs would be required. Parallel simulations of REEs and iron under a common boundary source around coastal regions may provide a clearer picture of the iron cycle in the ocean.