Including an ocean carbon cycle model into i LOVECLIM ( v 1 . 0 )

The atmospheric carbon dioxide concentration plays a crucial role in the radiative balance and as such has a strong influence on the evolution of climate. Because of the numerous interactions between climate and the carbon cycle, it is necessary to include a model of the carbon cycle 5 within a climate model to understand and simulate past and future changes of the carbon cycle. In particular, natural variations of atmospheric CO2 have happened in the past, while anthropogenic carbon emissions are likely to continue in the future. To study changes of the carbon cycle and climate on 10 timescales of a few hundred to a few thousand years, we have included a simple carbon cycle model into the iLOVECLIM Earth System Model. In this study, we describe the ocean and terrestrial biosphere carbon cycle models and their performance relative to observational data. We focus on the 15 main carbon cycle variables including the carbon isotope ratios δ13C and the ∆14C. We show that the model results are in good agreement with modern observations both at the surface and in the deep ocean for the main variables, in particular phosphates, DIC and the carbon isotopes. 20


Introduction
The carbon cycle is a key component of climate and environmental sciences, both because CO 2 is a greenhouse gas (Tyndall, 1861) and has a direct impact on climate, but also because it plays an important role in ocean acidification (Orr et al., 2005) which directly impacts marine life. The three main carbon reservoirs involved on the timescale of a few thousand years are the atmosphere, the ocean and the land biosphere. The ocean is the biggest of the three reser-voirs with around 39 000 GtC, while the atmosphere contains around 589 GtC and the terrestrial biosphere between 1950 and 3050 GtC for the pre-industrial . The climate also impacts the carbon cycle and hence the concentration of atmospheric CO 2 through various dynamical, chemical and biological processes. For example, changes in the ocean temperature will modify the solubility of CO 2 : the warmer the ocean the less soluble CO 2 becomes, which decreases the carbon stock in the ocean and increases atmospheric CO 2 . Temperature, as well as humidity, also influences the development of the terrestrial biosphere and decomposition of terrestrial organic matter. Low temperature and dry conditions tend to favour lower rates of decomposition. The various climate-carbon interactions involve all three carbon reservoirs. Therefore it is necessary to include a model of the carbon cycle within a climate model to understand past changes and anticipate the future evolution of the carbon cycle and climate.
Such models have been developed during the last decades (Cox et al., 2000;Friedlingstein et al., 2001) and a subset of coupled models used in CMIP5 (Coupled Model Intercomparison Project 5) now include a complete description of the ocean and land carbon cycles. Eleven models have been compared within the framework of the Fourth Coupled Carbon Cycle Climate Model Intercomparison Project (C4MIP) (Friedlingstein et al., 2006). They include models of both the ocean and the land carbon cycle.
Climate models range from simple box models to global climate models (GCMs). The carbon models have gradually become more complex by including more types of plankton in the ocean and more plant functional types on land, as well as more nutrients, such as iron in the ocean or nitrogen on land (Anav et al., 2013). The number of addi-tional tracers directly impacts the computing time, therefore such complex models are well suited to study the climatecarbon evolution on timescales of a few decades to hundreds of years, but are too computationally expensive for longer simulations. Simpler carbon models such as the ocean carbon models based on NPZD (nutrient-phytoplanktonzooplankton-detritus) ecosystems, and simple terrestrial biosphere models with a few plant functional types, associated with intermediate-complexity climate models, are thus more convenient for the study of long timescales of more than a few thousand years.
Intermediate-complexity models are well suited for long term studies of a few thousand to hundred of thousand years, and in particular the glacial-interglacial cycles. The carbon cycle varies greatly during the glacial and interglacial periods, with atmospheric CO 2 concentrations of around 190 ppm during the relatively colder glacials periods and around 280 ppm during the warmer interglacials (EPICA community members, 2004). Although such periods have already been studied with intermediate-complexity models (Brovkin et al., 2007;d'Orgeville et al., 2010;Bouttes et al., 2010;Tschumi et al., 2011;Menviel et al., 2012), large uncertainties remain concerning the processes responsible for the changes of the carbon cycle.
Besides understanding and simulating CO 2 concentrations in the past and future, the carbon cycle also provides indirect yet valuable information about changes of the ocean dynamics and biology, as well as the land vegetation, through carbon isotopes changes (Duplessy et al., 1988;Crowley, 1995). Indeed, there are no direct data of ocean circulation changes in the past (except for the last decades, see for example Mielke et al., 2013), but the measurement of δ 13 C and 14 C in sediment cores can help constrain the ocean and land vegetation changes. Moreover, the measure of atmospheric δ 13 C in ice cores (Lourantou et al., 2010;Schmitt et al., 2012) and the calibration curves of atmospheric 14 C (Reimer et al., 2013(Reimer et al., , 2009 provide additional data and constraints. By explicitly simulating the carbon isotopes within the carbon cycle model, as we have done in the iLOVECLIM model, it is possible to directly compare model results with data to calibrate model simulations and improve our understanding. Our long-term objective is to study past and future carbon cycle changes over timescales of a few thousand to hundred of thousand years, typical of glacial-interglacial changes. The iLOVECLIM model is perfectly suited to such studies since it includes the relevant physical and dynamical components of the ocean, atmosphere and terrestrial biosphere while running fast enough to simulate thousands of years in a reasonable amount of time (500 simulated years per day). To avoid increasing the computing time excessively, the ocean carbon cycle that we included in iLOVECLIM is based on a NPZD ecosystem which provides the main mechanisms relevant on the timescales of hundreds to thousands years, and includes the carbon isotopes. Sedimentary processes would also be relevant to such timescales. However, the introduction of a sediment model is beyond the scope of this study and remains to be done in future work. The terrestrial biosphere already included in iLOVECLIM has been further developed to add the carbon pools and carbon isotopes. Here, we evaluate the results obtained by including the model of ocean carbon in iLOVECLIM. We focus on the main variables from the carbon cycle and on the ocean carbon isotopes (δ 13 C and 14 C).

iLOVECLIM
The iLOVECLIM model is a new development branch (code fork) of the LOVECLIM model in its version 1.2, as presented in Goosse et al. (2010). It is identical to the latter with respect to its base components: atmosphere, ocean and vegetation (AOV). It has been modified in a number of aspects to include water oxygen isotopes (Roche, 2013) and an interactive ice sheet model (Roche et al., 2014). The general goal of the new developments within iLOVECLIM is to include the suite of processes needed for climate simulations on the Milankovic timescale. We summarize in the following the main characteristics of the AOV components as described in Roche et al. (2007) and Goosse et al. (2010). The following paragraph is taken from Roche et al. (2014).
The atmospheric component ECBilt was developed at the Dutch Royal Meteorological Institute (KNMI) (Opsteegh et al., 1998). Its dynamical core is based on quasi-geostrophic approximation with additional ageostrophic terms added to improve the representation of the Hadley cell dynamics. It is run on a spectral grid with a T21 truncation ( 5.6 • in latitude/longitude in the physical space). ECBilt has three vertical layers at 800, 500 and 200 hPa. Only the first layer contains humidity as a prognostic variable. The time step of integration of ECBilt is 4 h. The oceanic component (CLIO) is a 3-D oceanic general circulation model (Goosse and Fichefet, 1999) based on the Navier-Stokes equations. It is discretized on an Arakawa B-grid at approximately 3 • × 3 • resolution. The vertical discretization follows a "z coordinate" on 20 levels. It has a free surface that allows the use of real freshwater fluxes, a parametrization of downslope currents (Campin and Goosse, 1999) and a realistic bathymetry. CLIO includes a dynamicalthermodynamical sea-ice component that is an updated version of Morales Maqueda (1997, 1999). The dynamic land vegetation model (VECODE) was specifically designed for longterm computation and coupling to coarse resolution models (Brovkin et al., 1997). VECODE con- sists of three sub-models: (1) a model of vegetation structure (bioclimatic classification) calculates plant functional type (PFT) fractions in equilibrium with climate; (2) a biogeochemical model computes net primary productivity (NPP), allocation of NPP, and carbon pool dynamics (leaves, trunks, litter and soil carbon pools); (3) a vegetation dynamics model. The latter computes two plant functional types (PFTs: trees and grass) and a dummy type (bare soil). The vegetation model is resolved on the atmospheric grid (hence at T21 resolution) and allows fractional allocation of PFTs in the same grid cell to account for the small spatial scale needed by vegetation. The different modules exchange heat, stress and water.
For the sake of clarity, it shall be reminded that the carbon cycle model described here does not have any relationship with the LOCH model as described in Goosse et al. (2010).

Carbon cycle in the ocean
The ocean carbon cycle model is originally based on the NPZD ecosystem model described in Six and Maier-Reimer (1996) (Fig. 1). It is the same model as the one included in the CLIMBER-2 model of intermediate complexity (Brovkin et al., 2002a(Brovkin et al., , b, 2007 using the same parameter values, except for the remineralization profile and the atmospheric 14 C, which are described below. The carbon cycle is divided into inorganic and organic parts. The inorganic carbon is simulated as dissolved inorganic carbon (DIC) and alkalinity (ALK). Both tracers are advected and mixed in the ocean by the advectiondiffusion scheme of iLOVECLIM. As in Brovkin et al. (2002a), the flux of carbon at the air-sea surface is computed from the difference between the partial pressure of CO 2 in the atmosphere and ocean (with a gas exchange coefficient of 0.06 mol m −2 yr −1 ). The sea surface pCO 2 is computed from temperature, salinity, DIC and ALK following Millero (1995). The O 2 concentration is prescribed to saturation in the surface cell of the ocean. The organic carbon pool includes six additional tracers on top of inorganic carbon pool, O 2 and the nutrients (phosphate and nitrate, which is diagnostically deduced from phosphate by the Redfield ratio): phytoplankton, zooplankton, dissolved organic carbon (DOC), slow dissolved organic carbon (DOCs), particulate organic carbon (POC) and calcium carbonate (CaCO 3 ). The phytoplankton synthesizes carbon using the light and nutrients available in the first 100 m of the ocean (euphotic zone). It then either dies and sinks or is grazed by zooplankton. Part of the plankton is remineralized to DIC, while part of it is exuded to DOC (and DOCs) and the rest is allocated to POC. The CaCO 3 production is linearly dependent on the organic carbon production with a fixed coefficient. Both POC and CaCO 3 are heavy enough to sink and are instantly remineralized and dissolved at depth. All POC and CaCO 3 are remineralized and dissolved in the water column and there is no riverine input. The remineralization profile follows an exponential law as in Brovkin et al. (2002a), but this profile has been slightly modified to have less remineralization in the upper levels and more below (Fig. 2). All the tracers (except for the particulate pools CaCO 3 and POC) are also transported by the advection-diffusion scheme of iLOVECLIM.

Carbon cycle in the terrestrial biosphere
The VECODE terrestrial biosphere model (Brovkin et al., 1997) was already included in iLOVECLIM (Goosse et al., 2010). The model simulates two types of plants -trees and grass -as well as desert. The plants are divided into four compartments that exchange carbon: leaves, wood, litter and soil. Photosynthesis depends on the local climate (precipitation and temperature) and on the atmospheric CO 2 (CO 2 fertilization). We have added the isotopes of carbon to this pre-existing version of VECODE in every carbon compartment as was done in CLIMBER-2.

Carbon isotopes
Following the original CLIMBER-2 version of the carbon cycle model (Brovkin et al., 2002a(Brovkin et al., , b, 2007, the carbon isotopes 13 C and 14 C are simulated in the ocean and terrestrial biosphere. The 13 C is modelled as in Brovkin et al. (2007), while the numerical code has been modified for the 14 C which is now interactively dependent on cosmogenic production and carbon cycling in the atmosphere instead of having a fixed atmospheric value (Mariotti et al., 2013).
The 13 C simulated in the model is then used as the ratio of 13 C on 12 C to compare to the δ 13 C data from sediment cores. The δ 13 C is defined as follows: R ref is the PDB (Pee Dee Belemnite) carbon isotope standard, which corresponds approximately to average limestone (Craig, 1957). The 13 C distribution in the ocean depends on the airsea exchange, the transport by the ocean circulation (by the advection-diffusion scheme), and the marine biology fractionation. In the terrestrial biosphere, it only depends on the exchange with the atmosphere and the biological fractionation. Indeed, both the marine and terrestrial organisms preferentially use the lighter 12 C over 13 C during photosynthesis, which tends to increase the δ 13 C in the surrounding environment. When the remineralization occurs, the 12 C-rich carbon is released, which decreases the δ 13 C in the atmosphere or ocean.
The 14 C is defined as follows (Stuiver and Polach, 1977): In the model, the simulated 14 C is not subject to any isotopic fractionation (neither biological nor through air-sea exchanges). This formulation allows comparison directly with observations and reconstruction data from the sediment cores that are expressed in 14 C without performing a fractionation correction. The content of 14 C in a reservoir reflects the time since when this reservoir has been in direct contact with the atmosphere. Thus, ocean 14 C gives a good estimate of the age of water masses, which provides useful indications on ocean circulation pathways. This is particularly interesting in palaeoceanography in order to reconstruct past ocean circulation changes. Moreover, the 14 C representation in the model can take into account temporal changes in atmospheric 14 C, which has been the case for example during the historical bomb period or the last deglaciation characterized by changes in the production rate. This aspect of the 14 C representation will thus be particularly useful on future palaeosimulations.

Reference simulation
The model is run under control boundary conditions set to the pre-industrial values for the orbital parameters, ice sheet reconstruction and atmospheric gas concentrations (CO 2 = 280 ppm, CH 4 = 760 ppb and N 2 O = 270 ppb). There are indeed two different CO 2 variables in the model: the CO 2 used for the radiative code and set to 280 ppm, and the one computed by the carbon model. The CO 2 used for the radiative code is set to 280 ppm for simplicity and to make sure that the climate is correctly simulated by avoiding feedbacks arising from the wrong CO 2 computed by the carbon cycle model. For the reference simulation, as the CO 2 concentration simulated by the model is close to 280 ppm it is possible to set the radiative CO 2 equal to the CO 2 computed in the carbon cycle module, but it would be important to keep them separate for other boundary conditions such as the Last Glacial Maximum as long as the computed CO 2 concentration is not equal to the data value of the period studied. Hence the two variables are considered separately in this version of the model, but they could be the same value in future studies. The cosmogenic production of 14 C is set to 2.19 atom 14 C cm −2 s −1 , which is in the pre-industrial data error bar (Masarik and Beer, 2009). This production flux is then integrated over the Earth surface and added to the 14 C concentration of the atmosphere box. The simulation starts from an equilibrium run for the climate, and uniform distribution of tracers in the ocean. The total amount of carbon has been adjusted in iterative runs to reach a value close to the pre-industrial CO 2 level of 280 ppm in the atmosphere. The simulation was run until it reached an equilibrium for deep ocean variables ( 10 000 years), and the mean of the last 100 years is used to compare the results with existing data.

Data
We compare the model results with existing observations and CMIP5 model simulations. We use temperature, salinity, phosphate and oxygen data from the World Ocean Atlas 2009 Antonov et al., 2010;Garcia et al., 2010a, b). For the DIC, alkalinity and 14 C we compare results with data from GLODAP . The pCO 2 data come from Takahashi et al. (2009) and the δ 13 C data from Schmittner et al. (2013).
The global climate models considered from the Coupled Model Intercomparison Phase 5 (CMIP5) are CESM1-BGC, CMCC-CESM, GFDL-ESM2G, GFDL-ESM2M,HadGEM2-ES, IPSL-CM5A-LR, IPSL-CM5A-MR, MPI-ESM-LR, MPI-ESM-MR and NorESM1-ME. For each variable, the models for which the data were available are listed in Table 1. For more detailed information on the models see Bopp et al. (2013). The results are averaged over the period 1890-1899 from the "historical" simulation. The end of the 19th century is chosen because it is more   Table 1. CMIP5 models considered for each variable ("×" for yes, and "−" for no).
similar to the iLOVECLIM simulation. It can be noted that very similar results for the ocean interior are obtained when considering the end of the 20th century instead, due to the long timescale of the deep ocean (a few hundred years).

Results
After equilibrium, the atmospheric CO 2 concentration is 287 ppm, the atmospheric δ 13 C value −6.4‰ and the atmospheric 14 C value 1.5 ‰, close to the pre-industrial values of respectively 279 ppm, −6.4 ‰ (Elsig et al., 2009) and 0 ‰ (Reimer et al., 2009). In the case of 14 C, the simulated −1.5 ‰ is a particularly good estimate of the observed 0 ‰, because uncertainty on pre-industrial 14 C values is of the order of 10 ‰ (Reimer et al., 2009). The ocean contains 39 019 GtC and the terrestrial biosphere 2142 GtC.
The total vegetation cover simulated by the model (Fig. 3) is in agreement with the one from another version of LOVE-CLIM (Goose et al., 2010, Fig. 14). Likewise, it is similar to the data but with an overestimation of the cover in the Tropics because of too much precipitation. In terms of carbon content, iLOVECLIM simulates low carbon contents in the regions of low vegetation cover, and particularly high carbon contents in the southern and eastern parts of North America, the northeastern part of South America, the southeastern part of Africa and on the maritime continent. This results in 2142 GtC globally, corresponding to 863 GtC for vegetation and 1279 GtC for soils (and litter). This is in the range of other model estimates which vary between around 320 and 930 GtC for vegetation and between around 500 and 3100 GtC for soils (Anav et al., 2013), as well as close to data estimates although with an overestimation of vegetation carbon content and underestimation of soil carbon content (respectively 450 to 650 GtC for vegetation, Prentice et al., 2013, and 1500 to 2400 GtC for litter and soils, Batjes, 1996).
Because the objective of this coupling is to study the climate and carbon cycle on a timescale of more than thousands of years, and because the terrestrial biosphere has already been studied (apart from the isotopes) (Goosse et al., 2010), we focus mainly on the distribution of the tracers in the ocean, both at the surface and in the interior. We also compare the carbon isotope results with data as they constitute an important constraint for past climates.

Ocean dynamics
The ocean dynamics, which depend on temperature and salinity gradients, play an important role for the carbon cycle because they partly determine the distribution of the tracers that are transported. The iLOVECLIM model simulates relatively well the distribution of temperature and salinity both at the surface and in the ocean interior.
At the surface, the simulated temperature field is similar to the observations (Fig. 4), with higher temperatures at the low latitudes and lower at high latitudes. Some local discrepancies can be observed in the boundary currents which are not well represented in the model due to its low resolution. The salinity distribution is in agreement with the data in most places (Fig. 5), except in the two bands of higher salinity in the Pacific and Indian oceans around 30 • N and 30 • S and in the northwestern part of the Indian Ocean where the simulated salinity is too low compared to observations.
In the ocean interior, the major oceanic water masses display similar features as in the data (Figs. 6 and 7). The thermocline is well represented in both the Atlantic and Pacific oceans. The Antarctic Bottom Water (AABW), which forms around Antarctica and sinks to the bottom of the ocean, is characterized by very cold temperature and low salinity in the model as in the observations. The North Atlantic Deep Water (NADW) which forms in the North Atlantic high latitudes, has relatively warmer and saltier water, in agreement with data. The low-salinity tongue of the Antarctic Intermediate Water (AAIW), which spreads northward at interme-   Ocean, AABW is too cold, so that most of the bottom ocean is slightly too cold compared to the data. In the North Atlantic the water that sinks with NADW is too salty because the surface water is also slightly too salty (Fig. 5).
The simulated streamfunction (Fig. 8) is in the range of other models, with a maximum Atlantic Meridional Overturning Circulation (AMOC) value of 21 Sv, compared to values between 14 and 31 Sv for CMIP5 models (Weaver et al., 2012). Comparing to observation of the AMOC strength (e.g. Srokosz et al., 2012, and references therein), we find an upper limb transport at 26 • N of about 15 Sv, lower that the 17 to 22 Sv estimates (Kanzow et al., 2010;Srokosz et al., 2012) from direct measurements. At 16 • N, we obtain a lower limb of about 19 Sv, in good agreement with observations (Send et al., 2011;Srokosz et al., 2012) that infer a transport of 17 ± 3.5 Sv.

Nutrients and oxygen
The distribution of nutrients depends on the transport by the diffusion-advection scheme of the ocean model, their use by marine biota (net productivity) and remineralization at depth. In the euphotic zone in the first 100 m below the surface, nutrients are consumed by phytoplankton during photosynthesis, while oxygen is produced. There are thus less nutrients at the surface than in the deep ocean, which can be seen in simulated phosphate concentrations, in agreement with data (Figs. 9 and 10). The surface distribution of simulated phosphates tends to lead to an underestimate of the intensity of boundary currents and upwellings as already seen in the surface temperature field; nonetheless, the low-tohigh latitudes gradient observed in data is well represented (Fig. 9). At the surface the oxygen is set to the saturation level (Figs. 11 and 12). The simulated surface distribution of oxygen tends to be underestimated in the Northwest Atlantic and in the Benguela upwelling, as well as in parts of the Southern Ocean (Fig. 11) but this is due to the too warm temperatures in these areas compared to data (Fig. 4), which decreases the solubility of atmospheric oxygen in the surface water. In the North Atlantic, this error then propagates in the interior resulting in too low oxygen values in the deep North Atlantic. In the ocean interior, the remineralization of plankton consumes oxygen and releases nutrients. This explains the minimum of oxygen and maximum of nutrients around 500-1000 m which is relatively well represented in the model compared to data (Figs. 10 and 12). The differences between the Atlantic and Pacific basins are also well represented. In the North Atlantic, the NADW sinks with lower phosphate values ( Fig. 10a and b) and higher O 2 values ( Fig. 12a and b) from the surface where the waters are enriched in O 2 and where nutrients are consumed for photosynthesis. The O 2 values in the ocean interior where NADW penetrates are slightly too small in the model because the surface values are too low. In the Pacific, the water is progressively enriched in PO 4 ( Fig. 10c and d) while it becomes depleted in O 2 (Fig. 12c and d) during its transport from the south to the north, because of the constant remineralization which enriches the water in PO 4 and uses O 2 .

Carbon
The simulated distribution of DIC and alkalinity is in relative agreement with the data in the oceans. At the surface, DIC is higher at high latitudes and lower at low latitudes like in the data (Fig. 13), although the DIC levels in the Tropics are slightly too low compared to the data. The alkalinity values are similar to the data, but with some small differences especially in the Atlantic where the data display two zones of higher values in the middle of the tropical gyres which are not very well represented by the model (Fig. 14). This could be due to the dissolution profile of CaCO 3 which is a function of depth, as for POC but with different values, and could be improved to be more realistic.
In the ocean interior, NADW is characterized by relatively low DIC values in the model as in the data, although the model values are slightly too high (Fig. 14). In the Pacific, the water becomes progressively enriched in DIC and alkalinity as it goes from the south to the north because of remineralization (Figs. 15 and 16). This is well represented in the model for DIC, however the alkalinity distribution is less well represented in the model, which could be due to the simple linear relation between the production of CaCO 3 and the production of organic matter, or the fixed vertical profile of remineralization.
The regions of high and low pCO 2 are generally well represented in the model compared to the data (Fig. 17). In particular, the pCO 2 values are higher around the equator, where the upwelling brings water with a high carbon content that is lost to the atmosphere, even if the model underestimates these high values. At high latitudes, especially in the North Atlantic and Arctic regions, the pCO 2 values are low where the ocean takes up carbon from the atmosphere. However, in the Southern Ocean the data indicate low values, even if they are sparse, which are not shown by the model, but the cause of this mismatch is unknown.

Carbon isotopes
During photosynthesis, the organisms preferentially use the relatively light 12 C over 13 C. This leads to higher δ 13 C values in the surface and lower values deeper in the ocean where remineralization takes place and 12 C is released. This is well represented in the model (Fig. 18), as well as the minimum value in the subsurface equatorial Atlantic due to higher remineralization in that region. The δ 13 C also depends on circulation, so that NADW is characterized by relatively high values and AABW by lower values, in agreement with data. In the Pacific, the water is progressively enriched in 12 C from remineralization from south to north, resulting in the low δ 13 C values. However, the high δ 13 C values in the North Atlantic do not penetrate far enough south, which could be due to too much diffusion.
As opposed to simulated δ 13 C, simulated 14 C does not depend on biology effects, so it allows us to separate the biological and circulation effects registered by δ 13 C. The general structure of oceanic 14 C is well simulated by the model (Fig. 19) and reflects the penetration of water masses in the interior of the ocean: from north to south in the Atlantic Ocean and from south to north in the Pacific Ocean. The model performs well compared to other ocean GCMs (Mariotti et al., 2013;Tschumi et al., 2011;Franke et al., 2008;Matsumoto et al., 2004), especially in the intermediate to deep Pacific Ocean. The model values seem nonetheless to decrease too rapidly following the penetration of NADW in the North Atlantic, similarly to PO 4 , which could indicate that the diffusion is too strong in that region. In the Pacific, the water becomes too depleted in 14 C in the northern part, possibly due to an underestimate of the mixing in that region.

Discussion
Because the main feature added to iLOVECLIM for the carbon cycle concerns the ocean, we only discuss the results for the oceanic variables. The terrestrial biosphere has only been slightly modified to include the carbon reservoirs, but could benefit from further improvements such as more plant functional types, as well as additional modules such as permafrost, which is work in progress (Kitover et al., 2013).

Model-data comparison
The iLOVECLIM model simulates most of the variables in agreement with data, especially the main characteristics of the water masses. However, a number of discrepancies exist. Some are due to errors in the simulation of surface regional features which then propagate in the ocean interior, such as the North Atlantic where the high salinity from the Tropics is transported too much northward compared to the data. This could be partly due to the resolution of the model which limits the representation of small-scale features. The misrepresentation of temperature has a direct impact on oxygen, for example again in the North Atlantic where the temperatures are too high, which leads to too small values of oxygen in the surface and in the ocean interior. Another source of error   could come from the diffusion which seems too strong in the North Atlantic where the characteristic values of NADW for salinity, PO 4 , DIC and carbon isotopes decrease too rapidly while it penetrates southward. This highlights the crucial role of a correct representation of temperature and salinity and the associated ocean circulation in setting the distribution of the biogeochemical variables. The distribution of the variables strongly depends on salinity and temperature distribution: if it is improved it should also improve the carbon cycle.

Inter-model comparison
We compare the iLOVECLIM results with other models using the data from the Coupled Model Intercomparison Phase 5 (CMIP5). We focus on three key variables (dissolved inorganic carbon, alkalinity and oxygen) using the average over years 1890-1899 of the "historical" simulation (see Sect. 2.4). The data are zonally averaged for the Atlantic and Pacific basins (including the Southern Ocean). Note that the simulations that are compared are not exactly the same: the iLOVECLIM simulation is a long simulation of a few thousand years under pre-industrial conditions, whereas the CMIP5 simulations are run under evolving boundary conditions of the historical period since 1850 starting from spin-up simulations of a few hundred to one thousand years. Additionally, the spatial resolution is higher in the CMIP5 models which are fully coupled GCMs. Nevertheless, we show here that the skill scores of iLOVECLIM are similar to those of more complex Earth system models used in CMIP5. For most variables, iLOVECLIM is in the range of other models performance. For DIC the models that statistically perform best in both the Atlantic and Pacific are the IPSL-CM5A-LR and IPSL-CM5A-MR models (Fig. 20). iLOVECLIM is less accurate than the IPSL models, but still reproduces most of the pattern and gives better results than other models such as NorESM1-ME, CMCC-CESM, GFDL-ESM2G or MPI-ESM-LR in terms of correlation and root mean square error. For alkalinity, most models simulate poorly the distribution especially in the Atlantic basin, where iLOVECLIM is performing particularly poorly (Fig. 21). In the Pacific, which represents a larger volume, the models yields better results and so does iLOVECLIM, which lies in the middle of the ensemble. This highlights the need of better understanding the processes responsible for the change of alkalinity to improve its distribution in models. For the oxygen, iLOVECLIM lies behind most models in the Atlantic but is in the middle of the range in the Pacific (Fig. 22). In the Atlantic basin, this is partly due to the representation of the high O 2 values penetrating in the North Atlantic with NADW that is not well reproduced in iLOVECLIM because the O 2 values are too low at the surface. Future work will focus on understanding the causes of the mismatch to improve the O 2 distribution. In the Pacific basin iLOVECLIM has a good correlation at around 0.8 like most models. This is not as good as a few models with correlations higher than 0.9 such as CESM1-BGC, MPI-ESM-MR and MPI-ESM-LR, but relatively good and better than NorESM1-ME with a correlation of only 0.5.

Future developments
Overall, iLOVECLIM does a relatively good job compared to the data and other models and usually lies in the middle of the CMIP5 range. This is a good performance given that iLOVECLIM is an EMIC and has a less complex and comprehensive representation of the different processes than the CMIP5 GCMs. The GCMs usually simulate better the ocean circulation which yields better distribution of the geochemical variables. There are however a few points that need to be improved in iLOVECLIM, namely the O 2 representation in the Atlantic and the alkalinity distribution (like in all other models). Some limitations arise from the simplicity of the NPZD model which does not include iron nor silicate. This could be added in future work. The air-sea flux of oxygen has not yet been parametrized depending on the difference between the atmosphere and surface water values and the wind, but this will be explored in future studies. It could improve the regional distribution of oxygen values, and would also modify the temporal evolution of oxygen values in transient simulations. Work has been done in other models showing the importance of remineralization on the carbon cycle (Schneider et al., 2008;Kwon et al., 2009). The profile, which depends on depth, is currently fixed, but the effect of changing the values depending on the temperature or other variables should be evaluated. The production and dissolution of CaCO 3 could also be improved, which would yield better results for the alkalinity distribution. In particular, CaCO 3 production is currently proportional to the production of organic matter, which could be modified, and the vertical dissolution profile is fixed, which could be changed to take into account the saturation state.

Conclusions
We have described the implementation of a carbon cycle module in the iLOVECLIM model, including the carbon isotopes 13 C and 14 C. Comparison with modern data show that the model performs well for the main carbon cycle variables, and reproduces the most important features of the different water masses. In particular, the good representation of the 13 C and 14 C in the ocean interior paves the way for past studies for which they represent most of the available data. Therefore the iLOVECLIM model with the carbon cycle is well suited for long-term simulations of a few thousand years in the past but also in the future. Some improvements will be considered in future work, such as the inclusion of iron and silicate, a better parametrization of the O 2 air-sea exchange with wind and better parametrization of the remineralization and dissolution profiles. Finally, a sediment model remains to be coupled to include all relevant oceanic components of the carbon cycle on timescales of a few thousand years.

Code availability
The iLOVECLIM source code is based on the LOVE-CLIM model version 1.2 whose code is accessible at http:// www.elic.ucl.ac.be/modx/elic/index.php?id=289. The developments on the iLOVECLIM source code are hosted at https: //forge.ipsl.jussieu.fr/ludus but are not publicly available due to copyright restrictions. Access can be granted on demand by request to D. M. Roche (didier.roche@lsce.ipsl.fr).