Simulating stable carbon isotopes in the ocean component of the FAMOUS general circulation model with MOSES1 (XOAVI)

. Ocean circulation and the marine carbon cycle can be indirectly inferred from stable and radiogenic carbon isotope ratios ( δ 13 C and (cid:49) 14 C, respectively), measured directly in the water column, or recorded in geological archives such as sedimentary microfossils and corals. However, interpreting these records is non-trivial because they reﬂect a complex interplay between physical and biogeochemical processes. By directly simulating multiple isotopic tracer ﬁelds within numerical models, we can improve our understanding of the processes that control large-scale isotope distributions and interpolate the spatiotemporal gaps in both modern and palaeo datasets. We have added the stable isotope 13 C to the ocean component of the FAMOUS coupled atmosphere– ocean general circulation model, which is a valuable tool for simulating complex feedbacks between different Earth system processes on decadal to multi-millennial timescales. We tested three different biological fractionation parameterisations to account for the uncertainty associated with equilibrium fractionation during photosynthesis and used sensitivity experiments to quantify the effects of fractionation during air–sea gas exchange and primary productivity on the simulated δ 13 C DIC distributions. Following a 10 000-year preindustrial spin-up, we simulated the Suess effect (the isotopic imprint of anthropogenic fossil fuel burning) to assess the performance of the model in replicating modern observations. Our implementation captures the large-scale structure and range of δ 13 C DIC observations in the surface ocean, but the simulated values are too high at all depths, which we infer is due to biases in the biological pump. In the ﬁrst in-stance, the new 13 C tracer will therefore be useful for recalibrating both the physical and biogeochemical components of FAMOUS.

Abstract. Ocean circulation and the marine carbon cycle can be indirectly inferred from stable and radiogenic carbon isotope ratios (δ 13 C and 14 C, respectively), measured directly in the water column, or recorded in geological archives such as sedimentary microfossils and corals. However, interpreting these records is non-trivial because they reflect a complex interplay between physical and biogeochemical processes. By directly simulating multiple isotopic tracer fields within numerical models, we can improve our understanding of the processes that control large-scale isotope distributions and interpolate the spatiotemporal gaps in both modern and palaeo datasets. We have added the stable isotope 13 C to the ocean component of the FAMOUS coupled atmosphereocean general circulation model, which is a valuable tool for simulating complex feedbacks between different Earth system processes on decadal to multi-millennial timescales. We tested three different biological fractionation parameterisations to account for the uncertainty associated with equilibrium fractionation during photosynthesis and used sensitivity experiments to quantify the effects of fractionation during air-sea gas exchange and primary productivity on the simulated δ 13 C DIC distributions. Following a 10 000-year preindustrial spin-up, we simulated the Suess effect (the isotopic imprint of anthropogenic fossil fuel burning) to assess the performance of the model in replicating modern observations. Our implementation captures the large-scale structure and range of δ 13 C DIC observations in the surface ocean, but the simulated values are too high at all depths, which we infer is due to biases in the biological pump. In the first instance, the new 13 C tracer will therefore be useful for recalibrating both the physical and biogeochemical components of FAMOUS.

Introduction
Carbon isotopes are often used as proxies for ocean circulation and the marine carbon cycle. There are three naturally occurring carbon isotopes: the stable isotopes 12 C (98.9 %) and 13 C (1.1 %), and the radioactive isotope 14 C (1.2 × 10 −10 %), which is also known as radiocarbon (Key, 2001). In this study, we focus on the stable isotopes, with 14 C being discussed in detail elsewhere (Dentith et al., 2019a). The relative proportions of 12 C and 13 C in a given oceanic pool (e.g. dissolved inorganic carbon, DIC, or particulate organic carbon, POC) are controlled by ocean circulation and mixing, and mass-dependent fractionation during biogeochemical processes such as air-sea gas exchange (Lynch-Stieglitz et al., 1995;Zhang et al., 1995), photosynthesis (e.g. Sackett et al., 1965;Rau et al., 1989;Hollander and McKenzie, 1991;Keller and Morel, 1999), and calcium carbonate formation (Emrich et al., 1970;Turner, 1982;Ziveri et al., 2003). This is typically reported in delta (δ) notation, which is the heavy to light isotope ratio (R) of a sample relative to a standard in per mille (‰) units: δ 13 C = 13 C/ 12 C sample 13 C/ 12 C standard − 1 × 1000. (1) Oceanic δ 13 C is primarily used to track individual water masses (Curry and Oppo, 2005), study past changes in the carbon cycle (e.g. de la Fuente et al., 2017), and investigate changes in ocean circulation on glacial-interglacial timescales (e.g. Spero and Lea, 2002;Campos et al., 2017). It has also been used to constrain air-sea gas exchange rates (Gruber and Keeling, 2001) and to estimate the uptake of an-thropogenic carbon by the global oceans (Quay et al., 1992(Quay et al., , 2003. Oceanographic surveys conducted since the 1970s, such as the World Ocean Circulation Experiment (WOCE; Orsi and Whitworth, 2005;Talley, 2007Talley, , 2013Koltermann et al., 2011), and synthesis projects such as Carbon dioxide in the Atlantic Ocean (CARINA; Key et al., 2010), Pacific Ocean Interior Carbon (PACIFICA; Suzuki et al., 2013), and the Global Ocean Data Analysis Project (GLODAP; Key et al., 2004Key et al., , 2015Olsen et al., 2016), provide an indication of large-scale carbon isotope distributions in the modern oceans. The two main drawbacks of these surveys are that they include relatively few measurements from the subsurface ocean and that there were only a limited number of repeat measurements at fixed locations, which were often taken decades apart. These datasets are therefore insufficient for studying transient changes in carbon isotope distributions at subdecadal resolution.
Geological archives such as corals (e.g. Guilderson et al., 2013) and sediment cores (e.g. Oliver et al., 2010) are used to extend the record further back in time. However, interpreting isotopic ratios in geological archives is non-trivial because they result from a complex interplay between physical processes and biogeochemical processes, both in the water column itself and during biomineralisation, which can be difficult to disentangle.
By including carbon isotopes in climate models, we can fill in the spatiotemporal gaps in both modern and palaeo datasets, and improve our understanding of the processes that control their large-scale distributions (Tagliabue and Bopp, 2008;Schmittner et al., 2013;Menviel et al., 2017). The Ocean Carbon-Cycle Modelling Intercomparison Project (OCMIP) was initiated in 1995 with the aim of evaluating the major differences between global ocean carbon cycle models and advancing our understanding of the ocean as a long-term CO 2 reservoir (Orr, 1999). Carbon isotopes are not routinely incorporated into climate models because of the computational expense associated with the long equilibration between the deep ocean and the atmosphere (Bardin et al., 2014). However, since OCMIP produced a legacy of standard input fields (Orr, 1999;Orr et al., 2000Orr et al., , 2017, carbon isotopes have increasingly been implemented into models of varying complexities to validate physical and biogeochemical schemes, to investigate the spatiotemporal variability in isotope distributions, and to reconcile the interpretation of ocean proxy data. As outlined in Table 1, the community of 13 C-enabled models currently includes HAMOCC3.1 (Hofmann et al., 2000), the Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model (MOM; Murnane and Sarmiento, 2000), CLIMBER-2 (Brovkin et al., 2002), MoBidiC (Crucifix, 2005), GENIE (Ridgwell et al., 2007), Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) (Tagliabue and Bopp, 2008), LOVECLIM (Mouchet, 2011;Menviel et al., 2015), Bern3D+C (Tschumi et al., 2011), the University of Victoria (UVic) Earth system model (ESM) (Schmittner et al., 2013), iLOVECLIM (Bouttes et al., 2015), the Community Earth System Model (CESM) (Jahn et al., 2015), and the Commonwealth Scientific and Industrial Research Organisation Mark 3L climate system model with the Carbon of the Ocean, Atmosphere and Land (CSIRO Mk3L-COAL) biogeochemical model (Buchanan et al., 2019). Most of these are lowresolution (3 to 5 • ), intermediate-complexity models that are valuable tools for studying changes in ocean biogeochemistry on multi-millennial timescales. The more complex models (e.g. PISCES and CESM) provide a more sophisticated representation of physical and biogeochemical processes because of increased spatial resolution and/or the inclusion of more carbon pools. However, the higher-complexity models are computationally expensive; for example, at the time of their study, a 6010-year spin-up simulation with CESM took over 7 months to run (Jahn et al., 2015). Without employing offline or accelerated spin-up techniques (e.g. Lindsay, 2017), these models are therefore less practical for running the long simulations required to fully spin up the components of the Earth system that evolve on millennial timescales, such as deep ocean circulation (England, 1995) and ocean biogeochemical cycles (Falkowski et al., 2000;Key et al., 2004).
Here, we describe the implementation of 13 C in the ocean component of the FAMOUS general circulation model (GCM). FAMOUS is well suited to studying complex interactions between different components of the Earth system on decadal to multi-millennial timescales, due to its reduced spatial resolution and increased time step relative to the latest generation of state-of-the-art GCMs (Sect. 2.1). We use sensitivity experiments to quantify the effects of isotopic fractionation during air-sea gas exchange and primary productivity on the simulated δ 13 C DIC distributions (Sects. 2.3.3 and 3.1) and test three different parameterisations for photosynthetic fractionation to account for the uncertainty associated with the relative influence of ambient conditions, physiological effects and transport mechanism on the fractionation of carbon isotopes during photosynthetic CO 2 fixation (Sects. 2.2.2 and 3.3). We evaluate the overall performance of the model in simulating large-scale δ 13 C DIC distributions by comparing to modern observations (Sect. 3.2) and discuss the potential of the new 13 C tracer as a tuning target for future recalibration work (Sect. 3.4).

Model description
FAMOUS is a coupled atmosphere-ocean GCM (Jones et al., 2005;Smith et al., 2008;Smith, 2012;Williams et al., 2013) based on HadCM3 (Gordon et al., 2000;Pope et al., 2000). Both are configurations of the UK Met Office Unified Model (UM) version 4.5 (Valdes et al., 2017). The quasi-hydrostatic primitive equation atmospheric model is Maier-Reimer (1993), Popp et al. (1989), Rau et al. (1996)  in latitude by 7.5 • in longitude, with 11 vertical levels on a hybrid sigma-pressure coordinate system. The rigidlid ocean model has a horizontal resolution of 2.5 • × 3.75 • and 20 unevenly spaced vertical levels, which are approximately 10 m thick in the near-surface ocean and 600 m thick in the deep ocean. The atmosphere and ocean operate on 1 and 12 h time steps, respectively, and are coupled once per model day. The model currently includes oxygen (Williams et al., 2014) and chlorofluorocarbons (Pope et al., 2000) as optional tracers. At the time of this study, FA-MOUS is capable of simulating 400 to 500 model years per wall-clock day on tier 2 (regional) and tier 3 (local) highperformance computers at the University of Leeds, which is more than 5 times the run speed of HadCM3. This makes FAMOUS ideal for running long (multi-millennial) simulations (Smith and Gregory, 2012;Gregoire et al., 2012Gregoire et al., , 2015 or large (hundred-member) ensembles (Gregoire et al., 2011;Sagoo et al., 2013). Further technical documentation can be found in an ongoing special issue in Geoscientific Model Development (http://www.geosci-model-dev.net/ special_issue15.html, last access: 23 December 2019).
We added 13 C as an optional passive tracer into the ocean component of FAMOUS, using the Met Office Surface Exchange Scheme (MOSES) version 1 (Cox et al., 1999) generation of the model to evaluate our scheme. Although a newer version of the land surface model exists, which includes the terrestrial carbon cycle and interactive vegetation (MOSES2.2; Essery et al., 2001Essery et al., , 2003Williams et al., 2013;Valdes et al., 2017), problems have been identified with its representation of meridional overturning circulation (MOC) in multi-millennial simulations with constant pre-industrial boundary conditions (Dentith et al., 2019b). Specifically, FAMOUS-MOSES2.2 simulates a collapsed Atlantic MOC (AMOC) and a strong, deep Pacific MOC when the run length exceeds 6000 years, resulting in spurious ocean tracer distributions. However, our code is directly transferable between the different generations of the model, meaning that the isotope system can be extended into the terrestrial carbon cycle following additional tuning to improve the physical ocean circulation in FAMOUS-MOSES2.2.

Hadley Centre Ocean Carbon Cycle Model (HadOCC)
The marine carbon cycle in FAMOUS is modelled by HadOCC, a coupled physical-biogeochemical model that simulates air-sea gas exchange, the circulation of DIC, and the cycling of carbon by marine biota (Palmer, 1998;Palmer and Totterdell, 2001). The ecosystem model provides a simplified representation of the ocean biological system, with a single (nitrogenous) nutrient, a single class of phytoplankton, a single class of (non-migrating) zooplankton, and detritus. Changes in the size of these pools are calculated through a series of coupled differential equations that describe primary production, respiration, mortality, grazing, excretion, and the sinking and remineralisation of detritus. The system is nitrogen limited and carbon flows are coupled to the nitrogen flows by stoichiometric ratios that are fixed for each pool of organic matter. In addition to the four biological components, HadOCC also explicitly simulates DIC and alkalinity. Modelled DIC concentrations depend upon phytoplankton growth and biological breakdown. Alkalinity is similarly affected by biological processes and is used to calculate the proportion of DIC that is in the form of CO 2 in the surface waters and consequently the air-sea CO 2 flux. All six tracers are advected, diffused, and mixed across all levels, although phytoplankton and zooplankton concentrations are negligible below the euphotic zone (approximately the uppermost 100 m of the ocean). Detritus is the only biological tracer that is subject to sinking, which is parameterised at a constant rate of 10 m d −1 . However, there is no representation of sediments: any detrital material that reaches the ocean floor is therefore immediately refluxed back to the top layer of the ocean to conserve carbon and nitrogen. Calcium carbonate (CaCO 3 ) production is represented as an instantaneous redistribution of DIC and alkalinity below the lysocline, the depth of which is spatially and temporally constant (approximately 2500 m below sea level). HadOCC accurately simulates low primary production in the subtropical gyres and high production in the regions with the greatest nutrient supply: the subpolar North Pacific and North Atlantic oceans, and around the Antarctic Convergence Zone (Fig. 1). However, primary production is higher than observed in the equatorial Pacific, which is attributed to excessive upwelling in the eastern equatorial Pacific (Palmer and Totterdell, 2001). Production is lower than observed northwards of 50 • N in the Atlantic and Pacific basins because sea ice formation and melt do not affect salinity distributions (an area for future development), although the model does include an iceberg meltwater flux (Smith et al., 2008). Consequently, stably stratified, low-salinity layers of meltwater, which promote phytoplankton growth, are not represented in the model (Palmer and Totterdell, 2001). Furthermore, the simulated production in coastal regions is lower than observed. There are three main reasons for this: (1) HadOCC does not simulate riverine input of nutrients, which are a significant source of coastal nutrients; (2) most of the coastlines in FAMOUS are directly adjacent to ocean grid cells that are more than 1 km deep, which dilutes nearsurface nutrient concentrations; and (3) upwelling is spread out over several grid points, which causes production to be more diffuse than observed (Palmer and Totterdell, 2001).
The level of representation of ecosystem processes in HadOCC is of intermediate complexity, making it computationally faster than more sophisticated ecosystem models that include additional POC species and/or multiple nutrients (e.g. PISCES). Previous studies have found that errors in biogeochemical simulations are largely driven by biases in the physical ocean circulation (i.e. inaccuracies in the climate or ocean model to which the ecosystem model has been cou-  (Behrenfeld and Falkowski, 1997), (b) the std simulation in the 1990s, and (c) simulated minus observed. Monthly mean primary productivity data were obtained from the Oregon State University Ocean Productivity website (http://www.science.oregonstate.edu/ ocean.productivity, last access: 23 December 2019). pled; Doney, 1999;Doney et al., 2004;Najjar et al., 2007). Thus, simulating carbon isotopes in a more complex ecosystem model would not necessarily yield substantially better results.

Carbon isotope implementation
We added 13 C to the four carbon pools in HadOCC: DIC, phytoplankton, zooplankton, and detritus (Fig. 2). We assume that modelled DIC is 12 C and carry 13 C as a ratio (DI 13 C/DI 12 C); therefore, virtual fluxes are not required to account for the dilution or concentration effects of surface freshwater fluxes (Appendix A). We also use model units to minimise the error associated with carrying small numbers: where 13 C/ 12 C std = 1.12372 × 10 −2 (Craig, 1957). We account for isotopic fractionation during air-sea gas exchange (Sect. 2.2.1 and Appendix B) and photosynthesis (Sect. 2.2.2 and Appendix C). Observational estimates suggest that isotopic fractionation during CaCO 3 formation is between +3 ‰ and −2 ‰ (Ziveri et al., 2003), which is small compared to the other fractionation effects (Turner, 1982). Previous 13 C isotope implementation studies have therefore assumed either no isotopic fractionation during CaCO 3 production (Schmittner et al., 2013) or prescribed constant values, for example, +1 ‰ (Tagliabue and Bopp, 2008) or +2 ‰ (Jahn et al., 2015). We conducted sensitivity tests where fractionation during CaCO 3 formation was included at constant rates of −2 ‰, 0 ‰, and +2 ‰, respectively. After 10 000 years, there was 0.001 ‰ difference in both the mean surface ocean δ 13 C DIC values and the surface standard deviations between all three simulations, and 0.02 ‰ difference between the three global volume-weighted integrals. Since these differences are small, we proceeded with the equivalent of no fractionation during CaCO 3 production (α CaCO 3 = 1.0).

Air-sea gas exchange
The air-sea gas flux of DI 12 C (F ) is calculated as where C sat is the saturation concentration of atmospheric CO 2 (in mol m −3 ), C surf is the surface aqueous concentration of CO 2 (in mol m −3 ), and PV is the piston velocity (in cm h −1 ), which is calculated as where a is a tuneable coefficient, u is the wind speed (in m s −1 ), a ice is the fractional ice cover, and Sc is the Schmidt number for CO 2 , calculated as a function of sea surface temperature (SST, in • C): The air-sea gas flux of DI 13 C/DI 12 C F 13 Figure 2. Schematic overview of the 13 C implementation in FAMOUS. Blue boxes represent permanent carbon pools. Grey boxes represent temporary carbon pools (note that CaCO 3 is a temporary carbon pool because the export of CaCO 3 in FAMOUS is represented as an instantaneous redistribution of alkalinity and carbon at depth). The orange box represents the prescribed atmospheric carbon pool. The dashed lines represent fluxes of 13 C/ 12 C. However, note that the outgassed 13 C/ 12 C has no effect on δ 13 C atm because FAMOUS does not currently have a fully interactive carbon cycle. Solid lines represent fluxes of 13 C. Dot-dashed lines represent processes that occur below the lysocline (≈ 2500 m below sea level). The dotted line represents the reflux of detrital material from the seafloor to the surface layer. Red lines represent fractionation effects. The orange line represents isotopic fractionation during calcium carbonate formation (α CaCO 3 ), which is included in the code as a user-specified constant. Note that all simulations presented in this study were run without fractionation during calcium carbonate formation (i.e. α CaCO 3 = 1.0, which is equivalent to a fractionation effect of 0 ‰).
where 13 A/ 12 A and 13 C/ 12 C are the 13 C/ 12 C ratios of the atmosphere and DIC, respectively, α k is the constant kinetic fractionation factor (0.99919), α aq←g is the temperaturedependent fractionation during gas dissolution: and α DIC←g is the temperature-dependent fractionation between aqueous CO 2 and DIC: All three fractionation factors are based on the equations of Zhang et al. (1995). However, following Schmittner et al. (2013), we neglect the effect that the carbonate fraction (f CO 3 ) has on α DIC←g because this is much smaller (0.05 ‰) than the temperature effect (3 ‰). Currently, atmospheric CO 2 and δ 13 C concentrations can either be held con-stant or prescribed from a file that contains a single global weighted-average value per year.

Photosynthesis
Isotopic fractionation during photosynthesis (α POC←DIC , herein α p ) is calculated as where α POC←aq is the equilibrium fractionation factor between aqueous CO 2 and POC. Empirical relationships for the different biogeochemical fractionation effects (α aq←g , α DIC←g , and α POC←aq ) have been established from laboratory experiments, modern oceans and lakes, and the sedimentary record. However, there are still uncertainties associated with the parameterisation of α POC←aq . Early studies investigated a potential temperature dependence of the carbon isotope composition of marine phytoplankton. For example, Sackett et al. (1965) proposed that photosynthetic fractionation is higher at lower temperatures (0.23 ‰ per • C) after observing that phytoplankton in the Drake Passage had more negative δ 13 C values than those in the tropics. Wong and Sackett (1978) also recorded small temperature effects (−0.13 ‰ to +0.36 ‰ per • C) in 17 species of marine phytoplankton; however, the authors concluded that the 15 ‰ range observed in their samples was primarily related to different metabolic pathways within the organisms. Numerous studies have suggested that the fractionation of carbon isotopes during photosynthetic CO 2 fixation relates to aqueous CO 2 concentrations (CO * 2 ) in the ambient environment (Popp et al., 1989;Rau et al., 1989;Jasper and Hayes, 1990;Hollander and McKenzie, 1991;Freeman and Hayes, 1992). However, these studies assumed that CO 2 only enters the phytoplankton by passive diffusion and neglected physiological effects, such as phytoplankton growth rate, cell size and geometry, and cell membrane permeability. Taking into consideration that physiological factors may modify, weaken, or eliminate the relationship between CO * 2 and photosynthetic fractionation, Rau et al. (1996) proposed a model that accounted for the isotopic composition of the ambient aqueous CO 2 , isotopic fractionation associated with diffusive transport into the cell, and isotopic fractionation associated with enzymatic, intracellular fixation. Laws et al. (1995) identified a linear relationship between phytoplankton growth rate, CO * 2 , and isotopic fractionation during photosynthesis under the assumption that the growth rate is proportional to the net transport of CO 2 into the cell. A later study by Laws et al. (1997), which analysed the same species of marine diatom over a larger range of CO * 2 , revised this to a non-linear relationship. Burkhardt et al. (1999) and Keller and Morel (1999) additionally included active bicarbonate transport in their calculations, recognising that aqueous CO 2 is not the only substrate for photosynthetic fixation and that processes other than diffusion can contribute to inor-ganic carbon acquisition. This has been a relatively inactive research area in the last 20 years, but there remains no single accepted model for fractionation during photosynthesis.
Consequently, previous carbon isotope implementation studies have used a number of different parameterisations for biological fractionation (Table 1), with the choice of scheme largely reflecting the complexity of the simulated biogeochemical and ecosystem processes. It is difficult to compare the success of the different parameterisations used by individual modelling groups because inter-model differences in the simulated isotopic distributions predominantly relate to resolution, complexity, and biases in the physical ocean circulation and ocean biogeochemistry, as opposed to the choice of fractionation scheme. However, Hofmann et al. (2000) tested three different fractionation schemes within a single model. In their study, the oversimplified assumption of constant biological fractionation, taken from Maier-Reimer (1993), failed to reproduce the observed latitudinal gradients in δ 13 C POC . Calculating the fractionation as a function of CO * 2 , as per Popp et al. (1989), successfully replicated the interhemispheric asymmetry in δ 13 C POC , but a growthrate-dependent fractionation (e.g. Rau et al., 1996) was required to additionally capture the seasonal variations. Jahn et al. (2015) also demonstrated differences between three different fractionation schemes within a single model. In their study, the simple scheme of Rau et al. (1989) produced lower δ 13 C DIC values in the surface ocean and higher δ 13 C DIC values below 150 m compared to the more complex parameterisations of Laws et al. (1995) and Keller and Morel (1999). The differences between the intermediate-complexity formulation (Laws et al., 1995) and the most complex formulation (Keller and Morel, 1999) were small, and the Laws et al. (1995) equation was chosen as the default scheme.
To account for the uncertainty associated with biological fractionation in FAMOUS, we tested three different parameterisations for α POC←aq . In the standard simulation (std), we calculated α POC←aq according to Popp et al. (1989): where CO * 2 is the aqueous CO 2 concentration (in µmol L −1 ). Both of the alternative parameterisations calculated α POC←aq as a function of the phytoplankton-specific growth rate (µ) and CO * 2 , representing an increase in complexity relative to the standard scheme. The first was a linear relationship derived from the experimental results of Laws et al. (1995): The second was a non-linear relationship derived from the experimental results of Laws et al. (1997): Because HadOCC is a relatively simple ecological model, with only a single representation of phytoplankton, we did not test more complex schemes, such as those that use phytoplankton-type-specific cell parameters (e.g. Burkhardt et al., 1999;Keller and Morel, 1999).

Advection
The default advection scheme in FAMOUS is Quadratic Upstream Interpolation for Convective Kinematics (QUICK) with flux limiter (Leonard et al., 1993). This scheme is used to compute the transport of tracers such as temperature, salinity, nutrients, and DIC throughout the ocean. For consistency, we use the same advection scheme to calculate 13 C concentrations in the ocean interior. For greater numerical stability, δ 13 C DIC is fixed at 0 ‰ in the Hudson Bay and Baltic Sea. With the model's standard pre-industrial land-sea mask, these inland bodies of water are isolated from the global oceans; therefore, their isotope concentrations will not affect large-scale tracer distributions.

Spin-up simulation
Carbon isotope simulations must be spun up over multiple millennia (5000 to 15 000 years; Orr et al., 2000) to reach steady state because of the long timescale of deep ocean ventilation (Bardin et al., 2014). We therefore ran our spinup simulation for 10 000 years with constant pre-industrial boundary conditions, where δ 13 C atm was fixed at −6.5 ‰ (Francey et al., 1999) and δ 13 C ocn was initialised at a globally uniform value of 0 ‰. The global volume-weighted integral of δ 13 C DIC started to stabilise after 7000 years, and at the end of the spin-up simulation, the drift was less than 0.001 ‰ yr −1 (Fig. S1 in the Supplement).

Historical simulation
A transient simulation for the period 1765 to 2000 CE was initialised from the end of the spin-up simulation to generate model output that is directly comparable to modern observations (Fig. 3). Atmospheric CO 2 concentrations were prescribed from the OCMIP-2 files (Orr et al., 2000) and δ 13 C atm was prescribed from the Law Dome and South Pole ice core records (Rubino et al., 2013). The decrease in δ 13 C atm from −6.5 ‰ in 1750 to approximately −8.0 ‰ in 2000 is due to the Suess effect. First observed in tree ring records of atmospheric composition, the Suess effect refers to the dilution of 13 C in any carbon pool due to fossil fuel burning (Suess, 1955;Keeling, 1979). Fossil fuels formed millions of years ago from organic matter, which is relatively 13 C depleted due to isotopic fractionation during photosynthesis. Their isotopic signature is therefore approximately 20 ‰ lower than that of the ambient atmosphere (Andres et al., 1994(Andres et al., , 1996. To act as a control, the spin-up simulation  (Rubino et al., 2013) and prescribed atmospheric CO 2 values (dashed) from the OCMIP-2 files (Orr et al., 2000).
was continued for an additional 235 years with constant CO 2 and δ 13 C atm .

Sensitivity experiments
Five further simulations were conducted to quantify the effects of fractionation during air-sea gas exchange and primary productivity on the simulated δ 13 C DIC distributions. All five simulations were run for 10 000 years with constant preindustrial boundary conditions. In each of the simulations, δ 13 C atm was fixed at −6.5 ‰ and δ 13 C ocn was initialised at 0 ‰. At the end of each of the spin-up simulations, the global volume-weighted δ 13 C DIC integral was drifting by less than 0.001 ‰ yr −1 . Three of the simulations were designed to quantify the effects of the individual processes outlined in Sect. 2.2 (Table 2). In the ki-fract-only simulation, α aq←g , α DIC←g , and α p were all set to 1; therefore, only kinetic fractionation effects were calculated. In the no-asgx-fract simulation, α k , α aq←g , and α DIC←g were all set to 1 to eliminate the effect of fractionation during air-sea gas exchange. Fractionation during photosynthesis continued to be calculated using the std biological fractionation scheme, as per Eqs. (9)-(10). In the no-bio-fract simulation, α p was set to 1 to remove the effect of fractionation during photosynthesis, but fractionation during air-sea gas exchange continued to be calculated as per Eqs. (6)-(8).
The other two simulations were designed to assess the sensitivity of the simulated δ 13 C DIC distributions to the choice of biological fractionation scheme (Sect. 2.2.2). In the L95 simulation, α POC←aq was calculated using Eq. (11), whilst in the L97 simulation, α POC←aq was calculated using Eq. (12). As with the std simulation, we initialised a 235-year transient simulation (with the 13 C-Suess effect) from the end of both of these spin-ups to allow the output from all three photo- synthetic fractionation schemes to be compared directly to observations.
3 Results and discussion

Validating the isotope scheme
Isolating the different fractionation effects allows us to assess the relative contribution of air-sea gas exchange and biology to the simulated δ 13 C DIC distributions, and validate that the new isotope scheme is responding to physical and biogeochemical processes as expected. If there is no fractionation during either air-sea gas exchange or photosynthesis, the ocean equilibrates at a uniform value of −6.5 ‰, in line with the atmosphere (simulation not shown). Kinetic fractionation has only a minor effect on surface ocean δ 13 C DIC distributions, with simulated δ 13 C DIC values in the ki-fract-only simulation ranging between −6.57 ‰ in the Labrador Sea and −6.42 ‰ in the eastern equatorial Pacific (Fig. 4a). This represents a −0.07 ‰ to +0.08 ‰ shift relative to no isotopic fractionation. Specifically, there is 13 C depletion (low δ 13 C DIC ) in areas of net CO 2 invasion, such as the extratropics and high latitudes, and 13 C enrichment (high δ 13 C DIC ) in the equatorial upwelling zones and the deep water formation regions where CO 2 is being outgassed. Kinetic fractionation has a negligible effect on the δ 13 C DIC depth profile, with globally averaged δ 13 C DIC values of −6.4955 ‰ in the surface ocean and −6.5011 ‰ in the abyssal ocean (Fig. 5). When both the equilibrium and kinetic fractionation effects are included during air-sea gas exchange (no-bio-fract), the large-scale δ 13 C DIC distributions are closely related to the SST patterns because of the temperature dependence of α aq←g and α DIC←g (Fig. 4b). In the absence of biological fractionation, relatively high δ 13 C DIC values (> +2.5 ‰) are simulated in the Southern Ocean due to the combined effect of CO 2 outgassing and low SSTs, both of which cause 13 C enrichment. The δ 13 C DIC values in the Arctic Ocean are comparably low because the model has more extensive sea ice in the Northern Hemisphere than in the Southern Hemisphere, which inhibits air-sea gas exchange. The highest values (+3.00 ‰) are simulated in the eastern equatorial Pacific where there are high rates of net CO 2 outgassing, and southern-sourced waters, which have a high δ 13 C DIC signature in this simulation because there is no biological fractionation, are upwelled. Low δ 13 C DIC values are simulated in the Indian Ocean, with the lowest values (+1.1 ‰) in southeast Asia, because the sea surface is warmer than at the equivalent latitudes in the Atlantic and Pacific oceans. The globally averaged δ 13 C DIC values in this simulation range between +2.03 ‰ in the surface ocean and +2.16 ‰ in the deep ocean, with a minimum value of +2.00 ‰ at a depth of approximately 200 m (Fig. 5). Below approximately 1500 m, the globally averaged δ 13 C DIC is near constant with depth, matching the simulated temperature profile.
When only biological fractionation effects are included (no-asgx-fract), δ 13 C DIC values in the surface ocean range between −7.65 ‰ in the eastern equatorial Pacific and −3.89 ‰ in the eastern equatorial Atlantic (Fig. 4c), representing a shift of −1.15 ‰ to +2.61 ‰ relative to no isotopic fractionation. The asymmetry between these two upwelling zones occurs because the waters that are being upwelled from the deep Pacific Ocean are approximately 600 years older than the equivalent waters in the Atlantic Ocean. They therefore contain a higher percentage of remineralised organic matter, which is enriched in 12 C. Relatively low δ 13 C DIC values are also simulated in the Southern Ocean and northeast North Atlantic Ocean where older water is mixed upwards from the abyssal ocean to the surface ocean at sites of deep water formation. The globally averaged δ 13 C DIC values in this simulation range between −5.85 ‰ in the productive surface ocean and −7.56 ‰ in the abyssal ocean, with a minimum value of −7.86 ‰ at a depth of approximately 1000 m, which corresponds to the depth of maximum remineralisation in the model (Fig. 5). The values change from greater than −6.5 ‰ (enriched in 13 C relative to no fractionation) to less than −6.5 ‰ (depleted in 13 C relative to no fractionation) at a depth of approximately 100 m, which corresponds to the photic zone.
The spatial patterns in the std simulation and the no-asgxfract simulation are closely matched, both in the surface ocean (Fig. 6) and at depth (Fig. 5), demonstrating the importance of biology to the large-scale δ 13 C DIC distributions. However, in the surface layer, air-sea gas exchange partly compensates for the biological effects in the Southern Ocean, the Northern Hemisphere deep water formation region, and the equatorial upwelling zones, as inferred from the peak surface zonal mean δ 13 C DIC values at 60 • S, 55 • N and 0 • in the no-bio-fract simulation, which correspond with reduced amplitude troughs in the std simulation relative to the noasgx-fract simulation. Similar results pertaining to the relative influence of air-sea gas exchange and biology were presented by Schmittner et al. (2013), who concluded that air-sea gas exchange and temperature-dependent fractionation reduce the spatial δ 13 C DIC gradients that are created by biology. Earlier work by Murnane and Sarmiento (2000) and  Schmittner et al. (2013) also supports the notion that biology is the dominant factor controlling δ 13 C DIC distributions in the interior ocean. Overall, the sensitivity experiments demonstrate that the new carbon isotope scheme is accurately responding to physical and biogeochemical processes in the model, such as temperature, air-sea gas exchange, and the biological pump.

Comparison to observations
To assess the model performance in representing modern large-scale 13 C distributions, we compare the simulated mean δ 13 C DIC values for the 1990s with observations from GLODAP version 2 (v2; Key et al., 2015;Olsen et al., 2016) and the gridded global ocean climatology of Eide et al. (2017). The δ 13 C DIC values in the std simulation are, on average, 0.97 ‰ higher than the GLODAPv2 observations in the surface ocean (Fig. 7) and 0.64 ‰ higher globally, with root mean square error (RMSE) values of 1.03 ‰ and 0.91 ‰, respectively. However, the simulated range in the surface ocean (3.2 ‰) is in excellent agreement with the observed range (3.3 ‰). Specifically, the simulated surface Figure 6. Zonally averaged mean annual surface δ 13 C DIC at the end of the sensitivity experiment spin-up simulations (years 9900 to 10 000). The std (black) and no-bio-fract (purple) simulations use the left-hand axis, whilst the ki-fract-only (red) and no-asgx-fract (blue) simulations use the right-hand axis. The dotted lines are the equivalent simulations conducted by Schmittner et al. (2013) with the UVic ESM: std (black) and no-bio (purple) on the left-hand axis; ki-only (red) and const-gasx (blue) on the right-hand axis.
Re-examining the results of the sensitivity experiments allows us to ascertain the underlying causes of the model-data discrepancy. Schmittner et al. (2013; herein S13) conducted a similar set of simulations with the UVic ESM to elucidate the relative influence of biology and air-sea gas exchange on the distribution of oceanic δ 13 C DIC (see Table 1 in S13). Overall, there is good agreement between our ki-fract-only and nobio-fract simulations and the equivalent simulations in S13 (ki-fract and no-bio, respectively), both in the surface ocean (Fig. 6) and at depth (Fig. 5). However, there is a clear difference between the results of our no-asgx-fract simulation and the equivalent simulation in S13 (const-gasx). Specifically, the surface ocean zonal mean δ 13 C DIC values in our no-asgx-fract simulation range between −6.6 ‰ at 60 • S and −5.5 ‰ in the subtropics, with a local minimum of −5.8 ‰ at the Equator (solid blue line in Fig. 6). For comparison, the surface ocean zonal mean values in const-gasx range between −8.0 ‰ in the Southern Ocean and −5.75 ‰ in the Southern Hemisphere subtropics, with a localised minimum of −6.25 ‰ at the Equator (dotted blue line in Fig. 6). Similarly, whilst the globally averaged deep ocean δ 13 C DIC values in our no-asgx-fract simulation have a comparable range (2.01 ‰) to the deep ocean values in const-gasx, there is an offset of approximately 1 ‰, with S13 simulating δ 13 C DIC values of −6.4 ‰ in the surface ocean, −8.4 ‰ in the deep ocean, and near constant values below 1000 m (blue lines in Fig. 5). Overall, the δ 13 C DIC values in the standard simulation with the UVic ESM are in good agreement with observations, with a global linear regression r 2 value of 0.91 and a global RMSE of 0.33 ‰ (Schmittner et al., 2013;Buchanan et al., 2019). We therefore postulate that the offset in the simulated δ 13 C DIC values in FAMOUS relates to biases in the biological carbon cycle.
Elucidating the exact cause of δ 13 C DIC model-data discrepancy is difficult. There are a number of fluxes going into and out of the DI 13 C pool (Fig. 2), each of which could have biases that are compounding or reducing the overall δ 13 C DIC bias. For example, if any of the rates of phytoplankton respiration, phytoplankton mortality, or zooplankton mortality are too low, the input of 12 C-enriched material back into the DIC pool would be insufficient. Similarly, if the model is not simulating enough remineralisation, either as a direct consequence of the parameterised remineralisation rate or as a result of insufficient POC export, the input of 12 C-enriched material back into the DIC pool would again be too low.
Primary producers preferentially take up 12 C during photosynthesis; therefore, higher-than-observed rates of net primary production in the photic zone would increase δ 13 C DIC . However, if the δ 13 C DIC discrepancy in FAMOUS was a simple function of the biases in net primary production, δ 13 C DIC would be lower than observed in the subtropical gyres, the Indian Ocean, and the northern North Atlantic and North Pacific oceans, and higher than observed in the equatorial upwelling zones and the Southern Ocean (Fig. 1). Thus, whilst the differences in net primary production could be contributing towards the δ 13 C DIC bias, particularly in the equatorial upwelling zones, they alone cannot explain the unidirectional offset.
Alternatively, the fractionation during photosynthesis could be too strong as a result of imbalances in the carbonate chemistry (Fig. S2). The global mean alkalinity in FAMOUS is 81 µmol kg −1 higher than observed and the mean alkalinity in the uppermost 50 m of the ocean is 107 µmol kg −1 too high (Key et al., 2004;Sarmiento and Gruber, 2006). In addition, the simulated global mean DIC concentration is 54 µmol kg −1 higher than observed and the mean DIC concentration in the uppermost 50 m of the ocean is 96 µmol kg −1 too high (Key et al., 2004;Sarmiento and Gruber, 2006). Furthermore, the mean ocean temperatures in FAMOUS are warmer than observed, both globally (2.2 • C) and in the uppermost 50 m of the ocean (1 • C; Sarmiento and Gruber, 2006;Locarnini et al., 2013). Increasing alkalinity increases CO * 2 , whilst increasing the temperature and DIC concentrations decreases CO * 2 . Hence, the overall effect of the carbonate chemistry biases in FAMOUS result in the global mean CO * 2 being 3.03 µmol L −1 too low and the mean CO * 2 in the uppermost 50 m of the ocean being 0.58 µmol L −1 too high. In the photic zone, this translates to a simulated α p of 0.97378 compared to an observed α p of 0.97415 using the std fractionation parameterisation. Thus, we postulate that imbalances in the carbonate chemistry, and the consequent differences in α p , are contributing towards the δ 13 C DIC bias, but the overall effect is small.
The smallest model-data discrepancies in the surface layer are in the Southern Ocean and the northeast North Atlantic Ocean where deep convection mixes 12 C-enriched waters upwards (Fig. 7). In contrast, in the equatorial upwelling zones, the effect of higher than observed primary productivity (increasing δ 13 C DIC ) outweighs the effect of vertical mixing (reducing δ 13 C DIC ); therefore, the overall model-data biases are higher in these regions. Despite the global offset, the model correctly simulates lower δ 13 C DIC values in the Indian Ocean compared to the Atlantic and Pacific oceans, because the Indian Ocean is relatively nutrient poor, both in the model and reality (Fig. S3), which limits primary productivity (Fig. 1). Similar to previous 13 C modelling studies (e.g. Hofmann et al., 2000;Tagliabue and Bopp, 2008;Schmittner et al., 2013), FAMOUS also accurately simulates the observed latitudinal gradient in mixed layer δ 13 C POC , with relatively high val-ues (≈ −20 ‰) in the low latitudes and relatively low values (≈ −27 ‰) at high latitudes (Fig. 8).
As observed, δ 13 C DIC decreases with depth in all basins due to the remineralisation of isotopically light organic matter (Fig. 9). The maximum remineralisation depth in the model is approximately 1000 m, which is 200 to 500 m shallower than observed. In the deep ocean, the highest δ 13 C DIC values are in the Atlantic basin, with intermediate values in the Indian basin, and the lowest values in the Pacific basin, where the waters are older and therefore contain more remineralised organic material (enriched in 12 C). However, there are notable structural differences in the zonal means (Fig. 10), which arise due to inaccuracies in the physical ocean circulation in FAMOUS. Specifically, FAMOUS does not capture the observed structure in the Atlantic basin because, in this generation of the model, the AMOC is characterised by an overdeep North Atlantic Deep Water (NADW) cell and insufficient Antarctic Bottom Water formation (Dentith et al., 2019b). FAMOUS also simulates weak (less than 3 Sv) ventilation to depths of 2000 m in the North Pacific Ocean (Dentith et al., 2019b), which prevents the accumulation of old, 12 C-enriched (low δ 13 C DIC ) waters at inter- Figure 8. Zonally averaged mean annual mixed layer δ 13 C POC during the 1990s: observations (Goericke and Fry, 1994; red), the std simulation (black), the L95 simulation (grey), and the L97 simulation (blue). mediate depths in the Northern Hemisphere. Instead, the oldest carbon in the model is in the eastern equatorial Pacific. In addition, the surface winds in the model are weaker than observed (Kalnay et al., 1996), resulting in a relatively shallow mixed layer. This promotes the excessive accumulation of high δ 13 C DIC values in the surface ocean, which is particularly notable in the Southern Hemisphere subtropical gyres. These physical model biases are also clearly visible in the zonal mean profiles of other tracers, such as nutrients (Fig. S4) and DIC (Fig. S5). The overall shape of the simulated depth profile reaffirms the notion that there are inaccuracies in both the physical and biogeochemical components of the model (Fig. 9). Below approximately 1000 m, the simulated δ 13 C DIC values increase with depth in each ocean basin, whilst the observed basin averages are near constant with depth. The offset between the simulated and observed values is greatest in the deep Atlantic Ocean, where too much 13 C-enriched water from the shallow ocean is being circulated into the abyssal ocean. However, the trend towards increasing δ 13 C DIC with depth could also be in part explained by insufficient remineralisation in the model. This is supported by lower-than-observed nutrient concentrations in the deep ocean (Fig. S4). HadOCC's global export ratio at 2000 m is within the observed range, but a lack of spatial variation means that the geographic distributions are partially incorrect (Palmer and Totterdell, 2001). Hence, we postulate that localised inaccuracies in the export ratio, together with deficiencies in the parameterisation of the remineralisation rate, are contributing towards the δ 13 C DIC offset. The basin-averaged δ 13 C DIC bias is smallest in the Pacific Ocean, where the waters are old and therefore have had more time to remineralise, thereby partially compensating for the biogeochemical biases. Indeed, the shape of the simulated and observed basin-averaged depth profiles are in good agreement below approximately 2000 m in the Pacific Ocean, despite the structural differences in the zonal mean.
As outlined in Sect. 3.1, our carbon isotope implementation is sensitive to physical and biogeochemical processes in the model. Thus, whilst biases in the overturning circulation and the biological pump are currently limiting the model's ability to accurately represent modern large-scale 13 C distributions, the model-data agreement could be improved if the physical and ecological components of FAMOUS were recalibrated. This will be discussed further in Sect. 3.4.

Biological fractionation parameterisations
Given the uncertainty associated with biological fractionation (Sect. 2.2.2), we tested three different parameterisations for equilibrium fractionation during photosynthesis. For all three parameterisations, the total fractionation during photosynthesis is greatest in the high latitudes (where SSTs are relatively low and CO * 2 is relatively high) and lowest in the equatorial regions (where SSTs are relatively high and CO * 2 is relatively low; Fig. 11). The std parameterisation produces the largest range of α p values (between approximately 0.97 and 0.98), whilst the L95 parameterisation produces the smallest range (between approximately 0.964 and 0.970). The total fractionation during photosynthesis increases with the complexity of the parameterisation, with L97 producing the largest overall effect (with a minimum α p of 0.9635). For all three parameterisations, α p decreases (i.e. the strength of fractionation increases) with depth in the photic zone, with the largest gradient produced by the std parameterisation (Fig. S6).
The large-scale δ 13 C DIC patterns are very similar for all three photosynthetic fractionation schemes, but the parameterisations that take the phytoplankton growth rate in account simulate higher surface ocean δ 13 C DIC values everywhere except in the Southern Ocean, the Nordic Seas, and the eastern equatorial regions, where older 13 C-depleted waters are mixed upwards from the abyssal ocean during deep water formation and upwelling (Fig. S7). The differences are amplified when using the L97 parameterisation (RMSE = 1.24 ‰, bias = 1.15 ‰), which specifies a non-linear relationship between µ and CO * 2 , compared to the L95 parameterisation (RMSE = 1.21 ‰, bias = 1.13 ‰), which specifies a linear relationship (Fig. S8). Conversely, the alternative parameterisations decrease δ 13 C DIC at depth compared to the std simulation, bringing the simulated values closer to the observations (Fig. 9). Below approximately 500 m depth, the δ 13 C DIC values are consistently lower when using the L97 parameterisation compared to the L95 parameterisation. This is due to the preconditioning of δ 13 C DIC and δ 13 C POC as  Eide et al. (2017). The simulated values have also been subsampled at the locations where there is a corresponding observation in the GLODAPv2 dataset (Key et al., 2015;Olsen et al., 2016;dashed). The red shading shows the estimated uncertainty in δ 13 C DIC observations due to unresolved intercalibration between different laboratories (±0.2 ‰; Schmittner et al., 2013;Eide et al., 2017). a result of fractionation during photosynthesis in the photic zone. In the L95 and L97 simulations, δ 13 C POC is lower than in the std simulation due to increased uptake of 12 C during primary production (lower α p ). The latitudinal δ 13 C POC gradients in the mixed layer in these simulations are lower than observed, with zonal mean values ranging between approximately −30 ‰ at the Equator and −33 ‰ at 60 • N/S (Fig. 8). When the POC is remineralised, a relatively low δ 13 C signal is therefore being released back into the DIC pool, which causes the δ 13 C DIC in the deep ocean to be lower than in the std simulation. Thus, although the rates of biological exchange and overturning circulation are the same in all three simulations, the preconditioning of δ 13 C DIC and δ 13 C POC in the photic zone creates differences between the three simulations at depth. Whilst the global RMSE compared to the GLODAPv2 dataset is lower in the L95 and L97 simulations (0.86 ‰ and 0.87 ‰, respectively), it is still almost double the RMSE in other models (Buchanan et al., 2019). Overall, increasing the complexity of the fractionation scheme does not significantly improve the model-data agreement because of the aforementioned physical and biogeochemical biases.

A new tuning target
In contrast with earlier studies, we have demonstrated that the new carbon isotope scheme in FAMOUS is sensitive to both physical and biogeochemical processes. The simulated δ 13 C DIC distributions reflect known physical inaccuracies (such as overdeep NADW and weak convection in the subpolar North Pacific Ocean) and have allowed us to identify previously undisclosed biogeochemical biases (e.g. in the representation of remineralisation). The new tracer there-fore offers excellent potential as a holistic tuning target for recalibrating FAMOUS in the future.
FAMOUS has previously been tuned both systematically (Jones et al., 2005;Gregoire et al., 2011;Williams et al., 2013) and manually (Smith et al., 2008). Most recently, Williams et al. (2013) tuned the 20 structural parameters in HadOCC (coupled to FAMOUS-MOSES2.2) using an objective hypercube technique. Specifically, the parameter set included the C : N ratios for the different carbon pools, phytoplankton-specific parameters (e.g. maximum rate of photosynthesis), zooplankton-specific parameters (e.g. linear and quadratic zooplankton mortality rates), detritus-specific parameters (e.g. shallow and deep remineralisation rates), and carbonate-specific parameters (e.g. calcite export ratio). The main diagnostics used to evaluate the performance of the ensemble members were December-January-February and June-July-August surface air temperatures, annual mean total precipitation rate, annual mean nitrate concentrations, and primary productivity. Crucially, this study only ran each perturbed parameter simulation for 200 years and neglected to evaluate the strength and structure of the AMOC. The optimal parameter set therefore had small but important imbalances in the surface climate, which caused the AMOC to collapse over longer (multi-millennial) timescales (Dentith et al., 2019b).
HadOCC has not yet been tuned for the configuration of the model used in our study (FAMOUS-MOSES1). Simultaneously recalibrating HadOCC and the physical ocean circulation in FAMOUS-MOSES1 could therefore improve the simulated δ 13 C DIC distributions. In the first instance, further sensitivity studies would provide more insight into the extent to which our results could be improved by small adjustments Figure 10. Zonal mean δ 13 C DIC during the 1990s in the Atlantic Ocean (left), Pacific Ocean (centre), and Indian Ocean (right): (a-c) gridded observations (Eide et al., 2017), (d-f) the std simulation, (g-i) the L95 simulation, and (j-l) the L97 simulation.
to the model's biogeochemistry (e.g. modifying the remineralisation rate and/or the export ratio). Longer term, we propose that the addition of carbon isotope ratios as tuning targets (both the δ 13 C presented here and the new 14 C tracer for FAMOUS discussed by Dentith et al., 2019a) would improve the work of Williams et al. (2013) because they provide an objective and straightforward way of assessing whether the balance between all of the ecological processes in the model is correct.

Summary
We have added the stable isotope 13 C to the ocean component of the FAMOUS GCM, using the MOSES1 generation of the model to validate our scheme. We account for fractionation during air-sea gas exchange and photosynthesis, and have tested three different parameterisations for the latter. The model captures the range of observed δ 13 C DIC values in the surface ocean, but the simulated values are approximately 1 ‰ too high at all depths. The differences between the three fractionation schemes are relatively minor, but when fractionation during photosynthesis accounts for phytoplankton growth rates as opposed to just aqueous CO 2 concentrations the discrepancies between the model and observations are further increased in the surface ocean and reduced at depth. The sensitivity experiments suggest that the simulated values are too high because of underlying biases in the biological carbon cycle; therefore, retuning HadOCC could improve the model-data agreement. Biases in the large-scale ocean circulation also inhibit the model's ability to accurately sim-ulate the large-scale distribution of tracers in the deep ocean. Retuning the ocean circulation to improve the representation of the AMOC, in particular, would further reduce the model-data discrepancies. Thus, our results emphasise the utility of implementing carbon isotopes in GCMs; the simulated isotope distributions provide an additional measure against which the physical and biogeochemical model performance can be evaluated and offer an extra tuning metric for prospective development work. In the future, we intend to use the isotope-enabled model to study ocean circulation and the marine carbon cycle in both a modern and palaeo context, for example, at the Last Glacial Maximum (21 000 years ago) and during the last deglaciation (21 000 to 11 000 years ago). The standard equation for calculating the virtual flux to account for the dilution or concentration effect of surface freshwater fluxes is where E is evaporation, P is precipitation, and dz is layer depth.
As we carry 13 C as a ratio ( 13 C/ 12 C), virtual fluxes are not required: Appendix B: Air-sea gas exchange equations The standard equation for calculating the change in DI 13 C due to air-sea gas exchange is where PV is the piston velocity (Eq. 4), C sat is the saturation concentration of atmospheric CO 2 (in mol m −3 ), C surf is the surface aqueous concentration of CO 2 (in mol m −3 ), α k is the constant kinetic fractionation factor, α aq←g is the temperature-dependent fractionation during gas dissolution (Eq. 7), α DIC←g is the temperature-dependent fractionation between aqueous CO 2 and DIC (Eq. 8), and 13 A/ 12 A and 13 C/ 12 C are the 13 C/ 12 C ratios of the atmosphere and DIC, respectively. The equation for calculating the change in DI 13 C/DI 12 C due to air-sea gas exchange is

Appendix C: Biological equations
For consistency with the standard biological tracers, the 13 C contents of phytoplankton ( 13 P), zooplankton ( 13 Z), and detritus ( 13 D) are carried in mmol N m −3 , with the carbon concentrations and fluxes calculated using fixed stoichiometric ratios. The DI 13 C/DI 12 C values are therefore converted from a ratio in model units (Eq. 2) to normalised DI 13 C concentrations before entering the soft tissue pump. The conversion is reversed at the end of each time step.

C1 Phytoplankton (P)
The standard equation for calculating the change in phytoplankton ( 12 P) is where R P is the specific growth rate of phytoplankton, G P represents grazing by zooplankton, m P represents phytoplankton mortality due to overpopulation, and η P represents phytoplankton respiration. The equation for calculating the change in 13 P is where α p is the isotopic fractionation that occurs during photosynthesis (Eq. 9), 13 C/ 12 C is the 13 C/ 12 C ratio of DIC, and 13 P/ 12 P is the 13 C/ 12 C ratio of phytoplankton.
J. E. Dentith et al.: Simulating stable carbon isotopes with FAMOUS The 13 P tracer is updated using the forward Euler method: The standard equation for calculating the change in zooplankton ( 12 Z) is where β P and β D are the assimilation efficiencies associated with zooplankton grazing on phytoplankton (G P ) and detritus (G D ), respectively, and m Z represents zooplankton mortality due to predation and natural causes. The equation for calculating the change in 13 Z is where 13 P/ 12 P, 13 D/ 12 D, and 13 Z/ 12 Z are the isotopic ratios of phytoplankton, detritus and zooplankton, respectively. The 13 Z tracer is updated using the forward Euler method: (C6)

C3 Dissolved inorganic carbon (DIC, C)
The standard equation for calculating the change in DI 12 C is where R P is the specific growth rate of phytoplankton, λ D is detrital remineralisation, which is specified at a constant rate (0.1 d −1 in the uppermost 250 m of the ocean and 0.02 d −1 at all other depths), β P and β D are the assimilation efficiencies associated with zooplankton grazing on phytoplankton (G P ) and detritus (G D ), respectively, m Z represents zooplankton mortality due to predation and natural causes, m P represents phytoplankton mortality due to overpopulation, and η P represents phytoplankton respiration. The equation for calculating the change in DI 13 C is d 13 C dt = −R P × 13 C 12 C × α p + λ D × 13 D 12 D + (1 − β P ) × G P × 13 P 12 P + (1 − β D ) × G D × 13 D 12 D + m Z × 13 Z 12 Z + m P × 13 P 12 P + η P × 13 P 12 P , where α p is the isotopic fractionation that occurs during photosynthesis (Eq. 9), and 13 C/ 12 C, 13 D/ 12 D, 13 P/ 12 P, and 13 Z/ 12 Z are the isotopic ratios of DIC, detritus, phytoplankton, and zooplankton, respectively. The DI 13 C tracer is updated using the forward Euler method: 13 C (t+ t) = 13 C (t) + t × −R P (t) × 13 C 12 C (t) × α p (t) + λ D (t) × 13 D 12 D (t) + 1 − β P (t) × G P (t) × 13 P 12 P (t) +m P (t) × 13 P 12 P (t) + η P (t) × 13 P 12 P (t) . (C9)

C4.1 Biological effects
The standard equation for calculating the change in detritus ( 12 D) due to biology is where m Z represents zooplankton mortality due to predation and natural causes, m P represents phytoplankton mortality due to overpopulation, λ D is detrital remineralisation, which is specified at a constant rate (0.1 d −1 in the uppermost 250 m of the ocean and 0.02 d −1 at all other depths), and β P and β D are the assimilation efficiencies associated with zooplankton grazing on phytoplankton (G P ) and detritus (G D ), respectively.
The equation for calculating the change in 13 D due to biology is where 13 Z/ 12 Z, 13 P/ 12 P, and 13 D/ 12 D are the isotopic ratios of zooplankton, phytoplankton, and detritus, respectively.

C4.2 Reflux
The small amount of detritus that reaches the ocean floor is immediately refluxed back to the surface layer to conserve nitrogen and carbon. dD dt sink_in (k=1) = γ dz/100 × D (k=KMT) , where k is the model level, γ is the sinking rate, which is parameterised at 10 m d −1 , dz is the depth of the layer (in cm), D is the detritus concentration, and KMT is the maximum depth of the ocean. The same equation applies for 13 D.
Code availability. FAMOUS can be obtained from http://cms. ncas.ac.uk/wiki/UmFamous (last access: 23 December 2019). The code detailing the advances described in this paper is available via the Research Data Leeds Repository (Dentith, 2019, https://doi.org/10.5518/621) under a Creative Commons Attribution 4.0 International (CCBY 4.0) license. These files are known as code modification ("mod") files and should be applied to the original model code, which can be viewed online at http://cms. ncas.ac.uk/code_browsers/UM4.5/UMbrowser/index.html (last access: 23 December 2019). All of the additional modification files that are required to run the simulations described in this paper are available in the Supplement. These standard FAMOUS updatessome of which have been described by Smith et al. (2008), Smith (2012), and Valdes et al. (2017) -and the original model code are protected under UK Crown Copyright. The UM configuration ("basis") files for the simulations described in this paper are also available in the Supplement.
The simulations described in this study, as denoted by their unique five-letter Met Office UM identifiers and the notation used within this paper, are as follows: Supplement. The supplement related to this article is available online at: https://doi.org/10.5194/gmd-13-3529-2020-supplement.
Author contributions. RFI designed and supervised the project. JED wrote and implemented the code with input from JCT, LJG, and RFI. JED ran the simulations, analysed the results, and prepared the manuscript with input from all co-authors.
Competing interests. The authors declare that they have no conflict of interest.