Incorporating the stable carbon isotope ¹³C in the ocean biogeochemical component of the Max Planck Institute Earth System Model

The stable carbon isotopic composition (δ¹³C) is an important variable to study the ocean carbon cycle across different timescales. We include a new representation of the stable carbon isotope ¹³C into the HAMburg Ocean Carbon Cycle model (HAMOCC), the ocean biogeochemical component of the Max Planck Institute Earth System Model (MPI-ESM). ¹³C is explicitly resolved for all oceanic carbon pools considered. We account for fractionation during air–sea gas exchange and for biological fractionation ϵ_p associated with photosynthetic carbon fixation during phytoplankton growth. We examine two ϵ_p parameterisations of different complexity: ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ varies with surface dissolved CO₂ concentration (Popp et al., 1989), while ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ additionally depends on local phytoplankton growth rates (Laws et al., 1995). When compared to observations of δ¹³C of dissolved inorganic carbon (DIC), both parameterisations yield similar performance. However, with regard to δ¹³C in particulate organic carbon (POC) ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ shows a considerably improved performance compared to ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$. This is because ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ produces too strong a preference for ¹²C, resulting in δ¹³C_POC that is too low in our model. The model also well reproduces the global oceanic anthropogenic CO₂ sink and the oceanic ¹³C Suess effect, i.e. the intrusion and distribution of the isotopically light anthropogenic CO₂ in the ocean.

The satisfactory model performance of the present-day oceanic δ¹³C distribution using ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ and of the anthropogenic CO₂ uptake allows us to further investigate the potential sources of uncertainty of the Eide et al. (2017a) approach for estimating the oceanic ¹³C Suess effect. Eide et al. (2017a) derived the first global oceanic ¹³C Suess effect estimate based on observations. They have noted a potential underestimation, but their approach does not provide any insight about the cause. By applying the Eide et al. (2017a) approach to the model data we are able to investigate in detail potential sources of underestimation of the ¹³C Suess effect. Based on our model we find underestimations of the ¹³C Suess effect at 200 m by 0.24 ‰ in the Indian Ocean, 0.21 ‰ in the North Pacific, 0.26 ‰ in the South Pacific, 0.1 ‰ in the North Atlantic and 0.14 ‰ in the South Atlantic. We attribute the major sources of underestimation to two assumptions in the Eide et al. (2017a) approach: the spatially uniform preformed component of δ¹³C_DIC in year 1940 and the neglect of processes that are not directly linked to the oceanic uptake and transport of chlorofluorocarbon-12 (CFC-12) such as the decrease in δ¹³C_POC over the industrial period.

The new ¹³C module in the ocean biogeochemical component of MPI-ESM shows satisfying performance. It is a useful tool to study the ocean carbon sink under the anthropogenic influences, and it will be applied to investigating variations of ocean carbon cycle in the past.

1 Introduction

The stable carbon isotopic composition (δ¹³C) measured in carbonate shells of fossil foraminifera is one of the most widely used properties in paleoceanographic research (Schmittner et al., 2017). It is defined as a normalised ratio between the stable carbon isotopes ¹³C and ¹²C:

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\left(\mathrm{‰}\right)=\left({\displaystyle \frac{{}^{\mathrm{13}}\mathrm{C}/{}^{\mathrm{12}}\mathrm{C}}{R_{\mathrm{std}}}}-\mathrm{1}\right)\cdot \mathrm{1000},\end{array}$$

where R_std is an arbitrary standard ratio. In observational studies, the ratio ${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ in Pee Dee Belemnite (PDB; Craig, 1957) is conventionally used for R_std.

δ¹³C provides information on past changes in water mass distribution and properties (e.g. Curry and Oppo, 2005; Peterson et al., 2014). Direct comparison between paleo-δ¹³C measurements and simulated δ¹³C facilitates evaluating the ability of Earth system models (ESMs) to simulate paleo-ocean states. For this reason, we present a new implementation of ¹³C in the HAMburg Ocean Carbon Cycle model (HAMOCC6), the ocean biogeochemical component of the Max Planck Institute Earth System Model (MPI-ESM). A comprehensive representation of δ¹³C is a timely extension of MPI-ESM in support of planned simulations of a complete last glacial cycle within the German climate modelling initiative PalMod (Latif et al., 2016). Before applying the new ¹³C module to paleo-simulations, we evaluate it by comparison to observational data in the present-day ocean.

Earlier versions of HAMOCC already featured a ¹³C module, for instance HAMOCC2s (Heinze and Maier-Reimer, 1999) and HAMOCC3 (Maier-Reimer, 1993). HAMOCC3 included prognostic ¹³C variables for dissolved inorganic carbon (DIC), particulate organic matter and calcium carbonate. HAMOCC3 also accounted for temperature-dependent isotopic fractionation during air–sea gas exchange (higher δ¹³C of surface DIC in colder water) and biological fractionation during carbon fixation. Due to the simplified representation of marine biological production in HAMOCC3, biological fractionation was based on fixation of inorganic carbon into non-living particulate organic matter and was parameterised by a spatially and temporally uniform factor. This approach for biological fractionation of ¹³C, however, could not reproduce the observed large meridional gradient of δ¹³C in particulate organic matter (Goericke and Fry, 1994). Since then, HAMOCC3 was refined in particular with regard to its representation of plankton dynamics. The current version HAMOCC6 resolves bulk phytoplankton, zooplankton, detritus, dissolved organic carbon (Six and Maier-Reimer, 1996) and nitrogen-fixing cyanobacteria (Paulsen et al., 2017). We thus develop an updated ¹³C module that considers the refined ecosystem representation and test different non-uniform parameterisations for biological fractionation during phytoplankton growth.

To choose a suitable biological fractionation parameterisation for our model, we test the parameterisations of Popp et al. (1989) and Laws et al. (1995). These parameterisations are selected for two reasons. First, they are of different complexities. The parameterisation of Popp et al. (1989) empirically relates ¹³C biological fractionation to the concentration of dissolved CO₂ in seawater, whereas that of Laws et al. (1995) considers dissolved CO₂ concentration and phytoplankton growth rate. Second, input variables in these two parameterisations are explicitly computed in the model. We omit more complex parameterisations that include effects of cell membrane permeability of molecular CO₂ diffusion, cell size and shape (e.g. Rau et al., 1996; Keller and Morel, 1999), as HAMOCC6 does not resolve these features of plankton cells.

Oceanic δ¹³C measurements were mostly carried out in the late 20th century. In the upper ocean δ¹³C in dissolved inorganic carbon (δ¹³C_DIC) has been observed to noticeably decrease in response to the intrusion of anthropogenic CO₂ from fossil fuel combustion which carries a lower ${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ signal (Gruber et al., 1999; Quay et al., 2003). Such δ¹³C_DIC decrease is referred to as the oceanic ¹³C Suess effect (Keeling, 1979). Recently, Eide et al. (2017a) derived an observation-based estimate of the global ocean ¹³C Suess effect since pre-industrial times. Such an observation-based estimate is valuable as it is the basis of an almost independent estimate of the global ocean anthropogenic carbon uptake. And it could be used for evaluating models at pre-industrial states (Buchanan et al., 2019; Tjiputra et al., 2020) and for setting up paleo-simulations (O'Neill et al., 2019).

Yet, Eide et al. (2017a) have noted that their approach might underestimate the oceanic ¹³C Suess effect. They conjectured an underestimation of the ¹³C Suess effect between 0.15 ‰–0.24 ‰ at 200 m depth in 1994. However, the quantitative spatial distribution of this underestimation is unclear. Moreover, although Eide et al. (2017a) have related the underestimation to several assumptions in the approach they applied, the quantitative impact of these assumptions is still unclear as the measurements are too limited in space and time to perform in-depth investigation.

Our model data include all parameters needed to apply the Eide et al. (2017a) procedure, which relies on regressional relationships between preformed δ¹³C_DIC (related to the transport of surface waters with specific DIC and DI¹³C) and CFC-12 (chlorofluorocarbon-12) partial pressure. Thus, our consistent model framework, with the complete spatio-temporal information of the hydrological and biogeochemical variables, enables us to investigate the spatial distribution of the above-mentioned potential underestimation of the oceanic ¹³C Suess effect. Moreover, our model framework also allows for the attribution of the underestimation to the assumptions of the procedure Eide et al. (2017a) applied.

In the following sections, we first provide a brief introduction to the global ocean biogeochemical model HAMOCC6, followed by a description of the new ¹³C module including the experimental setup (Sect. 2). Section 3 presents the model evaluation against observations in the late 20th century, and Sect. 4 evaluates the simulated oceanic ¹³C Suess effect. Section 5 addresses our findings on testing the Eide et al. (2017a) approach for estimating the oceanic ¹³C Suess effect. Summary and conclusions are given in Sect. 6.

2 Model description

2.1 The global ocean biogeochemical model (HAMOCC6)

HAMOCC6 (Ilyina et al., 2013; Paulsen et al., 2017; Mauritsen et al., 2019) includes biogeochemical processes in the water column and in the sediment. In the water column, the following biogeochemical tracers are simulated: dissolved inorganic carbon (DIC), total alkalinity (TA), phosphate (PO₄), nitrate (NO₃), nitrous oxide (N₂O), dissolved nitrogen gas (N₂), silicate (SiO₄), dissolved bioavailable iron (Fe), dissolved oxygen (O₂), bulk phytoplankton (Phy), cyanobacteria (Cya), zooplankton (Zoo), dissolved organic matter (DOM), particulate organic matter (POM), opal shells, calcium carbonate shells (CaCO₃), terrigenous material (Dust) and hydrogen sulfide (H₂S). Below the model-defined export depth (100 m), the sinking speed of POM linearly increases with depth. Theoretically, this leads to a power-law-like attenuation of POM fluxes as observations (Martin et al., 1987; Kriest and Oschlies, 2008). Constant sinking speeds are set for opal, CaCO₃ and Dust. Except for CaCO₃ and opal, whose sinking speeds (30 and 25 m d⁻¹, respectively) are considerably faster than the horizontal velocities of ocean flow, the water-column biogeochemical tracers are transported by the hydrodynamical fields in the same manner as salinity.

The sediment module is based on Heinze et al. (1999). It simulates remineralisation and dissolution processes as in the water column concerning dissolved tracers (PO₄, NO₃, N₂, O₂, SiO₄, Fe, H₂S, DIC and TA) in the pore water and the solid sediment constituents (POM, opal, CaCO₃). The tracers in the pore water are exchanged with the overlying water column by diffusion. Pelagic sedimentation fluxes of POM, CaCO₃ and opal are added to the solid components of the sediment. Below the active sediment there is one layer containing only solid sediment components and representing burial. To balance the loss of nutrients, TA, DIC and SiO₄ in the water column, constant input fluxes of DOM, CO${}_{\mathrm{3}}^{\mathrm{2}-}$ and SiO₄ are added uniformly at the ocean surface, whose rates are derived from a linear regression of the long-term (approximately 100 years) temporal evolution of the sediment (active and burial) inventory.

A detailed description of HAMOCC6 is provided in Mauritsen et al. (2019) and the references therein. Different to the HAMOCC6 version in Mauritsen et al. (2019), we allow DOM degradation in low-oxygen conditions until all available O₂ is consumed.

2.2 The stable carbon isotope ¹³C in HAMOCC6

HAMOCC6 simulates total carbon C, which is the sum of the three natural isotopes ¹²C, ¹³C and ¹⁴C. Because in nature ¹²C constitutes about 98.9 % of the total carbon and ¹³C only constitutes about 1.1 % (Lide, 2002), in HAMOCC6 we assume ¹²C=C. We include a ¹³C counterpart for each ¹²C prognostic variable; that is, we introduce seven new tracers for the water column and three for the sediment. ¹³C only mimics the ¹²C biogeochemical fluxes, modified by the corresponding isotopic fractionation. We assume ¹³C inventory to be as large as the inventory of ¹²C to reduce numerical errors. Consequently, the reference standard of the stable carbon isotope ratio R_std is set to 1 in Eq. (1). In this section, we describe the implementation of ¹³C fractionation during air–sea exchange and carbon uptake by bulk phytoplankton and by cyanobacteria. Because the isotopic fractionation during the production of calcium carbonate is small (Turner, 1982) and uncertain (Zeebe and Wolf-Gladrow, 2001), it is not considered in this study, following the model studies of e.g. Lynch-Stieglitz et al. (1995), Schmittner et al. (2013) and Tjiputra et al. (2020).

2.2.1 Fractionation during air–sea gas exchange

The net air–sea CO₂ gas exchange flux F reads

$$\begin{array}{}\text{(2)}& F=-k_{{\mathrm{CO}}_{\mathrm{2}}}{\mathit{\gamma }}_{{\mathrm{CO}}_{\mathrm{2}}}\left(p{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{surf}}-p{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{atm}}\right).\end{array}$$

Here, $p{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{surf}}$ and $p{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{atm}}$ are the partial pressures of CO₂ in the surface seawater and in the atmosphere, respectively. The piston velocity $k_{{\mathrm{CO}}_{\mathrm{2}}}$ (m s⁻¹) for CO₂ and the solubility ${\mathit{\gamma }}_{{\mathrm{CO}}_{\mathrm{2}}}$ (mol L⁻¹ atm⁻¹) of CO₂ are calculated following Wanninkhof (2014) and Weiss (1974), respectively.

Similar to the air–sea flux of total carbon in Eq. (2), the net air–sea ¹³CO₂ exchange flux ¹³F reads

$$\begin{array}{}\text{(3)}& {}^{\mathrm{13}}F={-}^{\mathrm{13}}{k}_{{\mathrm{CO}}_{\mathrm{2}}}{}^{\mathrm{13}}{\mathit{\gamma }}_{{\mathrm{CO}}_{\mathrm{2}}}\left(p{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{surf}}{R}_g-p{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{atm}}{R}_{\mathrm{atm}}\right),\end{array}$$

in which R_g and R_atm are the ratios of ${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ in surface pCO₂ and in atmospheric CO₂, respectively. Following Zhang et al. (1995), we can rewrite Eq. (3) as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{}^{\mathrm{13}}F& =-k_{{\mathrm{CO}}_{\mathrm{2}}}{\mathit{\alpha }}_{\mathrm{k}}{\mathit{\gamma }}_{{\mathrm{CO}}_{\mathrm{2}}}{\mathit{\alpha }}_{\mathrm{aq}\leftarrow g}\\ & \left(p{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{surf}}{\displaystyle \frac{R_{\mathrm{DIC}}}{{\mathit{\alpha }}_{\mathrm{DIC}\leftarrow g}}}-p{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{atm}}{R}_{\mathrm{atm}}\right).\end{array}\end{array}$$

Here, ${\mathit{\alpha }}_{\mathrm{k}}{=}^{\mathrm{13}}{k}_{{\mathrm{CO}}_{\mathrm{2}}}/k_{{\mathrm{CO}}_{\mathrm{2}}}$ is the kinetic fractionation factor, ${\mathit{\alpha }}_{\mathrm{aq}\leftarrow g}{=}^{\mathrm{13}}{\mathit{\gamma }}_{{\mathrm{CO}}_{\mathrm{2}}}/{\mathit{\gamma }}_{{\mathrm{CO}}_{\mathrm{2}}}$ is the equilibrium isotopic fractionation factor for gas dissolution (from gaseous to aqueous CO₂), ${\mathit{\alpha }}_{\mathrm{DIC}\leftarrow g}=R_{\mathrm{DIC}}/R_g$ is the equilibrium isotopic fractionation factor from gaseous CO₂ to DIC and $R_{\mathrm{DIC}}={{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{DIC}}{/}^{\mathrm{12}}{\mathrm{C}}_{\mathrm{DIC}}$. Parameters α_k, α_aq←g and α_DIC←g are temperature-dependent, and they are obtained from laboratory experiments (Zhang et al., 1995), often expressed in terms of a per mil fractionation factor $\mathit{ϵ}\left(\mathrm{‰}\right)=(\mathit{\alpha }-\mathrm{1})\times {\mathrm{10}}^{\mathrm{3}}$:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathit{ϵ}}_k=-\mathrm{0.85},\text{(6)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{aq}\leftarrow g}=\mathrm{0.0049}{T}_{\mathrm{C}}-\mathrm{1.31},\text{(7)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{DIC}\leftarrow g}=\mathrm{0.014}{T}_{\mathrm{C}}{f}_{{\mathrm{CO}}_{\mathrm{3}}}-\mathrm{0.105}{T}_{\mathrm{C}}+\mathrm{10.53}.\end{array}$$

Here, T_C is the seawater temperature in ^∘C, and $f_{{\mathrm{CO}}_{\mathrm{3}}}={\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}/\mathrm{DIC}$ is the fraction of carbonate ions in DIC. Because in Eq. (6) the temperature dependency is weak, we use a constant ${\mathit{ϵ}}_{\mathrm{aq}\leftarrow g}=-\mathrm{1.24}$, obtained at T_C=15 ^∘C in the model, following Schmittner et al. (2013). In Eq. (7) we neglect the first term $\mathrm{0.014}{T}_{\mathrm{C}}{f}_{{\mathrm{CO}}_{\mathrm{3}}}$, because $f_{{\mathrm{CO}}_{\mathrm{3}}}$ is generally smaller than 0.1 and because the constant factor is 1 order of magnitude smaller than that of the second term 0.105T_C.

Note that Eq. (5) $({\mathit{ϵ}}_k=-\mathrm{0.85})$ and the simplified Eq. (7) $({\mathit{ϵ}}_{\mathrm{DIC}\leftarrow g}=-\mathrm{0.105}{T}_{\mathrm{C}}+\mathrm{10.53})$ in this study, adopting those of Schmittner et al. (2013), are slightly different from the OMIP protocol (Orr et al., 2017; ${\mathit{ϵ}}_k=-\mathrm{0.88}$ and ${\mathit{ϵ}}_{\mathrm{DIC}\leftarrow g}=\mathrm{0.014}{T}_{\mathrm{C}}{f}_{{\mathrm{CO}}_{\mathrm{3}}}-\mathrm{0.107}{T}_{\mathrm{C}}+\mathrm{10.53})$. Results of a short pre-industrial simulation with ϵ_k and ϵ_DIC←g from OMIP protocol yield a negligible difference (not shown). In our future simulations ϵ_k and ϵ_DIC←g suggested by the OMIP protocol will be used.

2.2.2 Fractionation during phytoplankton growth

The lighter stable carbon isotope ¹²C is preferentially utilised over ¹³C during photosynthesis (O'Leary, 1988). Following Schmittner et al. (2013), we formulate this isotopic fractionation during net growth of the bulk phytoplankton and cyanobacteria as

$$\begin{array}{}\text{(8)}& {}^{\mathrm{13}}G=R_{\mathrm{DIC}}{\mathit{\alpha }}_{\mathrm{Phy}\leftarrow \mathrm{DIC}}G,\end{array}$$

with

$$\begin{array}{}\text{(9)}& {\mathit{\alpha }}_{\mathrm{Phy}\leftarrow \mathrm{DIC}}={\mathit{\alpha }}_{\mathrm{aq}\leftarrow \mathrm{DIC}}{\mathit{\alpha }}_{\mathrm{Phy}\leftarrow \mathrm{aq}}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{aq}\leftarrow g}}{{\mathit{\alpha }}_{\mathrm{DIC}\leftarrow g}}}{\mathit{\alpha }}_{\mathrm{Phy}\leftarrow \mathrm{aq}}.\end{array}$$

Here G (µmol C L⁻¹ d⁻¹) denotes the growth of bulk phytoplankton or cyanobacteria. α_Phy←DIC is the isotopic fractionation factor for DIC fixation, which is determined by the equilibrium fractionation factor α_aq←DIC from DIC to aqueous CO₂(aq) and by the biological fractionation factor ${\mathit{ϵ}}_{\mathrm{p}}=({\mathit{\alpha }}_{\mathrm{Phy}\leftarrow \mathrm{aq}}-\mathrm{1})\times {\mathrm{10}}^{\mathrm{3}}$ related to the fixation of CO₂(aq). Here the subscript “Phy” denotes either the bulk phytoplankton or cyanobacteria.

We test the parameterisations for biological fractionation from Popp et al. (1989) and from Laws et al. (1995), i.e.

$$\begin{array}{}\text{(10)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}=-\mathrm{17}\log \left({\mathrm{CO}}_{\mathrm{2}}\right(\mathrm{aq}\left)\right)+\mathrm{3.4},\text{(11)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}=\left({\displaystyle \frac{\mathit{\mu }}{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)/{\mathit{\rho }}_{\mathrm{sea}}}}-\mathrm{0.371}\right)/\mathrm{0.015}.\end{array}$$

Here, CO₂(aq) (µmol L⁻¹) is aqueous CO₂ in surface water, and variable μ (d⁻¹) is the specific growth rate of bulk phytoplankton or of cyanobacteria. Note that Laws et al. (1995) measured ϵ_aq←Phy. Because α_Phy←aq is close to unity, ${\mathit{ϵ}}_{\mathrm{p}}\approx -{\mathit{ϵ}}_{\mathrm{aq}\leftarrow \mathrm{Phy}}$ (Zeebe and Wolf-Gladrow, 2001). In Eq. (11), we set the seawater density ρ_sea a constant value of 1.025 kg L⁻¹. Then Eq. (11) is simplified to

$$\begin{array}{}\text{(12)}& {\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}=\mathrm{68.3}{\displaystyle \frac{\mathit{\mu }}{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)}}-\mathrm{24.7}.\end{array}$$

Both CO₂(aq) and μ (depending on local conditions of light, water temperature and nutrient availability) are determined in HAMOCC6. Figure 1 illustrates the values of ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ and ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ under typical ranges of CO₂(aq) and μ in the ocean. When μ≤1, ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ is generally more negative than ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$. For high μ values, e.g. μ=2, ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ is constantly less negative than ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$. Under high μ and low CO₂(aq), ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ becomes positive, which is unrealistic. However, our simulated ratios of phytoplankton growth rate to dissolved CO₂ concentration do not produce unrealistic positive ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ at any time step in this study.

Figure 1 The per mil biological fractionation factor ϵ_p against aqueous CO₂ concentration. The solid line illustrates ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$, in which the biological fractionation during phytoplankton growth is only a function of CO₂(aq). The dashed–dotted lines show ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$, which depends on $\mathit{\mu }/{\mathrm{CO}}_{\mathrm{2}}$, the ratio of phytoplankton growth rate to CO₂(aq), for μ=0.2 (blue), 0.6 (red), 1.2 (yellow) and 2.0 (purple) d⁻¹.

2.3 Model setup and experimental design

2.3.1 Setup

We conduct ocean-only simulations using the MPIOM-1.6.3p1 (Jungclaus et al., 2013; Notz et al., 2013; Mauritsen et al., 2019) with HAMOCC6. MPIOM is a free-surface ocean general circulation model. It uses a curvilinear grid with the grid poles located over Greenland and Antarctica. We use a low-resolution configuration with a nominal horizontal resolution of 1.5^∘. This configuration has a minimum grid spacing of 15 km around Greenland and a maximum grid spacing of 185 km in the tropical Pacific. There are 40 unevenly spaced vertical levels. The layer thickness increases from 10 m in the upper ocean to 600 m in the deep ocean. The upper 100 m of the water column is represented by nine levels. The time step is 1 h. In this setup, we additionally include the oceanic uptake and transport of CFC-12. CFC-12 is chemically inert and can therefore be treated as a conservative and passive tracer participating in all hydrodynamical processes within the ocean identical to e.g. salinity. The implementation of the air–sea gas exchange of CFC-12 follows the OMIP protocol (Orr et al., 2017).

2.3.2 Experimental design

For the pre-industrial spin-up simulations we cyclically apply the 1905–1929 sea-surface boundary conditions from ERA20C (Poli et al., 2016, covering 1901–2010). The atmospheric CO₂ mixing ratio is set to 280 ppmv. A spin-up run is first conducted without ¹³C tracers until the long-term averaged global net air–sea CO₂ flux is smaller than 0.05 Pg C yr⁻¹ (adequate to the C4MIP criterion for steady-state conditions of <0.1 Pg C yr⁻¹; Jones et al., 2016). This model state is the starting point for the two spin-up runs including ¹³C tracers, PI_Popp and PI_Laws, which are based on the biological fractionation parameterisation ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ (Eq. 10) and ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ (Eq. 12), respectively.

The ¹³C tracers are initialised as follows. The mean δ¹³C of the marine organic matter is about −20 ‰ (Degens et al., 1968). Therefore, we set the initial concentrations of ¹³C in the bulk phytoplankton, cyanobacteria, zooplankton, dissolved organic carbon, particulate organic carbon in the water column and particulate organic carbon in the sediment to 0.98 (according to Eq. 1) of their ¹²C counterparts. The initial ¹³C_DIC in the water column is calculated using the relation between δ¹³C_DIC and PO₄ (Lynch-Stieglitz et al., 1995),

$$\begin{array}{}\text{(13)}& {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{DIC}}=\mathrm{2.7}-\mathrm{1.1}{\mathrm{PO}}_{\mathrm{4}},\end{array}$$

and Eq. (1). Here PO₄ and DIC are from the quasi-equilibrium state of the spin-up run without ¹³C tracers. The initial concentrations of ¹³C${}_{{\mathrm{CaCO}}_{\mathrm{3}}}$ in the water column and in the sediment and the initial concentration of ¹³C_DIC in pore water are set identical to their ¹²C counterparts.

The pre-industrial stable carbon isotope ratio δ¹³CO₂ of atmospheric CO₂ is fixed at −6.5 ‰. The inputs of dissolved organic ¹³C (DO¹³C) and ¹³CO${}_{\mathrm{3}}^{\mathrm{2}-}$ are uniformly added at the ocean surface. The input rate of DO¹³C is calculated as the product of the input rate of DOC and the sea-surface ${\mathrm{DO}}^{\mathrm{13}}\mathrm{C}/\mathrm{DOC}$ ratio; the input rate of ¹³CO${}_{\mathrm{3}}^{\mathrm{2}-}$ is the product of the input rate of CO${}_{\mathrm{3}}^{\mathrm{2}-}$ and the sea-surface ${}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}/{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ ratio. This approach to determine ¹³C input rates results in a small drift in the water-column ¹³C inventory, but it only has minor impact on the simulation results (see Appendix A).

PI_Popp and PI_Laws are spun up for 2500 simulation years. Equilibrium states are reached with 98 % of the ocean volume having a δ¹³C_DIC drift of less than 0.001 ‰ yr⁻¹ (employing the same criteria as for ¹⁴C in OMIP protocol, Orr et al., 2017). An equilibrium of the sediment is, however, not achieved for either ¹³C or other biogeochemical tracers.

In the transient simulations for the historical period 1850–2010, Hist_Popp and Hist_Laws, we prescribe increasing atmospheric CO₂ mixing ratios (Meinshausen et al., 2017) due to anthropogenic activities and decreasing atmospheric δ¹³CO₂ following OMIP and C4MIP protocols (Jones et al., 2016) (Fig. 2a). For the period 1850–1900, when forcing data are absent, we continue applying the 1905–1929 ERA20C cyclic forcing. From 1901 to 2010, we use the transient ERA20C forcing. The evolution of the atmospheric CFC-12 concentration (Fig. 2b) follows Bullister (2017). Because the atmospheric CFC-12 is slightly higher in the Northern Hemisphere, we prescribe a linear transition between 10^∘ S and 10^∘ N. Input rates of DO¹³C, DOC, ¹³CO${}_{\mathrm{3}}^{\mathrm{2}-}$, CO${}_{\mathrm{3}}^{\mathrm{2}-}$ and SiO₄ are kept constant and are the same as those of pre-industrial simulations.

Figure 2 (a) The evolution of atmospheric CO₂ (blue, Meinshausen et al., 2017) and δ¹³CO₂ (red, Jones et al., 2016) during 1850–2010. (b) The evolution of atmospheric CFC-12 concentrations (Bullister, 2017). The solid blue line indicates the Northern Hemisphere, and the dashed red line indicates Southern Hemisphere.

3 Model results and observations in the late 20th century

Our model generally simulates the physical and biogeochemical state for the present-day ocean well. The detailed model–observation comparisons for the ocean physical variables (e.g. seawater temperature and salinity, Atlantic Meridional Overturning Circulation stream function, CFC-12) and for the ocean biogeochemical tracers (e.g. primary production, nutrients, DIC) are summarised in Appendix Sects. B and C.

In this section, we compare simulated ¹³C between the two simulations Hist_Popp and Hist_Laws and evaluate the two experiments by comparison to observed δ¹³C_POC and δ¹³C_DIC. The observations used here are the surface δ¹³C_POC measurements assembled by Goericke and Fry (1994) and the observed δ¹³C_DIC, for both the surface and the interior ocean, compiled by Schmittner et al. (2013). For the model–observation comparison, we first grid the observed δ¹³C_POC and δ¹³C_DIC horizontally onto a 1 × 1^∘ grid and vertically (only for δ¹³C_DIC) onto the 40 depth layers of the model. Multiple data points in the same grid cell in the same month and year are averaged. Then we bilinearly interpolate the simulated monthly-mean δ¹³C_POC and δ¹³C_DIC over a 1 × 1^∘ grid. To quantitatively compare the performance between Hist_Popp and Hist_Laws and to other ¹³C models, we calculate the spatial correlation coefficient r and the normalised root mean squared error (NRMSE, normalised by the standard deviation that is calculated using all the available measurements of δ¹³C_POC or δ¹³C_DIC during the observational periods) between model results and observation.

A global ocean climatology of pre-industrial δ¹³C_DIC has recently been derived by first estimating the oceanic ¹³C Suess effect (Eide et al., 2017a) and then removing it from the observed δ¹³C_DIC (Eide et al., 2017b). This pre-industrial δ¹³C_DIC estimate has been used to evaluate model performance (Tjiputra et al., 2020). We do not include a δ¹³C_DIC evaluation for the pre-industrial ocean because the historical simulations in this study facilitate the direct comparison to observations in the late 20th century, which is different from Tjiputra et al. (2020), who only include pre-industrial simulations with ¹³C tracers. Moreover, as has already been discussed by Eide et al. (2017a) and is discussed in Sect. 5 of this study, the ¹³C Suess effect is possibly underestimated by the Eide et al. (2017a) approach. This suggests Eide et al. (2017b) likely overestimate the pre-industrial δ¹³C_DIC.

3.1 Isotopic signature of particular organic carbon in the surface ocean

For comparison between Hist_Popp and Hist_Laws, the climatological mean state of δ¹³C_POC is derived by averaging over 1960-1991, the period when most δ¹³C_POC measurements were collected. In Hist_Popp, the climatological annual-mean surface δ¹³C_POC has a global mean value of −22.5 ‰, and it shows a distinct horizontal pattern (Fig. 3a). Less negative values up to −19.3 ‰ are found in the subtropical regions, where alkalinity is typically high and CO₂(aq) is consequently low. This low CO₂(aq) results in a smaller isotope fractionation during carbon fixation by phytoplankton (Eq. 10, Fig. 1) with a biological fractionation factor ${\mathit{ϵ}}_{\mathrm{p}}>-\mathrm{13}$ ‰ (Fig. 3c). Poleward of the subtropical regions, δ¹³C_POC gradually decreases. The reason for this is twofold. First, ϵ_p decreases from −13 ‰ to about −20 ‰ following the increase in CO₂(aq). Second, the thermal effect of equilibrium fractionation causes about 3 ‰ more fractionation in the polar regions than in the tropical and subtropical regions (according to Eqs. 7 and 9). The lowest δ¹³C_POC of about −30 ‰ occurs close to Antarctica where highest surface DIC concentrations are typically found because of the upwelling of deep waters and the reduced air–sea gas exchange by ice cover (Takahashi et al., 2014). The annual range of δ¹³C_POC (Fig. 3e), i.e. the difference between the minimum and the maximum of its climatological monthly-mean annual cycle, is low (<0.5 ‰) in the subtropical regions, and it increases polewards up to ∼9 ‰ in the Southern Ocean, mirroring meridional changes in the annual range of CO₂(aq).

Figure 3 The climatological (1960–1991) annual-mean surface values for Hist_Popp (a, c, e) and Hist_Laws (b, d, f) for δ¹³C_POC (a, b), ϵ_p (c, d), and for the annual range of δ¹³C_POC (e, f). All values are given in per mil (‰).

Compared to Hist_Popp, Hist_Laws shows lower annual-mean surface δ¹³C_POC (Fig. 3b), with a global-mean value of −29.9 ‰ due to more negative ϵ_p (Fig. 3d). This is because ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ (Fig. 1) is always more negative than ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ when the simulated mean growth rates (Fig. C1a and b) are lower than 1 d⁻¹. As ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ increases with growth rate (Eq. 12), we find less negative δ¹³C_POC (up to −24.1 ‰) in the central tropical Pacific, where highest growth rates are simulated (Fig. C1a and b). The lowest δ¹³C_POC of −33 ‰ occurs in the Arctic Ocean and around Antarctica due to the combination of low growth rate, high CO₂(aq) and low seawater temperature. The meridional range of the annual-mean δ¹³C_POC in Hist_Laws (∼9 ‰) is smaller than that of Hist_Popp (∼11 ‰) because for low growth rates ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ is generally less sensitive to CO₂(aq) changes compared to ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ (Fig. 1). This also results in a smaller annual range of δ¹³C_POC in high latitudes (Fig. 3f) than Hist_Popp (Fig. 3e). In the low and mid-latitudes, Hist_Laws shows larger annual range of δ¹³C_POC because in these regions CO₂(aq) concentrations are relatively stable but growth rates shows noticeable seasonal variability.

Hist_Popp captures major features of the observed δ¹³C_POC (Fig. 4a, c and e). The meridional gradient, with less negative values in the low latitudes and minimal values around 60^∘ S, is well reproduced. In contrast, Hist_Laws shows generally lower δ¹³C_POC than the observations (a global mean bias of −8 ‰) and a smaller δ¹³C_POC difference between low and high latitudes (Fig. 4b, d and f). This is also seen in a recent study by Dentith et al. (2020), who tested ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ and ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ with the FAst Met Office/UK Universities Simulator (FAMOUS). The underestimation in the global mean and in the meridional gradient of δ¹³C_POC in Hist_Laws suggests that the parameters of the linear fit in Eq. (12) (slope and intercept) would need to be increased to gain a better performance. Around 60^∘ S of the Atlantic Ocean (Fig. 4b), Hist_Laws simulates a smaller range of δ¹³C_POC than the observations. This is also a result of the small δ¹³C_POC annual range produced by ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ (Fig. 3f). Between 40^∘ S and 40^∘ N in the Atlantic Ocean, Hist_Laws simulates δ¹³C_POC peaks in the region of high growth rates south of the Equator, whereas the observed high δ¹³C_POC values are located between the Equator and 20^∘ N.

Figure 4 Comparison of surface δ¹³C_POC (‰) observations (blue triangle) from Goericke and Fry (1994) to model data (red circle) in Hist_Popp (a, c, e) and Hist_Laws (b, d, f) for the Atlantic, Pacific and Indian Ocean, respectively. Inserted maps show cruise tracks of the measuring campaigns. (g) Comparison of simulated CO₂(aq) (red star) to observations (blue diamond) in the South Indian Ocean (Francois et al., 1993; measurement locations indicated by black triangles in the inset map for the Indian Ocean). Panel (h) is as panel (g) but for particulate organic matter, represented by total POC in Francois et al. (1993) and by phytoplankton biomass in the model. The measurement precision is ±0.17 ‰ for δ¹³C_POC and 2 % for CO₂(aq) and particulate organic matter, according to Francois et al. (1993).

In the Indian Ocean around 45^∘ S, Hist_Popp does not capture the prominent δ¹³C_POC peak in the field data (Fig. 4e), despite the fact that the simulated CO₂(aq), the controlling factor in the parameterisation ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ (Eq. 10), well reproduces the meridional variation of the contemporaneous CO₂(aq) measurements (Fig. 4g). Although the empirical correlation between ϵ_p and CO₂(aq), such as Eq. (10), holds true to the first order over large areas of the global ocean, other factors, such as growth rate, affect the local variability in ϵ_p (Popp et al., 1998; Hansman and Sessions, 2016; Tuerena et al., 2019). Hist_Laws captures the δ¹³C_POC peak around 45^∘ S in the observations (Fig. 4f), owing to the dependency of ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ on phytoplankton growth rate and to the model successfully reproducing the high productivity in this region (illustrated by phytoplankton biomass, Fig. 4h). This is in alignment with the field study by Francois et al. (1993) and the model study by Hofmann et al. (2000), who ascribed this observed δ¹³C_POC peak to a local high phytoplankton production during the measurement period.

Overall, Hist_Popp (r=0.84 and NRMSE=0.57) better reproduces the observed δ¹³C_POC than Hist_Laws (r=0.71, NRMSE=2.5). Here a higher NRMSE indicates the model captures a smaller fraction of the variation in observations. The performance of Hist_Popp regarding δ¹³C_POC compares well to that of the FAMOUS model (Dentith et al., 2020; comparing their Fig. 8 and Fig. 4 in this study) and the University of Victoria (UVic) Earth System Model of intermediate complexity (with r=0.74 and NRMSE=0.92; Schmittner et al., 2013). Note that Schmittner et al. (2013) compared climatological annual-mean model output to the δ¹³C_POC measurements from Goericke and Fry (1994), whereas our study uses model results of the corresponding month and year of the measurements. This difference leads to a better comparison of Hist_Popp to the observed δ¹³C_POC in high latitudes, particularly in the South Atlantic Ocean around 60^∘ S, and therefore it is one reason for the slight better performance of Hist_Popp compared to Schmittner et al. (2013), aside from the underlying differences between the two models.

Hist_Popp also well reproduces the temporal changes of the biological fractionation factor ϵ_p when compared to the observation-based estimates of Young et al. (2013). In Hist_Popp, the change rate of ϵ_p has a global-mean value of −0.026 ‰ yr⁻¹ for the period 1960–2009 (Fig. C7a), similar to an estimate of −0.022 ‰ yr⁻¹ in Young et al. (2013). Modest ϵ_p changes are found in eastern tropical Pacific and south of 60^∘ S, in good agreement with Young et al. (2013). Hist_Laws, on the other hand, shows a very small global-mean ϵ_p change rate of −0.005 ‰ yr⁻¹ (Fig. C7b) as ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ is less sensitive to the increase in CO₂(aq) than ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$.

3.2 Isotopic signature of dissolved inorganic carbon δ¹³C_DIC

3.2.1 Comparison between Hist_Popp and Hist_Laws and to observations

Figures 5a–b and 6a–f compare the climatological annual mean of δ¹³C_DIC (averaged over 1990–2005, when most δ¹³C_DIC measurements were collected) between Hist_Popp and Hist_Laws. The two simulations exhibit very similar δ¹³C_DIC patterns for both the surface and interior ocean. The surface seawater DIC is enriched in ¹³C due to the preferential uptake of the light isotope ¹²C by phytoplankton during primary production. As particulate organic matter sinks and is remineralised at depth, the negative δ¹³C_POC signal is released. Consequently, in both Hist_Popp and Hist_Laws, δ¹³C_DIC at the surface is generally higher than in the ocean interior. At the surface of the equatorial central Pacific, the eastern boundary upwelling systems and the Southern Ocean south of 60^∘ S, lower δ¹³C_DIC (<1.6 ‰) is seen due to the upward transport of the ¹³C-depleted water (Fig. 5a and b). In the interior ocean, we find higher δ¹³C_DIC (>1 ‰) in well-ventilated water masses, in particular the North Atlantic Deep Water (NADW) (Fig. 6a and d). The lowest δ¹³C_DIC values $(<-\mathrm{0.5}$ ‰) occur at depth in tropical and subtropical regions (Fig. 6a–f), where a large amount of organic matter is remineralised.

Figure 5 Climatological (averaged over 1990–2005) annual-mean surface δ¹³C_DIC for Hist_Popp (a) and Hist_Laws (b), respectively. Panels (c) and (d) show the difference in the climatological annual-mean δ¹³C_DIC between Hist_Laws and Hist_Popp, and the difference in the climatological annual range of δ¹³C_DIC between the two simulations, respectively.

The global-mean surface δ¹³C_DIC of the two experiments only differs marginally (1.64 ‰ for Hist_Popp and 1.7 ‰ for Hist_Laws), which is expected as they are run using the same prescribed atmospheric δ¹³CO₂ (Schmittner et al., 2013). Given very similar mean surface DI¹³C, the larger vertical DI¹³C gradients in Hist_Laws, established by more negative δ¹³C_POC (Fig. 3a and b), yield lower DI¹³C concentration at depth. This adjustment of DI¹³C content in the ocean interior takes place during the pre-industrial spin-up phase of the simulations via air–sea ¹³CO₂ exchange (Appendix A). At the end of the 2500-year spin-up, the water-column DI¹³C inventory in PI_Laws is 1.1×10¹² kmol lower than PI_Popp, yielding a global mean δ¹³C_DIC difference of 0.25 ‰ (Fig. 6g–i). Such interior-ocean δ¹³C_DIC difference caused by using different parameterisation for biological fractionation is also seen in Jahn et al. (2015) and Dentith et al. (2020). The seasonal upward transport of the lower deep-ocean δ¹³C_DIC in Hist_Laws leads to lower annual-mean surface δ¹³C_DIC and larger δ¹³C_DIC annual range in regions of upwelling (Fig. 5c and d).

Figure 6 Zonal-mean δ¹³C_DIC of the Atlantic Ocean (a, d, g), the Pacific Ocean (b, e, h) and the Indian Ocean (c, f, i) for Hist_Popp (a–c), Hist_Laws (d–f) and for the difference between Hist_Laws and Hist_Popp (g–i).

When compared to the observed δ¹³C_DIC, Hist_Popp (r=0.81, NRMSE=0.7) has a slightly better performance than Hist_Laws (r=0.80, NRMSE=1.1). Hist_Laws generally shows vertical gradients of δ¹³C_DIC that are too strong and therefore δ¹³C_DIC values that are too low in the ocean interior, as is seen in the depth profiles of horizontally averaged δ¹³C_DIC (Fig. 7). This points to too strong a preference for the isotopically light carbon simulated by ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ as is already discussed in Sect. 3.1. Given the slightly better performance of Hist_Popp than Hist_Laws regarding δ¹³C_DIC, we focus in the following on the comparison between Hist_Popp and observed δ¹³C_DIC.

Figure 7 Depth profiles of horizontally averaged δ¹³C_DIC of Hist_Popp (solid blue line), Hist_Laws (dashed red line) and the observational data from Schmittner et al. (2013) (solid black line) for the global ocean (a), the Atlantic Ocean (b), the Pacific Ocean (c) and for the Indian Ocean (d). The grey shading indicates observation uncertainty of ±0.15 ‰, which relates to the estimated accuracy due to unresolved intercalibration issues between laboratories (0.1 ‰–0.2 ‰; Schmittner et al., 2013).

3.2.2 Source of surface δ¹³C_DIC biases in Hist_Popp

Figure 8 contains model–observation comparison for the surface δ¹³C_DIC. Overall, the magnitude and spatial distribution of the observed δ¹³C_DIC is well-captured by Hist_Popp. In the surface ocean, the mean δ¹³C_DIC is slightly overestimated by Hist_Popp (1.7 ‰ compared to 1.5 ‰ in observation). Positive biases are widely seen in the Indian and Pacific Ocean, and the negative biases are mostly found in the Atlantic Ocean (Fig. 8c). To better understand the source of differences between model and observations, we follow the method of Broecker and Maier-Reimer (1992) to decompose δ¹³C_DIC into a biological component δ¹³C${}_{\mathrm{DIC}}^{\mathrm{bio}}$ and a residual component δ¹³C${}_{\mathrm{DIC}}^{\mathrm{resi}}$, driven by air–sea exchange and ocean circulation:

$$\begin{array}{}\text{(14)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}^{\mathrm{bio}}={\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}{|}_{\mathrm{M}.\mathrm{O}.}+{\displaystyle \frac{{\mathrm{\Delta }}_{\mathrm{photo}}}{{\mathrm{DIC}}_{\mathrm{M}.\mathrm{O}.}}}{\mathrm{R}}_{\mathrm{C}:\mathrm{P}}({\mathrm{PO}}_{\mathrm{4}}-{\mathrm{PO}}_{\mathrm{4}}{|}_{\mathrm{M}.\mathrm{O}.}).\end{array}$$

Here the subscript M.O. refers to mean ocean values, Δ_photo is the carbon isotope fractionation during marine photosynthesis and R_C:P is the C:P ratio of marine organic matter. We use ${\mathrm{\Delta }}_{\mathrm{photo}}=-\mathrm{19}$ ‰ (Eide et al., 2017b) and ${\mathrm{R}}_{\mathrm{C}:\mathrm{P}}=\mathrm{122}$ (Takahashi et al., 1985) for both model and observational data. In reality Δ_photo shows spatial variability due to the variations of CO₂(aq) (Fig. 3c) and temperature (Eq. 7) at the sea surface. However, using a constant Δ_photo only has limited quantitative impact on the model–observation comparison of the two components. To calculate δ¹³C${}_{\mathrm{DIC}}^{\mathrm{bio}}$ from observations, we employ δ¹³C${}_{\mathrm{DIC}}{|}_{\mathrm{M}.\mathrm{O}.}=\mathrm{0.5}$ ‰, ${\mathrm{DIC}}_{\mathrm{M}.\mathrm{O}.}=\mathrm{2255}$ mmol m⁻³ (Eide et al., 2017b) and PO₄ from the World Ocean Atlas (WOA13; Garcia et al., 2013a). Considering the strong seasonality in PO₄ in the surface ocean, we select the phosphate concentration from the climatological monthly WOA data (available only for the upper 500 m of the water column) and the climatological monthly-mean model data for the same month as the δ¹³C_DIC observations. The observed mean ocean phosphate concentration PO${}_{\mathrm{4}}{|}_{\mathrm{M}.\mathrm{O}.}=\mathrm{1.7}$ mmol m⁻³ is obtained by first merging the time mean of the PO₄ monthly WOA data in the upper 500 m and the PO₄ annual-mean WOA data below 500 m and then mapping the combined data to the vertical grid of our model. For simulated δ¹³C${}_{\mathrm{DIC}}^{\mathrm{bio}}$, the model data of δ¹³C${}_{\mathrm{DIC}}{|}_{\mathrm{M}.\mathrm{O}.}=\mathrm{0.67}$ ‰, ${\mathrm{DIC}}_{\mathrm{M}.\mathrm{O}.}=\mathrm{2197}$ mmol m⁻³, ${\mathrm{PO}}_{\mathrm{4}}{|}_{\mathrm{M}.\mathrm{O}.}=\mathrm{1.5}$ mmol m⁻³ and PO₄ are used. The model–observation δ¹³C${}_{\mathrm{DIC}}^{\mathrm{resi}}$ difference is calculated by subtracting the model–observation δ¹³C${}_{\mathrm{DIC}}^{\mathrm{bio}}$ difference from the model–observation δ¹³C_DIC difference.

Figure 8 Observed surface δ¹³C_DIC (Schmittner et al., 2013) (a) and simulated δ¹³C_DIC in Hist_Popp sampled at the location, month and year of the observation (b). (c) The difference in δ¹³C_DIC between Hist_Popp and observations.

The model captures the major features of the observed δ¹³C${}_{\mathrm{DIC}}^{\mathrm{bio}}$ at the surface; that is, higher values are seen in the subtropical regions and lower values in the high latitudes (Fig. C8a and b). Nevertheless, noticeable quantitative differences exist (Fig. 9a), which resemble the distribution of $({\mathrm{PO}}_{\mathrm{4}}-{\mathrm{PO}}_{\mathrm{4}}{|}_{\mathrm{M}.\mathrm{O}.}$) bias (Fig. 9b). Between 30^∘ N and 30^∘ S in the surface ocean, we find a mean negative bias of about −0.1 ‰. This is caused by the underestimation of primary production in the subtropical gyres (due to the underestimation of phytoplankton growth rates; see Appendix C1) and the consequently reduced enrichment of ¹³C in surface DIC. A strong positive δ¹³C${}_{\mathrm{DIC}}^{\mathrm{bio}}$ bias of 0.6 ‰ to 1 ‰ is seen in the North Pacific, where in the model iron is not a limiting nutrient (Fig. C3), in contrast to observations (Moore et al., 2013). In the equatorial central Pacific, a weak positive δ¹³C${}_{\mathrm{DIC}}^{\mathrm{bio}}$ bias <0.2 ‰ is caused by a primary production that is too high. Specifically, the simulated phytoplankton growth rates in this region compare well to observations, whereas the simulated phytoplankton biomass is too high (Appendix C1). The latter is mainly induced by an upwelling that is too strong. The observed mean upward vertical velocity at 0, 140^∘ W and 60 m depth during May 1990–June 1991 is $\mathrm{2.3}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ (Weisberg and Qiao, 2000), whereas the model simulates $\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ for the same location and period.

Figure 9 Model–observation difference in the biological component δ¹³C${}_{\mathrm{DIC}}^{\mathrm{bio}}$ (a), $({\mathrm{PO}}_{\mathrm{4}}-{\mathrm{PO}}_{\mathrm{4}}{|}_{\mathrm{M}.\mathrm{O}.})$ (b), and the residual component δ¹³C${}_{\mathrm{DIC}}^{\mathrm{resi}}$ (c) at the ocean surface. (d) The net air–sea CO₂ flux (positive into the air, averaged over 1990–2005) difference between the model data and observation-based data product from Landschützer et al. (2015).

In the Southern Ocean, a strong positive δ¹³C${}_{\mathrm{DIC}}^{\mathrm{bio}}$ bias of 0.6 to 1 ‰ (Fig. 9a) results from a primary production that is too high under surface iron concentrations that are too high (0.2–0.4 nmol L⁻¹ compared to generally <0.25 nmol L⁻¹ from data of the GEOTRACES programme (http://www.geotraces.org, last access: ​​​​​​​15 April 2021, not shown). Primary production is limited by iron only south of 50^∘ S in the model compared to south of 40^∘ S from observation (Moore et al., 2013). One cause for the high surface iron concentration is that in HAMOCC6 organic matter is remineralised at depths that are too shallow. This can been seen from the positive apparent oxygen utilisation (AOU) biases above 500 m south of 45^∘ S (Fig. 10j–l).

Figure 10 Zonal-mean distribution in the Atlantic Ocean (left column), the Pacific Ocean (middle column) and the Indian Ocean (right column) for the δ¹³C_DIC observations from Schmittner et al. (2013) (a–c), for the difference between Hist_Popp (sampled at the same location, year and month of the observations) and δ¹³C_DIC measurement (d–f), for the $({\mathrm{PO}}_{\mathrm{4}}-{\mathrm{PO}}_{\mathrm{4}}{|}_{\mathrm{M}.\mathrm{O}.})$ difference between model and WOA data (WOA13; Garcia et al., 2013a) (g–i), and for the apparent oxygen utilisation (AOU) difference between model and WOA data (WOA13; Garcia et al., 2013b) (j–l). Here the climatological annual mean values of PO₄ and AOU are used for both model and WOA data because seasonal variation is negligible in the interior ocean and WOA only provides monthly data above 500 m.

Another reason for the high surface iron concentration in the Southern Ocean is that MPIOM simulates an upwelling that is too strong. In particular, below 1000 m, the simulated upward velocity shows noticeably larger magnitude $(>\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹, Fig. B4) than that of a dynamically consistent and data-constrained ocean state estimate (see Fig. 1 in Liang et al., 2017). The upwelling that is too strong in the model is consistent with the volume transport that is too large across the Drake Passage of 192 Sv compared to 134–173 Sv from observations (Nowlin Jr. and Klinck, 1986; Cunningham et al., 2003; Meredith et al., 2011; Donohue et al., 2016). Our model also features larger downward velocities than the estimate from Liang et al. (2017), which correspond to mixed layer depths that are too deep in the Southern Ocean (up to 3000 m, Fig. B5) than observations (<700 m; de Boyer Montégut et al., 2004; Holte et al., 2017).

We find strong δ¹³C${}_{\mathrm{DIC}}^{\mathrm{resi}}$ negative biases of −0.5 ‰ to −1 ‰ (Fig. 9c) in the North Pacific and the Southern Ocean, which partially compensate for the positive biases of δ¹³C${}_{\mathrm{DIC}}^{\mathrm{bio}}$ (Fig. 9a) in these regions. One major cause for the negative δ¹³C${}_{\mathrm{DIC}}^{\mathrm{resi}}$ bias in these two regions is our model overestimating the uptake of anthropogenic carbon, as is illustrated by the net air–sea CO₂ difference between the model and the observation (Fig. 9d). Consequently, the decreased atmospheric ${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ ratio over the industrial period further lowers δ¹³C_DIC in the two ocean regions in the model. In the Southern Ocean, the upward transport of ¹³C-depleted water is too large, which also contributes to a negative δ¹³C${}_{\mathrm{DIC}}^{\mathrm{resi}}$ bias.

3.2.3 Source of δ¹³C_DIC biases in the interior ocean of Hist_Popp

Figure 10 contains the model–observation comparison for zonal-mean δ¹³C_DIC in the Atlantic, Pacific and Indian Ocean. In the interior ocean, δ¹³C_DIC is controlled by remineralisation of ¹³C-depleted organic matter and by ocean circulation (Broecker and Peng, 1993; Lynch-Stieglitz et al., 1995; Schmittner et al., 2013). Low δ¹³C_DIC is often found in waters of high nutrient concentration and vice versa. Thus, we find positive (negative) δ¹³C_DIC biases coincide with negative (positive) phosphate biases (Fig. 10d–i). In the Atlantic Ocean between 1000 and 3000 m, the North Pacific above 1500 m and the Indian Ocean below 1000 m, positive δ¹³C_DIC biases and negative phosphate biases are mainly caused by a remineralisation that is too low, as is shown by the negative AOU biases (Fig. 10j–l).

North of 30^∘ S in the Atlantic Ocean, the negative δ¹³C_DIC biases below 3000 m, together with the positive δ¹³C_DIC biases between 1000 and 3000 m, suggest δ¹³C_DIC vertical gradients that are too strong in the model (Fig. 10d). This results from a lower boundary of the NADW cell that is too shallow, constantly located above 2800 m (Fig. B3), compared to an estimated NADW lower boundary of about 4300 m deep at 26^∘ N (Msadek et al., 2013; Smeed et al., 2017). A possible reason for the shallow NADW in the model is that the Lower North Atlantic Deep Water (LNADW), forming from the Denmark Strait Overflow Water and the Iceland–Scotland Overflow Water, is not dense enough to flow further southward. This can be seen from the CFC-12 distribution along the zonal Sect. A5 at 24^∘ N (Fig. B7). The observed deeper CFC-12 maximum (3000–4500 m west of 60^∘ W) indicates the presence of LNADW (Dutay et al., 2002), which is not represented in our model.

We find the strongest negative δ¹³C_DIC bias in the deep eastern equatorial Pacific (Fig. 10e). The cause is the “nutrient trapping” problem in the model, characterised by nutrient concentrations that are too high in the deep eastern equatorial Pacific (Fig. 10h), which is a persistent problem in many ESMs (Aumont et al., 1999; Dietze and Loeptien, 2013). Based on sensitivity experiments with the Geophysical Fluid Dynamics Laboratory model and the UVic model, Dietze and Loeptien (2013) concluded the primary cause of the nutrient trapping problem is likely model biases in physical ocean state – in particular, the poor representation of the Equatorial Intermediate Current system and equatorial deep jets. The latter two current systems are indeed poorly represented in our model as well. Specifically, the zonal current at 1000 m depth (typical depth for the Equatorial Intermediate Current system) shows too little spatial variability and too low speeds of ∼0.2 cm s⁻¹ (Fig. B6), compared to the observed alternating jets with a meridional scale of 1.5^∘ and speeds of ∼5 cm s⁻¹ (see Fig. 2 from Cravatte et al., 2012).

The performances of both Hist_Popp and Hist_Laws regarding δ¹³C_DIC are comparable with the Norwegian Earth System Model version 2 (NorESM2, Tjiputra et al., 2020; comparing their Fig. 21), the Commonwealth Scientific and Industrial Research Organisation Mark 3L climate system model with the Carbon of the Ocean, Atmosphere and Land (CSIRO Mk3L-COAL), Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) and LOch-Vecode-Ecbilt-CLio-agIsm Model (LOVECLIM) (see Table 2 and Figs. 3, S2 and S3 of Buchanan et al., 2019, and references therein), the Community Earth System Model (CESM, Jahn et al., 2015; comparing their Figs. 5 and 6 to our Figs. 7 and 6, respectively), and the UVic Earth System Model (Schmittner et al., 2013). The latter two studies used the same δ¹³C_DIC data set for model evaluation. Schmittner et al. (2013) reported a better performance (r=0.88 and NRMSE=0.5) than ours (r=0.81 and NRMSE=0.7 in Hist_Popp). One main reason is that the nutrient trapping problem in HAMOCC6 does not occur in the simulations of Schmittner et al. (2013).

4 Evaluation of the simulated oceanic ¹³C Suess effect

The oceanic δ¹³C measurements taken during the late 20th century already include a signal that originates from burning of isotopically light fossil fuel over the industrial period. The associated decrease in atmospheric δ¹³C (Fig. 2) affects oceanic δ¹³C via air–sea gas exchange, leading to a general decrease in δ¹³C_DIC. The distribution of this δ¹³C_DIC change, i.e. the oceanic ¹³C Suess effect, could serve as benchmark for ocean models to evaluate the uptake and re-distribution of the anthropogenic CO₂ emissions in the ocean.

The model is able to reproduce the size of the global oceanic anthropogenic CO₂ sink, though some local biases in the net air–sea CO₂ flux exist (Fig. 9d). The simulated sink by year 1994 is 99 Pg C, which compares well to the observation-based estimate of 118±19 Pg C from Sabine et al. (2004) and to other model estimates (e.g. 94 Pg C in Tagliabue and Bopp, 2008). For a direct comparison to published studies, we calculate the oceanic δ¹³C Suess effect, δ¹³C_SE, as the difference between the 1990s-averaged δ¹³C_DIC from Hist_Popp and the pre-industrial climatological (50-year mean) δ¹³C_DIC from PI_Popp. δ¹³C_SE calculated using the results of Hist_Laws and PI_Laws only shows a marginal difference (global mean <0.04 ‰) and is therefore not presented.

The surface mean δ¹³C_SE in this study is −0.66 ‰, similar to the model study of Schmittner et al. (2013) (−0.67 ‰) and to the estimate by Sonnerup et al. (2007) $(-\mathrm{0.76}\pm \mathrm{0.12}$ ‰), who used an observation-based approach. The strongest oceanic ¹³C Suess effect is found in the subtropical gyres in the model (Fig. 11a), where water masses have long residence times at the ocean surface and therefore receive a strong anthropogenic imprint (Quay et al., 2003). In the subtropical gyres, the simulated surface δ¹³C_SE generally varies between −0.8 ‰ and −1.1 ‰, which compares well to the surface ocean δ¹³C decrease of $-\mathrm{0.9}\pm \mathrm{0.1}$ ‰ recorded by coral and sclerosponges (Wörheide, 1998; Böhm et al., 1996, 2000; Swart et al., 2002, 2010) and to the estimates of $-\mathrm{1.0}\pm \mathrm{0.09}$ ‰ extracted from GLODAPv2 (Olsen et al., 2016; Eide et al., 2017a).

Figure 11 The simulated oceanic Suess effect δ¹³C_SE from the pre-industrial period to the 1990s at sea surface (a) and at 200 m (b).

Along the vertical sections A16, P19 and I8S9N, δ¹³C_SE is mainly confined to the upper 1000 m depth in the subtropical gyres of the South Atlantic, the Pacific Ocean and the Indian Ocean (Fig. 12a–c). In the North Atlantic, δ¹³C_SE penetrates deeper than the other ocean regions, due to the intensive ventilation related to the formation of NADW. The simulated δ¹³C_SE distributions show similar features to those of CFC-12 (Fig. B8). This is because both the decrease in δ¹³C_DIC and increase in CFC-12 in the ocean are predominantly caused by the uptake of atmospheric anthropogenic signals and the subsequent transport by ocean circulation. Since changes in δ¹³C_DIC are also induced by changes in marine biological activity, we separate δ¹³C_DIC into a component depicting changes due to the transport of the surface ¹³C signal, i.e. the “preformed” δ¹³C_DIC, and to a regenerated component δ¹³C^reg, following Sonnerup et al. (1999):

$$\begin{array}{}\text{(15)}& {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}^{\mathrm{pref}}={\displaystyle \frac{{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{DIC}}\cdot \mathrm{DIC}-\mathrm{AOU}\cdot {\left(\frac{\mathrm{C}}{{-\mathrm{O}}_{\mathrm{2}}}\right)}_{\mathrm{org}}\cdot {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{org}}}{\mathrm{DIC}-\mathrm{AOU}\cdot {\left(\frac{\mathrm{C}}{{-\mathrm{O}}_{\mathrm{2}}}\right)}_{\mathrm{org}}}}.\end{array}$$

The ${\left(\frac{\mathrm{C}}{{-\mathrm{O}}_{\mathrm{2}}}\right)}_{\mathrm{org}}$ ratio is 122:172 in HAMOCC6, and we use the simulated δ¹³C_POC for δ¹³C_org. Clearly, the change in the preformed component ${{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}}^{\mathrm{pref}}={{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1990}\mathrm{s}}^{\mathrm{pref}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{PI}}^{\mathrm{pref}}$ dominates δ¹³C_SE (comparing Fig. 12a–c to Fig. 12d–f). A major difference between δ¹³C${}_{\mathrm{SE}}^{\mathrm{pref}}$ and δ¹³C_SE is that positive δ¹³C${}_{\mathrm{SE}}^{\mathrm{pref}}$ is widely seen below 1000 m, particularly in the Pacific Ocean (Fig. 12e). These positive δ¹³C${}_{\mathrm{SE}}^{\mathrm{pref}}$ values relate to changes in the regenerated component δ¹³C^reg (see Appendix D).

Figure 12 The simulated oceanic Suess effect δ¹³C_SE since pre-industrial times for vertical sections A16 in the Atlantic Ocean (a), P16 in the Pacific Ocean (b) and I8S9N in the Indian Ocean (c). Panels (d)–(f) and (g)–(i) are as panels (a)–(c) but for the change in the preformed component ${{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}}^{\mathrm{pref}}={{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1990}\mathrm{s}}^{\mathrm{pref}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{PI}}^{\mathrm{pref}}$ and for the observation-based estimate of the oceanic Suess effect from Eide et al. (2017a), respectively. Inserted maps show the location of the vertical sections. The horizontal dashed black lines in panels (a)–(f) indicate 200 m depth, below which the Eide et al. (2017a) estimate is available. Note the bathymetry is different between the model and Eide et al. (2017a).

5 Potential sources of uncertainties in an observation-based global oceanic ¹³C Suess effect estimate

Eide et al. (2017a) (hereafter E17) derived the first observation-based estimate of the global ocean ¹³C Suess effect since pre-industrial times. E17's approach uses the concept of the similarity between the oceanic uptake of the anthropogenically produced CFC-12 and isotopically light CO₂ (see details in Appendix E1). Due to method- and data-specific limitations E17 stated that they potentially underestimate the oceanic ¹³C Suess effect. However, based on observations alone it is not possible to gain insight into the spatial distribution of this uncertainty or into its origin.

Our model simulations, particularly PI_Popp and Hist_Popp, provide an opportunity to learn more about the source of this uncertainty because the oceanic δ¹³C in the late 20th century (Sect. 3), the oceanic anthropogenic CO₂ sink (Sect. 4) and the invasion of CFC-12 into the ocean (Fig. B8) are well represented. Moreover, our simulated δ¹³C_SE qualitatively resembles the oceanic ¹³C Suess effect estimate of E17 (see comparison between Fig. 11b and E17's Fig. 7, as well as comparison between Fig. 12a–c and g–i).

Based on the similarity between the oceanic uptake of the atmospheric CFC-12 and δ¹³CO₂ signal, E17 link the ¹³C Suess effect since 1940 (when CFC-12 becomes detectable in the ocean) to CFC-12 partial pressure (pCFC-12) with a proportionality factor. Under the assumption of a temporally constant regenerated fraction δ¹³C^reg, this proportionality factor is considered equivalent to the slope of a linear regression relationship between the preformed component δ¹³C^pref and pCFC-12 at any time after 1940. Thus, this slope a can be obtained by performing linear regression for field measurements of δ¹³C^pref and pCFC-12. Multiplying a and pCFC-12 data yields the ¹³C Suess effect since 1940, which is then scaled to the full industrial period by a constant factor f_atm (Eq. E7) related to changes in the atmospheric δ¹³C signature:

$$\begin{array}{}\text{(16)}& {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}(t-\mathrm{PI})}=f_{\mathrm{atm}}\cdot a\cdot {\text{pCFC-12}}_t.\end{array}$$

Here a is the regression slope for the linear relationship between δ¹³C${}_t^{\mathrm{pref}}$ and pCFC-12_t (Eq. E5). The value of a is determined for each ventilation region defined in E17 (i.e. the Indian Ocean, North Pacific, South Pacific, North Atlantic and South Atlantic). Details of the E17 approach are given in Appendix E1.

By applying E17's approach to our model data that are sampled at the same geographical locations as observations used in E17, we obtain the regression slopes, hereafter referred to as a_pref, for each ventilation region. Taking year t=1994 we obtain the estimated oceanic ¹³C Suess effect, SE_pref, for the period from the pre-industrial to 1994 following Eq. (16). The detailed calculation of SE_pref is given in Appendix E2.

To quantify if SE_pref under- or overestimates the oceanic ¹³C Suess effect, we compare SE_pref to the simulated oceanic ¹³C Suess effect ${\mathrm{SE}}_{\mathrm{Mod}}={{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{DIC},\mathrm{1994}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{DIC},\mathrm{PI}}$. Figure 13a presents (SE_pref−SE_Mod) for 200 m depth. Positive values of(SE_pref−SE_Mod) indicate underestimation of the oceanic ¹³C Suess effect.

Figure 13 Distribution at 200 m depth for SE_pref−SE_Mod (a), SE_pref−SE_total (b) and SE_total−SE_Mod (c). The isoline increment is 0.2 ‰. In panels (b) and (c), the South Pacific Ocean is not presented because the relationship between the total oceanic ¹³C Suess effect δ¹³C_{SE(1994–1940)} and pCFC-12₁₉₉₄ is too weak (r²=0.07), and therefore SE_total can not be estimated (see Appendix E4).

At 200 m SE_pref mostly underestimates SE_Mod (Fig. 13a). The region-mean underestimation is 0.24 ‰ for the Indian Ocean, 0.21 ‰ for the North Pacific, 0.26 ‰ for the South Pacific, 0.1 ‰ for the North Atlantic and 0.14 ‰ for the South Atlantic (Table 1). Our model findings are very similar to the underestimation range discussed by E17. They determined an uncertainty range of 0.15 ‰ to 0.24 ‰ by comparing their global-mean estimate (−0.4 ‰ at 200 m depth) to an estimate (−0.55 ‰ to −0.64 ‰ at 200 m) which they deduced from previous model studies. Specifically, based on Broecker and Peng (1993) and Bacastow et al. (1996) E17 assumed an ocean-to-atmosphere ratio of the ¹³C Suess effect of 0.65 and the 200 m to surface ratio of the ¹³C Suess effect of 0.6–0.7. Multiplying the above two ratios with the atmospheric δ¹³CO₂ decrease of −1.4 ‰ by year 1994 yields the global-mean ¹³C Suess effect estimate of −0.55 ‰ to −0.64 ‰ at 200 m. In our model, the global-mean surface-ocean–atmosphere ratio of the ¹³C Suess effect is only 0.46, significantly lower than the five-box model of Broecker and Peng (1993). The 200 m to surface ratio of the ¹³C Suess effect is 0.75 in our model, and it is slightly higher than that of Bacastow et al. (1996), who employed an ocean general circulation model with coarse vertical resolution (four layers for the upper 200 m).

Table 1 Region-mean (SE_pref−SE_Mod), (SE_pref−SE_total) and (SE_total−SE_Mod) for five ventilation regions defined by E17, i.e. the Indian Ocean, North Pacific, South Pacific, North Atlantic and South Atlantic. The unit is per mil. (SE_pref−SE_total) is further decomposed into the two contributions $f_{\mathrm{atm}}\cdot (a_{\mathrm{pref}}-a_{\mathrm{total}})\cdot \text{pCFC-12}$ and $-f_{\mathrm{atm}}\cdot b_{\mathrm{total}}$ according to Eq. (20).

5.1 Source of underestimation attributed to data coverage

E17 have speculated that the major cause of the underestimation of the oceanic ¹³C Suess effect is that the available observations are mostly from the intermediate and deep waters. The ocean–atmosphere equilibration timescale for δ¹³C (10 years, Broecker and Peng, 1974) is significantly longer than that of pCFC-12 (1 month, Gammon et al., 1982). Thus, waters that have a shorter surface residence time, such as the deep waters ventilated in the South Hemisphere, would show less negative regression slope a_pref (for the linear relationship between δ¹³C^pref and pCFC-12, Eq. E5) than waters that have a longer surface residence time, e.g. subtropical gyres. In other words, a_pref for the subtropical gyre water should be more negative than a_pref for the entire corresponding ventilation region (the North Pacific, South Pacific, North Atlantic, South Atlantic or the Indian Ocean).

We test this potential explanation for the Indian Ocean and North Pacific. We are able to span regressional relationships for the subtropical gyres only because we have a larger data base. Specifically, we consider only model data at the geographical location of observations, but we use all model levels between 200 m and the pCFC-12 penetration depth (see Appendix E3). For the Indian Ocean, we combine model data from Subtropical Gyre Water (STGW) and Sub-Antarctic Mode Water (SAMW) as both water masses have a strong ¹³C Suess effect (Eide et al., 2017a). We find for this combined water mass (STGW) a_pref ($-\mathrm{0.65}\times {\mathrm{10}}^{-\mathrm{3}}$, r²=0.49) is more negative than that for the whole ventilation region $(-\mathrm{0.47}\times {\mathrm{10}}^{-\mathrm{3}}$, Fig. E3a). So indeed, with additional observations in the subtropical gyre we would receive a stronger ¹³C Suess effect estimate for the Indian Ocean. However, this difference in a_pref only corresponds to an underestimation of about 0.12 ‰ at 200 m for the Indian subtropical region (see calculation in Appendix E3), which does not explain the total underestimation of 0.24 ‰ in the Indian Ocean (Table 1). In the North Pacific a_pref for the Subtropical Gyre Water $(-\mathrm{0.44}\times {\mathrm{10}}^{-\mathrm{3}}$, r²=0.26) is even less negative than that for the whole ventilation region $(-\mathrm{0.71}\times {\mathrm{10}}^{-\mathrm{3}})$ in the model, which is in contrast to the conjecture of E17.

5.2 Source of underestimation attributed to assumptions of E17's approach

A potential under-representation of data from subtropical gyres does not fully explain the underestimation of the ¹³C Suess effect found in our model. Instead, we argue that the source of uncertainty mainly relates to different assumptions that have been made in the E17 approach. Specifically, in the expression of the preformed component δ¹³C${}_{\mathrm{1994}}^{\mathrm{pref}}$ (following Eq. E3),

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1994}}^{\mathrm{pref}}& ={{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}(\mathrm{1994}-\mathrm{1940})}+{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{pref}}\\ & -({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1994}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}}),\end{array}\end{array}$$

E17 assume that the regenerated component is constant in time, i.e. $-({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1994}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}})=\mathrm{0}$. Consequently, Eq. (17) is reduced to

$$\begin{array}{}\text{(18)}& {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1994}}^{\mathrm{pref}}={{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}(\mathrm{1994}-\mathrm{1940})}+{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{pref}}.\end{array}$$

Furthermore, they assume that the regression slope a_pref for δ¹³C${}_{\mathrm{1994}}^{\mathrm{pref}}$ and pCFC-12₁₉₉₄ is equivalent to the regression slope for the total ¹³C Suess effect δ¹³C_{SE(1994−1940)} and pCFC-12₁₉₉₄ (see Eqs. E1, E4 and E5). This implies that the preformed component δ¹³C${}_{\mathrm{1940}}^{\mathrm{pref}}$ of 1940 has to be spatially uniform.

However, we find a specific vertical structure in the simulated δ¹³C${}_{\mathrm{1940}}^{\mathrm{pref}}$ (Fig. 14a–c). Over large regions of the ocean, δ¹³C${}_{\mathrm{1940}}^{\mathrm{pref}}$ generally decreases with increasing depth. This vertical distribution of δ¹³C^pref is already present in pre-industrial times. High surface δ¹³C_DIC caused by biological fractionation is transported into the ocean interior. Therefore, the preformed component generally decreases with increasing water depth. From pre-industrial times to 1940, the decrease in the atmospheric ${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ ratio is relatively small (0.4 ‰, Fig. 2a), and therefore also the impact on the oceanic δ¹³C_DIC is small. Thus, δ¹³C${}_{\mathrm{1940}}^{\mathrm{pref}}$ has the similar vertical structure as that of the pre-industrial ocean.

Figure 14 (a–c) The zonal mean of the simulated δ¹³C${}_{\mathrm{1940}}^{\mathrm{pref}}$ for the locations where both observed CFC-12 and δ¹³C_DIC are available. The thick grey line is the pCFC-12₁₉₉₄=20 patm isoline, above which model data are used to perform linear regression. The thick black lines outline the Subtropical Gyre Water in the Atlantic and North Pacific Ocean, as well as the Subtropical Gyre Water and Sub-Antarctic Mode Water in the Indian Ocean and South Pacific Ocean (definition of water masses in Table E1). Panels (d)–(f), (g)–(i) and (j)–(l) are as panels (a)–(c) but for pCFC-12₁₉₉₄, for $-({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1994}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}})$ and for AOU changes between year 1940 and 1994, respectively. Note that for the Atlantic Ocean the upper 3 km is shown, whereas for the Pacific and Indian Ocean the upper 1.5 km is presented.

Both the total δ¹³C_{SE(1994−1940)} (mostly negative, similar to the distribution of δ¹³C_{SE(1990s−PI)} in Fig. 12a–c) and pCFC-12 (Fig. B8a–c) show larger absolute values at the surface than in the interior ocean. As δ¹³C${}_{\mathrm{1940}}^{\mathrm{pref}}$ is more positive in the upper ocean than the deep ocean, δ¹³C${}_{\mathrm{1994}}^{\mathrm{pref}}$ has a smaller vertical gradient than δ¹³C_{SE(1994−1940)} (see Eq. 18). Thus, a linear regression for δ¹³C${}_{\mathrm{1994}}^{\mathrm{pref}}$ and pCFC-12 results in a less negative slope than a slope obtained with a spatially uniform δ¹³C${}_{\mathrm{1940}}^{\mathrm{pref}}$, which implicates a contribution to an underestimation of the oceanic ¹³C Suess effect.

We also find that $-({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1994}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}})$ is non-zero, and it shows considerable spatial variability (Fig. 14g–i). Most prominently, in the North Atlantic $-({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1994}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}})$ is mostly negative above 500 m, and it is mostly positive below 500 m. This vertical structure of $-({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1994}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}})$ in the North Atlantic leads to stronger vertical gradient in δ¹³C${}_{\mathrm{1994}}^{\mathrm{pref}}$ and therefore a more negative regression slope than that obtained with $-({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1994}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}})=\mathrm{0}$. This implies the overestimation of the ¹³C Suess effect in the North Atlantic.

To evaluate the impact of assuming a spatially uniform δ¹³C${}_{\mathrm{1940}}^{\mathrm{pref}}$ and $-({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1994}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}})=\mathrm{0}$, we calculate an estimated ¹³C Suess effect from pre-industrial times to 1994, SE_total, based on a linear regression for the simulated total oceanic ¹³C Suess effect δ¹³C_{SE(1994−1940)} and pCFC-12:

$$\begin{array}{}\text{(19)}& {\mathrm{SE}}_{\mathrm{total}}=f_{\mathrm{atm}}\cdot (a_{\mathrm{total}}\cdot {\text{pCFC-12}}_{\mathrm{1994}}+b_{\mathrm{total}}).\end{array}$$

Here a_total and b_total are regression coefficients for δ¹³C_{SE(1994−1940)} and pCFC-12 (more details in Appendix E4). With Eqs. (16) and (19) we get

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\mathrm{SE}}_{\mathrm{pref}}-{\mathrm{SE}}_{\mathrm{total}}& =f_{\mathrm{atm}}\cdot (a_{\mathrm{pref}}-a_{\mathrm{total}})\\ & \cdot {\text{pCFC-12}}_{\mathrm{1994}}-f_{\mathrm{atm}}\cdot b_{\mathrm{total}}.\end{array}\end{array}$$

Comparison between the regressional slope a_pref (obtained for δ¹³C${}_{\mathrm{1994}}^{\mathrm{pref}}$ and pCFC-12) and a_total facilitates the quantification of the under- or overestimation of the ¹³C Suess effect linked to the above two assumptions.

In the Indian Ocean $a_{\mathrm{pref}}=-\mathrm{0.47}\times {\mathrm{10}}^{-\mathrm{3}}$ (Fig. E2a) is less negative than $a_{\mathrm{total}}=-\mathrm{0.74}\times {\mathrm{10}}^{-\mathrm{3}}$ (Fig. E3a). This results in an underestimation of 0.12 ‰ according to Eq. (20). Similarly, for the North Pacific $a_{\mathrm{pref}}=-\mathrm{0.71}\times {\mathrm{10}}^{-\mathrm{3}}$ (Fig. E2b) is less negative than $a_{\mathrm{total}}=-\mathrm{0.83}\times {\mathrm{10}}^{-\mathrm{3}}$ (Fig. E3b), which leads to an underestimation of 0.06 ‰. For the South Atlantic $a_{\mathrm{pref}}=-\mathrm{0.6}\times {\mathrm{10}}^{-\mathrm{3}}$ (Fig. E2e) and $a_{\mathrm{total}}=-\mathrm{0.7}\times {\mathrm{10}}^{-\mathrm{3}}$ (Fig. E3e), which yields an underestimation of 0.04 ‰. Such underestimation is mainly due to the decreasing δ¹³C${}_{\mathrm{1940}}^{\mathrm{pref}}$ with increasing depth in these regions.

Different from these three ventilation regions, in the North Atlantic $a_{\mathrm{pref}}=-\mathrm{0.81}\times {\mathrm{10}}^{-\mathrm{3}}$ (Fig. E2d) is more negative than $a_{\mathrm{total}}=-\mathrm{0.62}\times {\mathrm{10}}^{-\mathrm{3}}$ (Fig. E3d). This is due to the specific vertical structure of $-({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1994}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}})$ as previously discussed.

Another major difference between SE_pref and SE_total is the non-negligible negative intercept b_total (Eq. 20). This reveals the underestimation of SE_pref related to E17's assumption that the ¹³C Suess effect is directly proportional to pCFC-12. The intercept b_total emerges possibly due to the different atmospheric time history of the ¹³C Suess effect compared to CFC-12, as is discussed by E17 for the deep ocean with very low or zero CFC-12. The decreasing of δ¹³C_POC under increasing surface CO₂(aq) (Appendix D) also contributes to an non-negligible b_total as lower δ¹³C_POC leads to lower δ¹³C_DIC in the ocean interior. In the South Atlantic and Indian Ocean, $b_{\mathrm{total}}=-\mathrm{0.07}$ ‰ corresponds to an underestimation of 0.12 and 0.11 ‰ (Table 1), respectively.

Table 1 summaries of the contributions from (SE_pref−SE_total) for different ventilation regions. The comparison to the total underestimation given by (SE_pref−SE_Mod) shows that this underestimation, which is attributed to the assumption of E17's approach, is the largest contributor for the Indian Ocean and the South Atlantic.

The residual under-/overestimation of SE_pref given by $({\mathrm{SE}}_{\mathrm{total}}-{\mathrm{SE}}_{\mathrm{Mod}})=({\mathrm{SE}}_{\mathrm{pref}}-{\mathrm{SE}}_{\mathrm{Mod}})-({\mathrm{SE}}_{\mathrm{pref}}-{\mathrm{SE}}_{\mathrm{total}})$ shows how well a method based on linear regression relationships between δ¹³C_SE and pCFC-12₁₉₉₄ can estimate the global ocean Suess effect. (SE_total−SE_Mod) at 200 m generally show positive values, i.e. underestimation, in low latitudes (between 40^∘ S and 40^∘ N), and it is rather negative poleward of 40^∘ (Fig. 13c). This pattern results from pooling data from different water masses to generate one regression relationship for a large ventilation region. The waters ventilated in lower latitudes typically have a stronger ¹³C Suess effect than those ventilated in high latitudes. This is clearly reflected in the linear regression relationships between δ¹³C_{SE(1994–1940)} and pCFC-12₁₉₉₄ for the North Atlantic (Fig. E3d), which shows that the regression slope a_total for the Subtropical Gyre Water is noticeably steeper than that of the deep waters. Accordingly in the interior ocean, the water masses ventilated in the low latitudes generally show an underestimation of the ¹³C Suess effect (positive values of SE_total−SE_Mod), and the water masses ventilated in the high latitudes show an overestimation (Fig. E1g–i). In the North Atlantic Ocean, the region-mean underestimation $({\mathrm{SE}}_{\mathrm{pref}}-{\mathrm{SE}}_{\mathrm{Mod}})=\mathrm{0.1}$ ‰ is predominantly contributed by $({\mathrm{SE}}_{\mathrm{total}}-{\mathrm{SE}}_{\mathrm{Mod}})=\mathrm{0.09}$ ‰. In the North Pacific Ocean $({\mathrm{SE}}_{\mathrm{total}}-{\mathrm{SE}}_{\mathrm{Mod}})=\mathrm{0.13}$ ‰ accounts for more than half of the total underestimation 0.21 ‰. In the Indian and South Atlantic Ocean, however, (SE_total−SE_Mod) has hardly any influence to the region-mean underestimation.

In summary, our analysis points out two major causes for the underestimation of ¹³C in E17's approach. The first is the assumption of a spatially uniform preformed δ¹³C component in 1940. The second cause is the neglect of processes not directly linked to the oceanic uptake and transport of CFC-12, e.g. the uptake of anthropogenically light CO₂ in the times prior to the emission of CFC-12 and the decrease in δ¹³C_DIC due to the decrease in δ¹³C_POC over the industrial period.

6 Summary and conclusions

We present results of the new ¹³C module in the ocean biogeochemical model HAMOCC6 for the historical period forced by reanalyses data (ERA20C). We test two parameterisations of different complexity for the biological fractionation factor: ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ depends on dissolved CO₂ (Popp et al., 1989); ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ is a function of dissolved CO₂ and phytoplankton growth rate (Laws et al., 1995). Furthermore, we use our consistent model framework to assess the approach by Eide et al. (2017a), which yields the first global oceanic ¹³C Suess effect estimate based on a correlation between preformed δ¹³C_DIC and CFC-12 partial pressure.

The comparison between simulated and observed isotopic ratio of organic matter δ¹³C_POC reveals that ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ (r=0.84 and NRMSE=0.57) has a better performance than ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ (r=0.71 and NRMSE=2.5). Using ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ results in noticeably lower δ¹³C_POC values and smaller δ¹³C_POC gradients between low and high latitudes compared to observations. The parameterisation of Laws et al. (1995), obtained based on cultures of marine diatom Phaeodactylum tricornutum, results in too strong a preference of isotopically light carbon in our global ocean biogeochemical model.

Regarding δ¹³C_DIC, ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ also yields slightly better agreement with observations than ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{laws}}$ (r=0.81 and NRMSE=0.7 versus r=0.80 and NRMSE=1.1), because ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Laws}}$ produces lower δ¹³C_POC and therefore lower interior-ocean δ¹³C_DIC than those found in observations. ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ performs well considering the uncertainties in observed δ¹³C_DIC (0.1 ‰–0.2 ‰; Schmittner et al., 2013). Our model slightly overestimates surface δ¹³C_DIC. By decomposing δ¹³C_DIC into a biological component and a residual component, we find the overestimation in the high-latitude ocean is dominated by biases in the biological component caused by e.g. surface iron concentration that is too high. In the interior-ocean δ¹³C_DIC biases are mainly due to biases in the physical state (for instance, a boundary that is too shallow between the NADW cell and the Antarctic Bottom Water cell in MPIOM).

Our model represents well the temporal evolution of the oceanic δ¹³C_DIC since pre-industrial times, i.e. the oceanic ¹³C Suess effect due to the intrusion of isotopically light carbon into the ocean. With the complete information on the spatial and temporal ¹³C evolution in the ocean, together with the simulated evolution of CFC-12, we identify the sources for the potential uncertainties in the framework of Eide et al. (2017a) for deriving an observation-based oceanic ¹³C Suess effect. Based on our model, we find underestimations of the ¹³C Suess effect at 200 m by 0.24 ‰ in the Indian Ocean, 0.21 ‰ in the North Pacific Ocean, 0.26 ‰ in the South Pacific Ocean, 0.1 ‰ in the North Atlantic Ocean and 0.14 ‰ in the South Atlantic Ocean. These numbers are in line with the underestimation range 0.15 ‰ to 0.24 ‰ conjectured by Eide et al. (2017a). They speculated this underestimation is due to the under-representation of the water masses with a stronger ¹³C Suess effect, such as the Subtropical Gyre Water and Sub-Antarctic Mode Water, in the observational data. Our analysis shows that their hypothesis only explain half of the underestimation in the Indian Ocean. For the North Atlantic Ocean this hypothesis is not supported by the model data . We identify two major causes for the underestimation of the ¹³C Suess effect by the applied method. The first relates to the assumption of a spatially uniform preformed component of δ¹³C_DIC in year 1940. In our model this preformed component is generally more positive in the upper ocean than in the interior ocean, which contributes to the underestimation of the δ¹³C Suess effect. The second cause relates to the neglect of processes that are not directly linked to the oceanic uptake and transport of CFC-12 – for instance, the ¹³C Suess effect prior to the emission of CFC-12 and the decrease in δ¹³C_POC over the industrial period.

We conclude that the new ¹³C module with biological fractionation factor ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ from Popp et al. (1989) has a satisfactory performance. We are aware that the parameterisation ${\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{Popp}}$ omits any potential changes, e.g. in ecosystem structure, which might have occurred in the paleo-ocean. Our new ¹³C module will serve as a useful tool to evaluate the performance of MPI-ESM in paleoclimate and to investigate the past changes in the ocean, for instance within the ongoing research project PalMod (Latif et al., 2016).

Appendix A: Governing factors for the water-column DI¹³C inventory changes

The water-column DI¹³C inventory difference is primarily a result of the difference in the net air–sea ¹³CO₂ flux between PI_Popp and PI_Laws. This is demonstrated by the comparison of the contributions of the governing factors for the water-column DI¹³C inventory changes (Table A1), including air–sea gas exchange, loss of POC and CaCO₃ to marine sediment, diffusion of the remineralised DIC from sediment into the water column, input of DOC and CO${}_{\mathrm{3}}^{\mathrm{2}-}$, and the exchange with other marine carbon pools (phytoplankton, CaCO₃, etc.). Table A1 also reveals that the current method to determine the ¹³C input (see Sect. 2.3.2) only has a small contribution to the change in the water-column DI¹³C inventory.

Table A1 Contributions to the rate of the water-column DI¹³C inventory change (in Gmol yr⁻¹), averaged in the last 50 years in the corresponding pre-industrial spin-up simulations. Positive values denote contributions to the increase in the water-column DI¹³C inventory. The last column gives the relative contribution to the total rate difference with relative contribution = (PI_Laws-PI_Popp) / total rate difference.

Appendix B: Model–observation comparison of ocean physics

Sea surface temperature (SST) and salinity (SSS) generally show good performance (Fig. B1 and Table B1). The most striking bias is seen for SSS (2–3 psu) in the Arctic Ocean. In the ocean interior, the performance of temperature and salinity is similar to other ocean general circulation models, e.g. Tjiputra et al. (2020) (comparing our Table B1 to their Fig. 2). The model biases shown here (for instance, the surface layers are too cold, whereas the water between 500 and 2500 m is too warm and saline; see Fig. B2) are typically seen in MPIOM; see Jungclaus et al. (2013) for a detailed discussion.

Figure B1 Biases in sea surface temperature (SST, panel a) and salinity (SSS, panel b). Both model and observational data (EN4 version 4.2.0; Good et al., 2013) are averaged for 1960–1999.

Figure B2 Zonal-mean biases of seawater temperature (a–c) and salinity (d–f) with respect to observations (EN4 version 4.2.0; Good et al., 2013) for the Atlantic Ocean (a, d), Pacific Ocean (b, e) and Indian Ocean (c, f).

Table B1 Summary of the spatial correlation coefficient r and normalised root mean square error (NRMSE) between model data and observations from EN4 (version 4.2.0; Good et al., 2013).

Figure B3 Atlantic Meridional Overturning Circulation (AMOC) stream function (Sv).

Figure B4 1990–2009 mean vertical velocity (m s⁻¹) in the model at 1020 m depth (a) and 2920 m depth (b).

Figure B5 The mean of the annual maximum of the monthly mixed layer depth (m) for the period 1970–1999 in the model. The mixed layer depth is defined as the depth at which a 0.03 kg m⁻³ change in potential density with respect to the surface has occurred. Contour intervals are 50 for 0–500 and 500 for 500–3000.

Figure B6 The simulated zonal current (cm s⁻¹) at 960 m depth in the equatorial Pacific (averaged over January 2003–August 2009). Positive values indicate eastward flow.

Figure B7 CFC-12 concentration (pmol kg ⁻¹) in February 1998 along the A5 section in the Atlantic Ocean (see right panel) of the model (a) and of observations from GLODAPv1 database (panel b; Key et al., 2004). Contour intervals are 0.01, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.6, 0.8, 1.2 and 2 pmol kg ⁻¹.

Figure B8 (a–c) CFC-12 concentration (pmol kg ⁻¹) for the section A16 (a), P16 (b) and I8S9N (c). Panels (d)–(f) and (g)–(i) are as panels (a)–(c) but for the observed CFC-12 (GLODAPv1; Key et al., 2004) and for the difference between model and observation, respectively. The isolines in panels (a)–(f) are 0.01, 0.1, 0.4, 0.7, 1.0, 1.3, 1.6, 1.9 and 2.2 pmol kg⁻¹. The isoline increment in panels (g)–(i) is 0.2 pmol kg⁻¹.

Appendix C: Model–observation comparison of ocean biogeochemistry

C1 Net primary production, growth rate, biomass and limiting nutrients

The simulated net primary production, 48.7 Gt yr⁻¹ for bulk phytoplankton and 3 Gt yr⁻¹ for cyanobacteria, compares well with the satellite-based estimate of ∼52 Gt yr⁻¹ (Westberry et al., 2008; Silsbe et al., 2016).

The simulated growth rate μ (Fig. C1a and b, only shown for bulk phytoplankton because cyanobacteria has a much lower primary production) is broadly consistent with the large-scale patterns of the satellite-based μ estimates from Westberry et al. (2008) (Figs. C1c and C1d) and with field observations. In the central equatorial Pacific the simulated μ well reproduces the observed range (0.55–0.7 d⁻¹, Chavez et al., 1996; note the satellite-based estimates overestimate μ due to excluding iron limitation). In the subtropical gyres, the simulated μ (annual-mean 0.1–0.25 d⁻¹) is at the lower side of both the observations (annual mean 0.3–0.53 d⁻¹ in the North Pacific subtropical gyre, Letelier et al., 1996; annual mean 0.13–0.62 d⁻¹ in the North Atlantic subtropical gyre, Marañón, 2005) and the satellite-based μ estimates. In the Pacific sector of the Southern Ocean, the simulated μ (0.3–0.4 d⁻¹) in the austral summer is higher than the observations (about 0.1–0.2 d⁻¹; Boyd et al., 2000) and the satellite-based estimates. The simulated phytoplankton biomass is too high in the equatorial Pacific (>100 mg C m⁻³) and the Southern Ocean (>50 mg C m⁻³; Fig. C2) compared to the satellite-based estimates (<30 mg C m⁻³ for both regions; Westberry et al., 2008).

Figure C1 The 1999–2004 climatological-mean surface phytoplankton growth rates (d⁻¹) of the model (a, b, for bulk phytoplankton) and of the satellite-based estimates from Westberry et al. (2008) (c, d) for the boreal summer (a, c) and winter (b, d)​​​​​​​. The growth rate is identical between Hist_Popp and Hist_Laws.

Figure C2 The 1999–2004 averaged annual-mean surface phytoplankton biomass (mg C m⁻³) of the model.

Figure C3 Limiting nutrients for primary production in the model.

C2 Additional model–observation comparison for oceanic biogeochemical variables

The model captures the major features of the observed phosphate, DIC, oxygen and nitrate distribution. The biases of the above four variables are shown in Figs. 9b, 10g–i, C4, C5 and C6. We slightly underestimate the global mean phosphate by 0.2 mmol m⁻³, DIC by 41.3 mmol m⁻³, oxygen by 15 mmol m⁻³ and nitrate by 4.7 mmol m⁻³.

Figure C4 (a) DIC biases with respect to observation (GLODAPv1; Key et al., 2004) at the sea surface. (b–d) Zonal-mean DIC biases for the Atlantic, Pacific and Indian Ocean, respectively. Model data are averaged for 1990–1999.

Figure C5 As Fig. C4 but for simulated oxygen and observation from WOA13 (Garcia et al., 2013b).

Figure C6 As Fig. C4 but for simulated nitrate and observation from WOA13 (Garcia et al., 2013a).

Figure C7 The change rate of biological fractionation ϵ_p from 1960 to 2009.

Figure C8 The biological component δ¹³C${}_{\mathrm{DIC}}^{\mathrm{bio}}$ at the ocean surface for the model Hist_Popp (a) and observation (b). Panels (c–d) are as panels (a–b) but for the residual component δ¹³C${}_{\mathrm{DIC}}^{\mathrm{resi}}$.

Appendix D: The regenerated component of δ¹³C_DIC

The regenerated component of δ¹³C_DIC, δ¹³C^reg, relates to organic matter remineralisation and calcium carbonate dissolution. We neglect the dissolution of CaCO₃ following Sonnerup et al. (1999), who argued that this simplification only results in a small offset (<2 %). δ¹³C^reg is calculated as

$$\begin{array}{}\text{(D1)}& {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}^{\mathrm{reg}}={{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{DIC}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}^{\mathrm{pref}},\end{array}$$

with δ¹³C^pref given in Eq. (15). Note that the calculation of δ¹³C^pref in Eq. (15) only applies below the 200 m, which is roughly the euphotic zone depth (Eide et al., 2017a).

The temporal change of the regenerated component ${{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}}^{\mathrm{reg}}={{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1990}\mathrm{s}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{PI}}^{\mathrm{reg}}$ (Fig. D1a–c) generally shows a much smaller magnitude than δ¹³C${}_{\mathrm{SE}}^{\mathrm{pref}}$ (Fig. 12d–f). Above 1500 m, the δ¹³C${}_{\mathrm{SE}}^{\mathrm{reg}}$ is mainly caused by the change in remineralisation, as is illustrated by the change in AOU (Fig. D1d–f). Below 1500 m, the ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{SE}}^{\mathrm{reg}}$ is generally negative because δ¹³C_POC decreases by 2.2 ‰ from the pre-industrial period to the 1990s, mainly due to the decline of the biological fractionation factor ϵ_p under increasing surface CO₂(aq) (Fig. C7a).

Figure D1 The simulated change in the regenerated component ${{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}}^{\mathrm{reg}}={{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1990}\mathrm{s}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{PI}}^{\mathrm{reg}}$ for vertical sections A16 in the Atlantic Ocean (a), P16 in the Pacific Ocean (b) and I8S9N in the Indian Ocean (c). The locations of the vertical sections are shown in Fig. 12. Panels (d–f) are as panels (a–c) but for the change in AOU from pre-industrial times to the 1990s.

Appendix E: Applying the Eide et al. (2017a) approach to the model data

E1 Description of the Eide et al. (2017a) approach

To derive the global oceanic ¹³C Suess effect, Eide et al. (2017a) (hereafter E17) first applied the two-stage back-calculation method developed by Olsen and Ninnemann (2010) to calculate the ¹³C Suess effect using data from the World Ocean Circulation Experiment sections. The steps and assumptions of this stage are explained below. Next E17 mapped these ¹³C Suess effect estimates onto a 1 × 1^∘ grid with 24 vertical layers and obtained the three-dimension distribution of the ¹³C Suess effect in the global ocean. For simplicity, hereafter the above procedure is collectively referred to as E17's approach.

E17 first assume that any oceanic CFC-12 signal before 1940 is negligible and the oceanic ¹³C Suess effect at any time t after 1940, δ¹³C_SE(t−1940)​​​​​​​, is proportional to CFC-12 partial pressure at time t:

$$\begin{array}{}\text{(E1)}& {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}(\mathrm{t}-\mathrm{1940})}\sim a\cdot {\text{pCFC-12}}_{\mathrm{t}}.\end{array}$$

Here the proportionality factor a is time-invariant. δ¹³C_DIC at any time t after year 1940 is decomposed as

$$\begin{array}{}\text{(E2)}& {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{t}}={{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}(\mathrm{t}-\mathrm{1940})}+{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{pref}}+{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}}.\end{array}$$

The calculation of δ¹³C^pref is given in Eq. (15) and δ¹³C^reg in Eq. (D1). E17 include two additional terms on the right-hand side of the above equation Δδ¹³C^reg and Δδ¹³C^pref (see their Eq. 4), which represent any changes not related to the ¹³C Suess effect, e.g. changes in ocean carbon cycle. We do not explicitly write these two terms as they are set to zero by E17.

Decomposing the left-hand side of Eq. (E2) into a preformed component and a regenerated component gives

$$\begin{array}{}\text{(E3)}& \begin{array}{rl}{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{t}}^{\mathrm{pref}}& ={{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}(\mathrm{t}-\mathrm{1940})}+{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{pref}}\\ & -({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{t}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}}).\end{array}\end{array}$$

Following Gruber et al. (1996), E17 assume a steady-state ocean over the period of interest and set $({{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{t}}^{\mathrm{reg}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{reg}})$ to zero, and this gives

$$\begin{array}{}\text{(E4)}& {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{t}}^{\mathrm{pref}}={{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}(\mathrm{t}-\mathrm{1940})}+{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{1940}}^{\mathrm{pref}}.\end{array}$$

Combining Eqs. (E1) and (E4) yields linear relationship between δ¹³C${}_{\mathrm{t}}^{\mathrm{pref}}$ and pCFC-12_t:

$$\begin{array}{}\text{(E5)}& {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{t}}^{\mathrm{pref}}\sim a\cdot {\text{pCFC-12}}_{\mathrm{t}}+b,\end{array}$$

where b contains term δ¹³C${}_{\mathrm{1940}}^{\mathrm{pref}}$. Thus, the proportionality factor a can be determined with δ¹³C${}_{\mathrm{t}}^{\mathrm{pref}}$ and pCFC-12_t at time t, and δ¹³C_SE(t−1940) can be obtained with Eq. (E1).

To scale δ¹³C_SE(t−1940) to δ¹³C_SE(t−PI) for the full industrial period, the assumption is used that the oceanic δ¹³C_DIC change scales with the atmospheric δ¹³CO₂ change, i.e.:

$$\begin{array}{}\text{(E6)}& {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}(\mathrm{t}-\mathrm{PI})}=f_{\mathrm{atm}}\cdot {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}(\mathrm{t}-\mathrm{1940})}=f_{\mathrm{atm}}\cdot a\cdot {\text{pCFC-12}}_{\mathrm{t}},\end{array}$$

with

$$\begin{array}{}\text{(E7)}& f_{\mathrm{atm}}={\displaystyle \frac{{{\mathit{\delta }}^{\mathrm{13}}\mathrm{CO}}_{\mathrm{2},\mathrm{t}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{CO}}_{\mathrm{2},\mathrm{PI}}}{{{\mathit{\delta }}^{\mathrm{13}}\mathrm{CO}}_{\mathrm{2},\mathrm{t}}-{{\mathit{\delta }}^{\mathrm{13}}\mathrm{CO}}_{\mathrm{2},\mathrm{1940}}}}.\end{array}$$

E2 Calculation of SE_pref, the oceanic ¹³C Suess effect estimate using E17's approach and model data

To achieve a result comparable to E17, we select the model data at the geographic locations for which both CFC-12 and δ¹³C_DIC measurements are available. The observational data set of E17 has data from one cruise in the South Atlantic (A13.5) in 2010. We do not include this cruise data because the applied ERA20C forcing, and, thus, our simulations ends in 2009. Here we use the observations compiled by Schmittner et al. (2013) because δ¹³C_DIC in this data set has been quality controlled and is publicly available. Following E17, we use data at the model layers between 200 m and the simulated CFC-12 penetration depth (defined as pCFC-12 = 20 patm (pico-atmosphere); see the thick grey lines in Fig. 14). We take model data of year t=1994. By performing a linear regression (Eq. E5) for five ventilation regions (the North Atlantic, South Atlantic, North Pacific, South Pacific and Indian Ocean) we obtain the regression parameters, hereafter referred to as a_pref and b_pref. Applying Eq. (E6) to the three-dimension model data of pCFC-12 for t=1994, regression slope a_pref and f_atm=1.5 (determined with Eq. E7 for year 1994), we obtain the estimate of the global oceanic ¹³C Suess effect, SE_pref, in year 1994 (Eq. 16).

The regressional relationships between ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{1994}}^{\mathrm{pref}}$ and pCFC-12₁₉₉₄ and the regression coefficients a_pref and b_pref are shown in Fig. E2 (the water masses in this figure are defined in Table E1). The coefficient of determination r², the percentage of the variance in the data explained by the regressional relationship, ranges between 0.33 and 0.66. The strength of these linear relationships is acceptable considering the lowest r²=0.22 in E17.

The regression relationships between δ¹³C^pref and pCFC-12 in our model (Fig. E2) show some quantitative differences compared to those of E17 (see their Fig. 3). These differences originate from model biases in the distribution and properties of water masses. These mismatches do not affect the analysis and conclusions in Sect. 5. Nevertheless, we briefly discuss their causes for better understanding of the model behaviour.

First, the definitions of several water masses in the model are slightly different from those of E17 (comparing our Table E1 with their Table 2).

Second, our simulated δ¹³C${}_{\mathrm{t}}^{\mathrm{pref}}$ in the deep and bottom waters (Antarctic Bottom Water, Circumpolar Deep Water, Pacific Deep Water and Indian Deep Water) in the Southern Hemisphere (Fig. E2c and e and Fig. E3c) is higher than that in E17 (see their Fig. 3a, c and e). The possible reasons for this difference are related to mixing and primary production in the Southern Ocean. Here, the simulated deep convection, which primarily occurs in the open ocean rather than the along continental shelf, is too strong in the model. This can be seen by the large mixed layer depth (Fig. B5) and by the CFC-12 bias along selected vertical sections (Fig. B8), which feature persistent positive biases off the Antarctic continental shelf in the Atlantic, Pacific and Indian sectors of the Southern Ocean. Furthermore, the Southern Ocean has a primary production that is too high in the model (about a factor of 1.5 of the satellite-based net primary production estimates from Westberry et al., 2008). The high primary production causes higher surface δ¹³C_DIC than observations (see the South Pacific Ocean in Fig. 8c). Consequently, the simulated preformed component δ¹³C${}_{\mathrm{t}}^{\mathrm{pref}}$ in the bottom and deep water masses of the Southern Ocean is higher than observed values in E17.

Third, the lowest values of δ¹³C${}_{\mathrm{t}}^{\mathrm{pref}}$ (<1.4 ‰) are often found in the upwelling regions in the model. This is due to the upward transport of water from the ocean interior that has lower δ¹³C_DIC than observations (Fig. 10e and f).

Figure E1 The difference (SE_pref−SE_Mod) for the vertical sections A16 in the Atlantic Ocean (a), P16 in the Pacific Ocean (b) and I8S9N in the Indian Ocean (c). Panels (d–f) and (g–i) are as panels (a–c) but for (SE_total−SE_Mod) and (SE_pref−SE_total), respectively. The isoline increment is 0.05 ‰. The thick grey line is the pCFC-12₁₉₉₄=20 patm isoline, below which SE_pref is generally very small (<0.05 ‰).

E3 Linear regression for subregions in the Indian Ocean

We can span regressional relationships for the subtropical gyres of the Indian Ocean and North Pacific Ocean because we use all model levels between 200 m and the pCFC-12 = 20 patm isoline at a given geographical location, and therefore we have more data points than field measurements. In the Indian Ocean, performing linear regression for δ¹³C${}_{\mathrm{1994}}^{\mathrm{pref}}$ and pCFC-12₁₉₉₄ in the Subtropical Gyre Water and Sub-Antarctic Mode Water yields regression parameters $a_{\mathrm{pref}}^{\mathrm{STGW}}=-\mathrm{0.65}\times {\mathrm{10}}^{-\mathrm{3}}$, $b_{\mathrm{pref}}^{\mathrm{STGW}}=\mathrm{1.98}$ and r²=0.49. The more negative $a_{\mathrm{pref}}^{\mathrm{STGW}}$ compared to regression slope $a_{\mathrm{pref}}=-\mathrm{0.47}\times {\mathrm{10}}^{-\mathrm{3}}$ obtained for the whole Indian Ocean suggests an underestimation of the ¹³C Suess effect. The mean pCFC-12 in the Indian subtropical region at 200 m ${\text{pCFC-12}}_{\mathrm{1994}}^{\mathrm{STGW}}=\mathrm{440}$ patm. Following Eq. (E6), we can calculate the mean underestimation for the subtropical Indian Ocean as $f_{\mathrm{atm}}\cdot (a_{\mathrm{pref}}-a_{\mathrm{pref}}^{\mathrm{STGW}})\cdot {\text{pCFC-12}}_{\mathrm{1994}}^{\mathrm{STGW}}=\mathrm{0.12}$ ‰.

E4 Calculation of SE_total

To calculate SE_total we perform a linear regression for the total oceanic ¹³C Suess effect δ¹³C_{SE(1994–1940)} and pCFC-12₁₉₉₄:

$$\begin{array}{}\text{(E8)}& {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}\left(\text{1994–1940}\right)}\sim a_{\mathrm{total}}\cdot {\text{pCFC-12}}_{\mathrm{1994}}+b_{\mathrm{total}}.\end{array}$$

Here the model data are subsampled in the same manner as in Sect. E2. Next, applying a correction for the period prior to 1940 (in analogy to Eq. E6) we obtain the expression of SE_total in Eq. (19).

The regression relationships in Eq. (E8) and regression coefficients are given in Fig. E3. For the Indian, North Pacific, North Atlantic and South Atlantic Ocean, r² lies between 0.34 and 0.67, which suggests an acceptable strength of the relationships. In the South Pacific Ocean we find low r²=0.07. This low r² is a result of the high variability in the change in the regenerated component (Fig. 14h) which corrupts the regression. Therefore we omit the South Pacific in the calculation of SE_total.

Table E1 Water masses and their definitions in the model.

^a Water masses are combined together rather than separately defined as in Eide et al. (2017a).
^b A different σ_θ threshold is used here compared to Eide et al. (2017a).

Figure E2 Regressional relationships δ¹³C${}_{\mathrm{1994}}^{\mathrm{pref}}\sim a_{\mathrm{pref}}\cdot {\text{pCFC-12}}_{\mathrm{1994}}+b_{\mathrm{pref}}$ for the Indian Ocean (a), the North Pacific (b), the South Pacific (c), the North Atlantic (d) and the South Atlantic (e). Different colours and symbols indicate different water masses. The full names, as well as the definitions, of the water masses are listed in Table E1. The regression slopes a_pref are used to calculate SE_pref in Eq. (16). In the Indian Ocean the regression relationship for the Subtropical Gyre Water and Sub-Antarctic Mode Water (red upward triangle in panel a) is $y=-\mathrm{0.65}\times {\mathrm{10}}^{-\mathrm{3}}x+\mathrm{1.98},r^{\mathrm{2}}=\mathrm{0.49}$. In the North Pacific the regression relationship for the Subtropical Gyre Water (red upward triangle in panel b) is $y=-\mathrm{0.44}\times {\mathrm{10}}^{-\mathrm{3}}x+\mathrm{1.66},r^{\mathrm{2}}=\mathrm{0.26}$.

Figure E3 As Fig. E2 but for the regression relationships ${{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}_{\mathrm{SE}(\mathrm{1994}-\mathrm{1940})}\sim a_{\mathrm{total}}\cdot {\text{pCFC-12}}_{\mathrm{1994}}+b_{\mathrm{total}}$. The regression coefficients a_total and b_total are used to calculate SE_total following Eq. (19).