The Earth system model CLIMBER-X v1.0 – Part 2: The global carbon cycle

. The carbon cycle component of the newly developed Earth System Model of intermediate complexity CLIMBER-X is presented. The model represents the cycling of carbon through atmosphere, vegetation, soils, seawater and marine sediments. Exchanges of carbon with geological reservoirs occur through sediment burial, rock weathering and volcanic degassing. The state-of-the-art HAMOCC6 model is employed to simulate ocean biogeochemistry and marine sediments processes. The land model PALADYN simulates the processes related to vegetation and soil carbon dynamics, including permafrost and peatlands. 5 The dust cycle in the model allows for an interactive determination of the input of the micro-nutrient iron into the ocean. A rock weathering scheme is implemented into the model, with the weathering rate depending on lithology, runoff and soil temperature. CLIMBER-X includes a simple representation of the methane cycle, with explicitly modelled natural emissions from land and the assumption of a constant residence time of CH 4 in the atmosphere. Carbon isotopes 13 C and 14 C are tracked through all model compartments and provide a useful diagnostic for model-data comparison.


Introduction
Atmospheric CO 2 exerts a profound control on the state of the Earth system. Although it is present only in tiny concentrations 20 in the present-day atmosphere, by absorbing radiation in the longwave spectral range it has a substantial effect on the energy balance of the Earth. In the present day atmosphere, CO 2 is the second most important greenhouse gas, after water vapor. CO 2 is also a fundamental molecule for life on Earth, as it serves as 'food' in the photosynthesis process. The atmospheric CO 2 concentration is hence a main control on the growth rate of plants on land.
From ice core data it is well known that atmospheric CO 2 concentrations showed pronounced variations over the last million 25 years (e.g. Petit et al., 1999;Augustin et al., 2004) that played an important role for the climate evolution over the Pleistocene (last ∼2.6 million years) by amplifiying the variations associated with glacial-interglacial cycles (e.g. Ganopolski and Calov, 2011;Abe-Ouchi et al., 2013). Furthermore, on even longer time scales, a secular decrease in CO 2 is thought to have been the main driver of the gradual cooling over the Cenozoic (last 66 million years) (e.g. Raymo M.E. and Ruddiman W.F., 1992).
Over the last few centuries, human activities have strongly disrupted the natural CO 2 balance, by directly emitting CO 2 30 from fossil sources into the atmosphere. The resulting increase in atmospheric CO 2 has been the main factor for the observed rapid climate warming since the preindustrial period (e.g. Gulev et al., 2021).
Modelling the atmospheric CO 2 concentration is thus fundamental both for understanding past climate changes and for predicting the future evolution of the Earth system under different anthropogenic emission scenarios. However, it is far from trivial, because atmospheric CO 2 is the result of complex biogeochemical processes on land, in the ocean, in marine sediments 35 and in the lithosphere. Additionally, because of the long time scales involved in some of the carbon cycle processes, the interactive simulation of atmospheric CO 2 has been, and still is, a challenge for state-of-the-art Earth system models. Fast Earth system models of intermediate complexity have therefore been extensively employed for investigating carbon cycleclimate feedbacks (e.g. Bern3D (Müller et al., 2008;Tschumi et al., 2011;Stocker et al., 2013), cGENIE (Ridgwell et al., 2007;Cao et al., 2009), CLIMBER-2 (Brovkin et al., 2002(Brovkin et al., , 2007(Brovkin et al., , 2012, iLOVECLIM (Bouttes et al., 2015), LOVECLIM 40 (Goosse et al., 2010) and Uvic Zickfeld et al., 2011;Mengis et al., 2020)). Among these, CLIMBER-2 has successfully reproduced glacial-interglacial variations in CO 2 (Ganopolski and Brovkin, 2017;Willeit et al., 2019), but some of the processes involved remain uncertain. CLIMBER-X builds on the past experience in modelling the global carbon cycle with CLIMBER-2, but adds an improved and more detailed representation of carbon cycle processes both on land and in the ocean. Improvements include a generally higher spatial resolution, a 3D ocean model, a state-of-the-art ocean biogeochemistry 45 and marine sediment model, a more comprehensive description of vegetation and soil carbon processes, including permafrost and peatlands, and a new chemical weathering scheme.
In the following, the biogeochemistry components of CLIMBER-X are presented. The climate core of CLIMBER-X is described in detail in Willeit et al. (2022). = 4×10 −6 :1).
The marine sediment module, which is part of HAMOCC, is based on Heinze et al. (1999). It essentially simulates the same processes between dissolved tracers (DIC, TA, PO 4 , NO 3 , O 2 , Fe, SiO 2 , H 2 S and N 2 ) in pore water and solid sediment constituents (POC, opal, CaCO 3 and dust) as in the water column. Pore water tracers are exchanged with the overlying water 85 column via diffusion. Sedimentation fluxes of POC, CaCO 3 , opal and dust are added to the solid components of the sediment.
Accumulation of solid sediment material will lead to active sediment layer content being shifted to the burial layer and back if boundary conditions changes lead to chemical erosion of previously buried sediment.
Next we describe the changes introduced into HAMOCC as part of its implementation into CLIMBER-X.
N 2 fixation is represented by a diagnostic formulation, whereby the nitrate influx into the surface layer is a function of 90 the nitrate deficit relative to phosphate, multiplied by a constant fixation rate . Prognostic N 2 fixers have recently been included in HAMOCC (Paulsen et al., 2017), based on the physiological characteristics of the cyanobacterium Trichodesmium. However, for simplicity and because uncertainties in nitrogen fixation remain large (e.g. Zehr and Capone, 2020), in CLIMBER-X cyanobacteria are disabled by default.
Following Heinemann et al. (2019), we have implemented a representation of aggregates in the model. Particulate organic 95 carbon is assumed to form aggregates with the denser calcite and opal built during phytoplankton and zooplankton growth, as well as with dust particles. The sinking speed of these aggregates depends on their excess density (Gehlen et al., 2006;Heinemann et al., 2019). Note that this approach neglects the effects of, e.g., aggregate size distribution and porosity on the sinking speed (Maerz et al., 2020), and it does not, like other numerically more expensive schemes (e.g., Kriest and Evans, 2000) explicitly resolve the biological and physical aggregation and disaggregation processes. Following recent evidence that the remineralisation of organic carbon depends on temperature (e.g. Laufkötter et al., 2017), we have introduced a Q10 temperature dependence for the remineralisation of POC and DOC (Segschneider and Bendtsen, 2013;Crichton et al., 2021), with a default Q10 value of 2.
In the original HAMOCC, iron complexation by organic substances is assumed when the iron concentration exceeds a given threshold and dissolved iron is then removed from the water column at a fixed rate. In CLIMBER-X, we explicitly model iron 105 complexation differentiating between free and complexed iron forms following Archer and Johnson (2000) and Parekh et al. (2004). The complexed iron is associated with an organic ligand and only the free iron is available for scavenging. The ligand concentration is assumed to be constant at 1 nmolkg −1 with a ligand stability constant of 1×10 11 kgmol −1 . The speciation of iron is then determined by equilibrium kinetics. The scavenging rate of free iron is a combination of a minimum scavenging rate and a scavenging rate that is proportional to the POC, calcite and opal concentrations following Aumont et al. (2015) and 110 Hauck et al. (2013). Compared to HAMOCC we have also increased the stochiometric iron ratio in organic compounds from Fe:C = 3×10 −6 :1 to Fe:C = 4×10 −6 :1.
The carbon 13 isotope has been recently implemented in HAMOCC by Liu et al. (2021). In CLIMBER-X we extended this approach to also include radiocarbon.
Since the ocean model in CLIMBER-X is a rigid lid model, following the OMIP protocol (Orr et al., 2017), we explicitly take 115 into account the local concentration-dilution effect of the net surface freshwater flux, which changes surface DIC concentration and alkalinity.
Based on scale analysis, we have excluded fast sinking tracers (CaCO 3 , opal, POC and dust) from advection, as these particles have sinking speeds which are large enough so that vertical transfer between different grid cells is more rapid than horizontal transfer by advection would be, considering the relatively coarse resolution of the ocean model. Following a similar 120 line of thought also short-lived tracers like phytoplankton and zooplankton are excluded from oceanic transport. However, convection and wind-driven surface vertical mixing are applied to all biogeochemical tracers.
In CLIMBER-X, HAMOCC is integrated with a time step of one day, which is also the time step of the physical ocean model.

125
PALADYN is a comprehensive land surface-vegetation-carbon cycle model designed specifically for the use in CLIMBER-X (Willeit and Ganopolski, 2016). It includes a detailed representation of the land carbon cycle. Photosynthesis is computed following the Farquhar model (Farquhar et al., 1980;Collatz et al., 1991). Carbon assimilation by vegetation is coupled to the transpiration of water through stomatal conductance. The model includes a dynamic vegetation module with 5 plant functional types (PFTs) competing for the gridcell share based on their respective net primary productivity. The model distinguishes 130 between mineral soil carbon, peat carbon, buried carbon and shelf carbon. Each soil carbon 'type' has its own soil carbon pools generally represented by a litter, a fast and a slow carbon pool in each of the five soil layers. Carbon can be redistributed between the layers by vertical diffusion. For the vegetated macro surface type, decomposition is a function of soil temperature and soil moisture. Carbon in permanently frozen layers is assigned a long turnover time which effectively locks carbon in 5 https://doi.org/10.5194/gmd-2022-307 Preprint. Discussion started: 6 January 2023 c Author(s) 2023. CC BY 4.0 License.
permafrost. Carbon buried below ice sheets and on ocean shelfs is treated separately. The land model also includes a dynamic 135 peat module. PALADYN includes carbon isotopes 13 C and 14 C, which are tracked through all carbon pools in vegetation and soil. Isotopic discrimination is modelled only during the photosynthetic process. A simple methane module is implemented to represent methane emissions from anaerobic carbon decomposition in wetlands and peatlands. The integration of PALADYN into the coupled CLIMBER-X framework and subsequent sensitivity analyses of the land carbon cycle feedbacks, which were not performed with the offline PALADYN setup in Willeit and Ganopolski (2016), highlighted the need to improve certain 140 aspects of the model. These improvements are described next.
The rubisco-limited photosynthesis rate the version of PALADYN model described in Willeit and Ganopolski (2016) was based on the 'strong optimality' hypothesis of Haxeltine and Prentice (1996), which assumes that rubisco activity and the nitrogen content of leaves vary with canopy position and seasonally so as to maximize net assimilation at the leaf level (Schaphoff et al., 2018). However, we found that this formulation led to a relatively small increase in gross primary production over the 145 historical period, which resulted into an overestimation of atmospheric CO 2 in coupled historical simulations. We therefore introduced a new formulation for the maximum rubisco capacity, with dependencies on PFT-specific, constant foliage nitrogen concentration, specific leaf area and leaf temperature following Thornton and Zimmermann (2007) as implemented in CLM4.5 (Oleson et al., 2010).
In the original PALADYN formulation, the internal leaf CO 2 concentration used for photosynthesis was computed based 150 on the Cowan-Farquhar optimality hypothesis (Medlyn et al., 2011). In the new model version, for C3 plants, we have implemented an alternative scheme following the more general least-cost optimality model (Prentice et al., 2014;Lavergne et al., 2019) with the moisture dependence proposed by Lavergne et al. (2020).
In the isotopic discrimination during photosynthesis (∆) we included an explicit fractionation term for photorespiration as recommended by several recent studies (Ubierna and Farquhar, 2014;Schubert and Jahren, 2018;Lavergne et al., 2019): where c a and c i are the ambient and leaf internal CO 2 concentrations, p a is the ambient partial pressure of CO 2 and Γ is the CO 2 compensation point.
In the dynamic vegetation model a parameter (λ) is used to partition the net primary production (NPP) between local growth of existing vegetation and lateral expansion ('spreading') of vegetation coverage within the grid cell, with all of the NPP being 160 used for growth for small leaf area index (LAI) values, and all the NPP being used for 'spreading' for large LAI values. λ is assumed to be a piecewise linear function of the leaf area index between a minimum and maxium LAI. For small leaf area indices, all of the NPP is used for local growth (λ = 0); for LAI above a critical value LAI min , a fraction (λ > 0) is used for 'spreading': 165 However, since the simulated leaf area index depends strongly on NPP, which in turn has a pronounced dependence on atmospheric CO 2 , this formulation results in a strong dependence of λ on CO 2 , with an increasingly larger fraction of NPP being 6 https://doi.org/10.5194/gmd-2022-307 Preprint. Discussion started: 6 January 2023 c Author(s) 2023. CC BY 4.0 License.
used for 'spreading' as CO 2 increases. We have therefore implemented a CO 2 dependence in the maximum leaf area index to reduce this effect: (3)

170
The fraction of decomposed litter respired directly as CO 2 to the atmosphere has been reduced from 0.7 to 0.6 and the fraction of decomposed litter transferred to the slow soil carbon pool has been doubled from 0.015 to 0.03. Together these changes result in more carbon being accumulated into the soil.
A simple representation of land use change has been introduced into the model following Burton et al. (2019) as described in Willeit et al. (2022). A fraction of each grid cell is prescribed as being used for agriculture and land use is then represented as a 175 limitation to the space available for the woody PFTs to expand into. When forests and shrubs are affected by land use change, an additional disturbance rate of 1 yr −1 is prescribed on top of the standard background disturbance, leading to vegetation dying. The resulting dead vegetation carbon is then added as litter to the soil carbon pools, and a large part will be respired directly to the atmosphere within a few years. Soil carbon is assumed to not be directly affected by land use practices.
The partitioning of the soil carbon decomposed under anaerobic conditions into CO 2 and CH 4 was a prescribed constant 180 in Willeit and Ganopolski (2016). We modified this by making the fraction released as CH 4 dependent on temperature with a Q10 of 1.8, following Riley et al. (2011) and Kleinen et al. (2020).
We implemented a chemical weathering model to compute the riverine fluxes of bicarbonate ions (HCO − 3 ), (and therefore dissolved inorganic carbon and alkalinity) to the ocean and the consumption of atmospheric CO 2 . The weathering rate depends on the lithology and on the climate variables temperature and runoff. The lithological map of Hartmann and Moosdorf (2012) 185 distinguishing 16 different lithologies is used to describe the spatial distribution of rocks. The parameters for the chemical weathering equations for all lithologies, except for carbonate sedimentary rocks and loess, are based on a spatially explicit runoff-dependent model of chemical weathering, which was calibrated for 381 catchments in Japan (Hartmann, 2009), with the additional temperature dependence of Hartmann et al. (2014). The effect of soil shielding on the weathering rate suggested by Hartmann et al. (2014) has not been considered since information on soil shielding is not readily available for periods 190 beyond the recent past. For carbonate sedimentary rocks, the weathering rate follows the approach of Amiotte Suchet and Probst (1995) with a dependence on runoff. Alternatively, the temperature dependent formulation of Romero-Mujalli et al.
(2019) is available for use in the model. The weathering rate for loess sediments depends on runoff following Börker et al. (2020). The global distribution of loess cover for present-day and for the last glacial maximum, as well as the lithologies of the continental shelves that were exposed at the last glacial maximum, is taken from Börker et al. (2020).

195
The carbon isotopes fluxes from chemical weathering are computed assuming a δ 13 C of 1.8 permil for carbon originating from carbonate minerals (Derry and France-Lanord, 1996).
Equations describing silicate and phosphorus weathering fluxes are also available as part of the weathering model. However, silicate and phosphorus riverine fluxes are not considered in the default model setup, as they would result in further complications related to the conservation of nutrients in the ocean. Instead, as discussed in section. 2.1, the silicate and phosphorus 200 budgets are closed by assuming that the sediment burial flux is prescribed as input at the ocean surface.

Atmospheric CO 2
The atmospheric CO 2 concentration in CLIMBER-X is a globally uniform value. It can either be prescribed (as constant or time-dependent) or interactively computed by the model from the following prognostic equation for the total carbon content stored as CO 2 in the atmosphere (C atm ): The source and sink terms on the right hand side represent, from left to right, the net sea-air carbon flux, the global net land to atmosphere carbon flux, the anthropogenic carbon emissions (excluding land-use change), the CO 2 consumption by silicate and carbonate weathering, the volcanic degassing flux and the CO 2 flux from the oxidation of atmospheric methane originating from non-agricultural sources. The CO 2 consumption by weathering is computed assuming that all carbon in the The constant volcanic degassing rate is set to half the silicate weathering rate (e.g. Munhoven and François, 1994) as determined by an equilibrium spinup simulation: The flux from the oxidation of methane, F CH4ox , is computed by the CH 4 model as described in Sec. 2.4 below. The atmospheric CO 2 concentration is then computed from C atm using a conversion factor of 2.12 PgCppm −1 (Denman et al., 2007).
Equations similar to eq. 4 are used also for the carbon isotopes 13 C and 14 C. The prognostic equation for the stable isotope 13 C in atmospheric CO 2 is: The 13 C fluxes from land and ocean are explicitly computed by the land and ocean carbon cycle models as described in detail in Willeit and Ganopolski (2016) and Liu et al. (2021). The δ 13 C of anthropogenic carbon emissions is taken to be -26 permil and the 13 C flux from CO 2 consumption by weathering, assuming no fractionation, is simply computed as: The 13 C of volcanic degassing is computed assuming a δ 13 C of -5 permil.
The prognostic equation for radiocarbon 14 C in atmospheric CO 2 reads: Carbon sources originating from geological reservoirs, i.e. volcanic degassing, are assumed to contain no radiocarbon. Similarly, radiocarbon is assumed to be absent in anthropogenic carbon emissions from fossil fuel burning, because the age of fossils far exceeds the half-life of 14 C. The production rate of radiocarbon in the atmosphere (F 14 prod ) is prescribed in the model and the radiocarbon decay time is τ14 C = 8267 yr.

Atmospheric CH 4
Similarly to CO 2 , atmospheric CH 4 is also considered to be well-mixed in the atmosphere and is therefore represented as a globally uniform value. The atmospheric CH 4 concentration can be prescribed, or it can be interactively computed by the 235 model from: Methane sources include natural emissions from wetlands and peatlands (F emis lnd ), which are explicitly simulated by the model as originating from anaerobic decomposition processes of carbon in soils (Willeit and Ganopolski, 2016). Other natural sources of methane are generally smaller (e.g. Saunois et al., 2020;Kleinen et al., 2020) and are neglected here for simplicity. Anthro-240 pogenic methane emissions (F emis anth ) are prescribed in the model. The sink of methane from oxidation in the atmosphere is computed using a constant residence time of CH 4 , τ CH4 = 9.5 years, which is a reasonable first approximation at least for climate conditions ranging between the last glacial maximum and present-day (Kleinen et al., 2020;Levine et al., 2011;Hopcroft et al., 2017). The first (and simplest) setup consists of ocean, land and atmosphere carbon cycle components only. In this setup marine sediments are disabled and particulate fluxes that reach the ocean floor are completely remineralised/dissolved in the bottom ocean grid cell. Rock weathering from land is also switched off, so that the carbon exchange between ocean, land and atmosphere occurs only through air-sea fluxes and through land-atmosphere exchanges. In this setup the carbon system is closed, 250 in the sense that there are no natural sources and sinks from and to geological reservoirs. As a response to an external climate perturbation, carbon is then simply redistributed between atmosphere, ocean and land, with the total carbon in the system being conserved. This setup is equivalent to what is used in many state-of-the-art Earth System models for climate change projections on centennial time scales (e.g. Séférian et al., 2020). The model spinup for this simple setup is straightforward and requires only to run the model to steady state with a prescribed atmospheric CO 2 concentration for ≈ 10,000 years. The slowest time 255 scale in this setup is given by the slow decomposition rate of organic carbon in frozen soils, which is limited to a maximum value set by default to 5000 years. The initial state for the spinup run is given by observed present-day 3D concentrations of different tracers in the ocean (Lauvset et al., 2016;Olsen et al., 2016;Garcia et al., 2013b), while the land surface is assumed to be covered by bare soil and with no carbon stored on land.
The closed carbon cycle setup is applicable to simulations of up to 1000 yr. On longer time scales, sediment and weathering 260 processes become important and need to be accounted for when performing long-term transient simulations with interactive 9 https://doi.org/10.5194/gmd-2022-307 Preprint. Discussion started: 6 January 2023 c Author(s) 2023. CC BY 4.0 License.
Although it is unlikely that in reality the slow carbon cycle processes related to marine sediments, peatlands and permafrost carbon are in equilibrium at any specific point in time, for practical reasons we assume that such an equilibrium is a reasonable first approximation. Assuming that the preindustrial is an equilibrium state of the climate-carbon cycle system allows to run perturbation experiments with interactive carbon cycle without having to deal with possible long-term drifts in 265 atmospheric CO 2 . However, the long time scale of ∼100,000 years involved in ocean sediment processes represents a challenge in running the model into equilibrium, even for a high-throughput model like CLIMBER-X. We therefore implemented a scheme to run the physical ocean and ocean biogeochemistry models in an offline setup with prescribed climatological daily input fields at the ocean surface. This setup results in a speedup of a factor >2 relative to running the fully coupled climatecarbon cycle model, meaning that ocean carbon cycle and marine sediments can be run into equilibrium in about a week of 270 computing time on a high performance computer. In detail, the spinup procedure of the full carbon cycle configuration comprises two different stages. Atmospheric CO 2 is prescribed to a constant value throughout the process, at 280 ppm for the pre-industrial case. The first spinup stage aims at spinning up the sediment model. For that the full carbon-cycle climate model is run for 5,000 years and every 300 years the sediment model is run in stand-alone for 1000 years. During this stage all net fluxes into the sediments are compensated for and prescribed as input at the ocean surface in order to approximately conserve 275 water column tracer inventories while the sediments are filling up. In the second stage we switch to simulated DIC and alkalinity weathering fluxes from land and at the same time also switch to the more efficient offline ocean-biogeochemistry setup described above and run the model until an approximate equilibrium is reached after ∼100,000 years (Fig. 2). A simplification that is made in the open carbon cycle setup is that organic carbon and opal that are buried into the sediments, and are therefore effectively leaving the system, are returned in remineralized form to the surface ocean so that phosphorus and silica inventories 280 of the ocean-sediment system are conserved throughout the simulation.
The carbon fluxes among the different model components in the open setup for equilibrium pre-industrial conditions are schematically illustrated in Fig. 3. The volcanic degassing rate is equal to half the atmospheric CO 2 consumption by silicate weathering, in accordance with theory (Munhoven and François, 1994). Note that not only the carbon budget of the different compartments (atmosphere, ocean, lithosphere) is well balanced, but also the ocean alkalinity budget is.   give an overview of how CLIMBER-X compares to state-of-the-art Earth system models based on general circulation models, we also include results from model simulations from the recent coupled model intercomparison project CMIP6 (Eyring et al., 2016). In particular, for ocean biogeochemistry, we highlight how the model compares with results from the MPI-ESM1-2-LR 300 employing the original marine carbon cycle model HAMOCC6.

Ocean biogeochemistry and marine sediments
An overview of simulated global variables characterising the ocean carbon cycle are presented and compared to observationbased estimates in Table 1, providing a summary of model performance for the present day.
The spatial pattern of air-sea CO 2 exchange is well captured by the model (Fig. 4), with outgassing generally taking place 305 in the tropics and CO 2 being taken up in mid-to high northern latitudes and in mid-latitudes of the Southern Hemisphere.
The main difference compared to other models is observed around the Equator, with a less pronounced peak in CO 2 release simulated by CLIMBER-X (Fig. 5), which is likely related to deficiencies in the simulated ocean circulation close to the equator, where the geostrophic approximation employed in CLIMBER-X reaches its limit of applicability. In the Southern Ocean, most CMIP6 models tend to overestimate the CO 2 uptake compared to observations (e.g. Gruber et al., 2009) (Fig. 5), 310 while CLIMBER-X is apparently more consistent with recent estimates, although with substantial differences in the spatial distribution of the CO 2 flux (Fig. 4). Notably, in the Southern Ocean the CLIMBER-X air-sea CO 2 exchange diverges from  Opal burial 6.3 2.7-9.9 TmolSiyr −1 Tréguer and De La Rocha (2013) that simulated by the MPI-ESM1-2-LR model (Fig. 5), which employs the original HAMOCC6 ocean biogeochemistry model. This is possibly related to the lower simulated net primary production in the Southern Ocean in CLIMBER-X compared to MPI-ESM (Fig. 6a).

315
The export of particulate organic carbon from the euphotic layer drives the biological pump and generally follows the primary productivity pattern, with modifications due to varying sinking speeds and remineralisation rates of POC in the water column. While the net primary productivity in CLIMBER-X is in line with CMIP6 models (Fig. 6a) and the globally integrated value agrees well with observations (Tab. 1), the export production in the model is generally at the lower end of the CMIP6 model range (Fig. 6b). CaCO 3 and opal export are compared to CMIP6 models in Fig. 6c,d.

320
Primary production in the ocean is limited by the availability of nutrients. Over large parts of the surface ocean nitrogen concentrations constitute the main limiting factor for photosynthesis in CLIMBER-X (Fig. 7). However, over the Southern    Ocean, in the equatorial Pacific and in the North Pacific production is limited by the availability of iron (Fig. 7). This is in accordance with observations showing that iron limitation is usually important where subsurface nutrient supply is enhanced, such as in oceanic upwelling regions (e.g. Moore et al., 2013). Since one of the main iron sources in the ocean is from mineral 325 dust deposited at the ocean surface (e.g. Tagliabue et al., 2016), iron limitation is confined to regions with low dust deposition.
The dust cycle is an integral part of CLIMBER-X, and the dust deposition is therefore explicitly modelled. The simulated dust deposition compares reasonably well with estimates from complex ESMs for the present-day (Fig. 8), although they are relatively poorly constrained. A comparison of dust deposition fluxes with observations over land further indicates that the model is able to capture the general pattern of dust deposition rate (Fig. 9).

330
The simulated dissolved iron concentration in surface water is closely related to the dust deposition shown in Fig. 8. It is therefore high in the Atlantic and Indian oceans, lower in the Southern Ocean and very small over large parts of the Pacific   North Pacific and low values elsewhere (Fig. 12, Fig. 11a). The most pronounced model biases are found in too high nitrate concentrations in the Arctic and too low values in the North Pacific. The simulated basin-wide vertical distribution of nitrate is 340 in very good agreement with observations (Fig. 13). The simulated basin-wide vertical distribution of different ocean tracers compares well with observations (Fig. 13), and generally lies within the relatively large range of values returned by CMIP6 models.
The 3D phosphate distribution in the global ocean is nicely captured by the model (Fig. 14, Fig. 11b, Fig. 13), except for too low concentrations simulated in the surface ocean of the North Pacific and the northern Indian Ocean. The negative bias in 345 the North Pacific is consistent with the too low simulated surface nitrate concentrations, both originating from a too vigorous ventilation of water masses in the upper kilometer in the physical ocean model.
As a result of reduced primary productivity in the Southern Ocean in CLIMBER-X compared to MPI-ESM1-2-LR, both surface nitrate and phosphate concentrations are consistently higher in CLIMBER-X (Fig. 11a,b), as less nutrients are assimilated during photosynthesis.

350
Silicate concentration is generally overestimated in the sub-surface ocean and underestimated in the deep North Pacific and North Indian oceans (Fig. 15). A positive surface silicate concentration bias in the Southern Ocean is seen in most CMIP6 models (Fig. 11c).
The large scale patterns of oxygen concentration in ocean waters simulated by CLIMBER-X is largely consistent with observations ( Fig. 16), but the extent and depth of the oxygen minimum zones, in particular in the Eastern equatorial Pacific, is 355 overestimated. This bias is common to many CMIP models (e.g. Cabré et al., 2015). Other biases include a too oxygen depleted Southern Ocean and too high oxygen concentrations in the upper North Pacific and North Indian oceans, again resulting from the excessive water mass ventilation in those regions as discussed above. Both DIC and alkalinity are generally overestimated in the upper ocean (Fig. 13), particularly in Antarctic intermediate water masses (Fig. 17,18). The simulated DIC concentration is too low in the upper North Pacific and Indian Oceans, while 360 alkalinity is underestimated in the deep Pacific.
Carbon isotopes in the ocean help to track the distribution of different water masses. The higher δ 13 C values in the Atlantic compared to the Pacific Ocean, originating mainly from the pronounced overturning circulation in the Atlantic, which is absent in the Pacific, are generally captured by the model (Fig. 19). The main bias is a too low δ 13 C associated with Antarctic intermediate waters and too high δ 13 C values in the upper North Pacific and Indian oceans. δ 13 C is also overestimated at 365 all depths in the Arctic. The radiocarbon ventilation age of the deep ocean is nicely reproduced by CLIMBER-X, while radiocarbon age is systematically overestimated in the upper kilometer across all ocean basins (Fig. 20). The too old (in terms of radiocarbon age) sub-surface waters could be a result of the model not explicitly resolving synoptic processes in the atmosphere, leading to an underestimation of the wind-driven mixing of passive carbon tracers in the upper ocean.
In the Atlantic and Indian Ocean CaCO 3 dominates the sediment composition, in accordance with observations (Fig. 21a,d).

Land carbon cycle
A detailed evaluation of the land carbon cycle component has already been presented in the original PALADYN description paper (Willeit and Ganopolski, 2016). However, here we partly repeat the analysis to show the model performance in the 380 coupled climate model setup and with the additional modifications to the model described above.
A selection of simulated global variables characterising the land carbon cycle is presented and compared to observationbased estimates in Table 2, providing a summary of model performance for the present day.
Photosynthesis is the basic process by which carbon enters the land domain. The simulated gross primary production (GPP), which quantifies this process, is in good agreement with observational estimates, both in terms of global integral (Tab. 2) and 385 in terms of spatial distribution (Fig. 22a,b,c).   The total carbon stored in vegetation, both above ground and below ground, is slightly overestimated in the model (Tab. 2), but the meridional distribution, mainly originating from large-scale differences in precipitation, is well reproduced (Fig. 22d,e,f).
Most of soil carbon in CLIMBER-X is stored in cold soils of the NH high-latitudes, in agreement with observations ( Fig. 22g,h,i).
However, compared to estimates from Carvalhais et al. (2014) the soil carbon distribution is too skewed towards high northern 390 latitudes and there is too little carbon in the tropics. Most CMIP6 models underestimate soil carbon in the tropics as well (Fig. 22j).
Total soil carbon estimates vary considerably, in particular when considering the full soil depth (Fan et al., 2020, e.g.). In CLIMBER-X, three quarters (∼1500 PgC) of the total ∼2200 PgC in soils are stored in the top soil meter, in good agreement with different estimates (Tab. 2). The carbon contained in areas affected by permafrost is ∼860 PgC, a bit lower than suggested 395 by observations (Tab. 2). The simulated peatland extent is lower than estimated (Yu et al., 2010), and consistently also the peat carbon is underestimated accordingly (Tab. 2).
The turnover time of terrestrial ecosystem carbon is an integrated quantitative measure of the residence time of carbon on land, from the time it is fixed by photosynthesis to the time it is returned to the atmosphere through respiration processes. It is computed as the ratio between land carbon stocks (vegetation+soil) and gross primary production. The ecosystem carbon    based estimates from Fan et al. (2020) (Fig. 22j,k,l). However, it should be noted that the large uncertainties in soil carbon content result in a rather uncertain estimated ecosystem carbon turnover time (Fan et al., 2020).
The global maximum monthly wetland extent in CLIMBER-X agrees well with observations (Tab. 2), although with substantial differences in the geographic distribution (Fig. 23). Compared to the multi-satellite product from GIEMS (Global

405
Inundation Extent from Multi-Satellites) (Prigent et al., 2007;Papa et al., 2010) the model simulates larger wetland extent in tropical forest areas. However, if compared to other wetland products based on data other than from satellite, GIEMS is underestimating wetlands below dense forests (e.g. the Amazon forest) (e.g. Melack and Hess, 2010) In south-east Asia, the GIEMS wetland extent also includes extensive rice cultivation areas, which are not represented in the model.
In CLIMBER-X methane is emitted exclusively from wetlands. However, because of the dependence of methane emissions 410 on soil carbon decomposition rates and because of the temperature dependence of the fraction of wetland carbon respired as methane, wetland methane emissions are dominated by tropical sources (Tab. 2, Fig. 24), in agreement with observations (e.g. Saunois et al., 2020). The total CH 4 emissions from wetlands are at the high end of recent estimates, which is a result of tuning the emissions in the model to match the observed emissions from all natural sources.
Chemical weathering fluxes are generally high where runoff is high, with the separation between silicate and carbonate 415 weathering being modulated by lithological properties (Fig. 25). The global CO 2 consumption rate by weathering and the alkalinity flux to the ocean in form of bicarbonate produced by rock weathering are in good agreement with observational estimates (Tab. 2), while the partitioning between carbonate and silicate weathering is skewed toward carbonate weathering (Tab. 2).

Historical period 420
As shown by Willeit et al. (2022), the historical climate evolution is well simulated by CLIMBER-X. Here we extend this analysis by focusing on the carbon cycle response.
The historical atmospheric CO 2 concentration is well reproduced by the model, with CO 2 at the year 2015 being within a few ppm of direct measurements (Fig. 26). Biases in simulated CO 2 of ∼10 ppm are quite common in state-of-the-art ESMs (e.g. Hoffman et al., 2014;Friedlingstein et al., 2014).

425
The partitioning of the anthropogenic carbon emitted over the historical period among the different spheres is compared with recent estimates of the Global Carbon Budget (GCB) (Friedlingstein et al., 2022) by the Global Carbon Project in Fig. 27. The amount of fossil carbon emitted from anthropogenic activities is prescribed from empirical data and therefore by definition matches with estimates from Friedlingstein et al. (2022). The carbon emissions resulting from land use change practices are underestimated in CLIMBER-X compared to the GCB, although the actual values remain uncertain (e.g. Gasser et al., 2020).  A substantial fraction of this anthropogenic CO 2 emission is absorbed by the ocean and the land, while the rest remains in the atmosphere. In CLIMBER-X, the ocean carbon uptake is a bit lower and the land carbon uptake a bit higher than GCB estimates, but the net effect is a realistic airborne fraction of carbon remaining in the atmosphere. The ocean carbon uptake is driven by the chemical disequilibrium between surface air CO 2 concentrations and the concentration of dissolved CO 2 in the surface ocean water and is relatively well understood, as also indicated by the narrow uncertainty range obtained from 435 different CMIP6 models (Fig. 28a). The CLIMBER-X ocean carbon uptake falls within this narrow range although it tends to be at the lower end. The land carbon uptake is largely driven by an increase in gross primary productivity as a response to increasing atmospheric CO 2 . The net primary productivity increase simulated by CLIMBER-X over the historical period is in agreement with what is shown by most CMIP6 models (Fig. 28b). However, the effect of this NPP increase on vegetation carbon varies widely among models (Fig. 28c), also because of the confounding factor of land use change. In CLIMBER-X the 440 net effect is a vegetation carbon stock decrease by ∼50 PgC. The historical evolution of soil carbon additionally depends on the response of microbial decomposition to changing environmental conditions, particularly soil temperatures. The increasing NPP, and consequently larger input of litter carbon into the soil, dominates over the negative effect of increasing temperatures in CLIMBER-X, leading to an increase in soil carbon by ∼50 PgC (Fig. 28d).    (2002) Alkalinity flux to ocean 33.9 30-40 Tmolyr −1 Amiotte Suchet et al. (2003); Gaillardet et al. (1999) Since CLIMBER-X is enabled with carbon isotopes, it also allows for a comparison of isotopic signatures to observations, 445 thereby providing additional constraints on processes involved in carbon cycle exchanges. As an example, the historical δ 13 C of atmospheric CO 2 is compared to observations in Fig. 29.
The general historical trend in atmospheric CH 4 is captured by the model (Fig. 30a). Prescribed anthropogenic methane emissions are the dominant source for the increase of the atmospheric methane burden, but natural emissions from land are also increasing due to the increase in NPP and soil temperature (Fig. 30b).

Carbon cycle feedbacks
Carbon cycle processes both on land and in the ocean are sensitive to changes in climate and atmospheric CO 2 . This implies that carbon fluxes between ocean and atmosphere and between land and atmosphere will respond to changes in climate and CO 2 concentration, which will in turn affect CO 2 and consequently climate. Quantifying the strength of these feedbacks is important to understand how the climate will respond to anthropogenic CO 2 emissions. For that, a linear feedback decomposition analysis 455 has been proposed by Friedlingstein et al. (2006) to estimate these feedbacks in Earth system models. The analysis relies on a set of model simulations under idealized 1 %yr −1 CO 2 increase experiments, whereby in one simulation the CO 2 increase is seen only by the radiative code in the atmosphere (radiatively coupled), in another one the CO 2 increase is seen only by the carbon cycle (biogeochemically coupled) and in a final one both radiation and carbon cycle see the increase in atmospheric CO 2 (fully coupled). This set of simulations allows to roughly separate the effect of changes in climate and changes in CO 2 on 460 land and ocean carbon fluxes. To a first approximation, the carbon cycle feedback to climate is usually quantified using global temperature changes, while in reality climate obviously influences the carbon cycle in more complex ways than just through (global) temperature.
The carbon cycle feedback parameters have been estimated for several phases of the coupled model intercomparison projects (Friedlingstein et al., 2006;Arora et al., 2013Arora et al., , 2020. In Figure 31 the carbon cycle-climate (γ) and the carbon cycle-465 concentration (β) feedbacks in CLIMBER-X are compared to CMIP6 model results separately for land and ocean. An increase in CO 2 has a positive effect on the uptake of carbon by both land and ocean, resulting in a lowering of CO 2 and therefore a negative feedback on climate (implying positive β, Fig.31a,c). Conversely, a climate warming has a generally negative impact on the ability of ocean and land to store carbon, leading to a positive feedback loop (implying negative γ, Fig.31b,d). The β and γ are generally in good agreement with estimates from CMIP6 models (Arora et al., 2020), although in CLIMBER-X the 470 CO 2 fertilisation effect on land is rather high (Fig.31a) and the ocean carbon uptake as a response to an increase in atmospheric CO 2 is at the lower end of the CMIP6 models (Fig.31c).

The zero emissions commitment
The zero emissions commitment (ZEC) is a measure of the amount of additional future temperature change following a complete cessation of CO 2 emissions (e.g. Matthews and Solomon, 2013). A model intercomparison project has been established 475 in an effort to analyze and compare the ZEC in different Earth system models (Jones et al., 2019). Here we use this standardized and idealized experimental setup to compare the carbon cycle response in CLIMBER-X with results from the ZECMIP models for the 1000 PgC emission pulse (MacDougall et al., 2020). The experiment branches off from a 1 % per year CO 2 increase run with CO 2 emissions set to zero at the point of 1000 PgC of total carbon emissions. We performed this experiment with both the open and closed carbon cycle setups.

480
The results of the CLIMBER-X simulations are generally well within the large range of results from state-of-the-art ESMs and EMICs participating in ZECMIP (Fig. 32). Atmospheric CO 2 concentration is rapidly decreasing after stopping the carbon emissions (Fig. 32b), while global temperature shows a more modest decrease (Fig. 32a). The ocean continues to take up carbon throughout the simulation (Fig. 32c), while the land turns from a sink to a source of carbon roughly at the time of the peak in CO 2 (Fig. 32d). CLIMBER-X initially shows a relatively weak ocean carbon uptake compared to ZECMIP models, while the 485 land carbon uptake is at the high end of the ZECMIP model ensemble.
The difference between the experiments with open and closed carbon cycle setup are negligible for the first few centuries but continue to grow over time, with CO 2 decreasing faster in the open setup (Fig. 32b).

Conclusions
We have described the major features of the carbon-cycle components of the newly developed CLIMBER-X Earth System  the coupled climate-carbon cycle evolution over timescales ranging from decades up to ∼100,000 years, while also allowing a quantification of related uncertainties.
Code and data availability.