Carbon isotopes in the marine biogeochemistry model FESOM2.1-REcoM3

. In this paper we describe the implementation of the carbon isotopes 13 C and 14 C (radiocarbon) 10 into the marine biogeochemistry model REcoM3. The implementation is tested in long-term equilibrium simulations where REcoM3 is coupled with the ocean general circulation model FESOM2.1, applying a low-resolution configuration and idealized climate forcing. Focusing on the carbon-isotopic composition of dissolved inorganic carbon (δ 13 C DIC and Δ 14 C DIC ), our model results are largely consistent with reconstructions for the pre-anthropogenic period. Our simulations also exhibit discrepancies, e.g., in 15 upwelling regions and the interior of the North Pacific. Some of these differences are due to the limitations of our ocean circulation model setup which results in a rather shallow meridional overturning circulation. We additionally study the accuracy of two simplified modelling approaches for dissolved inorganic 14 C, which are faster (15 % and about a factor of five, respectively) than the complete consideration of the marine radiocarbon cycle. The accuracy of both simplified approaches is better than 5 % which should 20 be sufficient for most studies of Δ 14 C DIC .


Introduction
Carbon isotopes are powerful tools for tracing present and past biogeochemical cycles and water masses.The stable isotope carbon-13 ( 13 C) can be used to study the anthropogenic perturbation of the global carbon cycle owing to the combustion of isotopically depleted fossil fuels via the so-called 13 C Suess effect (e.g.Keeling, 1979;Quay et al., 1992;Köhler, 2016;Graven et al., 2021).Carbon-13 may also help to decipher the exchange between atmospheric, marine, and terrestrial carbon reservoirs in the past, for example, during the last glacial cycle (Köhler et al., 2006;Ciais et al., 2012;Broecker and McGee, 2013;Jeltsch-Thömmes et al., 2019;Menking et al., 2022).Furthermore, 13 C is a proxy for oceanic nutrients and can be employed to infer past marine biological productivity, export production, and water mass distributions assuming that calcareous tests of foraminifera are faithful recorders of dissolved inorganic 13 C (e.g.Broecker and Peng, 1982;Lynch-Stieglitz et al., 2007;Hesse et al., 2011;Schmittner et al., 2017).The radioactive isotope carbon-14 (radiocarbon, 14 C) is the most important geochemical chronometer for dating organic matter and for assessing ocean ventilation rates and pathways over the last 55 000 years (Heaton et al., 2021;Rafter et al., 2022;Skinner andBard, 2022, Skinner et al., 2023).In addition, the penetration of bomb-produced 14 C into the oceans has provided a benchmark for ocean circulation models (Matsumoto et al., 2004).Although the potential of carbon isotopes has been known for a long time, ocean general circulation models have only recently been equipped with carbon isotopes and applied in Earth system modelling studies (e.g.Tschumi et al., 2011;Holden et al., 2013;Schmittner et al., 2013;Jahn et al., 2015;Menviel et al., 2015;Buchanan et al., 2019;Jeltsch-Thömmes et al., 2019;Dentith et al., 2020;Tjiputra et al., 2020;Claret et al., 2021;Liu et al., 2021;Morée et al., 2021).
2 Model description

Short overview of REcoM3
The Regulated Ecosystem Model, version 3 (REcoM3) is described in detail by Gürses et al. (2023).Here, we only give a brief summary of the common model features and describe the differences from the configuration that we use in this study.REcoM considers the marine biogeochemical cycles of carbon, nitrogen, silicon, iron, and oxygen.The ecosystem component of REcoM includes nutrients, two phytoplankton functional types (distinguishing between small phytoplankton and diatoms), one zooplankton functional type, one detritus type, and dissolved organic matter.Different from most other marine biogeochemistry models, REcoM does not rely on a fixed internal stoichiometry of phytoplankton.Instead, the composition of organic soft tissue is regulated (i.e.calculated) in response to light, temperature, and nutrient supply, which enables stoichiometric shifts between the present and past to be assessed.The model includes a sediment layer in which sinking detritus (consisting of particulate organic matter, calcite, and opal) is fully remineralized and where solutes are returned to the bottom water layer within days so that there is no accumulation of sedimentary matter.Alternatively, REcoM3 can be run with the comprehensive sediment model MEDUSA2 (Munhoven, 2021), which is described in another paper (Ye et al., 2023).Apart from the implementation of carbon isotopes, the main difference from the standard REcoM3 configuration by Gürses et al. ( 2023) is that RE-coM3 considers two zooplankton and two detritus groups instead of one, which would require inclusion of six additional carbon-isotopic tracers at the expense of model speed.We refer to the configuration presented here as REcoM3p as an initial set-up for paleo studies.To simulate biogeochemical tracer circulation, REcoM3 needs a transport model.Here, we utilize the ocean general circulation model FESOM2.1, which is an update of FESOM2.0 (Danilov et al., 2017).The coupling of FESOM2.1 with REcoM3 is also described by Gürses et al. (2023), where all biogeochemical model equations except for carbon isotopes can be found.

Implementation of carbon isotopes
Carbon-13 and 14 C are implemented as additional passive tracers, tripling the number of carbon-containing tracers in REcoM from 8 to 24.The model considers carbon isotopes for dissolved inorganic and organic carbon, small phytoplankton, diatoms, zooplankton, detritus, and calcite.As the abundances of 13 C and 14 C are small ( 13 R std ≈ 0.0113, Assonov et al., 2020, and 14 R std = 1.170 × 10 −12 ; Orr et al., 2017), we approximate total C with 12 C and neglect 13 C and 14 C in the kinetic calculations of the carbonate system.For the same reason, and to ensure numerical stability in the other model calculations involving 13 C and 14 C, the model considers scaled concentrations 13 C = 13 C/ 13 R std and 14 C = 14 C/ 14 R std , which have the same magnitude as 12 C, following Jahn et al. (2015).As a consequence, 13 C and 14 C have to be rescaled when compared with observed concentrations.However, as the scaling factors cancel out when 13 C and 14 C are converted to δ 13 C and 14 C, Eqs. ( 1)-( 3) can be directly evaluated with scaled concentration ratios and 13 R std = 14 R std = 1.We consider isotopic fractionation during air-sea gas exchange, dissolution of CO 2 in seawater, and photosynthesis by phytoplankton.In addition, the model accounts for radioactive decay and, optionally, cosmogenic production of 14 C.The details are explained in the following subsections.

Air-sea exchange
The isotopic fractionation during uptake and dissolution of 13 CO 2 is calculated according to Zhang et al. (1995), similar to the biogeochemical protocol of the CMIP6 Ocean Model Intercomparison Project (OMIP-BGC protocol, Orr et al., 2017).That is, the air-sea exchange flux 13 F is proportional to the difference between the saturation and in situ concentrations of aqueous 13 CO 2 : where k w is the CO 2 gas transfer velocity (calculated according to Wanninkhof, 2014, additionally considering seaice cover), and [ 13 CO * 2 ] sat and [ 13 CO * 2 ] are the saturation and in situ concentrations of aqueous 13 CO 2 . 13R atm and 13 R DIC are the 13 C/ 12 C concentration ratios of atmospheric CO 2 and dissolved inorganic carbon (DIC).Isotopic fractionation factors α denote kinetic fractionation during CO 2 gas transfer ( 13 α k ), equilibrium fractionation during gas dissolution ( 13 α aq−g ), and equilibrium fractionation between DIC and gaseous CO 2 ( 13 α DIC−g ).Numerical values were taken from Zhang et al. (1995), who measured kinetic and equilibrium fractionation of 13 C in acidified freshwater.For kinetic fractionation we employ a constant value of 13 α k = 0.99912, which is the average between 5 °C ( 13 α k = 0.99919) and 21 °C ( 13 α k = 0.99905).Equilibrium fractionation between aqueous and atmospheric 13 CO 2 is expressed as 13 α aq−g = 1 + 0.001(0.0049T− 1.31), (5) where T (°C) is the temperature of surface water.The numerical values of 13 α k and 13 α aq−g were determined for fresh water.To account for enhanced 13 C fractionation in seawater associated with hydration reactions, Eq. ( 4) includes a correction factor of −0.0002 following Zhang et al. (1995), which is not considered in the OMIP-BGC protocol.Fractionation between DIC and gaseous CO 2 is calculated as where f CO 3 = [CO 2− 3 ] / DIC is the carbonate fraction of DIC.

Biogenic fractionation
Photosynthesis of phytoplankton leads to isotopic depletion of particulate organic carbon (POC), which is expressed following Freeman and Hayes (1992): where [ 13 C POC ] and [ 12 C POC ] are the isotopic POC concentrations in phytoplankton, 13 R POC is the corresponding isotopic ratio, 13 R aq is the 13 C/ 12 C concentration ratio of aqueous CO 2 , and 13 α p is the isotopic fractionation factor associated with photosynthesis.Various experimental and modelling studies have determined 13 α p for certain phytoplankton species (Laws et al., 1995;Rau et al., 1996;Bidigare et al., 1997;Laws et al., 1997;Rau et al., 1997;Popp et al., 1998;Keller and Morel, 1999).However, it is uncertain to what extent these studies are globally representative and can be transferred into a single global modelling framework (see also the review by Brandenburg et al., 2022).Therefore, we pursue a less sophisticated but robust approach to calculating 13 α p , which does not rely on species-specific assumptions and has been inferred from a global dataset of 525 δ 13 C POC field measurements spanning the period 1962-2010 CE (Young et al., 2013): where [CO * 2 ] is in µmol L −1 .As no distinction is made between different phytoplankton species in Eq. ( 8), the carbonisotopic composition of small phytoplankton and diatoms in our model is the same.Similar to other models (e.g.Schmittner et al., 2013;Menviel et al., 2015;Buchanan et al., 2019;Tjiputra et al., 2020;Liu et al., 2021), we do not consider carbon-isotopic fractionation during the formation and dissolution of biogenic calcite because the effect is small and varies between species (α ∼ 0.999-1.003according to Ziveri et al., 2003).

Radiocarbon
Radiocarbon is subject to radioactive decay, cosmogenic production, and isotopic fractionation.The radioactive decay (applying a half-life of 5700 years, Audi et al., 2003;Bé and Chechev, 2012;Kutschera, 2013) is balanced in the model by cosmogenic 14 C production fluxes or, alternatively, by prescribed atmospheric 14 CO 2 concentrations corresponding to atmospheric 14 C values.
Fractionation factors are calculated analogously to 13 α, with 14 α = 2 13 α − 1 (e.g.Craig, 1954; see also the discussion by Fahrni et al., 2017).The slow equilibration of DI 14 C in the deep ocean requires spinup simulations of several thousands of years, which may be computationally too expensive for high-resolution simulations.Therefore, we have implemented 14 C in various ways that differ in their level of complexity and computational efficiency.The first implementation ("CC") considers the complete 14 C cycle parallel to 13 C.The second implementation ("IC") disregards the isotopic fractionation of 14 C and radioactive decay of organic matter.In turn, DI 14 C and DIC have identical sources and sinks, except for atmospheric CO 2 and radioactive decay.This "inorganic" 14 C approach is an approximation that is conceptually similar, but not identical, to the "abiotic" 14 C modelling approach described in the OMIP biogeochemical protocol (Orr et al., 2017).In our IC approach, DIC and DI 14 C include biological sources and sinks.This does not apply to the OMIPabiotic approach, which also considers alkalinity in a simplified way.In Sect.3.3 we scrutinize the accuracy of the IC approximation.The third approach takes advantage of the fact that marine 14 C values of DIC ( 14 C DIC ) are primarily governed by transport and radioactive decay (Fiadeiro, 1982; see also Mouchet, 2013).This implies that 14 C DIC can be implemented into ocean general circulation models as a single prognostic tracer without a full carbon cycle model (see Toggweiler et al., 1989, and numerous other studies later on).We evaluate this 14 C approximation ("DA") in an additional simulation and compare it with REcoM approaches CC and IC in Sect.3.3 too.
A list of the various model experiments and their key features can be found in Table 1.

Simulated ocean climatology
FESOM employs unstructured meshes with variable horizontal resolution.The default mesh of FESOM2.1-REcoM3 has about 127 000 horizontal surface nodes (Gürses et al., 2023).Here, our model resolution is radically reduced, considering 3140 surface nodes corresponding to a median horizontal resolution of 260 km (the so-called pi-mesh, see Fig. A1 in Appendix A).This configuration requires fewer computational resources and allows simulations over the time scale of mahttps://doi.org/10.5194/gmd-17-1709-2024 Geosci.Model Dev., 17, 1709-1727, 2024 rine carbon isotope equilibration to be performed (i.e. over several thousand years) within a few weeks.Vertical resolution is 47 layers using z* coordinates with non-linear free surface (for further details see Scholz et al., 2019).
In a first step we run FESOM without REcoM over 1000 years to spinup the overturning circulation and thermohaline fields.FESOM was initialized with seasonal winter temperatures and salinities by Steele et al. (2001), and driven with annually repeated atmospheric fields using Corrected Normal Year Forcing Version 2.0 (Large and Yeager, 2009; for an overview see also Griffies et al., 2012).As FE-SOM2.1 had originally been adjusted for higher resolution and different forcing, we retuned the model in our setup by reducing the maximum Gent-McWilliams thickness diffusivity from 2000 m 2 s −1 originally to 1000 m 2 s −1 .After 1000 simulated years there was no significant drift of thermohaline fields and the meridional overturning circulation (MOC) had stabilized, REcoM3p was turned on, and both models were run over an additional 5000 years.At this moment, temperature and salinity drifted by about 10 −6 K per year and 10 −6 PSU per year in the global average, and the major biogeochemical tracer inventories drifted by less than 10 −3 % per year (less than 10 −3 ‰ per year regarding DI 13 C and DI 14 C).REcoM3p was initialized with gridded concentration fields of total alkalinity, DIC, and nitrate from version 2 of the Global Ocean Data Analysis Project (GLODAPv2, Key et al., 2015;Lauvset et al., 2016), oxygen from the World Ocean Atlas 2018 (WOA18, Garcia et al., 2019), silicate from Ocean Atlas 2013 (WOA13, Garcia et al., 2014), and dissolved iron according to Aumont et al. (2003) and Tagliabue et al. (2012).Dissolved inorganic 13 C and 14 C were initialized with DIC, assuming initial fractionation values of δ 13 C DIC = 0 ‰ and 14 C DIC = −150 ‰ respectively.In the sediment layer all initial concentrations were close to zero.REcoM3p was forced with constant atmospheric CO 2 concentrations ([ 12 CO 2 ] atm = 284.3ppmv; [ 13 CO 2 ] atm and [ 14 CO 2 ] atm corresponding to δ 13 C atm = −6.61‰ and 14 C atm = 0 ‰ respectively, following Orr et al., 2017), and with climatological mean monthly fluxes of aeolian iron deposition (Albani et al., 2014).The CO 2 concentration forcing implies that carbon-isotopic mass is only conserved in the atmosphere-ocean system when the marine isotopic inventories have reached a corresponding steady state.In our simulations this is the case after a few thousand years (see further below).
As discussed in the following, our low-resolution test setup sufficiently captures the basic thermohaline circulation structures obtained with FESOM2.0 in higher-resolution simulations (Scholz et al., 2019(Scholz et al., , 2022;;compare their Figs. 14, 15, and 17 with our Figs.A2-A4 in Appendix A).Compared with observations, our simulations exhibit a warm bias for thermocline and intermediate water as well as for North Atlantic Deep Water (NADW; see Fig. A2 in Appendix A).The most notable exception is the region of the North Atlantic Gulf stream, which is not properly resolved, and where upper ocean temperatures are considerably underestimated.The temperature biases covary with salinity biases.That is, simulated salinities are also higher than observations where simulated temperatures are high, and salinities are too low where simulated temperatures tend to be low (Fig. A3).In the Atlantic, FESOM arrives at a MOC of about 16 Sv (1 Sv = 1 × 10 6 m 3 s −1 , Fig. A4), which is at the lower bound of observations, whereas the simulated overturning cell is too shallow compared with observational estimates (see Buckley and Marshall, 2016, and further references therein).For the purposes of this study, our simulations also reasonably reflect the observed large-scale pattern of DIC (Key et al., 2015;Lauvset et al., 2016).That is, low concentrations are found at upper levels of the subtropical oceans and in the freshly ventilated interior of the North Atlantic, whereas DIC concentrations progressively increase in the South Atlantic and in the deep North Pacific (Fig. A5).

Carbon-13
We now focus on the carbon-isotopic composition of DIC near the sea surface and along zonal-mean sections through the Atlantic and Pacific.Regarding δ 13 C DIC , we compare our model results with the gridded preindustrial (PI) δ 13 C DIC climatology by Eide et al. (2017b).Their reconstruction does not consider the upper 200 m to exclude the 13 C Suess effect.For the sake of completeness, we briefly note that at the sea surface our δ 13 C DIC results are largely in line with ungridded preindustrial values determined by Kwon et al. (2022; shown in Fig. A6).
In wide areas, our simulated near-surface δ 13 C DIC values at 200 m are within the range 1 ‰ to 2 ‰ (Fig. 1a).Higher δ 13 C DIC values are simulated in the subtropical oceans, particularly in the Atlantic, southeast Pacific, and southern Indian Ocean.Isotopic depletion of up to ∼ −1 ‰ is found in the eastern equatorial Pacific, the subpolar North Pacific, the Bay of Bengal, and in the Angola Basin.Although our simulation captures the reconstructed spatial pattern (shown in Fig. 1b), the simulated range of δ 13 C DIC variations is larger than in the reconstruction, that is, the model results are too high in the lower latitudes and too low in the abovementioned depletion regions (Fig. 2a).
Considering the interior of the Atlantic Ocean, our simulation yields higher δ 13 C DIC in the North Atlantic than in the South Atlantic, and small vertical δ 13 C DIC gradients in the high latitudes of both hemispheres (Fig. 1c).In the Pacific, our simulation displays a reversed meridional δ 13 C DIC gradient with the strongest isotopic depletion at intermediate depths in the North Pacific.Our results are roughly in line with the reconstruction by Eide et al. (2017b, see Fig. 1d) but overestimate δ 13 C DIC in the upper Pacific at low latitudes as well as in the North Pacific at a depth of between 1.5 and 3 km (Fig. 2b).Similar inaccuracies are seen in other mod-els (Tagliabue an Bopp, 2008;Schmittner et al., 2013;Jahn et al., 2015;Buchanan et al., 2019;Jeltsch-Thömmes et al., 2019;Dentith et al., 2020;Tjiputra et al., 2020;Liu et al., 2021).
The differences between simulated and reconstructed δ 13 C DIC (the δ 13 C DIC bias for short) can be attributed to various factors.First, the reconstructed δ 13 C DIC values may be too low in the upper ocean if the 13 C Suess effect is not completely removed from the observations.In fact, Eide et al. (2017a) presume that they may have underestimated the 13 C Suess effect at a depth of 200 m by about 0.1 ‰ to 0.2 ‰, which has been corroborated in recent modelling studies (Liu et al., 2021;Kwon et al., 2022).This reconstruction bias partly corresponds to the apparent simulation bias at a depth of 200 m, where our model results are 0.2 ‰ higher on global average than the reconstructed results.
To some extent, the δ 13 C DIC bias can be further attributed by decomposing δ 13 C DIC into biologic and thermodynamic sources: where δ 13 C BIO specifies the imprint of isotopic fractionation, respiration, and remineralization of organic matter in the absence of air-sea exchange.Following Broecker and Maier-Reimer (1992), δ 13 C BIO is frequently determined from the covariation of δ 13 C DIC with marine phosphate (PO where PO 3− 4 and DIN are in µmol kg −1 .Equation ( 10) does not distinguish between preformed and remineralized nutrient concentrations and cannot be used for further separation of δ 13 C BIO into preformed and remineralized components.The thermodynamic component δ 13 C AS describes the effects of air-sea exchange and ocean circulation and is the residual of observed δ 13 C DIC and reconstructed δ 13 C BIO .
In the following we validate simulated δ 13 C BIO and δ 13 C AS with the corresponding values reconstructed by Eide et al. (2017b), and we compare the δ 13 C DIC bias with the differences between the simulated and reconstructed δ 13 C DIC components, i.e. the δ 13 C BIO bias and δ 13 C AS bias.Simulation CC overestimates δ 13 C BIO in wide areas, except for subtropical and upwelling regions at 200 m and the North Atlantic below 3 km (Fig. 3).The simulation also overestimates δ 13 C AS in the thermocline but underestimates δ 13 C AS in the interior of the North Atlantic above about 3 km (Fig. 4).Comparing the biases of δ 13 C DIC , δ 13 C BIO and δ 13 C AS , it appears that the δ 13 C DIC bias corresponds to the δ 13 C BIO bias in the low latitudes, upwelling regions, and the interior of the Atlantic (see Figs. 2 and 3).This points to model deficiencies in describing the sinking and regeneration of 13 C-depleted organic carbon.On the other hand, the δ 13 C DIC bias corresponds to the δ 13 C AS bias in the upper thermocline of the open oceans in the Southern Hemisphere (shown in Fig. 4).Although our results may generally suffer from the coarse model resolution and simplified climate forcing, the specific reasons for this correspondence are not obvious.As a residual term δ 13 C AS may also reflect effects of biological 13 C cycling, which are not captured by Eq. (10).For example, δ 13 C BIO is estimated for constant isotopic fractionation of marine organic matter of −19 ‰ whereas δ 13 C POC varies by about 10 ‰ according to field data (Verwega et al., 2021; see also Fig. 6a).
Therefore, we explore the sensitivity of δ 13 C DIC to photosynthetic fractionation in an additional experiment ("NP") in which biogenic fractionation is disabled (i.e. 13 α p = 1 in Eq. 8).Compared with the default simulation, δ 13 C DIC|NP decreases by up to 1 ‰ in the pelagic euphotic zone whereas δ 13 C DIC|NP increases in the disphotic zone below, which is particularly the case in highly productive regions (Fig. 5a,  b).In the aphotic interior of the ocean δ 13 C DIC|NP progressively increases from the North Atlantic towards the North Pacific by up to 2.4 ‰ (Fig. 5c).These findings are in line with similar sensitivity studies (Schmittner et al., 2013;Dentith et al., 2020) as well as with simulations comparing different parametrizations for 13 α p (Jahn et al., 2015;Buchanan et al., 2019;Dentith et al., 2020;Liu et al., 2021).6a).The model overestimates δ 13 C POC at most locations (Fig. 6b; global RMS difference is 2.6 ‰).According to the sensitivity experiment NP, the overestimation of δ 13 C POC should result in overly enriched δ 13 C DIC in the twilight and dark zones of highly productive regions (Fig. 5b), whereas Fig. 2a indicates that the opposite is the case.However, this is only an apparent contradiction because the δ 13 C POC observations are biased by the 13 C Suess effect.Moreover, they exhibit a negative trend (by −3 ‰ between 1960 and 2010), which is about twice as high as the known 13 C Suess effect on aqueous CO 2 (Verwega et al., 2021).It has been presumed that this trend also reflects a shift in the composition of phytoplankton species (Lorrain et al., 2020;Verwega et al., 2021).Neither effect is considered in our simulation.A conclusive analysis would require transient simulations, including historical values of atmospheric 13 CO 2 .

Radiocarbon
We consider 14 C DIC and compare our model results with gridded fields of pre-bomb 14 C DIC provided by the Global Ocean Data Analysis Project, version 1.1 (GLODAPv1.1,Key et al., 2004).
In the comprehensive radiocarbon cycle simulation (CC) 14 C DIC|CC is within the range −40 ‰ to −140 ‰ (average value: −65 ‰) near the surface (at 50 m), with the highest values in the subtropical gyres and the lowest values in the Southern Ocean (Fig. 7a).In the interior of the Atlantic, 14 C DIC|CC ranges between −70 ‰ and −170 ‰ and decreases from the surface to the bottom, with small vertical gradients in the high latitudes (Fig. 8a).This is superimposed by a southward decline of 14 C DIC values.Conversely, 14 C DIC decreases from south to north in the interior of the Pacific.Different from the northern North Atlantic, there is no evidence of deep-sea ventilation in the North Pacific, where our model arrives at minimum 14 C DIC|CC https://doi.org/10.5194/gmd-17-1709-2024 Geosci.Model Dev., 17, 1709-1727, 2024 values exceeding −290 ‰ (Fig. 8a).This depletion corresponds to a 14 C age τ of about 2800 years with respect to the atmosphere (calculated as τ = −8033 × ln[(1 + 0.001 14 C atm )/(1 + 0.001 14 C DIC )]).
Overall, simulation CC captures the large-scale distribution of pre-nuclear 14 C DIC as reconstructed by b,8a,b).However, at a depth of 50 m the simulated 14 C DIC is too high by 10 ‰-30 ‰ in the lowand mid-latitudes, and too low by about the same amount in the high latitudes (see Fig. 9a; according to GLODAP 14 C DIC ranges from −50 ‰ to −170 ‰ with −71 ‰ on average).In the interior of the oceans, our model results are typically 20 ‰-60 ‰ too low (Fig. 10).The excessive depletion reaches −70 ‰ in the abyssal North Atlantic and in the North Pacific at 3 km depth (Fig. 10a).In the upper layers of the oceans, the GLODAPv1.1 data probably reflect the 14 C Suess effect (Suess, 1955), which is not considered in our simulations.The excessive depletion in the deep sea indicates that the simulated MOC leads to overly shallow and weak ocean ventilation, which is particularly the case in the North Pacific.This effect is enhanced by the progressive ra-dioactive decay of DI 14 C along its passage through the interior of the ocean.
Different from δ 14 C DIC , 14 C DIC is corrected for isotopic fractionation.In practice (as well as in the comprehensive DI 14 C cycle modelling approach CC), the correction is made after the simultaneous determination of δ 13 C DIC and δ 14 C DIC .In the inorganic 14 C (IC) modelling approach the fractionation of 14 C is omitted beforehand, so that posterior corrections are not necessary.That is, δ 14 C DIC|IC equals 14 C DIC|IC , which should equal 14 C DIC|CC .As the IC approach also neglects the radioactive decay of organic carbon, it considers seven tracers less than the CC approach, which is accompanied by an increase in model speed (simulated years per day) of about 15 %.
At a depth of 50 m 14 C DIC|IC ranges from −50 ‰ to −160 ‰ with an average value of 72 ‰ (Fig. 7c).Similar to experiment CC, the highest values of 14 C DIC|IC are found in subtropical surface waters.The lowest surface water values of 14 C DIC|IC are also found in the Southern Ocean.In the interior of the oceans, the isotopic depletion with respect to the atmosphere ranges from −90 ‰ in the North Atlantic to −310 ‰ in the North Pacific (Fig. 8c).Comparing IC with the GLODAPv1.1 reconstruction, we find that the enrichment of subtropical surface values is less pronounced in IC, whereas the depletion in the high latitudes increases (Fig. 9b).The latter is also the case in the interior of the oceans (Fig. 10b).Most notably, the outcomes of experiment IC are everywhere lower (by up to 30 ‰) than the results of simulation CC (Fig. 11a, c) which is explained as follows.
Analogous to sensitivity experiment NP, the IC approach disregards photosynthetic fractionation, which leads to lower DI 14 C concentrations in the euphotic zone than in simulation CC.In addition, the IC approach disregards the DI 14 C enrichment of the mixed layer associated with air-sea exchange.Furthermore, as the IC approach disregards the radioactive decay of phytoplankton, the loss of DI 14 C due to photosynthesis is overestimated in the mixed layer.Therefore, preformed DI 14 C is systematically lower than in simulation CC and becomes further depleted through radioactive decay in the deep sea.This bias is similar to the lower values of "abiotic" 14 C than "biotic" 14 C simulated by Frischknecht et al. (2022).
Instead of computing absolute concentrations of DIC, DI 13 C, and DI 14 C and converting them to 14 C DIC a posteriori, the 14 C approximation (DA) simulates 14 R DIC|DA = 1+0.001 14C DIC|DA as a single prognostic tracer that is considered abiotic and conservative except for its radioactive decay, and that is connected to the carbon cycle only through the timescale of 14 CO 2 air-sea exchange.This approach is about five times faster than the CC and IC approaches.The first results of the implementation of 14 R DIC|DA into FE-SOM2 were shown by Lohmann et al. (2020) for the default FESOM mesh with 127 000 horizontal surface nodes.Here, we repeated the experiment, now using the low-resolution mesh of experiments CC and IC to be able to compare the re-   sults of all approaches directly.For the same reason, we discuss model results after 5000 simulated years but point out that experiment DA has been run over 17 000 years in total.
The maximum drift of 14 C DIC|DA between 5000 and 17 000 simulated years is about −3.5 ‰ in North Pacific Deep Water, which is much smaller than the 14 C DIC differences between the various modelling approaches in this study after 5000 years shown below.
At at depth of 50 m 14 C DIC|DA ranges from −40 ‰ to −130 ‰, with 14 C DIC|DA = −58 ‰ on average (Fig. 7d).In intermediate and deep water 14 C DIC|DA declines from −60 ‰ in the North Atlantic to −280 ‰ in the North Pacific (Fig. 8d).At upper levels, 14 C DIC|DA is almost everywhere higher than 14 C DIC according to GLODAPv1.1 (Figs. 9c, 10c).In the deep sea, 14 C DIC|DA is still too low, but the depletion is less pronounced than in simulations CC and IC (Fig. 10c).
Experiment DA yields the highest 14 C DIC values of all three modelling approaches (Fig. 11).The DA approach assumes that DIC concentrations are constant and homogeneous (for a rigorous treatise see Mouchet, 2013).Following Toggweiler et al. (1989), our calculation of 14 CO 2 air-sea exchange fluxes in simulation DA assumes a DIC concentration of 2000 mmol m −3 in the mixed layer, which is somewhat lower than that observed and simulated in most areas, most notably in the high latitudes (Fig. A5a, b).This leads to faster and hence higher 14 C uptake in DA than in CC and IC be- cause the 14 CO 2 invasion flux is inversely proportional to the DIC concentration in the mixed layer.The absolute 14 C DIC differences between DA and CC are largely less than 10 ‰ (Fig. 11b, d).Moreover, the relative uncertainty of the DA approach with respect to the correct DI 14 C implementation (CC) is less than 5 % (Fig. 12) and actually smaller than the error of 10 % originally estimated by Fiadeiro (1982).

Summary
We have added the carbon isotopes 13 C and 14 C to the marine biogeochemistry model REcoM3 and tested the implementation in long-term equilibrium simulations in which the configuration REcoM3p was coupled with the ocean general circulation model FESOM2.1.Regarding the carbonisotopic composition of DIC (δ 13 C DIC and 14 C DIC ), our model results are largely consistent with marine δ 13 C DIC and 14 C DIC fields reconstructed for the pre-anthropogenic period.The simulations also exhibit discrepancies, such as overly depleted δ 13 C DIC values in upwelling regions, overly steep vertical gradients of δ 13 C DIC and 14 C DIC in the interior of the North Atlantic, and excessive depletion of 14 C DIC in the interior of the North Pacific.To some extent, the inaccuracies of δ 13 C DIC indicate shortcomings in modelled organic carbon cycling.The radiocarbon results ( 14 C DIC ) are insensitive to carbon cycle changes and reflect the rather shallow overturning circulation provided by our low-resolution ocean general circulation model test configuration with idealized repeat year climate forcing.As future simulations with a scientific focus will be carried out with considerably higher horizontal resolution and more realistic climate forcing, we expect some of the biases discussed in this study to decrease.These simulations should also consider transient boundary conditions of 13 C and 14 C, which provide additional benchmarks for the model.For these reasons we did not attempt to further tune REcoM3p here, e.g. by adjusting semi-empirical biogeochemical parameters such as gas transfer velocity or biogenic fractionation coefficients.
As 14 C DIC is dominated by radioactive decay and transport processes, we have additionally explored the accuracy of two simplified modelling approaches that are more efficient than the complete consideration (CC) of the DI 14 C cycle.One approach (IC) neglects isotopic fractionation but still considers biological processes.The IC approach yields lower 14 C DIC values than reconstructed for high-latitude surface waters and for deep and bottom waters.Another approach (DA) only considers the DI 14 C / DIC ratio for constant and homogeneous DIC concentrations and further disregards the marine carbon cycle.It yields higher 14 C DIC values than reconstructed for surface water, which mitigates https://doi.org/10.5194/gmd-17-1709-2024 Geosci.Model Dev., 17, 1709-1727, 2024 the isotopic depletion in the deep sea, where the DA approach therefore better agrees with the reconstruction than the other approaches.The relative uncertainty between the comprehensive and simplified approaches is less than 5 %.Therefore, the simplified 14 C DIC modelling approaches should be sufficiently accurate for radiocarbon dating of marine cliarchives.
Appendix A

Figure 1 .
Figure 1.Preindustrial δ 13 C of DIC, (a, c) this study (simulation CC), (b, d) reconstruction (Eide et al., 2017b).(a, b) values at 200 m depth, (c, d) zonal-mean values in the Atlantic and Pacific.Map projection here and in other figures is area preserving (Equal Earth projection, Šavrič et al., 2019).

Figure 2 .
Figure 2. Differences between simulated (this study; CC) and reconstructed (Eide et al., 2017b; PI) δ 13 C of DIC for the preindustrial period; (a) at a depth of 200 m, (b) zonal-mean values in the Atlantic and Pacific.

Figure 3 .
Figure 3. Differences between simulated (this study; CC) and estimated (Eide et al., 2017b; PI) δ 13 C BIO (the biological component of δ 13 C DIC in the absence of air-sea exchange) during the preindustrial period; (a) at a depth of 200 m, (b) zonal-mean values in the Atlantic and Pacific.

Figure 4 .
Figure 4. Differences between simulated (this study; CC) and estimated (Eide et al., 2017b; PI) δ 13 C AS = δ 13 C DIC − δ 13 C BIO during the preindustrial period; (a) at a depth of 200 m, (b) zonal-mean values in the Atlantic and Pacific.

Figure 5 .
Figure 5. Changes in preindustrial δ 13 C of DIC if isotopic fractionation during photosynthesis is disabled (sensitivity experiment NP versus simulation CC), (a) at a depth of 50 m, (b) at a depth of 200 m, (c) zonal-mean values in the Atlantic and Pacific.

Figure 6 .
Figure 6.(a) δ 13 C of POC averaged over the upper 200 m.Shaded areas: simulation results (this study; simulation CC), dots: bulk matter observations for the period 1964-2015 (compilation by Verwega et al., 2021; see further references therein).(b) Differences between simulated and observed δ 13 C of POC.

Figure 7 .
Figure 7. Preindustrial 14 C of DIC at a depth of 50 m, (a) simulation CC considering the complete marine 14 C cycle, (b) reconstruction (Key et al., 2004), (c) simulation IC applying the inorganic 14 C approximation, (d) simulation DA applying the 14 C approach.See the main text for further explanations.

Figure 8 .
Figure 8. Preindustrial 14 C of DIC in the Atlantic and Pacific.(a) Simulation CC, (b) reconstruction (Key et al., 2004), (c) simulation IC, (d) simulation DA.See the main text for further simulation explanations.Note the different colour scale ranges in Figs.7 and 8.

Figure 9 .
Figure 9. Differences between simulated and reconstructed preindustrial 14 C of DIC at a depth of 50 m, (a) simulation CC minus reconstruction (GLODAP; Key et al., 2004), (b) simulation IC minus reconstruction, (c) simulation DA minus reconstruction.

Figure 10 .
Figure 10.Differences between simulated and reconstructed preindustrial 14 C of DIC in the Atlantic and Pacific: zonal-mean values are shown, (a) simulation CC minus reconstruction (GLODAP; Key et al., 2004), (b) simulation IC minus reconstruction, (c) simulation DA minus reconstruction.Note the different colour scale ranges in Figs. 9 and 10.

Figure 11 .
Figure 11.Absolute differences in preindustrial 14 C DIC between the various simulation approaches: results at a depth of 50 m (a, b) and in the Atlantic and Pacific (c, d) are shown.(a, c) Inorganic 14 C (IC) modelling approach versus complete 14 C cycle (CC), (b, d) 14 C approximation (DA) versus simulation CC.

Figure 12 .
Figure 12.Relative differences in preindustrial 14 C DIC between the various simulation approaches: absolute values at a depth of 50 m depth (a, b) and in the Atlantic and Pacific (c, d) are shown.(a, c) No-fractionation (IC) approach versus complete 14 C cycle (CC), (b, d) 14 C approximation (DA) versus simulation CC.

Figure A2 .
Figure A2.Differences between simulated and observed (WOA09, Locarnini et al., 2010) temperatures, (a) averaged over the upper 200 m, (b) in the Atlantic and Pacific.See Scholz et al. (2019) for comparison with higher-resolution simulations.

Figure A3 .
Figure A3.Differences between simulated and observed (WOA09, Antonov et al., 2010) salinities, (a) averaged over the upper 200 m, (b) in the Atlantic and Pacific.See Scholz et al. (2019) for comparison with higher-resolution simulations.

Figure A6 .
Figure A6.Preindustrial δ 13 C DIC of surface water at a depth of about 18 m.Shaded areas: simulation CC; filled circles: reconstructed values by Kwon et al. (2022).Note that their reconstruction yields higher preindustrial δ 13 C DIC values than Eide et al. (2017b) in thermocline and intermediate waters of the Southern Hemisphere (see the discussion by Kwon et al., 2022).

Table 1 .
List of model experiments discussed in this paper.