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

Direct comparison between paleo oceanic δC records and model results facilitates assessing simulated distributions and properties of water masses in the past. To accomplish this, 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 existing oceanic carbon pools. We account for 5 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: Popp p varies with surface dissolved CO2 concentration (Popp et al., 1989), while Laws p additionally depends on local phytoplankton growth rates (Laws et al., 1995). When compared to observations of δC in dissolved inorganic carbon (DIC), both parameterisations yield similar performance. However, with regard to δC in particulate organic carbon Popp p shows a considerably improved performance 10 than Laws p , because the latter results in a too strong preference for C. The model also well reproduces the oceanic C Suess effect, i.e. the intrusion of the isotopically light anthropogenic CO2 into the ocean, based on comparison to other existing C models and to observation-based oceanic carbon uptake estimates over the industrial period. We further apply the approach of Eide et al. (2017a), who derived the first global oceanic C Suess effect estimate based on observations, to our model data that has ample spatial and temporal coverage. With this we are able to analyse in detail the 15 underestimation of C Suess effect by this approach as it has been noted by Eide et al. (2017a). Based on our model we find underestimations of 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 the underestimation to two assumptions in Eide et al. (2017a)’s approach: a spatially-constant preformed component of δCDIC in year 1940 and neglecting C Suess effect in CFC-12 free water. 20 1 https://doi.org/10.5194/bg-2021-32 Preprint. Discussion started: 12 February 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
The stable carbon isotopic composition δ 13 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 13 C and 12 C: where R std is an arbitrary standard ratio. In observational studies, the ratio 13 C/ 12 C in Pee Dee Belemnite (PDB; Craig, 1957) is conventionally used for R std .
δ 13 C provides information on past changes of water mass distribution and properties (e.g. Curry and Oppo, 2005;Peterson et al., 2014). Direct comparison between paleo δ 13 C measurements and simulated δ 13 C facilitates evaluating the ability of simulations of a complete last glacial cycle within the German climate modeling initiative PalMod (Latif et al., 2016). Before applying the new 13 C module to paleo simulations, we evaluate it by comparison to observational data in the present day ocean.

35
Earlier versions of HAMOCC already featured a 13 C module. HAMOCC3 included prognostic 13 C variables for dissolved inorganic carbon (DIC), particulate organic matter and calcium carbonate (Maier-Reimer, 1993). HAMOCC3 also accounted for temperature-dependent isotopic fractionation during air-sea gas exchange (higher δ 13 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 40 was parameterised by a spatially and temporally uniform factor. This approach for biological fractionation of 13 C, however, could not reproduce the observed large meridional gradient of δ 13 C in particulate organic matter (Goericke and Fry, 1994).
Since then, HAMOCC was refined in particular with regard to its representation of plankton dynamics, which currently 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 13 C module that considers the refined ecosystem representa-45 tion and test different non-uniform parameterisations for biological fractionation during phytoplankton growth.
We choose two parameterisations for biological fractionation that suit the complexity of our model and were successfully applied in previous modelling studies (Hofmann et al., 2000;Schmittner et al., 2013;Tagliabue and Bopp, 2008;Jahn et al., 2015;Dentith et al., 2020): 1) the parameterisation of Popp et al. (1989), which empirically relates biological fractionation to the concentration of dissolved CO 2 in seawater; 2) the parameterisation of Laws et al. (1995), which considers dissolved CO 2 50 concentration and phytoplankton growth rate. We omit more complex parameterisations that include effects of cell membrane permeability of molecular CO 2 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. To assess the model's performance, we run pre-industrial and present-day HAMOCC6 version in Mauritsen et al. (2019), we allow DOM degradation in low oxygen conditions until all available O 2 is consumed.

The stable carbon isotope 13 C in HAMOCC6
HAMOCC6 simulates total carbon C, which is the sum of the three natural isotopes 12 C, 13 C and 14 C. Because in nature 12 C constitutes about 98.9% of the total carbon and 13 C only constitutes about 1.1 % (Lide, 2002), in HAMOCC we assume 100 12 C = C. We include a 13 C counterpart for each 12 C prognostic variable, that is, we introduce seven new tracers for the water column and three for the sediment. 13 C only mimics the 12 C biogeochemical fluxes, modified by the corresponding isotopic fractionation. We assume 13 C inventory to be as large as the inventory of 12 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 13 C fractionation during air-sea exchange and carbon uptake by bulk phytoplankton and by cyanobacteria. Because the isotopic 105 fractionation during the production of calcium carbonate is small (Turner, 1982), it is not considered in this study.

Fractionation during air-sea gas exchange
13 C isotopic fractionation during air-sea gas exchange is temperature-dependent. We adopt the calculation of 13 C air-sea gas exchange recommended by the OMIP protocol (Orr et al., 2017).
The net air-sea CO 2 gas exchange flux F reads Here, pCO 2 surf and pCO 2 atm are the partial pressures of CO 2 in the surface seawater and in the atmosphere, respectively. The piston velocity k CO2 (ms −1 ) for CO 2 and the solubility γ CO2 (mol L −1 atm −1 ) of CO 2 are calculated following Wanninkhof (2014) and Weiss (1974), respectively.
in which, R g and R atm are the ratios of 13 C/ 12 C in surface pCO 2 and in atmospheric CO 2 , respectively. Following Zhang et al. (1995), we can re-write Eq.
Here, T C is the seawater temperature in • C and f CO3 = CO 2− 3 /DIC is the fraction of carbonate ions in DIC. Because in Eq. (6) the temperature dependency is weak, we use a constant aq←g = −1.24, obtained at T C = 15 • C in the model, following Schmittner et al. (2013). In Eq. (7) we neglect the first term 0.014 T C f CO3 , because f CO3 is generally smaller than 0.1 and because the constant factor is one order of magnitude smaller than that of the second term 0.105 T C .

Fractionation during phytoplankton growth
The lighter stable carbon isotope 12 C is preferentially utilised than 13 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 Here G (µmol C L −1 day −1 ) 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 2 (aq) and by the biological fractionation factor p = (α Phy←aq − 1) × 10 3 related to the fixation of CO 2 (aq). Here the subscript "Phy" 140 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. Popp p = −17 log(CO 2 (aq)) + 3.4, Here, CO 2 (aq) (µmol L −1 ) is aqueous CO 2 in surface water, µ (day −1 ) is the specific growth rate of bulk phytoplankton or 145 of cyanobacteria. Note that Laws et al. (1995) measured aq←Phy . Because α Phy←aq is close to unity, p ≈ − aq←Phy (Zeeb and Wolf-Gladrow, 2001). In Eq. (11), we set the seawater density ρ sea a constant value of 1.025 kg L −1 . Then Eq. (11) is simplified to Both CO 2 (aq) and µ (depending on local conditions of light, water temperature and nutrient availability) are determined in 150 HAMOCC. Figure 1 illustrates the values of Popp p and Laws p under typical ranges of CO 2 (aq) and µ in the ocean. When µ ≤ 1,

Setup
We conduct ocean-only simulations using the MPIOM-1.6.3p1 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 160 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 are represented by nine levels. The time step is 1 hour. In this set-up, 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 hydro-dynamical processes within the ocean identical to e.g.

Experimental design
We carry out pre-industrial spin-up simulations followed by historical (1850-2010) simulations. We force the model with the sea-surface boundary conditions from ERA20C (Poli et al., 2016), which covers the period 1901-2010. For the pre-industrial period, we cyclically apply the forcing of 1905-1929 and set the atmospheric CO 2 mixing ratio to 280 ppmv. We first conduct a spin-up run without 13 C tracers until the long-term averaged global net air-sea 13 CO 2 flux is smaller than 0.05 Pg C yr −1 (ade-170 quate to the C4MIP criterion for steady state conditions of <0.1 Pg C yr −1 ; Jones et al., 2016). This model state is the starting point for the two spin-up runs including 13 C tracers, PI_Popp and PI_Laws, which are based on the biological fractionation parametrisation Popp p (Eq. 10) and Laws p (Eq. 12), respectively.
The 13 C tracers are initialised as follows. The mean δ 13 C of the marine organic matter is about −20‰ (Degens et al., 1968).
Therefore, we set the initial concentrations of 13 C in the bulk phytoplankton, cyanobacteria, zooplankton, dissolved organic 175 carbon, particulate organic carbon in the water column and particulate organic carbon in the sediment to 0.98 (according to Eq. 1) of their 12 C counterparts. The initial 13 C DIC in the water column is calculated following the relation between δ 13 C DIC and PO 4 (Lynch- Stieglitz et al., 1995): The initial concentrations of 13 C CaCO3 in the water column and in the sediment, and the initial concentration of 13 C DIC in pore 180 water are set identical to their 12 C counterparts.  The pre-industrial stable carbon isotope ratio δ 13 CO 2 of atmospheric CO 2 is fixed at −6.5‰. The input rate of dissolved organic 13 C (DO 13 C) is calculated as the product of the input rate of DOC and the sea-surface DO 13 C/DOC ratio; the input rate of 13 CO 2− 3 is the product of the input rate of CO 2− 3 and the sea-surface 13 CO 2− 3 /CO 2− 3 ratio. This approach to determine 13 C input rates results in a small drift in the water-column 13 C inventory but it only has minor impact on the simulation results (see 185 Appendix A). PI_Popp and PI_Laws are spun up for 2500 simulation years such that 13 C inventory adjusts to be consistent with the simulated biogeochemical and hydrodynamical processes. Equilibrium states are reached with 98% of the ocean volume having a drift of less than 0.001‰ year −1 (employing the same criteria as for 14 C in OMIP protocol, Orr et al., 2017).
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 rates of DO 13 C, DOC, 13 CO 2− 3 , CO 2− 3 and SiO 4 are kept constant, and are the same as those in the pre-industrial simulations.

195
3 Model results and observations in the late 20th century In this section, we compare simulated 13 C between the two simulations Hist_Popp and Hist_Laws and evaluate the two experiments by comparison to observed δ 13 C POC and δ 13 C DIC . The observations used here are the surface δ 13 C POC measurements assembled by Goericke and Fry (1994) and the observed δ 13 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 observational data sets horizontally onto a 200 1x1 degree grid and vertically 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 δ 13 C POC and δ 13 C DIC over a 1x1 de- gree grid. To quantitatively compare the performance between Hist_Popp and Hist_Laws and to other 13 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 δ 13 C POC or δ 13 C DIC during the observational periods) between model 205 results and observation.

Isotopic signature of particular organic carbon in the surface ocean
For comparison between Hist_Popp and Hist_Laws, the climatological mean state of δ 13 C POC is derived by averaging over 1960-1991, the period when most δ 13 C POC measurements were collected. In Hist_Popp, the climatological annual-mean surface δ 13 C POC has a global mean value of −22.5‰ and it shows a distict horizontal pattern (Fig. 3a). Less negative values up 210 to −19.3‰ are found in the subtropical regions, where alkalinity is typically high and CO 2 (aq) is consequently low. This low CO 2 (aq) results in a smaller isotope fractionation during carbon fixation by phytoplankton (Eq. 10, Fig. 1) with a biological fractionation factor p > −13‰ (Fig. 3c). Poleward of the subtropical regions, δ 13 C POC gradually decreases. The reason for this is twofold. First, p decreases from −13 to about −20‰ following the increase of CO 2 (aq). Second, the thermal effect of equilibrium fractionation causes about 3‰ more fractionation in the polar regions than in the tropical and subtropical 215 regions (according to Eqs. 7 and 9). The lowest δ 13 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 δ 13 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 2 (aq).

220
Compared to Hist_Popp, Hist_Laws shows lower annual-mean surface δ 13 C POC (Fig. 3b), with a global-mean value of −29.9‰ due to more negative p (Fig. 3d). This is because Laws p ( Fig. 1) is always more negative than Popp p when the simulated mean growth rates (Figs. B1a and B1b) are lower than 1 day −1 . As Laws p increases with growth rate (Eq. 12), we find less negative δ 13 C POC (up to −24.1‰) in the central tropical Pacific, where highest growth rates are simulated (Figs. B1a and B1b).
The lowest δ 13 C POC of −33‰ occurs in the Arctic Ocean and around Antarctica due to the combination of low growth rate, high 225 CO 2 (aq) and low seawater temperature. The meridional range of the annual-mean δ 13 C POC in Hist_Laws (∼ 9‰) is smaller than that of Hist_Popp (∼ 11‰) because for low growth rates Laws p is generally less sensitive to CO 2 (aq) changes compared to Popp p ( Fig. 1). This also results in a smaller annual range of δ 13 C POC in high latitudes (Fig. 3f) than Hist_Popp (Fig. 3e). In the low and mid latitudes, Hist_Laws show larger annual range of δ 13 C POC because in these regions CO 2 (aq) concentrations are relatively stable but growth rates shows noticeable seasonal variability.

230
Hist_Popp captures major features of the observed δ 13 C POC (Figs. 4a, 4c and 4e). 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 δ 13 C POC than the observations (a global mean bias of −8‰) and smaller δ 13 C POC difference between low and high latitudes (Figs. 4b,4d and 4f). This is also seen in a recent study by Dentith et al. (2020), who tested Popp  to gain a better performance. Around 60 • S of the Atlantic Ocean ( Fig. 4b), Hist_Laws simulates a smaller range of δ 13 C POC than the observations. This is also a result of the small δ 13 C POC annual range produced by Laws p (Fig. 3f). Between 40 • S and 40 • N in the Atlantic Ocean, Hist_Laws simulates δ 13 C POC peaks in the region of high growth rates south of the Equator, whereas the observed high δ 13 C POC values locate between the Equator and 20 • N.

240
In the Indian Ocean around 45 • S, Hist_Popp does not capture the prominent δ 13 C POC peak in the field data ( Fig. 4e), although the simulated CO 2 (aq), the controlling factor in the parameterisation Popp p (Eq. 10), well reproduces the meridional variation of the observations (Fig. 4g). This is because although the empirical correlation between p and CO 2 (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 δ 13 C POC  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 δ 13 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 δ 13 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 δ 13 C POC compares well to that of the FAMOUS model ( 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.
3.2 Isotopic signature of dissolved inorganic carbon δ 13 C DIC Figures 5a, 5b and 6a -6f compare the climatological annual mean of δ 13 C DIC (averaged over 1990 -2005, when most δ 13 C DIC 260 measurements were collected) between Hist_Popp and Hist_Laws. The two simulations exhibit very similar δ 13 C DIC patterns for both surface and interior ocean. The surface seawater DIC is enriched in 13 C due to the preferential uptake of the light isotope 12 C by phytoplankton during primary production. As particulate organic matter sinks and is remineralised at depth, the negative δ 13 C POC signal is released. Consequently, in both Hist_Popp and Hist_Laws, δ 13 C DIC at the surface is generally Hist_Laws), which is expected as they are run using the same prescribed atmospheric δ 13 CO 2 (Schmittner et al., 2013).
Given very similar mean surface DI 13 C, the larger vertical DI 13 C gradients in Hist_Laws, established by more negative δ 13 C POC (Figs 3a and 3b), yields lower DI 13 C concentration at depth. This adjustment of DI 13 C content in the ocean interior takes place during the pre-industrial spin-up phase of the simulations via air-sea 13 CO 2 exchange (Appendix A). At the end of the 2500-year spin-up, the water-column DI 13 C inventory in PI_Laws is 1.1 × 10 12 kmol lower than PI_Popp, yielding 275 a global mean δ 13 C DIC difference of 0.25‰ (Figs. 6g -6i). Such interior-ocean δ 13 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 δ 13 C DIC in Hist_Laws leads to lower annual-mean surface δ 13 C DIC and larger δ 13 C DIC annual range in regions of upwelling (Figs. 5c and 5d).
When compared to the observed δ 13 C DIC , Hist_Popp (r = 0.81, NRMSE = 0.7) has a slightly better performance than 280 Hist_Laws (r = 0.80, NRMSE = 1.1). Hist_Laws generally shows too strong vertical gradients of δ 13 C DIC and too low δ 13 C DIC values in the ocean interior, as is seen in the depth profiles of horizontally-averaged δ 13 C DIC (Fig. 7). This points to too strong preference for the isotopically light carbon simulated by Laws p as is already discussed in Section 3.1. Given the slightly better performance of Hist_Popp than Hist_Laws regarding δ 13 C DIC , we focus in the following on the comparison between Hist_Popp and observed δ 13 C DIC .

285
Figures 8 and 9 contain model-observation comparison for the surface and interior-ocean δ 13 C DIC , respectively. Overall, the magnitude and spatial distribution of the observed δ 13 C DIC is well-captured by Hist_Popp. In the surface ocean, the mean δ 13 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 290 decompose δ 13 C DIC into a biological component δ 13 C bio DIC and a residual component δ 13 C resi DIC , driven by air-sea exchange and ocean circulation: 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 ∆ photo = −19‰ (Eide et al., 2017b) and R C:P = 122 (Takahashi 295 et al., 1985). To calculate δ 13 C bio DIC from observations, we employ δ 13 C DIC | M.O. = 0.5‰, DIC M.O. = 2200 µmol kg −1 (Eide et al., 2017b), and PO 4 from the World Ocean Atlas (WOA13; Garcia et al., 2013a). Considering the strong seasonality in PO 4 Between 30 • N and 30 • S in the surface ocean, the simulated δ 13 C bio DIC is generally lower than the observation-based δ 13 C bio DIC 305 with a mean negative bias of about −0.1‰ (Fig. 8d). This is caused by the underestimation of primary production in the subtropical gyres (due to the underestimation of phytoplankton growth rates, see Appendix B) and the consequently reduced enrichment of 13 C in surface DIC. A strong positive δ 13 C bio DIC bias of 0.6 to 1‰ is seen in the North Pacific, where in the model iron is not a limiting nutrient, in contrast to observations (Moore et al., 2013). In the equatorial central Pacific, a weak positive δ 13 C bio DIC bias < 0.2‰ is caused by a too high primary production. Specifically, the simulated phytoplankton growth rates in In the Southern Ocean, a strong positive δ 13 C bio DIC bias of 0.6 to 1‰ (Fig. 8d) results from a too large nutrient supply from 315 the interior ocean to the surface. The cause for this too large nutrient supply is two fold. First, organic matter is remineralised at too shallow depths in HAMOCC, as is shown by the positive apparent oxygen utilisation (AOU) biases above 500 m south of 45 • S (Figs. 9j -9l). Second, MPIOM simulates a too large upward transport due to too strong upwelling. In particular, below 1000 m, the simulated upward velocity shows noticeably larger magnitude (> 5 × 10 −6 m s −1 , Fig. D1) than that of a dynamically consistent and data-constrained ocean state estimate (see Figure 1 in Liang et al., 2017). Our model also features We find strong δ 13 C resi DIC negative biases of −0.5 to −1‰ (Fig. 8e) in the North Pacific and the Southern Ocean, which partially compensate the positive biases of δ 13 C bio DIC (Fig. 8d) in these regions. One major cause for the negative δ 13 C resi DIC bias in these two regions is our model overestimating the uptake of anthropogenic carbon, as is illustrated by the net air-sea CO 2 325 difference between the model and the observation (Fig. 8f). Consequently, the decreased atmospheric 13 C/ 12 C ratio over the industrial period further lowers δ 13 C DIC in the two ocean regions in the model. In the Southern Ocean, a too large upward transport of 13 C-depleted water at depth to the surface also contributes to a negative δ 13 C resi DIC bias. In the interior ocean, δ 13 C DIC is controlled by remineralisation of 13 C-depleted organic matter and by ocean circulation (Broecker and Peng, 1993;Lynch-Stieglitz et al., 1995;Schmittner et al., 2013). Low δ 13 C DIC is often found in waters of 330 high nutrient concentration and vice versa. Thus, we find positive (negative) δ 13 C DIC biases coincide with negative (positive) phosphate biases (Figs. 9d -9i). In the Atlantic Ocean between 1000 and 3000m, the North Pacific above 1500 m and the Indian Ocean below 1000 m, positive δ 13 C DIC biases and negative phosphate biases are mainly caused by a too low remineralisation, as is shown by the negative AOU biases (Figs 9j -9l). North of 30 • S in the Atlantic Ocean, the negative δ 13 C DIC biases below 3000 m, together with the negative δ 13 C DIC biases between 1000 and 3000 m, suggest too strong δ 13 C DIC vertical gradients 335 in the model (Fig. 9d). This results from a too shallow lower boundary of the NADW cell, constantly located above 2800 m, This is can be seen from the CFC-12 distribution along the zonal Section A5 at 24 • N (Fig. D3). The observed deeper CFC-12 340 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 δ 13 C DIC bias in the deep eastern equatorial Pacific (Fig. 9e). The cause is the 'nutrient trapping' problem in the model, characterised by too high nutrient concentrations in the deep eastern equatorial Pacific (Fig. 9h), 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 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 the Equatorial Intermediate Current System) shows too little spatial variability and too low speeds of ∼ 0.2 cm s −1 (Fig. D4), compared to the observed alternating 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 HAMOCC does not occur in the simulations of Schmittner et al. (2013). Our model shows noticeable better performance than that of Dentith et al. (2020). The latter study simulates too high δ 13 C DIC over all depth levels, which the authors ascribe to underlying biases in the biological carbon cycle.
4 Oceanic 13 C Suess effect The oceanic δ 13 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 δ 13 C (Fig. 2) affects oceanic δ 13 C via air-sea gas exchange, leading to a generally decrease of δ 13 C DIC . The distribution of this δ 13 C DIC change, i.e. the oceanic 13 C Suess effect, could serve as benchmark for ocean models to evaluate the uptake and re-distribution of the anthropogenic 365 CO 2 emissions in the ocean.
The model is able to reproduce the size of the global oceanic anthropogenic CO 2 sink. 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 δ 13 C Suess effect, δ 13 C SE , as the difference between the 1990s-averaged δ 13 C DIC from Hist_Popp and the pre-industrial cli-370 matological (50-year mean) δ 13 C DIC from PI_Popp. δ 13 C SE calculated using the results of Hist_Laws and PI_Laws only shows marginal difference (global-mean < 0.04‰), and is therefore not presented. The surface mean δ 13 C SE in this study is −0.  (Figs. 10d -10f). The 375 strongest oceanic 13 C Suess effect is found in the subtropical gyres, where water masses have long residence times at the ocean surface and therefore receive a strong anthropogenic imprint (Quay et al., 2003). At the surface of the subtropical gyres, our simulated δ 13 C SE generally varies between −0.8 and −1.1‰, which compares well to the the surface ocean δ 13 C decrease of −0.9 ± 0.1‰ recorded by coral and sclerosponges (Wörheide, 1998;Böhm et al., 1996Böhm et al., , 2000Swart et al., 2002Swart et al., , 2010 and to the estimates of −1.0 ± 0.09‰ extracted from GLODAPv2 (Olsen et al., 2016;Eide et al., 2017a). In the subtropical gyres of 380 the South Atlantic, the Pacific and the Indian Ocean, δ 13 C SE is mainly confined to upper 1000 m depth. In the North Atlantic, δ 13 C SE penetrates deeper than the other ocean regions, due to the intensive ventilation related to the formation of NADW. One noticeable discrepancy between the simulated δ 13 C SE and the estimates of Eide et al. (2017a) is some local positive δ 13 C SE values occur in our model. This difference will be discussed in Section 4.2.  Olsen and Ninnemann (2010) to calculate 13 C Suess effect using data from the World Ocean Circulation Experiment sections. Next they mapped these 13 C Suess effect estimates over a 1x1 degree grid with 24 vertical layers and obtained the three-dimension distribution of 13 C Suess effect in the global ocean. For simplicity, hereafter the above procedure is collectively referred to as E17's approach. E17 have noted their outcome is likely to underestimate the oceanic 13 C Suess effect by 0.15 to 0.24‰ 200 m, globally. However, they can not provide a quantitative explanation for the sources of the underestimation.
As our model reasonably reproduces the anthropogenic CO 2 uptake and δ 13 C SE distribution in the ocean and it includes all necessary variables (such as DI 13 C, CFC-12) required in E17's approach, we are able to investigate such potential underestimation by applying E17's approach to our model data. Specifically, we aim to extract information on the spatial distribution 395 of the potential underestimation and to quantitatively explain the causes for the underestimation. Below we briefly present the key assumptions and equations of E17's approach. A detailed description for E17's approach is given in Appendix C.
E17 assumed the oceanic 13 C Suess effect at any time t after 1940 is proportional to CFC-12 partial pressure (pCFC-12): Here the proportionality factor a is time-invariant. By decomposing δ 13 C DIC into a preformed component δ 13 C pref arising from 400 the transport of the surface water with specific DIC and DI 13 C and a regenerated component δ 13 C reg due to organic matter remineralization and calcium carbonate dissolution, the following equation is derived: The calculation of δ 13 C pref t and δ 13 C reg t are detailed in Appendix C. Then E17 assumed the regenerated component is constant in time, i.e. δ 13 C reg t = δ 13 C reg 1940 , which gives: Combining Eq. (15) and Eq. (17), together with regarding the preformed component for year 1940, δ 13 C pref 1940 , as a term independent of pCFC-12, yields a linear relationship between the preformed component δ 13 C pref t and pCFC-12 t : The regression coefficients a and b are determined with δ 13 C pref t and pCFC-12 t from measurement deeper than 200 m depth (be-410 low which the approach applies). With a and observed pCFC-12 t , δ 13 C SE(t−1940) on the ocean observation sections is obtained using Eq. (15). To scale δ 13 C SE(t−1940) to δ 13 C SE(t−PI) for the full industrial period, the assumption is used that the oceanic δ 13 C DIC change scales with the atmospheric δ 13 CO 2 change, i.e.: δ 13 C SE(t−PI) = f atm · δ 13 C SE(t−1940) = f atm · a · pCFC-12 t , with f atm = δ 13 CO 2,t − δ 13 CO 2,PI δ 13 CO 2,t − δ 13 CO 2,1940 .
To achieve a result comparable to E17, we select the model data at the locations for which both CFC-12 and δ 13 C DIC measurements are available. Here we use the observations compiled by Schmittner et al. (2013) because δ 13 C DIC in this data set has been quality controlled and is publicly available. The only difference with respect to the observational data used in E17 is that Schmittner et al. (2013) do not include the data at one South Atlantic section (A13.5) measured in 2010. However, this 420 difference does not affect our analysis. Vertically, we use data at the model layers above the simulated pCFC-12 penetration depth (set at 20 patm, following E17). We take t = 1994 and perform the linear regression of Eq. (18) for five ventilation regions (the North Atlantic, South Atlantic, North Pacific, South Pacific and Indian Ocean), respectively.
The regressional relationships between δ 13 C pref 1994 and pCFC-12 1994 (Eq. 18) and the regression coefficients, hereafter referred to as a pref and b pref , are shown in Fig. D5 (the water masses in this figure are defined in Table D1). The coefficient of 425 determination r 2 , 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 2 = 0.22 in E17. Applying Eq. (19) to the three-dimension model data of pCFC-12 1994 for t = 1994, a pref and f atm = 1.5, we obtain the estimate of the global oceanic 13 C Suess effect in year 1994, which we refer to as SE pref .
To quantify if SE pref under-or overestimate the oceanic 13 C Suess effect, we compare SE pref to the simulated oceanic 13 C At 200 m SE pref mostly underestimates SE Mod (Fig. 11a). The ventilation 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). These findings confirm the underestimation range discussed by E17. Note that E17 deduced the range of  Broecker and Peng (1993) and Bacastow et al. (1996) E17 assumed an ocean-to-atmosphere ratio of the 13 C Suess effect of 0.65 and the 200 m-to-surface ratio of the 13 C Suess effect of 0.6-0.7. Multiplying the above two ratios with the atmospheric δ 13 CO 2 decrease of −1.4‰ by year 1994 yields global-mean 13 C Suess effect -0.55 to −0.64‰ at 200 m. In our model, the global-mean ocean-to-atmosphere ratio of the 13 C Suess effect is 0.46, significantly lower than the five-box model of Broecker 440 and Peng (1993). On the other hand, our model shows a slightly higher 200 m-to-surface ratio of the 13 C Suess effect (0.75) than Bacastow et al. (1996) who employed an ocean general circulation model with coarse vertical resolution (4 layers for the upper 200 m).
E17 have speculated that the major cause of the underestimation of oceanic 13 C Suess effect is that the available observations are mostly from the intermediate and deep waters. The ocean-atmosphere equilibration timescale for δ 13 C (10 years, Broecker 445 and Peng, 1974) is significantly longer than that of pCFC-12 (1 month, Gammon et al., 1982). Thus, waters that have shorter surface residence time, such as the deep waters ventilated in the South Hemisphere, would show less negative slope a pref b c a Figure 11. Distribution at 200 m depth for SEpref − SEMod (a), SEtotal − SEMod (b) and SEpref − SEtotal (c). The isoline increment is 0.1‰. In panels b and c, the South Pacific Ocean is not presented because the relationship between the total oceanic 13 C Suess effect δ 13 C SE(1994SE( −1940 and pCFC-12 1994 is too weak (r 2 = 0.07) and therefore SEtotal can not be estimated. ventilation regions because our model has higher vertical resolution in the upper ocean and therefore has more data points than field measurements. For the Indian Ocean, we combine the model data from Subtropical Gyre Water and Sub-Antarctic Mode Water as both water masses have a strong 13 C Suess effect (E17). We find in the Indian Ocean a pref for the Subtropical Gyre Water and Sub-Antarctic Mode Water (−0.65 × 10 −3 , r 2 = 0.49) is more negative than that for the whole ventilation  To reveal the source for the underestimation, we divide (SE pref −SE Mod ) into two components (SE pref −SE total ) and (SE total − SE Mod ). Here SE total is derived similarly to SE pref and it is based on linear regression relationships between δ 13 C SE(1994SE( −1940 and pCFC-12 1994 . Positive values of (SE pref − SE total ) show the underestimation of oceanic 13 C Suess effect induced by using a pref · pCFC-12 1994 to approximate the relationship between the δ 13 C SE(1994SE( −1940 and pCFC-12. The difference (SE total − SE Mod ) shows how well a method based on linear relationships between the δ 13 C SE(1994SE( −1940 and pCFC-12 1994 can estimate 465 the global ocean 13 C Suess effect. To calculate SE total we perform a linear regression for the total oceanic 13 C Suess effect δ 13 C SE(1994SE( −1940 and pCFC-12 1994 with the subsampled model data: SE(1994SE( −1940 ∼ a total · pCFC-12 1994 + b total . (21) This is followed by scaling from the period 1940 -1994 to the full industrial period in analogy to Eq. (19): The regression relationships in Eq. (21) and regression coefficients are given in Fig. D6. For the Indian, North Pacific, North Atlantic and South Atlantic Ocean, r 2 lies between 0.34 and 0.67, which suggests acceptable strength of the relationships. In the South Pacific Ocean we find low r 2 = 0.07, and therefore we don't compute SE total for the South Pacific Ocean. The causes for this low r 2 will be discussed later in this section.

475
With Eqs. (19) and (22) we get and it is rather negative poleward of 40 • (Fig. 11b). This pattern results from lumping together data from different water masses to generate one regression relationship for a large ventilation region. The waters ventilated in lower latitudes typically have 480 stronger 13 C Suess effect than those ventilated in high latitudes. This is clearly reflected in the linear regression relationships between δ 13 C SE(1994SE( −1940 and pCFC-12 1994 for the North Atlantic Ocean (Fig. D6d), which shows 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 underestimation of 13 C Suess effect (positive values of SE total − SE Mod ) and the water masses ventilated in the high latitudes show overestimation (Figs. 12d -12f). In the North Atlantic Ocean, the 485 region-mean underestimation SE pref − SE Mod = 0.1‰ is predominantly contributed by SE total − SE Mod = 0.09‰. In the North Pacific Ocean SE total − SE Mod = 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 the South Atlantic, North Pacific and Indian Ocean, (SE pref − SE total ) is always positive and it decreases with increasing depth (Figs. 11c, 12g -12i) because pCFC-12 decreases towards the interior ocean (see Eq. 23). In the South Atlantic and 490 Indian Ocean, (SE pref − SE total ) determines the region-mean underestimation (Table 1). In the North Pacific Ocean, it contributes to less than half of the underestimation (Table 1). We find two main causes of the underestimation from the component (SE pref − SE total ) in the above three regions. The first arises from the assumption that δ 13 C pref 1940 is a constant in the regression equation Eq. (18). As is shown for the zonal-averaged vertical sections, δ 13 C pref 1940 exhibits noticeable spatial variations (Figs. 13a -13c). Over a considerable fraction of ocean regions (e.g. north of 40 • S in the South Atlantic Ocean, south of 35 • N in the North Pacific Ocean, north of 40 • S in the Indian Ocean) δ 13 C pref 1940 generally decreases with increasing depth. This vertical distribution is similar to that of pCFC-12 1994 (Figs. 13d -13f). Such a seemingly positive correlation between δ 13 C pref 1940 and pCFC-12 1994 exists because during the pre-industrial times, in our model the preformed component δ 13 C pref generally decreases with increasing water depth, which has also been reported by Schmittner et al. (2013) (see their Figs. 5 and 6). In the industrial period prior to 1940, the decrease of the atmospheric 13 C/ 12 C ratio is relatively slow (Fig. 2). Thus, by the year 1940 500 the oceanic uptake of isopotically light CO 2 only partly offsets vertical gradient of the pre-industrial δ 13 C pref . Consequently, the seemingly positive correlation between δ 13 C pref 1940 and pCFC-12 1994 results in less negative a pref than a total . This therefore generates the underestimation of 13 C Suess effect in the South Atlantic, North Pacific and Indian Ocean in the model.
history of 13 C Suess effect than CFC-12, as is already discussed by E17. This point is supported by non-negligible values of the regression intercept b total in our study. In the South Atlantic and Indian Ocean, b total = −0.07‰ corresponds to an underestimation of 0.11‰ (= −f atm · b total , see Eq. 23). Thus, neglecting 13 C Suess effect in CFC-12 water contributes to almost half of the total underestimation at 200 m for the Indian Ocean and it contributes to about 80% for the South Atlantic Ocean (Table 1).

510
Different from the previously discussed three ventilation regions, in the North Atlantic Ocean negative (SE pref − SE total ) is found at 200 m (Fig. 11c), which becomes positive below about 250 m (Fig. 12g). The reason is a pref = −0.81 × 10 −3 is more negative than a total = −0.62 × 10 −3 (Figs. D5d and D6d). This is related to the assumption that the regenerated component can be mainly attributed to the change of organic matter remineralization in the ocean interior, as is illustrated by the temporal change of AOU (Figs. 13j -13l). Below 1500 m, the δ 13 C reg changes are generally negative (Figs. D7d -D7f) because δ 13 C POC decreases globally by about 1.3‰ during 1940-1994 (Fig. D8), mainly due to the decline of the biological fractionation factor p under increasing surface CO 2 (aq). In the North Atlantic Ocean, −(δ 13 C reg 1994 − δ 13 C reg 1940 ) is mostly negative above 500 m, where pCFC-12 1994 is relatively high (Fig. 13g). Below 500 m, where pCFC-12 1994 is relatively 520 low, −(δ 13 C reg 1994 − δ 13 C reg 1940 ) is mostly positive. Thus, an apparent negative correlation between the spatial distributions of −(δ 13 C reg 1994 − δ 13 C reg 1940 ) and pCFC-12 1994 leads to a more negative a pref than a total , according to Eq. (16). The consequential overestimation of δ 13 C SE(1994−PI) by 0.09‰ (= f atm · (a pref − a total ) · pCFC-12 1994 ) is compensated by a underestimation of 0.12‰ (= −f atm · b total ) due to the negative linear regression intercept b total = −0.08‰. Overall we find a negligible underestimation of mean (SE pref − SE total ) = 0.02‰ at 200 m depth in the North Atlantic Ocean.

525
The temporal change of δ 13 C reg also causes the positive values of δ 13 C SE at depth, for instance, below 1000 m on the vertical section P16 in the South Pacific Ocean (Fig. 10b). Here the positive change of δ 13 C reg is due to a decrease of remineralisation, as is shown by the change of AOU in this region (Fig. D8h). Hence, less negative δ 13 C POC signal is released in this region and δ 13 C DIC slightly increases.
We don't compute SE total for the South Pacific Ocean because of a low r 2 = 0.07 obtained for the linear regression between 530 δ 13 C SE(1994SE( −1940 and pCFC-12 1994 , suggesting no linear relation between the two variables (Fig. D6c). The model data for the Subtropical Gyre Water and the Antarctic Intermediate Water in the South Pacific Ocean are particularly scattered (Fig. D5f) because −(δ 13 C reg 1994 − δ 13 C reg 1940 ) shows significant spatial variability within each of these two water masses (Fig. D7e). These changes of δ 13 C reg are mainly caused by the changes of ocean carbon cycle in our model, as is illustrated by the AOU changes (Fig. D7h).

535
Although we can not directly evaluate (SE pref − SE total ) for the South Pacific Ocean, we can try to understand the total underestimation (0.26‰ at 200 m, being the largest among the five ventilation regions, see Table 1) by analysing the spatial distribution of the terms in Eq. (16). A seemingly positive correlation clearly exists between δ 13 C pref 1940 and pCFC-12 1994 , which both decreases with increasing depth (Fig. 13b). This contributes to the underestimation of 13 C Suess effect (according to Eq. 23), similar to the case for the Indian, South Atlantic and North Pacific Ocean. Above the pCFC-12=20 patm isoline, −(δ 13 C reg 1994 − δ 13 C reg 1940 ) mostly decrease with increasing depth (Fig. 13h), similar to the Atlantic Ocean, which contributes to an overestimation.
Note that due to inevitable model biases, our model does not perfectly reproduce the distribution and properties of the observed water masses (see more discussion in Appendix C). Thus our regression relationships between δ 13 C pref and pCFC-12 (Fig. D5) show some quantitative differences to those of E17 (see their Fig. 3). Nevertheless, our analysis provides a possible spatial 545 distribution of the underestimation of 13 C in E17's product. More importantly, we uncover two major causes for the underestimation: the assumption of a spatially constant δ 13 C pref 1940 and the neglect of 13 C Suess effect in CFC-12 free water.

Summary and conclusions
We present results of the new 13 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: Popp p depends on dissolved CO 2 (Popp et al., 1989); Laws p is a function of dissolved CO 2 and phytoplankton growth rate (Laws et al., 1995). Futhermore, we used our consistent model framework to assess the approach by Eide et al. (2017a), who used a correlation between preformed δ 13 C DIC and CFC-12 partial pressure to obtain an estimate of the global oceanic 13 C Suess effect.
The comparison between simulated and observed isotopic ratio of organic matter δ 13 C POC reveals that Popp p (r = 0.84 and 555 NRMSE = 0.57) has a better performance than Laws p (r = 0.71 and NRMSE = 2.5). Using Laws p results in noticeably lower δ 13 C POC values and smaller δ 13 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 a too strong preference of the isotopically light carbon. Therefore it is not a good representative for 13 C biological fractionation in our global ocean biogeochemical model.

580
We conclude that the new 13 C module with biological fractionation factor Popp p from Popp et al. (1989) has a satisfactory performance. Thus, our new 13 C module will serve as a useful tool to evaluate the performance of MPI-ESM in paleo-climate 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 13 C inventory changes The water-column DI 13 C inventory difference is primarily a result of the difference of the net air-sea 13 CO 2 flux between 585 PI_Popp and PI_Laws. This is demonstrated by the comparison of the contributions of the governing factors for the watercolumn DI 13 C inventory changes (Table A1), including air-sea gas exchange, loss of POC and CaCO 3 to marine sediment, diffusion of the remineralised DIC from sediment into the water column, input of DOC and CO 2− 3 , and the exchange with other marine carbon pools (phytoplankton, CaCO 3 , etc.). Table A1 also reveals that the current method to determine the 13 C input (see Section 2.3.2) only has a small contribution to the change of the water-column DI 13 C inventory. 590 Table A1. Contributions to the rate of the water-column DI 13 C inventory change (in Gmol yr −1 ), averaged in the last 50 years in the corresponding pre-industrial spin-up simulations. Positive values denote contributions to the increase of the water-column DI 13 C inventory.
Last column gives relative contribution to the total rate difference with relative contribution = (PI_Laws-PI_Popp) / total rate difference. B1d) and with field observations. In the central equatorial Pacific the simulated µ well reproduces the observed range (0.55-0.7 day −1 , 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 day −1 ) is at the lower side of both the observations (annual mean 0.3-0.53 day −1 in the North Pacific subtropical gyre, Letelier et al., 1996; annual mean 0.13-0.62 day −1 in the North Atlantic subtropical gyre, Marañón, 2005) and the satellite-based µ estimates. In the Pacific sector of the Southern Ocean, the simulated   E17's procedure first assumes that any oceanic CFC-12 signal before 1940 is negligible and the oceanic 13 C Suess effect at any time t after 1940 is proportional to CFC-12 partial pressure at time t: Here the proportionality factor a is time-invariant. δ 13 C DIC at any time t after year 1940 is decomposed as: δ 13 C t = δ 13 C pref 1940 + δ 13 C reg 1940 + δ 13 C SE(t-1940) + ∆δ 13 C reg + ∆δ 13 C pref . (C2)

610
The superscript "pref" represents the preformed component, which arises from the transport of the surface water with specific DIC and DI 13 C. Superscript "reg" denotes the regenerated component δ 13 C reg SE due to organic matter remineralisation and calcium carbonate dissolution. The two last terms contain any changes not related to the 13 C Suess effect, e.g. changes in ocean carbon cycle. Decomposing the left-hand side of Eq. (C2) gives Note here −(δ 13 C reg t −δ 13 C reg 1940 )+∆δ 13 C reg represents any change of the regenerated component and is equivalent to −(δ 13 C reg t − δ 13 C reg 1940 ) in Eq. (16) in Section 4.2. In E17, the terms (δ 13 C reg t − δ 13 C reg 1940 ), ∆δ 13 C reg and ∆δ 13 C pref are assumed zero. Combining Eq. (C1) and Eq. (C3) yields where b contains terms δ 13 C pref 1940 . Thus, the proportionality factor a can be determined with δ 13 C pref t and ·pCFC-12 t at time t, 620 and δ 13 C SE(t-1940) is obtained with Eq. (C1).
The preformed component is calculated following Sonnerup et al. (1999) and Eide et al. (2017a): This equation is only valid below the 200 m, which is roughly the euphotic zone depth (E17). The C -O2 org ratio is 122:172 in HAMOCC6, and we use the simulated δ 13 C POC for δ 13 C org . The dissolution of CaCO 3 is neglected following Sonnerup et al.
For mapping E17 performed another linear regression for δ 13 C SE(t−PI) and pCFC-12 t . This step is not need here as we take t = 1994 rather than various years from observations. To obtain three-dimention distribution of 13 C Suess effect for the global ocean, we apply Eq. (C6) to the simulated pCFC-12 1994 .
The regression relationships between δ 13 C pref and pCFC-12 in our model (Fig. D5) show some quantitative differences to 635 those of E17 (see their Fig. 3). The reason is our model does not perfectly reproduce the distribution and properties of the observed water masses and this can be seen in the following aspects. First, the definitions of several water masses in the model are slightly different from those of E17 (comparing our for this difference are two fold. First, in the Southern Ocean 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. D2), and by the zonal-mean CFC-12 bias distribution (Fig. D9), which features persistent positive biases off the Antarctic continental shelf in the Atlantic, Pacific and Indian sectors of the Southern Ocean. The second cause is the Southern Ocean has a too high primary production (about a factor of 1.5 of the satellite-based net primary production estimates from Westberry et al., 2008). The high 645 primary production causes higher surface δ 13 C DIC than observations (see the South Pacific Ocean in Fig. 8c). Consequently, the simulated preformed component δ 13 C pref t in the bottom and deep water masses of the Southern Ocean is higher than observed values in E17. Third, the lowest values of δ 13 C pref t (< 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 δ 13 C DIC than observations (Figs. 9e and 9f).     * Water masses are combined together rather than separately defined as in Eide et al. (2017a).
** A different σ θ threshold is used here compared to Eide et al. (2017a). c) e) a) Figure D5. Regressional relationships between pCFC-12 1994 and δ 13 C pref 1994 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 indicates different water masses. The full names, as well as the definitions, of the water masses are listed in Table D1. c) e) a) Figure D6. As Fig. D5, but for the relationships between pCFC-12 1994 and δ 13 CSE .