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

In this paper we set forth a 3-D ocean model of the radioactive trace isotopes Th-230 and Pa-231. The interest arises from the fact that these isotopes are extensively used in paleoceanography to reconstruct ocean circulation. The tracers are reversibly scavenged by biogenic and lithogenic particles. In most previous models the complexity of particle dynamics is oversimplified and fails to do rigorous inferences about ocean circulation and other aspects of the past ocean. Our model is more complex, without underdetermining it more than needed, so that it can be used with more confidence than many other models. Our simulations of Th-230 and Pa-231 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 were 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 up to today, there are no non-trivial, prognostic models available on a global scale. Our simulations reproduce Th-230 and Pa-231 dissolved concentrations: they compare well with recent GEOTRACES observations in many parts of the ocean. Particulate Th-230 and Pa-231 concentrations are significantly improved compared to previous studies, but they are still too low because of missing particles from nepheloid layers.


Introduction
Oceanic circulation and 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 is transfer into the deep ocean through particle sinking. Because of this, the strength of the MOC and particle scavenging 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 unique large data set 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 from sediment cores are used as proxies to infer past ocean circulation. Several examples include carbon isotopes, the cadmium/calcium ratio (Rosenthal et al., 1997), the ratio between protactinium-231 and thorium-230 ( 230 Th/ 231 Pa), and the neodymium isotope ratio 143 Nd/ 144 Nd (Lacan and Jeandel, 2005;Jones et al., 2008;Arsouze et al., 2009;Ayache et al., 2016). These proxies are affected by dynamical and biogeochemical processes. Including these proxies in a climate model is a way to better understand the signal they register.
We will focus on 230 Th and 231 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 arXiv:1708.04157v1 [physics.ao-ph] 14 Aug 2017 2 M.M.P. van Hulten et al.: Global ocean model of 230 Th and 231 Pa ocean is perfectly known: radioactive decay of uranium isotopes. The ratio 230 Th and 231 Pa is also 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 the underlying particle dynamics is essential in order to better simulate these tracers.
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, 231 Pa and 230 Th are produced at a relatively constant rate (β Pa = 2.33 × 10 −3 dpm m −3 a −1 and β Th = 2.52 × 10 −2 dpm m −3 a −1 , dpm = disintegrations per minute) (Henderson and Anderson, 2003). They are both scavenged rapidly by the many particles that reside in the ocean. 231 Pa is less sensitive to particle scavenging than 230 Th, which is reflected in the longer residence time of 231 Pa (80-200 a) compared to that of 230 Th (20-40 a) (Yu et al., 1996). 231 Pa and 230 Th are radioactive, decaying to radium isotopes and having a half-life of 32.76 ka and 75.40 ka, respectively. Each combination of particleradionuclide adsorption has a different reactivity. The vertical distributions of natural radionuclides, such as 230 Th and 231 Pa, are hence sensitive to the distribution and mixture of particles. As a consequence of the different particle reactivities of 230 Th and 231 Pa, the [ 231 Pa] D /[ 230 Th] D deviates from the production activity ratio of 0.093 (Anderson, 2003;Rutgers van der Loeff et al., 2016).
As both particle dynamics and circulation of the ocean affect 230 Th and 231 Pa, numerical biogeochemical general circulation models are used to study the processes. The isotopes have been simulated in models of intermediate complexity, for instance by Marchal et al. (2000) (EMIC 2.5D), Siddall et al. (2007) Heinze et al. (2006) (HAMOCC) and Luo et al. (2010). To try to capture all important processes, more complex models have been developed as well, namely by Dutay et al. (2009) (NEMO-PISCES) and Rempfer et al. (2017) (Bern3D). Dutay et al. (2009) demonstrated that the simulated particle concentration in the deep ocean was much too low. This lead to overestimated radionuclide concentrations in the deep ocean. Therefore, it is crucial to improve the representation of the particles . Rempfer et al. (2017) showed that taking into account additional sinks at the seafloor and at the ocean margins yields an improved agreement with observations, especially for the dissolved phases of 230 Th and 231 Pa. Particulate ratios improved to a lesser extend, and the authors have not presented particulate concentrations, suggesting that the particle issue has not been satisfactory resolved.
In this study we use an ocean general circulation model to simulate the circulation of 230 Th and 231 Pa, and a biogeochemical model with improved particles for the scavenging. Aumont et al. (2017) showed that dissolution rates of Particulate Organic Carbon (POC) were overestimated, and they improved on this by introducing a spectrum of different la-bilities. This improved the simulation of both small and large POC significantly, so we use this same model for our simulations. Furthermore, we develop a model of lithogenic particles that also scavenge 230 Th and 231 Pa. The main objective of this paper is to improve on the simulation of 230 Th and 231 Pa (the dissolved and two particulate size classes of particles for both nuclides), based on a more realistic 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 230 Th and 231 Pa have been measured (Hayes et al., 2015b), but also their adsorbed forms (Hayes et al., 2015b) and biogenic and lithogenic particle concentrations in two size classes (Lam et al., 2015).

Model description
In order to simulate the biogenic particle dynamics and its interaction with the 230 Th and 231 Pa trace elements, we use the biogeochemical circulation model NEMO-PISCES (Madec, 2008;Aumont et al., 2015). This model has been employed for many other studies concerning trace metals, as well as large-scale ocean biogeochemistry (e.g. Gehlen et al., 2007;Arsouze et al., 2009;Dutay et al., 2009;Tagliabue et al., 2010;Van Hulten et al., 2013. We force PISCES by a climatological year of circulation fields (including turbulent diffusion) that was obtained from the dynamical component of the Nucleus for European Modelling of the Ocean (NEMO). Table 1 gives an overview of this and other components of the model and their relevant properties.
All model fields are defined on the ORCA2 discrete coordinate system, an irregular grid covering the whole world ocean with a nominal resolution of 2 • ×2 • , with an increased resolution in the 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. 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).

Circulation
The circulation was obtained by forcing NEMO in the ORCA2 configuration with air-sea boundary conditions, consisting of heat, fresh water 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  Table 1. Model components and properties essential for the radionuclide model. re-analysis data (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 simulation is used as our one-year climatology with a resolution of five days of the dynamics.

Particle dynamics
This version of PISCES includes two size classes of POC, both with differential remineralisation rates (Aumont et al., 2017), one class of biogenic silica (bSiO 2 ) and one class of calcium carbonate (CaCO 3 ). In the model, particles sink down with two velocities: w s = 2 m d −1 and w b = 50 m d −1 , corresponding with "small" and "big" particles ( Fig. 1).
For CaCO 3 we changed the standard first-order dissolution kinetics parameterisation to a fourth-order dissolution. For our simulation we used a calcite dissolution rate constant of k = 2.5 and a dissolution order of n = 3.9 (Subhas et al., 2015). Findings and discussion on this may be found in the report included in the Supplement of this publication, or to be downloaded at http://klimato.org/pub/report-calcite.pdf.
In addition to biogenic particles, we introduced lithogenic dust particles in the model. The yearly average dust flux is presented in Fig. 2, and is used by the model as the input of lithogenic ocean particles as well as for nutrient supply in PISCES.
Small and large lithogenic particles are added to the upper layer of the ocean according to: where P Litho is the small (large) lithogenic particle concentration, f is the fraction of the dust that gets partitioned into the small (large) lithogenic particles in the ocean, ∆z 1 = 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 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 where w is the settling velocity, set to the constant 2 m d −1 for the small lithogenic particles and to 50 m d −1 for the big particles. ∇ h is the horizontal divergence, and A and B are respectively the horizontal and vertical eddy diffusivity coefficients. The material derivative includes a term for eddyinduced velocity (Gent and McWilliams, 1990;Gent et al., 1995). Equation 2 is identical to the settling of small and big POC (Aumont et al., 2015). Of course there are biological and chemical sources and sinks as well 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 bathymetry, which means that they are buried in the sediment.

Radionuclides
Thorium-230 and protactinium-231 are produced throughout the ocean from the decay of 234 U and 235 U, respectively. The concentrations of these uranium isotopes are homogeneously distributed throughout the ocean, and do not change much over time. Therefore, the production decay rates are constant, both in space and time.
The 230 Th and 231 Pa radionuclides are reversibly scavenged by biogenic and lithogenic particles (Fig. 3). We will assume, as in previous studies (e.g. Siddall et al., 2005;Dutay et al., 2009), that the adsorption and desorption reaction rate is much faster than radionuclide production, decay, advection, mixing, change in particle distribution and settling of the adsorbed phases. Because of that we equilibrate between the dissolved and adsorbed phases instantly, which means we  must solve this set of equations for every nuclide i:  stands for the concentration of particle j ∈ {sPOC, bPOC, CaCO 3 , bSiO 2 , sLitho, bLitho}. D is the set of non-sinking phases (here only dissolved), S of small particles that sink with 2 m d −1 and B of big particles that sink with 50 m d −1 (in the same way as Eq. (2) for lithogenic particles). For any particle set S and any radionuclide i, by definition, A S i = j∈S A i,j . K i,j is the equilibrium partition coefficient of nuclide i for particle j.
Since we cannot solve this analytically, we solve this numerically by first assigning a new value to the dissolved nu-These terms are used interchangeably throughout this paper. Moreover, they are considered identical: the activity or concentration  Table 2. Partition coefficients for the different modelled particles. Between braces is the value for the sensitivity simulation discussed in Section 6.1. clide activity: 2 Then we calculate the activity of i that is adsorbed onto small and big particles by applying Eqs (3b) and (3c). With this approach the total activity of every i is conserved, i.e. Eqn (3a) holds, and the small and big adsorbed concentrations equilibrate instantly.
The 230 Th and 231 Pa that are adsorbed onto the particles follow the same law as the small and big (lithogenic) particles (Eq. 2). Of course, by definition, the radio isotope and the particle settle with the same speed, and thus we have implemented it.
The decay terms of 230 Th and 231 Pa are much smaller than the other sources and sinks, but they are included in the model: where β i and λ i are respectively the production rate and radioactive decay of isotope i.

Simulations
The model was spun up for 400 a, after which it was in a steady state (yearly drift of −0.01 % for total 230 Th and −0.09 % for total 231 Pa).
The partition coefficients K i,j for i ∈ {Th, Pa}, depend on the type of particle j and are given in Table 2. Uncertainties are large in the value of these coefficients, which can be seen in the large differences between reported values from different studies (e.g. Chase et al., 2002;Dutay et al., 2009;Geibert and Usbeck, 2004;Hayes et al., 2015b;Luo et al., 2 When summation bounds are not specified, the union of all phases, D ∪ S ∪ B, is assumed. 2010). Therefore we had quite some freedom in the choosing values that yielded good results and at the same time being consistent with the overall range given by the literature. The adsorption onto calcium carbonate is a factor two decreased from Chase et al. (2002). This brings the K value of Pa closer to that found by Hayes et al. (2015b) based on field measurements (0.9 ± 0.4 Mg g −1 ), but since we maintained the ratio of K values, the partition coefficient of Th is much smaller than found by Hayes et al. (2015b). However, it is consistent with the laboratory results of Geibert and Usbeck (2004). Small lithogenic particles K is set about a factor five larger than literature (Geibert and Usbeck, 2004;Hayes et al., 2015b), whereas large lithogenic particles have a smaller value than reported in the literature. The consequences of chosen values of the partition coefficients is discussed further in Section 6.

Observations
In this study we will focus on the GEOTRACES GA03 transect in the North Atlantic Ocean. Recently a large amount of measurements on both radionuclides and their carrier particles have been collected at 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).
For the carrier particles most of our data comes 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) Table 3.
To compare the sediment 231 Pa/ 230 Th flux of the model, we will compare with compilations of The Holocene (i.e. top core particulate concentrations) (Table 3).

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 big importance, we here present basic results thereof.  (through the basin) and vertically (from the surface downwards) integrated meridional current speed. We use this as a measure for the Atlantic Meridional Overturning Circulation (AMOC). The upper overturning cell has a strength of about 14 Sv (1 Sv = 10 6 m 3 s −1 ). Observations suggest higher values of about 19±5 Sv (Talley et al., 2003;Cunningham et al., 2007;Rayner et al., 2011). The lower cell has an overturning strength of about 6 Sv, which is high but mostly due to the fact that the AMOC is shallow (back flow is at around 2 km depth), leaving a lot of depth for AABW (and mixtures with northern Atlantic water masses). At around 30 • N the AMOC has weakened notably, which is consistent with some studies (e.g. Johnson, 2008), but the AntArctic Bottom Water (AABW) does not go as far north as other studies suggest (e.g. 45 • N in Van Aken, 2011). Figure 5 shows the surface concentrations of total POC, large calcium carbonate, large biogenic silica and total lithogenic particles. Observations are plotted as dots on top of the modelled concentrations. The modelled biogenic particle concentrations include living material. For the model we assume POC includes all phytoplankton and microzooplankton (not mesozooplankton since that may swim away), calcium carbonate contains the assumed fraction of CaCO 3 in the phytoplankton, and opal includes living diatoms. The modelled concentrations and patterns of total POC compare well with the observations (Fig. 5a). For instance, coastal and equatorial Pacific POC concentrations are both in the model and the observations elevated compared to other regions, albeit that the spatial extend of the oligotrophic regions appears to be too small in the model. We present large CaCO 3 particles, because we have data of that for both the Atlantic and Pacific Ocean, whereas we only have small particle data from the Atlantic Ocean (and none from the Indian Ocean). Furthermore, our model only includes large (fast-sinking) calcium carbonate particles. As the model tries to represent CaCO 3 with only one size class, neither comparing with only big particles nor comparing with total CaCO 3 is completely fair. Whereas the meridional patterns of large CaCO 3 particles are reproduced, its concentrations are generally overestimated, especially in the Gulf of Alaska (Fig. 5b). Contrarily, total observed [CaCO 3 ], which includes smaller particles, of which measurements are available only at GA03, is higher than the prediction of our modelled (only large) CaCO 3 particle concentration. The model produces realistic 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. 5c).

Particles
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 the biogenic particles ( Fig. 6a-d). However, for all four biogenic particles (small and big POC, calcium carbonate and biogenic silica) there are shortcomings in their distributions as well.
-The total POC concentration in Fig. 6a 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. 5a). In both the observations and the model the fraction of large POC varies from zero to 0.6 ( Fig. 6b) 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 3 concentrations from 70 • W 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) observational large CaCO 3 is again a factor of up to 10 larger than what the model predicts.
-Biogenic silica concentrations are generally reproduced except for the western margin where, while also high in the observations, the model overestimates the observed biogenic silica concentration (Fig. 6d).
Finally, the modelled total lithogenic particle concentration shows, as expected, a close resemblance with dust deposition patterns (Fig. 2), and mostly compares well with observations (Fig. 5d). Only near the coasts of the USA and of Portugal the model strongly underestimates lithogenic particle concentrations. Concerning the deep ocean, our model captures quite nicely the general distribution at the GA03 transect (Figs. 6e and 6f). However, at the western boundary concentrations are strongly underestimated, especially that of large lithogenic particles. Our model reproduces the observed magnitude of ∼0.1 to about ∼0.3 in most of the deep ocean, but it does not capture the much larger fraction of big particles in the upper 200 m of the ocean (0.1-0.2 in the model versus 0.3-0.9 in the observations).    Fig. 9a for the GA03 transect) that generates realistic values for the concentration of this element in the dissolved phase. 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 respectively nepheloid layers and manganese oxides from hydrothermal vents, which are not produced in the simulation (Figs 9a and 9b).

Thorium-230 and protactinium-231
In the intermediate and deep waters of all the oceans, modelled [ 231 Pa] D is generally of the correct order of magnitude but is strongly overestimated below 2500 m depth (Fig. 9b), especially in the Pacific Ocean (Fig. 8c,d).
The Contrarily, small 231 Pa particles are overestimated in our model below the surface at the GA03 transect (Fig. 9d). Large 231 Pa particles are simulated more realistically (Fig. 9f) phase. Table 4 gives two types of information on the adsorbed phases. One is the fraction of radionuclide carried by each phase ("global amount"). 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. The amount and the flux are not necessary the same, because the big particles have a different settling velocity from the small particles. Only a small amount of 230 Th is adsorbed onto biogenic silica (3 %), though the amount is larger for 231 Pa (11 %). However, the global settling flux of 231 Pa is primarily (54 %) determined by bSiO 2 . Most of the modelled particulate 231 Pa and 230 Th is on the lithogenic particles, but lithogenic particles account for only 14 % of the 231 Pa flux and 9 % of the 230 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 3 : only 24 % of the Th is on calcium carbonate, but CaCO 3 is for 76 % responsible for the 230 Th export.
The relative contributions of the different particulate phases also vary for different regions of the ocean. Opal is for 94 % responsible for Pa export in the Southern Ocean.
Compared with previous studies, specifically Dutay et al. (2009), the improvement of the 230 Th profile shape is only to a small extend due to the large improvement in the POC representation (Aumont et al., 2017). It is mostly due to the improvements in the addition of lithogenic dust particles and the adjusted CaCO 3 dissolution. Lithogenic particles and CaCO 3 hardly dissolve compared to POC and hence the downward particle flux remains constant, which is in agreement with the 1-D 230 Th linear profile model (e.g. Roy-Barman, 2009).

Sedimentation flux ratio of 231 Pa/ 230 Th
The sedimentation flux of the total adsorbed 231 Pa/ 230 Th ratio captures the general patterns of the upper sediment core 231 Pa/ 230 Th concentration ratio (Fig. 10). For instance, the model produces the more elevated 231 Pa/ 230 Th ratio values observed in the Southern Ocean, the North Indian and Pacific Ocean and along the coastal upwelling regions. However, at some places, like most of the Southern Ocean, our model tends to overestimate the ratio derived from the sediment core observations.

Discussion
We use a reversible scavenging model at equilibrium for our simulation, this approach succeeds to generate globally realistic results. However, we underestimate radioactivities at the surface (Figs 8a and 7a). This shortcoming may be due to the instant equilibration between the dissolved and the adsorbed phases. After production of 231 Pa and 230 Th, the tracers instantly adsorb onto the particles in the surface ocean and are exported from the mixed layer. For future model developments, it is of interest to test if a non-instant equilibration model (with a k + and k − ) may yield better surface values.  The reversible scavenging model used partition coefficient values K i,j for the different isotopes i and particles j. Many studies provide constraint for the values of these coefficients but estimations vary strongly between different studies. Small particles are the most important scavengers and largely control the shape of the vertical profile of the tracers ). This is explained by the specific surface of small particles that 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 to simulate realistic tracer concentration using reasonable K values. This is due to improvement in small particle concentrations, that were generated with the lability parametrisation of POC (Aumont et al., 2017) and the addition of lithogenic dust in our model.
The provided dust flux does not distinguish between different dust particle sizes or types. Our model adds 80 % of the dust deposition into the large lithogenic particle compartment, and only 20 % goes into the small compartment. This is necessary to reproduce the distributions of lithogenic particles. However, airborne dust particles typically have a diameter of 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 Does et al., 2016), but they do not reach far enough into the atmosphere above the open ocean to explain the high concentration of large lithogenic ocean particles. The explanation for the apparently required high fraction of large 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 are partly disaggregated again to result in the ∼0.2 big lithogenic fraction in the deep ocean (Fig. 6f). 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 230 Th and 231 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 [ 230 Th] D and [ 231 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, there are no prognostic 3-D nepheloid models developed.
As a consequence of not having nepheloids, we overestimate the radionuclide concentrations. Especially thorium adsorbs well onto lithogenic particles. This probably plays a role in the too strong build-up of 230 Th in the intermediate and deep Arctic Ocean (Fig. 7). There could also be a significant effect in the intermediate and deep Atlantic Ocean. In our simulation, [ 230 Th] D is not that much overestimated (outside the Arctic Ocean). However, [ 231 Pa] D is strongly overestimated in the deep Pacific and Atlantic Oceans.
Clearly the model does not remove 231 Pa efficiently enough from the deep ocean (Figs 8c,d,9b and Fig. 12upper panel). Two likely reasons may be that either the waters are not well ventilated in the model (older than in reality), or that there is not enough scavenging. Looking in Fig. 12 (upper panel), we see that the Atlantic OSF (Fig. 4) has weakened near the high-[ 231 Pa] D region, i.e. in the deep North Atlantic Ocean, so it is possible that the high deep [ 231 Pa] D is due to the sluggish ventilation of the physical model. However, it is not obvious whether the weakening is large enough to explain the discrepancy between the modelled and the observed [ 231 Pa] D . It 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 literature. Therefore, at least in the Atlantic Ocean, it is uncertain if unrealistic circulation is to blame for the overestimation of [ 231 Pa] D in the deepest waters. In the actual ocean, there are many regions in the North Atlantic Ocean where there is sediment resuspension and where there are nepheloid layers. This is consistent with the fact that our model underestimates lithogenic particle concentrations in the West Atlantic Ocean (Fig. 6e). The additional lithogenic particles would scavenge more 231 Pa, resulting in lower [ 231 Pa] D . Of course, the enforced scavenging occurs mostly near the western boundary of the ocean, but a strong flux of water passes through this region, and is transported through the North Atlantic Gyre and through the rest of the Atlantic, diluting high tracer concentrations. Figure 11. Sedimented 231 Pa/ 230 Th ratio when we set K Pa,bSiO 2 /K Th,bSiO 2 = 0.4. Top core measurements are presented as dots. Rempfer et al. (2017) confirmed that an additional bottom sink affects 231 Pa and 230 Th significantly. They used a homogeneous extra scavenger in their model grid's bottom grid cell. Therefore, it would still be worth tested 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 large-particle classes that are already in the model. Manganese oxides are available especially 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 231 Pa and 230 Th (Rutgers van der Loeff et al., 2016).
6.1 Effect of bSiO 2 scavenging on the 231 Pa/ 230 Th sedimentation flux ratio Our simulation yields a too high sedimentation 231 Pa/ 230 Th ratio compared to top-core sediment observations at diatomrich places such the Southern Ocean. In that region bSiO 2 largely controls the particle scavenging of the two tracers, especially that of 231 Pa (Table 4). In order to estimate how strongly the model depends on the value of the bSiO 2 partition coefficient we performed a sensitivity experiment with a reduced K Pa,bSiO2 and increased K Th,bSiO2 . Previously it was suggested that Pa has a stronger affinity for bSiO 2 than Th has (e.g. Chase et al., 2002;Dutay et al., 2009). We argue here that Th has a significantly stronger affinity to bSiO 2 than Pa, similarly (though not quite) like other particles. Figure 11 shows the 231 Pa/ 230 Th flux ratio in the simulation where K Pa,bSiO2 = 0.4 Mg g −1 and K Th,bSiO2 = 1.0 Mg g −1 (Pa/Th ratio of 0.4). As a result of this reduction, the high-diatom 231 Pa/ 230 Th flux ratio is not that much over- estimated anymore. This result suggests that protactinium does not have a stronger affinity to biogenic silica than thorium does, though the fractionation is closer to one 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. Also Rempfer et al. (2017)  We included lithogenic particles in the model from deposited dust. Since they have a strong relative affinity with 230 Th (K Pa,Lith /K Th,Lith = 0.2) and since the lithogenic particles are mostly prevalent in the (northern) Atlantic Ocean, much of the 230 Th is scavenged in the Atlantic Ocean but it 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.

Code availability
We have provided the model equations, source code and forcing files such that our results can be reproduced. Moreover, we encourage the reader to use our model to build extended and improved models. We have introduced two particle size classes of adsorbed nuclides, which is the minimum that is needed to produce good results. Instead of introducing a new tracer for every particle type, we only distinguish between "adsorbed onto small particles" and "adsorbed onto big particles", meaning that one is restricted, in this set-up, to include only particles that sink with the two respective settling speeds for scavenging. Whereas this is somewhat restricting, it is computationally effective, and it restricts the number of degrees of freedom (which is usually a good thing in complex models). ProThorP can also be used to set up a non-instant equilibration model (with a k + and k − ) that may yield better surface values.
NEMO 3.6 can be downloaded from http://www. nemo-ocean.eu/ after creating a login. The radionuclide and lithogenic particles model, as well as adjustments to the NEMO code, can be obtained at http://klimato.org/pub/. One is advised to download the latest version of ProThorP (prothorp_ * .tar.gz).
The NEMO model is available under the CeCILL free software licence, modelled after the GNU GPL. The lithogenic particles and 230 Th/ 231 Pa-specific code is licensed under the same terms, or, at your option, under the GNU General Public License version 3 or higher. The authors do appreciate if, besides the legal adherence to copyleft and attribution, they are informed about the reuse of the code in any (upcoming) publications.

Conclusion
The purpose of this study is two-fold. Firstly, we wanted to set out a model of 231 Pa and 230 Th, complete with all the necessary particle improvements and additions. The latter includes an improved underlying particle dynamics from NEMO-PISCES as well as 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 circulation in the past.
The second purpose of this study was to show that particle dynamics and composition are very important in modelling 231 Pa and 230 Th. We have indeed shown that by improving on the particles, we improved the radionuclides, and we have shown to what extend the different particle phases drive the scavenging of 231 Pa and 230 Th.
Dissolved 230 Th and 231 Pa concentrations are realistic in intermediate-depth ocean (big improvement compared to Dutay et al. (2009)) but not in the surface: a finite equilibration time between Th and Pa phases may help. This could easily be implemented in the current model framework for in future versions of the model. It would double the number of adsorption/desorption parameters whose values should be determined through a careful literature and model sensitivity study.
In the deep ocean, overestimation at several regions is likely to be caused, at least partly, by missing nepheloid particles, but also manganese oxides (from hydrothermal vents) may improve the distributions. Additionally, the circulation may be too weak, not freshening AABW fast enough, and hence [ 230 Th] D and [ 231 Pa] D become too high near the ocean floor.