Carbon–concentration and carbon–climate feedbacks in CMIP6 models and their comparison to CMIP5 models

Results from the fully and biogeochemically coupled simulations in which CO₂ increases at a rate of 1 % yr⁻¹ (1pctCO2) from its preindustrial value are analyzed to quantify the magnitude of carbon–concentration and carbon–climate feedback parameters which measure the response of ocean and terrestrial carbon pools to changes in atmospheric CO₂ concentration and the resulting change in global climate, respectively. The results are based on 11 comprehensive Earth system models from the most recent (sixth) Coupled Model Intercomparison Project (CMIP6) and compared with eight models from the fifth CMIP (CMIP5). The strength of the carbon–concentration feedback is of comparable magnitudes over land (mean ± standard deviation = 0.97 ± 0.40 PgC ppm⁻¹) and ocean (0.79 ± 0.07 PgC ppm⁻¹), while the carbon–climate feedback over land (−45.1 ± 50.6 PgC ^∘C⁻¹) is about 3 times larger than over ocean (−17.2 ± 5.0 PgC ^∘C⁻¹). The strength of both feedbacks is an order of magnitude more uncertain over land than over ocean as has been seen in existing studies. These values and their spread from 11 CMIP6 models have not changed significantly compared to CMIP5 models. The absolute values of feedback parameters are lower for land with models that include a representation of nitrogen cycle. The transient climate response to cumulative emissions (TCRE) from the 11 CMIP6 models considered here is 1.77 ± 0.37 ^∘C EgC⁻¹ and is similar to that found in CMIP5 models (1.63 ± 0.48 ^∘C EgC⁻¹) but with somewhat reduced model spread. The expressions for feedback parameters based on the fully and biogeochemically coupled configurations of the 1pctCO2 simulation are simplified when the small temperature change in the biogeochemically coupled simulation is ignored. Decomposition of the terms of these simplified expressions for the feedback parameters is used to gain insight into the reasons for differing responses among ocean and land carbon cycle models.

1 Introduction

The Earth system responds to the perturbation of atmospheric CO₂ concentration ([CO₂]), caused by anthropogenic emissions of CO₂ or any other forcing, via changes in its physical climate. The changes in the globally averaged temperature and the subsequent changes in other components of physical climate due to changes in radiative forcing associated with [CO₂] are larger than what would be expected from the blackbody response alone. The reason for this is that the positive feedbacks associated with various aspects of the climate system enhance the initial warming. These primarily include changes in atmospheric water vapour, tropospheric lapse rate, surface albedo resulting from ice and snow, and clouds (Hansen et al., 1984; Gregory et al., 2009; Ceppi and Gregory, 2017).

The biogeochemical cycling of carbon is also affected by changes in [CO₂] and the physical climate. In fact, changes in both the physical climate and the biogeochemical carbon cycle affect each other through multiple feedbacks. The response of the Earth's carbon cycle for both land and ocean components has been characterized in terms of carbon–concentration and carbon–climate feedback parameters which quantify their response to changes in [CO₂] and the physical climate, respectively (Friedlingstein et al., 2006; Arora et al., 2013). The carbon–concentration feedback (β) quantifies the response of the carbon cycle to changes in [CO₂] and is expressed in units of carbon uptake or release per unit change in [CO₂] (PgC ppm⁻¹). The carbon–climate feedback (γ) quantifies the response of the carbon cycle to changes in physical climate and is expressed in units of carbon uptake or release per unit change in global mean temperature (PgC ^∘C⁻¹). The changes in physical climate, in this framework, are expressed simply in terms of changes in the global mean near-surface air temperature although, of course, the carbon cycle also responds to other aspects of changes in climate (in particular precipitation over land and circulation changes in the ocean). The assumption is that the effect of other aspects of changes in climate on the carbon cycle can be broadly expressed in terms of changes in near-surface air temperature. These feedback parameters can be calculated from Earth system model (ESM) simulations globally, separately over land and ocean, regionally, or over individual grid cells which makes somewhat more sense over land than over ocean to investigate their geographical distribution (Yoshikawa et al., 2008; Boer and Arora, 2010; Tjiputra et al., 2010; Roy et al., 2011; Friedlingstein et al., 2006; Arora et al., 2013). The feedback analysis has shown that the carbon–concentration feedback is negative from the atmosphere's perspective. That is, an increase in [CO₂] leads to an increased carbon uptake by land and ocean which leads to a decrease in [CO₂], thereby slowing CO₂ accumulation in the atmosphere. The carbon–climate feedback, in contrast, has been shown to be positive in ESM simulations (at the global scale) from the atmosphere's perspective since an increase in temperature decreases the capacity of land and ocean to take up carbon, thereby contributing to a further increase in atmospheric CO₂.

The carbon–concentration and carbon–climate feedback parameters serve several purposes. First, these feedback parameters allow for the comparison of models in a simple and straightforward manner despite their underlying complexities and different model structures. Intermodel comparisons offer several benefits, including common standards and experiment protocol, coordination, and documentation that facilitate the distribution of model outputs and the characterization of the mean model response (Eyring et al., 2016), as has been shown for multiple model intercomparison projects (MIPs). Second, they allow for the quantification of the contribution of the two feedback processes to allowable anthropogenic emissions for a given CO₂ pathway. For example, Arora et al. (2013) and Gregory et al. (2009) showed that the contribution of the carbon–concentration feedback to allowable diagnosed emissions is about 4–4.5 times larger than that of the carbon–climate feedback. Third, they allow the comparison of feedbacks between climate and the carbon cycle to other feedbacks operating in the climate system as was performed by Gregory et al. (2009). Fourth, the feedback parameters can be considered as emergent properties of the coupled carbon cycle climate system which can potentially be constrained by observations, as Wenzel et al. (2014) attempted for the carbon–climate feedback parameter over land.

Here, we build on the work carried out in earlier studies that compared the strength of the carbon–concentration and carbon–climate feedbacks in coupled general circulation models with land and ocean carbon cycle components. Friedlingstein et al. (2006; hereafter F06) reported the first such results from the Coupled Climate–Carbon Cycle Model Intercomparison Project (C⁴MIP). Arora et al. (2013) (hereafter A13) compared the strength of the carbon–concentration and carbon–climate feedbacks from models participating in the fifth phase of the Coupled Model Intercomparison Project (CMIP5; https://pcmdi.llnl.gov/mips/cmip5/, last access: 15 December 2019; Taylor et al., 2012). The A13 study found that the strength of the two feedbacks was weaker and the spread between models was smaller in their study than in F06. However, the results from these two studies are not directly comparable because of several reasons. The results from the F06 study were based on the Special Report on Emissions Scenarios (SRES) A2 emissions scenario, while those in the A13 study were based on the 1 % yr⁻¹ increasing CO₂ experiment in which the atmospheric CO₂ concentration increases from its preindustrial value of around 285 ppm until it quadruples over a 140-year period (referred to as the 1pctCO2 experiment in the framework of the Coupled Model Intercomparison Project, CMIP). The absolute values of the feedback parameters are known to be dependent on the state of the system, the timescale of forcing (i.e. underlying emissions and concentration scenario), and the approach used to calculate them (Plattner et al., 2008; Zickfeld et al., 2011; Hajima et al., 2014; Gregory et al., 2009; Boer and Arora, 2010). The varying approaches employed over the past decade have made the cross comparison of feedbacks among the studies and different generations of Earth system models difficult.

In order to address the diversity of approaches to diagnose climate carbon cycle feedbacks and to promote a robust standard moving forward the C⁴MIP community has endorsed a framework of tiered experiments (Jones et al., 2016) that builds upon the core preindustrial control and 1pctCO2 experiments performed as part of the CMIP DECK (Diagnostic, Evaluation and Characterization of Klima) experiments (Eyring et al., 2016). Here, we compare carbon–concentration and carbon–climate feedbacks from models participating in the C⁴MIP (Jones et al., 2016) contribution to the sixth phase of CMIP (CMIP6; Eyring et al., 2016). To maintain continuity and consistency, feedback parameters are derived from the 1pctCO2 experiments as they were in A13. The 1pctCO2 experiment is a DECK experiment in the CMIP6 framework. All participating modelling groups are expected to perform DECK experiments to help document basic characteristics of models across different phases of CMIP (Eyring et al., 2016).

2 Feedbacks metrics in the coupled climate–carbon system

We largely follow the climate carbon cycle feedbacks framework presented in A13 (which in turn was built on F06) but with some additional modifications that are explained below. Only the primary equations are presented here, while the bulk of the framework is summarized in the Appendix for completeness. We also provide some history of how the carbon feedbacks analysis reached its current stage.

Carbon feedbacks analysis is traditionally based on simulations run with fully, radiatively, and biogeochemically coupled model configurations of an Earth system model. The objective of these simulations is to isolate feedbacks discussed above. In a biogeochemically coupled simulation (referred to here as the BGC simulation), biogeochemical processes over land and ocean respond to increasing atmospheric CO₂, while the radiative transfer calculations in the atmosphere use a CO₂ concentration that remains at its preindustrial value. Small climatic changes occur in the BGC simulation due to changes in evaporative (or latent heat) flux resulting from stomatal closure over land (associated with increasing [CO₂]), changes in vegetation structure, and changes in vegetation coverage and composition (in models which dynamically simulate competition between their plant functional types, PFTs), all of which affect latent and sensible heat fluxes at the land surface. In a radiatively coupled simulation (referred to here as the RAD simulation) increasing atmospheric CO₂ affects the radiative transfer processes in the atmosphere and hence climate but not the biogeochemical processes directly over land and ocean, for which the preindustrial value of atmospheric CO₂ concentration is prescribed. In a fully coupled simulation (referred to here as the COU simulation) both the biogeochemical and the radiative processes respond to increasing CO₂.

Following the F06 methodology which uses time-integrated fluxes (which are the same as the changes in carbon pool sizes), the changes in land (L) or ocean (O) carbon pools (${\mathrm{\Delta }\mathrm{C}}_X,X=\text{L,O}$) can be expressed using three equations corresponding to the BGC, RAD, and COU experiments, as shown in Eq. (1) (see also the Appendix).

$$\begin{array}{}\text{(1a)}& {\displaystyle }\begin{array}{rl}& \text{Radiatively coupled simulation}\\ & {\mathrm{\Delta }\mathrm{C}}_X^{+}=\int F_X^{+}\mathrm{d}t={\mathit{\gamma }}_X{T}^{+}\end{array}\end{array}$$

$$\begin{array}{}\text{(1b)}& {\displaystyle }\begin{array}{rl}& \text{Biogeochemically coupled simulation}\\ & {\mathrm{\Delta }\mathrm{C}}_X^{\ast }=\int F_X^{\ast }\mathrm{d}t={\mathit{\beta }}_X{c}^{\prime }+{\mathit{\gamma }}_X{T}^{\ast }\end{array}\end{array}$$

$$\begin{array}{}\text{(1c)}& {\displaystyle }\begin{array}{rl}& \text{Fully coupled simulation}\\ & {\mathrm{\Delta }{\mathrm{C}}^{\prime }}_X=\int F_X^{\prime }\mathrm{d}t={\mathit{\beta }}_X{c}^{\prime }+{\mathit{\gamma }}_X{T}^{\prime }\end{array}\end{array}$$

Here F⁺, F^∗, and F^′ are the CO₂ flux changes (PgC yr⁻¹); ${\mathrm{\Delta }\mathrm{C}}_X^{+}$, ${\mathrm{\Delta }\mathrm{C}}_X^{\ast }$, and ${\mathrm{\Delta }{\mathrm{C}}^{\prime }}_X$ the changes in global carbon pools (PgC); T⁺, T^∗, and T^′ the temperature changes (^∘C) in the RAD, BGC, and COU simulations, respectively; and the subscript X=LO refers to either the land or ocean model components. c^′ is the change in [CO₂]. Here and elsewhere uppercase C is used to denote pools and lowercase c is used to denote atmospheric CO₂ concentration, [CO₂]. All changes are defined relative to a preindustrial equilibrium state represented by the preindustrial control simulation. In the context of a specified-concentration simulation (the 1pctCO2 experiment in our case), c^′ is the same in BGC and COU simulations. There is no β_Xc^′ term in the RAD simulation since the biogeochemistry sees the preindustrial value of [CO₂] although T⁺ is a function of increasing c^′ that is seen only by the radiative transfer calculations.

These equations assume linearization of the globally integrated surface–atmosphere CO₂ flux (for land and ocean components) in terms of global mean temperature and [CO₂] change (compared to a preindustrial control run) and serve to define the carbon–concentration (β_X) and carbon–climate (γ_X) feedback parameters. A similar set of equations can be written that defines the instantaneous values of the feedback parameters and is based on fluxes rather than their time-integrated values (see Eqs. A4 and A5 in the Appendix). Both the time-integrated flux and the instantaneous flux-based versions of the feedback parameters evolve over time in an experiment as shown in A13.

There are several different ways in which the feedbacks (β_X and γ_X) in a coupled climate and carbon cycle system may be evaluated: (1) the experiments may use specified (concentration-driven) or freely evolving (emissions-driven) [CO₂], (2) any two of the three configurations of an experiment (COU, RAD, and BGC) may be used to calculate the two feedback parameters, and (3) the experiment may be based on an idealized scenario (like the 1pctCO2 experiment) or a more realistic emissions scenario. In addition, the small temperature change in the BGC simulation, T^∗, may be ignored, and other external forcings such as nitrogen (N) deposition or land use change, which directly affect carbon fluxes, may or may not be taken into account. The original framework proposed by F06 used COU and BGC versions (referred to as coupled and uncoupled in the F06 study) of an emissions-driven simulation for the SRES A2 scenario. The F06 framework assumed that the small temperature change in the BGC simulation can be ignored. A13 used BGC and RAD versions of the 1pctCO2 experiment in which the evolution of [CO₂] is specified and took into account the small global mean temperature change in the BGC simulation.

With regard to the use of concentration-driven versus emissions-driven simulations, Gregory et al. (2009) recommended the use of specified concentration simulations, which ensures consistency of [CO₂] across models, and this recommendation has been adopted since CMIP5. C⁴MIP has also adopted the use of the 1pctCO2 simulation; i.e. an idealized scenario is preferred over a more realistic scenario. The 1pctCO2 experiment provides an ideal experiment to compare carbon–climate interactions across models as the experiment does not include the confounding effects of other climate forcings (including land use change, non-CO₂ greenhouse gases, and aerosols) and is a CMIP DECK experiment, as mentioned earlier.

Table 1 The values of the carbon–concentration (β) and carbon–climate (γ) feedback parameters can be solved using results from any two combinations of the RAD, BGC, and COU versions of an experiment as shown in Eq. (1). In addition, when using results from the BGC and COU simulations, the effect of temperature change in the BGC simulation (T^∗) can be neglected, as it was in the F06 study, yielding approximate values for β_X and γ_X.

Using Eq. (1) as an example, Table 1 shows how any two combinations of the three configurations of an experiment can be used to calculate the values of the two feedback parameters. The A13 study showed that under the assumption of a linear system and if the conditions $F^{\prime }=F^{+}+F^{\ast }$and $T^{\prime }=T^{+}+T^{\ast }$ are met, i.e. if the sum of flux and temperature changes in the RAD and BGC simulations is the same as that in the COU simulation, then all approaches yield exactly the same solution. However, this is not the case because of the nonlinearities involved (Gregory et al., 2009; Zickfeld et al., 2011; Schwinger et al., 2014).

The use of BGC and RAD simulations that have only biogeochemistry or radiative forcing responding to increases in [CO₂] to find the feedback parameters is attractive since these simulations were designed to isolate the feedbacks. In the RAD simulation (whose purpose is to quantify the carbon–climate feedback, γ_X) the preindustrial global carbon pools for both land and ocean typically decrease in response to an increase in global temperature (hence the positive carbon–climate feedback and the negative value of γ_X). Consequently, negative values of γ_X (positive carbon–climate feedback) are obtained when using the RAD–BGC and RAD–COU approaches (see Table 1). If, however, γ_X is determined using the BGC and COU simulations in both of which the globally summed carbon pools for land and ocean are increasing in response to increasing [CO₂], the calculated value of γ_X is different than that obtained using the RAD–BGC and RAD–COU approaches. In the ocean, the RAD simulation mainly measures the loss of near-surface carbon owing to warming of the surface ocean layer (Schwinger et al., 2014). The RAD simulation misses the suppression of carbon drawdown to the deep ocean due to weakening ocean circulation, because there is no buildup of a strong carbon gradient from the surface to the deep ocean in contrast to the BGC and COU simulations. Therefore, the absolute value of γ_X for the ocean is smaller (less negative) when calculated using the RAD–BGC and RAD–COU approaches compared to the BGC–COU approach (Schwinger et al., 2014). Over land, in the RAD simulation carbon is lost in response to increasing temperatures primarily due to an increase in heterotrophic respiration. However, an increase in temperature also potentially increases photosynthesis at high latitudes, and this increase compensates for carbon lost due to increased heterotrophic respiratory losses, especially in the presence of continuously increasing [CO₂] seen in the COU configuration. Therefore γ_X value for land calculated using RAD–BGC and RAD–COU approaches may be higher or lower than that calculated using the BGC–COU approach. These are some mechanisms that lead to nonlinearities. Since the ongoing climate change (predominantly caused by increasing [CO₂]) is best characterized by the COU simulation, it can be argued that feedback parameters are more representative when calculated using the BGC–COU approach. Here, we propose to use the COU and BGC configurations of an experiment as the standard set from which to calculate the feedback parameters as recommended in the C⁴MIP protocol (Jones et al., 2016). However, we also quantify the values of feedback parameters when using the RAD simulation for comparison. The calculated values of the carbon–concentration feedback parameter (β_X), in contrast, are less sensitive to the approach used as shown in A13.

There is no broad consensus on whether temperature change in the BGC simulation should be assumed to be zero ($T^{\ast }=\mathrm{0}$) as standard practice when calculating the strengths of the feedbacks, as it was in F06. While the globally averaged value of T^∗ is an order of magnitude smaller than T^′, the spatial pattern of T^∗ is quite different from that of T^′. The spatial pattern of temperature change in the COU simulation (T^′) is dominated by radiative forcing of increased [CO₂] with greater warming at high latitudes and over land than over ocean. In contrast, the spatial pattern of temperature change in the BGC simulation (T^∗) is determined primarily by the reduction in latent heat flux associated with stomatal closure as [CO₂] increases which reduces transpiration from vegetation (Bounoua et al., 1999; Ainsworth and Long, 2005). This process leads to a much more spatially variable pattern of temperature change (than T^′) and the associated changes in precipitation patterns due to soil moisture–atmosphere feedbacks (Chadwick et al., 2017; Skinner et al., 2017). The difference in spatial patterns of temperature and precipitation change in the RAD versus the COU simulation is another reason that the values of the carbon–climate feedback (γ_X) depend on the simulation used, and this is another pathway through which nonlinearities can occur. A complete analysis of the effect of differences in spatial patterns of climate change and the carbon state on the calculated value of γ_X when using the RAD versus the COU simulation and determining whether or not the assumption of $T^{\ast }=\mathrm{0}$ should be a standard practice are beyond the scope of this study but remain topics for additional scientific investigation. In the interim, we report here values of β_X and γ_X both by considering T^∗ and by ignoring it (i.e. $T^{\ast }=\mathrm{0}$) when using the BGC–COU approach.

Following Table 1, when using results from the BGC and the COU configurations of a specified-concentration experiment, the values of the feedback parameters are written as

$$\begin{array}{}\text{(2)}& {\displaystyle }{\mathit{\beta }}_X={\displaystyle \frac{\mathrm{1}}{c^{\prime }}}\left({\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_X^{\ast }{T}^{\prime }-{\mathrm{\Delta }{\mathrm{C}}^{\prime }}_X{T}^{\ast }}{T^{\prime }-T^{\ast }}}\right),\text{(3)}& {\displaystyle }{\mathit{\gamma }}_X={\displaystyle \frac{{\mathrm{\Delta }{\mathrm{C}}^{\prime }}_X-{\mathrm{\Delta }\mathrm{C}}_X^{\ast }}{T^{\prime }-T^{\ast }}}.\end{array}$$

Equations (2) and (3) may be rearranged to explicitly calculate the effect of the $T^{\ast }=\mathrm{0}$ assumption on calculated values of feedback parameters, as shown in Eqs. (4) and (5). Here, the T^∗ term is retained only in the second part of the equations, whose contribution becomes zero when T^∗ is ignored.

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\beta }}_X={\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_X^{\ast }}{c^{\prime }}}+{\displaystyle \frac{\mathrm{1}}{c^{\prime }}}\left[{\displaystyle \frac{\left({\mathrm{\Delta }{\mathrm{C}}^{\prime }}_X-{\mathrm{\Delta }\mathrm{C}}_X^{\ast }\right)T^{\ast }}{\left(T^{\prime }-T^{\ast }\right)}}\right]\text{(5)}& {\displaystyle }{\mathit{\gamma }}_X={\displaystyle \frac{{\mathrm{\Delta }{\mathrm{C}}^{\prime }}_X-{\mathrm{\Delta }\mathrm{C}}_X^{\ast }}{T^{\prime }}}+{\displaystyle \frac{\left({\mathrm{\Delta }{\mathrm{C}}^{\prime }}_X-{\mathrm{\Delta }\mathrm{C}}_X^{\ast }\right)T^{\ast }}{T^{\prime }\left(T^{\prime }-T^{\ast }\right)}}\end{array}$$

Finally, in regard to other external forcings such as nitrogen (N) deposition that directly affect carbon fluxes, the C⁴MIP protocol for CMIP6 (Jones et al., 2016) recommended performing additional simulations for BGC and COU versions of the 1pctCO2 experiment with time-varying N deposition in addition to their standard versions which keep N deposition rates at their preindustrial level. Simulations with N deposition can only be performed for models that explicitly model the N cycle and its interactions with the carbon (C) cycle. The rationale for recommending increasing N deposition, in conjunction with temperature and CO₂ increases, is to be able to quantify the response of feedback parameters to this third forcing. However, here we restrict ourselves to the traditional analysis that considers the climate and CO₂ forcings only. We do highlight, however, which models include coupled C–N cycle interactions over land. The analysis of runs with N deposition forcing is left for future studies.

2.1 Reasons for differences in feedback parameters among models

As shown later in this paper, the contribution of the second term involving T^∗ in expressions for the carbon–concentration (β_X) and carbon–climate (γ_X) feedback parameters (in Eqs. 4 and 5, when using the BGC–COU approach) is around 1 % to 5 %. This allows for the reasons for differences in the feedback parameters to be investigated across models as the expressions for the feedback parameters can be simplified in terms of the changes in the sizes of carbon pools (${\mathrm{\Delta }{\mathrm{C}}^{\prime }}_X$ and ${\mathrm{\Delta }\mathrm{C}}_X^{\ast }$), the temperature change in the COU simulation (T^′), and the specified change in [CO₂] (c^′) as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\beta }}_X\approx {\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_X^{\ast }}{c^{\prime }}},\text{(7)}& {\displaystyle }{\mathit{\gamma }}_X\approx {\displaystyle \frac{{\mathrm{\Delta }{\mathrm{C}}^{\prime }}_X-{\mathrm{\Delta }\mathrm{C}}_X^{\ast }}{T^{\prime }}}.\end{array}$$

2.1.1 Land

Over land, Eqs. (6) and (7) can be expanded to investigate, firstly, the contributions, from changes in live vegetation pools (ΔC_V) and dead litter plus soil carbon pools (ΔC_S) to the strength of the feedback parameters, since ${\mathrm{\Delta }\mathrm{C}}_{\mathrm{L}}={\mathrm{\Delta }\mathrm{C}}_{\mathrm{V}}+{\mathrm{\Delta }\mathrm{C}}_{\mathrm{S}}$. Secondly, Eq. (6) can be further decomposed to gain insight into the reasons for differences across models, in a manner similar to Hajima et al. (2014).

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}& {\mathit{\beta }}_{\mathrm{L}}\approx {\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\mathrm{L}}^{\ast }}{c^{\prime }}}={\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\mathrm{V}}^{\ast }+{\mathrm{\Delta }\mathrm{C}}_{\mathrm{S}}^{\ast }}{c^{\prime }}}=\left({\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\mathrm{V}}^{\ast }}{\mathrm{\Delta }{\text{NPP}}^{\ast }}}{\displaystyle \frac{\mathrm{\Delta }{\text{NPP}}^{\ast }}{\mathrm{\Delta }{\text{GPP}}^{\ast }}}{\displaystyle \frac{\mathrm{\Delta }{\text{GPP}}^{\ast }}{c^{\prime }}}\right)\\ & +\left({\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\mathrm{S}}^{\ast }}{\mathrm{\Delta }{R}_{\mathrm{h}}^{\ast }}}{\displaystyle \frac{\mathrm{\Delta }{R}_{\mathrm{h}}^{\ast }}{{\mathrm{\Delta }\mathrm{LF}}^{\ast }}}{\displaystyle \frac{{\mathrm{\Delta }\mathrm{LF}}^{\ast }}{c^{\prime }}}\right)={\mathit{\tau }}_{\text{veg}\mathrm{\Delta }}{\text{CUE}}_{\mathrm{\Delta }}{\displaystyle \frac{\mathrm{\Delta }{\text{GPP}}^{\ast }}{c^{\prime }}}\\ & +{\mathit{\tau }}_{\mathrm{soil}\mathrm{\Delta }}{\displaystyle \frac{\mathrm{\Delta }{R}_{\mathrm{h}}^{\ast }}{{\mathrm{\Delta }\mathrm{LF}}^{\ast }}}{\displaystyle \frac{{\mathrm{\Delta }\text{LF}}^{\ast }}{c^{\prime }}}\end{array}\text{(9)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{L}}\approx {\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\mathrm{L}}^{\prime }-{\mathrm{\Delta }\mathrm{C}}_{\mathrm{L}}^{\ast }}{T^{\prime }}}={\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\mathrm{V}}^{\prime }-{\mathrm{\Delta }\mathrm{C}}_{\mathrm{V}}^{\ast }}{T^{\prime }}}+{\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\mathrm{S}}^{\prime }-{\mathrm{\Delta }\mathrm{C}}_{\mathrm{S}}^{\ast }}{T^{\prime }}}\end{array}$$

The superscript ^∗ in Eq. (8) implies that the terms are calculated here using the BGC version of the 1pctCO2 experiment. In Eq. (8), ΔNPP and ΔGPP represent the change in net primary productivity (NPP) and gross primary productivity (GPP), ΔLF the change in litterfall flux, and ΔR_h the change in heterotrophic respiration, compared to the preindustrial control experiment. The multiplicative terms in Eq. (8) do indeed have physical meaning although they are based on change in the magnitude of quantities as opposed to their absolute magnitudes. We note here explicitly that as such, these terms cannot be compared directly to the terms which are based on absolute magnitudes.

The term $\frac{\mathrm{\Delta }\mathrm{NPP}}{\mathrm{\Delta }\mathrm{GPP}}$ (fraction) is the fraction of GPP (above its preindustrial value) that is turned into NPP after autotrophic respiratory losses are taken into account. We use the term carbon use efficiency (CUE) but subscripted by Δ (CUE_Δ) to represent $\frac{\mathrm{\Delta }\mathrm{NPP}}{\mathrm{\Delta }\mathrm{GPP}}$. The subscripted Δ allows CUE_Δ to be differentiated from CUE as used in the existing literature (Choudhury, 2000), which represents the fraction of absolute GPP that is converted to NPP rather than its change over some time period, as well as the point that we consider globally integrated rather than locally derived quantities. Similarly, the term $\frac{{\mathrm{\Delta }\mathrm{C}}_{\mathrm{V}}}{\mathrm{\Delta }\mathrm{NPP}}$ represents a measure of turnover or residence timescale of carbon in the vegetation pool (τ_vegΔ, years). The term $\frac{\mathrm{\Delta }\mathrm{GPP}}{c^{\prime }}$ (PgC yr⁻¹ ppm⁻¹) is a measure of the strength of the globally integrated CO₂ fertilization effect. However, in the models that dynamically simulate changes in vegetation cover, the effect of changes in vegetation coverage is implicitly included in this term. The term $\frac{{\mathrm{\Delta }\mathrm{C}}_{\mathrm{S}}}{\mathrm{\Delta }{R}_{\mathrm{h}}}$ is a measure of the average residence time of carbon in the dead litter and soil carbon pools (τ_soilΔ, years). However, as with CUE, this quantity cannot be compared directly to the residence time of carbon in the litter plus soil carbon pool calculated using the absolute values of C_S and R_h. Nor can it be compared to the changes in carbon residence time due to the “false priming effect” associated with the increase in NPP inputs, as [CO₂] increases, into the dead carbon pools (Koven et al., 2015). $\frac{\mathrm{\Delta }{R}_{\mathrm{h}}}{\mathrm{\Delta }\mathrm{LF}}$ (fraction) is a measure of the increase in heterotrophic respiration per unit increase in litterfall rate, and $\frac{\mathrm{\Delta }\mathrm{LF}}{c^{\prime }}$ (PgC yr⁻¹ ppm⁻¹) indicates the global increase in litterfall rate per unit increase in CO₂, which in principle, should be close to the change in net primary productivity per unit increase in CO₂, $\left({\mathrm{CUE}}_{\mathrm{\Delta }}\frac{\mathrm{\Delta }\mathrm{GPP}}{c^{\prime }}\right)$. A comparison of these terms across models can potentially yield insight into the reasons for large differences in land carbon uptake across models.

2.1.2 Ocean

Assuming changes in the biological organic carbon inventory are small, the change in the ocean carbon inventory, ΔC_O, is defined by an integral of the change in the dissolved inorganic carbon, ΔDIC, and density over the ocean volume:

$$\begin{array}{}\text{(10)}& \mathrm{\Delta }{\mathrm{C}}_{\mathrm{O}}=\mathrm{12}\mathrm{gC}{\mathrm{mol}}^{-\mathrm{1}}\underset{V}{\int }\mathrm{\Delta }\mathrm{DIC}\mathrm{d}V\times {\mathrm{10}}^{-\mathrm{15}},\end{array}$$

where ΔC_O is in petagrams of carbon; the ocean dissolved inorganic carbon, DIC, is in moles per cubic meter; the ocean volume V is in cubic meters; and the multiplier 10⁻¹⁵ converts grams to petagrams of carbon.

To gain insight into how the ocean carbon distribution is controlled, the ocean dissolved inorganic carbon, DIC, may be defined in terms of separate carbon pools (Ito and Follows, 2005; Williams and Follows, 2011; Lauderdale et al., 2013; Schwinger and Tjiputra, 2018):

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\text{DIC}=& {\text{DIC}}_{\text{preformed}}+{\text{DIC}}_{\text{regenerated}}\\ =& {\text{DIC}}_{\text{sat}}+{\text{DIC}}_{\text{disequilib}}+{\text{DIC}}_{\text{regenerated}},\end{array}\end{array}$$

where the preformed carbon, DIC_preformed, is the amount of carbon in a water parcel when in the mixed layer at the time of subduction and the regenerated carbon, DIC_regenerated, is the amount of dissolved inorganic carbon accumulated below the mixed layer due to the biological regeneration of organic carbon. The preformed carbon is affected by the carbonate chemistry and ocean physics. To gain further insight into how close the ocean is to an equilibrium with the atmosphere, the preformed carbon, DIC_preformed, is further split into saturated, DIC_sat, and disequilibrium, DIC_disequilib, components. The saturated component represents the concentration in surface water fully equilibrated with the contemporary atmospheric CO₂ concentration. The disequilibrium component represents the extent that surface water is incompletely equilibrated before subduction, which is affected by the strength of the ocean circulation altering the residence time in the mixed layer and the ocean ventilation rate. Each of these components is affected by the increase in atmospheric CO₂ and the changes in climate.

The change in the global ocean carbon inventory, ΔC_O, relative to the preindustrial period may then be related to the global volume integral of the change in each of these DIC pools:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\mathrm{\Delta }\mathrm{C}}_{\mathrm{O}}=& {\mathrm{\Delta }\mathrm{C}}_{\text{preformed}}+{\mathrm{\Delta }\mathrm{C}}_{\text{regenerated}}\\ =& {\mathrm{\Delta }\mathrm{C}}_{\text{sat}}+{\mathrm{\Delta }\mathrm{C}}_{\text{disequilib}}+{\mathrm{\Delta }\mathrm{C}}_{\text{regenerated}},\end{array}\end{array}$$

where ΔC_preformed is the preformed carbon inventory, ΔC_sat is the saturated carbon inventory, ΔC_disequilib is the disequilibrium carbon inventory, and ΔC_regenerated is the regenerated carbon inventory.

The simplified expressions for carbon cycle feedback parameters in Eqs. (6) and (7) based on the air–sea flux changes to the ocean may then be approximated by the global ocean carbon inventory changes, which may be expressed in terms of these different global ocean carbon pools (Williams et al., 2019):

$$\begin{array}{}\text{(13)}& {\displaystyle }\begin{array}{rl}{\mathit{\beta }}_{\mathrm{O}}\approx & {\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\mathrm{O}}^{\ast }}{c^{\prime }}}={\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\text{preformed}}^{\ast }}{c^{\prime }}}+{\displaystyle \frac{\mathrm{\Delta }{C}_{\text{regenerated}}^{\ast }}{c^{\prime }}}\\ =& {\displaystyle \frac{\mathrm{\Delta }{\mathrm{C}}_{\text{sat}}^{\ast }}{c^{\prime }}}+{\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\text{disequilib}}^{\ast }}{c^{\prime }}}+{\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\text{regenerated}}^{\ast }}{c^{\prime }}},\end{array}\text{(14)}& {\displaystyle }\begin{array}{rl}{\mathit{\gamma }}_{\mathrm{O}}\approx & {\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\mathrm{O}}^{\prime }-{\mathrm{\Delta }\mathrm{C}}_{\mathrm{O}}^{\ast }}{T^{\prime }}}={\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\text{preformed}}^{\prime }-{\mathrm{\Delta }\mathrm{C}}_{\text{preformed}}^{\ast }}{T^{\prime }}}\\ +& {\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\text{regenerated}}^{\prime }-{\mathrm{\Delta }\mathrm{C}}_{\text{regenerated}}^{\ast }}{T^{\prime }}}\\ =& {\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\text{sat}}^{\prime }-{\mathrm{\Delta }\mathrm{C}}_{\text{sat}}^{\ast }}{T^{\prime }}}+{\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\text{disequilib}}^{\prime }-{\mathrm{\Delta }\mathrm{C}}_{\text{disequilib}}^{\ast }}{T^{\prime }}}\\ +& {\displaystyle \frac{{\mathrm{\Delta }\mathrm{C}}_{\text{regenerated}}^{\prime }-{\mathrm{\Delta }\mathrm{C}}_{\text{regenerated}}^{\ast }}{T^{\prime }}}.\end{array}\end{array}$$

The anomalies for each of these carbon pools are calculated as

$$\begin{array}{}\text{(15)}& {\displaystyle }\begin{array}{rl}& {\mathrm{\Delta }\mathrm{DIC}}_{\text{regenerated}}=-R_{\mathrm{CO}}\mathrm{\Delta }\mathrm{AOU}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{\Delta }\mathrm{Alk}\right.\\ & -{\mathrm{\Delta }\mathrm{Alk}}_{\text{pre}}-R_{\mathrm{NO}}\mathrm{\Delta }\mathrm{AOU}\left.\right),\end{array}\text{(16)}& {\displaystyle }\begin{array}{rl}& {\mathrm{\Delta }\mathrm{DIC}}_{\text{sat}}=f{\left(p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}},T_{\mathrm{o}},S_{\mathrm{o}},\text{P},\text{Si},{\text{Alk}}_{\text{pre}}\right)}_t\\ & -f{\left(p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}},T_{\mathrm{o}},S_{\mathrm{o}},\mathrm{P},\mathrm{Si},{\text{Alk}}_{\text{pre}}\right)}_{t=\mathrm{0}},\end{array}\text{(17)}& {\displaystyle }{\mathrm{\Delta }\mathrm{DIC}}_{\text{disequilib}}=\mathrm{\Delta }\mathrm{DIC}-{\mathrm{\Delta }\mathrm{DIC}}_{\text{regenerated}}-{\mathrm{\Delta }\mathrm{DIC}}_{\text{sat}},\end{array}$$

where R_CO and R_NO are constant stoichiometric ratios, ΔAOU is the change in apparent oxygen utilization from its preindustrial value (where preformed oxygen is assumed to be approximately saturated with respect to atmospheric oxygen); ΔAlk is the change in alkalinity; T_o and S_o are the ocean temperature and salinity, respectively; P and Si are the phosphate and silicate concentrations; and ΔAlk_pre is the change in preformed alkalinity (Ito and Follows, 2005; Williams and Follows, 2011; Appendix of Lauderdale et al., 2013). In Eq. (16), ΔDIC_sat is calculated using the partial pressure of CO₂ in the atmosphere (p${\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}$) and preformed alkalinity as represented by the function f following the iterative solution for the ocean carbon system of Follows et al. (2006) and by considering the small contribution of minor species (borate, phosphate, silicate) to the preformed alkalinity, at time t, and the preindustrial values, at time t=0. In Eq. (16), the calculation uses the preformed alkalinity, the alkalinity at the time of water subduction, instead of the total instantaneous alkalinity to remove the effect of CaCO₃ dissolution from the time the water parcel lost contact with the atmosphere. The preformed alkalinity is estimated from a multiple linear regression using salinity and the conservative tracer PO ($\mathrm{PO}={\mathrm{O}}_{\mathrm{2}}-{\mathrm{R}}_{\mathrm{o}\mathrm{2}:\mathrm{P}}\mathrm{P}$; Gruber et al., 1996), with the coefficients of this regression estimated based on the upper ocean (first 10 m) alkalinity, salinity, oxygen, and phosphate in each model. Our diagnostics of the ocean feedbacks and carbon pools depend primarily upon changes in DIC and upon the preformed and regenerated pools, relative to the preindustrial period, although differences in the preindustrial ocean in our suite of models do affect the saturated DIC changes relative to the preindustrial period by ∼5 % or less due to the nonlinearity of the carbonate chemistry.

Table 2 Primary features of the physical atmosphere and ocean components and land and ocean carbon cycle components of the 11 participating models in this study.

3 Model descriptions

Table 2 summarizes the primary features of the 11 comprehensive ESMs that contributed results to this study. Brief descriptions of land and ocean carbon cycle components of these ESMs are provided in the Appendix. The 11 ESMs, in alphabetical order, are the (1) Commonwealth Scientific and Industrial Research Organisation (CSIRO) ACCESS-ESM1.5, (2) Beijing Climate Center (BCC) BCC-CSM2-MR, (3) Canadian Centre for Climate Modelling and Analysis (CCCma) CanESM5, (4) Community Earth System Model version 2 (CESM2), (5) Centre National de Recherches Météorologiques (CNRM) CNRM-ESM2-1, (6) Institut Pierre Simon Laplace (IPSL) IPSL-CM6A-LR, (7) Japan Agency for Marine-Earth Science and Technology (JAMSTEC) in collaboration with the University of Tokyo and the National Institute for Environmental Studies (team MIROC) MIROC-ES2L, (8) Max Planck Institute for Meteorology (MPI) MPI-ESM1.2-LR, (9) Geophysical Fluid Dynamics Laboratory (GFDL) NOAA-GFDL-ESM4, (10) Norwegian Climate Centre (NCC) NorESM2-LM, and (11) United Kingdom (UK) UKESM1-0-LL.

In contrast to the A13 study where only two of the eight participating comprehensive ESMs had the terrestrial N cycle implemented and coupled to their C cycle, in this study 6 of the 11 participating ESMs represent coupling of terrestrial C and N cycles. These six models are the ACCESS-ESM1.5, CESM2, MIROC-ES2L, MPI-ESM1.2-LR, NorESM2-LM, and UKESM1-0-LL. Note that CESM2 and NorESM2-LM employ the same land surface component – version 5 of the Community Land Model (CLM5) – so we expect the land carbon cycle to respond very similarly in the two models. Three of the ESMs have land components that dynamically simulate vegetation cover and competition between their PFTs – NOAA-GFDL-ESM4, MPI-ESM1.2-LR, and UKESM1-0-LL.

Figure 1 Temperature changes in the fully, biogeochemically, and radiatively coupled configurations of the 1pctCO2 experiment across participating CMIP6 (a) and CMIP5 (b) comprehensive ESMs that participated in this and the A13 study, respectively. The model mean is indicated by the solid lines, and the range across the models is indicated by shading around the solid lines. Individual CMIP6 model results are shown in Fig. A1 in the Appendix.

4 Results

4.1 Global surface CO₂ fluxes and temperature change

Figure 1 shows the simulated changes in temperature in the three model configurations (COU, BGC, and RAD) of the 1pctCO2 experiment. Here and in subsequent figures, results are also shown for the eight comprehensive ESMs that participated in the A13 study. The eight models in the A13 study are a subset of 11 models considered in this study although they have been updated since CMIP5.

As expected, temperature change is higher in the COU and RAD simulations than in the BGC simulation, since the radiative forcing responds to increasing [CO₂] in these simulations. The small temperature change in the BGC simulation is due to a number of not only contributing but also compensating factors: (1) reduction in transpiration and hence latent heat flux due to stomatal closure in response to increasing [CO₂] (Cao et al., 2010); (2) increase in vegetation leaf area index (LAI), which decreases land surface albedo and hence increases absorbed solar radiation; and (3) increase in vegetation fraction in models that explicitly simulate competition between their plant functional types (PFTs) over land (NOAA-GFDL-ESM4, MPI-ESM1.2-LR, and UKESM1-0-LL), which also leads to reduced land surface albedo. As a result, temperature change in the COU simulation is higher than in the RAD simulation since these biogeochemical processes are active and contribute to a small additional warming. This is seen in Fig. 1a for CMIP6 models and Fig. 1b for CMIP5 models.

When comparing CMIP5 and CMIP6 models, the CMIP6 models are on average slightly warmer than CMIP5 models in the COU and RAD simulations. In Fig. 1a, the globally averaged near-surface temperature change at CO₂ quadrupling in the COU simulation is 4.87 ^∘C in CMIP6 models, compared to 4.74 ^∘C in CMIP5 models. The CMIP6 ensemble considered here includes some high-climate-sensitivity models including CanESM5, CESM2, CNRM-ESM2-1, IPSL-CM6A-LR, and UKESM1-0-LL. The globally averaged temperature change at CO₂ quadrupling in the COU simulation for the eight models that are common to this (CMIP6) and the A13 (CMIP5) studies is 4.97 and 4.74 ^∘C, respectively. The temperature change in the BGC simulation in CMIP6 models (0.21 ^∘C) is, however, slightly lower than in the CMIP5 models (0.26 ^∘C). The values in Fig. 1 for participating CMIP5 models are slightly different than those reported in the A13 study because those numbers also included the UVic Earth System Climate Model (an intermediate complexity model), which we have omitted here to keep the comparison consistent between comprehensive ESMs. In addition, in contrast to A13, the temperature at the end of the simulations in this study is calculated after fitting a fourth-order polynomial in R to the model mean values rather than using the actual model mean value at the end of the simulation which can be higher or lower than that calculated using the polynomial fit due to interannual variability. A fourth-order polynomial fit has been shown to yield a good estimate of the forced response of the global mean temperature response and to minimize the potential influence of internal variability (Hawkins and Sutton, 2009).

Figure 2 Model mean values and the range across models for annual simulated atmosphere–land CO₂ flux (a, b) and their cumulative values (c, d) for participating CMIP6 (a, c) and CMIP5 (b, d) models from the fully, biogeochemically, and radiatively coupled versions of the 1pctCO2 experiment. Individual CMIP6 model results are shown in Fig. A1.

Figures 2 and 3 show simulated model mean values and the range across models for annual simulated atmosphere–land and atmosphere–ocean CO₂ fluxes and their cumulative values for participating CMIP6 and CMIP5 models from the COU, BGC, and RAD configurations of the 1pctCO2 experiment. The general results from CMIP6 models are broadly similar in nature to those from CMIP5 models, as would be expected, with higher annual and cumulative values of atmosphere–land and atmosphere–ocean CO₂ fluxes in the BGC simulation compared to the COU simulation in which the radiative warming caused by increasing CO₂ weakens the land and ocean sinks. In the RAD simulation, where land and ocean carbon cycle components do not respond to increasing [CO₂], both components lose carbon, for reasons discussed below.

Figure 3 Model mean values and the range across models for annual simulated atmosphere–ocean CO₂ flux (a, b) and their cumulative values (c, d) for participating CMIP6 (a, c) and CMIP5 (b, d) models from the fully, biogeochemically, and radiatively coupled versions of the 1pctCO2 experiment. Individual CMIP6 model results are shown in Fig. A1.

Over land, the model mean rate of increase in atmosphere–land CO₂ flux declines and even becomes negative in the COU and BGC simulations as the terrestrial CO₂ fertilization effect saturates and the carbon pools build up, which increases the respiratory losses. The biggest difference between the CMIP5 and CMIP6 models is that the cumulative land carbon uptake in the COU simulation is about 25 % higher in CMIP6 (635 ± 258 PgC, mean ± standard deviation) models than in CMIP5 (505 ± 297 PgC) models, although this increase is not statistically significant across the model ensemble (Mann–Whitney test). Here and hereafter, we use the sample (not population) standard deviation. The cumulative value of carbon loss in the RAD simulation is similar in both CMIP6 and CMIP5 models: 239 ± 120 vs. 252 ± 158 PgC, respectively. This carbon loss occurs due to both increased heterotrophic respiration per unit carbon mass and reduced GPP (and consequently NPP) in the RAD simulation (not shown). While NPP declines globally in response to increases in temperature, mid- to high-latitude net primary production increases (Qian et al., 2010), so the reduction in global NPP comes largely from the reduction in the tropics. The large spread across CMIP6 land carbon cycle models, seen also in earlier F06 and A13 studies, has not changed significantly compared to CMIP5 models, and its implications will be discussed in more detail in Sect. 5. As discussed later in Sect. 4.3, the standard deviation of land carbon–climate feedback increases from CMIP5 to CMIP6 models, while it decreases somewhat for the land carbon–concentration feedback.

Over the ocean, the response to increasing [CO₂] and changing climate remains fairly similar across CMIP5 and CMIP6 models. The cumulative ocean carbon uptake in the COU simulation is 593 ± 54 and 611 ± 50 PgC in CMIP6 and CMIP5 models, respectively. Unlike the land uptake, however, the ocean carbon uptake does not saturate over the length of the simulation in the BGC simulation (Fig. 3a, b); it keeps on increasing albeit at a declining rate. The cumulative ocean carbon loss in the RAD simulation is 23 ± 19 and 37 ± 17 PgC in CMIP6 and CMIP5 models, respectively, and is primarily associated with warmer temperatures which reduce CO₂ solubility (Goodwin and Lenton, 2009; Schwinger et al., 2014).

As in F06 and A13, the range in cumulative atmosphere–land CO₂ fluxes among models at the end of the COU simulation, in response to changes in atmospheric CO₂ concentration and surface temperature, is 3 to 4 times larger than for the atmosphere–ocean CO₂ fluxes. Figure A1 shows results from individual CMIP6 models for which model means and ranges were shown in Figs. 1, 2, and 3 and allows for the identification of models which behave differently compared to the majority of models.

Figure 4 Components of the carbon budget terms in cumulative emissions from the 11 participating CMIP6 models based on Eq. (A6) in (a) and Eq. (A7) in (b) using results from the fully coupled 1pctCO2 simulation. The models are arranged in an ascending order based on their cumulative emissions values. Results from participating CMIP5 models in the A13 study are shown in panels (c) and (d). In addition, ESMs whose land component includes a representation of the N cycle are identified by a red font colour for cumulative land carbon uptake (a, c) and fractional emissions taken up by land (b, d). Model mean is shown not only for all models but also separately for models whose land components include or do not include a representation of the N cycle.

4.2 Carbon budget terms

Figure 4a shows the carbon budget components of the diagnosed cumulative fossil fuel emissions at the end of the 140-year period of the 1pctCO2 COU experiment when CO₂ concentration quadruples (${\stackrel{\mathrm{̃}}{E}}_{\mathrm{4}\times {\mathrm{CO}}_{\mathrm{2}}}$ or simply $\stackrel{\mathrm{̃}}{E}$), from CMIP6 models. Cumulative emissions can similarly also be calculated at 2×CO₂ (${\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$). The term “carbon budget” in this context refers to the accounting of carbon internal to individual ESMs. The sum of ocean ($\mathrm{\Delta }{\mathrm{C}}_{\mathrm{O}}^{\prime }$) and land ($\mathrm{\Delta }{\mathrm{C}}_{\mathrm{L}}^{\prime }$) sinks and the resulting change in atmospheric carbon burden ($\mathrm{\Delta }{\mathrm{C}}_{\mathrm{A}}^{\prime }$) yield cumulative fossil fuel emissions which are consistent with the specified CO₂ pathway (the 1pctCO2 scenario in this case) as indicated in the Appendix (Eq. A6). The corollary to this is that, in a specified emissions simulation, if the respective fossil fuel emissions were to be used in their models, each model would yield CO₂ concentrations that rise at a rate of 1 % yr⁻¹. The term “diagnosed” implies that the cumulative fossil fuel emissions are calculated from changes in atmosphere, land, and ocean carbon pools in the specified-concentration 1pctCO2 experiment. Figure 4b shows the terms of the budgets as fractional components for atmosphere (A), land (L), and ocean (O) based on Eq. (A7), where f_A is the airborne fraction of emissions and f_L and f_O are the fractions of emissions take up by land and ocean, respectively. More details are provided in Sect. A1 of the Appendix.

$$\begin{array}{}\text{(18)}& {\displaystyle }{\mathrm{\Delta }{\mathrm{C}}^{\prime }}_{\mathrm{A}}+{\mathrm{\Delta }{\mathrm{C}}^{\prime }}_{\mathrm{L}}+{\mathrm{\Delta }{\mathrm{C}}^{\prime }}_{\mathrm{O}}={\int }_{\mathrm{0}}^{t}E\mathrm{d}t=\stackrel{\mathrm{̃}}{E}\text{(19)}& {\displaystyle }{f}_{\mathrm{A}}+f_{\mathrm{L}}+f_{\mathrm{O}}=\mathrm{1}\end{array}$$

All panels in Fig. 4 identify models whose land component includes a representation of the N cycle – the cumulative land carbon uptake (panels a and c) and fractional emissions taken up by land (panels b and d) for these models are shown in red.

Figure 5 Carbon–concentration (a) and carbon–climate (b) feedback parameters over land from participating CMIP6 models calculated using the approaches summarized in Table 1. The boxes show the mean ± 1 standard deviation range, and the individual coloured dots represent individual models. Models which include a representation of the land nitrogen cycle are identified with a circle around their dot. The model mean ± 1 standard deviation range of feedback parameters is also separately shown for models which do and do not represent the land nitrogen cycle using the BGC–COU approach. Results from participating CMIP5 models in the A13 study are shown in (c) and (d).

Consistent with Figs. 2 and 3, and CMIP5 results reported in the A13 study, the differences among models are primarily due to the diverse response of the land carbon cycle components. While the model mean cumulative carbon uptake by the ocean is fairly similar between participating CMIP5 (611 ± 50 PgC) and CMIP6 (593 ± 54 PgC) models, the land uptake is higher in CMIP6 (635 ± 258 PgC) compared to CMIP5 (505 ± 297 PgC) models, as mentioned earlier. This is the case even when the CanESM5, the model with the largest land carbon uptake, is omitted from CMIP6 models (model mean land carbon uptake for the remaining 10 models is 578 ± 185 PgC). As a result, the model mean cumulative diagnosed emissions from CMIP6 models (3031 ± 242 PgC) are about 4 % higher than those from CMIP5 models (2927 ± 294 PgC). In Fig. 5a, the land carbon uptakes in the CESM2 (656 PgC) and NorESM2-LM (652 PgC) models are very similar; as noted above these models employ the same land component. Model mean estimates that are reported separately for models whose land component does and does not include a representation of the N cycle, for both CMIP5 and CMIP6 models, show that the model mean land carbon uptake is lower for models that explicitly represent the N cycle. As a consequence, the airborne fraction of emissions is also higher for models that represent the land N cycle and their diagnosed cumulative fossil fuel emissions are lower (Fig. 4). Of the 11 CMIP6 ESMs considered in this study, 6 represent the N cycle over land compared to only 2 of the 8 considered in the A13 study based on CMIP5 models. Yet, the model mean land carbon uptake over land is higher in this study than in the A13 study. This is partly because of the three models with the largest land carbon uptake (CNRM-ESM2-1, BCC-CSM2-MR, and CanESM5) which do not include land the N cycle (Fig. 4a). In addition, inclusion of the N cycle does not universally imply lower land C uptake. In Fig. 4a, IPSL-CM6A-LR and NOAA-GFDL-ESM4, both of which do not include the land N cycle, yield lower land carbon uptake than four of the models that do include the land N cycle.

Figure 4a and c allow for the direct comparison of models from the same modelling group. CanESM2, from the CCCma, which had below-average land carbon uptake among CMIP5 models, has evolved into CanESM5, a model with the largest land carbon uptake among CMIP6 models. The reason for this is an increase in the strength of its CO₂ fertilization effect following the retuning of its photosynthesis downregulation parameters, using carbon budget constraints over the historical period, as explained in Arora and Scinocca (2016). CESM1, which had one of the lowest land carbon uptakes among CMIP5 models, because of its apparently excessive nitrogen limitation effect in CLM4, has evolved into CESM2 (with the CLM5 land component) with near-average land carbon uptake among CMIP6 models. The transition of CLM from CLM4 to CLM5 on the one hand and the reduction in its nutrient constraints on photosynthesis and the parametric controls on fertilization responses on the other are discussed in Wieder et al. (2019) and Fisher et al. (2019), respectively. The land carbon uptake in MIROC-ESM increased from the lowest among CMIP5 models to near average for MIROC-ES2L among CMIP6 models, due to a new terrestrial biogeochemical component (Ito and Inatomi, 2012). Although the CO₂ fertilization effect in this new land model is weaker likely due to the incorporation of the nitrogen cycle, the model yields relatively higher NPP (Hajima et al., 2020), due to a higher CUE_Δ (as confirmed later in Sect. 4.4.1). The land carbon uptake in the IPSL-CM5A-LR model decreased from being the second largest in CMIP5 models to below average for the IPSL-CM6A-LR model due to the implementation of terrestrial photosynthesis downregulation, as a function of CO₂ concentration, which leads to a decrease in GPP across all latitudes, with the largest decrease in the tropics. For the MPI ESM, the decrease in land carbon uptake in MPI-ESM-LR for CMIP5 and MPI-ESM1.2-LR for CMIP6 is associated with the implementation of a nitrogen cycle model (Goll et al., 2017) and a new soil carbon model Yasso (Goll et al., 2015). Compared to its predecessor HadGEM2-ES, UKESM1-0-LL represents a prognostic treatment of terrestrial nitrogen including its impact on carbon storage in vegetation biomass and soil organic matter. A limitation on terrestrial productivity from available nitrogen is likely also the main reason for reduced land carbon storage in UKESM1-0-LL compared to in HadGEM2-ES.

The ocean carbon uptake in the IPSL model decreased from being the largest among CMIP5 models in IPSL-CM5A-LR to being lower than average for IPSL-CM6A-LR, and this change is attributed to a greater ocean stratification in the IPSL-CM6A-LR. The annual mean mixed-layer depth is 46.7 and 40.2 m in IPSL-CM5A-LR and IPSL-CM6A-LR, respectively. While NorESM1-ME was one of the CMIP5 models with the largest ocean carbon uptake, NorESM2-LM has an ocean carbon uptake close to the CMIP6 model mean. This change is a consequence of changes in the simulated (shallower depth and weaker strength) Atlantic meridional overturning circulation and reduced mixed-layer biases particularly at high latitudes (less deep winter mixing). Due to these modifications, the efficiency of carbon export below the mixed layer in NorESM2-LM is considerably reduced compared to in NorESM1-ME. This, in turn, leads to less excess carbon stored in the North Atlantic Deep Water (below 2000 m) as well as in the Antarctic Intermediate Water.

Figure A2 in the Appendix shows a version of Fig. 4 but at the time of CO₂ doubling (at year 70). Interestingly, the ordering of the models according to their diagnosed cumulative emissions at 2×CO₂ is different from that at 4×CO₂. As expected, however, the model mean fractional emissions taken up by land and ocean at 2×CO₂ are higher than at 4×CO₂, because both land and ocean carbon sinks relatively weaken as CO₂ continues to increase.

Figure 6 Carbon–concentration (a) and carbon–climate (b) feedback parameters over ocean from participating CMIP6 models calculated using the approaches summarized in Table 1. The boxes show the mean ± 1 standard deviation range. Results from participating CMIP5 models in the A13 study are shown in (c) and (d).

4.3 Feedback parameters

Figure 5a and b compare the carbon–concentration feedback (β_L) and carbon–climate feedback (γ_L) parameters over land from participating CMIP6 models, calculated using results at the end of the BGC, RAD, and COU simulations. The plots show not only feedback parameters from different models as coloured dots but also their mean ± 1 standard deviation as a box. Three primary observations can be made from Fig. 5. First and foremost, the spread in the magnitude of the carbon–concentration and carbon–climate feedback over land in CMIP6 models is of a similar magnitude to that in CMIP5 models (Fig. 5c and d). Second, the carbon–climate feedback (γ_L) is more sensitive to the approach used (and hence the type of simulations used) to derive its value than the carbon–concentration feedback (β_L). The absolute value of β_L varies by around 7 %, while γ_L varies by up to 26 %, depending on the approach used. Third, in the model mean sense, the absolute strength of the feedback parameters is weaker for models that include a representation of the N cycle, for both CMIP5 and CMIP6 models. Both the carbon gain due to an increase in atmospheric CO₂ concentration and the carbon loss due to an increase in the global average temperature in models with representation of the land N cycle are much lower than in models that do not include the N cycle. This response is most likely explained by the N limitation of photosynthesis as CO₂ increases and the additional release of N from dead organic matter as warming increases, which boosts productivity, thereby compensating for carbon lost due to increased respiratory losses, as also discussed in A13. The values of the feedback parameters, however, overlap between models that do and do not include a representation of the N cycle, given the wider spread in the feedback parameter values among models that do not include a representation of the land N cycle compared to models that do.

Figure 6a and b compare the carbon–concentration feedback (β_O) and carbon–climate feedback (γ_O) parameters over the ocean from participating CMIP6 models. For both CMIP5 and CMIP6 models, the absolute spread in the magnitude of the feedback parameters across the participating models is an order of magnitude smaller for the ocean C cycle component compared to the land C cycle component, as was also seen in F06 and A13. Similar to the land, the calculated values of the ocean carbon–climate feedback (γ_O) are more sensitive to the approach used (and hence the type of simulations used) than the ocean carbon–concentration feedback (β_O). In agreement with Schwinger et al. (2014), the absolute values of γ_O are 2–3 times larger when calculated using the COU and BGC simulations, compared to cases when RAD simulation is used, for reasons mentioned earlier. Figures 5 and 6 show also that while the strength of the carbon–concentration feedback is similar over land and ocean, the strength of the carbon–climate feedback parameter over ocean is much weaker than over land.

Figures 5 and 6 provide justification for using the BGC–COU approach, over the RAD–BGC and RAD–COU approaches, in calculating the feedback parameters as discussed below. In Fig. 6, the absolute magnitude of γ_O when using the BGC–COU approach is about twice as large in CMIP5 models and more than 3 times larger in CMIP6 models compared to its model mean value calculated using the RAD–BGC and RAD–COU approaches. The reason for this is that the RAD simulation misses the suppression (due to weakening of the ocean circulation) of carbon drawdown to the deep ocean due to a lack of buildup of a strong carbon gradient from the atmosphere to the deep ocean, as mentioned earlier. This process is important when climate change is forced by increasing atmospheric CO₂, and therefore feedback parameters calculated using the BGC–COU approach are more likely to include all processes relevant to application to realistic scenarios. In Fig. 6, the value of γ_O changes sign for the CNRM-ESM2-1 model from positive when calculated using the RAD–BGC or RAD–COU approaches to negative when calculated using the BGC–COU approach, and this further illustrates the sensitivity of feedback parameters to the approach used to calculate them. Section A2 discusses the reasons for this sensitivity in the CNRM-ESM2-1 model.

In Fig. 5, although the carbon–climate feedback parameter over land (γ_L) is larger in absolute terms, it is comparatively less sensitive to the approach used than that over ocean, because over land an increase in temperature not only increases the respiratory losses but also affects photosynthetic processes especially in conjunction with increasing CO₂. Warmer temperatures increase photosynthesis over mid- to high-latitude regions where photosynthesis is currently limited by temperature and more so with increasing CO₂ but decrease photosynthesis over tropical regions where the temperatures are already too warm for optimal photosynthesis. The net result of these compensating processes plays out very differently in different models, and in the model mean sense this results in less sensitivity over land than over ocean of the calculated value of the carbon–climate feedback parameter (γ_L and γ_O, respectively) to the different approaches. This is seen in both CMIP5 and CMIP6 models. Over land, photosynthesis is also affected by temperature (with widely varying responses between models) in addition to respiration, and the γ_L values vary widely between models between the RAD–BGC and RAD–COU approach and the BGC–COU approach. This is seen, for example, for the ACCESS-ESM1.5, IPSL-CM6A-LR, and CanESM5 models in Fig. 5b.

Figures 5 and 6 also show that the effect of assuming T^∗ (the temperature change in the BGC simulation) to be zero is around 1 % for the calculated value of the carbon–concentration feedback parameter (${\mathit{\beta }}_X;X=\text{LO}$) and around 5 % for the carbon–climate feedback parameter (${\mathit{\gamma }}_X;X=\text{LO}$). This small effect of T^∗ on the calculated global values of the feedback parameter allows for the investigation of the reasons for differences among model by using simplified forms of β_X and γ_X as presented in Eqs. (6) and (7).

For completeness, Table A1 in the Appendix summarizes the values of feedback parameters for both land and ocean from CMIP6 and CMIP5 models (corresponding to Figs. 5 and 6) not only at 4×CO₂ but also at 2×CO₂. Table A1 also shows the value of parameter α, the linear transient climate sensitivity to CO₂, following F06 (their Eq. 6), which is calculated using values at the end of the COU simulation as

$$\begin{array}{}\text{(20)}& T^{\prime }=\mathit{\alpha }{c}^{\prime }.\end{array}$$

Figure 7 Carbon uptake over land in the BGC simulation, used to calculate land carbon–concentration feedback (β_L) and its partitioning into vegetation and soil + litter carbon pools across the participating CMIP6 models (a). Panel (b) shows the fractional land carbon uptake by vegetation and soil + litter carbon pools in the BGC simulation. No partitioning is shown for the BCC-CSM2-MR model because total land carbon uptake in this model exceeded the sum of changes in the vegetation and soil + litter carbon pools by more than 10 %. Total land carbon uptake in models which include a representation of the N cycle is shown in red. The results from the BCC-CSM2-MR model are not used in calculating the model mean values.

4.4 Reasons for differences among models

4.4.1 Land

Equations (8) and (9) in Sect. 2.1.1 are used to gain insight into reasons for differing responses of land models. In the BGC–COU approach and assuming $T^{\ast }=\mathrm{0}$ (Eq. 8), the carbon uptake in the BGC simulation is used to calculate the carbon–concentration feedback parameter (β_L). Figure 7 shows how this carbon uptake over land is separated into vegetation and soil + litter components both in absolute (panel a) and fractional (panel b) terms. Figure 7b shows that models vary widely in terms of how the carbon uptake over land is split into vegetation and soil + litter components. The model mean values indicate that slightly more of the carbon sequestered is allocated to vegetation (55 %) than to the soil + litter pools (45 %).

Figure 8 Individual terms of Eq. (8) which contribute to changes in vegetation (ΔC_V) and litter + soil (ΔC_S) carbon pools. Values from the BCC-CSM2-MR model are not used in calculating the model mean.

Figure 8 shows the individual components of Eq. (8) which contribute to terms corresponding to changes in vegetation (ΔC_V) and soil + litter (ΔC_S) carbon pools. Figure 7a is repeated in Fig. 8 for the easy matching of individual components of Eq. (8) with their corresponding models. The model mean values of individual terms do not take into account the results from the BCC-CSM2-MR model as explained in the figure caption. In essence, the terms in Fig. 8 are emergent properties of the land models of the individual ESMs and result from their multiple interacting processes. The comparison of the individual terms of Eq. (8) provides additional insight into the reasons for differences in land models. For example, the CNRM-ESM2-1 model has the highest land carbon uptake among all models in the BGC simulation. However, this is not caused by a strong CO₂ fertilization effect (the $\frac{\mathrm{\Delta }\mathrm{GPP}}{{}^{c^{\prime }}}$ term) but rather by the relatively high τ_vegΔ and τ_soilΔ values. The CO₂ fertilization effect is strongest for the three models that simulate vegetation cover dynamically (NOAA-GFDL-ESM4, MPI-ESM1.2-LR, and UKESM1-0-LL) since the $\frac{\mathrm{\Delta }\mathrm{GPP}}{c^{\prime }}$ term also implicitly includes the effect of increasing vegetation cover as CO₂ increases. The tree cover in the NOAA-GFDL-ESM4 model, for example, increases in the BGC simulation – particularly in dry, high-latitude regions above 50^∘ N (not shown). However, these models do not simulate the largest land carbon uptake because of their lower-than-average τ_vegΔ and τ_soilΔ values. The $\frac{\mathrm{\Delta }\mathrm{GPP}}{c^{\prime }}$ term is unable to capture the CO₂ fertilization effect separately from increasing vegetation cover, and this illustrates the challenge in comparing models that do and do not simulate vegetation cover dynamically. The CanESM5 model exhibits higher-than-average land carbon uptake despite its near-average strength of the CO₂ fertilization effect and τ_vegΔ and τ_soilΔ values. However, its CUE_Δ is the highest, and therefore a much larger fraction of GPP is converted to NPP. Although CUE_Δ is not the same as CUE, we found that CUE_Δ and CUE (calculated at the end of the 1pctCO2 simulation at 4×CO₂) are strongly correlated with a correlation of around 0.90 (not shown). Similarly, τ_vegΔ is strongly correlated, with a correlation of 0.96, to τ_veg = C_V ∕ NPP calculated at the end of the simulation. The ACCESS-ESM1.5 model exhibits the lowest land carbon uptake because of its weak CO₂ fertilization effect and the lowest CUE_Δ of all models. Finally, the $\frac{\mathrm{\Delta }{R}_{\mathrm{h}}}{\mathrm{\Delta }\mathrm{LF}}$ term shows the least variability across models, which is reflective of the fact that the magnitude of the heterotrophic respiration flux is dominated by NPP inputs into the dead carbon pools (Koven et al., 2015). Several of these individual terms are strongly correlated. The $\frac{\mathrm{\Delta }\mathrm{GPP}}{c^{\prime }}$ and $\frac{\mathrm{\Delta }\mathrm{LF}}{c^{\prime }}$ terms have a correlation of 0.77, and CUE_Δ $\frac{\mathrm{\Delta }\mathrm{GPP}}{c^{\prime }}$ and $\frac{\mathrm{\Delta }\mathrm{LF}}{c^{\prime }}$ have a correlation of 0.94, since a stronger CO₂ fertilization effect also implies a larger litterfall flux per unit CO₂. Surprisingly, CUE_Δ and τ_vegΔ are negatively correlated (correlation = −0.49) across models, indicating that models which retain a higher fraction of GPP as NPP typically get rid of vegetation carbon sooner via litterfall as indicated by a faster turnover of vegetation (lower τ_vegΔ), thereby partially compensating for higher CUE_Δ.

Figure 9 The changes in vegetation and soil + litter carbon pools in the COU relative to the BGC simulation, as shown in Eq. (9), which contribute to the calculation of the carbon–climate feedback over land (γ_L) in the BGC–COU approach. The names of models which include the N cycle are shown in a red font colour.

Figure 9 investigates the reasons for varying magnitudes of the carbon–climate feedback over land (γ_L). In Eq. (9), γ_L is a function of change in land carbon (divided into vegetation and soil + litter components) in the COU relative to the BGC simulation and the temperature change in the COU simulation (T^′). Over land, the higher temperatures in the COU relative to the BGC simulation affect not only both autotrophic and heterotrophic respiratory fluxes, from live and dead vegetation pools, respectively, but also gross photosynthesis rates. The primary effect of this temperature change in COU versus the BGC simulation is the loss of carbon from the soil + litter carbon pool (hence the negative sign of γ_L for most models; Fig. 6b and d), but changes in the vegetation carbon pool also occur. Although γ_L also depends on T^′, Fig. 9 arranges models in order from the largest to smallest loss of land carbon in COU relative to the BGC simulation to illustrate the varying response of the models. This ordering of models changes slightly if the carbon loss (or gain in the CanESM5 model) is divided by the temperature change T^′ in the COU simulation (yielding the value of γ_L which assumes $T^{\ast }=\mathrm{0}$ as in Eq. 9).

As shown in Fig. 9, all models lose carbon from the soil + litter carbon pool but with widely varying magnitudes. Although typically smaller than the change in the soil + litter carbon pool, the change in the vegetation carbon pool in the COU relative to the BGC simulation is not of the same sign across models. Of the 11 participating models, 6 lose carbon in the vegetation pool in the COU relative to the BGC simulation, thereby contributing to increasing the absolute magnitude of γ_L, while the remaining 5 exhibit an increase in the vegetation carbon pool, thereby decreasing the absolute magnitude of γ_L. The largest increase in the vegetation carbon pool is seen in the CanESM5 model that more than compensates for the carbon loss from the soil + litter carbon pool, yielding a positive value of γ_L in contrast to other models. This case is one of the few times a positive value of γ_L is seen in an Earth system model. Thornton et al. (2009) reported positive γ_L after their first attempt to include N cycle in the CLM. Preliminary analysis of CanESM5 data shows the increase in vegetation carbon in the COU relative to the BGC simulation is caused by the increase in GPP and the resulting vegetation growth at middle to high latitudes in response to warming temperatures and increasing CO₂. Interestingly, this response is not seen at 2×CO₂ (see Table A1 in the Appendix), and γ_L is still negative for CanESM5.

Table 3 Correlation between carbon–concentration (β_X) and carbon–climate (γ_X) feedback parameters over land and ocean across comprehensive ESMs from the CMIP5 intercomparison in the A13 study and CMIP6 intercomparison in this study. For the land correlation it is also shown when CanESM5 is excluded from CMIP6 models.

The loss in land carbon in the COU relative to the BGC simulation (except for the CanESM5 model that gains carbon), indicated by the dark orange bars in Fig. 9, is strongly correlated with the carbon gain in the BGC simulation (Fig. A1e; correlation is 0.59 for all models and 0.87 when CanESM5 is excluded) but not with the absolute amount of total land carbon. Figure A3 in the Appendix shows the absolute amount of carbon in soil + litter and vegetation pools, and their change from the beginning, for the BGC simulation. The models vary widely in terms of the absolute size of the carbon pools, especially for the soil + litter pool. There are two implications of models losing more carbon in the COU relative to the BGC simulation when they take up more carbon in the BGC simulation alone. First, the transient behaviour of a model is determined primarily by its response to CO₂ and temperature perturbations and not by the absolute amount of land carbon. Second, carbon–concentration (β_X) and carbon–climate (γ_X) feedback parameters must be correlated as well. Indeed, this is not only the case over land for both CMIP5 and CMIP6 models but also true for ocean feedbacks although the correlations are somewhat weaker over the ocean. These correlations are shown in Table 3 and are negative since higher positive values of β_X are correlated with higher negative values of γ_X, indicating that models that take up more carbon with increasing CO₂ also release more carbon when they “see” the associated higher temperatures.

4.4.2 Ocean

The time-integrated air–sea flux of carbon provides the dominant contribution to the increase in global ocean carbon through changes in the DIC inventory. However, the global ocean carbon inventory is also affected by the land-to-ocean carbon flux from river runoff and the carbon burial in ocean sediments (see Table A2 in the Appendix).

Ocean carbon cycle feedbacks are defined in terms of ocean carbon inventory changes for the COU simulation and by the differences in COU relative to the BGC simulation. To fully understand the ocean carbon cycle feedbacks, it is necessary to understand the ocean carbon distributions for the preindustrial period and then analyze the carbon anomalies relative to the preindustrial period for these climate model experiments.

Figure 10 Meridional section of the dissolved inorganic carbon, DIC (mol m⁻³), and constituent carbon pools in UK-ESM1-0-LL for the zonally averaged Atlantic and Southern Ocean: (a) the preindustrial absolute concentrations and (b–d) the anomalies relative to the preindustrial state at year 140 for (b) the COU configuration, (c) the BGC configuration, and (d) the COU minus the BGC configuration. The DIC is separated into saturated carbon, DIC_sat, the disequilibrium carbon, DIC_disequilib, and the regenerated carbon, DIC_regenerated. The Atlantic and Southern Ocean domains are separated by a vertical black line.

Ocean carbon distribution

The ocean dissolved inorganic carbon, DIC, distribution is controlled by a combination of physical, chemical, and biological processes. For the preindustrial period, there is less DIC in the warmer waters of the upper ocean and more DIC in the colder mid-depth and bottom waters (Figs. 10a, 11a; illustrated here for UKESM1-0-LL as a representative example, and Figs. S1 to S7 show similar distributions for all the diagnosed Earth system models). The vertical extent of the low DIC follows the undulations of the thermocline, which is defined by strong vertical temperature and density gradients, and is deeper over the subtropical gyres at 30^∘ N and 30^∘ S and shallower in the equatorial zone and at high latitudes. The greater DIC at depth is a consequence of greater solubility in colder waters and the accumulation of DIC from the regeneration of organic matter.

Figure 11 Meridional section of the dissolved inorganic carbon, DIC (mol m⁻³), and constituent carbon pools in UK-ESM1-0-LL for the zonally averaged Pacific and Southern Ocean: (a) the preindustrial absolute concentrations and (b–d) the anomalies relative to the preindustrial state at year 140 for (b) the COU configuration, (c) the BGC configuration, and (d) the COU minus the BGC configuration. The DIC is separated into saturated carbon, DIC_sat, the disequilibrium carbon, DIC_disequilib, and the regenerated carbon, DIC_regenerated. The Pacific and Southern Ocean domains are separated by a vertical black line.

To gain insight into how the ocean carbon distribution is controlled, the DIC is separated into three pools, DIC_sat, DIC_disequilib, and DIC_regenerated, as defined earlier. The DIC distribution for both the preindustrial period and after 140 years in the 1pctCO2 simulation reveal the following key features for each of these carbon pools (Figs. 10a, b and 11a, b):

• The saturated carbon pool provides the dominant contribution to the DIC, holding more than 2.15 mol C m⁻³, particularly within cooler waters below the thermocline.

• The regenerated carbon pool enhances the carbon stored below the surface waters, typically providing an additional 0.2 mol C m⁻³ within the Southern Ocean and older waters spreading from the Southern Ocean into the Atlantic and below the thermocline in the Pacific.

• The disequilibrium carbon is small close to the surface, representing waters close to an equilibrium with the atmosphere. There is sometimes a positive disequilibrium of up to 0.05 mol C m⁻³ in some surface waters, which is associated with upwelling transferring carbon-rich deeper waters to the surface. The disequilibrium carbon is more strongly negative below the thermocline, typically reaching −0.1 mol C m⁻³ in the Atlantic and −0.02 mol C m⁻³ in the Southern Ocean and Pacific. In the preindustrial period, the undersaturation in carbon below the thermocline is due to the subduction of cold waters at high latitudes that have not equilibrated fully with the atmosphere, which then spread by advection along density surfaces. In the model integrations reaching year 140, the carbon below the thermocline becomes further undersaturated relative to the contemporary atmosphere due to the rapid rise in [CO₂].

Next we consider the anomalies in the DIC at year 140 in the COU configurations of the 1pctCO2 simulation calculated relative to the preindustrial period. The carbon anomaly, ΔDIC, in the COU configuration is positive over the upper thermocline over the Atlantic and Pacific basins, reaching +0.3 mol C m⁻³, coinciding with regions that are well ventilated. This gain in carbon is made up of an increase in the saturated carbon over all depths due to higher atmospheric CO₂. There is a dipole in the disequilibrium anomaly (Figs. 10b, c and 11b, c); it is generally weakly positive in the upper ocean and more strongly negative in deeper waters below the thermocline reaching up to −0.2 mol C m⁻³. This negative disequilibrium anomaly in deeper waters is smallest in the relatively well-ventilated mid-depth waters of the North Atlantic but extends over nearly all of the more poorly ventilated mid-depth waters of the Pacific (Figs. 10b and 11b).

The regenerated carbon anomaly is relatively small in magnitude, reaching less than 0.05 mol C m⁻³, and varies regionally, enhanced within much of the North Atlantic and the thermocline of the Pacific but with little change in the deep waters of the Pacific (Figs. 10b and 11b). The increase in regenerated carbon is due to a weakening of ocean overturning leading to an increase in residence time and an associated accumulation of DIC from the regeneration of biologically cycled carbon (Bernardello et al., 2014; Schwinger et al., 2014). The regenerated carbon signal does not change in the mid-depths and deep Pacific as 140 years is too short an integration timescale for any effect to be detected.

To diagnose the carbon cycle feedback parameters, the ocean carbon response needs to be considered for the BGC configuration where there is only limited warming from the increase in atmospheric CO₂ and therefore limited change in climate and ocean circulation. The resulting DIC anomalies are generally very similar to those for the COU configuration (Figs. 10b, c and 11b, c), which is to be expected as the dominant effect for the ocean carbon response is the enhanced ocean uptake of carbon in response to the increase in [CO₂]. There is a weakening in ventilation in the COU configuration due to the additional radiative forcing. In comparison, in the BGC configuration, there is no change in the circulation as there is no radiative warming effect, so there is slightly more carbon uptake in the northern North Atlantic, such as revealed at around 50^∘ N, compared with the COU configuration. For the BGC configuration, the saturated carbon pool is slightly greater at depth due to the water masses being cooler than in the COU configuration, the disequilibrium anomaly shows a less negative anomaly in the northern North Atlantic because there is little or no change in ventilation, and there are only slight differences in the regenerated pool.

The climate response to rising [CO₂] is now considered in terms of the difference in the COU and BGC configurations, which includes the combined effects of warming and circulation changes (Figs. 10d and 11d). The surface warming drives a decrease in solubility, an increase in stratification, and a reduction in ventilation, which leads to an overall decrease in carbon uptake over the Southern Ocean and Pacific basins and much of the Atlantic basin. There is a decrease in the saturated carbon pool associated with the warming acting to inhibit carbon uptake. The regenerated carbon anomaly is enhanced in the deep northern North Atlantic and in the Southern Ocean. The regenerated carbon anomaly for this climate response is very similar to that for the COU configuration, suggesting that the regenerated carbon anomaly is mainly due to circulation changes: the gain in the regenerated carbon anomaly is consistent with the expected longer residence time from weaker overturning and ventilation. There is a more negative disequilibrium anomaly in the deep waters of the North Atlantic, which is a consequence of weaker ventilation.

To gain more insight into the disequilibrium response, the ocean DIC response is also considered for the radiatively coupled integration (RAD), where the ocean biogeochemistry does not see the increase in [CO₂]. The warming in the RAD simulation leads to a weakening in the overturning, which enhances the residence time in the surface waters and so generally decreases the magnitude of the disequilibrium anomaly in the North Atlantic (Fig. S8), making the disequilibrium less negative relative to the preindustrial period and so forming a positive disequilibrium anomaly at year 140. In comparison the COU–BGC approach captures the effect of the warming under rising [CO₂], leading to the disequilibrium anomaly instead becoming more negative at depth, since the weakening in the ventilation leads to more of the anthropogenic carbon remaining at the surface rather than being transferred into the deeper ocean (Schwinger et al., 2014).

Figure 12 Carbon uptake over the ocean in the biogeochemically coupled simulation, used to calculate the ocean carbon–concentration feedback and its partitioning into saturated, disequilibrium, and regenerated carbon pools across the participating CMIP6 models (left) using Eq. (12). No partitioning is shown for models for which 3D ocean fields were not available, and the results of these models are not used in calculating the model mean values (right). The sum of the partitions does not exactly match the total ocean uptake diagnosed from the air–sea fluxes due to land–ocean interactions involving storage in sediments and river inputs.

Changes in ocean carbon pools for diagnosing feedback parameters

The ocean carbon–concentration feedback parameter, β_O, is diagnosed from the changes in the ocean carbon inventories for the BGC configuration, which does not include radiative warming due to increasing [CO₂] (Eq. 13). There is a consistent increase in ocean carbon storage across all models with a model mean value of around 670 PgC (Fig. 12, turquoise bars). This increase in ocean carbon storage is made up of an increase in the saturated carbon inventory, ΔC_sat, by about 3100 PgC from the increase in [CO₂] (Fig. 12, red bars). This increase is partly offset by a more negative disequilibrium carbon, ΔC_disequilib, of typically −2500 PgC (Fig. 12, blue bars), representing how the ocean carbon uptake cannot keep up with the rate of [CO₂] increase. There is relatively little change in the regenerated carbon inventory, ΔC_regenerated. The resulting β_O is positive and mainly explained by the chemical response involving the rise in ocean saturation with no significant biological changes, although the physical uptake of carbon within the ocean is unable to keep pace with the rise in [CO₂].

The ocean carbon–climate feedback parameter, γ_O, is diagnosed from the difference between the COU model configuration and the BGC configuration and so includes the effect of an increasing surface warming under rising [CO₂] (Eq. 14). There is a broadly consistent response across models, with a model mean decrease in carbon inventory of around 80 PgC due to the additional warming in the COU configuration relative to the BGC configuration (Fig. 13, turquoise bars). The effect of this additional warming and the associated climate change leads to a decrease in both the saturated carbon and disequilibrium carbon of typically −60 and −70 PgC (Fig. 13, orange and blue bars), representing the decrease in solubility and decreased ocean ventilation. There is an increase in the regenerated carbon of typically 50 PgC (Fig. 13, green bars), which is due to a weaker circulation leading to a longer residence time of thermocline and deep waters, so there is more time for the accumulation of regenerated carbon below the mixed layer. The resulting γ_O is negative, indicating that the ocean takes up less carbon in response to the combination of surface warming and a weakening in ocean ventilation. This response involves a combination of chemical, physical, and biological changes where the warming reduces the solubility of the carbon in the ocean and a weakening in the circulation decreases the disequilibrium pool but lengthens the residence time and so increases the regenerated pool.

Overall, the ocean carbon inventory increases in the BGC configuration by 666 ± 53 PgC (model mean ± standard deviation) and decreases in COU relative to BGC by −80 ± 15 PgC. The resulting β_O is very similar across all the models (0.78 ± 0.06 PgC ppm⁻¹), reflecting not only the strong control of carbonate chemistry by rising atmospheric CO₂ (Katavouta et al., 2018) but also the use of similar carbonate chemistry schemes and bulk parameterizations of air–sea CO₂ fluxes across marine biogeochemical models (Séférian et al., 2020). The dominant contributions are composed of a positive contribution from the saturated carbon (3.66 ± 0.16 PgC ppm⁻¹) and a negative contribution from the disequilibrium carbon (−2.98 ±  0.16 PgC ppm⁻¹; see Table A3 in the Appendix); these intermodel differences are relatively small with ratios of the standard deviation to model mean of only 0.05 and 0.06, respectively. The regenerated contribution is over 2 orders of magnitude smaller than the sum of the saturated and disequilibrium contributions and so may be neglected for evaluating β_O.

Figure 13 Change in saturated, disequilibrium, and regenerated carbon pools in the fully coupled minus the biogeochemical simulation using Eq. (14), which contributes to the calculation of the carbon–climate feedback over the ocean. The sum of the partitions does not exactly match the total ocean uptake diagnosed from the air–sea fluxes due to land–ocean interactions involving storage in sediments and river inputs.

The values of γ_O differ more strongly across the models (−16.95 ± 5.62 PgC ^∘C⁻¹) and arise from differences in the extent of the surface warming and the dynamical changes in the ocean circulation and resulting changes in ventilation, residence time, and biological regeneration (Table A3). The contributions to γ_O include negative contributions from the saturated (−12.78 ± 2.50 PgC ^∘C⁻¹) and disequilibrium (−16.36 ± 5.31 PgC ^∘C⁻¹) components, which are partly opposed by a positive contribution from the regenerated component (12.25 ± 8.53 PgC ^∘C⁻¹). The largest intermodel differences are in the regenerated and disequilibrium responses and a relatively small spread in the saturated response, with the ratios of the standard deviation to the model mean being 0.70, 0.33, and 0.20, respectively (Table A3).

4.5 Transient climate response (TCR) and transient climate response to cumulative emissions (TCRE)

The idealized 1pctCO2 simulation is also routinely used for calculating two other climate metrics. The first is the transient climate response (TCR), which is defined as the temperature change relative to the preindustrial state at the time of CO₂ doubling ($\mathrm{\Delta }{T}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$), which occurs at 70 years after the start of the simulation. The second is the transient climate response to cumulative carbon emissions (TCRE), which is defined as the ratio of the TCR to cumulative fossil fuel emissions also at the time of CO₂ doubling (${\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$; Matthews et al., 2009), typically expressed in units of degrees Celsius per exagram of carbon (EgC⁻¹; 1 EgC = 1000 PgC).

$$\begin{array}{}\text{(21)}& \text{TCRE}={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}{{\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}}\end{array}$$

It has been shown that the TCRE is approximately constant over a wide range of cumulative emissions and emission pathways (MacDougall, 2016). Although non-CO₂ greenhouse gases and other climate forcings (e.g. aerosols and land use change) also affect the realized warming, the TCRE is considered to be a straightforward measure of peak warming caused by anthropogenic CO₂ emissions. The TCRE metric has gained significant policy relevance (Frame et al., 2014; IPCC, 2014; Millar et al., 2016), and it is a central component of frameworks used to calculate the remaining allowable carbon emissions to reach a specified temperature change target above the preindustrial level (Millar et al., 2017; Rogelj et al., 2018, 2019).

Table 4 Transient climate response (TCR; $\mathrm{\Delta }{T}_{\mathrm{2}\times C{O}_{\mathrm{2}}}$), diagnosed cumulative emissions at 2×CO₂ (${\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$), and transient climate response to cumulative emissions (TCRE) for the 11 CMIP6 models considered in this study.

Table 4 lists the TCR, ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$, and TCRE from the 11 CMIP6 models considered in this study. We calculate the TCR, following the standard approach, as the average temperature of 20 years (years 60–79) centered on the year when CO₂ doubles (year 70). The mean ± standard deviation for the TCR, ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$, and the TCRE from the 11 CMIP6 models considered here is 1.97 ± 0.39 ^∘C, 1121 ± 73 PgC, and 1.77 ± 0.37 ^∘C EgC⁻¹, respectively. For 15 CMIP5 models, Gillett et al. (2013) calculated the TCRE to be 1.63 ± 0.48 ^∘C EgC⁻¹ and a 5 %–95 % range for its observationally constrained value as 0.7–2.0 ^∘C EgC⁻¹. The CMIP5 and CMIP6 mean values quoted are not statistically different given the small sample size of available models, some of which duplicate processes or components.

The uncertainties in the TCRE stem from uncertainties in both the TCR and ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ which is directly affected by land and ocean carbon uptake. A large fraction of uncertainty in ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ comes from the diverse response of land carbon cycle models. Inclusion of the N cycle helps to reduce this uncertainty in terms of the spread across the land models (both the strength of their feedback parameters and diagnosed cumulative emissions). Cumulative diagnosed emissions for models with the land N cycle total 1088 ± 34 PgC compared to 1160 ± 91 PgC for models without the N cycle (from Table 4). For the results reported here from 11 CMIP6 models, however, the uncertainty in the TCR (mean ± standard deviation = 1.97 ± 0.39 ^∘C), as indicated by standard deviation normalized by the mean, is 3 times as big as the uncertainty in ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ (1121 ± 73 PgC). As a result, the TCR contributes about 90 % of the total variance in the calculated TCRE value (1.77 ± 0.37 ^∘C EgC⁻¹; see Sect. A6 in the Appendix).

The TCRE may also be expressed in terms of a product of a thermal contribution from the dependence of surface warming on radiative forcing and a carbon contribution from the dependence of radiative forcing on cumulative carbon emissions (Williams et al., 2016; Katavouta et al., 2018), as

$$\begin{array}{}\text{(22)}& \text{TCRE}={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}{\mathrm{\Delta }{R}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}}{\displaystyle \frac{\mathrm{\Delta }{R}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}{{\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}},\end{array}$$

where $\mathrm{\Delta }{R}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ is the change in radiative forcing relative to the preindustrial period. For a suite of 10 CMIP5 models, Williams et al. (2017) show that the intermodel spread in the TCRE calculated from the 1pctCO2 experiment has a larger contribution from the intermodel differences in the thermal contribution, $\frac{\mathrm{\Delta }{T}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}{\mathrm{\Delta }{R}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}$, due to climate feedback and ocean heat uptake over the first few decades, but the intermodel differences in the carbon contribution, $\frac{\mathrm{\Delta }{R}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}{{\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}$, due to land and ocean carbon uptake become of comparable importance after 80 years.

Thus we conclude, as shown by Jones and Friedlingstein (2020), that the contributions to TCRE uncertainty have changed since CMIP5 from being of similar magnitudes due to carbon feedbacks and climate feedbacks to now being dominated by climate feedbacks. The reduction in the spread of land carbon model feedbacks which we are beginning to see in CMIP6 has led to a reduction in the spread of the TCRE implying that a large fraction of uncertainty in the TCRE is now contributed to by physical climate system processes that determine the TCR. More CMIP6 models include a complete treatment of important processes – notably the terrestrial nitrogen cycle – that determine the airborne fraction and hence ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$. Reducing the uncertainty in land and ocean carbon uptake across models remains a priority and will contribute to further reducing the uncertainty in the estimates of the TCRE on centennial timescales.

5 Summary and conclusions

Model intercomparison projects offer several benefits including the calculation of the model mean response, the quantification of the uncertainty based on the spread across models, and showing how this uncertainty changes over time. The carbon feedback analysis presented here based on the C⁴MIP protocol of experiments (Jones et al., 2016) allows us to investigate how feedback strengths have evolved since CMIP5 and also to attempt to understand the reasons behind the spread in models.

The carbon uptake over land and ocean, in response to increasing atmospheric CO₂ concentration, is well known to be dominated by the positive contribution from the carbon–concentration feedback (Gregory et al., 2009; Arora et al., 2013). The strength of this feedback is of comparable magnitudes over land (mean ± standard deviation = 0.97 ± 0.40 PgC ppm⁻¹) and ocean (0.79 ± 0.07 PgC ppm⁻¹) although the feedback is much more uncertain over land as indicated by the standard deviation across the 11 models considered here. This dominant positive contribution from the carbon–concentration feedback is, however, opposed by the weaker negative carbon–climate feedback that is associated with the climate change that results due to increasing atmospheric CO₂. The absolute magnitude of this weaker negative feedback is about 3 times larger and an order of magnitude more uncertain over land (−45.1 ± 50.6 PgC ^∘C⁻¹) than over ocean (−17.2 ± 5.0 PgC ^∘C⁻¹). Model estimates of the ocean carbon–concentration feedback are consistent with each other, reflecting the strong control of how carbonate chemistry alters with rising atmospheric CO₂. There is a relatively wider range in the model estimates of the ocean carbon–climate feedback, particularly in terms of how changes in ocean circulation alter the disequilibrium and regeneration terms. Over land, however, since the carbon–concentration and carbon–climate feedbacks are determined entirely by biological process, which are much less understood, the resulting uncertainty is much higher across the land models than across the ocean models. This uncertainty in the strength of carbon–concentration and carbon–climate feedbacks over land is well known (Friedlingstein et al., 2006; Arora et al., 2013). The inclusion of the N cycle results in lower absolute strength of the feedback parameters over land. In addition, the land models that include a representation of the N cycle exhibit a reduced spread in their feedback parameters, despite the additional complexity, compared to when all models are considered. This suggests that if all models were to include the N limitation of photosynthesis, the spread across them will potentially reduce.

The additional analyses that we have performed provide insight into the reasons for the diverse responses among models, especially for land models. Over land, the diverse response of models is found to be primarily due to the wide range of the strength of the CO₂ fertilization effect, the fraction of GPP that is converted to NPP, and the residence times of carbon in the live (vegetation) and dead (litter plus soil) carbon pools across models. There is more consistency in the response of the ocean models, although intermodel differences arise from differences in the ventilation and residence time, altering the ocean disequilibrium and regenerated carbon.

In regard to the TCRE, while its uncertainty is dominated by physical processes affecting the thermal response involving climate feedbacks and heat uptake on decadal timescales, a reduction in the uncertainty in land and ocean carbon uptake across models will reduce the uncertainty in the TCRE on centennial timescales.

Finally, the decision to use fully and biogeochemically coupled configurations of the 1pctCO2 experiment as the standard simulations from which to diagnose carbon cycle and climate system feedbacks should provide consistency and continuity for future versions of Earth system models that can be compared against their predecessors.

Appendix A: The climate carbon cycle feedbacks framework

The rate of change in carbon in the combined atmosphere–land–ocean system is written as

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{G}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{A}}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{L}}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{O}}}{\mathrm{d}t}}=E,\end{array}$$

where the global carbon pool C_G = C_A + C_L + C_O is the sum of carbon in the atmosphere, land, and ocean components (PgC) and E is the rate of anthropogenic CO₂ emissions (PgC yr⁻¹) into the atmosphere. The equations for the atmosphere, land, and ocean are

$$\begin{array}{}\text{(A2)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{A}}}{\mathrm{d}t}}=F_{\mathrm{A}}\left(T,c\right)+E,\\ & {\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{L}}}{\mathrm{d}t}}=F_{\mathrm{L}}\left(T,c\right),\\ & {\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{O}}}{\mathrm{d}t}}=F_{\mathrm{O}}\left(T,c\right),\end{array}\end{array}$$

where $(F_{\mathrm{L}}+F_{\mathrm{O}})=-F_{\mathrm{A}}$ is the fluxes between the atmosphere and the underlying land and ocean, taken to be positive into the components. The fluxes F are expressed as functions of surface temperature T and the surface atmospheric CO₂ concentration c. Here and subsequently, uppercase C denotes carbon pools and lowercase c denotes atmospheric CO₂ concentration.

Figure A1 Individual model values from CMIP6 models for globally averaged surface temperature change (a–c), cumulative atmosphere–land CO₂ flux (d–f), and cumulative atmosphere–ocean CO₂ flux (g–h) from the fully, biogeochemically, and radiatively coupled versions of the 1pctCO2 experiment.

Figure A2 Components of the carbon budget terms in cumulative emissions from the 11 participating CMIP6 models based on Eq. (A6) in (a) and Eq. (A7) in (b), using results from the fully coupled 1 % yr⁻¹ increasing CO₂ simulation at 2×CO₂ (year 70) in contrast to Fig. 4 which showed these results at 4×CO₂. The models are arranged in ascending order based on their cumulative emissions values. Results from participating CMIP5 models in the A13 study are shown in (c, d). In addition, ESMs whose land component includes a representation of the N cycle are identified by a red font colour for cumulative land carbon uptake (a, c) and fractional emissions taken up by land (b, d). The model mean is shown not only for all models but also separately for models whose land components include or do not include a representation of the N cycle.

In the fully, biogeochemically, and radiatively coupled versions of the 1pctCO2 experiments analyzed here, the rate of change in atmospheric carbon dC_A / dt is specified in Eqs. (A1) and (A2). The uptake or release of CO₂ by the underlying land and ocean yields an effective emission E which serves to maintain the budget.

Figure A3 Absolute amounts and the change from the beginning of the BGC simulation for carbon in soil + litter (a, b) and vegetation (c, d) pools.

The changes in atmosphere carbon budgets, from the preindustrial control simulation, in the differently coupled simulations are represented as

$$\begin{array}{}\text{(A3a)}& {\displaystyle }\begin{array}{rl}& \text{radiatively coupled}\\ & {\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{A}}^{\prime }}{\mathrm{d}t}}-E^{+}=F_{\mathrm{A}}^{+}=-F_{\mathrm{L}}^{+}-F_{\mathrm{O}}^{+}={\mathrm{\Gamma }}_{\mathrm{A}}{T}^{+},\end{array}\end{array}$$

$$\begin{array}{}\text{(A3b)}& {\displaystyle }\begin{array}{rl}& \text{biogeochemically coupled}\\ & {\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{A}}^{\prime }}{\mathrm{d}t}}-E^{\ast }=F_{\mathrm{A}}^{\ast }=-F_{\mathrm{L}}^{\ast }-F_{\mathrm{O}}^{\ast }={\mathrm{\Gamma }}_{\mathrm{A}}{T}^{\ast }+B_{\mathrm{A}}{c}^{\prime },\end{array}\end{array}$$

$$\begin{array}{}\text{(A3c)}& {\displaystyle }\begin{array}{rl}& \text{fully coupled}\\ & {\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{A}}^{\prime }}{\mathrm{d}t}}-E=F_{\mathrm{A}}^{\prime }=-F_{\mathrm{L}}^{\prime }-F_{\mathrm{O}}^{\prime }={\mathrm{\Gamma }}_{\mathrm{A}}{T}^{\prime }+B_{\mathrm{A}}{c}^{\prime },\end{array}\end{array}$$

which serve to define the instantaneous carbon–concentration (B_A) and carbon–climate (Γ_A) feedback parameters and assume linearization of the globally integrated surface–atmosphere CO₂ flux in terms of global mean temperature and concentration change. In Eq. (A3), F⁺, F^∗, and F^′ are the flux changes and T⁺, T^∗, and T^′ the temperature changes in the radiatively, biogeochemically, and fully coupled simulations, and E⁺, E^∗, and E are the resulting implicit emissions. c^′ is the specified CO₂ concentration change above its preindustrial level in the 1pctCO2 simulations. In the biogeochemically coupled simulation there is no radiative forcing due to increasing CO₂, so T^∗ is small, although not zero, and exhibits a distinct spatial pattern. The assumption made in Eq. (A3) is that the feedback parameters are the same in the three cases.

Carbon budget changes for the land component parallel to Eq. (A3) but without the emissions terms as

$$\begin{array}{}\text{(A4a)}& {\displaystyle }\begin{array}{rl}& \text{radiatively coupled}\\ & {\displaystyle \frac{{{\mathrm{dC}}^{\prime }}_{\mathrm{L}}}{\mathrm{d}t}}=F_{\mathrm{L}}^{+}={\mathrm{\Gamma }}_{\mathrm{L}}{T}^{+},\end{array}\end{array}$$

$$\begin{array}{}\text{(A4b)}& {\displaystyle }\begin{array}{rl}& \text{biogeochemically coupled}\\ & {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{L}}^{\ast }}{\mathrm{d}t}}=F_{\mathrm{L}}^{\ast }={\mathrm{\Gamma }}_{\mathrm{L}}{T}^{\ast }+B_{\mathrm{L}}{c}^{\prime },\end{array}\end{array}$$

$$\begin{array}{}\text{(A4c)}& {\displaystyle }\begin{array}{rl}& \text{fully coupled}\\ & {\displaystyle \frac{{\mathrm{dC}}_{\mathrm{L}}^{\ast }}{\mathrm{d}t}}=F_{\mathrm{L}}^{\prime }={\mathrm{\Gamma }}_{\mathrm{L}}{T}^{\prime }+B_{\mathrm{L}}{c}^{\prime },\end{array}\end{array}$$

and similarly for the ocean component. Since $F_{\mathrm{A}}=-\left(F_{\mathrm{L}}+F_{\mathrm{O}},\right)$ it follows that ${\mathrm{\Gamma }}_{\mathrm{A}}=-\left({\mathrm{\Gamma }}_{\mathrm{L}}+{\mathrm{\Gamma }}_{\mathrm{O}}\right)$ and $B_{\mathrm{A}}=-\left(B_{\mathrm{L}}+B_{\mathrm{O}}\right)$. There are no terms involving c^′ in the radiatively coupled simulation (Eqs. A3a and A4a) since the preindustrial value of atmospheric CO₂ concentration is prescribed for the biogeochemistry components, so $c^{\prime }=\mathrm{0}$ and does not affect the flux.

The instantaneous feedback parameters (B_L and Γ_L) differ from those in the integrated flux approach of Friedlingstein et al. (2006), who express time-integrated flux changes (i.e. change in pool or reservoir sizes) as functions of temperature and CO₂ concentration changes with

$$\begin{array}{}\text{(A5a)}& {\displaystyle }\begin{array}{rl}& \text{radiatively coupled}\\ & \int F_{\mathrm{L}}^{+}={\mathrm{\Delta }\mathrm{C}}_{\mathrm{L}}^{+}={\mathit{\gamma }}_{\mathrm{L}}{T}^{+},\end{array}\end{array}$$

$$\begin{array}{}\text{(A5b)}& {\displaystyle }\begin{array}{rl}& \text{biogeochemically coupled}\\ & \int F_{\mathrm{L}}^{\ast }={\mathrm{\Delta }\mathrm{C}}_{\mathrm{L}}^{\ast }={\mathit{\gamma }}_{\mathrm{L}}{T}^{\ast }+{\mathit{\beta }}_{\mathrm{L}}{c}^{\prime },\end{array}\end{array}$$

$$\begin{array}{}\text{(A5c)}& {\displaystyle }\begin{array}{rl}& \text{fully coupled}\\ & \int {F^{\prime }}_{\mathrm{L}}={\mathrm{\Delta }{\mathrm{C}}^{\prime }}_{\mathrm{L}}={\mathit{\gamma }}_{\mathrm{L}}{T}^{\prime }+{\mathit{\beta }}_{\mathrm{L}}{c}^{\prime },\end{array}\end{array}$$

and similarly for the ocean component, with the assumption that the ${\mathrm{\Delta }{\mathrm{C}}^{\prime }}_{\mathrm{O}}$ term includes changes in the carbon amount of ocean sediment as well.

The units of instantaneous and integrated flux-based parameters are different (Γ–PgC yr^−1 ∘C⁻¹, B– PgC yr⁻¹ ppm⁻¹ and γ–PgC ^∘C⁻¹, β–PgC ppm⁻¹). Arora et al. (2013) show how the instantaneous and integrated flux-based feedback parameters are related to each other

Integrating Eqs. (A1) and (A2) from initial time to t gives

$$\begin{array}{}\text{(A6)}& {\mathrm{\Delta }{\mathrm{C}}^{\prime }}_{\mathrm{A}}+{\mathrm{\Delta }{\mathrm{C}}^{\prime }}_{\mathrm{L}}+{\mathrm{\Delta }{\mathrm{C}}^{\prime }}_{\mathrm{O}}={\int }_{\mathrm{0}}^{t}E\mathrm{d}t=\stackrel{\mathrm{̃}}{E}.\end{array}$$

Here ${\mathrm{\Delta }{\mathrm{C}}^{\prime }}_{\mathrm{A}}=\mathrm{2.12}\left(c\left(t\right)-c\left(\mathrm{0}\right)\right)$ is the change in atmospheric carbon burden (the factor 2.12 converts atmospheric CO₂ concentration from parts per million to atmospheric burden in petagrams of carbon) and ${\mathrm{\Delta }{\mathrm{C}}^{\prime }}_X={\int }_{\mathrm{0}}^t{F^{\prime }}_X\mathrm{d}t,X=\text{L,\hspace{0.17em}O}$ is the cumulative flux equal to the change in the land or ocean carbon pool for the fully coupled simulation. The terms in Eq. (A6) indicate the contribution of changes in atmosphere, land, and ocean carbon pools to cumulative emissions $\stackrel{\mathrm{̃}}{E}$. Finally, division by the cumulative emissions term in Eq. (A6) gives all the terms in a fractional form as

$$\begin{array}{}\text{(A7)}& f_{\mathrm{A}}+f_{\mathrm{L}}+f_{\mathrm{O}}=\mathrm{1},\end{array}$$

where f_A is the airborne fraction of cumulative emissions and f_L and f_O are fractional emissions taken up by the land and ocean. These components are evaluated at the time of CO₂ quadrupling.

A1 Reasons for nonlinearity in the ocean C cycle response in the CNRM-ESM2-1 model

In Fig. 6 the value of γ_O changes sign for the CNRM-ESM2-1 model from positive when calculated using the RAD–BGC or RAD–COU approaches to negative when calculated using the BGC–COU approach. This nonlinear behaviour for a previous version of the CNRM model was documented in Schwinger et al. (2014) and is caused by the large increase in regenerated DIC in the RAD simulation, similar to the increase in the COU relative to the BGC simulation, as shown in Fig. 13 for the CNRM-ESM2-1 model. This nonlinear behaviour is stronger in CNRM-ESM2-1 compared to CNRM-ESM1, its previous version (Séférian et al., 2016), most likely due to a new parameterization for N fixation which increases ocean NPP and a revised parameterization for organic matter remineralization (in PISCESv2-gas). A contribution to a positive γ_O is also made by declining sea ice in the RAD simulation which leads to changes in the sign of the air–sea carbon exchange in the Southern Ocean. The vertical profile of dissolved inorganic carbon in the Southern Ocean in BGC and COU simulations (with rising [CO₂]) is different from that in the RAD simulation (for the preindustrial [CO₂]), and this leads to additional nonlinearities.

A2 Additional figures and discussion

Figure A1 shows results from individual CMIP6 models for which model means and ranges were shown in Figs. 1, 2, and 3 and allows for the identification of models which behave differently compared to the majority of models. In Fig. A1a and c, CanESM5 shows the largest temperature increase and NorESM2-LM and MIROC-ES2L the smallest in response to increase in [CO₂] for the COU and RAD simulations, respectively. For cumulative atmosphere–land CO₂ flux in the COU simulation (panel d), CanESM5 simulates the largest land carbon uptake and ACCESS-ESM1.5 the smallest. This is not the case for the BGC simulation (panel e) where land carbon uptakes from BCC-CSM2-MR and CNRM-ESM2.1 are the largest among all models, while land carbon uptake from ACCESS-ESM1.5 is the lowest. Finally, in the RAD simulation (panel f) the loss of carbon from land in response to increasing temperatures is lowest in MPI-ESM1.2-LR and largest in BCC-CSM2-MR. Over the ocean, while most models behave very similarly, the carbon uptakes in BCC-CSM2-MR, ACCESS-ESM1.5, and NOAA-GFDL-ESM4 are larger than in most models in the COU and BGC simulations. In the RAD simulation, almost all models simulate a loss of carbon from the ocean, but CNRM-ESM2.1 shows a small uptake.

A3 Additional tables

Table A1 Values of carbon–concentration and carbon–climate feedback parameters for land and ocean calculated using the BGC–COU approach (using results from the BGC and COU simulations) and the linear transient climate sensitivity to CO₂, from CMIP6 and CMIP5 models at 4×CO₂ (i.e. at the end of the 1pctCO2 simulation) and 2×CO₂.

A4 Model descriptions

A4.1 Commonwealth Scientific and Industrial Research Organisation (CSIRO) ACCESS-ESM1.5

The Australian Community Climate and Earth System Simulator ACCESS-ESM1.5 (Ziehn et al., 2020) is comprised of a number of component models. The atmospheric model is the UK Met Office Unified Model version 7.3 (Martin et al., 2010, 2011) with their land surface model replaced with the Community Atmosphere Biosphere Land Exchange (CABLE) model (Kowalczyk et al., 2013). The ocean component is the NOAA GFDL Modular Ocean Model (MOM) version 5 (Griffies, 2014) with the same configuration as the ocean model component of ACCESS1.0 and ACCESS1.3 (Bi et al., 2013). Sea ice is simulated using the LANL CICE4.1 model (Hunke and Lipscomb, 2010). Coupling of the ocean and sea ice to the atmosphere is through the OASIS-MCT coupler (Valcke, 2013). The physical climate model configuration used here is very similar to the version (ACCESS1.3) that contributed to the Coupled Model Intercomparison Project Phase 5 (CMIP5; Bi et al., 2013). The carbon cycle is included in ACCESS through the CABLE land surface model and its biogeochemistry module, CASA-CNP (Wang et al., 2010), and through the World Ocean Model of Biogeochemistry and Trophic-dynamics (WOMBAT; Oke et al., 2013).

The WOMBAT model is based on an NPZD (nutrient, phytoplankton, zooplankton, and detritus) model with the additions of bioavailable iron limitation, dissolved inorganic carbon, calcium carbonate, alkalinity, and oxygen. Productivity drives uptake and formation of carbon and oxygen that are exchanged with the atmosphere. The sinking and remineralization of detritus carries biogeochemical tracers to the deep ocean. Iron is supplied by dust deposition, continental shelves, and background ocean values.

Table A2 Estimate of the change in the ocean carbon inventory (PgC) expected from a time integral of the global air–sea carbon flux into the ocean versus the volume integral of the change in the dissolved inorganic carbon, together with the small residual. The time integral of the air–sea carbon flux provides the dominant contribution to the change in the ocean carbon inventory, although there is a small mismatch due to the land-to-ocean carbon flux from river runoff and the ocean-to-land carbon flux from carbon burial in ocean sediments.

The Australian community model CABLE simulates the fluxes of momentum, heat, water, and carbon at the surface. The biogeochemistry module CASA-CNP simulates the flow of carbon and nutrients such as nitrogen and phosphorus between three plant biomass pools (leaf, wood, root), three litter pools (metabolic, structural, coarse woody debris), and three organic soil pools (microbial, slow, passive), plus one inorganic soil mineral nitrogen pool and three phosphorus soil pools.

In the CABLE configuration applied here the land surface is represented by 10 vegetation and 3 nonvegetation land cover types. CABLE calculates gross primary production (GPP) and leaf respiration at every time step using a two-leaf canopy scheme (Wang and Leuning, 1998) as a function of the leaf area index (LAI). This set-up uses a simulated (prognostic) LAI based on the size of the leaf carbon pool and the specific leaf area. Daily mean GPP and leaf respiration values are then passed onto CASA-CNP to calculate daily respiration fluxes and the flow of carbon and nutrients between the pools. Similar to the previous version, ACCESS-ESM1 (Law et al., 2017; Ziehn et al., 2017), the model is run with nitrogen and phosphorus limitation enabled.

A4.2 Beijing Climate Center (BCC) Climate System Model version 2 with medium resolution (BCC-CSM2-MR)

BCC-CSM2-MR (Wu et al., 2019) is the second generation of the BCC model with medium resolution that was released to run CMIP6 simulations. It is a fully coupled global climate model and updated from its previous version of BCC-CSM1.1 used for CMIP5 (Wu et al., 2013). The atmospheric component of BCC-CSM2-MR is the BCC Atmospheric General Circulation Model version 3 (BCC-AGCM3-MR; Wu et al., 2019). The land component is the BCC Atmosphere and Vegetation Interaction Model version 2.0 (BCC-AVIM2; Li et al., 2019) with the terrestrial carbon cycle. The oceanic component is the Modular Ocean Model version 4 with 40 levels (hereafter MOM4-L40). The sea ice component is the NOAA GFDL Sea Ice Simulator (SIS). These components are physically coupled through fluxes of momentum, energy, water, and carbon at their interfaces. The coupling was realized with the flux coupler version 5 developed by the National Center for Atmosphere Research (NCAR).

Table A3 Carbon cycle feedback parameters for the ocean, β_O and γ_O, diagnosed from the air–sea carbonfluxes and separately diagnosed for the ocean carbon inventory and its separate ocean saturated, disequilibrium, and regenerated DIC pools for the subset of eight CMIP6 models for which 3D ocean data were available; their sum does not exactly match the diagnostics from the air–sea fluxes due to land–ocean interactions involving storage in sediments and river inputs.

The atmospheric component of BCC-CSM2-MR has a T106 horizontal resolution of approximately 1.125^∘ and 46 vertical levels in a hybrid sigma–pressure vertical coordinate system with the top level at 1.459 hPa. The ocean component resolution of BCC-CSM2-MR is 1^∘ longitude by $\mathrm{1}/\mathrm{3}{}^{\circ }$ latitude between 30^∘ S and 30^∘ N which increases to 1^∘ latitude at 60^∘ S and 60^∘ N and nominally 1^∘ polarward with tripolar coordinates, and there are 40 z levels in the vertical plane.

The atmospheric component model BCC-AGCM3-MR in BCC-CSM2-MR is developed from its previous CMIP5 version (Wu et al., 2008). The main updates include a modification of deep convection parameterization, a new scheme for cloud fraction, indirect effects of aerosols through clouds and precipitation, and the gravity wave drag generated by deep convection (Wu et al., 2019). Atmospheric CO₂ concentration in BCC-AGCM3-MR for this work is a prognostic variable and calculated through a budget equation which considers advective transport in the atmosphere, anthropogenic CO₂ emissions, and interactive CO₂ fluxes at the interfaces with land and ocean. But chemical processes are not taken into account. The terrestrial carbon cycle in BCC-AVIM2 (Li et al., 2019) operates through a series of biochemical and physiological processes acting on the photosynthesis and respiration of vegetation and takes into account carbon loss due to turnover and mortality of vegetation and CO₂ release into the atmosphere through soil respiration. The vegetation litter on the ground surface and in the soil is divided into eight terrestrial carbon pools (surface structural, surface metabolic, surface microbial, soil structural, soil metabolic, soil microbial, slow, and passive carbon pools) according to the timescale of the decomposition of carbon in each pool and the transfers between different pools. Allocation to and from the three vegetation biomass pools (leaf, stem, root) leads to dynamic vegetation that in turn produces litterfall and which ultimately transfers carbon to soil organic matter. The allocation of carbon to the three vegetation biomass pools is dependent on light availability, water stress, and the phenology stages of the canopy and follows the formulations of Arora and Boer (2005).

The biogeochemistry module to simulate the ocean carbon cycle in MOM4_L40 is based on the protocols from the Ocean Carbon-Cycle Model Intercomparison Project Phase 2 (OCMIP2; http://ocmip5.ipsl.jussieu.fr/OCMIP/phase2/, last access: 15 December 2019). The OCMIP biogeochemistry module parameterizes the process of marine biology in terms of geochemical fluxes, without explicit representation of the marine ecosystem and food web processes, and includes five prognostic variables: phosphate, dissolved organic phosphorus, dissolved oxygen, dissolved inorganic carbon, and alkalinity. Ocean carbon cycle processes in BCC-CSM2-MR follow OCMIP, except for in parameterizing the export of organic matter from surface waters to deep oceans (Wu et al., 2013).

A4.3 Canadian Centre for Climate Modelling and Analysis (CCCma) fifth-generation Earth System Model, CanESM5

CanESM5 has evolved from its predecessor CanESM2 (Arora et al., 2011) that was used in the Coupled Model Intercomparison Project phase 5 (CMIP5). CanESM5 represents a major update to CanESM2 and is described in detail in Swart et al. (2019). The major changes relative to CanESM2 are the implementation of completely new models for the ocean, sea ice, and marine ecosystems and a new coupler. The resolution of CanESM5 (T63 or $\sim \mathrm{2.8}{}^{\circ }$ in the atmosphere and $\sim \mathrm{1}{}^{\circ }$ in the ocean) remains similar to CanESM2 and is at the lower end of the spectrum of CMIP6 models. The atmospheric component of CanESM5 is represented by version 5 of the Canadian Atmospheric Model (CanAM5) and has several improvements relative to its predecessor, CanAM4 (von Salzen et al., 2013), including changes to aerosol, clouds, radiation, land surface, and lake processes. CanAM5 uses a triangular spectral truncation in the model dynamical core, with an approximate horizontal resolution of 2.8^∘ in latitude and longitude. It uses a hybrid vertical coordinate system with 49 levels between the surface and 1 hPa, with a vertical resolution of about 100 m near the surface. Relative to the 35 levels used in CanESM2 most of the additional 14 levels were added in the upper troposphere and stratosphere.

The land surface in CanESM5 is modelled using the Canadian Land Surface Scheme (CLASS; Verseghy, 2000) and the Canadian Terrestrial Ecosystem Model (CTEM; Arora and Boer, 2005, 2010), which together form the land component of CanESM5. CLASS and CTEM simulate the physical and biogeochemical land surface processes, respectively, and together they calculate fluxes of energy, water, CO₂, and wetland CH₄ emissions at the land–atmosphere boundary. Over land, three permeable soil layers are used with default thicknesses of 0.1, 0.25, and 3.75 m for which liquid and frozen soil moisture and temperature are prognostically calculated. The depth to bedrock is specified on the basis of the global dataset which reduces thicknesses of the permeable soil layers where soil depth is less than 4.1 m. Snow is represented using one layer, whose snow water equivalent and temperature are modelled prognostically. The introduction of dynamic wetlands and their methane emissions is a new biogeochemical process added since the CanESM2 (Arora et al., 2018). The nitrogen cycle over land is not represented but the parameterization of photosynthesis downregulation as CO₂ increases is included. The magnitude of the parameter representing this downregulation is increased in CanESM5 compared to CanESM2, following Arora and Scinocca (2016), who found the best value of this parameter that reproduced various aspects of the historical carbon budget for CanESM4.2 (a model version more similar to CanESM2 than CanESM5). Other than wetlands and the changes to the strength of the CO₂ fertilization effect, the remaining terrestrial ecosystem processes are represented in the same way as in CanESM2.

The physical ocean component of CanESM5 is based on Nucleus for European Modelling of the Ocean (NEMO) version 3.4.1. It is configured on the tripolar ORCA1 C grid with 45 z-coordinate vertical levels, varying in thickness from ∼6 m near the surface to ∼250 m in the abyssal ocean. The horizontal resolution is based on a 1^∘ Mercator grid, varying with the cosine of latitude, with a refinement of the meridional grid spacing to $\mathrm{1}/\mathrm{3}{}^{\circ }$ near the Equator. Two modifications have been introduced to the NEMO mesoscale and small-scale mixing physics in CanESM5, and these are detailed in Swart et al. (2019). Sea ice is represented using the LIM2 sea ice model (Fichefet and Morales Maqueda, 1997; Bouillon et al., 2009), which is run within the NEMO framework.

The ocean carbon cycle is represented using the Canadian Model of Ocean Carbon (CMOC), which was developed for earlier versions of CanESM (Christian et al., 2010; Arora et al., 2011), and includes carbon chemistry and biology. The biological component is a simple nutrient–phytoplankton–zooplankton–detritus (NPZD) model, with fixed Redfield stoichiometry, and simple parameterizations of iron limitation, nitrogen fixation, and export flux of calcium carbonate.

A4.4 Community Earth System Model version 2 (CESM2)

CESM2 (Danabasoglu et al., 2020) contains substantial improvements on CESM1. The resolution remains the same as in CESM1 (0.9^∘ latitude × 1.25^∘ longitude for the atmosphere and land with 32 vertical atmospheric levels and 25 ground levels and $\sim \mathrm{1}{}^{\circ }$ for the ocean). The Community Atmosphere Model version 6 includes many changes to the representation of physical processes with the primary change being the inclusion of the Cloud Layers Unified By Binormals (CLUBB) unified turbulence scheme.

The CESM2 ocean component (Parallel Ocean Program version 2, POP2) is largely the same as that used in CESM1 except with a new parameterization for mixing effects in estuaries along with several other numerical and physics improvements. The sea ice model is CICE version 5.1.2 (CICE5; Hunke et al., 2015) . Ocean biogeochemistry is represented by the Marine Biogeochemistry Library (MARBL). MARBL represents multiple nutrient colimitations (N, P, Si, and Fe). It includes three explicit phytoplankton functional groups (diatoms, diazotrophs, and picophytoplankton and nanophytoplankton), one implicit phytoplankton group (calcifiers), and one zooplankton group. MARBL includes prognostic carbonate chemistry and simulates sinking particulate organic matter. Major updates relative to CESM1 include a representation of subgrid-scale variations in light and variable C : P stoichiometry. The atmospheric deposition of iron is computed prognostically in CESM2 as a function of dust and black carbon deposition simulated by CAM6. Riverine nutrient, carbon, and alkalinity fluxes are supplied to the ocean from a dataset.

The land component is the Community Land Model version 5 (CLM5; Lawrence et al., 2019) which simulates land water, energy, momentum, and carbon and nitrogen cycling. CLM5 includes an extensive suite of new and updated processes and parameterizations that collectively improve the model's hydrological, biogeochemical, and ecological realism and enhance the representation of anthropogenic land use activities on climate and the carbon cycle. The primary updates are as follows with details, references, and additional updates described and listed in Lawrence et al. (2019): (1) updated parameterizations and structure for hydrology and snow (spatially explicit soil depth, dry surface layer, revised groundwater scheme, revised canopy interception and canopy snow processes, updated fresh snow density, and inclusion of the Model for Scale Adaptive River Transport); (2) a plant hydraulics scheme to more mechanistically represent plant water use and limitation; (3) vertically resolved soil biogeochemistry with base organic matter decomposition rates varying with depth and modified by soil temperature, water, and oxygen limitation and nitrification and denitrification updated as in the CENTURY model; (4) a methane production, oxidation, and emissions model; (5) improved representation of plant N dynamics to address deficiencies in CLM4 through the introduction of flexible plant carbon : nitrogen (C:N) stoichiometry which avoids the problematic CLM4 separation of potential and actual plant productivity, explicitly simulating the photosynthetic capacity response to environmental conditions through the Leaf Utilization of Nitrogen for Assimilation (LUNA) module and accounting for how N availability affects plant productivity through the Fixation and Uptake of Nitrogen (FUN) module which determines the C costs of N acquisition; methane emissions and oxidation from natural land processes; (6) a global active crop model with six crop types and time-evolving irrigated areas and industrial fertilization rates; (7) updated canopy processes including a revised canopy radiation scheme and canopy scaling of leaf processes, colimitations on photosynthesis, and updated stomatal conductance; (8) a new fire model that includes representation of natural and anthropogenic ignition sources and suppression along with agricultural, deforestation, and peat fires; and (9) inclusion of carbon isotopes.

A4.5 Centre National de Recherches Météorologiques (CNRM) CNRM-ESM2-1

CNRM-ESM2-1 is the second-generation Earth System model developed by CNRM and Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique (CERFACS) for CMIP6 (Séférian et al., 2019).

The atmosphere component of CNRM-ESM2-1 is based on version 6.3 of the global spectral model ARPEGE-Climat (ARPEGE-Climat_v6.3). ARPEGE-Climat resolves atmospheric dynamics and thermodynamics on a T127 triangular grid truncation that offers a spatial resolution of about 150 km in both longitude and latitude. CNRM-ESM2-1 employs a “high-top” configuration with 91 vertical levels that extend from the surface to 0.01 hPa in the mesosphere; 15 hybrid σ–pressure levels are available below 1500 m.

The surface state variables and fluxes at the surface–atmosphere interface are simulated by the SURFEX modelling platform version 8.0 over the same grid and with the same time step as the atmosphere model. SURFEX v8.0 encompasses several submodules for modelling the interactions between the atmosphere, the ocean, the lakes, and the land surface. Over the land surface, CNRM-ESM2-1 uses the ISBA-CTRIP land surface modelling system (http://www.umr-cnrm.fr/spip.php?article1092&lang=en, last access: 15 December 2019) to solve energy, carbon, and water budgets at the land surface (Decharme et al., 2019; Delire et al., 2019). Its physical core explicitly solves the one-dimensional Fourier's and Darcy's laws throughout the soil, accounting for the hydraulic and thermal properties of soil organic carbon. It uses a 12-layer snow model of intermediate complexity that allows for separate water and energy budgets for the soil and the snowpack. It accounts for a dynamic river flooding scheme in which floodplains interact with the soil and the atmosphere through free-water evaporation, infiltration, and precipitation interception and a two-dimensional diffusive groundwater scheme to represent unconfined aquifers and upward capillarity fluxes into the superficial soil. More details on these physical aspects can be found in Decharme et al. (2019).

To simulate the land carbon cycle and vegetation–climate interactions, ISBA-CTRIP simulates plant physiology, carbon allocation and turnover, and carbon cycling through vegetation, litter, and soil. It includes a module for wild fires, land use and land cover changes, and carbon leaching through the soil and transport of dissolved organic carbon to the ocean. Leaf photosynthesis is represented by the semiempirical model proposed by Goudriaan et al. (1985). Canopy-level assimilation is calculated using a 10-layer radiative transfer scheme including direct and diffuse radiation. Vegetation in ISBA is represented by four carbon pools for grasses and crops (leaves, stem, roots and a nonstructural carbohydrate storage pool) with two additional pools for trees (aboveground wood and coarse roots). Leaf phenology results directly from the carbon balance of the leaves. The model distinguishes 16 vegetation types (10 tree and shrub types, 3 grass types, and 3 crop types) alongside desert, rocks, and permanent snow. In the absence of nitrogen cycling within the vegetation, an implicit nitrogen limitation scheme that reduces specific leaf area with increasing CO₂ concentration was implemented in ISBA following the meta-analysis of Yin (2002). Additionally, there is an ad hoc representation of photosynthesis downregulation. The litter and soil organic matter module is based on the soil carbon part of the CENTURY model (Parton et al., 1988). The four litter and three soil carbon pools are defined based on their location aboveground or belowground and potential decomposition rates. The litter pools are supplied by the flux of dead biomass from each biomass reservoir (turnover). Decomposition of litter and soil carbon releases CO₂ (heterotrophic respiration). During the decomposition process, some carbon is dissolved by water slowly percolating through the soil column. This dissolved organic carbon is transported by the rivers to the ocean. A detailed description of the terrestrial carbon cycle can be found in Delire et al. (2019).

The ocean component of CNRM-ESM2-1 is the Nucleus for European Modelling of the Ocean (NEMO) version 3.6 (Madec et al., 2016) which is coupled to both the Global Experimental Leads and ice for ATmosphere and Ocean (GELATO) sea ice model (Salas Mélia, 2002) version 6 and also the marine biogeochemical model Pelagic Interactions Scheme for Carbon and Ecosystem Studies version 2 gas (PISCESv2-gas). NEMOv3.6 operates on the ORCA1L75 grid (Mathiot et al., 2017) which offers a nominal resolution of 1^∘ to which a latitudinal grid refinement of $\mathrm{1}/\mathrm{3}{}^{\circ }$ is added in the tropics; this grid describes 75 ocean vertical layers using a vertical z^∗ coordinate with partial step bathymetry formulation (Bernard et al., 2006).

The atmospheric chemistry scheme of CNRM-ESM2-1 is Reactive Processes Ruling the Ozone Budget in the Stratosphere version 2 (REPROBUS-C_v2). This scheme resolves the spatial distribution of 63 chemistry species but does not represent the low-troposphere ozone nonmethane hydrocarbon chemistry. CNRM-ESM2-1 also includes an interactive tropospheric aerosol scheme included in the atmospheric component ARPEGE-Climat. This aerosol scheme, named Tropospheric Aerosols for ClimaTe In CNRM (TACTIC_v2), represents the main anthropogenic and natural aerosol species of the troposphere.

The ocean biogeochemical component of CNRM-ESM2-1 uses the PISCESv2-gas model, which derives from PISCESv2 as described in Aumont et al. (2015). PISCESv2-gas simulates the distribution of five nutrients (from macronutrients nitrate, ammonium, phosphate, and silicate to the micronutrient iron) which regulate the growth of two explicit phytoplankton classes (nanophytoplankton and diatoms). Dissolved inorganic carbon (DIC) and alkalinity (Alk) are involved in the computation of the carbonate chemistry, which is resolved by routines to model the ocean carbonate system version 2 (mocsy 2.0; Orr and Epitalon, 2015) in PISCESv2-gas. The use of mocsy 2.0 enables a better and faster resolution of the ocean carbonate chemistry at thermodynamic equilibria. Oxygen is prognostically simulated using two different oxygen-to-carbon ratios, one when ammonium is converted to or mineralized from organic matter and the other when oxygen is consumed during nitrification. Their values have been set to 131∕122 and 32∕122, respectively.

At the ocean surface, PISCESv2-gas exchanges carbon, oxygen, dimethylsulfide (DMS), and nitrous oxide (N₂O) tracers with the atmosphere using the revised air–sea exchange bulk as published by Wanninkhof (2014). PISCESv2-gas uses several boundary conditions which represent the supply of nutrients from five different sources: atmospheric deposition, rivers, sediment mobilization, sea ice, and hydrothermal vents.

A4.6 Institut Pierre Simon Laplace (IPSL) IPSL-CM6A-LR

IPSL-CM6A-LR is the coupled climate model of the Institut Pierre Simon Laplace. It results from the integration of the following components: the LMDZ atmospheric general circulation model (version 6A-LR; Hourdin et al., 2019), the NEMO oceanic model (version 3.6; Vancoppenolle et al., 2009; Aumont et al., 2015; Rousset et al., 2015; Madec et al., 2016), and the ORCHIDEE land surface model (version 2.0).

The atmospheric general circulation model LMDZ6A-LR builds onto its previous version that has notably incorporated advances in the parameterization of turbulence, convection, and clouds. More specifically, LMDZ6A-LR includes a turbulent scheme based on the prognostic equation for the turbulent kinetic energy that follows Yamada (1983), a mass flux representation of the organized structures of the convective boundary layer called the thermal plume model (Hourdin et al., 2002; Rio and Hourdin, 2008; Rio et al., 2010), and a parameterization of the cold pools or wakes created below cumulonimbus by the evaporation of convective rainfall (Grandpeix and Lafore, 2010; Grandpeix et al., 2010). It is based on a regular horizontal grid with 144 grid points regularly spaced in longitude and 142 in latitude, corresponding to a resolution of $\mathrm{2.5}{}^{\circ }\times \mathrm{1.3}{}^{\circ }$, and 79 vertical layers.

IPSL-CM6A-LR further includes NEMO (Nucleus for European Models of the Ocean), which is itself composed of three major building blocks: the ocean physics NEMO-OPA (Madec et al., 2016), the sea ice dynamics and thermodynamics NEMO-LIM3 (Vancoppenolle et al., 2009; Rousset et al., 2015), and the ocean biogeochemistry NEMO-PISCES (Aumont et al., 2015). The grid used has a nominal resolution of 1^∘ in the zonal and meridional directions with a latitudinal grid refinement of $\mathrm{1}/\mathrm{3}{}^{\circ }$ in the tropics. Vertical discretization uses a partial step formulation (Bernard et al., 2006), which ensures a better representation of bottom bathymetry, with 75 levels. The initial layer thicknesses increase nonuniformly from 1 m at the surface to 10 m at 100 m depth and reach 200 m at the bottom and are subsequently time-dependent. NEMO-PISCES (Aumont et al., 2015) models the lower trophic levels of the marine ecosystem (phytoplankton, microzooplankton, and mesozooplankton) and the biogeochemical cycles of carbon and of the main nutrients (P, N, Fe, and Si). This model also computes air–sea carbon fluxes.

Finally, IPSL-CM6A-LR includes ORCHIDEE, a global process-based terrestrial biosphere model (Krinner et al., 2005) that calculates carbon, water, and energy fluxes between the land surface and the atmosphere. Photosynthesis and all components of the surface energy and water budgets are calculated at a 15 min resolution, while the dynamics of the carbon storage (including carbon allocation in plant reservoirs, soil carbon dynamics, and litter decomposition) are resolved on a daily basis. Photosynthesis depends on light availability and CO₂ concentration, soil moisture, and temperature and is parameterized based on Farquhar et al. (1980) and Collatz et al. (1992) for C₃ and C₄ plants, respectively. In the absence of an explicit representation of the nitrogen cycle, this latest version of ORCHIDEE includes a downregulation capability that models a reduction in the terrestrial photosynthesis rates as a function of CO₂ concentration. While the N and C cycles interact in multiple ways, a key aspect of these interactions is the limitation of carbon uptake by nitrogen availability, especially under increasing atmospheric CO₂ concentrations. The downregulation parameterization models a reduction in the maximum photosynthetic rate as the atmospheric CO₂ concentrations increase using a logarithmic function of the CO₂ concentration relative to 380 ppm (Sellers et al., 1996). The PFT-dependent parameters of this parameterization are chosen to broadly reproduce the change in GPP observed at the free-air CO₂ enrichment (FACE) experiment sites (Norby and Zak, 2011). In ORCHIDEE, the spatial distribution of vegetation is represented using 15 plant functional types (PFTs; Prentice et al., 1992; Cramer, 1997; Wullschleger et al., 2014). More precisely these PFTs are decomposed into three groups according to their physiological behaviour under similar climate conditions: tall vegetation (forests) represented by eight PFTs, short vegetation (grasses and crops) represented by six PFTs, and bare soil. The fractional coverage of these PFTs varies geographically. A soil type is associated with each one of these three PFT groups. This three-group partitioning allows for the dividing of each grid box into three tiles for which an independent hydrological budget is calculated, using the 11-layer physically based hydrology scheme. In ORCHIDEE the wood harvest product from the LUHv2h database is used in addition to the annual land cover maps.

A4.7 Team MIROC (Japan Agency for Marine-Earth Science and Technology, the University of Tokyo, the National Institute for Environmental Studies) MIROC-ES2L

MIROC-ES2L (Hajima et al., 2020) is based on the global climate model MIROC5.2 (Tatebe et al., 2018), which is a minor updated version of MIROC5 used for CMIP5 (Watanabe et al., 2010). The physical core shares almost the same structure and characteristics with the latest model MIROC6 (Tatebe et al., 2019), except for the atmospheric spatial resolution and treatment of cumulus clouds. This model interactively couples an atmospheric general circulation model (CCSR-NIES AGCM; Tatebe et al., 2019) including an on-line aerosol component (SPRINTARS; Takemura et al., 2000), an ocean general circulation model (GCM) with a sea ice component (COCO; Hasumi, 2015), and a land physical surface model (MATSIRO; Takata et al., 2003). The land and ocean biogeochemical components are represented by VISIT (Ito and Inatomi, 2012) and OECO2 (Hajima et al., 2020), respectively, which are interactively coupled to the atmospheric component. There exists another branched version that has an atmospheric chemistry component with a finer atmospheric grid (MIROC-ES2H), but it was not used in this study.

The atmospheric grid resolution is approximately 2.81^∘ with 40 vertical levels between the surface and about 3 hPa. For the ocean, the model employs a tripolar coordinate system with 62 vertical levels. To the south of 63^∘ N, the ocean model has longitudinal grid spacing of about 1^∘, while the meridional grid spacing varies from about 0.5^∘ near the Equator to 1^∘ in the midlatitudes. Over the Arctic ocean the grid resolution is even finer following the tripolar coordinate system. The physical terrestrial component resolves the vertical soil profile with six layers down to a 14 m depth, with two types of land use tiles (agriculture and nonagriculture). Terrestrial biogeochemical component considers two layers of soil organic matter (the upper litter layer and the lower humus layer), with five types of land use tiles (primary vegetation, secondary vegetation, urban, crop, and pasture).

The terrestrial biogeochemical component covers major processes relevant to the global carbon cycle, with vegetation (leaf, stem, and root), litter (leaf, stem, and root), and humus (active, intermediate, and passive) pools and with a static biome distribution. Details on carbon cycle processes in the model can been found in Ito and Oikawa (2002). The N cycle is simulated with N pools of vegetation (canopy and structural), organic soil (litter, humus, and microbe), and inorganic nitrogen (ammonium and nitrate). The model considers two major nitrogen influxes into the ecosystem (biological nitrogen fixation and external nitrogen inputs). Fluxes out of the land ecosystem in the model are N₂ or N₂O emissions, leaching, NH₃ emissions, and other emission-like volatilization from land use product pools. For installing into MIROC-ES2L, the terrestrial ecosystem processes were modified such that photosynthetic capacity is controlled by leaf N concentration. Processes associated with land use change are also modified to take full advantage of the CMIP6 LUC forcing dataset. Further details can be found in Hajima et al. (2020).

The new ocean biogeochemical component model, OECO2, is a NPZD-type model and modified from the previous model (Watanabe et al., 2011). The biogeochemical compartments of OECO2 are nitrate, phosphate, dissolved iron, dissolved oxygen, two types of phytoplankton (nondiazotroph and diazotroph), zooplankton, and particulate detritus. There exist other compartments of dissolved inorganic carbon (DIC), total alkalinity, calcium, calcium carbonate, and N₂O. All organic materials have an identical elemental stoichiometric ratio. The model considers external nutrient inputs (atmospheric N and Fe deposition, inorganic N and P from rivers, biological N fixation, Fe input from ocean bottom and shelf) and nutrient loss (denitrification for N and loss into sediment for N, P, and Fe). The emission, transportation, and deposition processes of iron are explicitly simulated by the atmospheric aerosol component.

A4.8 Max Planck Institute for Meteorology (MPI) MPI-ESM1.2-LR

The MPI-ESM1.2-LR model (Mauritsen et al., 2019) consists of ocean, atmosphere, land, and sea ice components which are connected via a coupler analogous to the predecessor MPI-ESM versions (Giorgetta et al., 2013). The atmosphere model, ECHAM6.3, at the LR resolution has a spectral truncation at T63 or approximately 200 km grid spacing with 47 vertical levels. It is directly coupled to the land model, JSBACH3.2, through surface exchange of mass, momentum, and heat. The ocean general circulation model, MPIOM1.6 in MPI-ESM1.2-LR, runs on a bipolar grid GR1.5 and has 40 unevenly placed levels. It computes transport of tracers of the ocean biogeochemistry model HAMOCC6 (Ilyina et al., 2013; Paulsen et al., 2017). The MPI-ESM-LR configuration computes 45–85 model years per physical day enabling new simulations which were not feasible previously, such as for instance, large ensemble simulations (Maher et al., 2019) or millennial-scale simulations with an interactive carbon cycle (Brovkin et al., 2019).

Terrestrial vegetation in JSBACH includes vegetation dynamics which interacts with land use changes (Reick et al., 2013), accounting for the latest changes in the land use harmonization dataset by Hurtt et al. (2006). The new SPITFIRE model simulates burned area and carbon emissions to the atmosphere due to wildfires and anthropogenic fires (Lasslop et al., 2014), replacing the old global fire parameterization used in the CMIP5 model. The soil carbon model Yasso simulates the dynamics of four fast soil carbon pools, which are different for leaf and woody litter types, plus a slow humus pool (Goll et al., 2015). Nitrogen and carbon pools are coupled based on CO₂-induced nitrogen limitation (Goll et al., 2017).

The ocean biogeochemistry model HAMOCC6 has been extended compared to the previous version described in Ilyina et al. (2013) to explicitly resolve nitrogen-fixing cyanobacteria as an additional prognostic phytoplankton class (Paulsen et al., 2017). This allows for the capture of the response of N₂ fixation and ocean biogeochemistry to changing climate conditions. Additionally, updates of existing processes have been performed. This includes for instance the addition of a vertically varying settling rate for detritus following the formulation by Martin et al. (1987). Finally some empirical relationships in the parameterized processes have been updated to follow recommendations of the C⁴MIP and OMIP protocols (Jones et al., 2016; Orr et al., 2017). The full overview of changes in HAMOCC is given in Mauritsen et al. (2019).

A4.9 Geophysical Fluid Dynamics Laboratory (GFDL) NOAA-GFDL-ESM4

GFDL-ESM4.1 is a comprehensive, fully coupled Earth system model developed by NOAA's Geophysical Fluid Dynamics Laboratory with a fully interactive carbon cycle and interactive atmospheric chemistry (Dunne et al., 2020) that builds on previous-generation modelling efforts of the carbon cycle (ESM2 series; Dunne et al., 2012, 2013) and atmospheric chemistry (CM3; Donner et al., 2011) along with increased resolution and improved numerics and physics, akin to the GFDL's fourth-generation coupled climate model (CM4.0; Held et al., 2019), and representation of additional Earth system processes.

The atmospheric component, GFDL AM4.1, is based on the third-generation finite-volume cube-sphere dynamical core (FV3; Lin, 2004) with a ∼100 km horizontal resolution and 49 vertical levels. The model top is located at ∼0.01 hPa to resolve the stratosphere. AM4.1 shares the critical developments in model physics with the AM4.0 model (Zhao et al., 2018), including radiation, convection, and clouds. AM4.1 differs from the AM4.0 model in its enhanced vertical resolution and its more explicit representation of atmospheric chemistry that motivated a separate radiative and gravity wave tuning.

AM4.1 includes interactive tropospheric and stratospheric gas-phase and aerosol chemistry represented through 77 prognostic (transported) tracers and 41 diagnostic (nontransported) chemical tracers. The tropospheric chemistry includes reactions for the oxidation of methane among other volatile organic compounds. The stratospheric chemistry accounts for the major ozone loss cycles and heterogeneous reactions on liquid and solid stratospheric aerosols.

Land hydrology and ecosystem dynamics are represented by the GFDL Land Model version 4.1 (LM4.1, Shevliakova et al., 2020) and build on the previous-generation LM3.1 model (Milly et al., 2014). Soil carbon dynamics and biogeochemistry are represented through the CORPSE model (Sulman et al., 2019) with an explicit treatment of soil microbes. LM4.1 also includes a new fire model FINAL (Rabin et al., 2018). Vegetation dynamics are represented by the second-generation age–height-structured approach, the perfect plasticity approximation (PPA; Weng et al., 2015; Martinez Cano et al., 2020). Allometric constraints and competition enable the simulation of size structure and carbon fluxes in dynamic vegetation. There are six carbon pools in LM4.1 representing leaves, fine roots, heartwood, sapwood, seeds, and nonstructural carbon. Litter is broken into leaf, fine roots, and coarse-wood categories. Soil has 20 vertical levels each with its own prognostic state for energy, water, and soil carbon variables. There are five types of vegetation in LM4.1 representing C₃ grass, C₄ grass, tropical trees, temperate deciduous trees, and cold evergreen trees. A combination of these vegetation types could coexist. The model also includes a new treatment of stomatal conductance and plant hydraulics. The vegetation state is used to drive a dust emission model that is coupled with the atmosphere for transport. The ESM4 implementation of LM4.1 does not include an interactive nitrogen cycle.

The ocean biogeochemical component of ESM4 is version 2 of the Carbon, Ocean Biogeochemistry and Lower Trophics (COBALTv2) model (Stock et al., 2014b). COBALTv2 uses 33 tracers to represent carbon, alkalinity, oxygen, nitrogen, phosphorus, iron, silica, calcite, and lithogenic mineral cycling within the ocean. Relative to previous-generation ocean biogeochemistry models developed at the GFDL, COBALTv2 includes an enhanced representation of plankton food web dynamics to resolve the flow of energy from phytoplankton to fish (Stock et al., 2014a) and enhance the model's capacity to resolve linkages between food webs and biogeochemical cycles. COBALTv2 explicitly includes small, large (split into diatoms and nondiatoms), and diazotrophic phytoplankton groups; three zooplankton groups; bacteria; and three labilities of dissolved organic matter. Other updates include a temperature dependence for sinking organic matter remineralization (Laufkötter et al., 2017), the addition of semilabile dissolved organic material and carbonate chemistry calculations based on the open-source mocsy 2.0 (Orr and Epitalon, 2015).

Data from the NOAA-GFDL-ESM4 model used in the analysis presented in this paper are accessible via the Earth System Grid Federation (ESGF) for the 1pctCO2 (Krasting et al., 2019b) simulation and for its radiatively and biogeochemically coupled configurations (Krasting et al., 2019a).

A4.10 Norwegian Climate Centre (NCC) NorESM2-LM

The NorESM2-LM (Seland et al., 2020) is based on the latest release of the Community Earth System Model (CESM2.1), whose development is supported by the National Center for Atmospheric Research in the United States. NorESM2 keeps the original land and sea ice components of CESM2.1 (i.e. CLM5 and CICE5, respectively). The atmospheric component is CAM6 (as in CESM) but with modifications regarding the energy and angular momentum conservation. Further, the atmospheric aerosol module of CAM6 has been replaced by the scheme developed by the Norwegian Meteorological Institute. The ocean physical and biogeochemical components of NorESM2 are the isopycnal ocean circulation and carbon cycle components updated from NorESM1 version (Tjiputra et al., 2013; Schwinger et al., 2016)

The CLM5 (Community Land Model version 5) prognostically simulates the carbon and nitrogen cycles, which include natural vegetation, crops, and soil biogeochemistry. The carbon and nitrogen budgets comprise leaf, live-stem, dead-stem, live-coarse-root, dead-coarse-root, fine-root, and grain pools. Each of these pools has short-term and long-term storage of nonstructural carbohydrates and labile nitrogen. In addition to the vegetation pools, CLM includes a series of decomposing carbon and nitrogen pools as vegetation successively breaks down into coarse woody debris and/or litter and subsequently into soil organic matter. Details on the CLM5 models are available in Lawrence et al. (2019).

Similar to the earlier version, the ocean carbon cycle component in NorESM2 is based on the Hamburg Ocean Carbon Cycle (HAMOCC; Maier-Reimer et al., 2005) model, which has been coupled to an isopycnic ocean general circulation model. The current version includes new processes and refined parameterizations, as well as new diagnostic tracers. The ecosystem model is based on an NPZD-type model with multinutrient limitation in its phytoplankton growth formulation. Riverine fluxes of inorganic and organic carbon as well as nutrients are now implemented. Unlike the earlier version, the sea-to-air dimethyl sulfate (DMS) fluxes alter the atmospheric radiative forcing and hence the climate–carbon cycle feedback. More details on the ocean carbon cycle of NorESM2 are available in Tjiputra et al. (2020).

A4.11 The United Kingdom Community Earth System Model, UKESM1-0-LL

UKESM1-0-LL (Sellar et al., 2019) is based upon the HadGEM3-GC3.1 (Williams et al., 2018) global climate model which includes coupled ocean, atmosphere, land, and sea ice components. The atmosphere component is the Unified Model with a resolution of 1.875^∘ by 1.25^∘ with 85 vertical levels up to a model top of 90 km (Walters et al., 2019) and includes a modal aerosol scheme (Mann et al., 2010). The ocean component uses the NEMO dynamical ocean at 1^∘ resolution with 75 vertical levels (Storkey et al., 2018). The sea ice component uses CICE on the same grid as the ocean with five ice thickness categories (Ridley et al., 2018). The land component uses the JULES land surface model (Wiltshire et al., 2020); however, the land surface configuration is substantially updated for UKESM. The primary differences between the physical and Earth system models are the inclusion of a terrestrial carbon and nitrogen cycle (Wiltshire et al., 2020), ocean biogeochemistry (Yool et al., 2013), and a tropospheric–stratospheric chemistry model. Atmospheric chemistry in UKESM1 is simulated by the UKCA chemistry and aerosol model with the specific configuration being a combination of tropospheric (O'Connor et al., 2014) and stratospheric chemistry (Morgenstern et al., 2009, 2017).

Terrestrial biogeochemistry is represented by the JULES-ES model (Wiltshire et al., 2020). The land surface is represented by 13 plant functional types (PFTs) including 4 managed crop and pasture land types. The height, leaf area index, and spatial distribution of the PFTs are dynamically simulated by the TRIFFID dynamic global vegetation model (DGVM; Cox, 2001). Soil carbon is represented by the four-pool RothC scheme (Coleman and Jenkinson, 1999). Terrestrial carbon uptake may be limited by the availability of nitrogen. Nitrogen does not directly affect photosynthetic capacity through leaf N concentrations but acts indirectly by controlling the biomass and leaf area index within the TRIFFID DGVM. A second mechanism acts through soil carbon by limiting the decomposition of litter into soil carbon in the RothC model. The vegetation model includes the retranslocation of nitrogen during the senescence of leaves and roots into a labile pool to supply nutrients for the following seasonal leaf-out. The soil model simulates mineralization and immobilization with mineralized nitrogen becoming available for plant uptake and ecosystem loss. Inorganic nitrogen is represented by a single grid box pool to which all PFTs have equal access. Nitrogen deposition is prescribed from ancillary data.

Land use change is represented by the application of time-varying fields of crop and pasture to the DGVM, which allocates space dynamically to C₃ and C₄, and crop and pasture types. Pasture is represented as natural grass, whereas crops include a harvest parameterization and are fertilized. Biogenic volatile organic compound (BVOC) emissions from vegetation are simulated and affect the formation of secondary organic aerosols. Mineral dust is emitted from bare soil and acts as both an aerosol and a fertilizer to the ocean.

Ocean biogeochemistry is represented by MEDUSA-2 (Yool et al., 2013), which resolves a dual size-structured ecosystem of small (nanophytoplankton and microzooplankton) and large (microphytoplankton and mesozooplankton) components. This explicitly includes the biogeochemical cycles of nitrogen, silicon, and iron nutrients as well as the cycles of carbon, alkalinity, and dissolved oxygen. Large phytoplankton are treated as diatoms and utilize silicic acid in addition to nitrogen, iron, and carbon. Like the living components, the detrital components are split into two size classes. At the seafloor, MEDUSA-2 resolves five reservoirs to temporarily store sinking organic material reaching the sediment. The model's nitrogen, silicon, and alkalinity cycles are closed and conservative (e.g. no riverine inputs), while the other three cycles (carbon, iron, oxygen) are open. The ocean's iron cycle includes aeolian (land-derived dust) and benthic sources, and is depleted by scavenging. The ocean's carbon cycle exchanges CO₂ with the atmosphere. The ocean's oxygen cycle exchanges oxygen with the atmosphere, and dissolved oxygen is additionally created by primary production and depleted by remineralization. Ocean biogeochemistry also feeds back on the atmosphere through the production of marine DMS and marine organic aerosols.

A5 Contribution of uncertainties in $\mathrm{\Delta }{T}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ and $E_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ to the TCRE.

The uncertainty in the TCRE, as indicated by its standard deviation (σ_TCRE), can be represented in terms of the standard deviation of $\mathrm{\Delta }{T}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ (σ_ΔT), the standard deviation of ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ (σ_E), and their means $\stackrel{\mathrm{‾}}{\mathrm{\Delta }T}$ and $\stackrel{\mathrm{‾}}{E}$ across the 11 CMIP6 models. Since $\mathrm{\Delta }{T}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ and ${\stackrel{\mathrm{̃}}{E}}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ are nearly independent (correlation between these two quantities is only 0.02 across the 11 CMIP6 models considered here), we can write

$$\begin{array}{}\text{(A8)}& {\mathit{\sigma }}_{\text{TCRE}}=\stackrel{\mathrm{‾}}{\text{TCRE}}\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{\Delta }T}}{\stackrel{\mathrm{‾}}{\mathrm{\Delta }T}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{E}}}{\stackrel{\mathrm{‾}}{E}}}\right)}^{\mathrm{2}}},\end{array}$$

which allows us to calculate contributions of ${\left(\frac{{\mathit{\sigma }}_{\mathrm{\Delta }T}}{\stackrel{\mathrm{‾}}{\mathrm{\Delta }T}}\right)}^{\mathrm{2}}$ and ${\left(\frac{{\mathit{\sigma }}_{\mathrm{E}}}{\stackrel{\mathrm{‾}}{E}}\right)}^{\mathrm{2}}$ to σ_TCRE.