2.1 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 REcoM3 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.

2.2 Implementation of carbon isotopes

Carbon-13 and ¹⁴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 ¹³C and ¹⁴C are small (¹³R_std≈0.0113, Assonov et al., 2020, and ${}^{\mathrm{14}}{R}_{\mathrm{std}}=\mathrm{1.170}\times {\mathrm{10}}^{-\mathrm{12}}$; Orr et al., 2017), we approximate total C with ¹²C and neglect ¹³C and ¹⁴C in the kinetic calculations of the carbonate system. For the same reason, and to ensure numerical stability in the other model calculations involving ¹³C and ¹⁴C, the model considers scaled concentrations ¹³C${}^{\prime }={}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{13}}{R}_{\mathrm{std}}$ and ¹⁴C${}^{\prime }={}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{14}}{R}_{\mathrm{std}}$, which have the same magnitude as ¹²C, following Jahn et al. (2015). As a consequence, ¹³C^′ and ¹⁴C^′ have to be rescaled when compared with observed concentrations. However, as the scaling factors cancel out when ¹³C^′ and ¹⁴C^′ are converted to δ¹³C and Δ¹⁴C, Eqs. (1)–(3) can be directly evaluated with scaled concentration ratios and ${}^{\mathrm{13}}{R}_{\mathrm{std}}{=}^{\mathrm{14}}{R}_{\mathrm{std}}=\mathrm{1}$.

We consider isotopic fractionation during air–sea gas exchange, dissolution of CO₂ in seawater, and photosynthesis by phytoplankton. In addition, the model accounts for radioactive decay and, optionally, cosmogenic production of ¹⁴C. The details are explained in the following subsections.

2.3 Carbon-13

2.4 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 ¹⁴C production fluxes or, alternatively, by prescribed atmospheric ¹⁴CO₂ concentrations corresponding to atmospheric Δ¹⁴C values.

Fractionation factors are calculated analogously to ¹³α, with ${}^{\mathrm{14}}\mathit{\alpha }={\mathrm{2}}^{\mathrm{13}}\mathit{\alpha }-\mathrm{1}$ (e.g. Craig, 1954; see also the discussion by Fahrni et al., 2017). The slow equilibration of DI¹⁴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 ¹⁴C in various ways that differ in their level of complexity and computational efficiency. The first implementation (“CC”) considers the complete ¹⁴C cycle parallel to ¹³C. The second implementation (“IC”) disregards the isotopic fractionation of ¹⁴C and radioactive decay of organic matter. In turn, DI¹⁴C and DIC have identical sources and sinks, except for atmospheric CO₂ and radioactive decay. This “inorganic” ¹⁴C approach is an approximation that is conceptually similar, but not identical, to the “abiotic” ¹⁴C modelling approach described in the OMIP biogeochemical protocol (Orr et al., 2017). In our IC approach, DIC and DI¹⁴C include biological sources and sinks. This does not apply to the OMIP-abiotic 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 Δ¹⁴C values of DIC (Δ¹⁴C_DIC) are primarily governed by transport and radioactive decay (Fiadeiro, 1982; see also Mouchet, 2013). This implies that Δ¹⁴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 Δ¹⁴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.

Table 1List of model experiments discussed in this paper.

Download Print Version | Download XLSX

3.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 marine 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 FESOM2.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² s⁻¹ originally to 1000 m² s⁻¹. 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⁻⁶ K per year and 10⁻⁶ PSU per year in the global average, and the major biogeochemical tracer inventories drifted by less than 10⁻³ % per year (less than 10⁻³ ‰ per year regarding DI¹³C and DI¹⁴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 ¹³C and ¹⁴C were initialized with DIC, assuming initial fractionation values of δ¹³C_DIC=0 ‰ and Δ¹⁴C${}_{\mathrm{DIC}}=-\mathrm{150}$ ‰ respectively. In the sediment layer all initial concentrations were close to zero. REcoM3p was forced with constant atmospheric CO₂ concentrations ([¹²CO₂]_atm=284.3 ppmv; [¹³CO₂]_atm and [¹⁴CO₂]_atm corresponding to δ¹³C${}_{\mathrm{atm}}=-\mathrm{6.61}$ ‰ and Δ¹⁴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₂ 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, 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⁶ m³ s⁻¹, 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).

Figure 1Preindustrial δ¹³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).

3.2 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 δ¹³C_DIC, we compare our model results with the gridded preindustrial (PI) δ¹³C_DIC climatology by Eide et al. (2017b). Their reconstruction does not consider the upper 200 m to exclude the ¹³C Suess effect. For the sake of completeness, we briefly note that at the sea surface our δ¹³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 δ¹³C_DIC values at 200 m are within the range 1 ‰ to 2 ‰ (Fig. 1a). Higher δ¹³C_DIC values are simulated in the subtropical oceans, particularly in the Atlantic, southeast Pacific, and southern Indian Ocean. Isotopic depletion of up to $\sim -\mathrm{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 δ¹³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).

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

Considering the interior of the Atlantic Ocean, our simulation yields higher δ¹³C_DIC in the North Atlantic than in the South Atlantic, and small vertical δ¹³C_DIC gradients in the high latitudes of both hemispheres (Fig. 1c). In the Pacific, our simulation displays a reversed meridional δ¹³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 δ¹³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 models (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 δ¹³C_DIC (the δ¹³C_DIC bias for short) can be attributed to various factors. First, the reconstructed δ¹³C_DIC values may be too low in the upper ocean if the ¹³C Suess effect is not completely removed from the observations. In fact, Eide et al. (2017a) presume that they may have underestimated the ¹³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 δ¹³C_DIC bias can be further attributed by decomposing δ¹³C_DIC into biologic and thermodynamic sources:

where δ¹³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), δ¹³C_BIO is frequently determined from the covariation of δ¹³C_DIC with marine phosphate (PO${}_{\mathrm{4}}^{\mathrm{3}-}$). We adopt this approach, employing revised parameter values by Eide et al. (2017b) and tentatively replacing PO${}_{\mathrm{4}}^{\mathrm{3}-}$ with dissolved inorganic nitrogen (DIN) divided by 16 because PO${}_{\mathrm{4}}^{\mathrm{3}-}$ is not considered by REcoM3p:

where PO${}_{\mathrm{4}}^{\mathrm{3}-}$ and DIN are in µmol kg⁻¹. Equation (10) does not distinguish between preformed and remineralized nutrient concentrations and cannot be used for further separation of δ¹³C_BIO into preformed and remineralized components. The thermodynamic component δ¹³C_AS describes the effects of air–sea exchange and ocean circulation and is the residual of observed δ¹³C_DIC and reconstructed δ¹³C_BIO.

Figure 3Differences between simulated (this study; CC) and estimated (Eide et al., 2017b; PI) δ¹³C_BIO (the biological component of δ¹³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 4Differences between simulated (this study; CC) and estimated (Eide et al., 2017b; PI) δ¹³C_AS=δ¹³C_DIC−δ¹³C_BIO during the preindustrial period; (a) at a depth of 200 m, (b) zonal-mean values in the Atlantic and Pacific.

In the following we validate simulated δ¹³C_BIO and δ¹³C_AS with the corresponding values reconstructed by Eide et al. (2017b), and we compare the δ¹³C_DIC bias with the differences between the simulated and reconstructed δ¹³C_DIC components, i.e. the δ¹³C_BIO bias and δ¹³C_AS bias. Simulation CC overestimates δ¹³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 δ¹³C_AS in the thermocline but underestimates δ¹³C_AS in the interior of the North Atlantic above about 3 km (Fig. 4). Comparing the biases of δ¹³C_DIC, δ¹³C_BIO and δ¹³C_AS, it appears that the δ¹³C_DIC bias corresponds to the δ¹³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 ¹³C-depleted organic carbon. On the other hand, the δ¹³C_DIC bias corresponds to the δ¹³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 δ¹³C_AS may also reflect effects of biological ¹³C cycling, which are not captured by Eq. (10). For example, δ¹³C_BIO is estimated for constant isotopic fractionation of marine organic matter of −19 ‰ whereas δ¹³C_POC varies by about 10 ‰ according to field data (Verwega et al., 2021; see also Fig. 6a).

Figure 5Changes in preindustrial δ¹³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(a) δ¹³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 δ¹³C of POC.

Therefore, we explore the sensitivity of δ¹³C_DIC to photosynthetic fractionation in an additional experiment (“NP”) in which biogenic fractionation is disabled (i.e. ¹³α_p=1 in Eq. 8). Compared with the default simulation, δ¹³C_DIC|NP decreases by up to 1 ‰ in the pelagic euphotic zone whereas δ¹³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 δ¹³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 ¹³α_p (Jahn et al., 2015; Buchanan et al., 2019; Dentith et al., 2020; Liu et al., 2021).

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

Figure 8Preindustrial Δ¹⁴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.

Download

Our default simulation with enabled photosynthetic fractionation yields δ¹³C_POC values between −18 ‰ and −25 ‰ whereas observations from the last decades range from −15 ‰ to −35 ‰ (Verwega et al., 2021; Fig. 6a). The model overestimates δ¹³C_POC at most locations (Fig. 6b; global RMS difference is 2.6 ‰). According to the sensitivity experiment NP, the overestimation of δ¹³C_POC should result in overly enriched δ¹³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 δ¹³C_POC observations are biased by the ¹³C Suess effect. Moreover, they exhibit a negative trend (by −3 ‰ between 1960 and 2010), which is about twice as high as the known ¹³C Suess effect on aqueous CO₂ (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 ¹³CO₂.

3.3 Radiocarbon

We consider Δ¹⁴C_DIC and compare our model results with gridded fields of pre-bomb Δ¹⁴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) Δ¹⁴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, Δ¹⁴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 Δ¹⁴C_DIC values. Conversely, Δ¹⁴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 Δ¹⁴C_DIC|CC values exceeding −290 ‰ (Fig. 8a). This depletion corresponds to a ¹⁴C age τ of about 2800 years with respect to the atmosphere (calculated as $\mathit{\tau }=-\mathrm{8033}\times \ln \left[\right(\mathrm{1}+\mathrm{0.001}{\mathrm{\Delta }}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{atm}})/(\mathrm{1}+\mathrm{0.001}{\mathrm{\Delta }}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{DIC}}\left)\right]$).

Figure 9Differences between simulated and reconstructed preindustrial Δ¹⁴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 10Differences between simulated and reconstructed preindustrial Δ¹⁴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.

Download

Overall, simulation CC captures the large-scale distribution of pre-nuclear Δ¹⁴C_DIC as reconstructed by GLODAPv1.1 (Figs. 7a, b, 8a, b). However, at a depth of 50 m the simulated Δ¹⁴C_DIC is too high by 10 ‰–30 ‰ in the low- and mid-latitudes, and too low by about the same amount in the high latitudes (see Fig. 9a; according to GLODAP Δ¹⁴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 ¹⁴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 radioactive decay of DI¹⁴C along its passage through the interior of the ocean.

Different from δ¹⁴C_DIC, Δ¹⁴C_DIC is corrected for isotopic fractionation. In practice (as well as in the comprehensive DI¹⁴C cycle modelling approach CC), the correction is made after the simultaneous determination of δ¹³C_DIC and δ¹⁴C_DIC. In the inorganic ¹⁴C (IC) modelling approach the fractionation of ¹⁴C is omitted beforehand, so that posterior corrections are not necessary. That is, δ¹⁴C_DIC|IC equals Δ¹⁴C_DIC|IC, which should equal Δ¹⁴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 Δ¹⁴C_DIC|IC ranges from −50 ‰ to −160 ‰ with an average value of 72 ‰ (Fig. 7c). Similar to experiment CC, the highest values of Δ¹⁴C_DIC|IC are found in subtropical surface waters. The lowest surface water values of Δ¹⁴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.

Figure 11Absolute differences in preindustrial Δ¹⁴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 ¹⁴C (IC) modelling approach versus complete ¹⁴C cycle (CC), (b, d) Δ¹⁴C approximation (DA) versus simulation CC.

Analogous to sensitivity experiment NP, the IC approach disregards photosynthetic fractionation, which leads to lower DI¹⁴C concentrations in the euphotic zone than in simulation CC. In addition, the IC approach disregards the DI¹⁴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¹⁴C due to photosynthesis is overestimated in the mixed layer. Therefore, preformed DI¹⁴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” Δ¹⁴C than “biotic” Δ¹⁴C simulated by Frischknecht et al. (2022).

Instead of computing absolute concentrations of DIC, DI¹³C, and DI¹⁴C and converting them to Δ¹⁴C_DIC a posteriori, the Δ¹⁴C approximation (DA) simulates ${}^{\mathrm{14}}{R}_{\mathrm{DIC}\mathrm{|}\mathrm{DA}}=\mathrm{1}+\mathrm{0.001}{\mathrm{\Delta }}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{DIC}\mathrm{|}\mathrm{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 ¹⁴CO₂ air–sea exchange. This approach is about five times faster than the CC and IC approaches. The first results of the implementation of ¹⁴R_DIC|DA into FESOM2 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 results 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 Δ¹⁴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 Δ¹⁴C_DIC differences between the various modelling approaches in this study after 5000 years shown below.

At at depth of 50 m Δ¹⁴C_DIC|DA ranges from −40 ‰ to −130 ‰, with Δ¹⁴C${}_{\mathrm{DIC}\mathrm{|}\mathrm{DA}}=-\mathrm{58}$ ‰ on average (Fig. 7d). In intermediate and deep water Δ¹⁴C_DIC|DA declines from −60 ‰ in the North Atlantic to −280 ‰ in the North Pacific (Fig. 8d). At upper levels, Δ¹⁴C_DIC|DA is almost everywhere higher than Δ¹⁴C_DIC according to GLODAPv1.1 (Figs. 9c, 10c). In the deep sea, Δ¹⁴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 Δ¹⁴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 ¹⁴CO₂ air–sea exchange fluxes in simulation DA assumes a DIC concentration of 2000 mmol m⁻³ 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 ¹⁴C uptake in DA than in CC and IC because the ¹⁴CO₂ invasion flux is inversely proportional to the DIC concentration in the mixed layer. The absolute Δ¹⁴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¹⁴C implementation (CC) is less than 5 % (Fig. 12) and actually smaller than the error of 10 % originally estimated by Fiadeiro (1982).

Figure 12Relative differences in preindustrial Δ¹⁴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 ¹⁴C cycle (CC), (b, d) Δ¹⁴C approximation (DA) versus simulation CC.