The potential of ²³⁰Th for detection of ocean acidification impacts on pelagic carbonate production

Concentrations of dissolved ²³⁰Th in the ocean water column increase with depth due to scavenging and downward particle flux. Due to the ²³⁰Th scavenging process, any change in the calcium carbonate (CaCO₃) fraction of the marine particle flux due to changes in biological CaCO₃ hard-shell production as a consequence of progressing ocean acidification would be reflected in the dissolved ²³⁰Th activity. Our prognostic simulations with a biogeochemical ocean general circulation model using different scenarios for the reduction of CaCO₃ production under ocean acidification and different greenhouse gas emission scenarios – the Representative Concentration Pathways (RCPs) 8.5 to 2.6 – reveal the potential for deep ²³⁰Th measurements to detect reduced CaCO₃ production at the sea surface. The time of emergence of an acidification-induced signal on dissolved ²³⁰Th is of the same order of magnitude as for alkalinity measurements. Interannual and decadal variability in factors other than a reduction in CaCO₃ hard-shell production may mask the ocean-acidification-induced signal in dissolved ²³⁰Th and make detection of the pure CaCO₃-induced signal more difficult so that only really strong changes in marine CaCO₃ export would be unambiguously identifiable soon. Nevertheless, the impacts of changes in CaCO₃ export production on marine ²³⁰Th are stronger than those for changes in POC (particulate organic carbon) or clay fluxes.

1 Introduction

Progressing ocean acidification is a fact. It can be directly seen from measurements at Eulerian time series stations (e.g. Bates, 2007; Dore et al., 2009; Santana-Casiano et al., 2007) and also at larger depth from high-quality deep hydrography data (Olafsson et al., 2009). Depending on the emission scenario for CO₂, the decrease in ocean pH and the decrease in carbonate saturation will continue and become more pressing during this century (e.g. Orr et al., 2005; Steinacher et al., 2009; Bopp et al., 2013). Temporal and regional patterns of changes in pH and carbonate saturation are relatively straightforwardly to project by Earth system models including marine inorganic carbon chemistry formulations (e.g. Bopp et al., 2013). They also can be monitored through long-term high-quality measurements of the inorganic carbon system. However, monitoring impacts of ocean acidification on biological processes remains challenging. While for some organisms, especially corals (Kleypas et al., 1999), detrimental effects due to decreasing pH and carbonate saturation have been determined, the various physiological impacts of ocean acidification on specific organisms and ecosystem functioning are still under investigation (e.g. Iglesias-Rodriguez et al., 2008; Kroeker et al., 2013; Meyer and Riebesell, 2015; Riebesell et al., 2007). At present, a series of possible pH dependencies of governing marine carbon cycle parameters (such as elemental stoichiometric ratios and nutrient uptake kinetics under biological particle production) are discussed. So far, the potential decrease in calcification due to the lowering of carbonate saturation under high pCO₂ is among the key changes which may be expected to occur, especially when it comes to organisms building aragonite shells (the metastable form of CaCO₃, calcite, has a lower solubility than aragonite) (Raven et al., 2005). Two important questions emerge: (1) if changes in biological calcification would indeed occur in the ocean during the coming years, how can they be detected and monitored by observational frameworks (through which methodology and through which variable)? (2) In which oceanic region could one observe early signals of these changes unambiguously at the earliest possible stage?

Development of such early warning systems is vital in order to check the validity of parameterisations of pH-dependent processes in ocean models and to take appropriate mitigation/adaptation measures to diminish the consequences of potential considerable ecosystem changes due to ocean acidification. Such changes could affect the marine food chain. In a global modelling study, Ilyina et al. (2009) quantified the detection thresholds for changes in alkalinity due to changes in a series of possible formulations for the reduction of calcification (with the term “calcification” we mean here the production of hard CaCO₃ parts by marine biota) with pH decrease. The logic behind the approach of Ilyina et al. (2009) is as follows: if ocean acidification leads to a decrease in calcification, more ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ ions would be retained in the surface water and not become incorporated into CaCO₃ shell material. Changes in ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ ion concentration would induce a change in ocean total alkalinity, which could eventually be measured against an earlier baseline. According to that study, the tropical Pacific (with high CaCO₃ production rates) would be the region for detecting such alkalinity changes first, as the anticipated change in CaCO₃ production would be largest there in absolute terms. For intermediate dependencies of CaCO₃ production on pH/carbonate saturation, a reduction in biological calcification could unequivocally only be diagnosed from approximately the year 2035 on. In the Arctic Ocean, where pH changes are expected to be most pronounced, large-scale changes in calcification would be detectable even later on the basis of alkalinity measurements due to the overall lower biological production rates. Likewise, detection may be additionally complicated due to signals from natural seasonal and interannual variations in surface total alkalinity (Carter et al., 2016). There is thus a need for the development of novel detection methods. Heinze et al. (2006) investigated the impact of varying rain ratios CaCO₃ : C_org on the distribution of the radionuclides thorium (²³⁰Th) and protactinium (²³¹Pa), as well as beryllium (¹⁰Be). As the rain ratio describes the average amount of carbon atoms incorporated into CaCO₃ shell material relative to the amount of carbon atoms incorporated into organic matter by plankton, a reduction in calcification would lead to a rain ratio decrease. Indeed, especially for ²³⁰Th, due to its affinity for scavenging to CaCO₃ particles, a considerable increase in the concentration of dissolved ²³⁰Th with depth and time was obtained in a sensitivity experiment with strong rain ratio reduction. In this paper, we explore the option of using radionuclides for diagnosing changes in calcification and respectively reduced downward marine particle flux in more detail.

Figure 1 Schematic illustration of the equilibration process between the dissolved and particle-attached phases of ²³⁰Th and the increasing concentrations downward in the water column.

2 The concept – radionuclides and particles

We focus here on ²³⁰Th, a long-lived radionuclide (half-life 7.5 × 10⁴ yr) and a highly particle reactive metal. ²³⁰Th is produced in the uranium decay series. As uranium has a very long residence time in the ocean and quasi-homogenous concentrations in seawater, the present natural marine ²³⁰Th source is constant everywhere in the ocean water column and thus ideally suited for using ²³⁰Th as a tracer in modelling studies. ²³⁰Th is strongly particle reactive and is removed from the water column quickly through scavenging by the vertical particle flux in the ocean (for a summary, please see Henderson et al., 1999). The majority of ²³⁰Th is removed from the water column and will not re-enter the bottom waters, e.g. through sediment pore water diffusion (though resuspension may potentially cause some re-release). In spite of its strong particle reactivity, concentrations of dissolved ²³⁰Th show horizontal as well as vertical gradients in the water column which are induced by the oceanic current field and differential particle fluxes as well as particle concentrations (e.g. Henderson et al., 1999; Yu et al., 1996). The distribution of particle-bound ²³⁰Th follows the particle concentrations and particle fluxes. Both the dissolved ²³⁰Th and the particle-bound ²³⁰Th show nutrient-type vertical profiles with low values at the surface and increasing concentrations with depth. Dissolved ²³⁰Th increases with depth as particles continuously carry ²³⁰Th downwards and new equilibria between the dissolved phase and the particle-attached phase establish as illustrated in Fig. 1. The equilibrium between the concentration of the dissolved phase [²³⁰Th_diss] and the concentration of the particle-attached phase [²³⁰Th_part] can be described in an analog way to a chemical reaction equation:

$${\displaystyle }{\displaystyle }{}^{\mathrm{230}}{\mathrm{Th}}_{\mathrm{diss}}\leftrightarrow {}^{\mathrm{230}}{\mathrm{Th}}_{\mathrm{part}}.$$

The respective analog for the mass action law constant describing to which extent the reaction from left to right is carried out is then given through the k_d value governing the equilibrium between the dissolved and particle-bound phases for the radionuclide:

$${\displaystyle }{\displaystyle }{k}_{\mathrm{d}}={\displaystyle \frac{\left[{}^{\mathrm{230}}{\mathrm{Th}}_{\mathrm{part}}\right]}{\left[{}^{\mathrm{230}}{\mathrm{Th}}_{\mathrm{diss}}\right]}}.$$

Often, k_d values are formulated to account for a shift of the equilibrium towards the particle-attached phase at low particle mass concentrations, assuming that low particulate concentrations occur in parallel to low particle sizes with correspondingly high reactive surface areas when compared to large particles (e.g. Honeyman et al., 1988). Respective formulations have been used successfully in ²³⁰Th modelling studies (Heinze et al., 2006; Henderson et al., 1999).

How would then the distribution of dissolved ²³⁰Th in the ocean reflect changes in marine calcification at the sea surface? There are indications that ²³⁰Th is attached first of all to CaCO₃ and clay particles in the water column as preferential carrier phases (Chase et al., 2002, 2003; Luo and Ku, 1999; Hayes et al., 2015b). Indeed, Heinze et al. (2006) reproduced the large-scale ²³⁰Th distribution in the world ocean with a coarse-resolution biogeochemical ocean general circulation model based on a formulation allowing ²³⁰Th scavenging only by CaCO₃ and clay particles. Other studies have indicated that preferential carrier phases may vary regionally in the ocean (Scholten et al., 2005). However, Heinze et al. (2006) could demonstrate that rain ratio changes (CaCO₃ : C_org) in marine biological particle export production could also be recorded, if ²³⁰Th became in addition attached to particulate organic carbon (POC). A weakening of CaCO₃ particle production would result in a less efficient ²³⁰Th scavenging as less particles (in terms of mass) would be available to carry ²³⁰Th to larger depths and finally the sediment. Changes in the strength of CaCO₃ production and the respective downward particle flux are reflected increasingly better in the ²³⁰Th distribution with increasing depth (see Heinze et al., 2006) due to two reasons. First, CaCO₃ particles get less well degraded as POC (which is remineralised mostly in the upper 1000 m of the water column) and thus reach larger depth; this is also corroborated from sediment trap measurements (e.g. Honjo, 1996). Second, due to the scavenging of ²³⁰Th by particles, the vertical downward particle flux, and the equilibration between dissolved and particle-bound phases, ²³⁰Th is continuously transferred from shallower layers to larger depths. Thus, temporal changes in ²³⁰Th scavenging in the upper ocean become enlarged in the deep ²³⁰Th distribution as in a kind of “magnifying glass” (see Fig. 1). We investigate here whether this feature can be exploited for an early detection method of large-scale reductions in calcification and correspondingly reduced rates in CaCO₃ particle export and CaCO₃ particle fluxes through the water column.

3 Model description

In this study, we use the Hamburg ocean carbon cycle circulation model, HAMOCC (Maier-Reimer, 1993), in its annually averaged version (time step of 1 year; Heinze and Maier-Reimer, 1999; Heinze et al., 2009, 2016) with a horizontal resolution of 3.5^∘ × 3.5^∘. This coarse-resolution model is computationally very efficient and useful when multiple integrations are needed, such as for the testing and adjusting of scavenging codes. An advantage of this fast model is that it includes a fully equilibrated early diagenesis module (10 layers) under each grid point and thus can account for alkalinity changes induced by dissolution of CaCO₃ from the seafloor. The model version employed here corresponds to the version used in Heinze et al. (2009, 2016) with a slightly updated scavenging module of Heinze et al. (2006) with respect to the formulation of the equilibrium coefficient governing the distribution of ²³⁰Th between the dissolved and particle-attached phase (see below). For details, please see those publications. The model uses a fixed ocean velocity field (and thus provides no dynamical computation of the ocean currents; velocities, temperature, salinity, and ice cover are rather read from an input file). The model includes an atmospheric compartment (“slab atmosphere”) which allows for prognostic computation of the atmospheric CO₂ concentration as well as meridional atmospheric CO₂ transport. We describe here only briefly features of specific relevance for this study. The water column is structured into 11 layers (centred at 25, 75, 150, 250, 450, 700, 1000, 2000, 3000, 4000, and 5000 m). The bioturbated top sediment zone of the early diagenesis module is divided into 10 layers which are separated by interfaces at 0, 0.3, 0.6, 1.1, 1.6, 2.1, 3.1, 4.1, 5.1, 7.55, and 10 cm “downcore”. We make the simplifying assumption that no pore water reactions take place below 10 cm depth in the sediment (see, e.g. Smith and Rabouille, 2002; Boudreau, 1997). The biogeochemical model includes the processes of air–sea gas exchange, biogenic particle export production out of the ocean surface layer, particle flux through the water column and particle degradation by dissolution as well as remineralisation, transport of dissolved substances with the ocean currents, deposition of particulate constituents on the ocean floor, pore water chemistry and diffusion, advection of solid sediment weight fractions (organic carbon, organic phosphorus, CaCO₃, opal, and clay), bioturbation, and sediment accumulation (export out of the sediment mixed layer). The model predicts the following tracer concentrations in the atmosphere, the ocean water column, and in the sediments. Atmospheric tracers include the concentrations of CO₂ (carbon dioxide) and O₂. In the water column, concentrations of DIC (dissolved inorganic carbon), POC (particulate organic carbon), POP (particulate organic phosphorus), DOC (dissolved organic carbon), CaCO₃ (calcium carbonate or particulate inorganic carbon), dissolved oxygen O₂, dissolved ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ as the biolimiting nutrient, silicic acid Si(OH)₄, and opal (biogenic particulate silica – BSi) are calculated. In the sediment pore waters, the same dissolved substances as in the water column, as well as solid sediment constituents such as clay, CaCO₃, opal, and organic carbon, are simulated. The inorganic carbon chemistry is computed following Dickson et al. (2007). In the advection scheme and for the other chemical reactions, DIC and total alkalinity (TAlk) are used as “master tracers” form which derived quantities such as the ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentration and the pH value are computed through a Newton–Raphson algorithm. In the annually averaged model as employed in this study, only export production of biogenic particles is modelled (and no explicit phytoplankton and zooplankton concentrations). Particle production takes place in the model surface layer representing the euphotic zone. Phosphate serves as the biolimiting nutrient. POC and opal export productions are simulated following Michaelis–Menten kinetics for nutrient uptake (e.g. Sarmiento and Gruber, 2006) (where the phytoplankton concentration is replaced by the phosphate concentration as ecosystem processes as such are not explicitly modelled):

$${\displaystyle }{\displaystyle }{P}_{\mathrm{POC}}={\displaystyle \frac{V_{\max }^{\mathrm{POC}}\cdot {\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}^{\mathrm{2}}\cdot \mathrm{Red}(\mathrm{C}:\mathrm{P})}{K_{\mathrm{s}}^{\mathrm{POC}}+\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}},$$

and

$${\displaystyle }{\displaystyle }{P}_{\mathrm{opal}}={\displaystyle \frac{V_{\max }^{\mathrm{opal}}\cdot {\left[\mathrm{Si}(\mathrm{OH}{)}_{\mathrm{4}}\right]}^{\mathrm{2}}}{K_{\mathrm{s}}^{\mathrm{PopalC}}+\left[\mathrm{Si}(\mathrm{OH}{)}_{\mathrm{4}}\right]}},$$

where P_POC and P_opal are the POC and opal export production rates (mol L⁻¹ yr⁻¹), Red(C : P) is the Redfield ratio C : P, $V_{\max }^{\mathrm{POC}}$ and $V_{\max }^{\mathrm{opal}}$ are the maximum uptake rate of phosphate and silicic acid from the water column (yr⁻¹), and $K_{\mathrm{s}}^{\mathrm{POC}}$ as well as $K_{\mathrm{s}}^{\mathrm{opal}}$ are the respective half-saturation constants. $V_{\max }^{\mathrm{POC}}$, $V_{\max }^{\mathrm{opal}}$, $K_{\mathrm{s}}^{\mathrm{POC}}$, and $K_{\mathrm{s}}^{\mathrm{opal}}$ are simulated as a function of sea surface temperature as described by Heinze et al. (2003). POP production follows POC production with a constant stoichiometry here. The export production of CaCO₃ is coupled to the local production ratio P_opal ∕ P_POC. It starts to increase gradually (parameter R; see below) if P_opal ∕ P_POC sinks below a threshold value S_opal, i.e. when not enough silicic acid is available in the ocean surface layer to fuel full diatom growth:

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle }{P}_{{\mathrm{CaCO}}_{\mathrm{3}}}=P_{\mathrm{POC}}\cdot R\cdot a\cdot \left(\mathrm{1}-{\displaystyle \frac{\frac{P_{\mathrm{opal}}}{P_{\mathrm{POC}}}}{S_{\mathrm{opal}}}}\right)\mathrm{for}{\displaystyle \frac{P_{\mathrm{opal}}}{P_{\mathrm{POC}}}}<S_{\mathrm{opal}};{\displaystyle }& {\displaystyle }{P}_{{\mathrm{CaCO}}_{\mathrm{3}}}=\mathrm{0}\mathrm{for}{\displaystyle \frac{P_{\mathrm{opal}}}{P_{\mathrm{POC}}}}\ge S_{\mathrm{opal}}.\end{array}$$

Parameter R is the maximum possible rain ratio C(CaCO₃) : C(POC), a is the CaCO₃ saturation-dependent factor to account for an ocean acidification impact (see Fig. 2, following Ilyina et al., 2009), and S_opal is the threshold value of P_opal ∕ P_POC for gradual onset of CaCO₃ production. Particle fluxes and particle degradation are simulated through balance equations for sinking particulate matter as in Heinze et al. (2009, 2016), taking the saturation state for CaCO₃ and biogenic silica into account.

Figure 2 Assumed sensitivities of CaCO₃ export production in response to changes in calcite saturation (Ω). The rain ratio coefficient is used as factor “a” applied to the right-hand side of the equation for CaCO₃ production (see Eq. 1). For the control run, a was fixed at each grid point according to the prevailing pre-industrial carbonate saturation (following the moderate sensitivity). Blue: moderate sensitivity. Green: linear sensitivity. Red: extreme sensitivity. The shape of the curves would look similar at each grid point; only the pre-industrial crossover point would be different.

Scavenging of ²³⁰Th is simulated according the reversible first-order scavenging reaction (Gehlen et al., 2003; Heinze et al., 2006):

$${\displaystyle \frac{\mathrm{d}{c}_{\mathrm{part}}}{\mathrm{d}t}}=K\cdot \left[c_{\mathrm{part}}^{\mathrm{EQ}}-c_{\mathrm{part}}\right];c_{\mathrm{part}}^{\mathrm{EQ}}=k_{\mathrm{d}}\cdot c_{\mathrm{diss}}\cdot M.$$

c_part is the concentration of particle-bound ²³⁰Th, c_diss the concentration of dissolved ²³⁰Th, and $c_{\mathrm{part}}^{\mathrm{EQ}}$ the equilibrium concentration of particle-bound ²³⁰Th. The first-order rate constant K (in yr⁻¹) is set here to 10⁴ yr⁻¹, thus assuming a quasi-instantaneous equilibration. Concentration of suspended particulate material is represented by M. For the partitioning coefficient k_d between the dissolved and particle-attached phases of ²³⁰Th, we follow the formulation of Honeyman et al. (1988), which accounts implicitly for the changing reactive surfaces of particles with particle size

$${\displaystyle }{\displaystyle }{\log }_{\mathrm{10}}{k}_{\mathrm{d}}^{\mathrm{first}\mathrm{guess}}=A+B\cdot {\log }_{\mathrm{10}}M,$$

where M is the particle concentration (here taken in mg particles per litre) and A as well as B are tuneable parameters. In addition to ²³⁰Th, we carry also ²³¹Pa as well as ¹⁰Be in our model (see Heinze et al., 2006) but focus here on ²³⁰Th only. For each of these radionuclides, the preferential carrier phase can be selected in a dedicated switchboard. For ²³⁰Th scavenging, we used here CaCO₃, POC, and clay as carrier phases. With introducing a particle-specific scavenging following Hayes et al. (2015b), the final formulation for the scavenging equilibrium coefficient then becomes

$${\displaystyle }{\displaystyle }{k}_{\mathrm{d}}=\left(C_{{\mathrm{CaCO}}_{\mathrm{3}}}\cdot F_{{\mathrm{CaCO}}_{\mathrm{3}}}+C_{\mathrm{POC}}\cdot F_{\mathrm{POC}}+C_{\mathrm{clay}}\cdot F_{\mathrm{clay}}\right)\cdot k_{\mathrm{d}}^{\mathrm{first}\mathrm{guess}},$$

where C and F are the weighting coefficients and fractions of total local particle mass. The A, B, and C values are included in Table 1. Atlantic and Pacific cross sections with the k_d values for the control simulation without CO₂ emissions are given in Fig. S1 (in the Supplement).

The clay flux is computed according to the modern dust deposition from Mahowald et al. (1999) assuming that respective clay particles are chemically quasi-inert in seawater. Particle-bound ²³⁰Th phases never get to zero, even if all biogenic particles are degraded, because at each grid point there is – at least a tiny – dust flux consisting of inert clay. Therefore, particle concentrations as such never go to zero (a minimum concentration could be prescribed, but this was not necessary in our case). Because the Atlantic dust deposition may be overestimated (Gehlen et al., 2003), we assumed for the ²³⁰Th scavenging only 25 % of the respective clay flux strength in order to avoid too strong scavenging in Atlantic surface waters (see also our sensitivity experiment assuming 100 % in comparison below). We allow biological production under sea ice scaled with the local sea ice thickness.

4 Control run, scenario experiments, and sensitivity experiments

The model was spun up restarting from an earlier integration and computed into quasi-equilibrium including the sediment over 40 000 years. The equilibrium coefficient k_d was determined through a visual fit to observations of dissolved ²³⁰Th taken from the literature (Bacon and Anderson, 1982; Bacon et al., 1989; Chase et al., 2002; Cochran et al., 1987, 1995; Colley et al., 1995; Guo et al., 1995; Huh and Beasley, 1987; Moore, 1981; Moran et al., 1995, 1997; Nozaki and Horibe, 1983; Nozaki and Yang, 1987; Nozaki et al., 1987; Roy-Barman et al., 1996; Scholten et al., 1995; Vanderloeff and Berger, 1993; Vogler et al., 1998) and combined with the data as given in the GEOTRACES intermediate data product version 3 (Mawji et al., 2015; Hayes et al., 2015a, b; Deng et al., 2014). Important global bulk numbers and the parameters for the partitioning coefficient k_d are listed in Table 1.

Table 1 Summary of control run results (run P0; see Table 2).

Figure 3 Meridional ²³⁰Th cross sections in dpm/1000 L (dpm indicates disintegration per minute) for the model control run under pre-industrial atmospheric CO₂. (a) Dissolved, Atlantic. (b) Dissolved, Pacific. (c) Particle attached, Atlantic. (d) Particle attached, Pacific.

Meridional sections of dissolved ²³⁰Th concentrations for the Atlantic and Pacific oceans are given in Fig. 3. The vertical distribution clearly shows the increase of concentrations with depth due to the downward transfer of ²³⁰Th with the marine particle flux. A comparison of simulated and observed dissolved ²³⁰Th values in the water column along a meridional Atlantic cross section (Fig. 4) indicates that the model dissolved ²³⁰Th is in fairly good agreement with observations, except for a strong overestimation in some bottom water locations. This could be due to deficiencies in the flow field and corresponding problems or due to the lack of explicit simulations of hydrothermal vents and sediment resuspension. We carried out a respective sensitivity experiment (see discussion further below). A scatter plot of modelled and observed ²³⁰Th values and the locations of stations with observed data are shown in Fig. 5. The standard-run export production rates for POC and CaCO₃ are given as maps in Fig. 6.

Figure 4 Meridional dissolved ²³⁰Th Atlantic Ocean cross sections in dpm/1000 L (dpm indicates disintegration per minute): (a) interpolation from observations (data from the GEOTRACES intermediate data product; for references, please see the text). (b) Difference between the model values from the control run without CO₂ emissions and the observations.

Figure 5 Comparison of dissolved ²³⁰Th model results with observations (for references of observations, see the text). (a) Correlation between observations and modelled values. (b) Location of all vertical profiles used in this work.

Figure 6 Model control run results for (a) biological export production of particulate organic carbon and (b) CaCO₃ export.

Table 2 Model runs' overview.

An overview concerning all model simulations is given in Table 2. We carried out model projections under four different scenarios for changes in calcification (including control runs with constant calcification but rising atmospheric CO₂) and four different scenarios of the future development of the atmospheric CO₂ concentration (C8.5–2.6, L8.5–2.6, M8.5–2.6, and E8.5–2.6; see Table 2). All scenarios were restarted from the same previously performed spin-up model run reflecting a pre-industrial biogeochemical state of the ocean. Throughout the experiments, the ocean circulation field was not changed. The model was integrated during the calendar years 1200–2300 spanning a period of 1100 years. In the model spin-up simulation, the atmospheric CO₂ concentration was a prognostic variable. For the computations from 1700 onwards, we prescribed the atmospheric CO₂ concentration according to the Representative Concentration Pathways (RCPs) including their extension to the year 2300 (van Vuuren et al., 2011) used in CMIP5 (Coupled Model Intercomparison Project phase 5). Our simulations follow scenarios RCP2.6, RCP4.5, RCP6.0, and RCP8.5 designed to span a range of radiative forcing between 2.6 and 8.5 W m⁻² by the year 2100. The respective atmospheric CO₂ concentrations are shown in Fig. S2. For the decrease in calcification with decreasing pH as well as carbonate saturation, functional relationships between CaCO₃ export production and carbonate saturation were chosen which correspond to the respective experiments in Ilyina et al. (2009). For each RCP, we carried out a control simulation with constant calcification, a moderate decrease of CaCO₃ production with decreasing saturation, a linear dependency simulation with respectively stronger decrease in CaCO₃ production), and an extreme scenario (Fig. 2). To date, no clear bulk formulation for the dependency of CaCO₃ export on CaCO₃ saturation exists. Therefore, the scenarios carried out here are only sensitivity experiments and not overall exhaustive simulations in order to reproduce the entire range of possible changes in ocean biogeochemistry.

As an indication for the robustness of our results, a suite of sensitivity experiments was carried out. In a first set of sensitivity experiments (P1–P5, Table 2), we investigated how much the pre-industrial control simulation depends on selected parameter choices. In P1, parameter A for $k_{\mathrm{d}}^{\mathrm{first}\mathrm{guess}}$ was slightly increased from 5.8 to 6.0 in order to test the dependency of the result from small changes in the still not very precisely known k_d value for ²³⁰Th. Given the effect of release of material from hydrothermal vents and resuspension of material from the sediments, we carried out experiment P2 with a strong increase in scavenging in the lowermost wet model grid cells (directly over the ocean floor). To this end, parameter A for $k_{\mathrm{d}}^{\mathrm{first}\mathrm{guess}}$ was increased by 0.5 units in these grid boxes. In our general reference run (P0), we applied scavenging to clay, as if the clay flux would only be 25 % of its prescribed value. Therefore, we also added the P3 run, where 100 % clay flux was assumed also for the ²³⁰Th scavenging. Finally, in experiments P4 and P5, we switched off the scavenging to POC and CaCO₃, respectively, in order to see the importance of these two biogenic particle species on the ²³⁰Th distribution.

In a second set of sensitivity experiments, we tested how sensitive the future projection of dissolved ²³⁰Th is in view of factors other than changes in CaCO₃ production (S1–S5; see Table 2). Because we use an annually averaged model without seasonal cycle and also employed fixed annual mean velocity field, we have to use approximations in order to see the effect of a change in the velocity field or other changes involving interannual/decadal variability. We use here the time series of the Atlantic meridional overturning from the fully fledged Earth system model NorESM for a historic ramp-up and subsequent RCP8.5 forcing (Tjiputra et al., 2016; see Fig. 2 therein). In experiment S1, we scaled the three-dimensional velocity field (and also the convective mixing representation) of our simplified model with this time series, leading to an overall reduction in circulation strength. The formulation for this scaling is

$${\displaystyle }{\displaystyle }{Y}^{\mathrm{scaled}}=\left(\mathrm{1}+\mathrm{\Delta }\mathrm{Overturning}\left(t\right)\right)\cdot Y^{\mathrm{reference}},$$

with Y being a scalar variable (such as the velocity components in the x, y, and, z directions, the convective adjustment, or other specific biogeochemical parameters), Y^scaled the value as updated for the sensitivity experiments, Y^reference the original value used in the reference run, and ΔOverturning the relative change in Atlantic meridional overturning between a specific year after the year 1850 and the year 1850 itself (ΔOverturning would be zero until 1850). Using this scaling, we also investigated separately from a circulation change the consequences of variations in the maximum nutrient uptake velocity V_max in our simple trophic model (a reduction in V_max could be seen as a deterioration of growth conditions for plankton under climate change) (run S2), and in the sinking velocity of particles (corresponding to, e.g. a loss in ballast material) (run S3). In experiments S4 and S5, finally, the simultaneous effects of changes in the circulation, V_max, and the particle sinking velocity were explored for the moderate (S4) and the extreme (S5) calcification scenarios. The list of potential factors with an influence on ²³⁰Th is much longer; however, our choice of factors covers at least the physical forcing through the velocity field, surface biogeochemistry through V_max, and three-dimensional biogeochemistry through the particle sinking speed.

Figure 7 Meridional dissolved ²³⁰Th Atlantic Ocean cross sections in dpm/1000 L (dpm indicates disintegration per minute) for sensitivity experiments. (a) K_d increased everywhere (by log (K_d) = +0.2). (b) Same as panel (a) but the difference between model and observations. (c) K_d increased only in the model grid cells directly over the ocean floor (see text). (d) Same as panel (c) but the difference between model and observations. (e) Clay flux increased to 100 % of control run value. (f) Same as panel (e) but the difference between model and observations. (g) No scavenging to POC. (h) Same as panel (g) but the difference between model and observations. (i) No scavenging to CaCO₃. (j) Same as panel (i) but the difference between model and observations. Please note changes in colour codes relative to panels (a)–(h). (Observed data are from the GEOTRACES intermediate data product; for references, please see the text.)

Figure 8 Temporal evolution of Atlantic (left column panels) and Pacific (right column panels) annual mean CaCO₃ export production under the different scenarios for the sensitivity of calcification under high CO₂ (unit: Gt C yr⁻¹). From top to bottom for greenhouse gas scenarios RCP8.5, RCP6.0, RCP4.5, and RCP2.6.

Figure 9 Time series for the evolution of mean Atlantic (left column panels) and Pacific (right column panels) dissolved ²³⁰Th concentrations at 2000 m under the different scenarios for reduction of calcification under high CO₂ (unit: dpm/1000 L). From top to bottom for greenhouse gas scenarios RCP8.5. RCP 6.0, RCP4.5, and RCP2.6. The orange line indicates the theoretical detection limit for changes with respect to the pre-industrial limit. (Respective diagrams for depth levels 700 and 4000 m are given in Figs. S3 and S4.)

5 Results and discussion

The results from the sensitivity experiments relative to the pre-industrial control run are summarised in Fig. 7. Increasing the scavenging equilibrium coefficient in general (run P1, Fig. 7a and b) or over the ocean bottom (run P2, Fig. 7c and d) corrects to a substantial degree the too high bottom water vales for dissolved ²³⁰Th in the southern Atlantic but leads to high relative errors in the upper few hundred metres (where ²³⁰Th concentrations are generally low). Too low upper ocean ²³⁰Th concentrations occur also for increasing the scavenging to 100 % of the clay flux (run P3, Fig. 7e and f). Omitting, the scavenging to POC has only minor influence on the ²³⁰Th distribution because of the relatively shallow recycling of POC (run P4, Fig. 7g and h). In contrast, cancelling the scavenging to CaCO₃ (and only retaining that to POC and clay) leads to a completely unrealistic dissolved ²³⁰Th distribution with at least a 5-fold overestimation (run P5, Fig. 7i and j; note change of colour code in Fig. 7i and j). These results demonstrate that the choice of critical parameters for our pre-industrial reference simulation (run P0) results in a reasonable dissolved ²³⁰Th distribution. The crucial role of CaCO₃ particles for determining water column ²³⁰Th has been confirmed.

We discuss now the results for the projections. According to the various scenarios applied, considerable changes in surface alkalinity (not shown) and corresponding changes in CaCO₃ export production evolve at high CO₂ (Fig. 8). The effect is on average larger in the Pacific than in the Atlantic due to the overall higher modelled biological production rates of CaCO₃ in the Pacific Ocean. The extreme calcification scenario would lead to vastly reduced CaCO₃ export production after the year 2100 in all RCPs except the most moderate (RCP2.6), where even a partial recovery occurs.

For dissolved ²³⁰Th, we show here time series for depth levels 700, 2000, and 4000 m (Figs. S3, 9, and S4, respectively). According to recent intercalibration experiments, sizable discrepancies still exist between the absolute values in ²³⁰Th measurements from different laboratories, though these measurements show a smaller scatter for replicates within a single laboratory. Intercalibration experiments reveal a standard deviation for measurements on one sample by different laboratories of about 0.07–0.08 dpm/1000 L (Anderson et al., 2012). We take 0.075 dpm/1000 L as an approximate indicative value for the detection level for changes in dissolved ²³⁰Th concentrations in seawater (Figs. S3, S4, and 9; see orange lines therein). This may be a somewhat optimistic estimate for samples of past decades, but future developments of measuring techniques could possibly reduce measurement errors and the spread across analyses from different laboratories and application of different measurement methods. The detection level is shown here relative to the modelled pre-industrial values. The time of emergence as indicated by the orange line in Figs. S3, S4, and 9 is therefore the earliest possible time of detection if pre-industrial ²³⁰Th would be known. Figures S3, S4, and 9 show the increase of dissolved ²³⁰Th activities with depth from 700 m (Fig. S3), over 2000 m (Fig. 9), to finally 4000 m (Fig. S4). At 700 m, the increase in ²³⁰Th activities due to the assumed ocean acidification effect is too small to unambiguously show an effect for RCP2.6 (Fig. S3). This picture changes when going down to 2000 m and below, where an effect would be detectable within this century at least for the stronger forcing scenarios (RCP8.5 and RCP6.0) (Figs. 9 and S4). For constant CaCO₃ production, the intermediate and deep-water ²³⁰Th activities start to rise around the year 2100 as well (see black curves in Figs. S3, S4, and 9). This effect is due to the increasing dissolution of CaCO₃ particles in the water column in parallel with downward mixing of waters that carry anthropogenic loads of dissolved organic carbon and hence subsurface and deep acidification. The effect is most important in areas where anthropogenic carbon is mixed down quickly and induces a significant shoaling of the CaCO₃ saturation level and CaCO₃ lysocline as well the carbonate compensation depth through deep-water acidification. Parts of the deep Pacific are not as much influenced by this as compared to the Atlantic within the 21st century. The control run in the figures does not represent the reference run with constant pre-industrial atmospheric pCO₂ but the run with constant CaCO₃ production and rising pCO₂ according to the RCP scenarios.

In order to look for most suitable regions for detection of ocean acidification effects on CaCO₃ particle fluxes using modern ²³⁰Th data, we plotted the year of emergence for different reference years and different assumed detection levels (Fig. 10). The year of emergence is defined here as the earliest possible year of a potential detection in the case that the ²³⁰Th signal (rate of change over time) is only influenced by a decrease in calcification and not any other processes (such as potential shifts in circulation and associated changes in biological particle production). We choose 2010 as reference year, to which the ²³⁰Th activity at a later stage would be compared to, and 0.075 dpm/1000 L as the limit for the detection of a signal. For the moderate calcification scenario, the time of emergence at least at 4000 m would be comparable to the time of emergence as potentially inferable through repeat surface alkalinity measurements (Ilyina et al., 2009) (Fig. 10a and b). For the extreme scenario, the signal would be identifiable considerably earlier. The signal emerges earliest in deep Pacific waters.

The results on the early warning indicator as shown in Fig. 10 are based on the simplifying assumption that only changes in calcification caused by ocean acidification influence the water column ²³⁰Th concentrations under progressing acidification and climate forcing. We therefore also mimicked the effect of interannual/decadal variability in circulation and biogeochemical parameters (runs S1–S5, Table 2). Comparing the projections under RCP8.5 and the moderate calcification scenario for the case without interannual/decadal variability (run M8.5, Table 2) with the one including variability (runs S1–S3, Table 2) for the middle of this century reveals differences in the order of the assumed detection level for changes in the velocity field and the particle sinking velocity and smaller changes for varying V_max (Fig. 11a–f). The effect on modulating water column ²³⁰Th can be quite variable for different oceanic domains. Deviations can add up or compensate for each other locally. For simultaneous changes of circulation, V_max, and particle sinking velocity, the overall deviations from the run without variability amount to 2–3 times the theoretical detection limit (Fig. 11g and h). We repeated the analysis of the time of emergence with the runs including interannual variability for the moderate (S4) and also the extreme (S5) scenario. For this analysis, we again chose 2010 as the reference year but increased the detection limit to 2.5 × 0.075 dpm/1000 L. Under these more realistic conditions, the time of emergence would be shifted to the end of the century for the moderate calcification scenario and to the middle of the century for the extreme scenario (Fig. 12).

Figure 10 Prospective year of emergence for an ocean-acidification-induced signal in dissolved ²³⁰Th activity as derived from the model runs. Shown are the calendar years of emergence for the depth levels at 2000 and 4000 m. All figures are shown for the strong RCP8.5 scenario concerning atmospheric CO₂ concentration and the moderate scenario of calcification decrease with saturation state, with the year 2010 as a reference year for ²³⁰Th activity and 0.075 dpm/1000 L as analytical threshold between different samples. (a) For the moderate calcification scenario, at 2000 m depth. (b) Same as panel (a) but for depth level 4000 m. (c) For the extreme calcification scenario, at 2000 m depth. (d) Same as panel (c) but for depth level 4000 m.

Figure 11 Meridional cross sections for sensitivity experiments concerning the projections including a representation of interannual/decadal variability. Shown are differences in dissolved ²³⁰Th for forcing relative to the normal run with RCP8.5 and the moderate calcification scenario for the year 2050. (a) Change in velocity field (for details, see the text). (b) As (a) but for the Pacific. (c) Change of V_max for phosphate uptake (for details, see the text), Atlantic. (d) As (c) but for the Pacific. (e) Change of the sinking velocity of particulate matter (for details, see the text). (f) As (e) but for the Pacific. (g) Simultaneous change of velocity field, V_max, and sinking velocity. (h) As (g) but for the Pacific.

Figure 12 Prospective year of emergence for an ocean-acidification-induced signal in dissolved ²³⁰Th activity as derived from the model runs. Shown are the calendar years of emergence for the depth levels at 2000 and 4000 m. All figures are shown for the strong RCP8.5 scenario concerning atmospheric CO₂ concentration and the moderate scenario of calcification decrease with saturation state, with the year 2010 as a reference year for ²³⁰Th activity and 2.5 ⋅ 0.075 dpm/1000 L as analytical threshold between different samples. In contrast to Fig. 11, here, an analysis including interannual/decadal variability in selected biogeochemical parameters is shown. (a) For the moderate calcification scenario, at 2000 m depth. (b) Same as panel (a) but for depth level 4000 m. (c) For the extreme calcification scenario, at 2000 m depth. (d) Same as panel (c) but for depth level 4000 m.

The question arises whether surface alkalinity and deep ²³⁰Th measurements possibly could be combined successfully in order to detect ocean acidification impacts on calcification. We could, so far, not detect any simple relationship between surface changes in CaCO₃ export production, surface alkalinity, and ²³⁰Th in spite of the conceptually clear interdependencies. In order to illustrate this, we plotted changes in surface alkalinity and 3000 m ²³⁰Th over time for the run with the moderate calcification scenario and RCP8.5 forcing between 2040 and 2010.

The distributions show a somewhat consistent picture without interannual/decadal variability (run M8.5, Fig. 13a and b) with general rises in both surface alkalinity and deep ²³⁰Th. For the case including interannual/decadal variability (run S4, Fig. 13c and d), the temporal gradients in both variables become smaller and the deep ²³⁰Th even may reverse. This does not mean that it would be impossible to construct a combined Alk-²³⁰Th tracer – only that such a tracer cannot be derived from this study.

As compared to the detection approach for ocean acidification impacts through total alkalinity measurements as pursued by Ilyina et al. (2009), the ²³⁰Th method presented here results in similar time of emergence of a signal reflecting changes in biocalcification if one omits interannual/decadal variability and other possible influence factors. Probably the most important limitation of the approach here is the fixed velocity field which does not vary realistically with seasons, interannual variability, and climatic variability modes (such as the El Niño–Southern Oscillation or North Atlantic Oscillation). Also the effect of boundary scavenging (Anderson et al., 1983; Roy-Barman et al., 2009), i.e. the transport of dissolved ²³⁰Th from areas of low particle concentrations to those with high particle concentrations (especially at ocean margins and shelf seas), is not spatially resolved in our model. Further uncertainties are associated with the choice of the k_d values and particle-specific reactivity.

The extreme scenario on pH-dependent CaCO₃ production is likely to be an overestimation. According to the results from this study and Ilyina et al. (2009), a respective large change in the real world would probably have already been detected through the alkalinity signal and other methods such as remote sensing and sediment trap measurements. In addition, the reaction of CaCO₃ shell producing organisms to sinking carbonate saturation levels is not a simple function and can vary between taxa (Kroeker et al., 2013).

Figure 13 Surface alkalinity changes and deep ²³⁰Th changes during 2040–2010 as for the RCP8.5 projection using moderate calcification changes. Left column without variability, (a) at the sea surface, (b) at 3000 m. Right column with decadal parameter variability, (c) at the sea surface, (b) at 3000 m.

The detection method for ocean acidification impacts through ²³⁰Th would work best if the pre-industrial values for dissolved oceanic ²³⁰Th activities could reliably be reconstructed. From core top samples and corals, one could possibly determine whether the particle-attached ²³⁰Th activities and hence those of the dissolved fraction would have undergone any variations over the last centuries and millennia. However, the overall uncertainties associated with the formation and analysis of the paleo-record values may be too high to provide an accurate enough baseline for comparisons with modern water column measurements.

6 Conclusions

In this study, we investigated the potential of the particle reactive radionuclide ²³⁰Th for detection of reduced calcification by biota due to progressing ocean acidification. The time of emergence of a dissolved ²³⁰Th activity signal with respect to a reference year and reference activity is about the same as for detection approaches employing alkalinity if one disregards effects of interannual/decadal variability in the flow field and influence factors other than ocean-acidification-induced changes in CaCO₃ export production. Taking into account interannual/decadal variability including ocean circulation changes may delay the signal emergence by several decades. Nevertheless, regular selected reoccupations of a series of deep stations in the Pacific Ocean and Southern Ocean with highest-quality ²³⁰Th measurements would be helpful to accompany alkalinity measurements that are easier to make in order to see whether both tracers give consistent results. Surface alkalinity measurements include signals of natural variability on seasonal and multiyear scales. Likewise, changes in ocean circulation and changes in biological particle production due to processes other than ocean acidification may lead to changes in the marine ²³⁰Th distribution (and alkalinity). As the deep ocean is less prone to effects of natural variability and the quality of observations does not change with depth, deep-ocean observations of ²³⁰Th could be advantageous for monitoring and detecting ocean acidification effects on calcification. In any case, ²³⁰Th represents a fascinating magnifying glass for changes in ocean surface processes seen through the deep-ocean signal. Its potential has not yet been fully exploited.