A global scavenging and circulation ocean model of thorium-230 and protactinium-231 with improved particle dynamics (NEMO–ProThorP 0.1)

In this paper we set forth a 3-D ocean model of the radioactive trace isotopes ²³⁰Th and ²³¹Pa. The interest arises from the fact that these isotopes are extensively used for investigating particle transport in the ocean and reconstructing past ocean circulation. The tracers are reversibly scavenged by biogenic and lithogenic particles.

Our simulations of ²³⁰Th and ²³¹Pa are based on the NEMO–PISCES ocean biogeochemistry general circulation model, which includes biogenic particles, namely small and big particulate organic carbon, calcium carbonate and biogenic silica. Small and big lithogenic particles from dust deposition are included in our model as well. Their distributions generally compare well with the small and big lithogenic particle concentrations from recent observations from the GEOTRACES programme, except for boundary nepheloid layers for which, as of today, there are no non-trivial prognostic models available on a global scale. Our simulations reproduce ²³⁰Th and ²³¹Pa dissolved concentrations: they compare well with recent GEOTRACES observations in many parts of the ocean. Particulate ²³⁰Th and ²³¹Pa concentrations are significantly improved compared to previous studies, but they are still too low because of missing particles from nepheloid layers. Our simulation reproduces the main characteristics of the ²³¹Pa∕²³⁰Th ratio observed in the sediments and supports a moderate affinity of ²³¹Pa to biogenic silica as suggested by recent observations relative to ²³⁰Th.

Future model development may further improve understanding, especially when this will include a more complete representation of all particles, including different size classes, manganese hydroxides and nepheloid layers. This can be done based on our model as its source code is readily available.

1 Introduction

Oceanic circulation and the carbon cycle play a major role in the regulation of the past and present climate. Heat and carbon dioxide in the atmosphere tend to equilibrate with the ocean surface and are transported down into the deep ocean through the Meridional Overturning Circulation (MOC). Biogeochemical cycling also generates organic carbon that transfers into the deep ocean through particle sinking. Because of this, the strength of the MOC and particle removal participate actively in the regulation of the climate on the Earth.

Trace elements are also affected by these mechanisms and represent useful tools to provide constraints on these processes. The GEOTRACES programme has generated a large unique dataset that can now be used to better understand biogeochemical oceanic processes. Modelling quantifies and provides more information on the processes that control the oceanic distribution of these new observations. In present day climate, it is difficult to measure the MOC strength, and for past climate there are no measurements available at all. Isotopes and trace elements from sediment cores are used as proxies to infer past ocean circulation. Several examples include carbon isotopes (Broecker et al., 1988), the cadmium∕calcium ratio (Rosenthal et al., 1997), the ratio between protactinium-231 and thorium-230 (²³¹Pa∕²³⁰Th) (Böhm et al., 2015), and the neodymium isotope ratio ¹⁴³Nd∕¹⁴⁴Nd (Piotrowski et al., 2005). These proxies are affected by dynamical and biogeochemical processes. Including these proxies in a climate model is a way to better understand the climatic signal they register.

We will focus on ²³⁰Th and ²³¹Pa because these isotopes are well documented by the international GEOTRACES programme, and they are particularly suitable to study the transfer of particulate matter since the isotopes' source in the ocean is perfectly known: radioactive decay of uranium isotopes. Others have modelled ¹⁴³Nd and ¹⁴⁴Nd (Arsouze et al., 2008, 2009; Ayache et al., 2016), ¹³C (Tagliabue et al., 2009), and ¹⁴C (Mouchet, 2013). The ratio between ²³⁰Th and ²³¹Pa is used as a proxy for past ocean conditions, but this signal is potentially affected by both circulation and biogeochemical changes. Therefore, a correct understanding of the scavenging and underlying particle dynamics is essential in order to better simulate these tracers (Dutay et al., 2015).

Protactinium-231 and thorium-230 are produced in the ocean by the α decay of uranium-235 and uranium-234, respectively. Because the activity of uranium is approximately uniform in the ocean, ²³¹Pa and ²³⁰Th are produced at a relatively constant rate (${\mathit{\beta }}_{\mathrm{Pa}}=\mathrm{2.33}\times {\mathrm{10}}^{-\mathrm{3}}$ $\mathrm{dpm}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ and ${\mathit{\beta }}_{\mathrm{Th}}=\mathrm{2.52}\times {\mathrm{10}}^{-\mathrm{2}}$ $\mathrm{dpm}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$; dpm: disintegrations per minute) (Henderson and Anderson, 2003). They are both scavenged rapidly by the many particles that reside in the ocean and settle towards the sea floor. ²³¹Pa is less sensitive to particle scavenging than ²³⁰Th, which is reflected in the longer residence time of ²³¹Pa (80–200 yr) compared to that of ²³⁰Th (20–40 yr) (Yu et al., 1996). ²³¹Pa and ²³⁰Th are radioactive, decaying to radium isotopes and having a half-life of 32.76 and 75.40 kyr, respectively. Each combination of particle–radionuclide adsorption has a different reactivity. Other factors may affect the adsorption of radionuclides onto particles. For example, smaller particles have larger surface area to volume ratios, which results in an increase in the number of these radionuclides adsorbed per particle. Additionally, while the adsorption rate is expected to increase with particle concentration, the adsorption rate may be partly limited by the coagulation of non-filtered particles (or colloids) to filtered particles (Honeyman and Santschi, 1988; Honeyman et al., 1988). The vertical distributions of natural radionuclides, such as ²³⁰Th and ²³¹Pa, are hence sensitive to the distribution and mixture of particles. As a consequence of the different particle reactivities of ²³⁰Th and ²³¹Pa, the dissolved concentration ratio ²³¹Pa∕²³⁰Th deviates from the production activity ratio of 0.093 (Anderson, 2003; Anderson et al., 1983; Rutgers van der Loeff et al., 2016).

As both particle dynamics and circulation of the ocean affect ²³⁰Th and ²³¹Pa, numerical biogeochemical general circulation models are used to study the relative contribution of these mechanisms. The isotopes have been simulated in models of intermediate complexity, for instance by Henderson et al. (1999) (LSG-OGCM), Marchal et al. (2000) (EMIC 2.5D), Heinze et al. (2006, 2018) (HAMOCC) and Luo et al. (2010). More complex ocean general circulation models have also been used to simulate these tracers, namely by Dutay et al. (2009) (NEMO–PISCES), Siddall et al. (2007), Rempfer et al. (2017) (Bern3D), and Gu and Liu (2017) (CESM). Dutay et al. (2009) demonstrated that the particle concentration simulated by the PISCES model in the deep ocean was too low by a factor of 2 in the mesopelagic zone to 50 in the deepest ocean. This led to overestimated radionuclide concentrations in the deep ocean. Therefore, it is crucial to improve the representation of the particles (Dutay et al., 2009). Rempfer et al. (2017) showed that taking into account additional sinks at the sea floor and at the ocean margins yields an improved agreement with observations, especially for the dissolved phases of ²³⁰Th and ²³¹Pa. Particulate ratios improved to a lesser extent; the authors have not presented an evaluation of their simulated particulate concentrations. Gu and Liu (2017) had similar goals as this study.

In this study, we try to improve on previous studies that simulate the distribution of ²³⁰Th and ²³¹Pa. Our approach is to improve the mechanistic description of the particle and radionuclide cycling, as well as on the side of system design and reusability of the model. We evaluate how well the model fits with observations, but tuning is not one of our main goals. Just like Dutay et al. (2009), we use the NEMO–OPA (Nucleus for European Modelling of the Ocean – Océan PAralelisé) ocean general circulation model (Madec and the NEMO team, 2016) and the PISCES biogeochemical model (Aumont et al., 2015). Aumont et al. (2017) showed that dissolution rates of particulate organic carbon (POC) in PISCES were overestimated, and they improved on this by introducing a spectrum of different labilities. This improved the simulation of both small and big POC significantly, so we use this same model for our simulations. There are several improvements since Dutay et al. (2009); namely, the new model

• includes lithogenic particles from dust deposition;

• has improved the biogeochemistry, affecting the biogenic particle distributions (Aumont et al., 2015, 2017);

• includes three different phases per nuclide (dissolved, and adsorbed onto big (and small) particles), whereas Dutay et al. (2009) include only a single compartment (total concentration) for each nuclide from which they calculated the respective phases based on chemical equilibrium;

• is more precise and explicit on the mathematical formalism;

• is written in Fortran 95 instead of Fortran 77, making the extension of its use to other modern Fortran models easier; and

• is part of a modern model framework, NEMO–TOP, facilitating its use with other models in this framework.

The main objective of this paper is to improve on the simulation of ²³⁰Th and ²³¹Pa (the dissolved phase and two particulate size classes of particles for both nuclides) based on an improved modelling of small and big particles.

New observations are available from the GEOTRACES programme. Especially the North Atlantic GA03 transects will be used for model validation because on this transect not only dissolved ²³⁰Th and ²³¹Pa have been measured (Hayes et al., 2015a), but also their adsorbed forms (Hayes et al., 2015b) and biogenic and lithogenic particle concentrations in two size classes (Lam et al., 2015).

2 Model description

In order to simulate the biogenic particle dynamics and their interaction with the ²³⁰Th and ²³¹Pa trace isotopes, we use the biogeochemical circulation model NEMO–PISCES (Aumont et al., 2015; Madec, 2008). This model has been employed for many other studies concerning trace metals, as well as large-scale ocean biogeochemistry (Arsouze et al., 2009; Dutay et al., 2009; Gehlen et al., 2007; Tagliabue et al., 2010; van Hulten et al., 2013, 2014, 2017b). We force PISCES by a climatological year of circulation fields (including turbulent diffusion) that was obtained from the dynamical component of NEMO. Table 1 gives an overview of this and other components of the model and their relevant properties.

(Aumont et al., 2015)(Aumont et al., 2017)

Table 1 Model components and properties essential for the radionuclide model. OPA stands for Océan PAralalisé.

All model fields are defined on the ORCA2 discrete coordinate system, an irregular grid covering the whole world ocean with a nominal resolution of ${\mathrm{2}}^{\circ }\times {\mathrm{2}}^{\circ }$, with an increased resolution in the meridional direction near the Equator and Antarctica and in both horizontal directions in the Mediterranean, Red, Black and Caspian seas. On the Northern Hemisphere, it has two coordinate singularities, one in Canada and the other in Russia, such that both singularities fall outside the computational domain. The vertical resolution of the ORCA2 grid is 10 m in the upper 100 m, increasing downwards to 500 m such that there are 30 layers in total and the ocean has a maximum depth of 5000 m (Madec and Imbard, 1996; Murray, 1996). The time step of the model is 6 h for the dynamics and the radionuclides and 1.5 h for the biogeochemistry (PISCES) and the lithogenic particles. When necessary, sub-time stepping is done for all sinking components in the model.

2.1 Circulation

The circulation was obtained by forcing NEMO in the ORCA2 configuration with climatological air–sea boundary conditions consisting of heat, freshwater and momentum fluxes that were derived from bulk formulae. They are functions of wind, sea surface temperature, air temperature, air humidity and evaporation minus precipitation. For the tropics, daily wind stress was used, which was based on European Remote Sensing (ERS) satellite data, and for the polar regions NCEP/NCAR reanalysis data were used (Kalnay et al., 1996; Kistler et al., 2001). Surface salinity was restored with a timescale of 60 days towards the seasonal Polar Science Center Hydrographic Climatology (PHC) dataset to avoid model drift (Timmermann et al., 2005). The last year of this 200-year simulation is used as our 1-year climatology with a resolution of 5 days of the dynamics.

2.2 Particle dynamics

This version of PISCES includes two size classes of POC, both with differential remineralisation rates (Aumont et al., 2017), one size class of biogenic silica (bSiO₂) and one class of calcium carbonate (CaCO₃) (Fig. 1). In the model, particles sink down with two velocities: w_s=2 m d⁻¹ and w_b=50 m d⁻¹. Aumont et al. (2015) mark the diameter boundary at 100 µm, but as we need to compare with observations, the size classes will be taken to correspond with “small” particles (∅≤51 µm) and “big” particles (∅>51 µm).

Figure 1 Conceptual model of the particle dynamics on which the numerical model is based. Nanophytoplankton and diatoms (in green) take up nutrients and CO₂, which are released again from respiration and remineralisation of bSiO₂, POC, dissolved organic matter and lithogenic particles. The nutrients, CO₂ and dissolved organic matter are not represented in the figure because only particles impact ²³¹Pa and ²³⁰Th. Zooplankton are denoted by the black box. All sinking particles are denoted by blue boxes. Effectively, this figure comprises the internal cycling of PISCES, minus details that are not of interest here, plus the lithogenic dust model. Small and big particulate organic matter is denoted by sPOM and bPOM, both subject to the differential lability scheme (Aumont et al., 2017), and bSiO₂ stands for biogenic silica. Similarly, sLith and bLith stand for small and big lithogenic particles. Sinking is denoted by the red arrows; triple arrows mean fast, and normal arrows slow.

The dissolution equation for calcium carbonate is given by

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\text{d}\left[{\mathrm{CaCO}}_{\mathrm{3}}\right]}{\text{d}t}}=-R\left[{\mathrm{CaCO}}_{\mathrm{3}}\right],\end{array}$$

where the rate variable for calcite dissolution is

$$\begin{array}{}\text{(2)}& R=k\cdot (\mathrm{1}-\mathrm{\Omega }{)}^n,\end{array}$$

where k and n are two somewhat restricted, though tunable, parameters that signify a dissolution rate constant and a reaction rate order, respectively. The calcite saturation state is given by

$$\begin{array}{}\text{(3)}& \mathrm{\Omega }={\displaystyle \frac{\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{K_{\mathrm{sp}}^{\prime }}},\end{array}$$

where $K_{\mathrm{sp}}^{\prime }$ is the product of the saturation concentrations of calcium and carbonate. This can be approximated by

$$\begin{array}{}\text{(4)}& \mathrm{\Omega }\approx {\displaystyle \frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}{]}_{\text{sat}}}}\end{array}$$

because relative variations in [Ca²⁺] are small. For CaCO₃ we changed the standard first-order dissolution kinetics parameterisation to a fourth-order dissolution based on evidence of Keir (1980) and subsequent studies. For our simulation we used a calcite dissolution rate constant of k=2.5 mo⁻¹ and a dissolution order of n=3.9 (Subhas et al., 2015).

In addition to biogenic particles, we introduced lithogenic dust particles in the model. The yearly average dust flux is derived from Hauglustaine et al. (2004) and is presented in Fig. 2. It is used by the model as the input of lithogenic ocean particles and for nutrient supply in PISCES. This dust deposition field has been tested in biogeochemical studies with this configuration of NEMO (van Hulten et al., 2013).

Figure 2 Dust deposition on a logarithmic scale ($\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). It is the integrated flux over the 12-month climatology based on Hauglustaine et al. (2004).

Small and big lithogenic particles are added to the upper layer of the ocean according to

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial P^{\mathrm{Litho}}}{\partial t}}{|}_{\mathrm{surface}}={\displaystyle \frac{f}{\mathrm{\Delta }{z}_{\mathrm{1}}}}\cdot {\mathrm{\Phi }}_{\mathrm{dust}},\end{array}$$

where P^Litho is the small (big) lithogenic particle concentration, f is the fraction of the dust that gets partitioned into the small (big) lithogenic particles in the ocean, Δz₁=10 m is the thickness of the upper model layer and Φ_dust is the dust flux.

Our model has two size classes for lithogenic particles, so this equation is applied for two different concentrations P^Litho and respective fractions f. We set the small lithogenic dust flux fraction to 20 % and the big one to 80 %. Once partitioned in the ocean, the lithogenic particles sink down, changing their concentrations throughout the ocean according to

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}{P}^{\mathrm{Litho}}}{\mathrm{d}t}}=-w\cdot {\displaystyle \frac{\partial P^{\mathrm{Litho}}}{\partial z}}+(\mathcal{A}{\mathrm{\nabla }}_{\mathrm{h}}^{\mathrm{2}}+\mathcal{B}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}})P^{\mathrm{Litho}},\end{array}$$

where w is the settling velocity set to the constant 2 m d⁻¹ for the small lithogenic particles and to 50 m d⁻¹ for the big particles. The depth, z, is positive upwards, ∇_h is the horizontal divergence, and 𝒜 and ℬ are respectively the horizontal and vertical eddy diffusivity coefficients. The material derivative includes a term for eddy-induced velocity (Gent and McWilliams, 1990; Gent et al., 1995). Equation (6) is identical to the settling of small and big POC (Aumont et al., 2015). Of course, there are also biological and chemical sources and sinks for POC. However, for lithogenic particles there are no such sources or sinks because the only source in our model is dust deposition and we assume the lithogenic particles are refractory. The lithogenic particles are removed from the model domain when arriving at the sea floor, which means that they are buried in the sediment.

2.3 Radionuclides

Thorium-230 and protactinium-231 are produced throughout the ocean from the decay of ²³⁴U and ²³⁵U, respectively. Because of long residence times of over 200 kyr, these uranium isotopes are approximately homogeneously distributed throughout the ocean and do not change much over time (Ku et al., 1977). The residence time of ²³⁴U is 3.2–5.6×10⁵ a, which is much longer than the full mixing time of the world ocean of about 10³ a (Dunk et al., 2002). The ²³⁴U concentration can vary about 10 %, depending mostly on the salinity (Owens et al., 2011), but that is smaller than uncertainties arising from other assumptions in our model. Therefore, we assume that the production rates of ²³⁰Th and ²³¹Pa are constant both in space and time.

Figure 3 The conceptual reversible scavenging model for the radionuclides. The radioisotopes, ²³⁰Th and ²³¹Pa, are depicted in orange when in the dissolved phase. Just like with the other PISCES tracers, the dissolved and particulate radionuclides are transported by circulation and eddy diffusion.

The ²³⁰Th and ²³¹Pa radionuclides are reversibly scavenged by biogenic and lithogenic particles (Fig. 3). We will assume, as in previous studies (Dutay et al., 2009; Siddall et al., 2005), that the adsorption and desorption reaction rates are sufficiently fast compared to radionuclide production, decay, advection, mixing, change in particle distribution and settling of the adsorbed phases (see discussion in Sect. 6). Because of that, we equilibrate between the dissolved and adsorbed phases instantly at each time step. This means, only considering adsorption and desorption processes at this point, that we must solve this set of equations for every nuclide i:

$$\begin{array}{}\text{(7a)}& {\displaystyle }\mathrm{0}& {\displaystyle }={\displaystyle \frac{\partial }{\partial t}}(A_{i,D}+A_{i,S}+A_{i,B}),\text{(7b)}& {\displaystyle }{A}_{i,S}& {\displaystyle }=A_{i,D}\sum _{j\in S}{K}_{ij}{P}^j,\text{(7c)}& {\displaystyle }{A}_{i,B}& {\displaystyle }=A_{i,D}\sum _{j\in B}{K}_{ij}{P}^j,\end{array}$$

where A_i stands for the activity of nuclide $i\in \mathit{\{}{}^{\mathrm{230}}\mathrm{Th},{}^{\mathrm{231}}\mathrm{Pa}\mathit{\}}$¹ and where P^j stands for the concentration of particle $j\in \mathit{\{}\mathrm{sPOC},\mathrm{bPOC},{\mathrm{CaCO}}_{\mathrm{3}},{\mathrm{bSiO}}_{\mathrm{2}},\mathrm{sLitho},\mathrm{bLitho}\mathit{\}}$ (in gram per gram of seawater, so strictly this is a mass fraction). D is the set of non-sinking phases (here only dissolved), S of small particles that sink with 2 m d⁻¹ and B of big particles that sink with 50 m d⁻¹ (in the same way as Eq. 6 for lithogenic particles). For any particle size class $J\in \mathit{\{}S,B\mathit{\}}$ and any radionuclide i, we define $A_{i,J}={\sum }_{j\in J}{A}_{ij}$. Finally, K_ij is the equilibrium partition coefficient of nuclide i for particle j. As with small and big POC, adsorbed ²³⁰Th and ²³¹Pa should aggregate and disaggregate. However, any transference of radionuclides between two particle types or size classes, in addition to the instant equilibration, would not have any effect as there would quickly be a new equilibrium of the radionuclides between the particles and the dissolved phase at the next instant re-equilibration (which happens each time step). Therefore, there is no equation for the radionuclides addressing aggregation.

Since we cannot solve this analytically, this will be done numerically by first assigning a new value to the dissolved nuclide activity:²

$$\begin{array}{}\text{(8)}& A_{i,D}:={\displaystyle \frac{A_{i,D}+A_{i,S}+A_{i,B}}{\mathrm{1}+{\sum }_j{K}_{ij}{P}^j}}.\end{array}$$

Then, we calculate the activity of i that is adsorbed onto small and big particles by applying Eqs. (7b) and (7c). With this approach, the small and big adsorbed concentrations equilibrate instantly. Assuming that the change in adsorption strength is much smaller than the relative change in tracer activity, ${\forall }_{i,J:}{\partial }_t{\sum }_{j\in J}{K}_{ij}{P}^j\ll {\partial }_t{A}_{i,D}/A_{i,D}$, the total activity of every i is conserved, i.e. Eq. (7a) holds.

Proof. Let the total adsorption strength for any isotope i be $Q_i={\sum }_j{K}_{ij}{P}^j$ and the total amount of the same isotope $T_i=A_{i,D}+A_{i,S}+A_{i,B}$, and let primes denote the updated concentrations. Assume that the adsorption strength for every isotope i is constant ($Q_i^{\prime }\equiv Q_i$).

$$\begin{array}{ll}{\displaystyle }{T}_i^{\prime }& {\displaystyle }={\displaystyle \frac{T_i}{\mathrm{1}+Q_i^{\prime }}}+A_{i,D}^{\prime }{Q}_i^{\prime }\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{T_i}{\mathrm{1}+Q_i}}+{\displaystyle \frac{T_i}{\mathrm{1}+Q_i}}{Q}_i\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{T_i}{\mathrm{1}+Q_i}}(\mathrm{1}+Q_i)=T_i\\ {\displaystyle }\Rightarrow {\partial }_t{T}_i^{\prime }& {\displaystyle }={\partial }_t{T}_i=\mathrm{0}.\mathit{\square }& {\displaystyle }& {\displaystyle }\end{array}$$

The ²³⁰Th and ²³¹Pa that are adsorbed onto the particles follow the same law as the small and big (lithogenic) particles (Eq. 6). Of course, by definition, the adsorbed radioisotope and the particle settle with the same speed, and thus we have implemented it.

The decay terms of ²³⁰Th and ²³¹Pa are much smaller than the other sources and sinks, but they are included in the model:

$$\begin{array}{}\text{(9a)}& {\displaystyle \frac{\mathrm{d}{A}_{i,D}}{\mathrm{d}t}}& {\displaystyle }={\mathit{\beta }}_i-{\mathit{\lambda }}_i{A}_{i,D}+(\mathcal{A}{\mathrm{\nabla }}_h^{\mathrm{2}}+\mathcal{B}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}})A_{i,D},\text{(9b)}& {\displaystyle \frac{\mathrm{d}{A}_{i,S}}{\mathrm{d}t}}& {\displaystyle }=-{\mathit{\lambda }}_i{A}_{i,S}-w_{\mathrm{s}}{\displaystyle \frac{\partial A_{i,S}}{\partial z}}+(\mathcal{A}{\mathrm{\nabla }}_h^{\mathrm{2}}+\mathcal{B}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}})A_{i,S},\text{(9c)}& {\displaystyle \frac{\mathrm{d}{A}_{i,B}}{\mathrm{d}t}}& {\displaystyle }=-{\mathit{\lambda }}_i{A}_{i,B}-w_{\mathrm{b}}{\displaystyle \frac{\partial A_{i,B}}{\partial z}}+(\mathcal{A}{\mathrm{\nabla }}_h^{\mathrm{2}}+\mathcal{B}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}})A_{i,B},\end{array}$$

where β_i and λ_i are respectively the production rate and radioactive decay of isotope i.

3 Simulations

For the numerical simulations, we forced the just described model, an adjusted version of PISCES with lithogenic particles and radiotracers, with the off-line circulation fields. The model was spun up for 500 years, after which it was in an approximate steady state (decadal drift of −0.002 % for total ²³⁰Th and +0.058 % for total ²³¹Pa). Protactinium-231 has a larger drift than thorium-230 because ²³⁰Th is more quickly removed everywhere in the ocean because of its high particle reactivity. The lithogenic particles are in a steady state, and the PISCES variables are in an approximate steady state (e.g. phosphate shows a drift of −0.005 % per decade).

The partition coefficients K_ij, with $i\in \mathit{\{}\mathrm{Th},\mathrm{Pa}\mathit{\}}$, depend on the type of particle j and are given in Table 2.

Table 2 Partition coefficients for the different modelled particles in 10⁶ g of seawater per gram of particles. In parentheses is the value for the sensitivity simulation discussed in Sect. 6.1.

These coefficients still have large uncertainties, but their values can be constrained by reported values from different experimental studies (Chase et al., 2002; Geibert and Usbeck, 2004; Hayes et al., 2015b). Therefore, we had quite some freedom in prescribing values. The adsorption onto calcium carbonate is a factor of 2 decreased from Chase et al. (2002). This brings the K value of Pa closer to that in Hayes et al. (2015b) based on field measurements, namely 0.9±0.4 Mg g⁻¹ (10⁶ g of seawater per gram of particles), but since we maintained the ratio of K values, the partition coefficient of Th is much smaller than found by Hayes et al. (2015b). It is consistent with the laboratory results of Geibert and Usbeck (2004), but they cleaned the biogenic silica aggressively, bought carbonate that should be organic-free and clay (smectite) was purchased as a commercial standard and not cleaned (Geibert and Usbeck, 2004; and Walter Geibert, personal communication, January 2018). In the real ocean, solid minerals are not clean but have organic matter around them, giving rise to increased scavenging (Chuang et al., 2014). Therefore, our K values for at least biogenic silica and calcium carbonate may be lower than what is typically considered realistic. The actual reason for the relatively low K values is that we only have fast-sinking biogenic silica and calcium carbonate particles in the model that also represent smaller particles that sink slower. Thus using higher K values would result in too much export of Pa and Th. For small lithogenic particles, K is set about a factor of 5 larger than in the literature (Geibert and Usbeck, 2004; Hayes et al., 2015b), whereas big lithogenic particles have a smaller value than reported in the literature. The Supplement contains a more exhaustive table of information, including an estimated literature range (Table S3). The consequences of the chosen values for the partition coefficients are discussed further in Sect. 6.

4 Observations

In this study, we will focus on the GEOTRACES GA03 transect in the North Atlantic Ocean. Recently, a large number of measurements on both radionuclides and their carrier particles have been collected on this transect. This unique combination makes it especially useful to evaluate a radionuclide scavenging model. We will also compare our global ocean model with several observational datasets throughout the ocean (Table 3). Observations obtained from the GEOTRACES programme are denoted with the respective GEOTRACES transect number (Mawji et al., 2015; Schlitzer et al., 2018).

Brun-Cottan et al. (1991)Druffel et al. (1992)Lam et al. (2011)Lam et al. (2015)Bacon and Anderson (1982)Roy-Barman et al. (1996)Henderson et al. (1999)Marchal et al. (2007)Venchiarutti et al. (2011)Venchiarutti et al. (2011)Hayes et al. (2013)Deng et al. (2014)Hayes et al. (2015a, b)Burckel et al. (2015, 2016)

Table 3 Observations used for comparison with the model simulations.

^* http://climotope.earth.ox.ac.uk/data_compilations (last access: 27 August 2018).

For the carrier particles most of our data come from Lam et al. (2015), which is a recent GEOTRACES dataset at the GA03 transect in the North Atlantic Ocean. An older compilation of particles is taken from Lam et al. (2011). We use some older data as well (Table 3). Particle concentrations were determined by filtering seawater in situ.

Concentrations of both dissolved and particulate phases of ²³⁰Th and ²³¹Pa were taken from Hayes et al. (2015a) (GA03) and Hayes et al. (2014) (Pacific Ocean). Other data are listed in Table 3.

To evaluate the sediment ²³¹Pa∕²³⁰Th flux of the model, we will compare with compilations of the Holocene (i.e. top core particulate concentrations) (Table 3).

5 Results

5.1 Circulation

As mentioned before, our model is part of a global general circulation model, but instead of solving the Navier–Stokes equations, our tracers are advected by the circulation fields of a previous simulation of the dynamical ocean model. Since the overturning circulation is of importance for the redistribution of the tracers, we present basic results here.

Figure S1 of the Supplement presents the overturning stream function (OSF) of the Atlantic Ocean. The OSF is defined as the zonally (through the basin) and vertically (from the surface downwards) integrated meridional component of the current velocity. We use this as a measure for the Atlantic Meridional Overturning Circulation (AMOC). The upper overturning cell transports 14 Sv (1 Sv =10⁶ m³ s⁻¹) on average. Observations suggest higher values of about 19±5 Sv (Cunningham et al., 2007; Rayner et al., 2011; Sinha et al., 2018). The lower cell has an overturning strength of about 6 Sv, which is stronger than the about 2 Sv from estimates (Talley et al., 2003). This overestimation is mostly due to the fact that the AMOC is shallow. The AMOC reaches about 2500 m of depth, whereas observational studies show the deep water sinks down to about 4000 m. At around 30^∘ N the AMOC has weakened notably, which is consistent with some studies (Johnson, 2008), but the Antarctic Bottom Water (AABW) does not go as far north as other studies suggest (Middag et al., 2015).

The Antarctic Circumpolar Current is also an important feature for large-scale ocean circulation. The throughflow at the Drake Passage is a good measure for that. In our model the flux through the Drake Passage is 200 Sv, which slightly overestimates recently reported values (Donohue et al., 2016).

5.2 Particles

Figure 4 Particle concentrations from the PISCES model at the surface ocean. Observations are represented as coloured discs (Brun-Cottan et al., 1991; Druffel et al., 1992; Lam et al., 2011, 2015). The unit of every modelled and measured particle concentration is mg m⁻³.

Figure 4 shows the surface concentrations of total POC, big calcium carbonate, big biogenic silica and total lithogenic particles. Observations are plotted as coloured discs on top of the modelled concentrations. The modelled biogenic particles include living matter. For the model, we assume POC includes all phytoplankton and microzooplankton (not mesozooplankton since it may swim away), calcium carbonate contains the assumed fraction of CaCO₃ in the phytoplankton, and biogenic silica includes living diatoms. The modelled concentrations and patterns of total POC compare well with the observations (Fig. 4a). For instance, coastal and equatorial Pacific POC concentrations are elevated in both the model and the observations compared to other regions, although the spatial extent of the oligotrophic regions appears to be too small in the model. We compare our modelled CaCO₃ (big, fast-sinking calcium carbonate particles), with the big calcium carbonate particles from observations. As the model tries to represent CaCO₃ with only one size class, neither comparing with only big particles nor comparing with total CaCO₃ is completely fair. Incidentally, we have data of big CaCO₃ for both the Atlantic and Pacific Ocean, whereas only for the Atlantic Ocean do we have small particle data. Even though the meridional patterns of big CaCO₃ particles are reproduced, the concentrations are generally overestimated, especially in the Gulf of Alaska (Fig. 4b). Contrarily, total observed [CaCO₃], which is only provided along section GA03 where small particulate CaCO₃ measurements are available, is higher than the prediction of our modelled (only big) CaCO₃ particle concentration. The model produces a reasonable spatial distribution of biogenic silica: there are more elevated values in the high latitudes, but concentrations are underestimated at the surface at low latitudes (Fig. 4c).

Figure 5 Particle fields at the GA03 North Atlantic GEOTRACES transect; observations (Lam et al., 2015) as coloured discs.

The recent large dataset obtained from the full-ocean depth GEOTRACES GA03 transect (Lam et al., 2015) provides an opportunity to analyse the performance of PISCES in more detail. PISCES generally produces the right order of magnitude for all four biogenic particles (small and big POC, calcium carbonate, and biogenic silica), but there still remain some shortcomings in their distributions (Fig. 5a–d).

• The total POC concentration in Fig. 5a is up to a factor of 4 underestimated in the deep oligotrophic ocean. In the upper 200 m of the ocean the model overestimates the observations, though at some points at the surface the POC concentration is underestimated (also Fig. 4a). In both the observations and the model the fraction of big POC varies from zero to 0.6 (Fig. 5), but the spatial distributions are very different. More detailed results and discussion on the simulation of POC are to be found in Aumont et al. (2017).

• The model underestimates the CaCO₃ concentrations from 70 to 25^∘ W by a factor of 2 to 10, but east of that up to Africa the prediction lies close to or overestimates the observations. In the Canary Basin and up to Portugal (meridional transect at the right), the model reproduces the right order of magnitude, but it does not reproduce the correct profiles everywhere (Fig. 5c).

• Biogenic silica concentrations are generally reproduced but the model overestimates the higher concentrations observed along the western margin (Fig. 5d).

Finally, the modelled total lithogenic particle concentration at the surface shows, as expected, a close resemblance with dust deposition patterns (Fig. 2) and mostly compares well with observations (Fig. 4d). Only near the coasts of the USA and Portugal does the model strongly underestimate lithogenic particle concentrations. Concerning the deep ocean, our model captures quite nicely the general distribution at the GA03 transect (Fig. 5e and f). However, the concentrations near the western boundary are strongly underestimated, especially those of big lithogenic particles. Our model reproduces the observed fraction of big lithogenic particles of ∼0.1 to ∼0.3 in most of the deep ocean, but it underestimates the much larger fraction of big particles observed in the upper 200 m of the ocean (0.1–0.2 in the model versus 0.3–0.9 in the observations).

5.3 Thorium-230 and protactinium-231

Figure 6 The dissolved thorium-230 activity at four depth levels (mBq m⁻³); observations as discs on the same colour scale.

Figure 7 The dissolved protactinium-231 activity at four depth levels (mBq m⁻³); observations as coloured discs.

[²³⁰Th]_D stands for the concentration of ²³⁰Th that is in the dissolved phases. Similarly, [²³⁰Th]_S is the concentration of ²³⁰Th that is adsorbed onto small particles S, [²³⁰Th]_B the concentration adsorbed onto big particles B and similarly for ²³¹Pa. The modelled dissolved distributions of ²³⁰Th and ²³¹Pa are underestimated at the surface (Figs. 6 and 7). Below the surface, [²³⁰Th]_D is much better captured by the model (also Fig. 8a for the GA03 transect). However, observations show lower values near the bottom in the western part of the GA03 section and at 45^∘ W associated with more intense scavenging related to respective nepheloid layers and manganese oxides from hydrothermal vents, which are not produced in the simulation (Fig. 8a and b).

In the intermediate and deep waters of all the oceans, modelled [²³¹Pa]_D is generally of the correct order of magnitude but is strongly overestimated below 2500 m of depth (Fig. 8b), especially in the Pacific Ocean (Fig. 7c, d).

The concentrations of the adsorbed phases of ²³⁰Th at GA03 are presented in Fig. 8c and e. In the deep ocean, modelled [²³⁰Th]_S and [²³⁰Th]_B have lower values than [²³⁰Th]_D, which is consistent with observations. Furthermore, compared with Dutay et al. (2009), we have notably improved [²³⁰Th]_S and [²³⁰Th]_B. Other global modelling studies have not reported adsorbed concentrations of ²³⁰Th and ²³¹Pa. However, the model still underestimates the observed [²³⁰Th]_S and [²³⁰Th]_B.

Figure 8 Radionuclide concentrations at the GA03 North Atlantic GEOTRACES transect. Panels (a) and (b) display dissolved concentrations, panels (c) and (d) show the amounts on small particles, and panels (e) and (f) those on big particles. Everything is in mBq m⁻³. Observations (Hayes et al., 2015a) are shown as coloured discs.

Contrarily, small ²³¹Pa particles are overestimated in our model below 1500 m at the GA03 transect (Fig. 8d). Big ²³¹Pa particles are simulated more realistically (Fig. 8f).

Figure 9 Radionuclide concentrations in mBq m⁻³ at stations in the North Pacific (around the ALOHA station, 158^∘ W, 23^∘ N) and North Atlantic Ocean (around the BATS station, 64^∘ W, 29^∘ N). The left two panels (a, b) display [²³⁰Th]_D, and the right two (c, d) [²³¹Pa]_D. The line is the model, and the speckles are the observations.

To clearly show the change in ²³⁰Th and ²³¹Pa with respect to depth, profiles are plotted in Fig. 9. [²³⁰Th]_D is linearly increasing downwards, which matches the observations, but below 3000 m near Bermuda, [²³⁰Th]_D increases where it should decrease according to the observations. The shape of [²³¹Pa]_D matches that of the observations for the upper 2 km, below which our model overestimates the observations by up to a factor of 4 at the bottom of the ocean.

Table 4 Relative global budget of ²³¹Pa and ²³⁰Th in the different particulate phases (“global stock”) in percentages. Globally averaged fluxes in the different phases (“global particle flx”), and the same but for the North Atlantic Ocean (44–24^∘ W, 25–56^∘ N) and the Weddell Sea (44–24^∘ W, 76–63^∘ S), all in percentages associated with each particle type.

Globally in the model, only 0.3 % of the ²³¹Pa sits in the particulate pool, and for ²³⁰Th this is 1.3 %. The rest is in the dissolved phase. Observational data show that more than 10 % of ²³⁰Th sits in the particulate pool, typically adsorbed onto small calcium carbonate particles. The underestimation by our model can be explained by the fact that we only have big, fast-sinking CaCO₃ particles in the model. Since the big particles are removed quickly, these cannot act as a realistic CaCO₃-associated ²³⁰Th pool without depleting ²³⁰Th, and thus we underestimate the particulate radionuclide concentrations.

Table 4 gives two types of information on the adsorbed phases. One is the amount fraction of radionuclide adsorbed onto each particle (“global stock”). The other is the fraction of the flux carried by each type of particle (the other rows) for the global ocean and two regions of the ocean. No distinction is made between big and small particles (they are summed). The stock and the flux are not necessarily the same because the big particles have a different settling velocity from the small particles. For example, only for POC and lithogenic particles do we have slowly sinking particles on which much of the radionuclides are adsorbed, whereas the big affinity of Pa to bSiO₂ mostly results in a strong export (and similarly for Th on CaCO₃).

Only a small amount of ²³⁰Th is adsorbed onto biogenic silica (3 %), though the amount is larger for ²³¹Pa (11 %). However, the global settling flux of ²³¹Pa is primarily (54 %) determined by bSiO₂, particularly due to the high bSiO₂ fluxes in the Southern Ocean. Most of the modelled particulate ²³¹Pa and ²³⁰Th is on the lithogenic particles, but lithogenic particles account for only 14 % of the ²³¹Pa flux and 9 % of the ²³⁰Th flux. If lithogenic particles from nepheloid layers were included, this fraction would be higher. The flux contribution is also different from what is in the pools when we look at CaCO₃: only 24 % of the ²³⁰Th is on calcium carbonate, but CaCO₃ is responsible for 76 % of the ²³⁰Th export. These discrepancies between the contribution of each particle type to the flux and to the stock of radionuclides arise from the different settling velocities for each type of particle.

The relative contributions of the different particulate phases also vary for different regions of the ocean. Most notably, biogenic silica is responsible for 94 % of ²³¹Pa export in the Southern Ocean.

Compared to Dutay et al. (2009), the ²³⁰Th profile shape is improved, but only to a small extent due to the large improvement in the POC representation (Aumont et al., 2017) because the POC accounts for only 6 % of the vertical transport of ²³⁰Th (Table 4). The improvement is mostly due to the addition of lithogenic dust particles and the adjusted CaCO₃ dissolution. Lithogenic particles and CaCO₃ hardly dissolve compared to POC and hence the downward particle flux remains constant, which is in agreement with 1-D models which show linear profiles of ²³⁰Th_D (Roy-Barman, 2009).

5.4 Sedimentation flux ratio of ²³¹Pa∕²³⁰Th

Figure 10 Sedimented ²³¹Pa∕²³⁰Th ratio. Top core measurements are presented as coloured discs. This is the standard simulation with $K_{\mathrm{Pa},{\mathrm{bSiO}}_{\mathrm{2}}}/K_{\mathrm{Th},{\mathrm{bSiO}}_{\mathrm{2}}}=\mathrm{1}$.

The sedimentation flux of the total adsorbed ²³¹Pa∕²³⁰Th ratio captures the general patterns of the upper sediment core ²³¹Pa∕²³⁰Th concentration ratio (Fig. 10). Consistent with observations, the model produces more elevated ²³¹Pa∕²³⁰Th ratio values in the Southern Ocean, the northern part of the Indian Ocean, the Pacific Ocean and along the coastal upwelling regions. However, at some places, like much of the Southern Ocean, our model tends to overestimate the ratio derived from the sediment core observations.

6 Discussion

We succeed in generating the global patterns of the dissolved nuclide concentrations. However, we underestimate radioactivities in the surface waters (Figs. 6a and 7a). This shortcoming may be due to the instant equilibration between the dissolved and the adsorbed phases (Henderson et al., 1999). After the decay of U into ²³¹Pa and ²³⁰Th, the tracers instantly adsorb onto the particles in the surface ocean and are exported from the mixed layer. The equilibrium scavenging approach is commonly used in Pa∕Th modelling (Dutay et al., 2009; Gu and Liu, 2017; Siddall et al., 2005). It allows the model to use partition coefficients that can be directly constrained by observations (although consensus values for these coefficients are still not available). A non-equilibrium approach to scavenging could be considered, but it would necessitate the use of adsorption and desorption coefficients (k₊, k₋) that are currently very poorly constrained by observations. Observations shows equilibrium partition coefficient K values varying (decreasing) with particle flux (Chase et al., 2002; Roy-Barman et al., 1996), and we apply these in our simulation (i.e. lower K values for large, fast-sinking particles). This parameterisation produces effects that are similar to those of a non-equilibrium system for fast-sinking particles.

The reversible scavenging model uses partition coefficients K_ij for the different isotopes i and particles j. Many studies provide constraints on these coefficients but estimates vary strongly between different studies. Small particles are the most important scavengers and largely control the shape of the vertical profile of the tracers (Dutay et al., 2009). This also shows from sedimentary isotope records (Kretschmer et al., 2011). This is explained by the specific surface of small particles, which is larger than that of big particles. Therefore, higher adsorption values are set for small particles compared to the bigger particles of the same type (namely, 5 times higher for POC and 10 times higher for lithogenic particles). Compared to Dutay et al. (2009), we succeed in simulating the tracer concentrations using reasonable K values. This is due to improvement in (especially small) particle concentrations, which were generated with the lability parameterisation of POC (Aumont et al., 2017), the change in the CaCO₃ dissolution parameterisation and the addition of lithogenic dust in our ocean model. This shows that at least two particle size classes are needed to simulate dissolved and particulate thorium and protactinium. Our model does not include small particles for all types (there are neither small CaCO₃ nor small bSiO₂ particles), even though we expect that better results can be reached upon adding those particles akin to small POC and small lithogenic particles.

Figure 11 Sedimented ²³¹Pa∕²³⁰Th ratio when we set $K_{\mathrm{Pa},{\mathrm{bSiO}}_{\mathrm{2}}}/K_{\mathrm{Th},{\mathrm{bSiO}}_{\mathrm{2}}}=\mathrm{0.4}$. Top core measurements are presented as discs.

Figure 12 Dissolved protactinium-231 (a) and thorium-230 (b) at the Drake Passage (GIPY5) and western Atlantic (GA02) transects (mBq m⁻³); observations as discs. The Atlantic OSF, in Sv, is superimposed as a contour.

On the GA02 transect, the correlation coefficient r of the data–model comparison of [²³⁰Th]_D is 0.80, which is close to what was reported by Rempfer et al. (2017) in their simulations without bottom-boundary scavenging. On the GA03 section, r=0.78. The correlations between the model and observations are presented as scatter plots (Supplement Fig. S2).

There are a couple of interesting features to note on the comparisons. Firstly, the hydrothermal plume over the ridge results in ²³⁰Th-poor seawater that is not reproduced by the model (red triangle in Fig. S2, GA03). Secondly, in the western boundary currents, very low dissolved ²³⁰Th concentrations are observed due to strong nepheloids and/or possibly strong deep ventilation from the north. Neither of these processes is represented in the model. The lack of nepheloids has an especially large impact on the GA02 transect, which has a relatively large number of measured profiles near the western boundary. Finally, between 32 and 20^∘ W, the model tends to underestimate ²³⁰Th. The value of the correlation index is not that important because the goal of this study is not to give a perfect model–data comparison.

The provided dust flux does not distinguish between different dust particle sizes or types. Furthermore, dust deposition is uncertain; therefore, we chose a dust flux that was available to the ORCA2 configuration of NEMO, and its biogeochemistry has been tested with that dust flux. However, given the significant impact of dust particles on thorium and protactinium concentrations, it would be useful to look into the effect of different dust deposition fields and using different dust particle size classes. A good model of ²³⁰Th (and other isotopes of thorium) could even constrain the sporadic and uncertain dust deposition, but all this falls outside the scope of the present study.

Our model adds 80 % of the dust deposition into the fast-sinking lithogenic particle compartment (corresponding to ∅>51 µm), and only 20 % goes into the slowly sinking compartment (corresponding to ∅≤51 µm). This is necessary to reproduce the distributions of lithogenic particles. However, airborne dust particles typically have a diameter smaller than 20 µm (Knippertz and Stuut, 2014) and thus fall clearly within the small size class of less than 50 µm. Recently, there have been observations of aerosols larger than 20 µm (van der Does et al., 2016), but they do not reach far enough into the atmosphere above the open ocean to explain the high concentration of big lithogenic ocean particles. The explanation for the apparently required high fraction of big lithogenic particles is that smaller aerosols aggregate in the upper layers of the ocean. Moreover, though more hypothetically, there may be strong aggregation with biogenic particles just below the surface, below which the aggregates partly disaggregate again to result in the ∼0.2 big lithogenic fraction in the deep ocean (Fig. 5f).

The underestimation of lithogenic particle concentrations at the western boundary is expected. The reason is that we only have dust deposition as a source of lithogenic particles, whereas nepheloid layers are not included. Nepheloid layers are at least locally an important source of lithogenic (and biogenic) particles (Lam et al., 2015). With the transport of ²³⁰Th and ²³¹Pa through the western boundary current, a significant portion may be scavenged. Moreover, the periodic transport through the North Atlantic Gyre would result in lower [²³⁰Th]_D and [²³¹Pa]_D throughout a large part of the North Atlantic Ocean, possibly improving the concentrations in the deep ocean. Therefore, it will be useful to include nepheloid layers in the future. We have not done this so far because we do not know how to model nepheloid layers. Except for trivial models, like the one of Rempfer et al. (2017), who forced an additional constant scavenging rate in the bottom box of the ocean model, no large-scale prognostic nepheloid models have been developed.

We overestimate the radionuclide activity in the deep ocean, which is partly because we underestimate particle concentrations. Thorium adsorbs especially well onto lithogenic particles, which are underestimated in the western Atlantic Ocean (in the subsurface, west of about 43^∘ W, Fig. 5e). Below 2 km of depth, POC is also underestimated (Fig. 5a), even though the lability parameterisation improved the distribution by over an order of magnitude compared to previous versions of PISCES (Aumont et al., 2017). Small POC and the small lithogenic particles are the only small particles. The overall underestimation of small particles in the deep ocean results in an overestimation of the radionuclides in the deep ocean. In our simulation, [²³⁰Th]_D is not that much overestimated (outside the Arctic Ocean). However, [²³¹Pa]_D is strongly overestimated in the deep Pacific and Atlantic Ocean.

Clearly the model does not remove ²³¹Pa efficiently enough from the deep ocean (Figs. 7c, d, 8b and 12a). Two likely reasons may be that the waters are not well ventilated in the model (older than in reality) or that there is not enough scavenging. The Atlantic OSF (Fig. S1) shows that the upper overturning cell is too shallow compared to observational studies, so there is too little ventilation. This means that ²³¹Pa can build up in the deep water due to sluggish ventilation of the physical model (Biastoch et al., 2008; Dutay et al., 2002). Moreover, the AABW has weakened in the deep North Atlantic Ocean near the high [²³¹Pa]_D region (Fig. 12a), which may also contribute to the high dissolved ²³¹Pa concentrations. The volume transport of the lower cell is about 6 Sv near 20^∘ N and still 2 Sv near 40^∘ N, which is in the order of magnitude of what is reported by the literature. Moreover, the relationship between sediment ²³¹Pa∕²³⁰Th at any given site and the overturning circulation is very complex. In general, any increased opportunity for protactinium to scavenge would increase the sedimentary ²³¹Pa∕²³⁰Th ratio. Specifically, a slower overturning would mean more Pa scavenging and a higher ratio “just downstream of the sites of deep water formation” (Luo et al., 2010). Similarly, relatively stronger scavenging of Pa means a higher sedimentary ²³¹Pa∕²³⁰Th, which we can see in the Southern Ocean and some coastal areas where the concentration of bSiO₂ is high. For these reasons, it is not obvious whether the weak overturning is large enough to explain the discrepancy between the modelled and the observed [²³¹Pa]_D.

In reality, there are many regions in the North Atlantic Ocean where there is sediment resuspension and where there are nepheloid layers (Gardner and Sullivan, 1981; Gardner et al., 2017; Lam et al., 2015). This is consistent with the fact that our model underestimates lithogenic particle concentrations in the western Atlantic Ocean (Fig. 5e). The additional lithogenic particles would scavenge more ²³¹Pa, resulting in lower [²³¹Pa]_D. Even though the enforced scavenging occurs near the floor and the western boundary of the ocean, a strong flux of water passes through the latter region and is transported through the North Atlantic Gyre and through the rest of the Atlantic, diluting high tracer concentrations. The small POC and lithogenic particles are now underestimated in much of the deep ocean. They would contribute to lower radio tracer concentrations as well. Rempfer et al. (2017) confirmed that an additional bottom sink affects ²³¹Pa and ²³⁰Th significantly. They used a homogeneous extra scavenger in their model grid's bottom grid cell. Therefore, it would be worth testing if a realistic distribution of nepheloids would result in the right amount of scavenging.

Other particles may include resuspended (nepheloidal) biogenic particles, but also manganese and iron hydroxides that are not included in our model, and a smaller class of calcium carbonate and biogenic silica particles in addition to the big particle classes that are already in the model. Manganese oxides are especially available near hydrothermal vents, but also throughout the Pacific Ocean in low concentrations but much higher than in the Atlantic Ocean. Indeed, it has been argued that hydrothermal inputs may provide additional removal of ²³¹Pa and ²³⁰Th (German et al., 2016; Lopez et al., 2015; Rutgers van der Loeff et al., 2016), and recently this has been confirmed based on an analysis of the ∼10^∘ S zonal GP16 transect in the Pacific Ocean (Pavia et al., 2017).

Finally, the particle concentration effect (Honeyman et al., 1988) should be taken into account in future models, particularly when regions of high particle concentration (e.g. ocean margins and nepheloid layers) are included.

6.1 Effect of scavenging by bSiO₂ on the ²³¹Pa∕²³⁰Th sedimentation flux ratio

Previously, it was suggested that Pa has a stronger affinity for bSiO₂ than Th, i.e. $F_{\mathrm{Pa}/\mathrm{Th},{\mathrm{bSiO}}_{\mathrm{2}}}>\mathrm{1}$ (Chase et al., 2002; Dutay et al., 2009; Henderson et al., 1999). This “standard view” has weakened over time, so it is often accepted that $F_{\mathrm{Pa}/\mathrm{Th},{\mathrm{bSiO}}_{\mathrm{2}}}\mathit{\gtrsim }\mathrm{1}$ (Rempfer et al., 2017; Rutgers van der Loeff et al., 2016; Venchiarutti et al., 2011). We argue here that it is Th which has the stronger affinity to bSiO₂ ($F_{\mathrm{Pa}/\mathrm{Th},{\mathrm{bSiO}}_{\mathrm{2}}}<\mathrm{1}$), similarly (though not quite) to other particles.

Our standard simulation (with $F_{\mathrm{Pa}/\mathrm{Th},{\mathrm{bSiO}}_{\mathrm{2}}}=\mathrm{1}$) yields a too-high sedimentation ²³¹Pa∕²³⁰Th ratio compared to top core sediment observations below diatom-rich areas such as the Southern Ocean. In that region, bSiO₂ largely controls the particle scavenging of the two tracers, especially that of ²³¹Pa (Table 4). In order to estimate how strongly the model depends on the value of the bSiO₂ partition coefficient we performed a sensitivity experiment in which we set $K_{\mathrm{Pa},{\mathrm{bSiO}}_{\mathrm{2}}}=\mathrm{0.4}\mathrm{Mg}{\mathrm{g}}^{-\mathrm{1}}$ and $K_{\mathrm{Th},{\mathrm{bSiO}}_{\mathrm{2}}}=\mathrm{1.0}\mathrm{Mg}{\mathrm{g}}^{-\mathrm{1}}$ such that now $F_{\mathrm{Pa}/\mathrm{Th},{\mathrm{bSiO}}_{\mathrm{2}}}=\mathrm{0.4}$ (Table 2).

As a result of this reduction, the ²³¹Pa∕²³⁰Th flux ratio is no longer that much overestimated in diatom-rich regions (Fig. 11). This result suggests that protactinium has a weaker affinity to biogenic silica than thorium, though the fractionation is closer to 1 than that of other particle types (Table 2). This result is consistent with Geibert and Usbeck (2004) whose laboratory study shows no definitive affinity of either Th or Pa to biogenic silica. Their average K values actually suggest that thorium has a higher affinity to biogenic silica.

We included lithogenic particles in the model from deposited dust. Since lithogenic particles have a strong relative affinity with ²³⁰Th ($K_{\mathrm{Pa},\mathrm{Lith}}/K_{\mathrm{Th},\mathrm{Lith}}=\mathrm{0.2}$) and since they are mostly prevalent in the (northern) Atlantic Ocean, much of the ²³⁰Th is scavenged in the Atlantic Ocean. Calcium carbonate K_Pa∕K_Th is even smaller (0.024) and is present throughout much of the Atlantic but not in the Southern Ocean. This leaves protactinium to be scavenged by biogenic silica in the Southern Ocean when it arrives there through the AMOC. Therefore, even with the smaller-than-conventional K ratio of 0.4 for biogenic silica, the familiar pattern of higher values of ²³¹Pa∕²³⁰Th in the Southern Ocean compared to most of the rest of the ocean is still reproduced.

7 Conclusions

The purpose of this study is two-fold. Firstly, we set out a model of ²³¹Pa and ²³⁰Th, complete with all the necessary particle improvements and additions. This includes improved underlying particle dynamics from NEMO–PISCES and the addition of lithogenic particles. We have succeeded in this, and the model has been presented and implemented such that it can be used for future studies, including for the modelling and validation of ocean circulations in the past, present and future.

Secondly, this model helps us to study the interplay between the particle and water transport that controls the ²³¹Pa and ²³⁰Th profiles. We have shown that, by improving on the particles, we improved the radionuclide distributions and fluxes. We have also shown to what extent the different particle phases drive the scavenging of ²³¹Pa and ²³⁰Th.

Dissolved ²³⁰Th and ²³¹Pa concentrations are realistic in the intermediate-depth ocean (big improvement compared to Dutay et al., 2009) but not in the upper part of the ocean. A finite equilibration time between the different radionuclide phases may help. The model can be extended to include this in future versions. It would double the number of adsorption and desorption parameters whose values should be determined through a careful literature and model sensitivity study.

In the deep ocean, the overestimation of ²³⁰Th and ²³¹Pa in several regions is likely to be caused, at least partly, by missing nepheloid particles, but manganese oxides (from hydrothermal vents) may also improve the distributions. Of course, reported adsorption and desorption parameters vary a lot between different studies, meaning that a different, yet still plausible, combination of partition coefficients could better describe the available data. Additionally, the circulation may be too weak, not freshening AABW fast enough, and hence ²³⁰Th and ²³¹Pa concentrations become too high near the ocean floor.