A New Method for Characterising Carbon Cycle Sensitivities in Earth System Models

1 Introduction

Since the pioneering work of Keeling (1960), published in Tellus, it has been well established that the concentration of atmospheric CO₂ has been increasing rapidly since the start of the Industrial Revolution, which has in turn led to a rise in global average temperatures. The growth of CO₂ is due to human activity, primarily the burning of fossil fuels, but with a significant contribution also from deforestation and other land-use change. As global populations, economies, and industries have expanded, their increasing energy and material demands have resulted in anthropogenic CO₂ emissions that increased almost exponentially over the past 200 years. However, the atmospheric CO₂ concentration has only increased by less than half of the anthropogenic emissions (Friedlingstein et al., 2024). This discrepancy is the result of significant natural carbon sinks on the land and in the ocean (Bolin, 1977; Bolin et al., 1989; Revelle and Suess, 1957).

However, land and ocean carbon sinks are not static. Terrestrial and marine carbon uptake varies over seasonal to millennial time-scales, and the magnitude of these variations depends on the responses of physical and biological processes to changes in both atmospheric composition and physical climate. These responses can either slow or accelerate the rate of climate change, leading to ‘climate-carbon cycle feedbacks’ which are a key source of uncertainty in future projections (Friedlingstein et al., 2014). The role of climate-carbon cycle feedbacks is also very evident in past climate variations. Most notably, during glacial–interglacial periods, changes in CO₂ drive a significant part of the glacial–interglacial changes in temperature (Petit et al., 1999; Sigman and Boyle, 2000).

The gross uptake of carbon by land and ocean depends on the concentration of CO₂ in the atmosphere. Terrestrial ecosystems sequester CO₂ through photosynthesis, which is stimulated by increased atmospheric CO₂ if nutrients are also available (Walker et al., 2021). Similarly, diffusion of CO₂ from the atmosphere into the surface layer of the ocean is controlled by the difference in partial pressures (pCO₂) between the ocean and atmosphere. As a result, since the start of the industrial revolution, land and ocean sinks have kept pace with the near-exponential increase in anthropogenic emissions, evidenced by the long-term stability of the airborne fraction of CO₂ emissions (Friedlingstein et al., 2024; Jones and Friedlingstein, 2020). However, the fraction of anthropogenic emissions removed by land and ocean carbon sinks is projected to decrease under climate change (Canadell et al., 2021).

On the land, CO₂ fertilisation relies heavily on sufficient nutrient supply, which, despite large anthropogenic production of reactive nitrogen (Zhu et al., 2025), is not always available (LeBauer and Treseder, 2008). Nutrient cycles in some areas remain relatively closed. Old-growth tropical forests, for example, are typically limited by phosphorus (Walker and Syers, 1976), which restricts their potential CO₂ uptake in the future. Changes in climate may also take plants outside their natural ranges. In the ocean, rising sea surface temperatures and increasing dissolved inorganic carbon (DIC) concentrations reduce the rate of CO₂ uptake (Egleston et al., 2010; Revelle and Suess, 1957). These effects contribute to ‘sink saturation’ and declining sink efficiency under increasing CO₂, which may reduce carbon uptake and accelerate global climate change in the future (Canadell et al., 2021).

Nevertheless, it was not until the early years of this millennium that general circulation models (GCMs) of the climate began to include interactions between the climate and carbon cycle (Cox et al., 2000; Friedlingstein et al., 2001). Indeed, Tellus published an intercomparison of the first two climate-carbon cycle GCMs in 2003 (Friedlingstein et al., 2003). These models were the forerunners of what are now called ‘Earth System Models (ESMs)’, and ultimately led to the ‘Coupled Climate-Carbon Cycle Model Intercomparison Project (C⁴MIP)’ (Arora et al., 2013, 2020; Friedlingstein et al., 2006). C⁴MIP compares results from multiple ESMs using a quasi-equilibrium linear formulation that was first proposed in that 2003 paper in Tellus. This ‘$\beta –\gamma$’ formulation has the advantage of being very simple to use, but is unable to capture vital aspects of sink saturation.

This paper for the Tellus 75^th Anniversary issue aims to summarise work on coupled-climate carbon cycle modelling (Section 2), and also to propose a long-overdue update to the ‘$\beta –\gamma$’ characterisation of carbon cycle sensitivities to CO₂ and climate, to capture sink saturation effects (Section 3).

2 The Evolution of Coupled Climate-Carbon Cycle Models

As described above, the importance of climate-carbon interactions had been revealed throughout the last 30 years of the 20th century. Although ‘full-form’ GCMs were beginning to include coupling between the climate and the ocean carbon cycle on long timescales (Sarmiento and Toggweiler, 1984), these had not yet included the land carbon cycle. GCM land surface schemes, which had traditionally focused on calculating land-atmosphere fluxes of energy, water, and momentum, were diversifying to also calculate plant photosynthesis (Cox et al., 1999; Sellers et al., 1996a). Initially, this was to improve the simulation of water fluxes over land, which is heavily influenced by stomatal pores on plant leaves (Sellers et al., 1997), by making use of empirical relationships between photosynthesis and stomatal conductance (Ball et al., 1987). These second-generation land surface models therefore coupled the water and carbon fluxes, but still required prescribed fixed maps of vegetation types and land-surface parameters (Pitman, 2003; Sellers et al., 1996b).

In parallel, a separate group of researchers had begun to develop dynamic global vegetation models (DGVMs), with the aim of simulating the changing geographical distribution of vegetation types under climate change (Cramer et al., 2001). These models were derived from biogeographical models developed for palaeoclimate reconstructions (Haxeltine and Prentice, 1996; Prentice et al., 1992), but were also increasingly based on the principles of carbon and water balance (Woodward, 1987). Before the end of the 20th century, a first attempt was made to asynchronously couple a DGVM to a GCM climate model (Betts et al., 1997). However, the pieces were now in place for a GCM of the fully coupled climate-carbon system, which are today referred to as ‘ESMs’ (Randall et al., 2019).

At the beginning of the 21st century, the first ESMs were published. A group based at the UK Met Office-Hadley Centre developed a variant of the HadCM3 climate model which included representations of the ocean carbon cycle (Palmer and Totterdell, 2001) and land carbon cycle (Cox et al., 1999), and also a simple representation of vegetation dynamics (Cox, 2001). A fully coupled climate-carbon cycle run was compared to a model run with prescribed CO₂ concentrations, to isolate the impact of climate-carbon feedbacks on the rate of CO₂ increase and global warming. Alarmingly, this study reported that climate-carbon cycle feedbacks increased projected atmospheric CO₂ by 270 ppmv, and global warming from 4$\mathrm{°}$C to 5.5$\mathrm{°}$C, by the end of the 21st century (Cox et al., 2000). The acceleration of climate change in this model was primarily due to the suppression of the global land carbon sink by climate change, which ultimately transitioned to a carbon source from about 2050 onwards, due to increased rates of soil respiration and Amazon forest dieback under global warming (Cox et al., 2004).

Independently, a group at the IPSL in France was also developing a climate-carbon cycle model (Friedlingstein et al., 2001), making use of the IPSL-OGCM1 ocean carbon cycle model (Le Quéré et al., 2000) and the SLAVE terrestrial carbon cycle model (Friedlingstein et al., 1995). Under a prescribed 1% per year increase in atmospheric CO₂, this model also simulated a suppression of the land carbon sink by climate change, with land carbon uptake at the time of doubling CO₂ reduced by about a quarter. These first two ESM simulations (Cox et al., 2000; Friedlingstein et al., 2001) suggested that climate change would suppress natural carbon sinks and therefore provide an amplifying positive feedback on climate change. However, the magnitude of the positive feedback differed markedly in the two models. In order to isolate the reasons for this, a collaborative study was undertaken between the French and UK groups. The resulting paper, published in Tellus (Friedlingstein et al., 2003), identified a stronger negative response of the land carbon cycle to climate change in the UK model as the key reason for the differences between them.

Prompted by these initial studies which suggested that carbon cycle feedbacks on future climate could be large but were uncertain, the C⁴MIP was set up to enable a cleaner comparison between the results from the growing number of coupled climate-carbon cycle models. The first phase of C⁴MIP used a common emissions scenario (SRES A2) and required global atmospheric CO₂ concentration to be simulated internally within each model. All models produced a positive climate-carbon feedback on climate change, but with a range that differed by an order of magnitude from an extra 20 to 200 ppmv of CO₂ by 2100 (Friedlingstein et al., 2006). Two further phases of C⁴MIP have since been carried-out, based, respectively, on the CMIP5 and CMIP6 generations of ESMs (Arora et al., 2013, 2020). These latter C⁴MIP studies used a simpler experimental design in which the atmospheric CO₂ concentration is prescribed to increase at 1% per year, and model runs are carried out with (‘COU’; coupled runs) and without (‘BGC’; biogeochemical runs) climate effects on the carbon cycle, to separate the direct-CO₂ effects (see Figure 1; Jones et al. (2016)). Each generation of C⁴MIP models has seen further developments. For example, additional nutrient limitations on land and ocean carbon sinks were implemented for CMIP6, and the range of carbon cycle projections has generally been reduced (Canadell et al., 2021). However, significant uncertainty remains, especially in the response of the land carbon cycle, which is also impacted by land-use change (Arora et al., 2020). As a result, methods to separate impacts on the carbon cycle remain important.

Figure 1

3 A More Robust Diagnosis of Carbon Cycle Feedbacks

The assumption that increases in land and ocean carbon storage are proportional to changes in atmospheric CO₂ concentration is inconsistent with our understanding and modelling of the processes that lead to sink saturation at high CO₂, as explained in Section 1. Also, the standard $\beta –\gamma$ formulation does not represent the timescales involved in carbon uptake by the land and ocean, instead assuming an instantaneous adjustment of carbon stores. It therefore predicts the same increase in land and ocean carbon storage as a function of the CO₂ concentration and global warming, regardless of the scenario used to get to that state. This is at odds with the outputs of complex ESMs, which show a dependence of carbon storage on the rate of increase of CO₂ (Canadell et al., 2021; Melnikova et al., 2021; Varney et al., 2023). In this paper, we seek to improve the standard $\beta –\gamma$, by representing sink saturation, and also the characteristic timescales of land and ocean carbon uptake.

We aim to improve the standard formulation by relaxing each of these assumptions in turn, in a way that minimises the additional complexity introduced. We design our new approach so that it requires no additional model runs in addition to those required by the standard approach. That is, we continue to deduce the carbon-concentration sensitivity from the BGC runs, and the carbon-climate sensitivity from the difference between the fully coupled and the biogeochemically coupled runs.

3.2 Results

Figure 2 shows global changes from the coupled runs (COU) of the CMIP6 generation of C⁴MIP (Arora et al., 2020). As for the previous CMIP5 generation of C⁴MIP (Arora et al., 2013), all modelling groups were required to prescribe a 1% per year increase in CO₂ from a preindustrial CO₂ concentration of about 284 ppmv to at least four times that value (panel (a)). As the models have rather different climate sensitivities to CO₂ (Nijsse et al., 2020), this common CO₂ scenario still leads to a wide range of global warming projections, ranging from about 4K (NorESM2-LM, MIROC-ES2L, GFDL-ESM4) to about 7K (CanESM5, UKESM1-0-LL, IPSL-CM6-LR) at 4×CO₂ (panel (b)). The CO₂ increase and global warming lead to significant changes in the land and ocean carbon cycles in the models. On the land, global NPP increases in all models, primarily due to CO₂-fertilisation of photosynthesis. However, the extent of the NPP increase varies hugely across the models from less than 10% (ACCESS-ESM1-5) to almost 150% (CanESM5), most likely due to differing levels of nutrient limitations. As a result, the change in global land carbon also varies hugely across the models, from an increase of about 200 PgC (ACCESS-ESM1-5) to an increase in excess of 1300 PgC (CanESM5), as shown in panel (d). The model projections for changes in ocean carbon storage are more clustered, but still range from a little more than 400 PgC (CNRM-ESM2-1) to over 750 PgC (BCC). In general, most models show evidence of saturation of CO₂ effects on NPP and the change in land carbon (the only exception being CanESM5). Sink saturation effects in the ocean carbon uptake are not immediately evident in panel (e). However, panel (e) does reveal a lagged (‘cold-start’) in ocean uptake, for all models. This is indicative of a non-negligible characteristic timescale for ocean carbon uptake, which hides the sink saturation effects, as we show below.

Figure 2

In order to reveal the sink saturation, we estimate the quasi-equilibrium state from Equation 5. Figure 3 shows the original transient BGC C⁴MIP runs (black lines), and this reconstructed quasi-equilibrium state (blue dots), with one panel for each of the C⁴MIP models. Now, ocean sink saturation is clearly seen in the reconstructed quasi-equilibrium of every model. We can now compare these quasi-equilibrium responses to the saturating CO₂ response of Equation 6 using the estimated values of ${\beta }_o$ and ${\kappa }_o$. These reconstructions are shown by the blue lines in Figure 3 (and we do likewise for land carbon, see Supplementary Material, Figure S.4).

Figure 3

Figure 4 shows a similar process for the climate response of land carbon. Here the climate response (shown by black dots) is diagnosed as the difference between the coupled runs (COU) and the BGC runs. We plot this residual against the temperature anomaly from the coupled run. In general, CO₂-only responses are smoother (see Figure 3 for the ocean and Figure S.4 for the land), while climate-only responses show more scatter, especially for the land (see Figure 4 for the land and Figure S.5 for the ocean). However, all models show continuous reductions in the rate of land carbon uptake due to global warming, apart from CanESM5 for which the diagnosed climate change effect is negative to start with and then switches to become positive beyond a global warming of approximately 4$\mathrm{°}$C. As for ocean carbon uptake, models show evidence of lags in their responses, which introduces a curvature in the line linking the black dots (see Supplementary Material, Figure S.5). Again, we can estimate the quasi-equilibrium response using Equation 5 and the characteristic timescales diagnosed using Equation 4. The red dots in Figure 4 show this reconstructed quasi-equilibrium response for each model, which is always more linear in global warming than the transient response (black dots). The fitted straight line for the quasi-equilibrium response, with gradient ${\gamma }_L$, is also shown. We do likewise for ${\gamma }_O$ for each model (Figure S.5).

Figure 4

Table 2 shows the best-fit parameters for each model using our new $\beta –\gamma –\kappa –\tau$ method, and Figure S.6 compares the new $\beta$ and $\gamma$ values with the old ones. It should be noted here that the $\beta$ parameters in the old and new approach are not directly comparable. In the new approach $\beta$ defines the initial sensitivity of the carbon store to CO₂, from the pre-industrial CO₂ concentration. In the old approach, $\beta$ relates to the mean impact of CO₂ on the carbon stores from pre-industrial to 2×CO₂ or 4×CO₂. As a result of sink saturation, which is neglected in the standard $\beta –\gamma$ approach, the new $\beta$ values are much larger. For the new approach, as with the old, we see more variation in the best-fit parameters for land carbon than for ocean carbon. Indeed, there is a good deal of model agreement on the ocean carbon responses to both CO₂ and climate. For both land and ocean, it is notable that the new fit implies more negative land and ocean responses to global warming (ensemble mean ${\gamma }_L$ moves from –45$\pm$51 PgC/K to –66$\pm$73 PgC/K; ensemble mean ${\gamma }_O$ moves from –17$\pm$5 PgC/K to –40$\pm$10 PgC/K). This is primarily due to fitting to diagnosed quasi-equilibrium states rather than the dynamical evolution of the carbon stores (see Figures 4 and S.5).

Figure 5 shows the fit of the new $\beta \mathrm{–}\gamma \mathrm{–}\kappa \mathrm{–}\tau$ approach (thick dashed lines) to the changes in land and ocean storage in the coupled C⁴MIP runs (thick continuous lines). For comparison, we also show the old $\beta \mathrm{–}\gamma$ fits (thin continuous lines). It is clear that the new approach picks up the sink saturation effects that are evident in the C⁴MIP outputs but which cannot be captured by the linear $\beta \mathrm{–}\gamma$ approach. This is further evidenced by Figure 6, which shows root mean square errors (RMSEs) for the old and new approaches relative to the C⁴MIP outputs. With the single exception of the change in land carbon in CanESM5 (which is unusually linear when plotted against CO₂), RMSEs are reduced significantly for all models and for both land and ocean, in the new approach. The improvement in the quality of the fit for land carbon uptake is especially impressive, with mean RMSE reduced by the new approach from about 70 PgC in the old approach, to about 17 PgC in the new.

Figure 5

Figure 6

4 Discussion

Our extension to the old $\beta –\gamma$ model rests on four key assumptions. First, we assume that the quasi-equilibrium response can be reconstructed from the transient response using a single timescale. This assumption is more justifiable for land carbon than for ocean carbon, as ocean uptake is known to respond on multiple timescales (Canadell et al., 2021; Sarmiento and Toggweiler, 1984; Strassmann and Joos, 2018). To demonstrate this, we examine how well our new fit predicts the land and ocean carbon uptake 120 years after an abrupt quadrupling of CO₂ (i.e., by examining the abrupt 4×CO₂ runs with the C⁴MIP models). Figure S.7 shows that our $\beta –\gamma –\kappa –\tau$ fit to the 1% year runs does surprisingly well in predicting the land uptake in the 4×CO₂ runs (see Figure S.7a). However, the assumption of a single timescale of the ocean response tends to overestimate carbon uptake 120 years after the abrupt CO₂ increase (see Figure S.7b).

Second, we assume that the single timescale for land is approximated by the ratio of the equilibrium land carbon to the equilibrium NPP. The quality of the fits to the 1% per year runs (see Figure S.4), and the ability to approximately reproduce the abrupt 4×CO₂ response (see Figure S.7a), suggest that this is a good assumption for the land. Thirdly, we assume that a single ocean timescale of 100 years is sufficient to identify the decadal-to-century timescale quasi-equilibrium response of the ocean. This is arguably the biggest weakness of our $\beta –\gamma –\kappa –\tau$ approach. It is an assumption that we have made in our attempt to develop a transparently simple extension to the existing $\beta –\gamma$ formulation, but one that warrants further analysis in the future.

Finally, we assume that the quasi-equilibrium response to CO₂ approximates a (Michaelis–Menton-like) hyperbolic response. The quality of the fit of the saturating hyperbolic to the diagnosed quasi-equilibrium from the CO₂-only (BGC) runs suggests that this is a good assumption at least up to 900 ppmv for the ocean (see Figure 3), and for most land models up to 4×CO₂ (see Figure S.4). The exceptions are ACCESS-ESM1.5 and IPSL-CM6A-LR, which show a loss of the quasi-equilibrium land carbon sink at high CO₂ even in the absence of climate change effects. This may be indicative of these models having characteristic timescales that are not well approximated by the ratio of the equilibrium land carbon to equilibrium NPP.

Despite these simplifications within our new $\beta –\gamma –\kappa –\tau$ approach, it reproduces the C⁴MIP 1% per year runs much better than the existing $\beta –\gamma$ approach without requiring any new runs to be analysed.

5 Conclusions

The carbon cycle and the climate system interact on timescales from seasonal to multi-millennial. Over the last 25 years, GCMs of the climate have started to include models of the land and ocean carbon cycle as interactive components. These models indicate that carbon cycle feedbacks on future climate change could be very important, but differences between models remain, especially in the response of the land carbon cycle to increases in CO₂ and changes in climate. Reducing the remaining uncertainties will require a combination of better diagnosis, improvements in model process representation, and better model calibration. Until now, the diagnosis of the sources of differences between ESMs has relied on an assumption that changes in land and ocean carbon storage respond approximately linearly to changes in atmospheric CO₂ and global warming (the ‘$\beta –\gamma$’ framework). In this paper, we have shown that the results from ESMs instead show strong evidence of sink saturation at high-CO₂, which cannot be captured by the standard $\beta –\gamma$ approach. We have therefore developed a simple extension to the existing approach to capture sink saturation and the characteristic timescales of carbon uptake by the land and ocean. The resulting $\beta –\gamma –\kappa –\tau$ provides a much improved fit to the outputs from complex coupled climate-carbon cycle GCMs. We hope that our new approach will be helpful in analysing the varying carbon cycle responses across different models, and between model generations. The representation of sink saturation effects within our approach should make it easier to identify the impact of new process representations such as nutrient limitations. In addition, it is vital to keep sink saturation in mind when considering future carbon sinks, and our new $\beta –\gamma –\kappa –\tau$ approach allows these effects to be much more clearly diagnosed.

Code availability

Code accompanying this paper can be found on GitHub at https://github.com/josephjclarke/ClimateCarbonCycleFeedbacks

Additional File

The additional file for this article can be found as follows:

Data availability

The CMIP6 Earth system model output data is available online: https://esgf-node.ipsl.upmc.fr/projects/cmip6-ipsl/ and the C⁴MIP data used in Arora et al. (2020) is available at: https://c4mip.net/results.

Author Contributions

PMC conceived of the new $\beta –\gamma –\kappa –\tau$ approach and designed the study. JJC processed the C⁴MIP models data and drafted Section 3. SJ drafted Section 1 and RMV drafted Section 2 of the paper. All authors contributed to the study and to the submitted paper.