GNOM v1.0: an optimized steady-state model of the modern marine neodymium cycle

Spatially distant sources of neodymium (Nd) to the ocean that carry different isotopic signatures (ε_Nd) have been shown to trace out major water masses and have thus been extensively used to study large-scale features of the ocean circulation both past and current. While the global marine Nd cycle is qualitatively well understood, a complete quantitative determination of all its components and mechanisms, such as the magnitude of its sources and the paradoxical conservative behavior of ε_Nd, remains elusive. To make sense of the increasing collection of observational Nd and ε_Nd data, in this model description paper we present and describe the Global Neodymium Ocean Model (GNOM) v1.0, the first inverse model of the global marine biogeochemical cycle of Nd. The GNOM is embedded in a data-constrained steady-state circulation that affords spectacular computational efficiency, which we leverage to perform systematic objective optimization, allowing us to make preliminary estimates of biogeochemical parameters. Owing to its matrix representation, the GNOM model is additionally amenable to novel diagnostics that allow us to investigate open questions about the Nd cycle with unprecedented accuracy. This model is open-source and freely accessible, is written in Julia, and its code is easily understandable and modifiable for further community developments, refinements, and experiments.

1 Introduction

Rare earth elements (REEs) have long been recognized to provide unique insight into ocean circulation and biogeochemical cycles (e.g., de Baar et al., 1983, 1985; Bertram and Elderfield, 1993; Elderfield, 1988; Elderfield and Greaves, 1982; German et al., 1995; Goldberg et al., 1963; Haley et al., 2014; Høgdahl et al., 1968; Lacan and Jeandel, 2001, 2004; Piepgras and Jacobsen, 1992; Piper, 1974; Sholkovitz and Schneider, 1991; Zheng et al., 2016). Isotopic variations in neodymium (Nd), in particular, have been extensively used as a tracer of ocean circulation, which plays a fundamental role in Earth's climate over a wide range of timescales, from millennia to millions of years (e.g., Adkins, 2013; van de Flierdt et al., 2016a; Frank, 2002; Goldstein and Hemming, 2003; Piepgras and Wasserburg, 1980; Sigman et al., 2010; Tachikawa et al., 2017).

Neodymium is part of a long-lived isotope system. Samarium-147 (¹⁴⁷Sm) decays to neodymium-143 (¹⁴³Nd) with a half-life of 106 Gyr. While the Sm:Nd ratio varies within the earth, these ratios are remarkably similar in most rocks in the continental crust and across the geological timescale and about 40 % lower than the bulk earth Sm:Nd (DePaolo and Wasserburg, 1976; McCulloch and Wasserburg, 1978; Goldstein et al., 1984), and as a result, the ε_Nd values in continental rocks generally directly reflect the average crustal age. Therefore, the Nd isotope ratio $R={}^{\mathrm{143}}\mathrm{Nd}{/}^{\mathrm{144}}\mathrm{Nd}$ is mainly a reflection of the amount of time the Nd in a rock has been a part of the continental crust, with lower values indicating older ages and longer crustal residence times (DePaolo and Wasserburg, 1976; McCulloch and Wasserburg, 1978; Goldstein and Hemming, 2003; Jeandel et al., 2007; van de Flierdt et al., 2016a; Robinson et al., 2021). Because R variations are typically small, Nd isotope signatures are usually defined as follows:

$$\begin{array}{}\text{(1)}& {\mathit{\epsilon }}_{\mathrm{Nd}}=R/R_{\mathrm{CHUR}}-\mathrm{1},\end{array}$$

expressed in parts per 10 000 (; DePaolo and Wasserburg, 1976), where R is the measured ${}^{\mathrm{143}}\mathrm{Nd}{/}^{\mathrm{144}}\mathrm{Nd}$ ratio, and the chondritic uniform reservoir (CHUR) represents an estimate of the average Nd isotope ratio of chondritic meteorites and the bulk earth. For consistency with previously published data, we use R_CHUR=0.512638 from Jacobsen and Wasserburg (1980), rather than the updated value from Bouvier et al. (2008).

Early measurements of ε_Nd in seawater (Piepgras and Wasserburg, 1980) and ferromanganese oxide crusts (Elderfield et al., 1981; Goldstein and O'Nions, 1981; O'Nions et al., 1978; Piepgras et al., 1979) showed systematic variation across the ocean basins, with the lowest ε_Nd values in the North Atlantic (−14  to −10 ), the highest values in the Pacific (−5  to 0 ), and intermediate values in the Southern Ocean (−11  to −8 ). The latter value broadly reflects mixing between waters from the North Atlantic, which are influenced by old continental terrains in northern Canada and Greenland, and the Pacific, which is influenced by mantle-derived volcanics (van de Flierdt et al., 2016a; Frank, 2002; Garcia-Solsona et al., 2014; Goldstein and Hemming, 2003; Goldstein and O'Nions, 1981; Lambelet et al., 2016; Piepgras and Wasserburg, 1980; Stichel et al., 2012a). These observations led to the recognition that ε_Nd values could be used to trace mixing between North Atlantic and Pacific waters over time, thus making ε_Nd a potentially powerful paleoceanographic tracer.

More recently, the GEOTRACES program was created to better understand the sources and cycling of trace elements and isotopes in the ocean and how they impact broader marine biogeochemical cycles. The GEOTRACES Science Plan (GEOTRACES Planning Group, 2006) identified Nd isotopes as a key trace element or isotope that is expected to be measured on all GEOTRACES cruises because of its use as a paleoceanographic proxy. Thanks to this international effort, considerable amounts of new Nd concentration and isotope data have been generated in recent years, collected notably in the GEOTRACES Intermediate Data product 2017 (IDP17; Schlitzer et al., 2018), which increased the Nd data inventory by about 50 %, and post-IDP17 GEOTRACES data yet to be released in future data products.

To gain the most useful and accurate information from these observations, however, it is paramount to understand their modern ocean biogeochemical cycles and tracer budgets. Neodymium and other REEs enter the ocean via rivers, submarine groundwater discharge, eolian deposition, pore waters, and/or interaction with sediments (Fig. 1). Once in the ocean, they are redistributed by the ocean circulation, scavenged by particulate matter, and exit the ocean via sedimentation and incorporation into authigenic ferromanganese oxides (Frank, 2002; Byrne and Kim, 1990; Elderfield, 1988; Elderfield et al., 1981; Elderfield and Sholkovitz, 1987; Sholkovitz et al., 1989, 1994, 1992; Haley et al., 2004; Blaser et al., 2016; Du et al., 2016). While most sources of REEs to the ocean have likely been identified, there are still large uncertainties associated with the magnitudes of these different fluxes due to the inherent challenges of measuring sources that are temporally variable and globally widespread.

Models of the marine Nd cycle, in conjunction with seawater measurements, offer a way to constrain the magnitudes and isotopic compositions of these various inputs to the ocean and identify the most important sources. Of the distinct types of models, the following four have been used to simulate the modern ocean Nd cycle: simple box models, ocean general circulation models (OGCMs), steady-state circulation models, and boundary propagation models, each with their strengths and weaknesses (Table 1). Some of these models have explicitly tracked the concentrations of each Nd isotope (¹⁴³Nd and ¹⁴⁴Nd, thus allowing for the estimation of the Nd concentration and its isotopic composition), while other models have simply tracked ε_Nd as a single conservative tracer.

Bertram and Elderfield (1993)Tachikawa et al. (2003)Du et al. (2020)Arsouze et al. (2007)Arsouze et al. (2009)Arsouze et al. (2010)Rempfer et al. (2011)Gu et al. (2019)Gu et al. (2020)Pöppelmeier et al. (2020)Jones et al. (2008)Siddall et al. (2008)Du et al. (2020)

Table 1 Previous modeling studies. Note: GCM is general circulation model, and TMI is total matrix intercomparison.

^a Only the Pacific is modeled. ^b With deep boundary condition. ^c With surface boundary condition. ^d With deep and/or surface boundary conditions. ^e Bulk residence time estimated from the total Nd source magnitude for an ocean volume of 1.32×10¹⁸ m³ and a mean Nd concentration of 22 pM. ^f Total exterior surface flux calculated from their model, 90 % of which is missing compared to their estimate based on observations.

Box models typically refer to models consisting of 10 or fewer well-mixed boxes that exchange tracer with each other through prescribed mixing and overturning rates. Owing to their small size, box model simulations are the fastest to run and require very little computational power. Thus, they facilitate parameter optimization and scientific exploration by allowing for quick experimentation. For example, box models have been successfully used to determine that Nd must exchange between seawater and particles in the water column or at the sediment–water interface (Bertram and Elderfield, 1993) and that riverine and eolian sources are not sufficient to explain regional ε_Nd variability (Tachikawa et al., 2003). However, very low spatial resolution prevents box models from capturing many important features of ocean circulation.

Ocean general circulation models sit on the other end of the spectrum of computational complexity, with better spatial resolution and resolved physics. Their computational costs generally prohibit systematic parameter space exploration or parameter optimization. These models have thus been used primarily to run well-defined experiments that target specific hypotheses, such as the importance of continental margin sources (boundary exchange) on ε_Nd distributions, either as the sole source of Nd to the ocean (Arsouze et al., 2007, 2010) or an additional source to rivers and eolian deposition (Arsouze et al., 2009; Rempfer et al., 2011; Gu et al., 2019, 2020). To our knowledge, only Rempfer et al. (2011) and Gu et al. (2019) have attempted to optimize a Nd cycling model, using the low-resolution Bern3D OGCM and 3^∘ resolution Community Earth System Model (CESM1.3), respectively, and each optimizing only two parameters.

More recently, a new class of steady-state models has emerged with unique potential to combine the advantages of OGCMs with the computational speed of box models. These models do not resolve the physics at runtime and, instead, rely on a prescribed, steady-state circulation. They can thus directly solve for the steady-state solution of the three-dimensional tracer equations, avoiding costly spin-ups, and drastically reducing simulation times. Thus far, to our knowledge, these models have only been used to test the top-down hypothesis by propagating a surface boundary condition into the ocean interior. Using the transport matrix method (TMM; Khatiwala et al., 2005; Khatiwala, 2007), Jones et al. (2008) showed that conservative mixing and advection from the surface alone cannot reproduce interior ε_Nd observations, while Siddall et al. (2008) showed that including reversible scavenging captures the observed decoupling between quasi-conservative ε_Nd and nutrient-like Nd concentration ([Nd]) distributions (the Nd paradox; Goldstein and Hemming, 2003).

The fourth class of models, which we have termed boundary propagation models, entirely bypasses expressing fluxes between model grid cells by connecting interior grid cells directly to the surface, using the total matrix intercomparison method (TMI; Gebbie and Huybers, 2010). Specifically, boundary propagation models estimate the fractional contribution of each surface grid cell to each interior grid cell. These models have been used to explicitly test the conservativeness of Nd isotopes as a tracer, since they do not incorporate external fluxes of Nd or internal cycling processes, and can thus only be used to simulate conservative transport. Indeed, similar to the experiment of Jones et al. (2008) referenced above, Du et al. (2020) used the TMI to inquire how well interior ε_Nd values can be explained by conservative mixing and advection alone.

Our goal is to fill the current gap in the marine Nd modeling landscape and leverage the largely unexplored benefits of steady-state circulation models. Hence, here, we present the Global Neodymium Ocean Model (GNOM) v1.0, a mechanistic model of the modern ocean Nd cycle embedded in a state-of-the-art, steady-state estimate of the modern ocean circulation from the Ocean Circulation Inverse Model version 2 (OCIM v2.0; DeVries and Primeau, 2011; DeVries, 2014; DeVries and Holzer, 2019). The computational efficiency afforded by the model allows us to objectively optimize the model's parameters, making GNOM v1.0 the first inverse model of the Nd cycle and producing a good match to observations.

The GNOM v1.0 thus provides the community with a realistic yet computationally affordable tool to model the marine Nd cycle that we hope will be used to further improve our understanding of Nd cycling in the ocean. The model code and its optimization script are available publicly on GitHub at https://github.com/MTEL-USC/GNOM (last access: 26 May 2022). We used the free and open-source Julia language (Bezanson et al., 2017) and its packages, AIBECS.jl (https://github.com/JuliaOcean/AIBECS.jl, last access: 26 May 2022) in particular (Pasquier, 2020a; Pasquier et al., 2022b), as our main development platform. Owing to its open-source design, simplicity, and computational speed, the GNOM v1.0 is ideal for Nd cycle investigations. Except for the GEOTRACES dataset which must be downloaded manually, the GNOM is self-contained and version controlled, making it easy to reproduce simulations.

Additionally, the steady-state formulation of the GNOM is amenable to novel Green-function-based diagnostics that can provide important new insights into major open questions on the marine Nd cycle. Green functions (sometimes spelled Green's functions) can be used for solving ordinary differential equations with an initial condition and/or boundary values (see, e.g., Morse et al., 1953). In our case, they can be thought of as the [Nd] responses to unit local sources of Nd and allow us to partition [Nd] or ε_Nd into components of interest, such as Nd from a particular source or location. Here, we introduce new partitions of Atlantic Nd and ε_Nd (following, e.g., Holzer et al., 2016; Pasquier and Holzer, 2017, 2018; Holzer et al., 2021) that are helpful for disentangling the neodymium paradox (Siddall et al., 2008). We show that we can accurately partition [Nd] and ε_Nd in the central Atlantic into contributions from northern- and southern-sourced waters. These preliminary diagnostics already reveal important information. They help quantify the conservativeness of ε_Nd along water pathways and unveil underlying mechanisms by evaluating the effects of local sources and sinks. Detailed investigations of these diagnostics are out of the scope of this study and will be carried out in future work using a subsequent version of the GNOM with more finalized parameter values. We invite paleoceanographers and modelers alike to use the GNOM v1.0 model, to improve its implementation, to explore its capabilities, and thus to contribute to quantitatively answering long-lasting questions on the Nd cycle.

2 The GNOM model

Neodymium concentrations are controlled by the interplay between circulation, external sources, and reversible scavenging and burial in the sediments (Fig. 1). These components completely define the state of the Nd cycle in our Global Neodymium Ocean Model (GNOM) v1.0. The three-dimensional partial differential equation for the Nd concentration tracer is discretized onto the grid of the Ocean Circulation Inverse Model (OCIM v2.0; DeVries and Holzer, 2019), yielding a system of 200 160 ordinary differential equations. Reorganizing the discretized three-dimensional arrays into column vectors, the steady-state tracer equation is recast in matrix form, as follows:

$$\begin{array}{}\text{(2)}& \left({\mathbf{T}}_{\mathrm{circ}}+{\mathbf{T}}_{\mathrm{scav}}\right){\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}=\sum _k{\mathit{s}}_k,\end{array}$$

where ${\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}$ is the modeled Nd concentration vector, T_circ is the OCIM v2.0 advection–diffusion operator or transport matrix, T_scav is the reversible-scavenging matrix, and the s_k are the external sources of neodymium. Note that ${\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}$ and s_k are 200 160-element column vectors and that T_circ and T_circ are sparse 200 160×200 160 matrices, such that the linear system represented by Eq. (2) can be solved in a few seconds on a modern laptop via lower–upper (LU) factorization and forward and backward substitution (often referred to as matrix inversion).

Figure 1 Diagram of the Nd cycle model as implemented in GNOM v1.0. External sources of dissolved Nd are represented by black arrows. Localized sources, rivers, groundwater, and hydrothermal vents are indicated by a small circle at the origin of their respective arrows. A fraction of Nd is reversibly scavenged and pumped downwards. A fraction of scavenged Nd that reaches the sediments is buried in the sediments and removed from the system. Nd is also continuously transported by the ocean circulation model (not represented in the schematic).

The global ε_Nd distribution is determined by both the distribution of ¹⁴³Nd and ¹⁴⁴Nd. Following, e.g., John et al. (2020), instead of explicitly simulating two additional tracers, we recover ε_Nd values by simulating a single additional fictitious tracer for R[Nd], which we denote by RNd (and its column vector by ${\mathit{\chi }}_{\mathrm{RNd}}^{\mathrm{mod}}$). This is equivalent to assuming that ¹⁴⁴Nd:Nd is constant, such that ${\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}$ and ${\mathit{\chi }}_{\mathrm{RNd}}^{\mathrm{mod}}$ nominally track [¹⁴⁴Nd] and [¹⁴³Nd], respectively, multiplied by this constant ¹⁴⁴Nd:Nd. We omit stable isotope fractionation during scavenging because its effect is negligible compared to the effect of radioactive decay from ¹⁴⁷Sm. Thus, in Eq. (2), only the external sources s_k differ in their isotopic composition. Thus, in practice, ${\mathit{\chi }}_{\mathrm{RNd}}^{\mathrm{mod}}$ is computed by solving Eq. (2) with the sources replaced by R_ks_k (element-wise multiplication), where R_k is the vector of the isotopic ratio of Nd injected by source k. The modeled ε_Nd values are then given by the vector ${\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{mod}}=\frac{{\mathit{\chi }}_{\mathrm{RNd}}^{\mathrm{mod}}/{\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}}{R_{\mathrm{CHUR}}}-\mathrm{1}$ (where all the operations are element-wise).

2.1 Ocean circulation

The ${\mathbf{T}}_{\mathrm{circ}}{\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}$ term in Eq. (2) captures the flux divergence of [Nd] as it is carried along the mean ocean currents of the model and mixed by subgrid-scale eddies. The advection–diffusion operator T_circ is represented as a 200 160×200 160 sparse matrix. (Most of the entries of T_circ are zero because water can only travel directly between neighboring grid cells.) It comes from the output of the OCIM v2.0 (DeVries and Holzer, 2019), which provides a state-of-the-art data-assimilated steady-state ocean circulation (DeVries and Primeau, 2011; DeVries, 2014; DeVries and Holzer, 2019). Physically, T_circ can be interpreted as the equivalent of $\mathrm{\nabla }\cdot (\mathit{u}-\mathbf{K}\mathrm{\nabla })$, where u is the climatological mean water velocity field, and K is an eddy-diffusivity matrix of which the horizontal component is slanted along isopycnals. The spatial resolution of its grid is fixed at a nominal 2^∘ × 2^∘ in the horizontal and consists of 24 vertical levels of increasing height with depth. We emphasize that the OCIM v2.0 is particularly suited to this type of model because it arguably provides the best available estimate of the current climate long-term large-scale ocean circulation, while it affords spectacular computational efficiency.

2.2 External sources

The GNOM v1.0 explicitly represents the following six sources of Nd into the ocean (Fig. 1): (i) atmospheric mineral dust deposition, s_dust, (ii) atmospheric volcanic ash deposition, s_volc, (iii) riverine discharge, s_river, (iv) groundwater discharge, s_gw, (v) sedimentary remobilization (including pore water fluxes), s_sed, and (vi) hydrothermal-vent release, s_hydro. The column vectors s_k summed together constitute the total source of Nd (Eq. 2). Each source term is detailed in the following sections. Their spatial patterns and isotopic signatures are shown on Fig. 2, and their magnitudes and contributions to the total inventory of Nd are collected in Table 3.

Figure 2 Vertically integrated Nd sources and corresponding vertical mean ε_Nd.

2.2.1 Eolian dust

We assume that atmospheric dust deposition injects Nd in the surface ocean only. That is, soluble Nd from dust is instantly released as dissolved Nd in the top layer of the model grid. Although it can vary with location and mineralogy (Goldstein et al., 1984), for simplicity, we assume a constant dust Nd content of $(\mathrm{Nd}:\mathrm{dust})=\mathrm{40}$ µg g⁻¹ (which is within the 11.93 to 45.76 ppm range of atmospheric dust observations of Goldstein et al., 1984). The spatial pattern of the dust source is prescribed by an atmospheric model output (Scanza et al., 2018) and is shown in Fig. 2a.

The isotopic signature of atmospheric mineral dust deposited on the ocean surface is not homogeneous (Goldstein et al., 1984). Instead, dust ε_Nd varies with composition and mineralogy, which derives from its land origin. It is also likely that Nd solubility varies with composition and mineralogy. Thus, the GNOM v1.0 uses nine separate annual mineral dust deposition fields (dataset available from Adebiyi et al., 2020) from nine different regions. These dust deposition fields were generated by Kok et al. (2021a) and Kok et al. (2021b), who partitioned dust emissions according to nine different regions of origin, using the global climate model of Scanza et al. (2018). The nine regions we use are northwestern Africa (NWAf), northeastern Africa (NEAf), southern Sahara and Sahel (Sahel), Middle East and Central Asia (MECA), East Asia (EAsia), North America (NAm), Australia (Aus), South America (SAm), and southern Africa (SAf). Figure 3 shows the extent of these regions.

Figure 3 Extent of the dust regions of origin (Kok et al., 2021a, b; Adebiyi et al., 2020) and their ε_Nd values as optimized in GNOM v1.0. See the text for region names.

We assign a distinct Nd solubility and isotopic signature to each region of origin, controlled by the 2×9 corresponding parameters (denoted β_r and ε_r for each region r; see the parameters in Table 2). The dust source of Nd into the ocean is hence given by the following:

$$\begin{array}{}\text{(3)}& {\mathit{s}}_{\mathrm{dust}}=\sum _r{\mathit{\beta }}_r{\displaystyle \frac{(\mathrm{Nd}:\mathrm{dust}){\mathit{\varphi }}_{\mathrm{dust},r}}{M_{\mathrm{Nd}}\mathrm{\Delta }{z}_{\mathrm{1}}}}(\mathit{z}=z_{\mathrm{1}}),\end{array}$$

where ϕ_dust,r is the dust deposition flux from region r taken from the Adebiyi et al. (2020) dataset and rearranged into a 200 160-element vector, Δz₁ and z₁ are the height and depth of the top layer of the model grid, z is the 200 160-element vector of depths, and M_Nd=144.24 g mol⁻¹ is the molar mass of Nd. (All the operations in Eq. 3 are element-wise, and (z=z₁) acts like a mask so that s_dust only injects Nd in the top layer of the model grid.)

Each isotopic signature parameter ε_r uniquely defines the isotopic ratio of each region via $R_r=R_{\mathrm{CHUR}}({\mathit{\epsilon }}_r+\mathrm{1})$, which is then used to compute the dust source for the RNd tracer via the following:

$$\begin{array}{}\text{(4)}& {\mathit{R}}_{\mathrm{dust}}{\mathit{s}}_{\mathrm{dust}}=\sum _r{R}_r{\mathit{\beta }}_r{\displaystyle \frac{(\mathrm{Nd}:\mathrm{dust}){\mathit{\varphi }}_{\mathrm{dust},r}}{M_{\mathrm{Nd}}\mathrm{\Delta }{z}_{\mathrm{1}}}}(\mathit{z}=z_{\mathrm{1}}).\end{array}$$

This allows for the eolian dust source to carry an elaborate and more realistic isotopic signature than previous models (Fig. 2g). Figure 3 also shows the optimized ε_r values of each region.

2.2.2 Volcanic ash

Despite a smaller atmospheric loading than mineral dust, we include volcanic ash as a separate, potentially important, eolian source of Nd because of its typically high reactivity and solubility compared to mineral dust. This reactivity partly reflects the high surface area of volcanic ash and the thermodynamic instability of volcanic glass (Gaillardet et al., 1999; Dessert et al., 2003). We use the geographic pattern of volcanic ash deposition, as used in the work of Chien et al. (2016) and Brahney et al. (2015), which provides estimates of the global deposition fields of dust and soluble iron from different aerosol types (mineral dust, volcanic ash, combustion fire, and so on). Assuming a constant neodymium content identical to dust, the volcanic ash source of Nd into the ocean is thus given by the following:

$$\begin{array}{}\text{(5)}& {\mathit{s}}_{\mathrm{volc}}={\mathit{\beta }}_{\mathrm{volc}}{\displaystyle \frac{(\mathrm{Nd}:\mathrm{dust}){\mathit{\varphi }}_{\mathrm{volc}}}{M_{\mathrm{Nd}}\mathrm{\Delta }{z}_{\mathrm{1}}}}(\mathit{z}=z_{\mathrm{1}}),\end{array}$$

where ϕ_volc is the column vector of the volcanic ash deposition flux from the Chien et al. (2016) dataset, and β_volc is the Nd solubility in volcanic ash. Similar to the dust source formulation, the magnitude of the volcanic ash source of RNd is controlled by the parameter ε_volc, as follows:

$$\begin{array}{}\text{(6)}& {\mathit{R}}_{\mathrm{volc}}{\mathit{s}}_{\mathrm{volc}}=R_{\mathrm{volc}}{\mathit{\beta }}_{\mathrm{volc}}{\displaystyle \frac{(\mathrm{Nd}:\mathrm{dust}){\mathit{\varphi }}_{\mathrm{volc}}}{M_{\mathrm{Nd}}\mathrm{\Delta }{z}_{\mathrm{1}}}}(\mathit{z}=z_{\mathrm{1}}),\end{array}$$

where $R_{\mathrm{volc}}=R_{\mathrm{CHUR}}({\mathit{\epsilon }}_{\mathrm{volc}}+\mathrm{1})$. (Note that R_volc=R_volc everywhere because the volcanic ash source comprises a single term, unlike the region-of-origin partitioned dust source.) The geographical patterns of the volcanic ash source and its uniform isotopic signature are shown in Fig. 2b and h.

2.2.3 Sediments

Sedimentary Nd is likely released via pore waters located in the upper few centimeters below the seafloor (e.g., Elderfield and Sholkovitz, 1987; Sholkovitz et al., 1989; Haley et al., 2004; Lacan and Jeandel, 2005; Wilson et al., 2013; Haley et al., 2017; Abbott et al., 2015a, b; Du et al., 2016, and references therein). The flux magnitude of this sedimentary release likely depends on sediment composition and reactivity (Lacan and Jeandel, 2005; Pearce et al., 2013; Wilson et al., 2013; Blaser et al., 2016, 2020). Other sedimentary environmental factors also likely play a role, such as oxygenation and organic matter flux (Elderfield and Sholkovitz, 1987; Sholkovitz et al., 1989, 1992; Haley et al., 2004; Lacan and Jeandel, 2005; Wilson et al., 2013). At high latitudes, mechanical glacial erosion likely increases sedimentary Nd fluxes by exposing fresh material and increasing surface area by producing fine particulates (Anderson, 2005; von Blanckenburg and Nägler, 2001), while increased bottom water eddy-kinetic energy may also enhance Nd release (Lacan and Jeandel, 2005; Gardner et al., 2018; Pöppelmeier et al., 2019).

There is no established quantitative flux model for sedimentary Nd release that works on the global scale, especially given the limited spatial coverage of direct sedimentary flux measurements (which are almost entirely restricted to the coastal northwestern Pacific). Therefore, the GNOM v1.0 implements the sedimentary Nd flux into the ocean as a flexible and optimizable function of depth z and local sedimentary ε_Nd. The base sedimentary Nd flux, ϕ(z), is modeled as an exponential function of depth, as follows:

$$\begin{array}{}\text{(7)}& {\mathit{\varphi }}_{\mathrm{sed}}\left(z\right)=({\mathit{\varphi }}_{\mathrm{0}}-{\mathit{\varphi }}_{\mathrm{\infty }})e^{-z/z_{\mathrm{0}}}+{\mathit{\varphi }}_{\mathrm{\infty }},\end{array}$$

where ϕ₀, ϕ_∞, and z₀ are optimizable parameters. The rationale behind the parameterization of Eq. (7) is versatility. For ϕ_∞<ϕ₀ and small z₀, the flux profile is larger near the surface and smaller in the deepest parts of the ocean, while for ϕ_∞>ϕ₀ and large z₀, the sedimentary flux increases quasi-linearly with depth (as in, e.g., Du et al., 2020, Fig. 1c). The optimization only enforces weak direct constraints on these parameters, allowing for any such profile shape. The optimized base Nd flux profile as a function of depth is shown in Fig. 4c.

Figure 4 (a) Sedimentary source enhancement as a quadratic function of observed ε_Nd. (b) Difference between effectively released ε_Nd and in situ ε_Nd as a function of in situ ε_Nd. (c) Sedimentary source flux profile as a function of depth. (d) Horizontally integrated sedimentary source. (e) Map of sedimentary reactivity/scaling factor α(ε_sed).

The base sedimentary flux is further scaled by a reactivity factor, α, which controls the effective sedimentary Nd release. Sediment reactivity is modeled via a simple parameterized quadratic function of local sedimentary ε_Nd, as follows:

$$\begin{array}{}\text{(8)}& \mathit{\alpha }\left({\mathit{\epsilon }}_{\mathrm{Nd}}\right)=a{\left({\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{Nd}}-c}{{\mathit{\epsilon }}_{\mathrm{10}}}}\right)}^{\mathrm{2}}+\mathrm{1},\end{array}$$

where the optimizable parameters a and c control the curvature and minimum of the quadratic, respectively, while ε₁₀=10  is a normalization constant. (The specific value of ε₁₀ is unimportant because it is absorbed by the optimizable parameter a during the optimization). Sedimentary ε_Nd is taken from a modified version of the interpolated global map of sedimentary ε_Nd of Robinson et al. (2021). (Our modification caps the central and northern Pacific ε_Nd values to a minimum of −5  because it appears that Robinson et al. (2021) artificially disconnected seafloor areas at the 180^∘ meridian, resulting in ε_Nd values that we flagged as too negative.) This quadratic parameterization is motivated by the fact that extreme sedimentary ε_Nd values are often associated with rather fresh, and thus reactive, detrital material. We emphasize that this enhancement can be turned on or off, depending on the choice of parameters (a=0 turns it off). However, maybe coincidentally, extremely high ε_Nd values are generally associated with relatively young volcanic Nd that is more reactive and readily soluble (Lacan and Jeandel, 2005; Pearce et al., 2013; Wilson et al., 2013; Blaser et al., 2016, 2020), and previous model studies have resorted to different enhanced Nd release parameterizations to achieve a similar effect (see, e.g., Pöppelmeier et al., 2020). The same is not necessarily true for rocks with extremely low ε_Nd values; however, it so happens that much of the region around the Labrador Sea (Greenland and northern Canada) is currently, or was previously, glaciated, which has resulted in a large amount of fine-grained crystalline (and thus labile) detritus with extremely negative ε_Nd (von Blanckenburg and Nägler, 2001). The quadratic function α(ε_Nd) is shown in Fig. 4a, and the resulting scaling factor for the global map is show in Fig. 4e.

Finally, to account for large glaciers that may produce fine-grained glacial flour from previously unexposed bedrock that likely contains reactive Nd (von Blanckenburg and Nägler, 2001), we additionally scale Nd release from the sedimentary source along the coast of Greenland by a factor α_GRL. (For simplicity, we did not account for potentially enhanced Nd release in Antarctic because we assume that extreme ε_Nd released by sediments in the Antarctic would be relatively rapidly mixed along the circumpolar current.) Combined with the reactivity α, the resulting sedimentary source is given by the following:

$$\begin{array}{}\text{(9)}& {\mathit{s}}_{\mathrm{sed}}={\mathit{\alpha }}_{\mathrm{GRL}}\mathit{\alpha }\left({\mathit{\epsilon }}_{\mathrm{sed}}\right){\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{sed}}\left({\mathit{z}}_{\mathrm{bot}}\right)}{\mathrm{\Delta }\mathit{z}}}{\mathit{f}}_{\mathrm{topo}},\end{array}$$

where α_GRL is the vector equal to α_GRL for grid cells against the coast of Greenland and equal to 1 otherwise, ε_sed is the column vector of the ε_Nd map of Robinson et al. (2021) repeated at all depths throughout the water column, z_bot is the vector of the depths of the bottom of each model grid cell, Δz is the vector of the height of each model grid cell, and f_topo is a mask equal to 1 for grid cells on the seafloor and equal to 0 otherwise. (Functions and operations in Eq. 9 are applied element-wise.) The horizontally integrated sedimentary source is shown Fig. 4d.

For the RNd sedimentary flux, we use the interpolated seafloor map of ε_Nd values from Robinson et al. (2021) (modified with a −5  minimum in the Pacific north of 40^∘ S). For grid cells with heterogeneous ε_Nd in sediments, our quadratic implementation of the reactivity α as a function of ε_Nd implies a statistical shift in the mean released ε_Nd towards extreme values because extremely light or heavy Nd is released more efficiently. Assuming that ε_sed represents the observed in situ mean of a normally distributed ε_Nd sediment composition in each grid cell and a uniform global standard deviation σ_ε within each grid cell, the ε_Nd that is effectively released at any location, denoted ${\mathit{\epsilon }}_{\mathrm{sed}}^{\mathrm{eff}}$, is given by the following:

$$\begin{array}{}\text{(10)}& {\mathit{\epsilon }}_{\mathrm{sed}}^{\mathrm{eff}}={\displaystyle \frac{a({\mathit{\epsilon }}_{\mathrm{sed}}-\mathrm{2}(c-{\mathit{\epsilon }}_{\mathrm{sed}}\left)\right){\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}+\left(a\right(c-{\mathit{\epsilon }}_{\mathrm{sed}}{)}^{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{10}}^{\mathrm{2}}){\mathit{\epsilon }}_{\mathrm{sed}}}{a{\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}+a(c-{\mathit{\epsilon }}_{\mathrm{sed}}{)}^{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{10}}^{\mathrm{2}}}}.\end{array}$$

The difference between ${\mathit{\epsilon }}_{\mathrm{sed}}^{\mathrm{eff}}$ and in situ ε_sed is shown in Fig. 4b to shift ε_Nd values by up to about ±0.03 . In other words, for in situ ε_Nd values lower than the minimum of α (dashed gray line in Fig. 4a and b), the released ε_Nd value is pushed toward even lower values, and for in situ ε_Nd values greater than the minimum of α, released ε_Nd is pushed toward even larger values. (The derivation of Eq. 10 is given in Appendix C). Applying ${\mathit{R}}_{\mathrm{sed}}^{\mathrm{eff}}=R_{\mathrm{CHUR}}({\mathit{\epsilon }}_{\mathrm{sed}}^{\mathrm{eff}}+\mathrm{1})$ gives the sedimentary source of RNd as ${\mathit{R}}_{\mathrm{sed}}^{\mathrm{eff}}{\mathit{s}}_{\mathrm{sed}}$.

2.2.4 Rivers

For riverine sources, we use the Global River Flow and Continental Discharge Dataset (Dai, 2017), originally described by Dai and Trenberth (2002) and later updated by Dai et al. (2009) and Dai (2016). This dataset provides an estimate of the volumetric flow rate of the 200 largest rivers on Earth. As a simplification, and to reduce the total number of free parameters in the model, we assume that all rivers share the same Nd concentration c_river, which is the parameter that controls the global riverine source magnitude (see Table 2). (Future improvements of the GNOM could include optimizable [Nd] parameters for each individual major river, constrained by ranges based on observations.) Because the GNOM v1.0 does not resolve estuary removal processes, our c_river is to be understood as an effective Nd concentration that implicitly accounts for Nd removal in estuaries and is thus the concentration that makes it into the ocean. Hence, the vector of the riverine source is given by the following:

$$\begin{array}{}\text{(11)}& {\mathit{s}}_{\mathrm{river}}=c_{\mathrm{river}}{\displaystyle \frac{{\mathit{Q}}_{\mathrm{river}}}{\mathit{v}}},\end{array}$$

where Q_river is the vector of the volumetric flow rates of the rivers from the dataset of Dai (2017) gridded (cumulatively) onto the OCIM v2.0 grid, and v is the vector of the volumes of the model grid cells. Note that in order to prevent numerical noise from the large gradients caused by the discrete nature of their distribution, we additionally artificially spatially smooth out the riverine sources by spreading it over neighboring grid boxes (see Fig. 2d).

Riverine ε_Nd values are taken from the global map of interpolated sedimentary ε_Nd by Robinson et al. (2021), which we also use for the sedimentary source, so that the RNd riverine source is given by ${\mathit{R}}_{\mathrm{sed}}^{\mathrm{eff}}{\mathit{s}}_{\mathrm{river}}$. We note that the sedimentary ε_Nd map of Robinson et al. (2021) overlays the nearest continental ε_Nd signal where sediment thickness is more than 1 km, such that the ε_Nd of the GNOM v1.0 riverine sources are mostly from continental measurements that lie within or close to the river drainage basins. Riverine ε_Nd values are shown in Fig. 2j.

2.2.5 Groundwater

Neodymium also enters the oceans via coastal groundwater (Johannesson and Burdige, 2007). We use the coastal submarine and terrestrial groundwater discharge dataset of Luijendijk et al. (2019), described by Luijendijk et al. (2020), which provides the location and volumetric flow rate of 40 082 coastal watersheds. Similar to the riverine sources, we assume that [Nd] is constant across river watersheds and implicitly accounts for local Nd removal processes. The single parameter c_gw is thus the effective groundwater concentration that makes it into the ocean and controls the global magnitude of the GNOM v1.0 groundwater Nd source (see Table 2). The groundwater Nd source is given by the following:

$$\begin{array}{}\text{(12)}& {\mathit{s}}_{\mathrm{gw}}=c_{\mathrm{gw}}{\displaystyle \frac{{\mathit{Q}}_{\mathrm{gw}}}{\mathit{v}}},\end{array}$$

where Q_gw is the groundwater volumetric flow rates from the dataset of Luijendijk et al. (2020), gridded cumulatively onto the GNOM grid. The pattern of s_gw is shown in Fig. 2e.

Following Jeandel et al. (2007), we assume that the ε_Nd of Nd released through groundwater is determined by the local lithogenic isotopic composition. However, instead of the dataset of Jeandel et al. (2007), we use the more recent Robinson et al. (2021) data, exactly like for the riverine source. These ε_Nd values, which are located near the coast, are likely adequately representing the local lithogenic composition because Robinson et al. (2021) assign continental values where sediment thickness is greater than 1 km. The groundwater ε_Nd values are shown on Fig. 2k.

2.2.6 Hydrothermal vents

A minor fraction of the marine neodymium budget presumably comes from hydrothermal vents, which deliver likely young Nd (high ε_Nd) along the mid-ocean ridges (Piepgras and Wasserburg, 1985; Stichel et al., 2018). Here, we assume that the release of hydrothermal Nd is proportional to that of helium. For consistency, we use the mantle helium source field that was used in the data assimilation of the OCIM v2.0 (DeVries and Holzer, 2019). The global magnitude and ε_Nd of the hydrothermal Nd source are set by parameters σ_hydro and ε_hydro, respectively (see Table 2), with the following:

$$\begin{array}{}\text{(13)}& {\mathit{s}}_{\mathrm{hydro}}={\mathit{\sigma }}_{\mathrm{hydro}}{\displaystyle \frac{{\mathit{s}}_{\mathrm{He}}}{{\mathit{v}}^{\mathsf{T}}{\mathit{s}}_{\mathrm{He}}}},\end{array}$$

where s_He is the vector of the ³He mantle source, and v^𝖳s_He is its global magnitude, i.e., its volume integral, used here for normalization. (One can easily check that v^𝖳s_hydro=σ_hydro.) The hydrothermal source of RNd is simply given by ${\mathit{R}}_{\mathrm{hydro}}{\mathit{s}}_{\mathrm{hydro}}=R_{\mathrm{CHUR}}(\mathrm{1}+{\mathit{\epsilon }}_{\mathrm{hydro}}){\mathit{s}}_{\mathrm{hyrdo}}$. Figure 2f and l show the spatial distribution of the hydrothermal Nd source and its ε_Nd, respectively.

Arguably, the hydrothermal system as a whole acts as a net sink of Nd in the ocean (Stichel et al., 2018). As described in Sect. 2.3, the GNOM v1.0 does not include a parameterization of scavenging due to hydrothermal particles. Future versions of the GNOM should attempt to include such a removal process in order to properly balance the hydrothermal source and allow the ε_Nd signature to be modified along hydrothermal vents without increasing the [Nd] concentration at the same time.

2.3 Reversible scavenging

Neodymium is removed from the system through scavenging onto particles. We follow, e.g., Bacon and Anderson (1982), Siddall et al. (2008), and Arsouze et al. (2009) and assume that dissolved and scavenged Nd are exchanged via a first-order kinetic reaction, as follows:

$$\begin{array}{}\text{(R1)}& \mathrm{Nd}+X\rightleftharpoons X\mathrm{Nd},\end{array}$$

where X is a given particle type. The rate of change of [Nd] in Reaction (R1) can be written as follows:

$$\begin{array}{}\text{(14)}& -{\displaystyle \frac{\partial \left[\mathrm{Nd}\right]}{\partial t}}={\displaystyle \frac{\partial \left[X\mathrm{Nd}\right]}{\partial t}}=k_X^{+}\left[\mathrm{Nd}\right]\left[X\right]-k_X^{-}\left[X\mathrm{Nd}\right],\end{array}$$

with the following equilibrium constant:

$$\begin{array}{}\text{(15)}& K_X={\displaystyle \frac{k_X^{+}}{k_X^{-}}}={\displaystyle \frac{\left[X\mathrm{Nd}\right]}{\left[\mathrm{Nd}\right]\left[X\right]}}.\end{array}$$

We further assume that each scavenging particle type X has a constant settling velocity w_X that dominates the transport rates of the ocean circulation. For each particle type X, we construct its flux divergence operator, denoted T_X, such that T_X x is the discrete equivalent of ∇⋅(w_X [X]) (where x is the particulate concentration vector, and w_X is the downward three-dimensional settling velocity vector). We use T_X to compute the rate at which reversible scavenging adds or removes Nd in each grid box.

However, we avoid the explicit simulation of scavenged neodymium, XNd, by having a fraction of the dissolved neodymium pool sink to the box below as if it were adsorbed onto a falling particle. To do this in practice, we take advantage of the direct relationship between free and scavenged Nd, Eq. (15), assuming that Reaction (R1) operates on shorter timescales than either vertical particulate transport or ocean transport. (This assumption is common in models that include scavenging and simpler than resolving the adsorption/desorption rates dynamically (e.g., van Hulten et al., 2018).) Since dissolved and scavenged Nd are in equilibrium, Eq. (15) uniquely determines [XNd]=K_X [X] [Nd], given the modeled [Nd] and the prescribed particle concentration [X] (from the four particle fields included in the GNOM v1.0, described below). Consequently, the corresponding partial downward flux of dissolved Nd is given by w_X [XNd], where w_X is the settling velocity of particle X. We further assume that a fraction f_X of the scavenged Nd that reaches the seafloor is removed from the system, providing a net sink for our model. (Note that this is the same implicit approach as in the AWESOME OCIM (John et al., 2020).)

We consider four different particle types for scavenging Nd. (i) Scavenging by dust particles is modeled using the dust deposition fields of Kok et al. (2021b), assuming a vertically constant concentration as dust particles settle with velocity w_dust=1 km yr⁻¹ through the water column. (ii) Scavenging by particulate organic carbon (POC) is modeled using the three-dimensional POC concentration field from the work of Weber et al. (2018) and following the AWESOME OCIM implementation of John et al. (2020), with a settling velocity w_POC=40 m d⁻¹. (iii) Scavenging by biogenic silica (bSi), or opal, is modeled using a simple, nutrient-restoring offline model of Si cycling described in Appendix A. (iv) A particle-independent scavenging is included to prevent the accumulation of Nd where the concentration fields of dust, POC, and opal are unrealistically low.

This fourth scavenging mechanism effectively behaves like spontaneous precipitation, and, as such, will be referred to as precipitation throughout this study (subscript prec). Precipitation is implemented by using a spatially uniform fictitious particle concentration of 1 mol m⁻³ that settles with velocity w_prec=0.7 km yr⁻¹. We note that, while this additional particle-independent scavenging sink could compensate for additional types of particles not currently implemented in the model, it is likely that more scavenging particle types are required for an accurate representation of the Nd cycle. These include hydrothermal particles (which should result in hydrothermal systems being a net sink; Stichel et al., 2018) and iron–manganese oxides (which are potentially the most important scavenging particles; Lagarde et al., 2020; Schijf et al., 2015; Sholkovitz et al., 1994). Overall, the scavenging transport operator is thus defined by summing the flux divergence for all particle types and using Eq. (15), as follows:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\mathbf{T}}_{\mathrm{scav}}& =K_{\mathrm{dust}}{\mathbf{T}}_{\mathrm{dust}}{\mathbf{D}}_{\mathrm{dust}}+K_{\mathrm{POC}}{\mathbf{T}}_{\mathrm{POC}}{\mathbf{D}}_{\mathrm{POC}}\\ & +K_{\mathrm{bSi}}{\mathbf{T}}_{\mathrm{bSi}}{\mathbf{D}}_{\mathrm{bSi}}+K_{\mathrm{prec}}{\mathbf{T}}_{\mathrm{prec}},\end{array}\end{array}$$

where D_X is a diagonal matrix with diagonal x, which is the vector of the concentrations of particle type X. Hence, for each scavenging particle type X, the corresponding scavenging rate and downward transport is controlled by the concentration [X], the equilibrium constant K_X, the settling velocity w_X, and the burial fraction f_X. Figure 5 shows maps and profiles of net scavenging rates for each particle type.

Figure 5 (a–d) Vertically integrated and (e–f) horizontally integrated scavenging by precipitation, dust, POC, and opal particles, where positive values represent Nd removal and negative values represent a source of pumped Nd. Note that in panels (e)–(f), net scavenging rates in the top layer are positive and largest by construction because Nd can only be removed there, as opposed to all the layers below which can receive Nd from superjacent layers.

2.4 Optimization

The output of our model is governed by a set of 43 free parameters that control the magnitude and ε_Nd of each source and the reversible scavenging and burial rates. The computational speed afforded by the model implementation allows us to jointly optimize almost all of these parameters by minimizing the mismatch of modeled and observed [Nd] and ε_Nd values. This is done in practice by minimizing an objective function $\widehat{f}\left(\mathit{p}\right)$ that quantifies the mismatch between model and observations, for a given set of parameters p.

2.4.1 Objective function

The mismatch with each observation is quantified by the square of the difference between the observed value and the modeled value from the closest grid cell. Because we use observations of [Nd] and ε_Nd, the mismatch function, denoted f, depends on the 3D fields of the two modeled tracers (${\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}$ and ${\mathit{\chi }}_{\mathrm{RNd}}^{\mathrm{mod}}$). We also include an additional cost for parameter values themselves. The mismatch function is defined by the following:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}f({\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}},{\mathit{\chi }}_{\mathrm{RNd}}^{\mathrm{mod}},\mathit{p})& ={\mathit{\omega }}_{\mathrm{Nd}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{{\sum }_{\mathit{r}\in {\mathcal{O}}_{\mathrm{Nd}}}{\left({\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}\left[\mathit{r}\right]-{\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{obs}}\left[\mathit{r}\right]\right)}^{\mathrm{2}}}{{\sum }_{\mathit{r}\in {\mathcal{O}}_{\mathrm{Nd}}}{\left({\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{obs}}\left[\mathit{r}\right]\right)}^{\mathrm{2}}}}\\ & +{\mathit{\omega }}_{{\mathit{\epsilon }}_{\mathrm{Nd}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{{\sum }_{\mathit{r}\in {\mathcal{O}}_{{\mathit{\epsilon }}_{\mathrm{Nd}}}}{\left({\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{mod}}\left[\mathit{r}\right]-{\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{obs}}\left[\mathit{r}\right]\right)}^{\mathrm{2}}}{{\sum }_{\mathit{r}\in {\mathcal{O}}_{{\mathit{\epsilon }}_{\mathrm{Nd}}}}{\left({\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{obs}}\left[\mathit{r}\right]\right)}^{\mathrm{2}}}}\\ & -{\mathit{\omega }}_{\mathit{p}}\sum _i\log \left(P\right({\mathit{p}}_i\left)\right),\end{array}\end{array}$$

where the first term represents the normalized mismatch between modeled and observed [Nd], the second term represents the normalized mismatch between modeled and observed ε_Nd, and the last term represents the inverse of the likelihood of the model parameters. We detail each term below.

In Eq. (17), ${\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{obs}}$ is a vector of all the [Nd] observations, and ${\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{obs}}\left[\mathit{r}\right]$ denotes the observed [Nd] at location r, which spans all the locations of observations, 𝒪_Nd. (One can think of r as indexing the vector of observations ${\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{obs}}$.) We compare each observation with the model output from the closest model grid cell, denoted ${\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}\left[\mathit{r}\right]$ for simplicity. (Technically, the location of the observed and modeled value being compared may not match, in which case we use the closest wet model grid cell using a nearest-neighbor algorithm.) We use the same approach for ε_Nd by comparing the modeled vector ${\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{mod}}$ to the observed vector ${\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{obs}}$ at the locations $\mathit{r}\in {\mathcal{O}}_{{\mathit{\epsilon }}_{\mathrm{Nd}}}$ of each ε_Nd observation. The ω_Nd and ${\mathit{\omega }}_{{\mathit{\epsilon }}_{\mathrm{Nd}}}$ values, fixed at 1 in this study, control the relative contributions of the mismatches in Nd concentrations and ε_Nd values. Given these, an error of 1  in ε_Nd weighs the same as an error of about 4.5 pM.

The third term of Eq. (17) adds a direct penalty constraint on the parameters to prevent them from reaching unrealistic values. If any parameter reaches a value close to the limits we impose in the model, then this third term will grow large; since the algorithm tries to minimize Eq. (17), it will push that parameter back to a more acceptable value. For each parameter p_i, we prescribe a realistic domain D_i and an initial guess (see Table 2) that we use to determine a reasonable prior probability distribution d_i and to randomize the initial parameter values. Specifically, each parameter with a semi-infinite range (from 0 to ∞) is given a lognormal prior of which the logarithm has mean equal to the logarithm of the initial guess and has variance equal to 1. Each parameter with a finite range is given a logit-normal prior that is scaled and shifted, such that its support matches the range D_i exactly, such that the initial guess equals the median of d_i. For example, the ϕ₀ parameter for the sedimentary flux at the surface is given the (0,∞) range and an initial guess of 20 $\mathrm{pmol}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. Taken as a random variable, the prior distribution of ϕ₀ is the lognormal distribution, such that $\log ({\mathit{\varphi }}_{\mathrm{0}}/(\mathrm{pmol}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}\left)\right)\sim \mathcal{N}\left(\log \right(\mathrm{20}),\mathrm{1})$. (𝒩(μ,σ²) denotes the normal distribution with mean μ and variance σ².)

For the performance and robustness of the optimization, we additionally perform a variable transform λ_i on each parameter using a bijection from the parameter domain D_i to $(-\mathrm{\infty },\mathrm{\infty })$. This variable transform prevents parameters from reaching beyond their prescribed ranges. We also carefully chose the bijection, such that the prior distribution is normally distributed in the transformed parameter space and thus incurs an inverse log-likelihood that is quadratic, a property that benefits the performance of the optimization. In the case of the parameter ϕ₀, which is transformed via the bijection ${\mathit{\lambda }}_i:{\mathit{\varphi }}_{\mathrm{0}}\mapsto \log ({\mathit{\varphi }}_{\mathrm{0}}/(\mathrm{pmol}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}\left)\right)$, the corresponding transformed random variable is normally distributed by construction. For bounded parameters, such as β_volc, a shifted and scaled logit transform is applied, which also yields a transformed random variable that is normally distributed by construction.

The ω_p value, fixed at 10⁻⁴ in this study, controls the relative size of the penalty for the parameters compared to the cost of Nd and ε_Nd. It is chosen so that the [Nd] and ε_Nd mismatch costs are generally about 2 to 3 orders of magnitude larger than the parameter penalty (although there is no bound on the parameter penalty for extreme parameter values). The primary role of this added parameter cost and the associated variable transform is to improve the convergence rate of the optimization and help prevent it from becoming stuck in valleys of parameter space (see, e.g., Nocedal and Wright, 1999).

The objective function depends on p only and is defined by the following:

$$\begin{array}{}\text{(18)}& \widehat{f}\left(\mathit{p}\right)\equiv f\left({\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}\right(\mathit{p}),{\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{mod}}(\mathit{p}),\mathit{p}),\end{array}$$

where we have explicitly marked ${\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}\left(\mathit{p}\right)$ and ${\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{mod}}\left(\mathit{p}\right)$ as functions of the parameters p. That is, for any choice of parameters p, before evaluating the model mismatch as quantified by the objective function, we must first compute the vectors ${\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}}$ and ${\mathit{\chi }}_{\mathrm{RNd}}^{\mathrm{mod}}$ by solving for the steady-state solution to Eq. (2). The gradient, $\mathrm{\nabla }\widehat{f}$, and Hessian, ${\mathrm{\nabla }}^{\mathrm{2}}\widehat{f}$, of the objective function are computed using a combination of autodifferentiation and adjoint techniques available from within the AIBECS.jl package (Pasquier, 2020a; Pasquier et al., 2022b) or specifically developed in parallel for computational efficiency (F1Method.jl; Pasquier, 2020b).

2.4.2 Dissolved neodymium and ε_Nd data

The [Nd] and ε_Nd observations used in this study consist of the following three datasets: (i) the pre-GEOTRACES compilation of Nd and ε_Nd data by van de Flierdt et al. (2016a), (ii) the GEOTRACES Intermediate Data Product 2017 (IDP17; Schlitzer et al., 2018) (including specifically Nd-linked publications; Stichel et al., 2012a, b, 2015; Garcia-Solsona et al., 2014; Basak et al., 2015; Fröllje et al., 2016; Lambelet et al., 2016, 2018; Behrens et al., 2018a, b), and (iii) our post-IDP17 compilation of data from the Indian Ocean (Amakawa et al., 2019), the Barents Sea (Laukert et al., 2018, 2019), the northern Iceland Basin (Morrison et al., 2019), the northwestern Pacific (Che and Zhang, 2018), the Kerguelen Plateau (Grenier et al., 2018), the southeastern Atlantic Ocean (GA08, Rahlf et al., 2020, 2019, 2021), the Bay of Biscay (Dausmann et al., 2020, 2019), the western North Atlantic (Stichel et al., 2020), the Arctic (Laukert et al., 2017a, d), and the Bermuda Atlantic Time-series Study (BATS; Laukert et al., 2017b, c). The spatial distribution of these observations as used in this study are shown in Fig. 6.

Figure 6 Locations of (a) Nd and (b) ε_Nd observations used in this study, labeled per the data source.

2.4.3 Minimization algorithm

We use the Newton trust region algorithm from the Optim.jl package (Mogensen and Riseth, 2018) to minimize the objective function $\widehat{f}\left(\mathit{p}\right)$. This requires, at every iteration, the objective function, its gradient, and its Hessian, which are evaluated using the F1Method.jl and AIBECS.jl packages (Pasquier, 2020a; Pasquier et al., 2022b; Pasquier, 2020b).

Thanks to the computationally efficient gradient optimization algorithm that leverages gradient and Hessian information, the entire optimization run takes a few hours on a modern laptop. In our experience, for comparison, using the more standard finite differences approach or an optimization algorithm that does not have access to derivatives would likely take multiple months.

Because the Newton trust region algorithm performs local rather than global optimization, we run multiple optimizations, starting from randomized initial parameter values sampled from the parameter distributions d_i (as shown in Fig. B1). Although not all optimization runs end up in the same state because of many local minima, we find that most of them converge towards similar solutions with a similar small objective function value, which we denote as our best estimate, and out of which all the figures in this work are created. Figure B1 also shows the initial values and final values of a dozen of optimization runs, with the best estimate shown in blue.

3 Results

3.1 Parameter values

Our best estimate of the state of the Nd cycle is given by the set of parameters that minimizes the objective function defined in Eqs. (17)–(18). We emphasize that our estimate is determined by a local minimum of a specific parameterization, such that the best here is somewhat subjective. In all likelihood, there exist other models and other parameter choices which produce a similar fit to global observations, though we expect all such models to capture the same key features of global Nd biogeochemical cycling. Initial guesses and optimizable ranges for each parameter were determined from the literature and the expertise of authors. Initial guesses and final parameter values, along with unit, prescribed range, and a brief description, are given in Table 2. Parameter prior distributions, their randomized initial values, and final optimized values are shown in Fig. B1.

Table 2 List of parameters. Realistic parameter ranges were prescribed based on the literature and the expertise of the authors. Final values have been rounded to three significant digits. Parameters without a range are not optimized (final value equals initial guess). Scavenging reaction constants, K_X, are reported in units of inverse concentration of the particle X.

In Table 2, parameters without a range indicate that they were not optimized and held fixed at the given previous model or literature values. For example, in the case of scavenging by each particle type X, we only optimized K_X and f_X (not w_X). As described in Sect. 2.3 above, the settling velocities for POC and bSi are not optimized and are instead fixed to match the values of their respective parent offline models. While there are no parent models for dust and precipitation, we do not optimize the corresponding settling velocities for these particle types either because K_X and w_X can perfectly compensate each other. For example, doubling K_X while halving w_X has no effect on Nd distributions and the objective function. Only their product, K_X w_X, which sets the strength of the scavenging pump through the operator matrix T_scav, appears in the tracer equations (see Eq. 16 or, e.g., John et al., 2020), such that these parameters cannot be easily optimized independently.

3.2 Fit to observations

The general fit to observations is illustrated in Fig. 7, which shows the percentiles of the cumulative joint probability distribution of the modeled and observed Nd concentrations and ε_Nd values. Despite the slightly visible spread, most of the modeled–observed [Nd] and ε_Nd values lie close to the 1:1 line, indicating a good match, with a root mean square error of about 6.80 pM and 2.09 , respectively.

Figure 7 Quantiles of the cumulative joint probability density functions of modeled and observed (a) Nd concentrations and (b) ε_Nd values. Darker colors indicate a high density of data, such that n % of the modeled and observed data lie outside of the nth percentile contour. The closer the darker contours are to the 1:1 black dashed line, the better the fit.

While statistics such as Fig. 7 provide important information at a quick glance, they do not retain any geographical information, so that a more detailed investigation is required to fully assess the model's skill. Indeed, the deviations shown by [Nd] and ε_Nd clusters slightly off the 1:1 line (Fig. 7) likely reflect groups of geographically proximate data points that may be symptomatic of systematic biases, which must be analyzed in further detail.

We explore the regional variations in the model's skill with depth in Fig. 8, which shows the basin-averaged profiles of modeled and observed [Nd] and ε_Nd for the Atlantic, Pacific, Indian, and Southern oceans. Simulated [Nd] fits the nutrient-like profiles of basin mean observations and captures the bulk of interbasin variance fairly well despite systematic biases of about −3 pM in the mid-depth Atlantic, +3 pM in the deep Atlantic, up to +6 pM in the mid-depth Indian Ocean, and −6 pM in the 4,250 m deep Southern Ocean (Fig. 8a–d). Similarly, for ε_Nd values, we find an overestimate of about +1  below 700 m in the Atlantic and an underestimate of up to −2  in the Pacific, particularly near the surface, while the modeled and observed basin-averaged Southern Ocean profiles are a tight fit (Fig. 8e–h).

Figure 8 Basin-averaged profiles of (a–d) Nd concentrations and (e–f) ε_Nd values versus depth. The basins (Atlantic, Pacific, Indian, and Southern Ocean) with the number of observations for each tracer are reported in the bottom right corner of each panel. The mean and standard deviation of observations are calculated at each vertical grid level of the OCIM v2.0 grid and represented by the thin line and error bars. The mean and standard deviation of the model are represented by the thick line and lighter-colored ribbon.

We further assess the model skill by looking at GEOTRACES transects individually. Out of the 3483 observations of [Nd] that we use to constrain our model, 1575 (∼45 %) come from the IDP17 (Schlitzer et al., 2018) and were collected along the GA02, GA03, GA10, GA11, GAc01, GIPY04, GIPY05, GIPY06, GPc02, GPpr04, and GPpr05 cruises (Fig. 9b–l). Similarly, out of the 2988 ε_Nd observations, 790 (∼26 %) come from IDP17 cruises GA02, GA03, GA11, GIPY04, GIPY05, GIPY06, GPc02, GPpr04, and GPpr05 (Fig. 10b–j).

Figure 9 reveals in detail how well GNOM output matches observational [Nd] data. The model captures the broad interbasin and intrabasin variations with high fidelity despite some slight mismatches. Specifically, Fig. 9c reveals overestimates of mid-depth and deep [Nd] in the west Pacific (cruise GPpr04) with an overestimated gradient with depth, potentially due to too large a sedimentary source or too strong scavenging. Mid-depth overestimates of [Nd] also appear in the Atlantic (GA03, GAc01, and GA10; Fig. 9h, j, k). However, the deepest [Nd] are underestimated in association with a generally underestimated vertical gradient, particularly along GA10. Hence, the [Nd] mismatches in the Pacific are suggestive of either too weak a deep sedimentary source or too efficient scavenging and burial in the deep. These systematically opposed mismatches between the Atlantic and Pacific are likely due to the optimization procedure, which balances out all the mismatches simultaneously. Future improvements in the GNOM that resolve these discrepancies could include different parameterizations of the sedimentary source and scavenging or the addition of a currently absent third mechanism, such as a nepheloid source or sink.

Figure 9 (a) GEOTRACES cruise tracks with [Nd] observations. The legend layout matches the layout of the other subplots of the figure. (b–l) Modeled and observed [Nd] along GEOTRACES transects. Modeled values are shown as filled contours, while observed values are overlaid as a scatterplot.

Figure 10 shows modeled and observed ε_Nd values along IDP17 cruise transects. While the observed interbasin variability is adequately represented, the GNOM does not perfectly capture the finer spatial details of observed ε_Nd, suggesting that there is still room for model improvement. The model agrees well with observations along Southern Ocean transects (Fig. 10e, i, j; GPc02, GIPY04, and GIPY05). However, in the North Atlantic (e.g., Fig. 10f, GA02), the GNOM does not entirely capture a strong negative ε_Nd plume along the North Atlantic Deep Water (NADW). Conversely, in the western Pacific (Fig. 10c; GPpr04), our model misses strongly positive surface ε_Nd observations and instead displays a deep plume of positive ε_Nd that is absent from the observational data. This is potentially due to missing mechanisms, sources, or sinks, the correct implementation of which would likely benefit from more observational ε_Nd data in the Pacific and Indian basins.

Figure 10 (a) GEOTRACES cruise tracks with ε_Nd observations. The legend layout matches the layout of the other subplots of the figure. (b–j) Modeled and observed ε_Nd along GEOTRACES transects. Modeled values are shown as filled contours, while observed values are overlaid as a scatterplot.

Figure 11 shows the model and observed [Nd] and ε_Nd for all observations averaged over different depth ranges. Contrary to Figs. 9 and 10, this includes all observations used to constrain the GNOM v1.0 (i.e., not just IDP17). As expected, the model broadly matches the observational data well, with some systematic mismatches in different locations. Figure 11d–e reveal an underestimate of deep [Nd] in the northern Indian Ocean in the Bay of Bengal, which is likely attributable to too strong scavenging or too weak sedimentary fluxes into the deeper layers of the model. Fig. 11c shows elevated Nd concentrations in the deepest parts of Baffin Bay, potentially due to too large sources, lack of data, or even circulation issues related to the resolution of the model grid in that region. Notably, Fig. 11c–e reveal discrepancies among observations, with a few [Nd] values near the GA02 transect that stand out compared to neighboring observations. Figure 11f–i show a substantial underestimate of equatorial Pacific ε_Nd values above 1500 m depth, from east Indonesia to the west coast of Ecuador and Peru. Strongly negative observed ε_Nd values in the northwestern Atlantic and near the west coast of Africa from the Congo River to Namibia are not well captured by the model. Conversely, rather positive ε_Nd values at the surface going north along the western coast of North Africa also seem to evade the capability of GNOM to represent observations. The lack of resolution or some important missing mechanisms are likely the cause of these larger mismatches. In future versions of GNOM, we intend to improve the model by targeting these regions of particularly pronounced misfits.

Figure 11 Model (heatmaps) and observations (markers) for (a–f) Nd concentrations and (g–k) ε_Nd values. Model values are averaged over the indicated depth range. Individual observed values are overlaid on top of the modeled heatmap and on top of each other, so that perfect model–observation matches and deeper observations can be hard to see or hidden.

3.3 Diagnostics

One of the biggest advances in the GNOM v1.0, compared to earlier models of the marine Nd cycle, is due to the steady-state matrix formulation of the model, which allows us to compute advanced and detailed diagnostics that can directly address fundamental questions about the distribution of tracers and better understand their cycle. In the following sections, we showcase a few such diagnostics.

3.3.1 Source magnitudes

The optimized parameters determine the magnitude of the sources, which are collected in Table 3. In our best estimate, about 66 Mmol of Nd (or about 9500 metric tonnes) are injected into the global ocean every year. This falls slightly above the 38–57 Mmol yr⁻¹ range of previous GCM models (see Table 1 for model references).

The eolian dust and sedimentary sources are the dominant ones contributing 24 and 32 Mmol yr⁻¹ (about 35 % and 50 % of to total source), respectively. While this falls within the 0.69–60 Mmol yr⁻¹ range for global eolian source magnitudes of previous modeling studies, our eolian sources are an order of magnitude larger than previous GCM-based modeling studies (0.69–3.5 Mmol yr⁻¹; Table 1) and than the 4.4 Mmol yr⁻¹ estimate of Greaves et al. (1994). This is likely due to our optimization procedure, during which Nd solubility is allowed to be adjusted within the whole 0 %–100 % range, compared to previous GCM-based studies that typically use a fixed 2 % solubility (Arsouze et al., 2009; Gu et al., 2019; Pöppelmeier et al., 2020). This is despite the initial guesses for β_r values for dust set at 5 %, which penalizes large solubilities more than low solubilities (see the increased probability densities for low solubilities in Figs. B1a and B1b). The high optimized solubility of volcanic ash β_volc is also likely unrealistic, although the total contribution of volcanic ash is much smaller than mineral dust. We note that generally worse fits to [Nd] and ε_Nd observations have been achieved with some of our optimization runs ending in distinct local minima with significantly smaller dust solubilities. (We do not show these worse mismatches, but we show the corresponding initial and final parameter values in Fig. B1.) Finally, we emphasize that it is not the goal of this work to establish estimates of the GNOM parameters and that we welcome future GNOM users to apply narrower ranges for those parameters for which they have better constraints (for example, restricting dust Nd solubilities to values below 10 %).

At 32 Mmol yr⁻¹, the GNOM sedimentary source falls right within the 0–78 Mmol yr⁻¹ range of previous models (28–55 Mmol yr⁻¹ for GCM-based studies; Table 1) and in agreement with the 18–110 Mmol yr⁻¹ range of Abbott et al. (2015b). The third-largest source is riverine, with about 10 Mmol yr⁻¹, also in accord with the published 1.8–12.4 Mmol yr⁻¹ range in previous models (Table 1) and similar to the 4.6–12 Mmol yr⁻¹ values from Goldstein and Jacobsen (1987) and Greaves et al. (1994). The GNOM optimization did not favor a large source from submarine groundwater discharge, which has been estimated between 29–81 Mmol yr⁻¹ by Johannesson and Burdige (2007), but this source was not included in any other modeling studies, so it is difficult to compare with other global model estimates. Hence, apart from a relatively elevated dust source, the GNOM optimization generally supports previous estimates of the magnitude of these sources.

Table 3 Source magnitudes, their Nd contributions, and corresponding bulk residence time. Values are reported with two significant digits.

3.3.2 Partition according to source type

We partition Nd concentration according to source type simply by removing all the other sources. (The superposition principle applies directly to our model because it is linear in Nd; see, e.g., Holzer et al., 2016, who partitioned dissolved iron concentrations, first requiring the construction of a linear equivalent model.). With Nd_k denoting neodymium that was injected by source type k, its corresponding column vector is thus computed by solving the following:

$$\begin{array}{}\text{(19)}& \mathbf{H}{\mathit{\chi }}_{{\mathrm{Nd}}_k}^{\mathrm{mod}}={\mathit{s}}_k,\end{array}$$

where $\mathbf{H}={\mathbf{T}}_{\mathrm{circ}}+{\mathbf{T}}_{\mathrm{scav}}$. Taking the global volume-weighted mean of each ${\mathit{\chi }}_{{\mathrm{Nd}}_k}^{\mathrm{mod}}$ gives the contribution of each source type to the total Nd inventory and are collected in Table 3.

As is the case for most global biogeochemical cycles, the relative source magnitudes and their relative contribution to the standing stock do not necessarily match. For instance, mineral dust, volcanic ash, and hydrothermal vents contribute more to the mean Nd concentration than their relative source magnitudes. These variations can be directly linked to the bulk residence time of Nd molecules, which vary with location of injection and consequently with source type.

3.3.3 Bulk residence times

The bulk residence time of Nd_k is given by taking the ratio of its inventory to its source magnitude. Total Nd (i.e., from all sources) has a bulk residence time of roughly 440 years, which is within the 350–2900 year range of previous Nd cycling models. (The bulk residence time for GNOM falls slightly below the 560–750 year range of GCM-based models; see Table 1).

Unlike in the real ocean, where each molecule of Nd is indistinguishable from the next (with no information about its initial source), in a model we can track Nd coming from different sources and calculate source-partitioned residence times. We find that sedimentary-sourced Nd has the shortest residence time at 370 years. This is because it is injected just above the seafloor and is thus buried more quickly (i.e., a molecule of Nd sourced from the sediments at 5000 m, which only has to fall a few meters to be scavenged back to the sediments, leaves the ocean quicker than a molecule near the surface sourced from dust, which has to fall thousands of meters). In comparison, volcanic ash Nd, most of which is deposited onto the surface of the Pacific, remains on average 650 years in the system (i.e., about 75 % longer than sedimentary-sourced Nd). Mineral dust deposited Nd and riverine and surface groundwater Nd all show a residence time of about 550 years.

3.3.4 Tracking Nd to investigate ε_Nd conservativeness

Sediment core records of ε_Nd are of considerable importance for paleoceanography because they serve as a fingerprint of past ocean circulation and have been used, in particular, to infer changes in the Atlantic meridional overturning circulation (AMOC; e.g., Palmer and Elderfield, 1985; Rutberg et al., 2000; Piotrowski et al., 2004, 2005; van de Flierdt et al., 2006; Piotrowski et al., 2008; Pena et al., 2013; Pena and Goldstein, 2014; Huang et al., 2014; Kim et al., 2021; Pöppelmeier et al., 2021; Hines et al., 2022). Observations of modern ocean ε_Nd have thus been extensively used to estimate Atlantic water mass fractions, typically inferred from north–south end-members and assuming quasi-conservative transport and mixing of ε_Nd as an ocean circulation tracer (Hartman, 2015; Wu, 2019; Wu et al., 2022).

Our GNOM model – or, more precisely, the steady-state matrix schema with which it is built – allows for exact computations of these water masses and the contributions of various sources and regions to the modern ocean Nd and ε_Nd distributions. We emphasize that by exact, here, we do not mean exact for the real ocean, but instead for our given choice of ocean circulation model (in this case the OCIM v2.0; DeVries and Holzer, 2019), which is arguably the best steady-state ocean circulation model available for such climatological estimates (see, e.g., John et al., 2020). Here, we merely showcase these diagnostics within GNOM, but in the future we plan to further explore the underlying scientific questions that these diagnostics can address.

We chose two simple regions that cover the entire ocean, except for the central Atlantic between 30^∘ S and 40^∘ N. We denote these regions by Ω_N and Ω_S, such that Ω_N borders the northern Atlantic and Ω_S the southern Atlantic. The regions are shown in Fig. 12.

Figure 12 Northern Atlantic (Ω_N; blue) and southern Atlantic (Ω_S; orange) regions used for the ε_Nd conservativeness and the water-tagged Nd diagnostics within the central Atlantic region (from 30^∘ S to 40^∘ N; light gray).

First, we track Nd concentrations from each of these regions. Neodymium concentrations are not conservative, in part due to reversible scavenging and in part due to external sources that inject new Nd along transport pathways. For example, we can track Nd that came from region Ω_N by tagging Nd that comes into contact with Ω_N and removing that tag when Nd enters Ω_S. Mathematically, we can perform this tagging/untagging by simulating a fictitious tracer, denoted Nd_N−tag, for which we enforce a concentration equal to simulated [Nd] in region Ω_N, a concentration equal to zero in Ω_S, and allowing reversible scavenging and burial to remove Nd_N−tag in between. In practice, we compute the corresponding column vector ${\mathit{\chi }}_{{\mathrm{Nd}}_{\mathrm{N}-\mathrm{tag}}}$ by solving the following:

$$\begin{array}{}\text{(20)}& \left(\mathbf{H}+{\mathbf{M}}_{\mathrm{N}}+{\mathbf{M}}_{\mathrm{S}}\right){\mathit{\chi }}_{{\mathrm{Nd}}_{\mathrm{N}-\mathrm{tag}}}={\mathbf{M}}_{\mathrm{N}}{\mathit{\chi }}_{\mathrm{Nd}}^{\mathrm{mod}},\end{array}$$

where the M_i are diagonal matrices of which the diagonals are 1 s⁻¹ for indices (i.e., coordinates) within the region Ω_i and 0 s⁻¹ otherwise. Because of the very short timescale of 1 s employed, Eq. (20) effectively enforces that $\left[{\mathrm{Nd}}_{\mathrm{N}-\mathrm{tag}}\right]=\mathrm{0}$ pM in the region Ω_S and that $\left[{\mathrm{Nd}}_{\mathrm{N}-\mathrm{tag}}\right]=\left[\mathrm{Nd}\right]$ in the region Ω_N (where [Nd] denotes the simulated concentration from the GNOM). We track Nd_S−tag in the same way by replacing M_N with M_S on the right-hand side of Eq. (20). The neodymium concentration that neither came from Ω_N nor Ω_S is then simply given by $\left[\mathrm{Nd}\right]-\left[{\mathrm{Nd}}_{\mathrm{N}-\mathrm{tag}}\right]-\left[{\mathrm{Nd}}_{\mathrm{S}-\mathrm{tag}}\right]$.

Figure 13 shows the Atlantic zonal averages of [Nd_N−tag], [Nd_S−tag], and non-tagged [Nd] (i.e., Nd that was injected in the central Atlantic between 30^∘ S and 40^∘ N). In the Ω_N region (north of 40^∘ N), 100 % of Nd is tagged as Nd_N−tag (Fig. 13a). Similarly, 100 % of Nd is Nd_S−tag in Ω_S (Fig. 13b). Figure 13c shows the remaining non-tagged Nd and Fig. 13d combines Fig. 13a–c by showing only the dominant fraction. A clear signal of the influence of surface sources is visible down to about 1500 m, confining the dominance of the N- and S-tagged fractions of Nd to deep waters and to high latitudes, close to their respective tagging regions. In the future, we intend to further explore the extent of the influence of high latitudes on the distribution of Nd.

Gu et al. (2020) performed a north–south end-member partitioning using the CESM model and its POP2 circulation. They quantified water mass fractions using dye injections at the surface and compared them with water mass reconstructions from deep ε_Nd values. Here we present more detailed partitions that have never been estimated in previous modeling studies to our knowledge.

Figure 13 (a) Atlantic zonal average of the fraction of Nd tagged in the Ω_N region (north of 40^∘ N; blue). (b) Same for Nd tagged in the Ω_S region (south of 30^∘ S; orange). (c) Same for Nd tagged injected within the midlatitude Atlantic (between 30^∘ S and 40^∘ N; green). (d) Atlantic zonal average showing only the dominant fraction of Nd coming from either Ω_N or Ω_S or from within 30^∘ S–40^∘ N. Contour lines are shown for each 10 % increment. Tagging regions are shown in Fig. 12.

We now track ε_Nd, as if it were a conservative tracer, from the Ω_N and Ω_S regions. Technically, this is done by tracking water itself as it leaves either region. Water mass fractions estimated with this approach can be used to directly propagate modeled ε_Nd values from the boundaries of the northern and southern Atlantic regions to provide an exact end-member mixing estimate that serves as a reference.

Water mass fractions have been estimated using a Green function boundary propagator in similar model contexts (e.g., Holzer and Hall, 2008; Primeau, 2005; Holzer and Primeau, 2008). They can be used to calculate the conservative transport of any tracer. As an illustrative example, here, we propagate modeled ε_Nd simultaneously from both the Ω_N and Ω_S regions into the central Atlantic. This theoretical conservative value is denoted by ${\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{\Omega }}$ (with Ω as a superscript to denote that its value is entirely determined by the ε_Nd inside the Ω_N and Ω_S regions), and the corresponding column vector ${\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{\Omega }}$ is computed by solving the following:

$$\begin{array}{}\text{(21)}& \left({\mathbf{T}}_{\mathrm{circ}}+{\mathbf{M}}_{\mathrm{N}}+{\mathbf{M}}_{\mathrm{S}}\right){\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{\Omega }}=\left({\mathbf{M}}_{\mathrm{N}}+{\mathbf{M}}_{\mathrm{S}}\right){\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{mod}}.\end{array}$$

Like Eq. (20), Eq. (21) effectively enforces that ${\mathit{\epsilon }}_{\mathrm{Nd}}^{\mathrm{\Omega }}$ matches modeled ε_Nd inside Ω_N and Ω_S but is conservatively propagated by the advective–diffusive transport operator T_circ outside of Ω. (Note the different operators on both the left-hand side and right-hand sides between Eqs. 20 and 21.) Atlantic zonal averages of modeled ε_Nd, conservatively propagated ε_Nd, and their differences (Δ(ε_Nd)) are shown in Fig. 14. Figure 14c, in particular, shows where ε_Nd behaves conservatively (white), and where it does not (green or pink), in our model. In accordance with Fig. 13, ε_Nd is the least conservative close to the surface, where most of the Nd was never in contact with Ω_S or Ω_N and was instead injected in the midlatitude Atlantic. This surface overestimate is most pronounced away from Ω_N and Ω_S, with Δ(ε_Nd) values of up to +10 . Going from 200 to 1000 m depth, we find Δ(ε_Nd) values decreasing from +5 to +1 . Conservative ε_Nd and true ε_Nd remain within 1  of each other below 1000 m depth, despite a slight Δ(ε_Nd) overestimate near the seafloor (likely due to the effect of local sedimentary flux) and a slight underestimate around 1500 m (potentially due to reversible scavenging). We intend to investigate the distinct conservativeness of Nd and ε_Nd (the neodymium paradox) further in future work.

Figure 14 (a) Atlantic zonal average of modeled ε_Nd. (b) Same for conservatively transported ε_Nd from regions Ω_S and Ω_N. (c) Difference between panels (a) and (b). Contour lines are shown for each 1  increment. (Contour lines for $\mathrm{\Delta }\left({\mathit{\epsilon }}_{\mathrm{Nd}}\right)>+\mathrm{5}$  are not shown to avoid clutter). Regions are shown in Fig. 12.

4 Conclusions

The most prominent caveat of GNOM v1.0 is the steady-state assumption, which we apply to both the circulation and the Nd cycle. Hence, by construction, daily, seasonal, decadal, or multi-decadal fluctuations that deviate from the climatological mean cannot be captured by our model. However, we trust that the circulation model used (the OCIM v2.0; DeVries and Holzer, 2019), which is data-assimilated with ventilation tracers, captures the predominant features and pathways of the modern ocean circulation and provides the most realistic steady-state transport to date (e.g., DeVries and Holzer, 2019; John et al., 2020).

We note that, compared to previous modeling studies, the GNOM does not represent scavenging by calcium carbonate (CaCO₃) because there is no publicly available particulate CaCO₃ field, to the best of our knowledge. Modeling scavenging is a challenging task that the GNOM model does not pretend to achieve with high accuracy. However, we deem the current implementation satisfactory, considering the quality of the overall model–observation fit. Future versions of the GNOM could include CaCO₃ particle scavenging or a generally improved scavenging parameterization.

Our model reveals some locations of particular interest for improving our understanding of the Nd cycle and ε_Nd patterns. For instance, there are only two GEOTRACES transects in the Pacific Ocean which cover the western Pacific and the Southern Ocean, the zonal transects in the Atlantic contain only a few ε_Nd measurements compared to [Nd] (e.g., GA03, GA10, GA11, and GAc01), and there are no published transects in the Indian Ocean, which may contribute a non-negligible fraction of Nd to both the Atlantic and Pacific. In the future, we hope that more data will be made available and improve the capabilities of data-constrained model estimates of ε_Nd and the Nd cycle.

While our model endeavors to use formulations and parameter constraints which have reasonable biogeochemical interpretations, there is always room for improvement. For example, our sedimentary source parameterization, which we plan to investigate further in future work, assumes an exponential profile for the sedimentary source flux. While this parameterization is flexible enough to reproduce most of the qualitative features of the sedimentary fluxes in previous models, one might argue that a more mechanistic model of the sedimentary budget would result in a more realistic overall Nd cycling model.

Despite the theoretical advantages they confer to the convergence rate of our optimization procedure, our specific choice of prior distributions for the parameters (Fig. B1) remains arbitrary, though the ranges chosen are of course informed by prior work. Different ranges and initial conditions would yield different optimal solutions and therefore different parameter values. Our choice of a local Newton trust region optimizer also exposes our strategy to the risk of getting stuck in local minima. To counter that risk, we have opted for the traditional ad hoc countermeasure of running the optimization from multiple, randomized initial guesses. While this offers the advantage of spanning a large number of cases, it offers no guarantees that the optimal set of parameters represents a global minimum. It is thus entirely likely that a better set of parameter values would reduce the objective function value and improve the model skill, though again we note that the most important features of Nd biogeochemical cycling are likely to converge to similar values regardless of initial conditions.

Our specific choice of objective function gives a measure of the skill of the model for reproducing [Nd] and ε_Nd observations. Roughly speaking, despite our arbitrary choices for the weights ω_Nd, ${\mathit{\omega }}_{{\mathit{\epsilon }}_{\mathrm{Nd}}}$, and ω_p involved in Eq. (17), the value of the objective function can also be interpreted as the negative log-likelihood of observing the measured ε_Nd and [Nd], given the model and its parameters. In the future, we plan to determine the values of hyperparameters, such as ω_Nd, ${\mathit{\omega }}_{{\mathit{\epsilon }}_{\mathrm{Nd}}}$, and ω_p through a more formal Bayesian approach.

Qualitatively, GNOM compares well to previous models that are embedded in ocean general circulation models and simulate two explicit Nd isotopes (Arsouze et al., 2009; Rempfer et al., 2011; Gu et al., 2019; Pöppelmeier et al., 2020). Given that these models were built with a different objective, less available data, and without a systematic optimization of all parameters, we expect GNOM v1.0 to perform better against the objective function used in this study. However, we emphasize that these previous models are more suited to specific experiments, including simulations of transient changes in circulation on millennial timescales. In other words, these previous models are not all restricted by the caveats of a coarse steady-state circulation. We also note that the underlying circulation used in GNOM can be swapped with any other circulation available through the AIBECS.jl framework (see Pasquier et al., 2022b), although we reiterate that the OCIM v2.0 likely offers the best available representation of the current ocean circulation.

The main advantages of the GNOM are skill and computational efficiency. The GNOM v1.0 owes its low computational cost and quick simulation time to the steady-state OCIM v2.0 circulation in which it is embedded and the linear representation that allows us to solve the system of tracer equations in a single matrix inversion. The model's skill comes from the optimization procedure and likely benefits from the quality of the OCIM v2.0 circulation. The GNOM is also versatile in many respects owing to its simplicity. Parameter values and acceptable ranges can be tuned, entire mechanisms can be turned off, eliminating free parameters, or added with a few changes of simple lines of code. This versatility is compounded by computational speed, which makes GNOM ideally suited for quick experimentation and further optimization. The GNOM model is also easily diagnosed, owing to the powerful tools of linear algebra. Novel diagnostics offer new insights by revealing features often hidden in standard model output. Finally, as it is available in a self-contained package (except for the GEOTRACES dataset, which is not programmatically accessible), the GNOM v1.0 offers unprecedented reproducibility, which is sorely lacking in advanced research (Peng, 2011; Irving, 2016). Thanks to the advantages listed above, the GNOM is well positioned for answering long-lasting questions and exploring new ideas about the Nd cycle and the ε_Nd distributions.

Appendix A: Particulate Si model

To represent scavenging by opal (particulate silica), we designed and optimize a simple Si cycling model in parallel to our Nd cycling model. Our Si cycling model is a simple nutrient restoring model embedded in the same OCIM v2.0 circulation. We emphasize that the goal here is only to generate a reasonable 3D field for particulate biogenic silica concentrations.

The Si cycling model considered here explicitly tracks two tracers, DSi and PSi. We thus denote the modeled column vectors for DSi and PSi concentrations by the following: ${\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}$, and ${\mathit{\chi }}_{\mathrm{PSi}}^{\mathrm{mod}}$. Biological uptake of silicate in the euphotic zone is essentially modeled after the simple nutrient restoring scheme of the OCIMP-2 protocol (Najjar et al., 2007) with some slight modification. Specifically, the uptake rate ${\mathit{J}}_{\mathrm{up}}\left({\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}\right)$ vector is defined by the following:

$$\begin{array}{}\text{(A1)}& {\mathit{J}}_{\mathrm{up}}\left({\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}\right)=(\mathit{z}<z_{\mathrm{0}}){\displaystyle \frac{{\left({\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}-{\mathit{\alpha }}_{\mathrm{up}}{\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{obs}}\right)}^{+}}{{\mathit{\tau }}_{\mathrm{up}}}}\end{array}$$

where (z<z₀) ensures that uptake only happens in the euphotic zone, and ${\left({\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}-{\mathit{\alpha }}_{\mathrm{up}}{\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{obs}}\right)}^{+}/{\mathit{\tau }}_{\mathrm{up}}$ is the modified nutrient-restoring scheme. Technically, (z<z₀) is an abuse of notation and is here meant to be equal to 1 in the euphotic zone and 0 otherwise (z is the column vector depths and z₀=80 m is the depth of the base of the euphotic zone). The column vector ${\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{obs}}$ represents World Ocean Atlas silicate observations (WOA18; Garcia et al., 2019) regridded to the OCIM v2.0 grid. Here, (x)⁺ is a shortcut notation for x(x>0), such that $({\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}-{\mathit{\alpha }}_{\mathrm{up}}{\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{obs}}{)}^{+}/{\mathit{\tau }}_{\mathrm{up}}$ only activates when ${\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}>{\mathit{\alpha }}_{\mathrm{up}}{\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{obs}}$. The difference with the standard Ocean Carbon-cycle Model Intercomparison Project (OCMIP-2) protocol lies in the addition of the α_up modifier, which is a scalar that scales the field of observed silicate, allowing the optimized model to better fit observations. With α_up close to 1, this parameterization essentially allows for the model to take up silicate in the euphotic zone when the simulated DSi concentration exceeds the observations.

All the silicate that is taken up in this model is converted to sinking particulate PSi, which gravitationally settles with optimizable velocity parameter w_Si. Particulate biogenic silica is assumed to remineralize and redissolve into DSi in the water column with an optimizable timescale τ_rem. The rate of remineralization, ${\mathit{J}}_{\mathrm{rem}}\left({\mathit{\chi }}_{\mathrm{PSi}}^{\mathrm{mod}}\right)$, is thus simply defined by the following:

$$\begin{array}{}\text{(A2)}& {\mathit{J}}_{\mathrm{rem}}\left({\mathit{\chi }}_{\mathrm{PSi}}^{\mathrm{mod}}\right)={\mathit{\chi }}_{\mathrm{PSi}}^{\mathrm{mod}}/{\mathit{\tau }}_{\mathrm{rem}},\end{array}$$

which essentially closes the Si cycle.

Figure A1 Quantiles of the cumulative joint probability density functions of modeled and observed dissolved Si concentrations for the parallel Si cycling model. Darker colors indicate high density of data, such that n % of the modeled and observed data lie outside of the nth percentile contour. The closer the darker contours are to the 1:1 black dashed line, the better the fit.

Table A1 Optimized Si cycling model parameters.

Hence, the steady-state tracer equation for dissolved silicate (DSi) is as follows:

$$\begin{array}{}\text{(A3)}& {\mathbf{T}}_{\mathrm{circ}}{\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}={\mathit{J}}_{\mathrm{up}}\left({\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}\right)-{\mathit{J}}_{\mathrm{rem}}\left({\mathit{\chi }}_{\mathrm{PSi}}^{\mathrm{mod}}\right)+{\mathit{J}}_{\mathrm{geo}}\left({\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}\right),\end{array}$$

where ${\mathit{J}}_{\mathrm{geo}}\left({\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}\right)$ is an added term used to constrain the global inventory of Si in the system, which is set on geological timescales. In practice, we use the following:

$$\begin{array}{}\text{(A4)}& {\mathit{J}}_{\mathrm{geo}}\left({\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}\right)=\left(\right[\mathrm{DSi}{]}_{\mathrm{geo}}-{\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}})/{\mathit{\tau }}_{\mathrm{geo}},\end{array}$$

where [DSi]_geo is the optimizable global mean DSi concentration to which DSi concentrations are restored to with the timescale τ_geo=1 Myr. Conversely, the connected steady-state tracer equation for particulate biogenic silica (PSi) is as follows:

$$\begin{array}{}\text{(A5)}& \left({\mathbf{T}}_{\mathrm{circ}}+{\mathbf{T}}_{\mathrm{grav}}\right){\mathit{\chi }}_{\mathrm{PSi}}^{\mathrm{mod}}={\mathit{J}}_{\mathrm{rem}}\left({\mathit{\chi }}_{\mathrm{PSi}}^{\mathrm{mod}}\right)-{\mathit{J}}_{\mathrm{up}}\left({\mathit{\chi }}_{\mathrm{DSi}}^{\mathrm{mod}}\right),\end{array}$$

where T_grav is the vertical (downward) transport operator representing the flux divergence of particles. Particulate Si reaching the seafloor is not buried and instead remains in the deepest grid cell until eventual remineralization and redissolution into DSi.

A similar optimization procedure as for the Nd cycle is applied to optimize the five parameters collected in Table A1. The joint probability distribution of the mismatch between modeled and observed silicate is shown in Fig. A1. While not perfect overall, this simple model achieves a good fit to observations with a root mean square error of about 11.7 Mmol m⁻³.

Appendix B: Parameters

Figure B1 shows the prior distributions and optimized values of the parameters listed in Table 2.

Figure B1 Parameter prior distributions (color-filled densities) and initial and optimized parameter values for a dozen of optimization runs (lines). Each line starts by showing the initial parameter value at the top and is connected to the final optimized value at the bottom. The thicker blue line represents the optimization run that was used as our best estimate. (a) Dust ε_Nd values (ε_r) and solubilities (β_r) for each region of origin r. (b) Volcanic ash ε_Nd value (ε_volc) and solubility (β_volc). (c) Enhanced Nd release parameters, α curve parameters (curvature α_a and center α_c), Greenland enhancement (α_GRL), and global standard deviation of per-grid-cell in situ sedimentary ε_Nd (σ_ε). (d) Rivers and groundwater Nd concentrations (c_river and c_gw). (e) Global hydrothermal source magnitude (σ_hydro) and ε_Nd (ε_hydro). (f) Sedimentary flux parameters (ϕ₀, ϕ_∞, and z₀). (g) Reversible scavenging reaction equilibrium constants, K_X, and fractions returned to water when reaching seafloor, f_X, for each particle type X. See also Table 2.

Appendix C: Shift in effectively released ε_Nd

This appendix describes the effect that enhanced Nd release with extreme ε_Nd values has on effectively released ε_Nd. Let X be a random variable with the normal distribution 𝒩(μ_ε,σ_ε) denoting the observable in situ ε_Nd value, and let Y denote the random variable of the effectively released ε_Nd value. The mean expected value for Y is given by the following:

$$\begin{array}{}\text{(C1)}& \text{E}\left[Y\right]={\displaystyle \frac{\text{E}\left[X\mathit{\alpha }\right(X\left)\right]}{\text{E}\left[\mathit{\alpha }\right(X\left)\right]}},\end{array}$$

where α is defined in Eq. (8).

From the moments of a normal distribution, one can show, in the following, that:

$$\begin{array}{}\text{(C2)}& \text{E}\left[\mathit{\alpha }\left(X\right)\right]={\displaystyle \frac{a{\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}+a(c-{\mathit{\mu }}_{\mathit{\epsilon }}{)}^{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{10}}^{\mathrm{2}}}{{\mathit{\epsilon }}_{\mathrm{10}}^{\mathrm{2}}}},\end{array}$$

and, in the following, that:

$$\begin{array}{}\text{(C3)}& \mathrm{\text{E}}\left[X\mathit{\alpha }\left(X\right)\right]={\displaystyle \frac{a({\mathit{\mu }}_{\mathit{\epsilon }}-\mathrm{2}(c-{\mathit{\mu }}_{\mathit{\epsilon }}\left)\right){\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}+\left(a\right(c-{\mathit{\mu }}_{\mathit{\epsilon }}{)}^{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{10}}^{\mathrm{2}}){\mathit{\mu }}_{\mathit{\epsilon }}}{{\mathit{\epsilon }}_{\mathrm{10}}^{\mathrm{2}}}},\end{array}$$

so that, in the following:

$$\begin{array}{}\text{(C4)}& \mathrm{\text{E}}\left[Y\right]={\displaystyle \frac{a({\mathit{\mu }}_{\mathit{\epsilon }}-\mathrm{2}(c-{\mathit{\mu }}_{\mathit{\epsilon }}\left)\right){\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}+\left(a\right(c-{\mathit{\mu }}_{\mathit{\epsilon }}{)}^{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{10}}^{\mathrm{2}}){\mathit{\mu }}_{\mathit{\epsilon }}}{a{\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}+a(c-{\mathit{\mu }}_{\mathit{\epsilon }}{)}^{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{10}}^{\mathrm{2}}}}.\end{array}$$