Carbon isotopes in the ocean model of the Community Earth System Model ( CESM 1 )

Carbon isotopes in the ocean are frequently used as paleoclimate proxies and as present-day geochemical ocean tracers. In order to allow a more direct comparison of climate model results with this large and currently underutilized data set, we added a carbon isotope module to the ocean model of the Community Earth System Model (CESM), containing the cycling of the stable isotope C and the radioactive isotope C. We implemented the C tracer in two ways: in the “abiotic” case, the C tracer is only subject to air– sea gas exchange, physical transport, and radioactive decay, while in the “biotic” version, the C additionally follows the C tracer through all biogeochemical and ecological processes. Thus, the abiotic C tracer can be run without the ecosystem module, requiring significantly fewer computational resources. The carbon isotope module calculates the carbon isotopic fractionation during gas exchange, photosynthesis, and calcium carbonate formation, while any subsequent biological process such as remineralization as well as any external inputs are assumed to occur without fractionation. Given the uncertainty associated with the biological fractionation during photosynthesis, we implemented and tested three parameterizations of different complexity. Compared to present-day observations, the model is able to simulate the oceanic C bomb uptake and the C Suess effect reasonably well compared to observations and other model studies. At the same time, the carbon isotopes reveal biases in the physical model, for example, too sluggish ventilation of the deep Pacific Ocean.


Introduction
A large fraction of paleoclimatic reconstructions are based on isotopic measurements (e.g., Petit et al., 1999;McDermott, 2004;Curry and Oppo, 2005;Polka et al., 2013), yet there are many uncertainties associated with the interpretation of these records in terms of physical climate variables such as temperature, precipitation, and ocean circulation rates.More direct comparisons of paleo data with climate models would therefore be beneficial, both to test the interpretation of the isotopic proxy data and to allow for better comparisons of model simulations with proxy data.Furthermore, many isotope tracers are currently being measured in the ocean, and including them in ocean models can help us better understand the ocean circulation and diagnose model biases (e.g., Matsumoto et al., 2004).For all of these reasons, we have added a carbon isotope module to the ocean model of the Community Earth System Model (CESM) (Hurrell et al., 2013).
Carbon has two stable isotopes, 12 C and 13 C.More than 98.9 % of carbon on earth is 12 C, while 13 C makes up most Published by Copernicus Publications on behalf of the European Geosciences Union.
of the remaining 1 %.The radioactive carbon isotope 14 C, also called radiocarbon, is present only in trace amounts (approximately 1 × 10 −10 % of all carbon) and has a half-life of 5730 years (Godwin, 1962).Radiocarbon is a useful tracer to evaluate the ventilation of the deep ocean because it acts as a clock, measuring the time since water was last in contact with the atmosphere (e.g., Toggweiler et al., 1989;Orr, 2002;Meissner et al., 2003;Waugh et al., 2003;Key et al., 2004;Doney et al., 2004;Matsumoto et al., 2004;Meissner, 2007;Bardin et al., 2014).Because of the atmospheric nuclear weapons tests in the 1950s and 1960s and the wellknown input function of radiocarbon during this time, radiocarbon is also useful to evaluate the recent penetration of anthropogenic carbon into the ocean (e.g., Graven et al., 2012).Furthermore, oceanic radiocarbon has been used to determine the mean gas exchange velocity used in ocean models (e.g., Wanninkhof, 1992;Sweeney et al., 2007;Naegler et al., 2006;Naegler, 2009).Oceanic δ 13 C has been used in paleoclimate studies as a tracer of the ocean circulation (e.g., Marchal et al., 1998;Curry and Oppo, 2005;Crucifix, 2005), to calculate the uptake of anthropogenic carbon dioxide (e.g., Keeling et al., 1980;Quay et al., 1992;Gruber et al., 1999;Sonnerup et al., 1999;Gruber and Keeling, 2001), and to diagnose biases in marine ecosystem models (e.g., Schmittner et al., 2013).
We added the carbon isotopes to the code so that they follow the cycling of total carbon through all ecosystem and physical/chemical processes.In this biotic formulation, a new 13 C and 14 C state variable was added to each carbonbearing state variable, resulting in a total of 14 new state variables.For 14 C, we also added the option of a simplified representation, where the isotope is only subject to the main chemical and physical processes during gas exchange and decay, but does not cycle through the ecosystem.This abiotic formulation of 14 C was implemented based on the Ocean Carbon Model Intercomparison Project Phase 2 (OCMIP-2) protocol (Orr et al., 2000).
Abiotic radiocarbon had been added previously to the NCAR ocean model (in NCOM1.4, Orr, 2002, andPOP1/CCSM3, Graven et al., 2012), and biotic 13 C was implemented in the ecosystem model of the CCSM3 by X. Giraud and N. Gruber in 2009-2010.However, neither development was ever added to the trunk of the ocean model code of the CESM, so it was not maintained as the model evolved over the years and consequently none of these implementations still work in the current ocean model of the CESM.In contrast, the addition of a biotic radiocarbon tracer is completely new in this implementation in the CESM.In order to increase the chances of maintaining these developments as the model continues to evolve, the current implementation has been added to the code trunk of the current ocean model of CESM.By including carbon isotopes in the ocean model of the CESM1, the CESM1 joins the community of other comprehensive ocean general circulation models that include abiotic radiocarbon and/or biotic 13 C in the ocean (e.g., Mo-BidiC, Crucifix, 2005; PISCES, Tagliabue and Bopp, 2008;CM2Mc ESM, Galbraith et al., 2011;, HAMOCC2s, Hesse et al., 2011;and UVic ESCM, e.g., Meissner et al., 2003;Schmittner et al., 2013).While the abiotic radiocarbon implementation tends to follow the OCMIP-2 protocol (Orr et al., 2000) in all models, the implementations of biotic 13 C differ between models, mainly due to the complexity of the ocean biogeochemistry model used in them, but also due to different choices in regards to the parameterization of the biological fractionation during photosynthesis and calcium carbonate (CaCO 3 ) formation.
As a reference for future studies using these new capabilities in the CESM, we describe the model used (Sect.2), describe the details of the implementation of the abiotic and biotic carbon isotopes (Sect.3), and compare the simulated carbon-isotope fields to observational data to show the general performance of the model (Sect.4).

Model
This work was done using the code base of the Community Earth System Model (CESM) (Hurrell et al., 2013), version 1.0.5.The isotope code has been updated to the current version of the CESM and is targeted for public release in 2016 as part of CESM2 (see the section on code availability at the end of the article).The CESM is a fully coupled climate model with components for the atmosphere, land, river runoff, sea ice, ocean and ice sheets, coupled by a coupler.Its components and simulations have been described in a large collection of articles, many of them contained in a special collection in the Journal of Climate (http://journals.ametsoc.org/page/CCSM4/CESM1).The simulations analyzed here were performed using the ocean model coupled to data models for the atmosphere, the land, the sea ice, and the river routing, using repeated normal year forcing from CORE-II (Large and Yeager, 2009).The ocean model was run at a nominal 3 • horizontal resolution with 60 vertical levels, which is the lowresolution configuration of the ocean model (Shields et al., 2012).This ocean-only model version with ocean biogeochemistry at 3 • model resolution is used as a low-cost test bed for model development, but is not scientifically validated or supported.Hence, future science applications of the carbon isotopes should use the scientifically supported 1 • horizontal resolution model version of the CESM.

Carbon isotope implementation
The carbon isotopes were added as optional passive tracers, with the biotic and abiotic implementations as two different options that can be set at the compilation and build time.The abiotic 14 C can be run with or without the ocean ecosystem model, while the biotic 13 C and 14 C require the ocean ecosystem model to be turned on.

Abiotic 14 C
In this implementation, DI 14 C is the model's normalized concentration of total dissolved inorganic 14 C, following the OCMIP2 protocol (Orr et al., 2000).DI 14 C is used as normalized concentration in order to minimize the numerical error of carrying very small numbers.The normalization is done by dividing the real DI 14 C by the standard ratio of 14 C / 12 C = 1.176 × 10 −12 (Karlen et al., 1968).To obtain comparable DI 14 C values as measured, we multiply the simulated DI 14 C by this scaling factor of 1.176 × 10 −12 .Since the abiotic radiocarbon is designed to be run without the ocean ecosystem active, we also carry an abiotic DI 12 C tracer to calculate the isotope ratio 14 R = DI 14 C / DI 12 C.For comparisons with observations, we calculate 14 C as a diagnostic variable: (1) By construction, the abiotic DI 12 C and DI 14 C tracers only depend on the solubility of carbon in seawater and neglect all biological activity.The error in 14 C due to neglecting biology activity has been estimated to be on the order of 10 % (Fiadiero, 1982).
Note that we do not multiply 14 R by 14 R std in Eq. ( 1), as we are using a normalized DI 14 C (following Orr et al., 2000).Given that this abiotic implementation does not account for the fractionation during gas exchange, we do not apply the correction for fractionation that is commonly applied to observational measurements of 14 C / 12 C ratios (as well as for the biotic 14 C implementation; see Eq. ( 27) in Sect.3.2.4).The simulated abiotic 14 C is therefore directly comparable to observed data reported as 14 C (see Toggweiler et al., 1989, for more details).

Surface fluxes
We follow the abiotic OCMIP-2 protocol (Orr et al., 2000) for most of the implementation of the abiotic radiocarbon surface fluxes, with the following notable differences.
-We use a coefficient a of 0.31 cm h −1 (Wanninkhof, 1992) instead of 0.337 cm h −1 as used in OCMIP-2.This is higher than what most recent estimates suggest (e.g., Sweeney et al., 2007;Naegler et al., 2006;Naegler, 2009;Graven et al., 2012), but makes it consistent with the gas-transfer formulation used in other parts of the CESM.
-We use the daily mean of the squared 10 m wind speed (either from the prescribed CORE-II forcing or from the coupled atmospheric model) instead of the climatology of the squared monthly average of the instantaneous SSMI velocity and its instantaneous variance as used in OCMIP-2.
-We use the daily mean of the ice fraction and atmospheric pressure (either from the data models or the coupled sea ice and atmosphere models) instead of the monthly averaged climatology used in OCMIP-2.
-We use a constant reference value (1944 µmol m −3 ) for the virtual fluxes of abiotic radiocarbon, rather than an annually updated average of the surface DI 14 C as suggested in OCMIP-2.This is done to conserve total 14 C in the model (in the absence of radioactive decay).
To compute the partial pressure of CO 2 from the abiotic DI 12 C, we require an estimate of surface alkalinity.We follow again OCMIP-2, i.e., we estimate surface alkalinity (Alk) by scaling the ocean mean alkalinity, Alkbar = 2310 microeq kg −1 with sea-surface salinity, SSS, i.e., with S Ref = 34.7 and ρ sw = 4.1/3.996g cm −3 (these two are constants in the CESM).We alter this calculation in the Baltic Sea and the Black Sea to avoid unrealistic alkalinity values, following the procedure developed by K. Lindsay for creating initial conditions for the marine ecosystem model: in the Black Sea, the surface alkalinity is independent of SSS: alkalinity = 3300 • ρ sw .In the Baltic Sea, we calculate alkalinity depending on the surface salinity, with alkalinity = 119 + 196 • SSS when SSS is equal to or below 7.3, and alkalinity = 1237 + 43 • SSS when the SSS is above 7.3.The computation of pCO 2 also requires an assumption about the surface ocean concentrations of silicic acid and phosphate, for which we use OCMIP-2's global constants, i.e., 7.5 µmol kg −1 for silicic acid, Si(OH) 4 , and 0.5 µmol kg −1 for phosphate, PO 4 .

Air-sea gas exchange
As in OCMIP-2, the air-sea gas exchange flux of 12 C is calculated as with PV being the CO 2 gas transfer velocity (called the piston velocity) in m s −1 , calculated as The coefficient a is taken as 0.31 cm h −1 as mentioned earlier, aice is the fraction of the ocean covered by sea ice, u 2 10 is the squared 10 m wind speed from the coupler, and Sc CO 2 is the Schmidt number of CO 2 .Sc CO 2 is calculated as in the ecosystem model, following Wanninkhof (1992): C surf in the gas flux calculation above is the surface aqueous CO 2 concentration in mol m −3 (also called CO * 2 , which is the aqueous CO 2 concentration in mol m −3 in the ocean in general).C sat is the saturation concentration in mol m −3 , with C sat = CO * 2 + DCO * 2 and DCO * 2 being the difference in CO 2 concentration between the surface ocean and the atmosphere.SST is the sea surface temperature.CO * 2 and DCO * 2 in turn are calculated by the carbonate solver from the ecosystem model, based on SST, SSS, ALK, PO 4 , Si(OH) 4 , pH, atmospheric pCO 2 , atmospheric pressure, and the abiotic DI 12 C and DI 14 C concentrations in the surface water.
As in OCMIP-2, we do not account for fractionation during gas exchange in this abiotic formulation, as the effect of isotopic fractionation is almost completely accounted for by the standard correction made when calculating 14 C from observations (see Toggweiler et al., 1989, for details).
The gas flux of the normalized abiotic DI 14 C is calculated as with and The values of the atmospheric pCO 2 and 14 C atm can be set to be constants or can be read in from a file.For atmospheric pCO 2 , it can also be taken from the coupler, to ensure the use of a consistent atmospheric pCO 2 value across model components.Currently the code is set up to read in three files of 14 C atm values, one each for the Northern Hemisphere, the equatorial region (20 • N-20 • S), and the Southern Hemisphere, in order to represent the spatial inhomogeneity of 14 C atm , for example, after the atmospheric nuclear bomb tests.

Virtual fluxes
The CESM ocean model is a volume-conserving model where water fluxes at the surface (from precipitation, evaporation, and river input) are added as virtual fluxes.These virtual fluxes represent the dilution or concentration effect from adding or removing freshwater.For the abiotic carbon isotope tracers, we have a virtual DI 12 C and DI 14 C flux.As for salinity and for DIC in the ecosystem model, we use a constant surface reference DI 12 C and DI 14 C for the calculation of virtual fluxes in order to conserve tracers.The reference values are 1944 µmol m −3 for both DI 12 C and normalized DI 14 C, the same as for DIC in the ecosystem model of CESM.

Interior processes
In the interior of the ocean, the only additional term to the transport of the tracers by the physical ocean model is the decay term for DI 14 C, following the OCMIP-2 protocol: with L being the 3-D transport operator and λ being the radioactive decay constant for 14 C in s −1 , using a half-life of 5730 years (Godwin, 1962): The radiocarbon age (relative to AD 1950 = 0 yr BP) is calculated from 14 C following 14 C age = −5730/ ln 2 × ln(1 + 14 C/1000). ( 5730 years / ln 2 = 8267 years is the mean life of 14 C, which differs from the often used mean life of 8033 years (e.g., Stuiver and Polach, 1977), which is based on the earlier Libby half-life of 5568 (Libby, 1955).

Biotic 13 C and 14 C
In the biotic implementation of 13 C and 14 C, we use the ocean ecosystem model (e.g., Moore et al., 2013) to compute the carbon pools as well as all other biological variables (like silicic acid and alkalinity).The ecosystem model currently has seven carbon pools: DIC, DOC (dissolved organic carbon), CaCO 3 , diazotrophs, diatoms, small phytoplankton, and zooplankton.We carry passive tracers for each of these in the isotope-enabled version of the code.As 12 C makes up over 98 % of the carbon earth and does not fractionate, we assume that the ecosystem carries 12 C.This means that the isotope ratio R can be calculated as the ratio of the new isotopic carbon pools to the ecosystem carbon pools.As for the abiotic radiocarbon, we use scaled variables for 13 C and 14 C in order to minimize the numerical error of carrying very small numbers (particularly for 14 C).The scaling factor is the commonly used standard iso C / 12 C for each isotope, i.e., 1.12372 × 10 −8 for iso = 13 C (Craig, 1957) and 1.176×10 −12 for iso = 14 C (Karlen et al., 1968).This means that we use 13 R Std = 1 and 14 R Std = 1 in the code, and that the model-simulated isotopic carbon pools are multiplied by the respective scaling factor to compare them with observations.
In the biotic formulation, we account for the fractionation of 13 C and 14 C during gas exchange and during biological processes.The fractionation ( ) of 14 C is always twice that of 13 C, as all relevant processes have a mass-dependent fractionation for carbon (Bigeleisen, 1952;Craig, 1954).The isotopic fractionation is related to the fractionation factor α through As diagnostic variables, we compute the δ iso C values by first computing the ratio iso R = DI iso C / DIC and then using As for the abiotic 14 C calculation in Eq. ( 1), we do not multiply by iso R Std in the calculation of δ iso C because we are using normalized DI iso C.

Air-sea gas exchange of 13 C
The air-sea flux of13 C is calculated based on Zhang et al. (1995): Here, C sat and C surf are obtained from the ecosystem model.α k = −0.99919 is the constant kinetic fractionation factor from Zhang et al. (1995) (with = −0.81 and α = /1000+1).α aq g is the temperature (TEMP, in • C) dependent isotopic fractionation factor during gas dissolution, based on the equation for aq g from Zhang et al. (1995): The temperature and carbonate fraction (f CO 3 ) dependent fractionation factor (α DIC g ) between total DIC and CO 2 is based on the empirical relationship for DIC g from Zhang et al. (1995): R 13 C atm is the 13 C to 12 C ratio in atmospheric CO 2 , calculated using the atmospheric δ 13 C atm record and R atm = 1 + δ 13 C atm /1000 (scaled by 13 R Std ).The values of δ 13 C atm can be set to be a constant or it can be read in from a file.Currently δ 13 C atm is assumed to be well mixed globally, so only one global value is read in.With small code modifications globally inhomogeneous δ 13 C atm values can easily be read in instead.R 13 C DIC is the 13 C to 12 C ratio of dissolved inorganic carbon, calculated from the simulated biotic DIC and DI 13 C.

Virtual fluxes of 13 C
As stated in Sect.3.1, we account for the dilution and concentration effect of surface freshwater fluxes in the model by adding a virtual flux, using a constant surface reference DI 13 C (and DI 14 C) of 1944 µmol m −3 for the calculation of virtual fluxes.

Biological fractionation of 13 C
The isotopic carbon fixation by photosynthesis (photo 13 C) is computed from the 12 C fixation during photosynthesis (pho-toC, from the ecosystem model), using Table 1.Parameters used in the parameterization of p for the implementation following Keller and Morel (1999).The values for small phytoplankton are based on E. huxleyi, the value for diatoms are based on P. tricornumtum, and the values for diatoms are based on Synechococcus sp.(Keller and Morel, 1999;Popp et al., 1998).and The strength of the biological fractionation of carbon during photosynthesis ( p ), as well as the key controlling parameters, are still being debated in the literature (e.g., Keller and Morel, 1999), and many of the existing 13 C implementations in models use different parameterizations.We therefore implemented three different parameterizations for p to test the sensitivity of our results to the choice of biological fractionation.
The simplest model for p by Rau et al. (1989) gives the same p value for all types of autotrophs: This relationship is based on the empirical relationship found by Rau et al. (1989) between the isotopic composition of the autotroph (δ C p ) and CO * 2 : limiting δ C p to values between −18 and −32 ‰ (Rau et al., 1989).Laws et al. (1995) assumed that CO 2 enters the cell by diffusion and that the fractionation depends on the rate of photosynthesis, and therefore parameterized p as a function of CO * 2 and the specific photosynthesis rate of each phytoplankton group (µ, in s −1 , calculated by the ecosystem model): Keller and Morel (1999) argued that only considering diffusive CO 2 transport into cells and assuming a linear relationship between p and CO * 2 concentration and the specific growth rate (µ) does not agree with laboratory and field data, citing work by Sikes et al. (1980), Tortell et al. (1997), and Laws et al. (1997).Keller and Morel (1999) therefore proposed to use phytoplankton type-specific (constant) cell parameters (see Table 1) to compute the fractionation during photosynthesis: where with Qc being the cell carbon content, cell permea being the cell wall permeability to CO 2 (aq), cell surf being the surface areas of cells, C up being the ratio of active carbon uptake to carbon fixation, fix being a constant phytoplankton typedependent fractionation effect of carbon fixation, diff = 0.7 representing the fractionation by diffusion (O' Leary, 1984), and δ d13C = −9.0being the difference between the isotopic compositions of the external CO 2 and the organic matter pools (Goericke et al., 1994).
While the fractionation during calcium carbonate formation is much smaller than the fractionation during photosynthesis (Turner, 1982), we include a small constant fractionation of 2 ‰ for calcium carbonate formation, based on work by Ziveri et al. (2003) that found a range of 3 ‰ to −2 ‰ for different species.Other implementations of 13 C in ocean models have used values of 1 ‰ (e.g., Sonnerup et al., 1999;Tagliabue and Bopp, 2008) or have assumed no isotopic fractionation for calcification (e.g., Marchal et al., 1998;Schmittner et al., 2013).However, as shown by Schmittner et al. (2013), the effect of the calcium carbonate pump on δ 13 C is small, so the choice of the value for the fractionation during calcium carbonate formation is not expected to have a big impact on the results in the current ecosystem model with one species of calcium carbonate.

Biotic 14 C
The 14 C air-sea flux is calculated in the same way as shown in Eq. (15) for 13 C, but with the fractionation for 14 C being twice as large as for 13 C ( 14 = 2 • 13 , Zeebe and Wolf-Gladrow, 2001) and with R 14 C atm and R 14 C DIC instead of R 13 C atm and R 13 C DIC .The biological fractionation is also the same as for 13 C, except that 14 = 2 • 13 everywhere in Sect.3.2.3.The surface reference value for DI 14 C for the virtual flux calculation is 1944 µmol m −3 , the same as for DI 13 C (and DI 12 C).
To compare the model-simulated δ 14 C values that we save as diagnostics (see Eq. 14) with published observations of 14 C, we apply the same fractionation correction to it that is used for observations to convert δ 14 C to 14 C: In the following we always show 14 C.As for the abiotic 14 C implementation, the value of 14 C atm can be set to be a constant or it can be read in from three files (one for the Northern Hemisphere, one for the equatorial region, and one for the Southern Hemisphere).

Ecosystem driver
We added an ecosystem driver (ecosys_driver) to the ocean model of the CESM in order to make it easier to expand the model to carry additional passive tracers that require variables from the ecosystem model, without adding these additional tracers to the ecosystem model itself.The ecosystem driver is structured similarly to the passive_tracers subroutine that calls all passive tracer modules, but it handles only the passive tracers that use the ecosystem model (see Fig. 1).It is called from the passive_tracers subroutine, and determines how many ecosystem-related passive tracers the model carries based on the namelist options set at buildtime.It then calls all subroutines in the ecosystem model and the related tracer modules, after being called by passive tracers with the corresponding tracer indices.Variables computed in the ecosystem model but used by other modules are shared via the new ecosys_share module.Only the ecosystem model changes the value of the variables in ecosys_share at this point.Other modules currently only read them from there, but do not modify them.With this infrastructure in place, additional tracers can be easily added without changing the ecosystem model too much.The only changes to the ecosystem model should be the copying of ecosystem variables to ecosys_share if they need to be shared with a new module as well as potentially the addition of new definitions and calculations of derived ecosystem variables that are needed but that are not currently computed in the ecosystem model (or not present in the required format, i.e., defined as local 2-D variables instead of a global 3-D variable).Nitrogen isotopes in the ocean model have already been added using this infrastructure (S.Yang, personal communication, 2014).

Simulations and spin-up
We have performed several simulations with the new carbonisotope-enabled model.As described in Sect.2, we used the ocean-only version of the CESM1.0.5, at a nominal 3 • horizontal resolution, forced by CORE-II climatological forcing (Large and Yeager, 2009).To spin up the carbon isotopes, we performed spin-up simulations that lasted several thousands of years.Radiocarbon takes a long time (5000-15 000 years, according to Orr et al., 2000) to equilibrate, due to the long timescale of deep ocean ventilation.
The abiotic radiocarbon has been spun up for 10 000 years using an atmospheric CO 2 concentration of 284.7 ppm and a 14 C value of 0 ‰.The abiotic DI 14 C and DI 12 C were started from the standard ecosystem initial conditions, scaled to yield a global initial state of 0 ‰ 14 C (following Orr et al., 2000), in order to simplify early interpretation and code verification.After 10 000 simulated years, the models satisfy the OCMIP2 surface CO 2 flux criterion of less than 0.01 Pg C year −1 .In terms of the drift in 14 C, 91 % of the ocean volume is spun up to the OCMIP2 criterion of a drift of less than 0.001 ‰ year −1 (compared to the required 98 % for OCMIP2).Compared to the fully spun-up solution (obtained using a new online Newton-Krylov method, manuscript in preparation by K. Lindsay, NCAR), differences are seen in the deep ocean only.
For the biotic carbon isotopes, we spun up the carbon isotopes for 6010 years, starting from the initial conditions of the ecosystem model, scaled to give a δ 13 C of 0 ‰ and a 14 C of −100 ‰.The atmospheric CO 2 concentration was set to 284.7 ppm, the atmospheric 14 C was set to 0 ‰, and the atmospheric δ 13 C was set to −6.379 ‰.In order to study the different biological fractionation parameterizations, two additional spin-up simulations were branched from the first spin-up simulation at year 2560 and run to year 6010.After 6010 years, the surface CO 2 flux is well below the OCMIP2 criterion of less than 0.01 Pg C year −1 , and over 99.99 % of the ocean volume shows a drift of less than 0.001 ‰ year −1 in δ 13 C.However, only 26 % of the ocean satisfies the OCMIP2 criterion of a drift of less than 0.001 ‰ year −1 for 14 C. Another 4000 years or more are likely required to get the biotic 14 C fully spun-up according to the OCMIP2 criterion.However, if we weaken the OCMIP2 criterion by an order of magnitude to less than 0.01 ‰ year −1 , 99.98 % of the ocean satisfies this new criterion for 14 C. Due to the long time required to run the ocean model with the ecosystem and the biotic carbon isotopes (the 6010 years took over 7 months of constant running on a supercomputer), we are currently not able to run the biotic radiocarbon to full equilibrium.In order to reach equilibrium in the future, a fast spin-up technique for the ecosystem model is currently in development by Keith Lindsay and will be applied to the biotic carbon isotopes when it is ready.We believe that for the purpose of this paper, which documents the implementation of the carbon isotopes in the model, the current spin-up is sufficient, in particular because future science applications will use the fully coupled 1 • CESM, rather than the 3 • ocean-only model version used here for model development.For these future science applications, the biotic radiocarbon should be spun up further in order to be fully trustworthy, and the carbon isotope simulation should be carefully validated in this model configuration.
After the long spin-up, we then performed experiments from 1765 to 2008, with the initial conditions from the end of the spin-up simulations in year 6010 for the biotic carbon isotopes and in year 10 000 for the abiotic radiocarbon.The atmospheric CO 2 , 14 C, and δ 13 C were prescribed based on the OCMIP-2 files (Orr et al., 2000) up to 1989, and H. Graven's formulation of the global average for 1990-2008 (personal communication, 2012).The atmospheric state was the same repeating climatological CORE-II forcing as used for the spin-up, so changes related to warming or changes in the wind forcing over the 20th century are not included.We continued with the climatological CORE-II forcing rather than use the interannually varying CORE-II forcing for 1948-2007 in order to avoid shocks to the ocean when switching the forcing and when the forcing jumps from 2007 back to 1948 every 60 years.This jump in the forcing impacts the simulation for 10 years or more (as described in Danabasoglu et. al., 2014), and would overlap with the start of the introduction of bomb radiocarbon into the atmosphere.
We also continued the spin-up simulations for 243 years, so that we could remove the influence of any continuing drift on the radiocarbon results shown in Sect.4.2.2.To investigate the influence of the net CO 2 uptake on the simulation results in the second part of the 20th century, we also performed sensitivity experiments where the atmospheric CO 2 was fixed at 1949 conditions, while 14 C atm and δ 13 C atm changed as usual.

Simulated distributions of 14 C
The radiocarbon simulation shows good agreement with the gridded GLODAP data for the 1990s (Key et al., 2004), reflecting the main features of the 14 C distribution: (i) at the surface (see Fig. 2) the model shows the observed M shape of 14 C distribution, with the highest values in the relatively stable subtropical waters, intermediate values in the equatorial upwelling zone, and low values in the polar regions, where the residence time is short and sea ice limits the uptake of atmospheric 14 C, with the overall lowest values in the Southern Ocean, where the upwelling of old, low 14 C waters further dilutes the surface waters.(ii) In the zonal mean (see Fig. 3), newly formed deepwater with high 14 C values can clearly be separated from old water masses with low 14 C values.(iii) In the depth profiles (see Fig. 4), it is obvious that the 14 C in the deep water decreases from the Atlantic Ocean over the Indian Ocean to the Pacific Ocean, which has the lowest 14 C values (i.e., oldest water).Consistently, the abiotic 14 C values are higher than the biotic 14 C values, but both show the same general features also shown in GLODAP (Key et al., 2004) and in the cruise data compiled by Schmittner et al. (2013) because their distribution is set mainly by the physical ocean simulation.The difference between the abiotic and biotic simulations due to biological effects has been estimated to be on the order of 10 % (Fiadiero, 1982), but since the biotic radiocarbon simulation is less spun-up than the abiotic simulation at this point, a detailed investigating of the impact of the biological effects is premature and will be the topic of a future study when we can spin up both radiocarbon implementations using a fast-spin up technique.
Above 1000 m, the depth structure of the simulated 14 C agrees reasonably well with observations, with the best agreement with the GLODAP 14 C in the upper 250 m of the Indian Ocean (see Fig. 4).The largest biases are found at depths below 1000 m (see Fig. 4), with the model showing 14 C values that are too negative (i.e., water that is too old).The largest bias is located in the deep Pacific, where the 14 C is up to 100 ‰ too negative (see Figs. 3 and 4).In terms of radiocarbon age, the maximum bias in the deep Pacific is 1000 years compared to GLODAP, revealing that the deep Pacific Ocean in the model is not ventilated as much as it should be.This is a well-known bias in the CESM, which was also present in the ocean model of a previous version of the CESM, the CCSM3 (Graven et al., 2012), as well as in the nominal 1 • resolution version of the current CESM1 ocean model (Bardin et al., 2014), and is related to too weak Antarctic Bottom Water formation in the CESM (Danabasoglu et. al., 2011) and too shallow mixed layers in the Southern Ocean (Moore et al., 2013).Currently radiocarbon is being used to test improvements to the ventilation in the Southern Ocean in the ocean model in the CESM, in order to improve this bias in future versions of the CESM (K.Lindsay, personal communication, 2014).

14 C bomb inventory
The excess oceanic radiocarbon inventory is frequently used to investigate the ocean uptake of anthropogenic carbon (e.g., Key et al., 2004;Graven et al., 2012) and to determine the mean gas exchange velocity used in ocean models (e.g., Wanninkhof, 1992;Sweeney et al., 2007;Naegler et al., 2006;Naegler, 2009).To establish how well the newly developed radiocarbon tracer compares to observations, we here compare the simulated excess radiocarbon inventory with observational estimates.The excess radiocarbon in the ocean includes change in the oceanic radiocarbon from the atmospheric nuclear bomb tests of the 1950s and 1960s, as well as from the Suess effect and changes in net CO 2 uptake, compared to the reference period of the 1940s, following Naegler (2009).In 1975, the excess radiocarbon inventory in the abiotic and biotic simulation is 297 × 10 26 atoms 14 C and 295 × 10 26 atoms 14 C, respectively.This lies within the range of observational estimates of the excess radiocarbon in 1975, which range from 225 ×10 26 atoms 14 C to 314 ± 35 × 10 26 atoms 14 C (see Table 2).It has been shown that the earlier estimates from Broecker et al. (1985Broecker et al. ( , 1995) ) were high by about 25 % (e.g., Hesshaimer et al., 1994;Peacock, 2004;Sweeney et al., 2007), which suggests that the simulated values are probably on the high end of the observational range.One reason for this could be the choice of the coefficient a = 0.31 cm h −1 in Eq. ( 3), which has been shown  2013) (top row), the gridded GLODAP data (Key et al., 2004) (second row), the 14 C from the biotic model (third row), and the abiotic model (bottom row).Note that due to the sparse observational data (see Fig. 2a for the coverage at the surface), the zonal average from the cruise data in the top row is more of a zonal composite than a zonal average.to be high (e.g., Sweeney et al., 2007;Naegler, 2009).Graven et al. (2012) showed that in the ocean model of the CCSM3, the simulated excess radiocarbon inventory was lower when a coefficient a = 0.23 cm h −1 rather than a = 0.31 cm h −1 was used in Eq. (3).However, since a = 0.31 cm h −1 is the parameter used in the CESM in general to compute air-sea gas fluxes, we did not change it here.For 1995, the excess ra-diocarbon inventories in the abiotic and biotic simulation are 389×10 26 atoms 14 C and 390×10 26 atoms 14 C, respectively, which is close to but slightly higher than the observational estimates of 313-383 ×10 26 atoms 14 C (see Table 2).
The natural radiocarbon inventory, before anthropogenic disturbances from the Suess effect and from increased oceanic net CO 2 uptake, has been estimated to be 19 000 ± Table 2. Excess oceanic radiocarbon inventory, measured in 10 26 atoms of 14 C, from various sources for 1975 (GEOSECS) and 1995 (WOCE).Corrections by Naegler et al. (2006) are for neglected ocean regions; corrections by Naegler (2009) are for neglected contributions from increasing DIC.The values from this study are listed at the bottom, for the abiotic and biotic implementation.The biotic excess radiocarbon inventories are the same for all biological fractionation choices tested.1200 × 10 26 atoms of 14 C (Naegler, 2009).In the simulation the inventory is just outside the error bar for the biotic model (17 763-17 770 × 10 26 atoms of 14 C, depending on the biological fractionation used), and slightly lower for the abiotic model (16 190×10 26 atoms of 14 C).The natural radiocarbon inventories are calculated for years 6185-6194 of the control simulations, which corresponds to the same total run time as years 1940-1949 in the 1765-2008 experiments, which were started from the control in year 6010.However, the biotic model estimate of the natural radiocarbon inventory might still not be the final value, as the biotic radiocarbon is still spinning up.In terms of the anthropogenic radiocarbon inventories presented next, this biases should not play any large role, however, as we remove any remaining drift.Specifically, to calculate the early anthropogenic radiocarbon inventory present in the 1940s, we take the difference between the natural radiocarbon inventory in simulation years 6185-6194 (with constant atmospheric CO 2 , 14 C, and δ 14 C) and the inventory in the 1940s (with changing atmospheric CO 2 , 14 C, and δ 14 C since 1765).By taking this difference between years of equal total run time, we remove the impact of any remaining drift in 14 C. We find an anthropogenic radiocarbon inventory of 20 × 10 26 atoms of 14 C for the abiotic model and 5 × 10 26 atoms of 14 C for the biotic model (independent of the biological fractionation used).Both of these anthropogenic radiocarbon inventories for the 1940s are within the error bar of the estimate of 4 ± 20 × 10 26 of 14 C from Naegler (2009), with the biotic model giving a very good match.
Using sensitivity experiments from 1950 to 2008 with atmospheric CO 2 held constant at 1949 levels but normally increasing atmospheric 14 C, we can calculate the impact of increased ocean uptake of anthropogenic CO 2 on the excess radiocarbon inventory: in 1975, the excess oceanic radiocarbon inventory relative to the 1940s due to atmospheric 14 C changes alone (from the atmospheric bomb tests and the Suess effect) is 282 × 10 26 atoms of 14 C for the abiotic model and 280 × 10 26 atoms of 14 C for the biotic model, while for 1995 the numbers are 353 × 10 26 atoms of 14 C and 354 × 10 26 atoms of 14 C, respectively.This means that the increase in net CO 2 uptake contributed 15 × 10 26 atoms of 14 C in 1975 and 36 × 10 26 atoms of 14 C in 1995 compared to the 1940s (for both the abiotic and biotic models), which is 5 and 9 % of the total radiocarbon excess in these years, respectively.These changes are in excellent agreement with calculations from Naegler (2009), which showed an excess radiocarbon inventory in 1995 of 346 ± 98 × 10 26 atoms 14 C due to atmospheric 14 C changes, and 27 ± 9 × 10 26 atoms 14 C due to net CO 2 uptake.The percentage contribution of the net CO 2 uptake to the total radiocarbon excess was given as 3 % in 19753 % in and 8 % in 19953 % in in Naegler (2009)), which again compares very well with our model simulation.

Simulated δ 13 C and the impact of different biological fractionation parameterizations
In the literature, models of biological fractionation are still under debate (e.g., Keller and Morel, 1999).We therefore tested three different parameterizations of biological fractionation, to investigate the impact on the simulated δ 13 C (as described in Sects.3.2.3 and 4.1).As shown in Fig. 5a, the simulated globally averaged p depth profiles differ when these different parameterizations are used, with p values ranging from 15 to 30.By design, p is the same for diatoms, diazotrophs, and small phytoplankton when using Rau et al. (1989), while p shows large variations between species for the method of Keller and Morel (1999), due to the dependence on species-specific cell parameters (see Table 1).The method of Laws et al. (1995) leads to small differences between species in the surface ocean only.Below 200 m, only the p following Rau et al. (1989) still changes with depth (see Fig. 5a), due to the sole dependence of p on CO * 2 and the export of organic carbon and carbonates to depth.
The impact of the different biological fractionation choices on δ 13 C DIC is noticeable (see Fig. 5b), with the globally averaged δ 13 C DIC based on p from Rau et al. (1989) being larger below 150 m compared to the δ 13 C DIC from Laws et al. (1995) and Keller and Morel (1999), but slightly smaller at the surface.Despite the more complex formulation of p in Keller and Morel (1999) compared to Laws et al. (1995) and the significantly different p profiles, the resulting δ 13 C DIC from both methods is very similar and only differs slightly at depth (most notably between 150 and 2000 m).To compare the simulated δ 13 C DIC to the cruise data of δ 13 C DIC compiled by Schmittner et al. (2013), we re-gridded the model output to subsample the model at the same points as covered by the cruise data.The resulting globally averaged depth profiles are shown in Fig. 5c, and are remarkably similar to the full globally averaged model results in Fig. 5b.Both show the expected increase in δ 13 C DIC directly below the surface, due to the preferential uptake of the light isotope during photosynthesis, followed by the expected decrease in δ 13 C DIC with depth due to the remineralization of the isotopically light organic material back into the water column.The modelsimulated global depth profile of δ 13 C DIC lies within the error range of ±0.2 ‰ around the cruise δ 13 C DIC data between the surface and 150 m and below 1000 m, but shows smaller δ 13 C DIC values than observed between 150 and 1000 m.
For individual basins, the model bias compared to the cruise data is smallest in the Atlantic, with the δ 13 C DIC based on the biological fractionation from Rau et al. (1989) almost entirely within the uncertainty range of the data (see Fig. 5d).All three basins contribute to the bias seen between 150 and 2000 m in the global average, with the Indian Ocean contributing the most to this bias in the upper ocean and the Pacific Ocean contributing the most at intermediate depths (see Fig. 5c-f).In general, the model simulated δ 13 C DIC tends to be smaller than the observed δ 13 C DIC .While the difference between the full global average in the model and the subset global average based on the cruise data locations is small, the difference between the total basin average (shown as dashed lines in Fig. 5d-f) and the subset basin averages (shown as solid lines) is larger for the individual basins.
At the surface, the simulated δ 13 C DIC values show a systematic bias in that they are generally larger than the observational data suggest, but the same general spatial pattern is visible (see Fig. 6).While both gas-exchange and biological processes are important for the surface ocean δ 13 C DIC pattern (Schmittner et al., 2013), the details of the biological fractionation parameterizations appear to have a very small impact at the surface, as shown in the almost identical surface distributions from the model (see Fig. 6c-e).The zonal means of δ 13 C DIC from the different biological fractionation parameterizations on the other hand do show some small differences (see Fig. 7), with the biological fractionation from Rau et al. (1989) leading to the largest δ 13 C DIC values in all three ocean basins, while the fractionation based on Keller and Morel (1999) shows the lowest δ 13 C DIC values.Overall, all three parameterizations lead to the expected pattern of high values of δ 13 C DIC in water that has recently been in contact with the surface (e.g., North Atlantic Deep Water) and low δ 13 C DIC values in water that has been out of contact with the atmosphere for a long period of time and has accumulated a large amount of remineralized (isotopically light) organic mater (e.g., in the deep Pacific).
We choose the biological formulation from Laws et al. (1995) as the default biological fractionation in our model, as it considers the growth rate of different species, but the differences in the simulated δ 13 C DIC compared to the more complex formulation from Keller and Morel (1999) are small.The other parameterizations of biological fractionation remain an option in the model that can be chosen at build time.

Oceanic surface 13 C Suess effect
The surface oceanic Suess effect, which is the decrease in the surface ocean δ 13 C due to the penetration of carbon originating from the burning of fossil fuels, has been calculated from observational data as well as from other models that include 13 C, and it is often used to derive the oceanic anthropogenic carbon uptake (e.g., McNeil et al., 2001;Tagliabue and Bopp, 2008).In our model simulation, the surface δ 13 C change between 1975 and 1995 is −0.164 to −0.167 ‰ decade −1 (the range is for the different biological fractionations used).This compares well with other estimates of −0.171 ‰ decade −1 (Bacastow et al., 1996), −0.18 ‰ decade −1 (Gruber et al., 1999), −0.15 ‰ decade −1 (Sonnerup et al., 1999), and −0.174 ‰ decade −1 (Tagliabue and Bopp, 2008).As already shown by Quay et al. (1992) and Gruber et al. (1999), the surface ocean Suess effect is not uniform (see Fig. 8), and the model simulation of the spatial   1989), (d) Laws et al. (1995), and (e) Keller and Morel (1999).
Suess effect agrees well with the model results of Tagliabue and Bopp (2008): the largest changes (i.e., most negative values in Fig. 8) occur in regions with little deep ventilation and therefore longer residence times of water at the surface (e.g., the subtropical gyres), while the smallest changes (i.e., least negative or zero in Fig. 8) occur in regions of reduced airsea gas exchange (e.g., under sea ice), in regions with active deep convection (and therefore short residence times at the surface, e.g., around the Antarctic), as well as in regions with upwelling (which dilutes the surface δ 13 C, for example, off the western coast of South America).
Compared to the pre-industrial ocean, the total surface ocean 13 C Suess effect is −0.064 to −0.065 ‰ decade −1 for 1860-2000 (depending on the different biological fractionations used), compared to −0.07 ‰ decade −1 found by Tagliabue and Bopp (2008).The fact that the simulated oceanic  13 C Suess effect calculated over different periods agrees reasonably well with other available estimates suggests that our model is able to simulate the change in the oceanic δ 13 C inventory correctly, despite some mean biases in the distribution of δ 13 C described and shown in Sect.4.3.1.

Summary
We have developed carbon isotope tracers in the ocean model of the CESM, including a biotic and an abiotic radiocarbon tracer and a biotic 13 C tracer.The details of the implementation are described here in order to serve as a reference for future users of these new model features and/or for model developers planning to modify the code or add carbon isotopes to other ocean models.In particular, we tested three different formulations for the fractionation during phytoplankton growth that have been discussed in the literature, and show that the effect on the simulated δ 13 C in the ocean is relatively minor.A comparison of the simulation results from the coarse nominal 3 • resolution ocean model forced with climatological CORE-II atmospheric forcing and with present-day data for 14 C and δ 13 C shows that the simulated carbon isotopes can represent the large-scale features of the observed distributions as well as the anthropogenic changes due to nuclear bomb tests and the burning of fossil fuels.The carbon isotopes also reflect some known model biases, for example too sluggish ventilation of the deep Pacific Ocean.Once a fast spin-up technique for the biotic carbon isotopes has been implemented, we are planning to further validate the carbon-isotope simulation in the fully coupled CESM framework at 1 • resolution.Ultimately, we plan to use the carbon isotopes for both present-day and paleo simulations in the fully coupled framework of the CESM at the standard nominal 1 • resolution in the ocean, in order to investigate details of changes in the ocean circulation over the 20th century, the last Millennium, and at the Last Glacial Maximum.

Figure 1 .
Figure 1.Schematic of the passive tracer modules with the new ecosystem driver and carbon isotope modules.Existing modules are shown in blue, new modules are shown in red, and edited modules are shown in blue with a red box.Dashed lines indicate future developments.This schematic shows how the ecosystem driver acts as an interface between the ecosystem-related modules and the passive tracers module that drives all tracer modules as well as how ecosys_share is used to share variables computed by the ecosystem model and used by other modules beside the ecosystem model.

Figure 3 .
Figure 3. Zonal averages of total 14 C for the Atlantic, Pacific, and Indian oceans for the 1990s, from cruise data compiled by Schmittner et al. (2013) (top row), the gridded GLODAP data(Key et al., 2004) (second row), the 14 C from the biotic model (third row), and the abiotic model (bottom row).Note that due to the sparse observational data (see Fig.2afor the coverage at the surface), the zonal average from the cruise data in the top row is more of a zonal composite than a zonal average.

Figure 4 .
Figure 4. Depth profiles of 14 C for (a) the global ocean, (b) the Atlantic Ocean, (c) the Pacific Ocean, and (d) the Indian Ocean.The simulated biotic (green) and abiotic (blue) 14 C is compared to the global gridded GLODAP 14 C (black) data set (Key et al., 2004).In addition, dashed lines show the cruise data compiled by Schmittner et al. (2013) (gray) and the model-simulated data subsampled at the same locations as these data (green and blue dashed lines).

Figure 5 .
Figure 5. (a) Depth profiles over the top 500 m (where p is important because of primary production) of the globally averaged values of p produced by the three tested parameterizations for biological fractionation for diatoms (solid line), diazotrophs (short dashes), and small phytoplancton (large dashes).The simulated globally averaged depth profile (0-6000 m) of δ 13 C DIC in the 1990s is shown in (b), and the global average depth profile of the subset model δ 13 C DIC for the same grid points as in the cruise data compiled by Schmittner et al. (2013) is shown in (c).Basin average depth profiles are shown in (d-f), with dashed lines showing the full basin average from the model and solid lines showing the subset averages for the same points as the cruise data compiled by Schmittner et al. (2013).The uncertainty for the cruise data is shown as gray shading in (c), and is ±0.2 ‰ (Schmittner et al., 2013).Note that the irregular y axis in (b-f) emphasizes the upper ocean.

Figure 7 .
Figure 7. Zonal ocean basin composites from the cruise data compiled by Schmittner et al. (2013) (top row) compared to 1990s zonal basin averages from the model simulation using the biological fractionation from Rau et al. (1989) (second row), Laws et al. (1995) (third row), and Keller and Morel (1999) (bottom row).