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 the 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 sediment processes. The land model PALADYN simulates the processes related to vegetation and soil carbon dynamics, including permafrost and peatlands. 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 in 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₄ in the atmosphere. Carbon isotopes ¹³C and ¹⁴C are tracked through all model compartments and provide a useful diagnostic for model–data comparison.

A comprehensive evaluation of the model performance for the present day and the historical period shows that CLIMBER-X is capable of realistically reproducing the historical evolution of atmospheric CO₂ and CH₄ but also the spatial distribution of carbon on land and the 3D structure of biogeochemical ocean tracers. The analysis of model performance is complemented by an assessment of carbon cycle feedbacks and model sensitivities compared to state-of-the-art Coupled Model Intercomparison Project Phase 6 (CMIP6) models.

Enabling an interactive carbon cycle in CLIMBER-X results in a relatively minor slow-down of model computational performance by ∼ 20 % compared to a throughput of ∼ 10 000 simulation years per day on a single node with 16 CPUs on a high-performance computer in a climate-only model set-up. CLIMBER-X is therefore well suited to investigating the feedbacks between climate and the carbon cycle on temporal scales ranging from decades to >100 000 years.

1 Introduction

Atmospheric CO₂ exerts a profound control on the state of the Earth system. Although it is present only in tiny concentrations 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₂ is the second most important greenhouse gas after water vapour. CO₂ is also a fundamental molecule for life on Earth, as it serves as “food” in the photosynthesis process. The atmospheric CO₂ 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₂ concentrations showed pronounced variations over the last million 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 amplifying the variations associated with glacial–interglacial cycles (e.g. Ganopolski and Calov, 2011; Abe-Ouchi et al., 2013). Furthermore, on even longer timescales, a secular decrease in CO₂ is thought to have been the main driver of the gradual cooling over the Cenozoic (last 66 million years) (e.g. Raymo and Ruddiman, 1992).

Over the last few centuries, human activities have strongly disrupted the natural CO₂ balance by directly emitting CO₂ from fossil sources into the atmosphere. The resulting increase in atmospheric CO₂ has been the main factor for the observed rapid climate warming since the pre-industrial period (e.g. Gulev et al., 2021).

Modelling the atmospheric CO₂ 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₂ is the result of complex biogeochemical processes on land, in the ocean, in marine sediments and in the lithosphere. Additionally, because of the long timescales involved in some of the carbon cycle processes, the interactive simulation of atmospheric CO₂ 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 cycle–climate 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, 2007, 2012), iLOVECLIM (Bouttes et al., 2015), LOVECLIM (Goosse et al., 2010) and Uvic (Eby et al., 2009; Zickfeld et al., 2011; Mengis et al., 2020). Among these, CLIMBER-2 has successfully reproduced glacial–interglacial variations in CO₂ (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 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).

2 Model description

CLIMBER-X represents the cycling of carbon through the atmosphere, vegetation, soils, seawater and marine sediments. Through sediment burial, chemical weathering of rocks and volcanic degassing, carbon is also exchanged with geological reservoirs. A schematic illustration of the carbon cycle in the model is shown in Fig. 1. The carbon cycle component of CLIMBER-X consists of the ocean biogeochemistry and marine sediment models from HAMOCC6 (Maier-Reimer and Hasselmann, 1987; Ilyina et al., 2013; Heinze et al., 1999; Mauritsen et al., 2019) and the land model PALADYN (Willeit and Ganopolski, 2016), which includes dynamic vegetation, a soil carbon model and the weathering model of Hartmann (2009) and Börker et al. (2020). The atmospheric CO₂ concentration is determined interactively by the exchange of carbon between the atmosphere, seawater, land and lithosphere. The model includes a representation of the dust cycle, with simulated dust deposition determining the input of the micro-nutrient iron into the ocean. CLIMBER-X also 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₄ in the atmosphere. The model is enabled with the carbon isotopes ¹³C and ¹⁴C, which are tracked through all model compartments.

The different model components are described in more detail in the following sections.

Figure 1 Schematic illustration of the natural biogeochemical cycles in the CLIMBER-X model.

2.1 Ocean biogeochemistry and marine sediments: HAMOCC

HAMOCC (Maier-Reimer and Hasselmann, 1987; Maier-Reimer et al., 1993; Ilyina et al., 2013) is a state-of-the-art ocean biogeochemistry model, which is part of the MPI-ESM, the Earth system model of the Max Planck Institute for Meteorology (MPI). The latest version (Mauritsen et al., 2019), which is the version employed by the MPI in the Coupled Model Intercomparison Project Phase 6 (CMIP6), has been the starting point for the implementation of the model in CLIMBER-X. As a first step, the original HAMOCC6 code has been adapted to the CLIMBER-X structure. Notably, for easier parallelization, it has been transformed from a 3D model into a 1D vertical column model in which each water column is independent of the others. This is possible because the biogeochemical processes in the model are restricted to local vertical interactions. The different columns interact only through horizontal advection by ocean currents, which takes place in the ocean model.

HAMOCC represents the biogeochemical processes in the water column, in the sediments and at the air–sea interface. Marine biology dynamics are based on an extended NPZD (nutrients, phytoplankton, zooplankton and detritus) approach (Six and Maier-Reimer, 1996). The carbonate chemistry in the model follows the latest OMIP protocol (Orr et al., 2017), which uses the robust and safe pH calculation routines from SolveSAPHE-r1 (Munhoven, 2013). In the water column, the following biogeochemical tracers are simulated: dissolved inorganic carbon (DIC), total alkalinity (TA), phosphate (PO₄), nitrate (NO₃), nitrous oxide (N₂O), dissolved nitrogen gas (N₂), silicate (SiO₂), dissolved bioavailable iron (Fe), dissolved oxygen (O₂), phytoplankton (Phy), zooplankton (Zoo), dissolved organic matter (DOC), particulate organic matter (POC), opal shells, calcium carbonate shells (CaCO₃), terrigenous material (dust) and hydrogen sulfide (H₂S). The composition of organic material follows a constant Redfield ratio ($\mathrm{C}:\mathrm{N}:\mathrm{P}:{\mathrm{O}}_{\mathrm{2}}\mathrm{\hspace{0.17em}}=\mathrm{\hspace{0.17em}}\mathrm{122}:\mathrm{16}:\mathrm{1}:\mathrm{-}\mathrm{172}$) after Takahashi et al. (1985) and for the micro-nutrient iron ($\mathrm{Fe}:\mathrm{C}=\mathrm{4}\times {\mathrm{10}}^{-\mathrm{6}}:\mathrm{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₄, NO₃, O₂, Fe, SiO₂, H₂S and N₂) in porewater and solid sediment constituents (POC, opal, CaCO₃ and dust) as in the water column. Porewater tracers are exchanged with the overlying water column via diffusion. Sedimentation fluxes of POC, CaCO₃, 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 condition changes lead to chemical erosion of previously buried sediment.

Next we describe the changes introduced into HAMOCC as part of its implementation in CLIMBER-X.

N₂ fixation is represented by a diagnostic formulation, whereby the nitrate influx into the surface layer is a function of the nitrate deficit relative to phosphate, multiplied by a constant fixation rate (Ilyina et al., 2013). Prognostic N₂ 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 carbon is assumed to form aggregates with the denser calcite and opal built during phytoplankton and zooplankton growth and 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. The introduction of the ballasting scheme required a re-tuning of the dissolution rates of calcite and opal as shown in Table 1.

Following recent evidence that the remineralization of organic carbon depends on temperature (e.g. Laufkötter et al., 2017), we have introduced a Q₁₀ temperature dependence for the remineralization of POC and DOC (Segschneider and Bendtsen, 2013; Crichton et al., 2021), with a default Q₁₀ value of 2. The complete set of remineralization parameters is listed in Table 1.

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 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 nmol kg⁻¹ with a ligand stability constant of 1×10¹¹ kg mol⁻¹. 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 Hauck et al. (2013). Compared to HAMOCC, we have also increased the stoichiometric iron ratio in organic compounds from $\mathrm{Fe}:\mathrm{C}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}:\mathrm{1}$ to $\mathrm{Fe}:\mathrm{C}=\mathrm{4}\times {\mathrm{10}}^{-\mathrm{6}}:\mathrm{1}$. The parameters related to the iron cycle are also reported in Table 1.

Table 1 Modified HAMOCC parameters used in CLIMBER-X compared to HAMOCC6 (i.e. Table 2 in Ilyina et al., 2013).

The carbon-13 isotope was 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 into account the local concentration-dilution effect of the net surface freshwater flux, which changes surface DIC concentration and alkalinity. Two options are available in the model to implement the dilution effect on DIC and alkalinity. The first one ensures that the net global surface tracer flux is zero by applying deviations from the global average freshwater flux to the global average surface tracer concentration. The second (default) option applies the actual local surface freshwater flux to compute a new virtual top ocean layer thickness and then dilutes the tracers accordingly. In this case, the conservation of tracer inventories is ensured by compensating for imbalances over the global ocean. Additionally, during times when ocean volume is changing because of build-up or melt of land ice, concentrations of all tracers are globally adjusted while conserving tracer inventories. This is a reasonable simplification, considering that land ice volume changes occur on multi-millennial timescales, over which the ocean can be considered well mixed.

Based on scale analysis, we have excluded fast sinking tracers (CaCO₃, 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 line of thought, short-lived tracers like phytoplankton and zooplankton are also 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 1 d, which is also the time step of the physical ocean model.

2.2 Land carbon cycle: PALADYN

PALADYN is a comprehensive land surface–vegetation–carbon cycle model designed specifically for 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) and depends on absorbed shortwave radiation, air temperature, vapour pressure deficit between leaf and ambient air, atmospheric CO₂ and soil moisture. Carbon assimilation by vegetation is coupled to the transpiration of water through stomatal conductance. The model includes a dynamic vegetation module with five plant functional types (PFTs) competing for the grid-cell share based on their respective net primary productivity. The model distinguishes 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 and fast and slow carbon pools 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 permafrost. Carbon buried below ice sheets and on ocean shelves is treated separately. The land model also includes a dynamic peat module. PALADYN includes carbon isotopes ¹³C and ¹⁴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 set-up in Willeit and Ganopolski (2016), highlighted the need to improve certain aspects of the model. These improvements are described next.

We have updated the parameterization of the roughness length for heat and moisture. Originally, it was simply taken to be proportional to the roughness length for momentum, but there is ample evidence from observations that the roughness length for scalars can be orders of magnitude lower than that for momentum when the surface roughness is large (e.g. Zilitinkevich, 1995; Chen and Zhang, 2009; Yang et al., 2008; Zheng et al., 2012). We have therefore implemented the parameterization from Zilitinkevich (1995), which includes a dependence of the surface roughness length for heat and moisture on the roughness Reynolds number. With this new parameterization, the exchange coefficient for the turbulent surface fluxes shows a much weaker dependence on the roughness of the surface, which has an impact on the vegetation feedback.

We have introduced a topographic erodibility factor for dust emissions following Ginoux et al. (2001). It assumes that a basin with pronounced topographic variations contains a large amount of sediments which have accumulated in the valleys and depressions and which can easily be mobilized by wind. The following topographic factor is then applied to scale dust emissions:

$$\begin{array}{}\text{(1)}& f_{\mathrm{topo}}={\left({\displaystyle \frac{max(\mathrm{0},z_{\max }-z)}{z_{\max }-z_{\min }}}\right)}^{\mathrm{5}},\end{array}$$

where z is the grid-cell mean elevation, and z_max and z_min are the maximum and minimum surface elevations computed from the high-resolution topography in the surrounding 15×15^∘. The exponent 5 is taken from Zender et al. (2003).

The RuBisCO-limited photosynthesis rate in the version of the 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 historical period, which resulted in an overestimation of atmospheric CO₂ 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₂ concentration used for photosynthesis was computed based on the Cowan–Farquhar optimality hypothesis (Medlyn et al., 2011). In the new model version, for C₃ 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):

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }=\mathrm{4.4}{\displaystyle \frac{c_{\mathrm{a}}-c_{\mathrm{i}}}{c_{\mathrm{a}}}}+\mathrm{27}\cdot {\displaystyle \frac{c_{\mathrm{i}}}{c_{\mathrm{a}}}}-\mathrm{12}{\displaystyle \frac{{\mathrm{\Gamma }}_{\ast }}{p_{\mathrm{a}}}},\end{array}$$

where c_a and c_i are the ambient and leaf-internal CO₂ concentrations, p_a is the ambient partial pressure of CO₂ and Γ_∗ is the CO₂ compensation point.

For the distinction between evergreen and summergreen trees, in addition to a threshold on the coldest month's temperature, we have introduced a PFT-specific threshold on the growing degree days above 5 ^∘C, which is set to 600 for needleleaf trees and 900 for shrubs following Sitch et al. (2003).

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 used for growth for small leaf area index (LAI) values and all of 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 maximum 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:

$$\begin{array}{}\text{(3)}& \mathit{\lambda }={\displaystyle \frac{\mathrm{LAI}-{\mathrm{LAI}}_{\min }}{{\mathrm{LAI}}_{\max }-{\mathrm{LAI}}_{\min }}}.\end{array}$$

However, since the simulated leaf area index depends strongly on NPP, which in turn has a pronounced dependence on atmospheric CO₂, this formulation results in a strong dependence of λ on CO₂, with an increasingly larger fraction of NPP being used for spreading as CO₂ increases. We have therefore implemented a CO₂ dependence in the maximum leaf area index to reduce this effect:

$$\begin{array}{}\text{(4)}& {\mathrm{LAI}}_{\max }={\mathrm{LAI}}_{\max }^{\mathrm{ref}}\cdot \left(\mathrm{1}+\mathrm{0.5}\cdot \log \left({\mathrm{CO}}_{\mathrm{2}}/{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ref}}\right)\right).\end{array}$$

The fraction of decomposed litter respired directly as CO₂ 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 accumulating 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 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⁻¹ 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. Storage of wood from deforestation in products such as paper or wood for construction is not accounted for in the model and soil carbon is assumed to not be directly affected by land use practices. Following deforestation, the model will grow C₃ or C₄ grasses, depending on climate conditions.

The partitioning of the soil carbon decomposed under anaerobic conditions into CO₂ and CH₄ used a prescribed constant ratio in Willeit and Ganopolski (2016). We modified this by making the fraction released as CH₄ dependent on temperature with a Q₁₀ 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 (${\mathrm{HCO}}_{\mathrm{3}}^{-}$) (and therefore dissolved inorganic carbon and alkalinity) to the ocean and the consumption of atmospheric CO₂. The weathering rate depends on the lithology and on the climate variables temperature and runoff. The lithological map of Hartmann and Moosdorf (2012) 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 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 the 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 are taken from Börker et al. (2020). The weathering fluxes are transferred from the land to the ocean in the same way as water runoff, following the runoff routing scheme.

The carbon isotope fluxes from chemical weathering are computed assuming a δ¹³C of 1.8 ‰ 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 set-up, as they would result in further complications related to the conservation of nutrients in the ocean. Instead, as discussed in Sect. 2.1, the silicate and phosphorus budgets are closed by assuming that the sediment burial flux is returned as input at the ocean surface.

2.3 Atmospheric CO₂

The atmospheric CO₂ 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₂ in the atmosphere (C_atm):

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{atm}}}{\mathrm{d}t}}=F_{\mathrm{ocn}}+F_{\mathrm{lnd}}+F_{\mathrm{anth}}-F_{\mathrm{weath}}+F_{\mathrm{volc}}+F_{{\mathrm{CH}}_{\mathrm{4}}\mathrm{ox}}.\end{array}$$

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₂ consumption by silicate and carbonate weathering, the volcanic degassing flux and the CO₂ flux from the oxidation of atmospheric methane originating from non-agricultural sources. The CO₂ consumption by weathering is computed assuming that all carbon in the ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ originating from the weathering of silicate rocks ($F_{{\mathrm{HCO}}_{\mathrm{3}}^{-}}^{\mathrm{sil}}$) comes from the atmosphere, while only half of the carbon in the ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ originating from the weathering of carbonate rocks and sediments ($F_{{\mathrm{HCO}}_{\mathrm{3}}^{-}}^{\mathrm{carb}}$) comes from the atmosphere:

$$\begin{array}{}\text{(6)}& F_{\mathrm{weath}}=F_{{\mathrm{HCO}}_{\mathrm{3}}^{-}}^{\mathrm{sil}}+\mathrm{0.5}\cdot F_{{\mathrm{HCO}}_{\mathrm{3}}^{-}}^{\mathrm{carb}}.\end{array}$$

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 spin-up simulation:

$$\begin{array}{}\text{(7)}& F_{\mathrm{volc}}=\mathrm{0.5}\cdot F_{{\mathrm{HCO}}_{\mathrm{3}}^{-}}^{\mathrm{sil}}.\end{array}$$

The flux from the oxidation of methane, $F_{{\mathrm{CH}}_{\mathrm{4}}\mathrm{ox}}$, is computed by the CH₄ model as described in Sect. 2.4 below. The atmospheric CO₂ concentration is then computed from C_atm using a conversion factor of 2.12 PgC ppm⁻¹ (Denman et al., 2007).

Equations similar to Eq. (5) are also used for the carbon isotopes ¹³C and ¹⁴C. The prognostic equation for the stable isotope ¹³C in atmospheric CO₂ is

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{13}}{C}_{\mathrm{atm}}}{\mathrm{d}t}}=F_{\mathrm{ocn}}^{\mathrm{13}}+F_{\mathrm{lnd}}^{\mathrm{13}}+F_{\mathrm{anth}}^{\mathrm{13}}-F_{\mathrm{weath}}^{\mathrm{13}}+F_{\mathrm{volc}}^{\mathrm{13}}.\end{array}$$

The ¹³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 δ¹³C of anthropogenic carbon emissions is prescribed as time-dependent from historical data of Andres et al. (2017), and the ¹³C flux from CO₂ consumption by weathering, assuming no fractionation, is simply computed as

$$\begin{array}{}\text{(9)}& F_{\mathrm{weath}}^{\mathrm{13}}=F_{\mathrm{weath}}{\displaystyle \frac{{}^{\mathrm{13}}{C}_{\mathrm{atm}}}{C_{\mathrm{atm}}}}.\end{array}$$

The ¹³C of volcanic degassing is computed assuming a δ¹³C of −5 ‰.

The prognostic equation for radiocarbon ¹⁴C in atmospheric CO₂ reads as

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{14}}{C}_{\mathrm{atm}}}{\mathrm{d}t}}=F_{\mathrm{ocn}}^{\mathrm{14}}+F_{\mathrm{lnd}}^{\mathrm{14}}-F_{\mathrm{weath}}^{\mathrm{14}}+F_{\mathrm{prod}}^{\mathrm{14}}-{\displaystyle \frac{{}^{\mathrm{14}}{C}_{\mathrm{atm}}}{{\mathit{\tau }}_{{}^{\mathrm{14}}\mathrm{C}}}}.\end{array}$$

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 ¹⁴C. The production rate of radiocarbon in the atmosphere ($F_{\mathrm{prod}}^{\mathrm{14}}$) is prescribed in the model and the radiocarbon decay time is ${\mathit{\tau }}_{{}^{\mathrm{14}}\mathrm{C}}=\mathrm{8267}$ years.

2.4 Atmospheric CH₄

Similarly to CO₂, atmospheric CH₄ is also considered to be well mixed in the atmosphere and is therefore represented as a globally uniform value. The atmospheric CH₄ concentration can be prescribed, or it can be interactively computed by the model from

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}{\mathrm{CH}}_{\mathrm{4}}}{\mathrm{d}t}}=F_{\mathrm{lnd}}^{\mathrm{emis}}+F_{\mathrm{anth}}^{\mathrm{emis}}-{\displaystyle \frac{{\mathrm{CH}}_{\mathrm{4}}}{{\mathit{\tau }}_{{\mathrm{CH}}_{\mathrm{4}}}}}.\end{array}$$

Methane sources include natural emissions from wetlands and peatlands ($F_{\mathrm{lnd}}^{\mathrm{emis}}$), 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. Anthropogenic methane emissions ($F_{\mathrm{anth}}^{\mathrm{emis}}$) are prescribed in the model. The sink of methane from oxidation in the atmosphere is computed using a constant residence time of CH₄, ${\mathit{\tau }}_{{\mathrm{CH}}_{\mathrm{4}}}=\mathrm{9.5}$ years, which is a reasonable first approximation, at least for climate conditions ranging between the Last Glacial Maximum and the present day (Kleinen et al., 2020; Levine et al., 2011; Hopcroft et al., 2017).

3 Closed and open carbon cycle model configurations and model spin-up

Two different configurations of the carbon cycle model are available and can be chosen according to the specific needs.

The first (and simplest) set-up consists of ocean, land and atmosphere carbon cycle components only. In this set-up marine sediments are disabled and particulate fluxes that reach the ocean floor are completely remineralized/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 set-up the carbon system is closed 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 set-up is equivalent to what is used in many state-of-the-art Earth system models for climate change projections on centennial timescales (e.g. Séférian et al., 2020). The model spin-up for this simple set-up is straightforward and requires only that the model is run to steady state with a prescribed atmospheric CO₂ concentration for ≈10 000 years. The slowest timescale in this set-up 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 spin-up 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 set-up is applicable to simulations of up to 1000 years. On longer timescales, sediment and weathering processes become important and need to be accounted for when performing long-term transient simulations with interactive CO₂. 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 pre-industrial is an equilibrium state of the climate–carbon cycle system allows us to run perturbation experiments with the interactive carbon cycle without having to deal with possible long-term drifts in atmospheric CO₂. However, the long timescale 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 set-up with prescribed climatological daily input fields at the ocean surface. This set-up results in a speed-up of a factor >2 relative to running the fully coupled climate–carbon cycle model, meaning that ocean carbon cycle and marine sediments can be run into equilibrium in about a week of computing time on a high-performance computer. In detail, the spin-up procedure of the full carbon cycle configuration comprises two different stages. Atmospheric CO₂ is prescribed to a constant value throughout the process, at 280 ppm for the pre-industrial case. The first stage aims at spinning up the sediment model. For this purpose the full carbon cycle–climate model is run for 5000 years, and every 300 years the sediment model is run offline for 1000 years. During this stage all net fluxes into the sediments are compensated for and returned as inputs at the ocean surface in order to approximately conserve 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 set-up 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 set-up is that organic carbon and opal that are buried in the sediments, and are therefore effectively leaving the system, are returned in remineralized form to the surface ocean, so that phosphorus and silica inventories of the ocean–sediment system are conserved throughout the simulation.

Figure 2 Open versus closed carbon cycle spin-up for pre-industrial conditions. The figure shows surface input of (a) DIC and (b) alkalinity, the evolution of (c) DIC and (d) alkalinity inventories in the ocean, and (e) the air–sea CO₂ flux. The grey vertical lines indicate the switch between the first and second spin-up phases, as described in the text.

The carbon fluxes among the different model components in the open set-up for equilibrium pre-industrial conditions are schematically illustrated in Fig. 3. The volcanic degassing rate is equal to half the atmospheric CO₂ 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 that the ocean alkalinity budget also is.

Figure 3 CLIMBER-X carbon fluxes and reservoirs in equilibrium with pre-industrial conditions for the open carbon cycle set-up.

4 Model evaluation for the historical period and present day

Here we present results from a CLIMBER-X simulation with interactive CO₂ and CH₄ in the open carbon cycle set-up for the historical period (1850–2015) and provide a comprehensive evaluation of model performance against various observational datasets. The forcings for this simulation include variations in solar radiation (Matthes et al., 2017), radiative forcing of volcanic eruptions (Prather et al., 2013), globally uniform N₂O concentrations from Köhler et al. (2017), globally uniform CFC11 and CFC12 concentrations from Meinshausen et al. (2017), 3D O₃ concentrations and 2D ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ load from the ensemble mean of CMIP6 models and land use change (pasture and cropland fractions) from Ma et al. (2020). The model is initialized from an 80 000-year equilibrium simulation with the open carbon cycle set-up for pre-industrial boundary conditions and a prescribed atmospheric CO₂ of 280 ppm, as described in Sect. 3 and shown in Fig. 2.

4.1 Present day

In the following, different simulated climatological characteristics are compared to observations to assess the model performance for the present day. Unless stated otherwise, the comparison with observations is for the time interval from 1981 to 2010. To 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 CMIP6 (Eyring et al., 2016). The following CMIP6 models are included for ocean biogeochemistry: CESM2, IPSL-CM6A-LR, MRI-ESM2-0, MIROC-ES2L, MPI-ESM1-2-LR, UKESM1-0-LL and CanESM5. For the land carbon cycle, the following models are used for comparison: ACCESS-ESM1-5, BCC-CSM2-MR, CanESM5, CNRM-ESM2-1, GFDL-ESM4, IPSL-CM6A-LR, MIROC-ES2L, MPI-ESM1-2-LR, MRI-ESM2-0, NorESM2-LM and UKESM1-0-LL. For ocean biogeochemistry, we highlight how the model compares with results from the MPI-ESM1-2-LR employing the original marine carbon cycle model HAMOCC6.

4.1.1 Ocean biogeochemistry and marine sediments

An overview of simulated global variables characterizing the ocean carbon cycle are presented and compared to observation-based estimates in Table 2, providing a summary of model performance for the present day.

Jacobson et al. (2007)Regnier et al. (2013)Buitenhuis et al. (2018)(Lauvset et al., 2016; Olsen et al., 2016)Garcia et al. (2013b)Garcia et al. (2013b)Garcia et al. (2013b)Johnson and Bif (2021)Carr et al. (2006)Karl et al. (2002)Großkopf et al. (2012)Dunne et al. (2007)Dunne et al. (2007)Tréguer and De La Rocha (2013)Dunne et al. (2007)Battaglia et al. (2016)Milliman and Droxler (1996)Tréguer and De La Rocha (2013)Tréguer et al. (2021)Cartapanis et al. (2018)Cartapanis et al. (2018)Tréguer and De La Rocha (2013)

Table 2 Global values of the main ocean biogeochemical variables for the present day.

Representing the ocean ventilation timescale reasonably well is a prerequisite for simulating biogeochemical tracers in the ocean. The ocean uptake of CFCs of anthropogenic origin over the historical period is often used to probe the ventilation of the ocean on decadal timescales, while the pre-industrial radiocarbon concentration in the ocean provides information on the age distribution of the water masses in an approximate equilibrium state. We therefore start by comparing how well the model reproduces the CFC11 and radiocarbon distributions in the ocean. The inventory of CFC11 in the ocean starts to increase after ≈1950 as a consequence of its increase in the atmosphere (Fig. 4). Estimates for the CFC11 inventory in the year ≈1994 are available from models from the OCMIP model intercomparison (Dutay et al., 2002) and from direct observations (Willey et al., 2004). CLIMBER-X results are generally consistent with these estimates (Fig. 4), indicating that, at least at the global scale, the decadal ventilation timescale in CLIMBER-X is well in line with observations and other models. In terms of spatial distribution, the CFC11 uptake is overestimated in the North Pacific Ocean, in the North Indian Ocean and around Antarctica, while too small CFC11 concentrations are simulated at mid latitudes in all basins at depths between 500 and 1000 m (Fig. 5).

Figure 4 Historical global cumulative ocean uptake of CFC11 in CLIMBER-X compared to observations (Willey et al., 2004) and OCMIP models (Dutay et al., 2002).

Figure 5 Zonally averaged CFC11 concentration for the year 1994 in CLIMBER-X (a, d, g) and GLODAP (Key et al., 2004) (b, e, h) for different basins: Atlantic (a–c), Pacific (d–f), and Indian (g–i) oceans. The model bias is shown in panels (c), (f), and (i).

The radiocarbon ventilation age in the pre-industrial gives additional insights into the ocean ventilation under quasi-equilibrium conditions, information which is complementary to CFC11. The radiocarbon ventilation age of the deep ocean is nicely reproduced by CLIMBER-X, while radiocarbon age is systematically overestimated in the upper kilometre across all ocean basins (Fig. 6). 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 and therefore not representing the non-linear effects of synoptic variability on vertical mixing of tracers. For instance, a one-time mixing down to 200 m depth by a wind storm could have a large effect on some tracers, which cannot be resolved by using climatological mean winds. We would expect this non-linear effect to be much more important for radiocarbon than for nutrients. The analyses of CFC11 and radiocarbon provide important insights into the ocean ventilation in the model and will be useful when discussing model biases in the distribution of other biogeochemical tracers below.

Figure 6 Zonally averaged pre-industrial radiocarbon ventilation age in CLIMBER-X (a, d, g) and GLODAP Key et al. (2004) (b, e, h) for different basins: Atlantic (a–c), Pacific (d–f), and Indian (g–i) oceans. The model bias is shown in panels (c), (f), and (i).

The spatial pattern of the air–sea CO₂ exchange is well captured by the model (Fig. 7), with outgassing generally taking place in the tropics and CO₂ being taken up at mid to high northern latitudes and at mid latitudes of the Southern Hemisphere. The main difference compared to the other models is observed around the Equator, with a less pronounced peak in CO₂ release simulated by CLIMBER-X (Fig. 8), 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₂ uptake compared to observations (e.g. Gruber et al., 2009) (Fig. 8), while CLIMBER-X is apparently more consistent with recent estimates, although with substantial differences in the spatial distribution of the CO₂ flux (Fig. 7). Notably, in the Southern Ocean the CLIMBER-X air–sea CO₂ exchange diverges from that simulated by the MPI-ESM1-2-LR model (Fig. 8), 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. 9a). However, the MPI-ESM seems to be an outlier in the simulated primary production in the Southern Ocean, possibly because of biases in climate, which are unrelated to the HAMOCC ocean carbon cycle model.

Figure 7 The 1985–2010 average air–sea CO₂ flux in (a) CLIMBER-X compared to (b) observations from Landschützer et al. (2016).

Figure 8 Zonal mean air–sea CO₂ flux (1985–2010 average) in CLIMBER-X compared to observations from Landschützer et al. (2016) and selected CMIP6 models.

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 remineralization rates of POC in the water column. While the net primary productivity in CLIMBER-X is in line with CMIP6 models (Fig. 9a) and the globally integrated value of 53 PgC yr⁻¹ agrees well with observations (Table 2), the export production in the model is generally at the lower end of the CMIP6 model range (Fig. 9b). CaCO₃ and opal export are compared to CMIP6 models in Fig. 9c, d.

Figure 9 The 1981–2010 average global zonal mean (a) net primary production, (b) particulate organic carbon export at 100 m depth, (c) CaCO₃, and (d) opal export at 100 m depth. Results from CLIMBER-X are compared to CMIP6 models.

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. 10). However, over the Southern Ocean, in the equatorial Pacific and in the North Pacific, production is limited by the availability of iron (Fig. 10). This is in accordance with observations showing that iron limitation is usually important where sub-surface 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 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. 11), 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 the dust deposition rate (Fig. 12).

Figure 10 Nutrient limitation of marine net primary productivity in CLIMBER-X.

Figure 11 (a) CLIMBER-X annual dust deposition flux compared to model-based products of (b) Albani et al. (2016), (c) Lambert et al. (2015), and (d) Hopcroft et al. (2015). The respective globally integrated deposition values are given in brackets in the panel titles.

Figure 12 Simulated versus observed dust deposition fluxes at different locations available from the AeroCom dataset (Huneeus et al., 2011, and references therein). The dashed lines indicate 1 order of magnitude deviation from the 1:1 line.

The simulated dissolved iron concentration in surface water is closely related to the dust deposition shown in Fig. 11. It is therefore high in the Atlantic and Indian oceans, lower in the Southern Ocean and very small over large parts of the Pacific (Fig. 13). This is broadly consistent with observations (e.g. Tagliabue et al., 2012), but measurements of iron concentration in ocean water are still relatively sparse.

Figure 13 Average surface dissolved iron concentration in CLIMBER-X over the period 1981–2010.

The main features of the surface nitrate concentration are well reproduced by CLIMBER-X, with large concentrations in the Southern Ocean, moderate values in the upwelling region of the eastern equatorial Pacific and in the North Atlantic and North Pacific and low values elsewhere (Figs. 15 and 14a). 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 in very good agreement with observations (Fig. 16).

Figure 14 The 1981–2010 average global zonal mean surface concentrations of the nutrients (a) nitrate, (b) phosphate, (c) silicate, and (d) dissolved iron. CLIMBER-X is compared to observations (Garcia et al., 2013b) and CMIP6 model results.

Figure 15 Surface NO₃ concentration in (a) CLIMBER-X (1981–2010 average) compared to (b) observations from the World Ocean Atlas 2013 (WOA13, Garcia et al., 2013b). The model bias is shown in panel (c).

Figure 16 Global and basin-wide average profiles of different biogeochemical tracers in the ocean, from top to bottom: DIC, alkalinity, nitrate, phosphate, silicate, and dissolved iron. CLIMBER-X results (black) are compared to observations (blue) (Lauvset et al., 2016; Olsen et al., 2016; Garcia et al., 2013a, b) and CMIP6 model results (grey). Results from the MPI-ESM1-2-LR are shown by the green dashed lines. The boundary of the Southern Ocean is set at 35 ^∘ S and the Southern Ocean section is not included in the profiles of the Atlantic, Pacific, and Indian oceans. CLIMBER-X and CMIP6 data are averages over the time period 1981–2010.

The 3D phosphate distribution in the global ocean is nicely captured by the model (Figs. 17, 14b, and 16), except for too low concentrations simulated in the surface ocean of the North Pacific and North Indian oceans. The negative bias in 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 kilometre 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. 14a, b), as fewer nutrients are assimilated during photosynthesis.

Figure 17 Zonally averaged PO₄ concentration in CLIMBER-X, 1981–2010 average (a, d, g), and WOA13 (Garcia et al., 2013b) (b, e, h) for different basins: Atlantic (a–c), Pacific (d–f), and Indian (g–i) oceans. The model bias is shown in panels (c), (f), and (i).

Silicate concentration is generally overestimated in the sub-surface ocean and is underestimated in the deep North Pacific and North Indian oceans (Fig. 18), similarly to other nutrients (Fig. 17).

Figure 18 Zonally averaged Si concentration in CLIMBER-X, 1981–2010 average (a, d, g), and WOA13 (Garcia et al., 2013b) (b, e, h) for different basins: Atlantic (a–c), Pacific (d–f), and Indian (g–i) oceans. The model bias is shown in panels (c), (f), and (i).

The large-scale patterns of oxygen concentration in ocean waters simulated by CLIMBER-X are largely consistent with observations (Fig. 19), but the extent and depth of the oxygen minimum zones, in particular in the eastern equatorial Pacific, are overestimated. This bias is common to many CMIP5 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.

Figure 19 Oxygen concentration in CLIMBER-X, 1981–2010 average (left column), and WOA13 (Garcia et al., 2013a) (middle column) at different ocean depths: from top to bottom, 100, 500, 1000, and 3000 m. The model bias is shown in the right column.

Both DIC and alkalinity are generally overestimated in the upper ocean (Fig. 16), particularly in Antarctic intermediate water masses, and underestimated in the deep ocean (Figs. 20 and 21). These biases in the simulated vertical distribution of DIC and alkalinity could be due to a relatively low CaCO3 export from the euphotic layer (Table 2), which leads to a too weak vertical redistribution. Additionally, the simulated DIC concentration is generally too low in the North Pacific and North Indian oceans.

Figure 20 Zonally averaged dissolved inorganic carbon in CLIMBER-X, 1981–2010 average (a, d, g), and GLODAPv2 (Lauvset et al., 2016; Olsen et al., 2016) (b, e, h) for different basins: Atlantic (a–c), Pacific (d–f), and Indian (g–i) oceans. The model bias is shown in panels (c), (f), and (i).

Figure 21 Zonally averaged total alkalinity in CLIMBER-X, 1981–2010 average (a, d, g), and GLODAPv2 (Lauvset et al., 2016; Olsen et al., 2016) (b, e, h) for different basins: Atlantic (a–c), Pacific (d–f), and Indian (g–i) oceans. The model bias is shown in panels (c), (f), and (i).

The carbon-13 isotope in the ocean helps to track the distribution of different water masses. The higher δ¹³C values in the Atlantic Ocean 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. 22). The negative biases at 500–1500 m depth are associated with the “nutrient trapping” problem (Aumont et al., 1999; Dietze and Loeptien, 2013) that is often seen in ESMs. This problem is characterized by high concentrations of remineralized nutrients and carbon and, therefore, low δ¹³C (Liu et al., 2021). The positive biases through the whole water column in the North Atlantic, North Pacific and North Indian oceans are possibly the result of too strong ventilation in these regions in the model.

Figure 22 Zonally averaged δ¹³C in CLIMBER-X, 1981–2010 average (a, d, g), and Eide et al. (2017) (b, e, h) for different basins: Atlantic (a–c), Pacific (d–f), and Indian (g–i) oceans. The model bias is shown in panels (c), (f), and (i).

Figure 23 Weight fraction of calcite, opal, and organic carbon in marine sediments as simulated by CLIMBER-X (a–c) compared to observations (Hayes et al., 2021) (d–f).

In the Atlantic and Indian oceans, CaCO₃ dominates the sediment composition, in accordance with observations (Fig. 23a, d). However, little CaCO₃ is simulated in large parts of the sediment in the eastern Pacific Ocean, where observations indicate widespread CaCO₃ content in the Southern Hemisphere (Fig. 23a, d). The underestimation of calcite weight fractions in sediments of the eastern South Pacific Ocean is caused by water being undersaturated with respect to calcite in this area. This leads to dissolution of most of the calcite produced at the surface before it can even reach the sediments. The strongly undersaturated water is ultimately a result of deficiencies in the simulated ocean circulation. Some other models show similar deficiencies in the simulated calcite fraction in Pacific sediments (e.g. Kurahashi-Nakamura et al., 2022). Global CaCO₃ sediment deposition and burial are in line with observational underestimates (Table 2), with around 25 % of the deposited CaCO₃ undergoing dissolution. The opal content in sediments in CLIMBER-X is overestimated (Fig. 23b, e), even though the global opal sedimentation and burial fluxes are fully consistent with observational estimates (Table 2). Opal is particularly abundant in the eastern equatorial Pacific simply as a result of missing CaCO₃ in the sediments in that area. Organic carbon is found mainly on the continental margins and in the equatorial eastern Pacific, in agreement with observations (Fig. 23c, f), although CLIMBER-X tends to underestimate the organic carbon content in sediments, possibly because of a too small sediment deposition flux of POC (Table 2).

4.1.2 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 coupled climate model set-up and with the additional modifications to the model described above.

A selection of simulated global variables characterizing the land carbon cycle is presented and compared to observation-based estimates in Table 3, providing a summary of model performance for the present day.

Beer et al. (2010)Ito (2011)Fan et al. (2020)Fan et al. (2020)Varney et al. (2022)Varney et al. (2022)Brown et al. (1998)Tarnocai et al. (2009)Hugelius et al. (2014)Yu et al. (2010)Yu et al. (2010)Prigent et al. (2007)Papa et al. (2010)Saunois et al. (2020)Saunois et al. (2020)Saunois et al. (2020)Munhoven (2002)Munhoven (2002)Munhoven (2002)Amiotte Suchet et al. (2003)Gaillardet et al. (1999)

Table 3 Global values for the main variables of the land carbon cycle.

Figure 24 (a) Simulated GPP compared to (b) observations (Jung et al., 2011). (c) Comparison of zonally integrated GPP. (d) Simulated vegetation carbon compared to (e) observations (Carvalhais et al., 2014). (f) Comparison of zonally integrated vegetation carbon. (g) Simulated soil carbon compared to (h) observations (Carvalhais et al., 2014). (i) Comparison of zonally integrated soil carbon. (j) Simulated ecosystem carbon turnover time compared to (k) observations (Fan et al., 2020). (l) Comparison of zonal mean ecosystem carbon turnover time. In panels (c), (f), (i), and (l), results from CLIMBER-X are shown in black, observations in blue, and CMIP6 models in grey. CLIMBER-X and CMIP6 data are averages over the time period 1981–2010.

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 integrals (Table 3) and in terms of spatial distribution (Fig. 24a, b, c).

The total carbon stored in the vegetation, both above ground and below ground, is slightly overestimated in the model (Table 3), but the meridional distribution, mainly originating from large-scale differences in precipitation, is well reproduced (Fig. 24d, e, f). Most of the soil carbon in CLIMBER-X is stored in cold soils of the Northern Hemisphere high latitudes, in agreement with observations (Fig. 24g, h, i). However, compared to estimates from Carvalhais et al. (2014), the soil carbon distribution is too skewed toward high northern latitudes, and there is too little carbon in the tropics. Most CMIP6 models underestimate soil carbon in the tropics as well (Fig. 24j).

In CLIMBER-X, ∼ 1500 PgC of carbon is stored in the top soil metre, in good agreement with different estimates (Table 3). However, with ∼ 2150 PgC, the total soil carbon content seems to be underestimated compared to observations, which suggest >3000 PgC. This indicates that too little carbon is simulated in soil below 1 m depth. However, total soil carbon content estimates vary widely between datasets (e.g. Fan et al., 2020), with e.g. 1952±198 PgC in WISE30sec (Batjes, 2016) and 3141±893 PgC in Sanderman et al. (2017) in the top 2 m of soil. Most of the carbon in mineral soil layers below 1 m is recalcitrant, and its response to changes in environmental conditions is uncertain. In Earth system models, total soil carbon storage is usually much lower (1206±445 PgC, Varney et al., 2022), as these models account for active carbon responding on centennial timescales. In CLIMBER-X, one possible explanation for underestimated carbon content in deeper soil layers is that the maximum turnover timescale of soil carbon is set to 5000 years in the model, which limits the amount of carbon that can be accumulated in cold, frozen soil layers. Other possible reasons include (i) a general underestimation of vertical carbon transport by diffusion, particularly into perennially frozen soil layers, and (ii) a possible depth dependence of soil carbon turnover due to processes other than temperature and moisture (e.g. Koven et al., 2013) and that are not included in the model. Consistently, the carbon contained in areas affected by permafrost is ∼ 800 PgC, which is also a bit lower than the ∼ 1100–1500 PgC suggested by observations (Table 3). Let us note that, even when models are initialized with the observed permafrost carbon stock of ≈1300 PgC, remapping at model resolution and accounting for differences in soil temperatures between models and observations generally leads to a reduction in permafrost carbon stocks (e.g. Kleinen and Brovkin, 2018). The CLIMBER-X-simulated peatland extent is lower than estimated (Yu et al., 2010), and the peat carbon is also consistently underestimated accordingly (Table 3).

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 turnover time simulated by CLIMBER-X is in line with CMIP6 models, while it is underestimated compared to observation-based estimates from Fan et al. (2020) (Fig. 24j, 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 (Table 3), although with substantial differences in the geographic distribution (Fig. 25). Compared to the multi-satellite product from GIEMS (Global Inundation Extent from Multi-Satellites) (Prigent et al., 2007; Papa et al., 2010), the model simulates a larger wetland extent in tropical forest areas. However, when compared to other wetland products based on data other than from satellites, GIEMS underestimates 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.

Figure 25 Maximum monthly wetland fraction (a) in CLIMBER-X compared to (b) the GIEMS dataset (Papa et al., 2010; Prigent et al., 2007).

In CLIMBER-X, methane is emitted exclusively from wetlands. Because of the dependence of methane emissions 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 (Table 3, Fig. 26), in agreement with observations (e.g. Saunois et al., 2020). The total CH₄ 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.

Figure 26 Natural methane emission simulated by CLIMBER-X for the present day (1981–2010 average).

Chemical weathering fluxes are generally high where runoff is high, with the separation between silicate and carbonate weathering being modulated by lithological properties (Fig. 27). The global CO₂ consumption rate by weathering and the alkalinity flux to the ocean in the form of bicarbonate produced by rock weathering are in good agreement with observational estimates (Table 3), while the partitioning between carbonate and silicate weathering is skewed toward carbonate weathering (Table 3).

Figure 27 CLIMBER-X (a) silicate and (b) carbonate weathering flux distribution for the present day (1981–2010 average).

4.2 Historical period

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₂ concentration is well reproduced by the model, with CO₂ at the year 2015 being within ∼ 5ppm of direct measurements (Fig. 28). Biases in simulated CO₂ of ∼ 10 ppm are quite common in state-of-the-art ESMs (e.g. Hoffman et al., 2014; Friedlingstein et al., 2014).

Figure 28 Historical atmospheric CO₂ concentration from a coupled CLIMBER-X simulation compared to observations (Köhler et al., 2017).

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. 29. 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 with the GCB, although the actual values remain uncertain (e.g. Gasser et al., 2020). A substantial fraction of this anthropogenic CO₂ 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₂ concentrations and the concentration of dissolved CO₂ in the surface ocean water and is relatively well understood, as also indicated by the narrow uncertainty range obtained from different CMIP6 models (Fig. 30a). 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₂. 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. 30b). However, the effect of this NPP increase on vegetation carbon varies widely among models (Fig. 30c), also because of the confounding factor of land use change. In CLIMBER-X the net effect is a vegetation carbon stock decrease by ∼ 10 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 dominate over the negative effect of increasing temperatures in CLIMBER-X, leading to an increase in soil carbon by ∼ 50 PgC (Fig. 30d).

Figure 29 Historical global carbon budget in CLIMBER-X. The dashed lines are estimates from Friedlingstein et al. (2022).

Figure 30 Historical anomalies of (a) air–sea CO₂ flux, (b) net primary production on land, (c) vegetation carbon, and (d) soil carbon in CLIMBER-X compared to CMIP6 models. The anomalies are computed relative to the time interval 1850–1880 CE.

Since CLIMBER-X is enabled with carbon isotopes, it also allows for a comparison of isotopic signatures to observations, thereby providing additional constraints on processes involved in carbon cycle exchanges. As an example, the historical δ¹³C of atmospheric CO₂ is compared to observations in Fig. 31.

Figure 31 Historical δ¹³C of atmospheric CO₂ in CLIMBER-X compared with observations (Eggleston et al., 2016).

The general historical trend in atmospheric CH₄ is captured by the model (Fig. 32a). Prescribed anthropogenic methane emissions are the dominant source for the increase in the atmospheric methane burden, but natural emissions from land are also increasing due to the increase in NPP and soil temperature (Fig. 32b).

Figure 32 Historical (a) atmospheric CH₄ concentration in CLIMBER-X compared to observations (Köhler et al., 2017) and (b) prescribed anthropogenic CH₄ emissions and natural land emissions as simulated in CLIMBER-X.

5 Carbon cycle feedbacks

Carbon cycle processes both on land and in the ocean are sensitive to changes in climate and atmospheric CO₂. This implies that carbon fluxes between ocean and atmosphere and between land and atmosphere will respond to changes in climate and CO₂ concentration, which will in turn affect CO₂ and consequently climate. Quantifying the strength of these feedbacks is important to understanding how the climate will respond to anthropogenic CO₂ emissions. For that, a linear feedback decomposition analysis was 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⁻¹ CO₂ increase experiments, whereby in one simulation the CO₂ increase is seen only by the radiative code in the atmosphere (radiatively coupled), in another one the CO₂ increase is seen only by the carbon cycle (biogeochemically coupled) and in a final one both the radiation and carbon cycle see the increase in atmospheric CO₂ (fully coupled). This set of simulations allows us to roughly separate the effect of changes in climate and changes in CO₂ on 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 the C⁴MIP, CMIP5 and CMIP6 models (Friedlingstein et al., 2006; Arora et al., 2013, 2020). In Fig. 33 the carbon cycle–climate (γ) and the carbon cycle–concentration (β) feedbacks in CLIMBER-X are compared to CMIP6 model results separately for land and ocean. An increase in CO₂ has a positive effect on the uptake of carbon by both land and ocean, resulting in a lowering of CO₂ and therefore a negative feedback on climate (implying positive β, Fig. 33a, c). Conversely, 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. 33b, d). The β and γ values obtained here are well within the range of estimates from CMIP6 models (Arora et al., 2020), although in CLIMBER-X the CO₂ fertilization effect on land is rather high (Fig. 33a) and the ocean carbon uptake as a response to an increase in atmospheric CO₂ is at the lower end of the CMIP6 models (Fig. 33c).

Figure 33 Carbon cycle feedbacks in CLIMBER-X compared to CMIP6 models (Arora et al., 2020). The (a, c) β and (b, d) γ parameters are shown at the time of CO₂ doubling (year ∼ 70) and at the time of CO₂ quadrupling (year ∼ 140) in idealized 1 % yr⁻¹ CO₂ increase experiments.

6 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₂ emissions (e.g. Matthews and Solomon, 2013). A model intercomparison project has been established in an effort to analyse and compare the ZEC in different Earth system models (Jones et al., 2019). Here we use this standardized and idealized experimental set-up to compare the carbon cycle response in CLIMBER-X with results from the Zero Emissions Commitment Model Intercomparison Project (ZECMIP) models for the 1000 PgC emission pulse (MacDougall et al., 2020). The experiment branches off from a 1 % per year CO₂ increase run with CO₂ 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 set-ups.

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. 34). Atmospheric CO₂ concentration rapidly decreases after stopping the carbon emissions (Fig. 34b), while global temperature shows a more modest decrease (Fig. 34a). The ocean continues to take up carbon throughout the simulation (Fig. 34c), while the land turns from a sink to a source of carbon roughly at the time of the peak in CO₂ (Fig. 34d). CLIMBER-X initially shows a relatively weak ocean carbon uptake compared to ZECMIP models, while the land carbon uptake is at the high end of the ZECMIP model ensemble.

The differences between the experiments with open and closed carbon cycle set-ups are negligible for the first few centuries but continue to grow over time, with CO₂ decreasing faster in the open set-up (Fig. 34b).

Figure 34 Comparison of CLIMBER-X simulations with ZECMIP model results in terms of (a) global temperature, (b) atmospheric CO₂ concentration, (c) cumulative ocean carbon uptake, and (d) cumulative land carbon uptake for the standard ZECMIP experiment with 1000 PgC cumulative carbon emissions. For CLIMBER-X, results with both the open (solid lines) and closed (dashed lines) carbon cycle set-ups are shown.

7 Conclusions

We have described the major features of the carbon cycle components of the newly developed CLIMBER-X Earth system model. The model includes a detailed representation of carbon cycle processes on land, in the ocean and in marine sediments, thus allowing the investigation of the complex interactions between climate and the carbon cycle. Two set-ups of the global carbon cycle, closed and open, are available in CLIMBER-X, allowing both a proper comparison with CMIP6-type models in centennial-scale and multi-millennium simulations. We have evaluated the model performance for the historical period against observations and state-of-the-art CMIP6 models, showing that many characteristics and feedbacks of the carbon cycle are reasonably well captured by the model. Biases in the simulated distribution of ocean biogeochemical tracers exist and can mostly be related to deficiencies in the simulated ocean circulation changes, which can at least partly be attributed to the comparatively coarse resolution of the ocean model and to the frictional–geostrophic approximation that it rests upon. On land, the carbon in soil layers below 1 m is underestimated, particularly in the permafrost zone, with possible implications for the land carbon cycle response to global warming.

Some possible directions for future model developments include the extension of the organic carbon cycle, allowing for burial on land and in sediments and fluxes from land to ocean, the refinement of the carbonate chemistry on shelves, including coral growth, and possibly the addition of the nitrogen cycle on land, which could be important for nitrogen limitation of photosynthesis and would allow for interactive atmospheric N₂O, considering that N₂O fluxes from the ocean are already available from the ocean biogeochemistry model HAMOCC.

The computational efficiency of CLIMBER-X will enable it to be used for systematic explorations of the coupled climate–carbon cycle evolution over timescales ranging from decades up to ∼ 100 000 years while also allowing a quantification of related uncertainties.