Assessing the lifetime of anthropogenic CO₂ and its sensitivity to different carbon cycle processes

Although it is well-established that anthropogenic CO₂ emitted into the atmosphere will persist for a long time, the duration of the anthropogenic climate perturbation will depend on how rapidly the excess CO₂ is removed from the climate system by different biogeochemical processes. The uncertainty around the long-term climate evolution is therefore linked to not only the future of anthropogenic CO₂ emissions but also our insufficient understanding of the long-term carbon cycle. Here, we use the fast Earth system model CLIMBER-X, which features a comprehensive carbon cycle, to examine the lifetime of anthropogenic CO₂ and its effects on the long-term evolution of atmospheric CO₂ concentration. This is done through an ensemble of 100 000-year-long simulations, each driven by idealized CO₂ emission pulses. Our findings indicate that depending on the magnitude of the emission, 75 % of anthropogenic CO₂ is removed within 197–1820 years of the peak CO₂ concentration (with larger cumulative emissions taking longer to remove). Approximately 4.3 % of anthropogenic CO₂ will remain beyond 100 kyr. We find that the uptake of carbon by land, which has only been considered to a small extent in previous studies, has a significant long-term effect, storing approximately 4 %–13 % of anthropogenic carbon by the end of the simulation. Higher-emission scenarios fall on the lower end of this range, as increased soil respiration leads to greater carbon loss. For the first time, we have quantified the effect of dynamically changing methane concentrations on the long-term carbon cycle, showing that its effects are likely negligible over long timescales. The timescale of carbon removal via silicate weathering is also reassessed here, providing an estimate (80–105 kyr) that is significantly shorter than some previous studies due to higher climate sensitivity, stronger weathering feedbacks, and the use of a spatially explicit weathering scheme, leading to faster removal of anthropogenic CO₂ in the long term. Furthermore, this timescale is shown to have a non-monotonic relationship with cumulative emissions. Our study highlights the importance of adding model complexity to the global carbon cycle in Earth system models to accurately represent the long-term future evolution of atmospheric CO₂.

1 Introduction

A large amount of research has been dedicated to studying the impact of anthropogenic CO₂ emissions on the climate, what oceanographer Roger Revelle stated was “man's great geophysical experiment” (Revelle, 1956; Revelle et al., 1965). The scope of such research is usually limited to centennial timescales due to their relevance for governance and policy. However, there is a growing societal need for and an increasing number of scientific inquiries into the long-term (> 10³ years) future. This concerns the site selection for nuclear waste disposal and post-closure safety assessments of deep geological repositories in particular, as a number of environmental factors, such as subterranean stress, permafrost, erosion, and subrosion, can potentially compromise long-term safety (Näslund et al., 2013; Lord et al., 2015a; Lindborg et al., 2018; Turner et al., 2023; Kurgyis et al., 2024).

Previous studies on long-term future climate evolution have primarily focused on one of two tightly related research areas. Some, for example, have investigated the length of the current interglacial period under natural conditions: Berger and Loutre (2002) showed that we were already positioned to experience an unusually prolonged interglacial period as we currently approach a minima in the 100 and 400 kyr eccentricity cycle. A second area of interest is the lifetime of anthropogenic CO₂ in the atmosphere¹: Archer et al. (1997) and Archer et al. (1998) improved previous estimates by including CaCO₃ dissolution kinetics and showed that the return of the atmospheric CO₂ concentration to pre-industrial conditions would take tens of thousands of years. These two overarching areas of focus converged in Archer (2005) and Archer and Ganopolski (2005), where they showed that (1) not only is anthropogenic CO₂ expected to survive in our atmosphere for several hundreds of thousands of years but (2) it could also significantly delay the next glacial period. This was later confirmed by other studies (Ganopolski et al., 2016; Talento and Ganopolski, 2021).

In the past, general circulation models (GCMs) and intermediate complexity models (EMICs) have been used to study the 1000–10 000 year carbon cycle response to anthropogenic CO₂ emissions (see Archer et al., 2009a; Joos et al., 2013). Modelling studies become more scarce when going even further forward in time due to the rising computational cost, and as such, experiments are based on either EMICs (Ridgwell and Hargreaves, 2007; Charbit et al., 2008; Eby et al., 2009; Shaffer, 2010; Goodwin and Ridgwell, 2010; Meissner et al., 2012; Ganopolski et al., 2016; Brault et al., 2017; Jeltsch-Thömmes and Joos, 2020; Duque-Villegas et al., 2022) or models of lower complexity (Sundquist, 1991; Walker and Kasting, 1992; Lenton and Britton, 2006; Tyrrell et al., 2007a; Uchikawa and Zeebe, 2008; Cawley, 2011; Herrero et al., 2013; Köhler, 2020; Talento and Ganopolski, 2021; Couplet et al., 2025). The results of these studies largely agree in some aspects. For example, it has been established that the magnitude of cumulative emissions rather than the emission pathway will predominantly control the long-term response (Eby et al., 2009; Matthews et al., 2009; Zickfeld et al., 2012; Herrington and Zickfeld, 2014). Contemporary predictions for deep-future CO₂ concentrations exhibit considerable diversity, however. For illustration, Fig. 1 provides a graphical overview of the current estimates of atmospheric CO₂ concentration 10 000 years after the present for different cumulative CO₂ emissions.

Figure 1 Estimates on the impact of different cumulative emissions on atmospheric CO₂ concentration 10 kyr after the start of the simulation. The data from previous studies shown here were acquired through visual inspection of graphs, meaning that the figure is semi-qualitative and mostly serves an illustrative purpose (small errors may be present). Experiments shorter than 10 000 years or those that considered emissions of 6000 Pg C or larger (i.e., unconventional fossil fuel resources) were omitted. Different experiments were chosen from these publications to highlight the fact that the large diversity is primarily due to how long-term carbon cycle processes are resolved. The reference experiment in this study (yellow) has been plotted for comparison. Further information can be found in the respective publications. ^† Data for Sundquist (1991) taken from Tyrrell et al. (2007a). ^‡ Data generated from the emulator based on cGENIE results. ${}^{**}$ Data calculated from the fraction of emissions remaining. ^‖ Best solution from the model ensemble. ^* Data calculated from the formula provided in the appendix of Duque-Villegas et al. (2022).

Quantifying the atmospheric lifetime of anthropogenic CO₂ is challenging because it is controlled by a number of biochemical and geological mechanisms operating on different timescales up to hundreds of thousands of years (Archer et al., 2009a). Modelling studies suggest that the terrestrial biosphere will initially take up anthropogenic CO₂ (owing to enhanced vegetation productivity from CO₂ fertilization), although increased cellular respiration, soil carbon decomposition, land use change, and wildfires could eventually turn the land into a carbon source in a warming world. Extensive fossil fuel emissions can therefore overpower the capacity of conventional land carbon reservoirs to absorb anthropogenic CO₂, and only the ocean remains as the CO₂ sink on long timescales through processes such as air–sea CO₂ exchange, ocean invasion, seafloor CaCO₃ reactions, and carbonate accumulation. The substantial influence that marine biogeochemistry has on the long-term carbon cycle has also been recognized by many studies, in particular regarding the reaction between dissolved anthropogenic CO₂ and calcium carbonate in deep-ocean sediments (Archer et al., 1997; Middelburg et al., 2020). On even longer timescales, carbon exchange with geological reservoirs becomes important through processes such as sediment burial, chemical weathering of rocks on land, and volcanic degassing. As a consequence of all of these processes, carbon cycle models predict a highly nonlinear removal of anthropogenic CO₂ over time. The decline in anthropogenic CO₂ is usually presented as a superposition of exponential decays (Maier-Reimer and Hasselmann, 1987; Archer et al., 1997; Archer and Brovkin, 2008; Colbourn et al., 2015; Lord et al., 2015b), with each function representing a different process in the carbon cycle that takes up carbon. This generally produces a long tail of anthropogenic CO₂ concentration, which has been shown to persist for hundreds of thousands of years.

Due to the complexity and variety of poorly understood processes, the long-term future climate evolution remains highly uncertain even when not considering the unpredictability of future anthropogenic emissions (which can be reconciled to some extent using an ensemble of emission scenarios). Given this complexity, it becomes essential to make assumptions based on the timescales involved, the specific research question being addressed, and the resolution of the model employed. Previous studies have often made significant simplifications to (e.g., Lenton and Britton, 2006) or completely neglected (e.g., Archer et al., 1997; Lord et al., 2015b; Köhler, 2020) certain aspects related to land carbon. Weathering processes have also been frequently simplified (e.g., Charbit et al., 2008; Uchikawa and Zeebe, 2008) or their feedbacks not considered (e.g., Montenegro et al., 2007; Ridgwell and Hargreaves, 2007; Eby et al., 2009). As computational expense is the biggest factor inhibiting high-resolution modelling on longer timescales, these studies also did not consider the interaction between CO₂ and ice sheets, methane deposits, and the effect of future glacial cycles (Archer, 2005).

Here, we provide a set of idealized long-term transient climate change scenarios for the next 100 000 years using the fast Earth system model CLIMBER-X with a comprehensive carbon cycle model. These experiments seek to quantify the degree of uncertainty in the anthropogenic factor by employing a wide range of idealized emission scenarios to evaluate the spread in possible climate responses. Considering the ongoing and substantial disagreement regarding future CO₂ evolution among different models, we provide an analysis of the sensitivity of the model response to climate sensitivity and different carbon cycle processes. This is a first step towards a fully coupled simulation of the next 100 kyr with CLIMBER-X, which includes interactive ice sheets. The paper is structured as follows: we first describe our experimental setup, model, and considerations (Sect. 2). The results are then contextualized in relation to other studies by analyzing the model response, atmospheric lifetime, and subsequent removal timescales of different carbon cycle processes (Sect. 3). This is followed by a corresponding sensitivity study (Sect. 4). We conclude with a brief discussion, summary of our findings, and an outlook on further investigations (Sect. 5).

2 Methods

2.1 Model description

We use CLIMBER-X v1.1 (Willeit et al., 2022, 2023), a fast Earth system model designed to simulate the climate evolution on timescales ranging from decades to > 100 000 years. CLIMBER-X includes the statistical–dynamical semi-empirical atmosphere model SESAM (Willeit et al., 2022), the 3D frictional–geostrophic ocean model GOLDSTEIN (Edwards et al., 1998; Edwards and Marsh, 2005), the dynamic–thermodynamic sea ice model SISIM (Willeit et al., 2022), the land surface and dynamic vegetation model PALADYN (Willeit and Ganopolski, 2016), and the ocean biogeochemistry and marine sediment model HAMOCC (Heinze and Maier-Reimer, 1999; Ilyina et al., 2013; Mauritsen et al., 2019). CLIMBER-X operates at a horizontal resolution of 5°×5°. It is computationally efficient and can simulate up to 10 000 years per day, allowing it to perform a large ensemble of long-term simulations. This is also possible because short-term processes (e.g., weather, diurnal cycle) are not resolved. CLIMBER-X can generally represent historical changes in climate (Willeit et al., 2022) and the carbon cycle well (Willeit et al., 2023) and is suited to provide credible simulations for very long timescales into the deep future. Some studies using CLIMBER-X to transiently simulate the next several thousand years with a focus on the stability of the Greenland ice sheet have already been performed (e.g., Höning et al., 2023, 2024).

2.2 Open carbon cycle

The CLIMBER-X carbon cycle model has already been described in detail (Willeit et al., 2023) and has been shown to effectively represent the cycling of carbon through the atmosphere, vegetation, soils, permafrost, seawater, and marine sediments. For our experiments, we run the simulations using the open carbon cycle setup (Willeit et al., 2023), with interactive marine sediments, weathering on land, and the resulting geological sources and sinks of carbon. Carbon is not conserved in this setup; it is removed from the system through sediment burial and introduced into the system via weathering, volcanic outgassing, and anthropogenic emissions. In the open carbon cycle setup, a simplification is made to ensure that the budgets for silicate and phosphorus within the ocean–sediment system are balanced. This was done by disabling riverine fluxes of phosphorus and silicate and by assuming that organic carbon (which includes P) and opal (which includes Si) that are buried in the sediments (and, therefore, removed from the system) are instead returned to the surface ocean in remineralized form. The conservation of these inventories in the ocean–sediment system removes challenges related to nutrient conservation that would otherwise complicate the analysis and the interpretation of model results. The globally uniform atmospheric CO₂ concentration in CLIMBER-X is computed interactively using source and sink terms in the following prognostic equation for the total carbon content stored as atmospheric CO₂ (C_atm in Pg C):

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

where F_lnd is the global net land-to-atmosphere carbon flux (Pg C yr⁻¹), F_ocn is the net sea–air carbon flux (Pg C yr⁻¹), F_anth are anthropogenic carbon emissions (Pg C yr⁻¹), F_volc is the volcanic outgassing flux (Pg C yr⁻¹), and F_weath is the CO₂ consumption by carbonate and silicate weathering (Pg C yr⁻¹) (Willeit et al., 2023). PALADYN includes a rock weathering scheme influenced by runoff and temperature (Hartmann, 2009a; Börker et al., 2020), accounting for 14 different lithologies, as described in Hartmann and Moosdorf (2012). The weathering module computes dissolved inorganic carbon (DIC) and alkalinity fluxes to the ocean based on the release of bicarbonate ions (${\mathrm{HCO}}_{\mathrm{3}}^{-}$) into rivers. Spatially explicit carbonate and silicate weathering fluxes take the form

$$\begin{array}{}\text{(2)}& {\displaystyle }\begin{array}{rl}& F_{\text{weath}}^{\text{carb}}=F_{\text{weath,0}}^{\text{carb}}\\ & \sum _{\mathrm{\ell }=\mathrm{1}}^{\mathrm{13}}\left\{\begin{array}{l}\mathit{\beta }\left(\mathrm{\ell }\right)\times \left(\mathrm{3.1692}\times {\mathrm{10}}^{-\mathrm{3}}\right)R_{\mathrm{off}}\\ \text{for}\mathrm{\ell }=\text{carbonate sedimentary rock},\\ \mathit{\beta }\left(\mathrm{\ell }\right)\times {\mathrm{10}}^{-\mathrm{6}}{\mathrm{10}}^{\mathrm{9}{\left(\mathrm{1}+e^{-\mathrm{0.63}\left({\log }_{\mathrm{10}}{R}_{\mathrm{off}}+\mathrm{0.76}\right)}\right)}^{-\mathrm{1}}-\mathrm{2}}\\ \text{for}\mathrm{\ell }=\text{loess,}\\ \mathit{\beta }\left(\mathrm{\ell }\right)\mathit{\alpha }\left(\mathrm{\ell }\right)\times b\left(\mathrm{\ell }\right)\times e^{\frac{\mathrm{1000}{E}_{\mathrm{a},\mathrm{carb}}}{R(T\cdot T_{\mathrm{0}})}\left(T-T_{\mathrm{0}}\right)}{R}_{\mathrm{off}}\\ \text{for}\mathrm{\ell }=\text{other lithologies (Table E1)},\end{array}\right.\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}{F}_{\text{weath}}^{\text{sil}}& =F_{\text{weath,0}}^{\text{sil}}\sum _{\mathrm{\ell }=\mathrm{1}}^{\mathrm{13}}\mathit{\beta }\left(\mathrm{\ell }\right)(\mathrm{1}-\mathit{\alpha }(\mathrm{\ell }\left)\right)\times b\left(\mathrm{\ell }\right)\\ & \times e^{\frac{\mathrm{1000}{E}_{\mathrm{a},\mathrm{sil}}\left(\mathrm{\ell }\right)}{R(T\cdot T_{\mathrm{0}})}\left(T-T_{\mathrm{0}}\right)}{R}_{\mathrm{off}}.\end{array}\end{array}$$

Weathering fluxes are largely described by the Arrhenius equation, with the exception of carbonate weathering for carbonate sedimentary rocks and loess, which is instead described by the runoff-dependent weathering equations from Amiotte Suchet and Probst (1995) and Börker et al. (2020), respectively. In Eqs. (2) and (3), ℓ represents the different lithologies, T is the annual mean near-surface air temperature (K), T₀ is 284.2 K, R_off is the annual runoff ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), b(ℓ) is the molality/weathering rate (mol C kg⁻¹ water), E_a,carb is the activation energy of carbonates (14 kJ mol⁻¹), E_a,sil(ℓ) is the activation energy of silicates (kJ mol⁻¹), R is the molar gas constant (8.3145 $\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), α(ℓ) is the fraction to weather as carbonate rocks, and β(ℓ) is the fraction of a given lithology in a grid cell (Fig. E1). The values for α, b, and E_a are all based on the calibrated runoff- and temperature-dependent models from Hartmann (2009a) and Hartmann et al. (2014) and are displayed in Table E1.

Although the Earth experiences a number of external forces, this study generally excludes factors like impact events, changes in tectonic configuration, and solar luminosity, as they are either unpredictable or operate on such long timescales that they are not relevant for our investigation. Similarly, we do not resolve the deep-carbon cycle (i.e., carbon subduction and recycling in the Earth's mantle), as its effects only become important on timescales longer than our experiment duration (i.e., 10⁶–10⁷ years) and involve processes that cannot be resolved with Earth system models. As such, we assume that the pre-industrial carbon cycle was in equilibrium, which implies that the constant volcanic outgassing is prescribed and is set to half the global silicate weathering flux at the pre-industrial time (0.0738 Pg C yr⁻¹ or 6.15 Tmol C yr⁻¹). This condition ensures that the atmospheric CO₂ is in equilibrium under pre-industrial conditions (Munhoven and François, 1994; Willeit et al., 2023), which has been a common assumption in previous investigations that looked at the long-term forced response of the climate (e.g., Colbourn et al., 2015; Lord et al., 2015b).

2.3 Experimental setup

Model simulations are started from a pre-industrial equilibrium state that has been obtained from a 100 000 year equilibrium spin-up of the carbon cycle model, as described in Willeit et al. (2023). The pre-industrial CO₂ concentration is taken to be 280 ppm, whereas the pre-industrial orbital parameters of eccentricity e is 0.0167, obliquity ϵ is 23.46°, and the perihelion ω is 100.33°. This corresponds to a maximum summer insolation at 65° N of 480.4 W m⁻². In the reference (REF) experiment, we set a constant methane concentration of 600 ppb, which is representative of Holocene conditions (Sapart et al., 2012; Mitchell et al., 2013; Beck et al., 2018). Although CLIMBER-X is capable of simulating ice sheets (e.g., Höning et al., 2023; Talento et al., 2024; Willeit et al., 2024), we conduct our simulations with prescribed present-day Greenland and Antarctica to ensure comparability with earlier studies. Interactive ice sheets will be enabled in a follow-up study. To simplify interpretation and ensure consistency with previous studies, orbital forcing was fixed at present-day values, with the combined effects of anthropogenic and orbital forcing to be explored in a future study. All simulations were run for 100 000 years, with constant orbital parameters and without any climate acceleration technique.

Table 1 Overview of the experimental configuration for the simulations performed in this study. For each experimental configuration, all emission forcings (0–5000 Pg C) were applied to produce an ensemble. Emissions pathways are shown in Fig. E2.

To investigate the lifetime of anthropogenic CO₂ and the effect it has on the long-term evolution of the climate, we introduce idealized CO₂ emission scenarios representing cumulative emissions from 0 to 5000 Pg C (Fig. E2) to cover a variety of different anthropogenic scenarios. The upper limit of ∼ 5000 Pg C has been used in many previous studies (e.g., Archer, 2005; Montenegro et al., 2007; Uchikawa and Zeebe, 2008; Archer et al., 2009a), as it just surpasses the current estimated maximum of conventional fossil fuel reserves (Lal, 2008; McGlade and Ekins, 2015). To simulate a variety of emissions pathways and explore the sensitivity of the results to the duration of the CO₂ emission scenarios, we designed three different Gaussian functions with different shapes (with an increasing mean and standard deviation with higher cumulative CO₂; Fig. E2) and one pulse-like perturbation where all emissions are released within the first year of the simulation. It should be noted here that the 100 kyr simulation duration includes the ramp-up and ramp-down of emissions. Following the emission function, emissions are set to 0 Pg C yr⁻¹ for the rest of the simulation.

A number of additional experiments were performed to isolate the effect of different carbon cycle processes on the long-term fate of anthropogenic CO₂. Specifically, we explored the role of the land carbon cycle response (noLAND), the impact of enabling interactive methane concentrations in the atmosphere (intCH4), and the impact of the emissions pathway (PULSE, PATH1, PATH2) on the lifetime of anthropogenic CO₂. To isolate the effect of the weathering feedback on the CO₂ lifetime, a dedicated set of experiments (noWEATH) was conducted in which weathering fluxes were forced to remain constant at their pre-industrial rate throughout the simulation (0.2376 Pg C yr⁻¹ or 19.8 Tmol C yr⁻¹ for carbonate weathering and 0.1476 Pg C yr⁻¹ or 12.3 Tmol C yr⁻¹ for silicate weathering). The sensitivity of the results to different equilibrium climate sensitivities (ECSs) between 2 and 4 °C was additionally investigated (ECS2, ECS4) by rescaling the equivalent CO₂ in the long-wave radiation scheme (further described in Appendix A). A summary of all experimental configurations is presented in Table 1.

Figure 2 Responses in (a) atmospheric CO₂ concentration, (b) global surface air temperature, (c) land carbon inventory, (d) ocean carbon inventory, and (e) sediment carbon inventory to the full ensemble of emission scenarios (0–5000 Pg C) in the REF experiment for 100 kyr.

3 Results

3.1 Long-term carbon cycle evolution

An overview of the 100 kyr atmospheric CO₂ concentration and global mean temperature responses across all emission scenarios in the reference (REF) experiment can be seen in Fig. 2. The initialization of our experimental setup is demonstrably robust, evidenced by the less than 2 ppm drift (from the pre-industrial value of 280 ppm) of atmospheric CO₂ concentration in the zero-emissions scenario over the span of 100 kyr (Fig. 2a). All simulations with anthropogenic emissions show a rapid increase in atmospheric CO₂ concentration during the ramp-up of emissions, with peak concentrations occurring within the first 200 years. In the REF ensemble, peak CO₂ concentrations fall between ∼ 360 and 1820 ppm across the different emission scenarios (Fig. 2a). These values are within the range of magnitudes reported by previous studies (e.g., Montenegro et al., 2007; Colbourn et al., 2015). We find that the increase in global mean surface temperature lags behind that of the CO₂ concentration (Fig. 2b) and that the size of the temporal lag tends to increase as the cumulative emission grows (Ricke and Caldeira, 2014; Zickfeld and Herrington, 2015). Peak temperature anomalies range between ∼ 0.6 and 6.4 °C across the different emission scenarios in the REF ensemble. However, the extended decline in the Atlantic Meridional Overturning Circulation (AMOC; Fig. 7e) in the 5000 Pg C scenario results in a cooling effect in the Northern Hemisphere that prevents global mean temperature from further rising, thereby truncating peak warming. After emissions cease (∼ 150–200 years after the start of the simulation; Fig. E2), both atmospheric CO₂ concentration and global temperatures decrease as carbon is taken up by the land and ocean pools (Fig. 2a–d). The only exception to this is the 5000 Pg C scenario, where there is a global temperature plateau, and global cooling is delayed until approximately ∼ 700 years after the start of the simulation, coinciding with a substantial recovery of the AMOC (Figs. 2b, 7e). Temperatures slowly decrease towards pre-industrial levels after the first millennium, and, after 10 kyr, global warming is reduced by more than half (Fig. 2b). After 100 kyr, atmospheric CO₂ concentration still remains ∼ 10–100 ppm above pre-industrial levels across the different emission scenarios, and temperatures are still between ∼ 0.1 and 1.4 °C higher than the pre-industrial level at the end of the simulation.

Figure 3 The relative partition of anthropogenic carbon in the REF experiment after peak atmospheric CO₂ concentration. The partition of the atmosphere (grey), land (green), ocean and sediments (blue), and silicate weathering (teal) is shown for the (a) 500, (b) 1000, and (c) 5000 Pg C scenarios. The calculations for these percentages is outlined in Appendix B using cumulative fluxes from the atmosphere. The carbon fraction attributed to the ocean and sediments is labelled this way, as the air–sea CO₂ flux may indirectly include carbon from sediments through dissolution. However, we do not account for sediment carbon fluxes directly since there is no direct air–sediment flux (see Appendix B for further details on this formulation). It should be noted that the ramp-up of emissions is excluded here.

The fate of the anthropogenic CO₂ emissions (and how this carbon is distributed among the different Earth system components) after the peak atmospheric CO₂ concentration is shown in Fig. 3. The land and ocean both rapidly absorb CO₂ during the ramp-up of emissions through enhanced vegetation productivity, air–sea CO₂ exchange, and its subsequent dissolution (Fig. 4). At peak CO₂ concentrations, the land and ocean are responsible for between 34 % and 61 % of carbon removal (19 %–39 % by land, 15 %–22 % by the ocean; Fig. 3), with the carbon removal fractions decreasing with increasing cumulative emissions. During the first 1 kyr after peak CO₂ concentrations, oceanic CO₂ uptake gradually slows as surface waters approach equilibrium with the atmosphere (Fig. 4b). Although the ocean carbon fraction largely increases over the first 1 kyr after peak CO₂ concentrations (accounting for 46 %–59 % 1 kyr after emissions cease; Fig. 3), the carbon fraction taken up by the ocean during the latter half of the millennium appears disproportionately larger than what is taken up during the first 500 years, as land carbon has already begun decreasing by this time (Figs. 2c, 4a). Beyond the first millennium, geological processes (e.g., marine sediments and weathering) largely control the evolution of atmospheric CO₂ concentration (Archer and Brovkin, 2008). Silicate weathering is responsible for 11 %–17 % of carbon removal 10 kyr after peak CO₂ concentrations. However, by 100 kyr, it is responsible for the majority of carbon removal (62 %–78 %) and becomes the dominant removal process on the timescale of hundreds of thousands of years.

Figure 4 Changes in (a) net atmosphere-to-land carbon flux, (b) net air–sea carbon flux, and (c) weathering consumption of CO₂ in different emission scenarios over 100 kyr. Colours correspond to the cumulative emission scenarios shown in Fig. 2.

3.1.1 Land carbon response

Two opposing mechanisms involving vegetation and soil carbon predominantly govern the land carbon response to climate change. The land carbon pool (as the sum of vegetation and soil carbon) is initially a carbon sink due to increases in vegetation carbon (Brovkin et al., 2013), as CO₂-driven fertilization enhances net primary productivity. These increases, however, can be partially or totally offset by warming-enhanced soil respiration, which strongly depends on model parameters; on whether permafrost, peatlands, and wetlands are included; and on how they are modelled in the carbon cycle model (Zickfeld et al., 2013; Eby et al., 2013). In our simulations, the peak increase in vegetation carbon significantly exceeds and more than compensates for the maximum loss in soil carbon (Fig. 5). Peak land carbon uptake in our simulations falls comfortably between the values given by previous studies (Mikolajewicz et al., 2006; Vakilifard et al., 2022).

Figure 5 Changes in (a) vegetation and (b) soil carbon inventories relative to the pre-industrial period in different emission scenarios over 100 kyr. Colours correspond to the cumulative emission scenarios shown in Fig. 2.

Figure 6 Mean annual (a–c) net primary productivity, (d–f) vegetation carbon, and (g–i) soil carbon in different emission scenarios at 1 kyr. The columns show (a, d, g) the absolute values in the zero-emissions scenario (equivalent to the pre-industrial values), (b, e, h) the 1000 Pg C scenario relative to the pre-industrial values, and (c, f, i) the 5000 Pg C scenario relative to the pre-industrial values. Absolute values for the zonal mean are additionally shown for (j) NPP, (k) vegetation, and (l) soil carbon for the zero-emissions (black), 1000 Pg C (purple), and 5000 Pg C (red) scenarios. Magenta lines (g–i) display the permafrost area.

In the latter half of the first millennium, vegetation carbon decreases as atmospheric CO₂ concentrations fall (Fig. 5a). Soil carbon exhibits non-monotonic behaviour with increasing emissions: low-emission scenarios (< 1000 Pg C) have gained up to 100 Pg C, and high-emission scenarios (> 4000 Pg C) have lost up to 150 Pg C by 1 kyr (Fig. 5b). This is because in the low-emission scenarios, the accumulation of carbon in the soil exceeds the carbon loss due to enhanced soil respiration (Fig. 6h), preventing the overall decline in soil carbon inventory (Fig. 5b). This is not the case in high-emission scenarios, as permafrost area decreases so significantly that soil respiration counteracts any potential increase when temperatures are sufficiently high (Fig. 6i). Most of the globe sees an increase in vegetation carbon with increasing emissions, mainly because of increases in net primary productivity (NPP; Fig. 6b, c, e, f). However, there are a few regions that see a decline in vegetation carbon due to a reduction in precipitation (Fig. 11h, i). Although the soil carbon inventory decreases in high-emission scenarios by 1 kyr, it recovers over the next 100 kyr as temperatures gradually decrease. In our simulations, the terrestrial storage of anthropogenic carbon is positive during the entire 100 kyr due to an increase in vegetation carbon and the recovery of soil carbon. At the end of the experiments, the land still stores approximately 7 % (range of 4 %–13 %) of the emitted anthropogenic carbon (Fig. 3).

Figure 7 Changes in (a) ocean buffering capacity (Revelle factor), (b) global mean surface pH, (c) relative changes in surface DIC concentration and global ocean DIC inventory, (d) relative changes in surface alkalinity concentration and global alkalinity inventory, and (e) strength of the Atlantic Meridional Overturning Circulation (AMOC) in different emission scenarios over 100 kyr. Colours correspond to the cumulative emission scenarios shown in Fig. 2. The surface ocean Revelle factor shown in (a) has been calculated via Eq. (C1) in Appendix C.

Figure 8 Changes in (a) marine net primary production, (b) CaCO₃ production, (c) sediment particulate organic carbon inventory, (d) sediment CaCO₃ inventory, and (e) sediment CaCO₃ dissolution rate in different emission scenarios over 100 kyr. Colours correspond to the cumulative emission scenarios shown in Fig. 2.

3.1.2 Ocean and sediment carbon response

The ocean absorbs and stores atmospheric CO₂ due to a sequence of chemical, biological, physical, and geological processes. Under anthropogenic emissions, the ocean initially takes up carbon at an accelerated rate from greater air–sea CO₂ exchange (Fig. 4b). At peak CO₂ concentrations, the ocean has a net cumulative uptake between ∼ 100 and 780 Pg C depending on the emission scenario (Fig. 2d). Oceanic CO₂ uptake yields bicarbonate and hydrogen ions from CO₂ dissolution. As a consequence, surface ocean pH decreases at a similar rate to that which carbon is taken up in the ocean, with larger pH decreases associated with higher cumulative emissions (Fig. 7b). The magnitude of the surface ocean pH decrease in our experiments matches well with that of other studies: for example, the peak pH anomaly of approximately −0.7 in the 5000 Pg C experiment falls within the range of −0.6 to −0.8 given by other models (Caldeira and Wickett, 2005; Montenegro et al., 2007; Uchikawa and Zeebe, 2008; Eby et al., 2009; Zickfeld et al., 2012; Colbourn et al., 2015). As the phytoplankton growth rate strongly depends on temperature in HAMOCC (Ilyina et al., 2013; Willeit et al., 2023), NPP generally increases with higher sea surface temperature (SST) under larger cumulative emissions (Fig. 8a). Warming also accelerates the remineralization of organic matter in the ocean, which increases nutrient availability in the surface ocean for primary production. Moreover, reduced sea ice cover at high latitudes increases light availability and, therefore, enhances local productivity. In our highest-emission scenario, the combined effects amount to an approximate ∼ 20 % increase in peak NPP under the 5000 Pg C scenario compared to the pre-industrial value (53.3 Pg C yr⁻¹; Fig. 8a).

As the amount of dissolved CO₂ in the surface ocean rises, its ability to absorb excess CO₂ decreases. This decrease in buffering capacity is illustrated by the Revelle factor (ξ, the sensitivity of pCO₂ to changes in DIC; Revelle and Suess, 1957), which increases with higher emissions (Fig. 7a). The pre-industrial Revelle factor of ∼ 9.4 matches reasonably well with that calculated by Jiang et al. (2019), and the peak Revelle factor in the 500 Pg C scenario (∼ 10.2) is close to that calculated using present-day GLODAPv2 observations (Terhaar et al., 2022). Between 100 years and 1 kyr, ξ slowly decreases owing to the gradual decrease in surface DIC (except for the 5000 Pg C scenario, Fig. 7c) and the increase in surface alkalinity (Fig. 7d). The decrease in surface DIC results from the combined effect of the downward transport of anthropogenic carbon in the surface ocean and the enhanced biological production (Fig. 8a), which is partly counteracted by the increase in the global DIC inventory (Fig. 7c). Both DIC and alkalinity inventories increase during this period due to high carbonate weathering fluxes (input as ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ into the ocean; see Fig. 9a) and enhanced refluxes from sediment. Although NPP increases with warming, the loss of organic carbon to the sediment (not shown) and the organic carbon content in the sediment (Fig. 8c) are reduced due to higher remineralization rates in both the water column and sediment. Deep convection, which transports surface waters with decreased pH and carbonate ion concentrations into the interior ocean, together with increased organic matter remineralization, causes carbonate ion concentrations in the interior ocean to fall with increasing emissions (Fig. E5). In response, the lysocline and carbonate compensation depth (CCD) shoal and sediment CaCO₃ dissolves (contributing more to the total loss of the sediment carbon inventory than sediment organic carbon does; Fig. 2e), driving the increase in the alkalinity inventory. In the 5000 Pg C scenario, surface DIC increases between the years 1100 and 1700, which is likely due to the extended period of the AMOC decline (Fig. 7e) that inhibits deep convection and downward carbon transport.

Surface DIC shows diverse variations between 1 and 10 kyr depending on the emission scenario, whereas surface alkalinity generally increases and is responsible for the recovery of the buffering capacity (shown by the continuous decline in ξ). The variation in surface alkalinity generally follows that of the global inventory, with the latter controlled by the weathering input and the net loss to the sediment. The sediment CaCO₃ continues to dissolve at a higher rate than in the first millennium (Fig. 8e). Beyond 10 kyr, both surface DIC and alkalinity decline following the global inventory as a response to decreased carbonate weathering (Fig. 9a) and increased loss to the sediment (Fig. 8c, d). Due to silicate weathering (Sect. 3.1.3), which reduces atmospheric CO₂, the ocean becomes a carbon source, contributing to the decrease in surface DIC and the recovery of the buffering capacity. The impact of carbonate compensation and seafloor CaCO₃ neutralization on the restoration of the buffering capacity on the millennial timescales shown in these simulations is in line with previous studies (Archer et al., 1997, 1998; Tyrrell et al., 2007a; Tyrrell, 2007b). The recovery of the sediment CaCO₃ inventory occurs roughly around the same time as that predicted by Lenton and Britton (2006). However, the timing of peak CaCO₃ dissolution in CLIMBER-X is earlier than in a CLIMBER-2 study by Montenegro et al. (2007), which may be due to increased biological production of CaCO₃ (as it is proportional to organic detritus production and, thus, NPP; Fig. 8a) leading to more CaCO₃ export to the interior ocean.

Figure 9 Changes in (a) carbonate and (b) silicate weathering fluxes in different emission scenarios over 100 kyr. Colours correspond to the cumulative emission scenarios shown in Fig. 2.

Figure 10 Relationship between (a) carbonate weathering, (b) silicate weathering, and temperature in previous studies and CLIMBER-X. Solid coloured lines corresponding to the cumulative emission scenarios seen in Fig. 2 represent the weathering response in the REF ensemble. The dashed black lines in (a) and (b) show an exponential fit of the trajectory of the weathering response. Calculated trajectories have additionally been plotted for the equations used by Sundquist (1991), Uchikawa and Zeebe (2008), Colbourn et al. (2015) (with runoff, R_off, and productivity, P, parametrized as a function of CO₂, as in Colbourn et al., 2013), and Lord et al. (2015b) specifically using CLIMBER-X variables for atmospheric CO₂ (CO₂) and land surface temperature (T_lnd). The fit of these trajectories has been plotted above using Numpy.Polyfit in Python. It should be noted that since this comparison is based on the simulated CLIMBER-X variables, it implicitly assumes an ECS of 3.1 °C.

Figure 11 Mean annual (a–c) carbonate weathering, (d–f) silicate weathering, and (g–i) precipitation in different emission scenarios at 1 kyr. The columns show (a, d, g) the absolute values in the zero-emissions scenario (equivalent to the pre-industrial values), (b, e, h) the 1000 Pg C scenario relative to the pre-industrial values, and (c, f, i) the 5000 Pg C scenario relative to the pre-industrial values. Absolute values for the zonal mean are additionally shown for (j) carbonate weathering, (k) silicate weathering, and (l) precipitation for the zero-emissions (black), 1000 Pg C (purple), and 5000 Pg C (red) scenarios.

3.1.3 Weathering response

While the pre-industrial weathering fluxes in our simulations (19.8 Tmol C yr⁻¹ for carbonate and 12.3 Tmol C yr⁻¹ for silicate; Fig. 9) fall within but toward the higher end of observational estimates (Meybeck, 1987; Gaillardet et al., 1999; Munhoven, 2002; Amiotte Suchet et al., 2003), they are significantly higher than previous studies investigating the long-term future CO₂ evolution (Uchikawa and Zeebe, 2008; Colbourn et al., 2015; Lord et al., 2015b). This is at least partly due to the spatially explicit weathering scheme employed in our model, as switching from a 0D to a 2D scheme has been shown to increase global weathering rates (Colbourn, 2011; Colbourn et al., 2013). Chemical weathering on land increases in high-emission scenarios due to warmer and wetter conditions, drawing down more carbon from the atmosphere (Fig. 9). Both carbonate and silicate weathering rapidly increase over the first 200 years as a response to elevated temperatures and increased runoff (Lenton and Britton, 2006). The carbonate weathering flux increases to 20.8–32.0 Tmol C yr⁻¹ in the different scenarios around the same time as the peak atmospheric CO₂ concentration (Fig. 9a). Peak silicate weathering occurs at the same time but has even larger differences from the pre-industrial value, nearly tripling under the high-emission scenario (Fig. 9b). After peak weathering fluxes are reached in the REF ensemble, they begin to decrease within the first millennium as atmospheric CO₂ decreases.

The sensitivity of carbonate and silicate weathering to climate change in CLIMBER-X is compared to estimates from earlier investigations in Fig. 10a, b. Generally, the CLIMBER-X response falls between those estimated by previous studies (Sundquist, 1991; Uchikawa and Zeebe, 2008; Colbourn et al., 2015; Lord et al., 2015b), with the response of carbonate weathering falling in the middle and silicate weathering tending to be on the stronger side. In both cases, however, the weathering response is higher than was simulated by the equations used in Lord et al. (2015b). In addition to this, the more sophisticated weathering scheme in CLIMBER-X produces a narrow hysteresis in the response to the change in global mean temperature (Fig. 10a, b), as weathering fluxes decrease more slowly than global mean temperature does. This hysteresis is primarily due to a lag between the evolution of precipitation and temperature. The combination of a rather large pre-industrial silicate weathering flux and high sensitivity of silicate weathering to climate in the model results in a weathering feedback that is stronger than in some previous studies (Lord et al., 2015b; Colbourn et al., 2015). This has implications for the lifetime of atmospheric CO₂ (Sect. 3.2).

Most previous studies on deep-future CO₂ concentrations did not spatially resolve weathering, instead relying on global mean values. The few studies that did use spatially explicit weathering employed either the GKWM (Bluth and Kump, 1991) or GEM-CO2 (Amiotte Suchet et al., 2003) lithological maps, which have quite coarse resolutions compared to the GLiM (Hartmann, 2009a) lithological map used by CLIMBER-X. The pre-industrial carbonate and silicate weathering distributions are shown in Fig. 11a, d. These distributions align quite closely with the calculated distribution of CO₂ consumption from Hartmann et al. (2009b). However, not all details can be reproduced, given the 5°×5° resolution of CLIMBER-X and biases in precipitation (Fig. 11g; see Willeit et al., 2022). Weathering rates rise with increasing temperature, particularly in regions where significant weathering is already taking place. Runoff, which indirectly depends on precipitation through soil infiltration and drainage (Willeit and Ganopolski, 2016), drives some of these changes (Figs. 11g, h, i, E6). As emissions grow, silicate weathering sees a significant increase in the tropical rain belt area due to higher precipitation. These include regions such as Southeast Asia, central Africa, and Brazil (Fig. 11k). Large changes in carbonate weathering are more globally distributed than silicate weathering (Fig. 11e, f). However, like silicate weathering, large changes in carbonate weathering (e.g., increases and decreases in carbonate weathering in Southeast Asia and central Asia compared to the pre-industrial values) can also be explained by increases and decreases in precipitation (Fig. 11b, c). Such spatial patterns in silicate and carbonate weathering under idealized 1000 and 5000 Pg C scenarios have also been seen by Brault et al. (2017). In our simulations, highly active weathering regions are dominant contributors to the global weathering flux. About 10 % of land area across the different emission scenarios is responsible for ∼ 13 %–14 % of the total carbonate weathering and ∼ 33 %–37 % of the total silicate weathering (Fig. 11b, c, e, f). Given that the simulated weathering rates for both carbonate and silicate rocks are quite spatially heterogeneous, simplifications used by previous studies to calculate weathering rates (e.g., weathering as a function of global mean land surface temperature) could be insufficient to describe the feedback between climate and weathering.

3.2 Atmospheric lifetime of anthropogenic CO₂

The atmospheric lifetime of anthropogenic CO₂ across the different scenarios is estimated using an impulse response function (IRF). This is a first-order approximation of the long-term effects of a short-term perturbation. It can be applied to changes in atmospheric CO₂ concentration (normalized by cumulative CO₂ emissions) such that it represents the fraction of emissions remaining (Maier-Reimer and Hasselmann, 1987; Joos et al., 2013; Jeltsch-Thömmes and Joos, 2020). In this formulation, it is functionally no different than the airborne fraction, but we distinguish between the two by assuming that the IRF is only valid for times after the peak CO₂ concentration (CO${}_{\mathrm{2}}\left(t_{\mathrm{0}}\right)={{\mathrm{CO}}_{\mathrm{2}}}^{\max }$). This can be written as the following equation:

$$\begin{array}{}\text{(4)}& {\mathrm{IRF}}_{{\mathrm{CO}}_{\mathrm{2}}}\left(t\right)={\displaystyle \frac{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}\left(t\right)\times \mathrm{2.124}}{E}},\end{array}$$

where E is cumulative emissions in Pg C, ΔCO₂(t) is the change in atmospheric CO₂ concentration at time t compared to the pre-industrial concentration (280 ppm), and 2.214 Pg C ppm⁻¹ is the conversion factor between concentration and mass (Friedlingstein et al., 2023). The IRF for the REF ensemble is shown in Fig. 12a, c. The ramp-up of emissions is not included in this analysis. Unlike previous studies (e.g., Archer et al., 2009a; Joos et al., 2013), we do not deploy an instantaneous pulse in the REF ensemble, and our fractions of emissions remaining do not start at 1, as carbon sequestration began during the ramp-up of emissions. While it is unlikely that the emissions pathway will be an important factor for the atmospheric lifetime on long timescales, the IRF for the PULSE, PATH1 and PATH2 ensembles has additionally been plotted in Fig. 12b to explicitly show the impact of the emissions pathway. Although the emissions pathway has a notable impact on the IRF within the first millennium, this effect becomes largely negligible after ∼ 3 kyr. This supports the long-term path independence found by Zickfeld et al. (2012) and Herrington and Zickfeld (2014), meaning that despite the potential complication from the ramp-up and ramp-down periods of emissions in the REF experiment, our results are fully comparable with those in the Long-Term Model Intercomparison Project (LTMIP) in Archer et al. (2009a).

Figure 12 Atmospheric lifetime of anthropogenic CO₂ described by IRF${}_{{\mathrm{CO}}_{\mathrm{2}}}$(t). Colours correspond to the cumulative emission scenarios shown in Fig. 2. (a) Fraction of emissions remaining in the atmosphere after peak CO₂ concentration for the different emission scenarios in the REF experiment. (b) Fraction of emissions remaining in the atmosphere over the first 10 kyr after peak CO₂ concentration. Here, the x axis is log-scaled and only shows the 1000 Pg C (purple) and 5000 Pg C (red) scenarios. The PULSE, PATH1, and PATH2 experiments have additionally been plotted to show the effect of the emissions pathway on the IRF. (c) Same as (a) but for the first 10 kyr after peak CO₂ concentration. The diamond markers and horizontal bars represent the mean time and spread in which the different emission scenarios reach 50 %, ≈ $\mathrm{1}/e$ % (mean lifetime), 25 %, 10 %, and 5 % of their original cumulative emissions.

Modelling studies have suggested that the majority (up to 80 %) of anthropogenic CO₂ is quickly taken up by fast carbon processes, leading to a mean lifetime of only a few centuries (Archer, 2005; Montenegro et al., 2007; Archer and Brovkin, 2008). In the REF ensemble, the mean lifetime is 390 (range of 0–955) years (Fig. 12c). This value depends on the cumulative emission load to the atmosphere (Fig. 12c) and, to some extent, the pathway of the emission (Fig. 12b). This explains why the lifetime determined by the timing of IRF${}_{{\mathrm{CO}}_{\mathrm{2}}}\left(t\right)=\mathrm{1}/e$ is close to zero in the 500 Pg C scenario, as the majority of the carbon is taken up from the atmosphere prior to reaching the peak atmospheric CO₂ concentration in our simulations (Fig. 3c). On average, it takes between 197 and 1820 years after peak CO₂ concentration to remove 75 % of anthropogenic emissions, depending on the emission scenario. Our results (Fig. 12c) generally agree with the most comprehensive assessment to date, which indicates that between 20 % and 35 % of CO₂ remains in the atmosphere after 200–2000 years (Archer et al., 2009a). After 10 kyr, it has been previously reported that CO₂ emissions in the atmosphere would remain somewhere between ∼ 5 % and 20 % for the 1000 Pg C scenario and ∼ 10 % and 35 % for 5000 Pg C scenario (Archer et al., 2009a). In all experiments, the atmospheric CO₂ concentration has decreased to approximately 10 % of its initial value after 10 kyr across the different emission scenarios (Fig. 12c). Our values are on the lower end of the range provided by Archer et al. (2009a) due to the combined effect of a higher pre-industrial weathering rate and stronger weathering feedbacks. After 100 kyr from the start of the REF experiments, the fraction of emissions remaining is on average 4.3 % (range of 4.0 %–4.4 %) across the different emission scenarios (Fig. 12a). This falls in between estimates by previous studies, which have reported that between 2 % and 7 % of emissions will remain in the atmosphere beyond 100 kyr (Archer, 2005; Lenton and Britton, 2006; Jeltsch-Thömmes and Joos, 2020).

3.3 Multi-exponential timescale analysis

We conducted a least-squares fit of a multi-exponential decay function to tease out different timescales that can otherwise be associated with long-term processes that take up excess CO₂. This analysis was exclusively done on the REF ensemble, the ensemble with a pulse-like perturbation of CO₂ (PULSE; Table 1), and the ensemble with land carbon disabled (noLAND; Table 1). The latter two experiments were included to (1) provide confidence in our estimated timescales by excluding the ramp-up period of emissions as a factor (aligning with the procedure used in similar studies, e.g., Archer et al., 1997) and (2) facilitate a direct comparison with the findings of Lord et al. (2015b). The Python function curve_fit in Scipy.Optimize was used, employing nonlinear least squares to fit the following function to our data:

$$\begin{array}{}\text{(5)}& {\mathrm{CO}}_{\mathrm{2}}\left(t\right)=\mathrm{280}+\mathrm{\Delta }{{\mathrm{CO}}_{\mathrm{2}}}^{\max }\sum _{i=\mathrm{1}}^n{A}_i\cdot \exp \left(-{\displaystyle \frac{t}{{\mathit{\tau }}_i}}\right),\end{array}$$

where CO₂(t) is the CO₂ concentration at time t (where CO${}_{\mathrm{2}}\left(t_{\mathrm{0}}\right)={{\mathrm{CO}}_{\mathrm{2}}}^{\max }$), and $\mathrm{\Delta }{{\mathrm{CO}}_{\mathrm{2}}}^{\max }$ is the peak CO₂ concentration compared to the pre-industrial concentration (280 ppm). This function represents the superposition of exponential decay functions each with their own time constant, τ_i, and amplitude, A_i (Maier-Reimer and Hasselmann, 1987). Here, the former represents the timescale of decay (i.e., the removal timescale) and the latter represents the fraction of emissions removed over that timescale. Since an instantaneous emission pulse was not used in the REF and noLAND experiments, the first ∼ 200 years of our simulations were dedicated to the ramp-up and ramp-down of atmospheric CO₂ concentrations (Fig. E2). Therefore, we decided to fit all data starting from the peak concentration in atmospheric CO₂ (as done in Sect. 3.2), thereby omitting the ramp-up period of emissions. Several superimposed exponential functions were used to fit the data (from n=1 to 6), but fast convergence times and significantly high (> 0.97) R² values were found for n=3.

Figure 13 Dependence of fitting parameters (A_i and τ_i) determined in the multi-exponential decay analysis of CO₂ concentration in the REF, PULSE, and noLAND experiments. Different markers have been used to convey different experiments. The arrows in (a) indicate how the ability of individual processes to remove anthropogenic CO₂ changes as total emissions increase. This is shown more explicitly in (b) and (c), which depict how the fraction of emissions remaining and the removal timescales for different processes change with increasing emissions. Scenarios with smaller cumulative emissions typically result in lower A_i values, as a larger proportion of emissions is taken up by short-term processes (sub-centennial timescales), which are not considered in our multi-exponential fitting procedure (see Lord et al., 2015b).

We should first highlight the fact that it is not surprising that n=3 yielded a good fit across the different experiments. Although the equilibrium timescales for fast carbon processes have been thoroughly examined in the literature previously (e.g., Eby et al., 2009), Archer and Brovkin (2008) showed that long-term carbon uptake is predominately controlled by the slow processes of (1) ocean invasion, (2) reactions with CaCO₃, and (3) reactions with igneous (i.e., silicate) rocks. In some studies, these are sometimes more broadly grouped into ocean and geological processes, with CaCO₃ reactions further split into sea floor and terrestrial CaCO₃ neutralization (corresponding to n=4, e.g., Archer et al., 1997). Other studies have used higher orders of n to account for (or further break down) processes related to land carbon, air–sea CO₂ exchange, ocean invasion, carbonate chemistry, mixed-layer mixing, and deep-ocean mixing (Colbourn et al., 2015; Lord et al., 2015b; Jeltsch-Thömmes and Joos, 2020). Short-term processes (i.e., sub-centennial timescales) cannot be fit for the REF and noLAND experiments, given potential complications with and the removal of the ramp-up period. However, n=3 is already sufficient to examine the long-term uptake of anthropogenic CO₂. It should be noted that due to this decision, the values for A_i in the REF, PULSE, and noLAND experiments do not sum to 1 (Fig. 12b). This is because for the REF and noLAND experiments, a fraction of emissions was already removed via short-term processes during the excluded ramp-up period. A more accurate fit could be achieved for the PULSE experiment by increasing the number of fitted exponentials; however, we kept n=3 to maintain consistency in our analysis. In the REF ensemble, removal timescales of 386–810 years for i=1, ∼ 4–17 kyr for i=2, and ∼ 108–159 kyr for i=3 were found across the different emission scenarios (Table E2). The estimated range of removal timescales marginally changes for the PULSE and noLAND experiments (Table E2).

The centennial magnitude of the first removal timescale (i=1) suggests the involvement of ocean circulation, which we interpret as the deep-ocean processes responsible for transporting CO₂ away from the surface. The range determined in the REF, PULSE, and noLAND experiments fits well within previous estimates: for example, Archer et al. (2009a) found a range of 250–450 years, while Lord et al. (2015b) found 230–880 years. The ability of this process to remove carbon with increasing emissions is nonlinear, as the fraction of emissions removed steadily increases until ∼ 3000–4000 Pg C of cumulative emissions (Fig. 13b). A major cause of this behaviour is the slowdown of ocean circulation during the first millennium, as seen by the weakening of the AMOC. In our experiments, both the duration and magnitude of the AMOC weakening increase with larger cumulative emissions (Fig. 7e), followed by a subsequent strengthening of the AMOC after emissions cease, which overshoots the pre-industrial value. This promotes the downward transport of anthropogenic carbon and therefore increases the oceanic carbon uptake capacity (Montenegro et al., 2007; Lord et al., 2015b). In the 5000 Pg C scenario, however, the behaviour of this timescale differs significantly between the three experiments (Fig. 13). While deep-ocean processes in the REF (and to some extent, in the PULSE) experiment stabilize at a similar removal capacity as in the 4000 Pg C scenario (Table E2), its ability to remove carbon is notably reduced in the noLAND experiment (∼ 37 %). In the REF experiment, the AMOC does not begin to recover until 200–300 years after emissions end in the 5000 Pg C scenario (Fig. 7e), leading to large amount of anthropogenic carbon accumulating in the upper ocean by 1 kyr (Fig. E3). However, the AMOC eventually recovers and follows a similar trajectory to that of the 4000 Pg C experiment. This is also true for the PULSE experiment, although not shown. Meanwhile, the AMOC fails to recover at 5000 Pg C in the noLAND experiment, shutting down for the majority of the simulation (not shown). This results in the decreased ability of deep-ocean processes to absorb additional anthropogenic CO₂.

Previous studies have given a range of 1–6 kyr for sediment dissolution and 6–14 kyr for carbonate weathering (Sundquist, 1991; Archer et al., 1997; Ridgwell and Hargreaves, 2007; Colbourn et al., 2015; Lord et al., 2015b; Jeltsch-Thömmes and Joos, 2020). However, it has been emphasized that these processes should not be considered separately, given their interdependence (Lord et al., 2015b). We attribute the second calculated timescale (i=2), which is similar in all three experiments and falls within the estimated range, to the combined influence of both processes. As emissions increase to 4000 Pg C, the removal timescale shortens (Fig. 13c), and a higher fraction of emissions is removed (Fig. 13b). The enhanced ability of this timescale to remove anthropogenic CO₂ with larger cumulative emissions is primarily driven by the resulting increase in carbonate weathering, which leads to more ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ in the ocean (Fig. 9a). In low-emission scenarios, the calculated removal timescale is closer to ∼ 17 kyr, suggesting that the terrestrial weathering of carbonate rocks is the dominant process of carbon removal. But considering the fact that this timescale generally has a limited capacity to remove additional anthropogenic CO₂ (∼ 3 %–4 % in the REF ensemble; Table E2), it suggests that the combined effect of enhanced carbonate weathering (Fig. 9a) and CaCO₃ dissolution (Fig. 8e) is relatively weak in low-emission scenarios.

In high-emission scenarios, however, the timescale shortens to about 4 kyr, closer to that estimated for marine sediment dissolution. Sediment dissolution from a shoaling lysocline and CCD in the 3000–5000 Pg C scenarios reaches a similar maximum of ∼ 0.17 Pg C yr⁻¹ shortly after the first millennium in the REF experiment (Fig. 8e). As this process restores buffering capacity to the ocean and is responsible for the continued removal of anthropogenic CO₂, it is perhaps not so surprising that there is a discernible threshold in the removal timescale at these levels of emissions. Again, the removal ability of this timescale in the 5000 Pg C scenario differs significantly between the REF, PULSE, and noLAND experiments. This scenario in the noLAND ensemble accounts for the removal of nearly twice as much anthropogenic CO₂ compared to the other two (41 % vs. 18 % and 24 %). This shows that a weakening ability of deep-ocean processes to take up CO₂ is compensated by an increasing ability of sediment dissolution, to some extent (similarly shown in Archer et al., 1997), and is because seafloor CaCO₃ neutralization will act to restore the buffering capacity that has been depleted from highly alkaline deep water (Fig. E4c, g, k).

The last timescale to be considered (i=3) pertains to silicate weathering. The multi-exponential analysis shown here provides an estimate of the silicate weathering timescale that is significantly shorter than what has been previously reported. For the first time, a non-monotonic dependence of the silicate weathering timescale and capacity for carbon removal of cumulative emissions were found. This trend is seen in the REF, PULSE, and noLAND experiments, although the spread of the timescale significantly changes between the different experiments. In the 500 Pg C scenario, silicate weathering is responsible for removing ∼ 4 %–5 % of anthropogenic CO₂ across the different experiments (Table E2). However, the removal timescale is approximately 159 kyr in the REF ensemble, 151 kyr in the PULSE ensemble, and 113 kyr in the noLAND ensemble – the latter nearly a factor 2 smaller than originally estimated by Archer et al. (1997), Colbourn et al. (2015), and Lord et al. (2015b). The ability of this process to remove carbon increases up to a maximum of 15 % (seen in the REF ensemble) as emissions increase. At the same time, the removal timescale of silicate weathering in the REF ensemble decreases to ∼107 kyr in the 1000 Pg C scenario before generally increasing to ∼ 129 kyr in the 5000 Pg C scenario (Fig. 13c). This trend can also be seen in the PULSE and noLAND experiments (Fig. 13c). The initial shortening of the silicate weathering timescale with increasing emissions can be partially explained by spatially explicit weathering (Colbourn et al., 2013) and strong weathering feedbacks (Sect. 3.1). This is not unexpected, given that differences in weathering parametrizations have been found to change the estimated timescale of silicate weathering by more than 100 kyr (Uchikawa and Zeebe, 2008; Colbourn et al., 2015).

The behaviour of the silicate weathering timescale mirrors that of CaCO₃ reactions, as the timescale generally decreases as emissions rise until approximately 3000–4000 Pg C (Fig. 13c). The threshold of ∼ 3000–4000 Pg C is also seen for deep-ocean processes, which increase in removal ability until that point (Fig. 13b). The non-monotonic behaviour of the silicate weathering timescale and carbon removal ability can be attributed to the same physical mechanisms responsible for the nonlinearity observed in deep-ocean processes and CaCO₃ reactions. If the removal by either of these processes is strong, it can diminish the effectiveness of silicate weathering by reducing the amount of atmospheric CO₂ available for it to act upon. As such, the AMOC likely plays a significant role in the removal timescale of silicate weathering, as experiments with a delayed AMOC recovery (Fig. 7e) or a complete shutdown (not shown) show an increase in the timescale. Since emission scenarios greater than 5000 Pg C were not explored, it is unclear whether this trend will continue with increasing emissions. Despite differences in our estimated timescale of silicate weathering, the ability of silicate weathering to remove carbon (∼ 4 %–15 % across all experiments and scenarios; Table E2) falls reasonably close to that of previous studies, which have reported values between 4 % and 10 % (Lenton and Britton, 2006; Ridgwell and Hargreaves, 2007; Archer, 2005; Colbourn et al., 2015; Lord et al., 2015b).

3.4 Continuum of removal timescales

The multi-exponential decomposition of atmospheric CO₂ (Sect. 3.3) is, to some extent, dependent on how well a discrete number of functions can describe the continuous decline in atmospheric CO₂. This drawback has been recognized and discussed by other studies (e.g., Jeltsch-Thömmes and Joos, 2020). Therefore, we have derived an alternative method to calculate the removal timescale of anthropogenic CO₂ emissions as a function of time. With this method, the following first-order exponential decay equation describing the decline in atmospheric CO₂ (Eq. 6) is solved to estimate the removal timescale, τ(t). This was done by estimating the derivative using a central-difference approach:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}\left(t\right)}{\mathrm{d}t}}={\displaystyle \frac{-\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}\left(t\right)}{\mathit{\tau }\left(t\right)}}\text{(7)}& {\displaystyle }\mathit{\tau }\left(t\right)={\displaystyle \frac{-\mathrm{2}\mathrm{\Delta }t\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}\left(t\right)}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}(t+\mathrm{\Delta }t)-\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}(t-\mathrm{\Delta }t)}},\end{array}$$

where ΔCO₂(t) is the change in atmospheric CO₂ concentration at time t (where CO${}_{\mathrm{2}}\left(t_{\mathrm{0}}\right)={{\mathrm{CO}}_{\mathrm{2}}}^{\max }$) compared to the pre-industrial concentration (280 ppm). In doing this, the time-dependence of the calculated removal timescales was examined, along with how it changes for different levels of cumulative CO₂ emissions. We employed a rolling mean with an increasing window size over time such that larger removal timescales are averaged over longer time periods. This has been more thoroughly described in Appendix D. As in the procedures in Sect. 3.2 and 3.3, the ramp-up period of emissions was excluded, meaning that this analysis was carried out for data after the peak CO₂ concentration had been reached (CO${}_{\mathrm{2}}\left(t_{\mathrm{0}}\right)={{\mathrm{CO}}_{\mathrm{2}}}^{\max }$). This method was conducted exclusively on the REF ensemble, as one of its objectives was to provide greater confidence in the results obtained from Sect. 3.3.

Figure 14 Estimated removal timescale τ at any given time, as determined from Eq. (7) for the different emission scenarios in the REF experiment. Data have been filtered with the variable moving window discussed in Appendix D for visibility, and this is why the estimated removal timescale is not available for the full 100 kyr.

The calculated removal timescale for the REF ensemble is shown in Fig. 14a. As expected, the removal timescale is generally short in the beginning. At peak CO₂ concentrations, the estimated removal timescale ranges between ∼ 0.3 and 1.6 kyr in the REF ensemble (Fig. 14a). Between 1 and 10 kyr, there is a period in time in which the removal timescale temporarily increases for the 500 and 1000 Pg C scenarios (Fig. 14a). This behaviour of the removal timescales is due to multi-millennial processes at this time, which temporarily stabilize atmospheric CO₂ concentrations (since $\mathit{\tau }\left(t\right)\propto {\left(\frac{{\mathrm{dCO}}_{\mathrm{2}}\left(t\right)}{\mathrm{d}t}\right)}^{-\mathrm{1}}$). As this occurs around the same time as when CaCO₃ reactions are particularly significant for the removal of anthropogenic CO₂ (Fig. 13), it is likely related to the weak combined effect of deep-sea carbonate dissolution and terrestrial weathering in these scenarios (discussed in Sect. 3.3). At 50 kyr, the effect of carbonate weathering becomes minimal, and the removal timescale primarily reflects silicate weathering. Interestingly, the removal timescales for the ∼ 1000–4000 Pg C scenarios converge to a value lower than those in the 500 and 5000 Pg C scenarios (Fig. 14b). This generally confirms that the non-monotonic behaviour of the silicate weathering timescale seen in the multi-exponential analysis (Fig. 13) is not just an artefact of the fit. Averaging over the last 5000 years of τ(t) yields a mean removal timescale of approximately 89 kyr in the REF ensemble (range of 80–105 kyr; Fig. 14b). This method generally provides estimates of the silicate weathering timescale that are slightly lower and exhibit less spread compared to those obtained in Sect. 3.3. The multi-exponential analysis is helpful for understanding the role and removal ability of each process individually. However, this second method provides more confidence in calculating the silicate weathering timescale, as it is not influenced by factors such as the number of superimposed exponential functions (n) or the goodness of fit.

4 Sensitivity experiments

To explore the sensitivity of our results to different carbon cycle processes and climate sensitivity, we have performed a number of additional experiments (Table 1). The results from these additional experiments are discussed below. To provide a linear analysis of our sensitivity experiments, we limit our emission scenarios to ≤ 3000 Pg C. This is largely because high levels of cumulative emissions make it difficult to ignore secondary and cascading effects (e.g., shutdown of the AMOC), which would introduce nonlinearities and provide additional uncertainties.

Figure 15 Difference in atmospheric CO₂ concentration in the different sensitivity experiments compared to the REF experiment. (a) The effect of disabling the land carbon pool. (b) The effect of constant weathering (i.e., removal of the weathering feedback). (c) The effect of a lower (2 °C) and a higher (4 °C) equilibrium climate sensitivity compared to the standard 3.1 °C in CLIMBER-X. (d) The effect of interactive methane. Negative anomaly values represent a larger atmospheric concentration in the REF experiment. Colours correspond to the cumulative emission scenarios shown in Fig. 2.

4.1 Land carbon sink

A complete disregard of the land carbon pool results in higher atmospheric CO₂ concentrations at the start of our simulations (Fig. 15a). Peak atmospheric CO₂ concentrations are between ∼ 65 and 340 ppm higher without the land carbon pool depending on the emission scenario. These peak magnitudes generally align with previous modelling work without a land carbon pool (Lord et al., 2015b; Köhler, 2020). After peak CO₂, the difference in atmospheric CO₂ concentrations between the noLAND and REF ensembles begins to decrease for all emission scenarios up to 3000 Pg C. Conversely, the carbon inventory of the ocean increases to approximately 2300 Pg C by 1 kyr, which is roughly 300 Pg C more than in the REF ensemble (Fig. 16d). The majority of this carbon uptake in the ocean occurs at peak CO₂ concentration from increased air–sea CO₂ flux. However, at least some of this carbon is from sediment dissolution (Fig. 16g), which contributes approximately 50 Pg C more by 1 kyr. At the beginning of the noLAND ensemble, less CO₂ is absorbed on the whole, as it is taken up by relatively slow ocean and weathering processes. This leads to a higher fraction of emissions remaining in the atmosphere compared to in the REF ensemble (Fig. 17a).

It has long been recognized that land plays a crucial role in modulating the effects of weathering (Nugent and Matthews, 2012; Brault et al., 2017). However, the influence of land carbon on the long-term removal of anthropogenic CO₂ through silicate weathering remains largely understudied. Since atmospheric CO₂ concentrations are higher in the noLAND ensemble (Fig. 15a), peak temperatures and runoff are also elevated. As weathering in CLIMBER-X is driven by changes in runoff and temperature (see Eqs. 2 and 3), the noLAND ensemble produces higher weathering fluxes compared to the REF ensemble. This leads to a faster decline in atmospheric CO₂ concentrations in the long term and a shorter silicate weathering timescale. Averaging over the last 5000 years of τ(t) in Fig. 17b yields an average value of approximately 79 kyr in the noLAND ensemble. By 100 kyr, when silicate weathering becomes the dominant carbon removal process, the higher silicate weathering in the noLAND ensemble brings atmospheric CO₂ concentrations closer to those in the REF ensemble (Fig. 15a). As a result, the average fractions of emissions remaining between the REF and noLAND experiments differ by only 0.2 % at this time (Fig. 17a). The higher weathering in the noLAND ensemble therefore partially compensates for both higher initial atmospheric CO₂ concentrations and a missing land carbon pool on long timescales.

Figure 16 Net cumulative uptake of carbon in the (a–c) land, (d–f) ocean, and (g–i) sediment carbon pools across the different sensitivity experiments and emission scenarios. The stacked bar plots are cumulative, meaning that the height of the bar (rather than the bar length) in each emission scenario reflects the magnitude of carbon uptake or loss. Positive values represents carbon into the land, ocean, and sediment carbon pools. The time slices of (a, d, g) 1, (b, e, h) 10, and (c, f, i) 100 kyr were chosen to capture the effects of the different timescales. It should be noted that the time slices here are measured from the start of the simulations, which includes the emissions ramp-up period.

4.2 Silicate weathering feedback

In the absence of a silicate weathering feedback, atmospheric CO₂ concentrations cannot be restored back to pre-industrial conditions (Fig. 17a). Initially, atmospheric CO₂ concentrations are quite similar between the REF and noWEATH ensembles, as weathering is only responsible for taking up approximately 3 % of anthropogenic CO₂ within the first 1 kyr after emissions end (Fig. 3). However, the two ensembles quickly diverge as atmospheric CO₂ is additionally consumed by enhanced weathering in the REF ensemble. The difference in atmospheric CO₂ concentration between the REF and noWEATH grows to ∼120 ppm over the course of 100 kyr (Fig. 15). This is because the value to which atmospheric CO₂ concentrations return in the noWEATH ensemble increases with larger emissions (with final values between ∼ 300 and 470 ppm compared to the pre-industrial value of 280 ppm). Across the different emission scenarios, the last ∼ 14 % (range of 13 %–15 %) of anthropogenic CO₂ in the noWEATH ensemble cannot be removed without the weathering feedback (Fig. 17). This estimate is larger than what was calculated for the fraction of emissions removed by the silicate weathering timescale in the noLAND ensemble (range of 5 %–12 %; Table E2) and range provided by previous studies (range of 4 %–10 %;  Lenton and Britton, 2006; Ridgwell and Hargreaves, 2007; Archer, 2005; Colbourn et al., 2015; Lord et al., 2015b; Jeltsch-Thömmes and Joos, 2020). As atmospheric CO₂ concentrations do not decrease beyond ∼ 300–470 ppm in the noWEATH ensemble, temperatures stay significantly elevated (and near constant) until the end of the simulation. Climate is therefore very different in the noWEATH ensemble, resulting in notably different behaviour in the land, ocean, and sediment carbon pools compared to the REF ensemble (Fig. 16).

4.3 Climate sensitivity

The best estimate of equilibrium climate sensitivity is reported as 3 °C by the Intergovernmental Panel on Climate Change (IPCC; Canadell et al., 2021), which is close to the estimated ECS of 3.1 °C in CLIMBER-X. However, the IPCC also reports that the very likely range is 2 to 5 °C, so the effect of a lower and higher ECS (corresponding to values of 2 and 4 °C; Appendix A) was also tested. As seen in Fig. 15c, simulations using a lower ECS produce lower atmospheric concentrations of CO₂ (and vice versa), and higher levels of cumulative emissions produce larger differences from the REF ensemble. The differences in atmospheric CO₂ concentration can largely be traced back to changes in the land carbon inventory (Fig. 16a). In the case of ECS4, terrestrial vegetation takes up approximately the same amount of carbon as in the REF ensemble (not shown). However, increased warming causes soil carbon to become a larger source of carbon, resulting in higher atmospheric CO₂ concentrations and lower land carbon across the different emission scenarios over 100 kyr. This behaviour is indicative of positive feedback as a result of warming-induced soil carbon release. For ECS2, the storage of anthropogenic CO₂ in the vegetation and soil carbon pools is positive for the entire duration of the simulation regardless of the emission scenario (not shown). Despite significant variations in the partitioning of carbon between the land and ocean carbon pools (Fig. 16), differences in atmospheric CO₂ concentrations between the ECS2 and ECS4 ensembles and the REF ensemble are within ∼ 5 ppm by 10 kyr and remain low for the rest of the simulation (Fig. 15c).

At first glance, this might suggest that the impact of ECS on atmospheric CO₂ concentration is minimal after 10 kyr regardless of the cumulative emission scenario. Similar behaviour was also seen in Shaffer et al. (2009). However, this conclusion is only valid for a short period. In Fig. 17c, d, higher ECS values lead to a faster decline in atmospheric concentration and a shorter removal timescale of silicate weathering. Although atmospheric CO₂ concentrations are initially higher in the ECS4 ensemble, they fall below those of the REF ensemble around 50 kyr due to substantial consumption by weathering (Fig. 15c). The reverse is true for ECS2. Averaging over the last 5000 years for τ(t) results in a mean silicate weathering timescale of 153 kyr for an ECS of 2 °C and 54 kyr for an ECS of 4 °C. This generally corroborates the findings by Colbourn et al. (2015), who found that a lower ECS would give a longer timescale for silicate weathering. With that in mind, the default ECS in Colbourn et al. (2015) and Lord et al. (2015b) was 2.64 °C, indicating that climate sensitivity could account for at least some of the differences in our estimated silicate weathering timescale (Fig. 13).

Figure 17 Comparison of (a, c) the atmospheric lifetime of anthropogenic CO₂ described by IRF${}_{{\mathrm{CO}}_{\mathrm{2}}}$(t) and (b, d) the removal timescale of anthropogenic CO₂ in the different sensitivity experiments. The mean of the 500–3000 Pg C scenarios is shown. As described in Sects. 3.2 and 3.4, the emissions ramp-up period was excluded from the calculation of IRF${}_{{\mathrm{CO}}_{\mathrm{2}}}$(t) and τ(t). Data in (b) and (d) have been filtered with the variable moving window discussed in Appendix D for visibility, and this is why the estimated removal timescale is not available for the full 100 kyr.

4.4 Interactive methane

Global warming from anthropogenic CO₂ emissions can cause additional methane emissions from natural sources that can result in an additional positive-feedback loop. Although methane has a relatively short atmospheric lifetime compared to our experiment duration (∼ 9.8 years; Voulgarakis et al., 2013), it has higher global warming potential compared to CO₂. An increase in methane concentration could therefore affect the capacity of different Earth system components to take up carbon over time. For the first time, we investigate the impact of interactive methane on the long-term carbon cycle. A simple methane cycle has been implemented in CLIMBER-X that explicitly simulates the natural methane emissions from land resulting from anaerobic respiration in wetlands and peatlands (Willeit and Ganopolski, 2016). Permafrost is treated in the methane model through the representation of anaerobic decomposition in saturated soils – a process that can take place in permafrost regions. Atmospheric methane concentration is then computed with the assumption of a constant atmospheric lifetime of methane (Willeit et al., 2023).

In the zero-emissions scenario, the CH₄ concentration stays largely around 600–700 ppb (close to the constant 600 ppb set in the REF ensemble; Fig. E7b). For the other scenarios (500–3000 Pg C), atmospheric CH₄ concentrations rise over the first ∼ 200 years, reaching a peak of approximately 1620 ppb in the 3000 Pg C scenario. With interactive methane, atmospheric CO₂ concentrations can be up to 20 ppm higher than in the REF ensemble (Fig. 15d). During much of the first 1 kyr, the differences in atmospheric CO₂ concentrations between the intCH4 and REF ensembles increase with higher emissions. However, this pattern changes as atmospheric CO₂ concentrations decline and does not hold beyond the first 1 kyr. The global mean temperature in the 3000 Pg C scenario is up to 0.5 °C higher in the intCH4 experiment compared to the REF experiment (Fig. E7a), but the difference remains largely negligible in the low-emission scenarios. While total vegetation carbon is similar between the REF and intCH4 ensembles, soil carbon in the intCH4 ensemble is notably lower across the different emission scenarios, causing the land carbon pool to take up less than the REF ensemble overall (Fig. 16a–c). Carbon stored in peatland is responsible for some of these differences, as carbon is lost from warming-induced soil respiration. Beyond the first 5 kyr, interactive methane has a ± 5 ppm effect on atmospheric CO₂ concentration (Fig. 15d). Temperature differences between the intCH4 and REF ensembles also tend to decrease around this time, and despite small variations in how temperature evolves, the experiments largely follow each other for the rest of the simulation (Fig. E7a). In this way, the lifetime of anthropogenic CO₂ and removal timescales in the intCH4 ensemble do not differ significantly from those in the REF ensemble over long timescales (Fig. 17a, b), suggesting that the indirect effect of CH₄ on the climate (through CO₂) is likely negligible over such periods.

5 Discussion and conclusions

The impacts of anthropogenic climate change will likely be felt for tens of thousands if not hundreds of thousands of years – long after the anthropogenic emission of CO₂ has ended. By simulating the long-term response to idealized emission scenarios, we highlighted how processes that remove carbon from the atmosphere change with different levels of cumulative emissions. With this, we have examined the capacity of these processes and extent to which they take up carbon on different timescales. To some extent, the results of this study align with previous investigations. For example, we have shown that the magnitude of cumulative emissions will predominantly control the long-term response, supporting the path independence found by previous work (Matthews et al., 2009; Zickfeld et al., 2012; Herrington and Zickfeld, 2014). The decline in anthropogenic CO₂ over time seen here generally falls within estimates given by the multi-model assessment in Archer et al. (2009a) and other studies: approximately 10 % of anthropogenic emissions will persist longer than 10 000 years and ∼ 4 % beyond 100 000 years. The response of deep-ocean processes and CaCO₃ reactions, as well as their ability to remove anthropogenic CO₂ at different levels of cumulative emissions, also supports earlier investigations (Lord et al., 2015b). However, the extensive sensitivity analysis, along with the higher resolution and complexity of CLIMBER-X, have lead to several new findings.

1. The removal timescale of silicate weathering is significantly influenced by a number of factors. In this study, we quantified the effects of climate sensitivity, the land carbon pool, weathering feedbacks, and dynamically changing methane concentrations on the atmospheric lifetime of anthropogenic CO₂. The sensitivity analysis revealed that climate sensitivity has the most significant influence on the silicate weathering timescale (Fig. 17). However, most factors were generally found to impact the silicate weathering timescale, albeit in opposing ways. Our results show that higher climate sensitivity, stronger silicate weathering feedbacks, a spatially explicit weathering scheme, and no land carbon pool reduce the removal timescale of silicate weathering, whereas lower climate sensitivity and a land carbon pool increase this timescale. This means that in theory, the same silicate weathering timescale could be obtained with no land carbon and either weaker weathering feedbacks or lower ECS. This emphasizes the importance of incorporating as many components of the carbon cycle as possible and increasing model complexity in Earth system models to more accurately represent the long-term future behaviour of atmospheric CO₂. Furthermore, the global heterogeneity of carbonate and silicate weathering with increasing emissions (Fig. 11) also showed the importance of highly active weathering regions for the global weathering flux (Sect. 3.1.3), highlighting the need for spatially explicit schemes. Improving the silicate weathering timescale would require more accurate constraints on both these factors and the global silicate weathering flux and would require a deeper understanding of how weathering responds to increasing global temperatures.

2. The removal timescale of silicate weathering is significantly shorter than previously estimated. The two methods used in this study provided significantly shorter removal timescales for silicate weathering compared to previous research, providing a best estimate of 80–105 kyr depending on the emission size. This discrepancy is likely due to differences in climate sensitivity, weathering feedbacks, the spatially explicit weathering scheme, and the land carbon pool. With this in mind, our estimates here do not address (potentially large) uncertainties related to the sensitivity of silicate weathering processes, such as erosion limitation or net primary production on land. Nevertheless, it is worth noting that ∼ 20–30 years has passed since the constant 200–400 kyr range was first established for the silicate weathering timescale using lower-complexity models (Sundquist, 1991; Archer et al., 1997, 1998), and our results indicate that it is time to re-examine this estimate. Considering the large uncertainty that exists in the literature concerning its value and the progress made in model resolution and process descriptions in our model and others since the last LTMIP, another intercomparison project could be beneficial.

3. The removal timescale of silicate weathering is dependent on the magnitude of cumulative emissions. A major finding in Lord et al. (2015b) was that the removal ability and timescale of silicate weathering were approximately invariant with respect to cumulative emissions. This is challenged here, as it was shown that both are non-monotonic with increasing emissions. In all experiments that were analyzed (REF, PULSE, and noLAND), there is an initial shortening of the silicate weathering timescale with increasing cumulative emissions until ∼ 3000–4000 Pg C. For the 5000 Pg C scenario, however, it was found (via two methods) that the timescale increases. Whether this trend persists for emissions beyond 5000 Pg C remains to be determined. Future research would benefit from exploring the connection between the extended AMOC decline and the long-term removal of anthropogenic CO₂ via silicate weathering.

4. Land carbon has a significant long-term influence on the future evolution of atmospheric CO₂. Despite variations in vegetation and soil carbon response with increasing emissions, the terrestrial storage of anthropogenic carbon is positive for the duration of our simulations and is responsible for storing approximately 4 %–13 % of anthropogenic carbon by the end of the simulation. The land's capacity to absorb anthropogenic CO₂ in the simulations buffers higher CO₂ concentrations and, therefore, peak global temperatures. As a consequence, the inclusion of land carbon cycle processes effectively reduces weathering fluxes, thereby increasing the removal timescale of silicate weathering. Experiments that exclude land carbon would show a shorter silicate weathering timescale than those that include land carbon. As silicate weathering is the ultimate control of the removal of anthropogenic CO₂, the extent to which the land buffers peak temperature in a model can significantly influence the atmospheric lifetime of anthropogenic CO₂.

5. The methane cycle will likely have a negligible effect on long timescales (> 10³ years). Rising global temperatures will increase methane levels, potentially amplifying positive feedbacks and impacting how Earth system components absorb carbon over time. These experiments show that a dynamic methane cycle in CLIMBER-X can contribute up to an additional ∼ 25 ppm in atmospheric CO₂ concentration. However, this effect is mainly limited to the first millennium after emissions end and is largely negligible on long timescales. As such, the methane cycle does not contribute to substantially different removal timescales, leading to the same rate of anthropogenic CO₂ removal over time as in the REF ensemble.

There are several assumptions in our study that should be considered when generalizing the long-term response of the climate–carbon system. Certain aspects of the carbon cycle in CLIMBER-X are simplified or not considered, as current computational restrictions limit the inclusion of all carbon cycle processes across all timescales. For example, CLIMBER-X does not consider deep-sea methane hydrates or nitrogen limitation in tropical or boreal forests, despite their potential impacts on future carbon uptake (Archer et al., 2009b; Norby et al., 2010). Some general uncertainties and poor empirical constraints should also be considered, as they can inhibit an accurate description of some processes. One such example is the strength of the biological pump in a warming world. Although state-of-the-art Earth system models do not agree on the overall response of marine primary production to increasing CO₂ levels on centennial timescales (Laufkötter et al., 2015), widespread ocean acidification is still expected to have some influence (Mora et al., 2013). This is not considered in our model, which relies on an extended NPZD (nutrients–phytoplankton–zooplankton–detritus) description, typical of models at this level of complexity. The weathering scheme used in CLIMBER-X is significantly more advanced than in other models of the same class; however, it does not explicitly account for changes in vegetation cover or net primary production. While this could be an important simplification, it is presently unclear to what extent it may affect long-term CO₂ dynamics. As 10 % of the land area in our simulations is responsible for more than a third of global silicate weathering fluxes, small changes in the controlling parameters in these areas can also significantly impact total weathering fluxes (Hartmann et al., 2009b). Like many other models, CLIMBER-X also does not account for the potentially significant effect of soil shielding, given the limited quantity of data and the ongoing uncertainties regarding their impact. However, there is some evidence to suggest that erosion limitation could weaken the strength of the weathering feedback (Goddéris et al., 2008; Maher and Chamberlain, 2014; Maffre et al., 2023), which may increase the estimated removal timescale of silicate weathering. These limitations highlight potential directions for further work.

A few studies (Lenton and Britton, 2006; Lord et al., 2015b; Köhler, 2020; Couplet et al., 2025) have looked at idealized emission scenarios with more than ∼6000 Pg C, corresponding to unconventional fossil fuels. This was not considered here, partially because the potential of unconventional fossil fuels is still speculative but also because such high levels of cumulative emissions make it difficult to ignore additional changes in the Earth system, such as sea level rise. Simulations in this study were carried out using prescribed present-day ice sheets to maintain consistency with previous research. However, the melting of the Greenland and West Antarctic ice sheets, inevitable in high-emission scenarios, can introduce additional feedback mechanisms into the carbon cycle (Wadham et al., 2019). These effects could be significant for the range of emissions considered here. Moreover, significant global warming is expected to trigger the crossing of one (or more) critical thresholds in the Earth system (Lenton et al., 2008; Armstrong McKay et al., 2022; Wunderling et al., 2024), which could affect the long-term capacity of different components to absorb carbon over time. In addition to this, it has been assumed that the pre-industrial global carbon cycle was in equilibrium, such that the silicate weathering feedback will eventually return the CO₂ concentrations to 280 ppm after an anthropogenic perturbation in the carbon cycle. This assumption can also be explicitly seen in the zero-emissions scenario in Fig. 2. While this is a standard assumption for investigations into the long-term response of the climate-carbon system to perturbations in the carbon cycle, it is perhaps not so appropriate for the very long timescales considered in this work. The implications of this assumption will be considered in a forthcoming paper. Future work with the CLIMBER-X model will investigate the long-term effects of anthropogenic CO₂ on the dynamic evolution of Northern Hemisphere ice sheets.

Appendix A: Emulation of different climate sensitivities

Previous studies have implemented straightforward ways to scale ECS (e.g., Colbourn et al., 2015). In CLIMBER-X, different climate sensitivities can be achieved by scaling the equivalent CO₂ in the long-wave radiation scheme in the following way:

$$\begin{array}{}\text{(A1)}& {\mathrm{CO}}_{\mathrm{2}}^{\mathrm{eq},\mathrm{ecs}}=e^{\left(\mathit{\alpha }\ln {\mathrm{CO}}_{\mathrm{2}}^{\mathrm{eq}}+(\mathrm{1}-\mathit{\alpha })\ln {\mathrm{CO}}_{\mathrm{2}}^{\mathrm{ref}}\right)},\end{array}$$

where α is the scaling coefficient and CO${}_{\mathrm{2}}^{\mathrm{ref}}=\mathrm{280}$ ppm. The equivalent CO₂ concentration for radiation is CO${}_{\mathrm{2}}^{\mathrm{eq}}$, which also includes the effect of other greenhouse gases (Willeit et al., 2023). It was determined a priori that α=0.66 yielded an ECS of 2 °C, whereas α=1.33 yielded an ECS of 4 °C. This method has been more extensively applied in Höning et al. (2024) and Kaufhold et al. (2025b). More information about this ECS scaling technique in CLIMBER-X, as well as a comparison with state-of-the-art Earth system models, can be found in Kaufhold et al. (2025b).

Appendix B: Fraction of cumulative fluxes from the atmosphere

The fraction of the different carbon reservoirs that contribute to the removal of anthropogenic CO₂ over a span of 100 kyr is assessed using Eq. (B3). This is similar to calculations for airborne, land-borne, ocean-borne, and sediment-borne fractions (Jones et al., 2013; Vakilifard et al., 2022), as the cumulative carbon fluxes of the atmosphere, land, ocean, and weathering are normalized by cumulative anthropogenic emissions, E (Pg C). However, we focus strictly on changes in cumulative carbon flux from the atmosphere relative to the pre-industrial state instead of changes in the carbon inventory, as our experiments were conducted with an open carbon cycle. From Eq. (1), we get

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

which, integrated in time, gives

$$\begin{array}{}\text{(B2)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{anth}}\mathrm{d}t& =\left({\mathrm{C}}_{\mathrm{atm}}\left(t\right)-{\mathrm{C}}_{\mathrm{atm}}\left(\mathrm{0}\right)\right)-\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{lnd}}\mathrm{d}t\\ & -\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{ocn}}\mathrm{d}t+\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{weath}}\mathrm{d}t-\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{volc}}\mathrm{d}t.\end{array}\end{array}$$

As we only consider times after the cessation of anthropogenic emissions (i.e., ${\int }_{\mathrm{0}}^t{F}_{\text{anth}}\mathrm{d}t=E$, total cumulative emissions), the equation above becomes

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}\mathrm{1}& ={\displaystyle \frac{\mathrm{1}}{E}}(\left({\mathrm{C}}_{\mathrm{atm}}\left(t\right)-{\mathrm{C}}_{\mathrm{atm}}\left(\mathrm{0}\right)\right)-\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{lnd}}\mathrm{d}t-\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{ocn}}\mathrm{d}t\\ & +\underset{\mathrm{0}}{\overset{t}{\int }}\left(F_{\text{weath}}^{\mathrm{carb}}+F_{\text{weath}}^{\mathrm{sil}}\right)\mathrm{d}t-\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{volc}}\mathrm{d}t).\end{array}\end{array}$$

The individual components seen in Fig. 3 are then computed as follows:

$$\begin{array}{}\text{(B4)}& {\displaystyle }{\mathrm{\Phi }}_{\text{atm}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{E}}\left({\mathrm{C}}_{\mathrm{atm}}\left(t\right)-{\mathrm{C}}_{\mathrm{atm}}\left(\mathrm{0}\right)\right)\equiv {\mathrm{IRF}}_{{\mathrm{CO}}_{\mathrm{2}}}\left(t\right),\text{(B5)}& {\displaystyle }{\mathrm{\Phi }}_{\text{lnd}}\left(t\right)=-{\displaystyle \frac{\mathrm{1}}{E}}\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{lnd}}\mathrm{d}t,\text{(B6)}& {\displaystyle }{\mathrm{\Phi }}_{\text{ocn}}\left(t\right)=-{\displaystyle \frac{\mathrm{1}}{E}}\left(\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{ocn}}\mathrm{d}t-\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{weath}}^{\text{carb}}\mathrm{d}t\right),\text{(B7)}& {\displaystyle }{\mathrm{\Phi }}_{\text{weath}}^{\text{sil}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{E}}\left(\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{weath}}^{\text{sil}}\mathrm{d}t-\underset{\mathrm{0}}{\overset{t}{\int }}{F}_{\text{volc}}\mathrm{d}t\right).\end{array}$$

Fluxes (in Pg C yr⁻¹) are relative to the average flux over the reference simulation to account for model drift. We add and subtract $F_{\text{weath}}^{\text{carb}}$ from F_ocn and F_weath, respectively, as carbonate weathering additionally transfers carbon into the ocean and causes ocean outgassing. Volcanic outgassing, F_volc, was additionally subtracted from F_weath to isolate silicate weathering consumption.

Appendix C: Calculation of the Revelle factor

The Revelle factor ξ quantifies the ocean's ability to absorb CO₂. It essentially compares the ratio of fractional changes in CO₂ to the fractional changes in DIC (Sarmiento and Gruber, 2006), making it a useful metric to assess the sensitivity of our simulations to anthropogenic CO₂ perturbations (Broecker et al., 1979). Given that there are competing definitions in the literature (see:  Egleston et al., 2010), we explicitly calculate the Revelle factor in the following way:

$$\begin{array}{}\text{(C1)}& \mathit{\xi }={\displaystyle \frac{\partial \ln p{\mathrm{CO}}_{\mathrm{2}}}{\partial \ln \mathrm{DIC}}}.\end{array}$$

Appendix D: Rolling mean with variable window size

One disadvantage of the formulation for τ(t) in Sect. 3.4 is that it is sensitive to small-scale variability and noise in CO₂(t). A simple rolling mean is not sufficient to smooth this, as the removal timescale is expected to increase over time, eventually converging to an 𝒪(10⁵⁻⁶) timescale that is representative of silicate weathering. For this reason, we use a forward-rolling mean to remove noise on both short and long timescales. While a rolling mean is generally defined with a fixed window size, we modified the standard equation to have a variable window size. This can be written as follows:

$$\begin{array}{}\text{(D1)}& {\mathrm{CO}}_{\mathrm{2}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{w\left(t\right)}}\sum _{i=t}^{t+w\left(t\right)}{\mathrm{CO}}_{\mathrm{2}}\left(i\right),\mathrm{where}{\mathrm{CO}}_{\mathrm{2}}\left(t_{\mathrm{0}}\right)={{\mathrm{CO}}_{\mathrm{2}}}^{\max },\end{array}$$

where w(t) is the window size as a function of time, t is time, CO₂(t) is the atmospheric CO₂ concentration at time t, and CO${}_{\mathrm{2}}^{\max }$ is the peak CO₂ concentration. Different functions were tested for w(t), but following a simple, linearly increasing window of $w\left(t\right)=\mathrm{0.5}t-\mathrm{70}$ provided the best results, as it still preserved the overall trend on multiple timescales. Given data limitations, the rolling window could not extend beyond 100 kyr due to the forward-rolling technique. For this reason, the data are cut around approximately t≈ 67 kyr, when the window size becomes approximately w(t)≈ 33 kyr, and this is why data for the full 100 kyr are not shown in Figs. 14 and 17b, d.

Appendix E: Supporting material

Figure E1 Distribution of the 14 lithological classes implemented in CLIMBER-X for the present day. These are based on the GLiM lithological map (Hartmann and Moosdorf, 2012). The dominant rock type in a given grid cell is shown in (a). The spatial distribution of grid cell lithological fractions, β(ℓ), is given for (b) metamorphic rocks, (c) acid plutonic rocks, (d) basic plutonic rocks, (e) intermediate plutonic rocks, (f) pyroclastics, (g) mixed sedimentary rocks, (h) siliciclastic sedimentary rocks, (i) unconsolidated sediments, (j) acid volcanic rocks, (k) basic volcanic rocks, (l) intermediate volcanic rocks, (m) loess, (n) carbonate sedimentary rocks, and (o) evaporites. Note that the colour bar range for different lithological fractions changes between lithological classes.

Table E1 Lithological classes in CLIMBER-X. These are based on the lithological map of GLiM (Hartmann and Moosdorf, 2012), and the activation energy of silicates is from Hartmann et al. (2014). The evaporites class (ev) is used only to compute phosphorus fluxes in CLIMBER-X and is therefore also not considered here.

Figure E2 Pathways of emitted carbon as described by three different functions. A fourth pathway called “pulse”, not shown here, assumes that all anthropogenic emissions are released in the first time step of the simulation.

Figure E3 Annual mean dissolved inorganic CO₂ in three oceanic basins at 1 kyr. The top row (a)–(c) shows the zonal mean cross section of the Atlantic Ocean, the middle row (d)–(f) shows the Pacific Ocean, and the bottom row (g)–(i) shows the Indian Ocean. The time slice here is measured from the start of the simulations, meaning that it includes the emissions ramp-up period.

Figure E4 Annual mean alkalinity in three oceanic basins at 1 kyr. The top row (a–c) shows the zonal mean cross section of the Atlantic Ocean, the middle row (d–f) shows the Pacific Ocean, and the bottom row (g–i) shows the Indian Ocean. The time slice here is measured from the start of the simulations, meaning that it includes the emissions ramp-up period.

Figure E5 Annual mean dissolved carbonate ion concentration in three oceanic basins at 1 kyr. The top row (a–c) shows the zonal mean cross section of the Atlantic Ocean, the middle row (d–f) shows the Pacific Ocean, and the bottom row (g–i) shows the Indian Ocean. The time slice here is measured from the start of the simulations, meaning that it includes the emissions ramp-up period.

Figure E6 Contribution of temperature and runoff to the silicate weathering response in CLIMBER-X. (a) Change in silicate weathering flux from the pre-industrial value (12.3 Tmol C yr⁻¹) for the 5000 Pg C scenario. The contributions from temperature and runoff were broken down for the 5000 Pg C scenario, and they do not linearly add to the calculated model response, as the weathering relation in Eq. (3) is nonlinear. (b) Response in silicate weathering as a function of global mean temperature. Like Fig. 10, the fit of these trajectories has been plotted above using Numpy.Polyfit in Python. The relative contributions of temperature and runoff to the silicate weathering response were determined from the 5000 Pg C scenario, i.e., the trajectories from (a).

Figure E7 Changes in (a) global temperature and (b) atmospheric CH₄ concentration in the intCH4 experiment compared to the REF experiment. The dotted grey line here in (b) shows the constant CH₄ concentration used in the REF experiment (600 ppb). Colours correspond to the cumulative emission scenarios shown in Fig. 2.

Table E2 Fitting parameters used for a multi-exponential decay function describing atmospheric CO₂ concentration in the REF, PULSE, and noLAND experiments. Here, the fraction of emissions (A_i) taken up over the removal timescales (τ_i) associated with ocean invasion (i=1), CaCO₃ reactions (i=2), and silicate weathering (i=3) are shown. We highlight the fact that this is the total fraction of emissions removed after peak CO₂ concentrations have been reached. Anthropogenic CO₂ removed during the ramp-up of emissions is not considered in this analysis. The fraction of emissions removed from short-term (sub-centennial timescales) processes can be broadly estimated from $\mathrm{1}-\sum A_i$. While n=3 was used for consistency across experiments, we acknowledge that a more accurate fit for the PULSE ensemble could theoretically be achieved by increasing the number of fitted exponentials (as in Colbourn et al., 2015; Lord et al., 2015b).