# A small climate-amplifying effect of climate-carbon cycle feedback

Xuanze Zhang, Ying-Ping Wang, Peter Rayner, Philippe Ciais, Kun Huang, Yiqi Luo, Shilong Piao, Zhonglei Wang, Jianyang Xia, Wei Zhao, et al.

► To cite this version:

Xuanze Zhang, Ying-Ping Wang, Peter Rayner, Philippe Ciais, Kun Huang, et al. A small climate-amplifying effect of climate-carbon cycle feedback. Nature Communications, 2021, 12 (1), 10.1038/s41467-021-22392-w. hal-03231668

HAL Id: hal-03231668

https://hal.science/hal-03231668

Submitted on 26 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-021-22392-w

<figure><img src="image_3.png" /></figure>

Check for updates

OPEN

# A small climate-amplifying effect of climate-carbon cycle feedback

Xuanze Zhang<sup>1,2✉</sup>, Ying-Ping Wang<sup>3,4✉</sup>, Peter J. Rayner<sup>5</sup>, Philippe Ciais<sup>6</sup>, Kun Huang<sup>2</sup>, Yiqi Luo<sup>7</sup>, Shilong Piao<sup>8</sup>, Zhonglei Wang<sup>9</sup>, Jianyang Xia<sup>2</sup>, Wei Zhao<sup>10</sup>, Xiaogu Zheng<sup>11</sup>, Jing Tian<sup>1</sup> & Yongqiang Zhang<sup>1✉</sup>

The climate-carbon cycle feedback is one of the most important climate-amplifying feedbacks of the Earth system, and is quantified as a function of carbon-concentration feedback parameter (β) and carbon-climate feedback parameter (γ). However, the global climate-amplifying effect from this feedback loop (determined by the gain factor, g) has not been quantified from observations. Here we apply a Fourier analysis-based carbon cycle feedback framework to the reconstructed records from 1850 to 2017 and 1000 to 1850 to estimate β and γ. We show that the β-feedback varies by less than 10% with an average of 3.22 ± 0.32 GtC ppm⁻<fnref n="1" /> for 1880–2017, whereas the γ-feedback increases from −33 ± 14 GtC K⁻<fnref n="1" /> on a decadal scale to −122 ± 60 GtC K⁻<fnref n="1" /> on a centennial scale for 1000–1850. Feedback analysis further reveals that the current amplification effect from the carbon cycle feedback is small (g is 0.01 ± 0.05), which is much lower than the estimates by the advanced Earth system models (g is 0.09 ± 0.04 for the historical period and is 0.15 ± 0.08 for the RCP8.5 scenario), implying that the future allowable CO₂ emissions could be 9 ± 7% more. Therefore, our findings provide new insights about the strength of climate-carbon cycle feedback and about observational constraints on models for projecting future climate.

<footnote n="1">¹Key Laboratory of Water Cycle and Related Land Surface Processes, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing, China. ²Research Center for Global Change and Ecological Forecasting, School of Ecological and Environmental Science, East China Normal University, Shanghai, China. ³Terrestrial Biogeochemistry Group, South China Botanical Garden, Chinese Academy of Sciences, Guangzhou, China. ⁴CSIRO Oceans and Atmosphere, Private Bag 1, Aspendale, Victoria, Australia. ⁵School of Earth Sciences, Climate and Energy College, University of Melbourne, Parkville, Victoria, Australia. ⁶Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France. ⁷Center for Ecosystem Science and Society, Northern Arizona University, Flagstaff, AZ, USA. ⁸Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University, Beijing, China. ⁹Wang Yanan Institute for Studies in Economics (WISE) and School of Economics, Xiamen University, Xiamen, China. ¹⁰National Meteorological Center, China Meteorological Administration, Beijing, China. ¹¹Key Laboratory of Regional Climate-Environment Research for East Asia, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing, China.</footnote>

<footnote>✉email: xuanzezhang@ignsnrr.ac.cn; yingping.wang@csiro.au; zhangyq@ignsnrr.ac.cn</footnote>
With increasing atmospheric CO₂ and a warming climate during the industrial era, land and ocean reservoirs have together absorbed >50% of anthropogenic CO₂ emissions¹, playing a significant role in reducing anthropogenic warming². Whether the ocean and land sinks will continue to take up a similar fraction in the future remains uncertain, and one of the major causes for that uncertainty is the feedback between the carbon cycle and the physical climate system³⁻⁵. Global climate change affects carbon uptake by land and oceans, which impacts the rate of increase in atmospheric CO₂ and, in turn, climate change. This feedback loop between the physical climate system and the global carbon cycle of the Earth system was quantified using a modeling approach⁶⁻⁹.

Previous studies quantified the climate-carbon cycle feedback as a function of the carbon-concentration feedback response parameter (β) and the carbon-climate feedback response parameter (γ)⁷,¹⁰. The β (GtC ppm⁻¹) and γ (GtC K⁻¹) are also defined as the rates of change in land and ocean carbon storages relative to a fixed reference time to atmospheric CO₂ concentration increase and to global climate change that is often quantified by the global-mean surface temperature change, respectively. From the perspective of the land and ocean reservoirs, β is positive, and γ is negative. Therefore, β-feedback reduces the impact of CO₂ emissions on atmospheric CO₂ concentrations and then global warming, while γ-feedback amplifies global warming. The combined effects of the β-feedback and the γ-feedback and the nonlinear interaction between them determine the strength of the climate-carbon cycle feedback loop which is known as the feedback gain factor (g).

Both β and γ were previously quantified using physical climate models coupled with the global carbon cycle under idealized experiments⁷,¹⁰,¹¹. Based on the experiments by the advanced Earth system models, the last two successive assessments by the International Panel on Climate Change (IPCC) found that the uncertainty of the climate-carbon cycle feedback was dominated by β¹²,¹³. These Earth system models under idealized experiments did not account for the biophysical effect of land-use change on climate and ecological effects on the residence time of carbon in the land biosphere. Modeling experiments showed that in a high CO₂-induced warming climate system, the nonlinear carbon-climate feedback can reduce the ocean carbon uptake by 3.6–10.6% based on the simulations of seven Earth system models¹⁴.

On the other hand, observations can be used to constrain the estimated climate-carbon cycle feedbacks. Applying an emerging constraint-based approach to instrumental records, Cox et al.¹⁵ estimated γ for tropical land to be −53 ± 17 GtC K⁻¹ by 2100. Using data from three ice cores and multiple temperature reconstructions from tree ring data over 1050–1800, Frank et al.¹⁶ reported a ΔCₐ/ΔTₐ (defined as η in this study) of 7.7 ppm K⁻¹ with a likely range of 1.7–21.4 ppm K⁻¹, which was much lower than previous estimates (e.g., ~40 ppm K⁻¹)¹⁷,¹⁸. However, the η is not truly comparable with γ, because their relationship also depends on β (see Eq. (3) below). In quantifying the response of climate to an increase in atmospheric CO₂ concentration or anthropogenic emissions, two other quantities are frequently used in the studies of carbon-climate interactions: the sensitivity of climate to atmospheric CO₂ (α ≡ ΔTₐ/ΔCₐ, in K GtC⁻¹, which is equivalent to 1/(2.12 η) in this study)⁷,¹⁹, and the transient climate response to cumulative CO₂ emission (TCRE ≡ ΔTₐ/ΔCₐ, in K GtC⁻¹)²⁰. The values of β, γ, α, and TCRE are related to each other through one equation (see Eq. (7))¹⁰,²¹, and that equation is applicable to any individual frequency (or timescale) (see “Methods”).

## Results

**Estimates of β and γ across timescales.** Based on the previous studies⁷,¹⁰,¹⁴,¹⁹, the climate-carbon cycle feedback framework with a nonlinear feedback term as a function of β and γ parameters, i.e., f(β, γ) in a unit of GtC ppm⁻¹ K⁻¹ or GtC GtC⁻¹ K⁻¹, in a CO₂ emission-driven coupled climate-carbon cycle system, can be expressed as:

$$ \Delta C_E = \Delta C_A + \beta \Delta C_A + \gamma \Delta T_A + f(\beta, \gamma) \Delta C_A \Delta T_A \quad (1) $$

where ΔCₐ, ΔCₐ, and ΔTₐ are changes in cumulative CO₂ emissions, atmospheric CO₂, and global surface temperature, respectively, during a time interval (Δt). Here, we applied a Fourier analysis-based spectral decomposition approach to observations or reconstructed records of global surface temperature and atmospheric CO₂ to quantify how the two feedback parameters vary over different periods of time, or across different timescales (1/frequencies) over the same time period. As shown in the “Methods” section, the relationship among the four parameters at any given timescale (k) is given by

$$ (TCRE^{-1})_k = (1 + \beta_k)(\alpha^{-1})_k + \gamma_k^* \quad (2) $$

where γ* = γ + f(β, γ)ΔCₐ represents the combined carbon-climate feedback from the linear γ-feedback and the nonlinear impact of atmospheric CO₂ on the carbon-climate feedback (see “Methods”). Equation (2) also shows that variations of TCRE⁻¹ and α⁻¹ for a given timescale are linearly related, with a slope of 1 + βₖ and an intercept of γₖ*. Our result from modeling experiments (see Supplementary Table 4) showed that the nonlinear feedback term f(β, γ) is relatively small, and its contribution (f(β, γ)ΔCₐ) to the γ* is 15 ± 23%, while its contribution (f(β, γ)ΔTₐ) to the β-feedback is negligible (3 ± 3%). Thus, we indicate that the nonlinear feedback term has negligible effect on the estimate of the slope of 1 + βₖ.

We estimated (TCRE⁻¹)ₖ from the ratio of the variability in Cₐ to the variability in Tₐ, and (α⁻¹)ₖ from the ratio of the variability in Cₐ to the variability in Tₐ on a given timescale (see “Methods”). In Fourier analysis, the variables of Cₐ, Cₐ, and Tₐ were represented by the amplitudes of their harmonics at a given frequency or period (Supplementary Fig. 1), which were then used to calculate TCRE⁻¹ and α⁻¹ across different timescales (see “Methods”). Different from the previous approaches based on a Taylor series expansion and modeling experiments (hereafter referred to as the FEA approach, see Eqs. (22–25) in “Methods”)⁷,¹⁰,¹⁴, the estimated β and γ using Fourier analysis-based approach do not depend on a reference time that was used to compute changes in ΔCₐ, ΔCₐ, and ΔTₐ as in the FEA approach.

In the following, we apply the Fourier analysis-based approach to observed time series of Cₐ, Cₐ, and Tₐ to estimate βₖ and γₖ for the industrial period (1850–2017). By ignoring the timescale dependence of β, and using the value of β for the industrial period, we estimated γₖ at different timescales using large ensembles based on combinations of ice-core atmospheric CO₂ records (Cₐ) and reconstructed surface air temperature (Tₐ) datasets (“Methods”) over the period 1000–1850. Finally, we compare the observation-based estimates of β and γ, and corresponding feedback gain factor (g) with those estimated from the simulations from multiple Earth system models.

**Estimates of β and γ* for the industrial period (1850–2017).** Figure 1 shows the observed increase in atmospheric CO₂ (ΔCₐ) from ~285 ppm in 1850 to 405 ppm in 2017 in response to the increase in cumulative CO₂ emissions (ΔCₐ), including fossil fuel combustion²² and land-use change²³,²⁴, which amounted to ~630 ± 42 GtC by 2017 since 1850, with respect to a concurrent
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 1** Variability in atmospheric CO₂ concentration and temperature over the past millennium. **a** Variations in atmospheric CO₂ concentrations over 1000–1850 from three Antarctic ice-core records at Law Dome, WAIS Divide, and Dronning Maud Land. The subplot in **a** shows changes in atmospheric CO₂ concentrations and cumulative anthropogenic CO₂ emissions over 1850–2017 with respect to values in 1850. The cumulative anthropogenic CO₂ emissions with an uncertainty of ±1σ confidence interval (red-shaded) include CO₂ release from both fossil fuel combustion and land-use change (unit was converted to ppm from GtC). **b** Variations in northern hemispheric-mean temperature anomalies from five reconstructions (PAGES2k, Frank2010, Mann2009, MannEIV, and Moberg2005). The subplot in **b** shows global-mean temperature anomalies calculated from land surface air temperature and sea surface temperature with respect to the average of 1961–1990 from datasets of HadCRUT4 and Berkeley Earth for 1850–2017, and GISTEMP and NOAA GlobalTemp for 1880–2017. The reconstructions of 1000–1850 in (a) and (b) were smoothed with 30-year splines. Temperature anomalies in 1850 (1880)–2017 in (b) were smoothed with 5-year splines.</figcaption></figure>

increase in global-mean annual temperature ($\Delta T_A$) of ~1.1 K over the same period based on the averages of the four global temperature datasets. Using Fourier analysis, we found that the amplitudes of $C_E$, $C_A$, and $T_A$ for the industrial period increased nonlinearly with timescale (Fig. 2a and b), but the variations in $(TCRE^{-1})_k$ and $(\alpha^{-1})_k$ with timescale were linearly related to each other ($R^2 = 0.99$, $P < 0.001$) with a slope (mean ± 1 standard deviation) of 2.52 ± 0.15 (see Fig. 2c and Supplementary Table 1 for individual estimates). Across different timescales from inter-annual to multi-decadal (e.g., $k$ from 2 to 90 years) for the 1880–2017 period, $(\alpha^{-1})_k$ varied from 50 to 250 GtC K⁻¹ and $(TCRE^{-1})_k$ from 125 to 600 GtC K⁻¹ (Fig. 2c).

The slope of the linear regression between $(TCRE^{-1})_k$ and $(\alpha^{-1})_k$ across different timescales can be used to estimate $(1 + \beta_k)$ based on Eq. (2). The result in Fig. 2c shows that the value of $\beta_k$ is approximately constant across different timescales, with an average value of 3.22 ± 0.32 GtC ppm⁻¹ (or 1.52 ± 0.15 GtC GtC⁻¹) throughout the industrial period in 1880–2017 (Supplementary Table 1). The uncertainty in $\beta$ results from the uncertainties in $C_E$ from fossil fuel combustion and land-use change estimates and

errors in the temperature datasets and the regression (see “Methods”). Using the estimated $\beta$ for the industrial period, we estimated the corresponding $\gamma^*$ to be −10.9 ± 3.6 GtC K⁻¹ for the period 1880–2017 based on Eq. (6).

To diagnose the constancy of $\beta$ estimated from the industrial period and the reliability of our Fourier analysis-based approach, we used a box model to predict changes in temperature and CO₂ concentration with CO₂ emissions as forcing input (see “Methods”, Eqs. (29)–(30)). We ran the box model using the estimate of two parameters based on observations (i.e., $\beta = 3.22$ GtC ppm⁻¹ and $\gamma^* = -10.9$ GtC K⁻¹) and the annual cumulative emissions over 1850–2017. The model predictions fitted very well to the observed values of atmospheric CO₂ concentration and surface temperature ($R^2 = 0.99$, RMSE = 3.5 ppm for CO₂, and $R^2 = 0.96$, RMSE = 0.17 K for temperature) (see Supplementary Fig. 2). The greater trend of the predicted global surface temperature relative to the observed after 1980s (Supplementary Fig. 2) is likely associated with $\gamma^*$ being treated as a constant, while in reality $\gamma^*$ could vary with internal climate variability and across different timescales. Using the box model predicted annual CO₂ and temperature, we in return, estimated a nearly constant $\beta$ (3.42 GtC ppm⁻¹) and a $\gamma^* (-12.2$ GtC K⁻¹) for the period 1880–2017 based on our approach of Eqs. (1) and (2). Both values of $\beta$ and $\gamma^*$ estimated from the predicted CO₂ and temperature time series by the box model fall within the uncertainties of observation-based $\beta$ and $\gamma^*$, suggesting that the $\beta$ for 1880–2017 being nearly constant is robust.

As the $\gamma^*$ for the industrial period consists of the $\gamma$-feedback and the nonlinear feedback effect of $f(\beta, \gamma)\Delta C_A$, which results from the multi-decadal climate variability and anthropogenic CO₂ emissions. Because of the relatively short records over the industrial period, we could not separate the $\gamma$-feedback parameter directly from the $\gamma^*$ using our analysis framework. The $\eta$ (or $\alpha^{-1}$) for the industrial period includes a possibly significant contribution from an emissions-driven increase in atmospheric CO₂ concentration to the climate-carbon feedback. As a result, the mean value of $\eta$ for 1850–2017 is 109 ppm K⁻¹, much higher than the estimates of 7.7–40 ppm K⁻¹ during the preindustrial period before 1850s¹⁶,¹⁷ that is considered to be at quasi-equilibrium with small variation in greenhouse gas forcing. The reason for this large discrepancy in $\eta$ between the preindustrial and industrial periods is likely a consequence of the nonlinear dependence of radiative forcing on atmospheric CO₂ concentrations²⁵, and the temperature change in response to the increase in atmospheric CO₂ during the industrial era has not reached steady state, as large part of atmospheric CO₂ increase during industrial period was driven by emissions not due to warming-induced CO₂ release from land and ocean reservoirs. It is also well known that equilibrium climate sensitivity is often considerably larger than the transient climate sensitivity²⁶. In the following, we used reconstructed records of atmospheric CO₂ and surface temperature during the preindustrial last millennium from 1000 to 1850¹⁶⁻¹⁸,²⁷,²⁸ to estimate the $\gamma$.

**Timescale dependence of $\gamma$ over the preindustrial last millennium (1000–1850).** The preindustrial last millennium (1000–1850) climate could be considered as a quasi-equilibrium state with very little change in CO₂ emissions, therefore the $\Delta C_A$ during 1000–1850 was largely driven by carbon-climate feedback without the complication of concurrent increase in $C_E$ as the industrial period. We ignored the possible influences of CO₂ emissions from the early land-use²⁹ on global surface temperature during 1000–1850. As the atmospheric CO₂ was relatively stable (280 ± 8 ppm) during the 1000–1850 period, i.e., $\Delta C_A$ was only ~3% of atmospheric CO₂ concentration (Fig. 1a), we assumed
<figure><img src="image_6.png" /><figcaption>**Fig. 2** Linear relationship between $\frac{1}{\alpha}$ and $\frac{1}{TCRE}$ across timescales. **a** Amplitude spectrum from Fourier analysis for annual atmospheric CO$_2$ ($\Delta C_A$) and cumulative anthropogenic CO$_2$ emissions ($\Delta C_E$) for 1880–2017. **b** Same as (a) but for annual global-mean temperature from four observational datasets (HadCRUT4, GISTEMP, Berkeley Earth, and NOAA GlobalTemp) for 1880–2017. **c** Estimates of $\frac{1}{\alpha}$ ($=\Delta C_A/\Delta T_A$) across timescales against $\frac{1}{TCRE}$ ($=\Delta C_E/\Delta T_E$) for four observational global-mean temperature datasets with atmospheric CO$_2$ records (or anthropogenic emissions) for 1880–2017. The solid line with the shaded area between the dashed lines is the linear regression of all datasets with a slope (p) of 2.52 ± 0.15 ppm ppm$^{-1}$, which indicates that the carbon-concentration feedback parameter $\beta$ is 3.22 ± 0.32 GtC ppm$^{-1}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 3** Ensemble estimates of climate-carbon cycle feedback parameters from observations and models. **a** Histogram of the $\beta$ estimates (error bars for mean ± 1σ) derived from four instrumental temperature and CO$_2$ datasets for 1880–2017 (in gray), $\beta^{BGC}$ derived from 11 C$^4$MIP models$^7$ for 1880–2017 (in dodger blue) and for 2018–2100 (in navy blue), and $\beta^{BGC}$ derived from nine CMIP5 models$^{11}$ for the 140-year 1pctCO$_2$ experiment (in deep pink). **b** Boxplot of the $\gamma^*$ derived from four instrumental temperature and CO$_2$ for 1880–2017 (in gray), and $\gamma_{100yr}$ estimates derived from an ensemble of 15 members (3 ice-core CO$_2$ records × 5 reconstructed temperature, in black) and from the EnOBS for 1000–1850 (in green), and $\gamma^{COU-BGC}$ derived from 11 C$^4$MIP models for 1880–2017 (in dodger blue) and for 2018–2100 (in navy blue), and $\gamma^{COU-BGC}$ derived from nine CMIP5 models (in deep pink). In **b**, the large ensemble estimates of the $\gamma_{100yr}$ from the EnOBS also provides probability distributions (right panel) of the cooler period of Little Ice Age (LIA, 1400–1700) and the warmer period of Medieval Climate Anomaly (MCA, 1000–1300), compared to those of the full preindustrial last millennium (PILM, 1000–1850). **c** Same as (a), but for the feedback gain factor (g) derived from the $\beta$ and AF estimates.</figcaption></figure>

that its effects on $\beta$-feedback and $\gamma$-feedback were small (e.g., $f(\beta, \gamma)\Delta C_A\Delta T_A \approx 0$). For the 1000–1850 period with $\Delta C_E \approx 0$, Eq. (2) can be further simplified as

$$ \gamma_k = -m(1 + \beta)\eta_k - f(\beta, \gamma)\Delta C_A \approx -m(1 + \beta)\eta_k \quad (3) $$

where $m = 2.12$ GtC ppm$^{-1}$ is a factor for converting atmospheric CO$_2$ in ppm to GtC. As the $\beta$ is found to be nearly constant during the industrial period, we assumed that the constant value is applicable to the 1000–1850 period. This is reasonable, as the error sensitivity analysis showed that the uncertainty in the estimated $\gamma$ for 1000–1850 induced from the uncertainty in $\beta$ is still smaller than those from the difference in the estimated $\gamma$ from ice-core CO$_2$ records and reconstructed temperatures (Supplementary Fig. 5). Thus, using 15 combinations of 3 ice-core CO$_2$ records$^{28,30,31}$ and 5 temperature reconstructions$^{16,32–35}$ as shown in Fig. 1 and Eq. (3) with $\beta$ of

3.22 ± 0.32 GtC ppm$^{-1}$, we estimated that the $\gamma$ at a 100-year timescale ($\gamma_{100yr}$) for 1000–1850 was −122.8 ± 60.2 GtC K$^{-1}$ (Fig. 3b and Supplementary Table 2). Furthermore, $\gamma_k$ increases with timescale, varying from −33 ± 14 GtC K$^{-1}$ on timescales of 10–70 years to −110 ± 40 GtC K$^{-1}$ over timescales of 200–800 years (Supplementary Fig. 6).

Using another set of >1500 combinations of 521 reconstructed temperature records from 1000 to 1850$^{16}$ and 3 ice-core CO$_2$ records (see “Methods”, Supplementary Fig. 4), we calculated much larger ensemble (>1500) estimates (EnOBS) of $\eta_k$ and $\gamma_k$ using the Fourier analysis-based approach. Our results confirmed that $\eta_k$ systematically increased with timescale, varying from 6 ± 6 ppm K$^{-1}$ on timescales of 10–70 years to 20 ± 8 ppm K$^{-1}$ over timescales of 200–800 years due to the increased $\gamma_k$ on magnitude with timescale (Supplementary Fig. 6). The range of the EnOBS-based $\eta_{100yr}$ for 1000–1850 was 7–23 ppm K$^{-1}$ with a median of
9 ppm K⁻¹, consistent with the estimate of η (1.7–21.4 ppm K⁻¹) by Frank et al.¹⁶. The resulting γ₁₀₀γr of EnOBS for 1000–1850 was −82 ± 42 GtC K⁻¹, varying from a minimum of −260 GtC K⁻¹ to a maximum of −15 GtC K⁻¹ (Fig. 3b). The EnOBS-based γ₁₀₀γr over timescales of 200–800 years was −80 ± 50 GtC K⁻¹ (Supplementary Fig. 6), which is 30% smaller than the estimate of γ₁₀₀γr from 3 ice-core CO₂ records × 5 temperature reconstructions. These results suggest that the timescale or temporal dependance of η over the 1000–1850¹⁶ is largely driven by the positive feedback of terrestrial and oceanic carbon pools to climate (i.e., the γ feedback), implying that on longer timescales, warming of the climate would cause more release of CO₂ into the atmosphere and in return, amplify warming.

The estimate of γ also depends on the mean climate state (quasi-equilibrium or transient). Here we compared estimates of γ between the warmer period of the Medieval Climate Anomaly (1000–1300) with a mean temperature anomaly of −0.18 ± 0.09 K and the cooler period of the Little Ice Age (1400–1700) with a mean temperature anomaly of −0.41 ± 0.03 K, and both anomalies were calculated relative to the mean global surface temperature from 1961 to 1990. We further analyzed probability distributions of the EnOBS-based γ₁₀₀γr for three different periods: 1000–1850, 1000–1300, and 1400–1700. The estimates of γ₁₀₀γr for the cooler period of 1400–1700 had a comparable probability distribution with the 1000–1850 period (Fig. 3b). However, the γ₁₀₀γr for the warmer period of 1000–1300 showed a much narrower distribution with a less negative mean (−38 ± 17 GtC K⁻¹, or 50% smaller) than that for 1400–1700 (−72 ± 56 GtC K⁻¹), suggesting that γ₁₀₀γr feedback is sensitive to the state of the mean climate. The higher γ₁₀₀γr for 1400–1700 resulted from a colder climate and some drastic fluctuations in CO₂ (e.g., larger ΔCₐ leads to more nonlinear feedback contribution from f(β, y)ΔCₐ), especially the strong dip at ~1600 AD, while the lower γ₁₀₀γr for 1000–1300 is associated with smaller variations in CO₂ and a warmer climate (Fig. 1). The CO₂ drop in the Little Ice Age (1400–1700) may have been driven not only by natural disturbances (e.g., volcanic eruptions) but also probably by human-induced land-use effects³⁶, which could have led to a more negative estimate of γ₁₀₀γr for 1400–1700.

**Comparisons of observation-based and model-based estimates of feedback parameters.** Previous studies have noted the timescale dependence of η (i.e., ΔCₐ/ΔTₐ) and the relationship between η and the climate-carbon cycle feedback parameters (β and γ)¹⁰,¹⁶,³⁷. Here we further quantified the different timescale dependence of β and γ. Using an observation-based constraint¹⁸, Cox and Jones estimated that β was between 3 and 5 GtC ppm⁻¹ and γ fell between −250 and −50 GtC K⁻¹. For the β-feedback of land, recent studies estimated that the CO₂ fertilization effect on global plant biomass carbon only was 25 ± 4 GtC year⁻¹ for a 100 ppm ΔCO₂ over 1980–2010³⁸, and for the global terrestrial C sink was 3.5 ± 1.9 GtC year⁻¹ per 100 ppm ΔCO₂ over 1959 to 2010³⁹, which is equivalent to a land β-feedback of 1.75 ± 0.95 GtC ppm⁻¹ over the same period. Our findings indicate that observation-based (land + ocean) β is approximately constant across different timescales at 3.22 ± 0.32 GtC ppm⁻¹ for 1880–2017. While magnitudes for γ increased with timescale from decadal to multi-centennial, and our estimate of γ₁₀₀γr had a narrower uncertainty than previous estimates, with an average of −122.82 ± 60.16 GtC K⁻¹ from the 3 ice-core CO₂ records × 5 temperature reconstructions and of −81.91 ± 41.90 GtC K⁻¹ from EnOBS for 1000–1850. Uncertainties related to the estimated γ₁₀₀γr from the two datasets (15 members versus >1500 members of EnOBS) did not overlap, which may suggest the uncertainty in γ₁₀₀γr was from −180 to −40 GtC K⁻¹.

We then compared our observation-based estimates of β and γ using Fourier analysis-based approach with those from Earth system models. Estimates of β and γ from Earth system models were calculated from three sets of model simulations: the biogeochemically-coupled (BGC), radiatively coupled (RAD), and fully coupled (COU) simulations based on the FEA approach⁷,¹⁴ (see “Methods”). Eleven first-generation coupled climate-carbon cycle models (C⁴MIP) were driven with the prescribed CO₂ emissions from the historical period (1860–2005) and the future period (2006–2100) under the IPCC SRES A2 scenario without land-use change⁷. Only COU and BGC simulations were conducted for the C⁴MIP models. The subsequent analysis used nine models from phase 5 of the Coupled Model Inter-comparison Project (CMIP5) used 1% year⁻¹ increasing CO₂ for 140 years, no land-use change (or the 1pctCO₂ experiment)¹¹. β and γ were estimated from three simulations (COU, BGC, and RAD) by each participating model. Because of the nonlinear feedback, the changes in size of carbon pools (land, ocean, and atmosphere) in the COU simulation were not equal to the sum of the simulated changes of those pool sizes in the RAD and BGC simulations¹⁴,⁴⁰,⁴¹. According to “Methods”, we calculated the direct β-feedback from the COU-BGC simulations (βGC ≈ ΔC<sup>BGC</sup><sub>B</sub>/ΔC<sup>BGC</sup><sub>A</sub>) and the direct γ-feedback from the COU-RAD simulations (γRAD = ΔC<sup>RAD</sup><sub>B</sub>/ΔT<sup>RAD</sup><sub>A</sub>) and the total γ-feedback (direct plus indirect) from the COU-BGC simulations (γCOU-BGC ≈ (ΔC<sup>COU</sup><sub>B</sub> − ΔC<sup>BGC</sup><sub>B</sub>)/ΔT<sup>COU</sup><sub>A</sub>), in theory this is the γ*) using the FEA approach, for the observation-overlapped period of 1880–2017 and the future emission scenario of 1880–2100 for the C⁴MIP models, and the 1pctCO₂ 140-year period for the CMIP5 models, respectively. We also estimated the nonlinear feedback term from the difference between COU simulations and the BGC and RAD simulations (f(β, y) ≈ [ΔC<sup>COU</sup><sub>B</sub> − (ΔC<sup>BGC</sup><sub>B</sub> + ΔC<sup>RAD</sup><sub>B</sub>)]/ΔC<sup>COU</sup><sub>A</sub>ΔT<sup>COU</sup><sub>A</sub>) and its contribution to γ-feedback (f(β, y)ΔC<sup>COU</sup><sub>A</sub>) for the CMIP5 models (“Methods”).

The estimated β<sup>BGC</sup> from 11 C⁴MIP models for 1880–2017 using the FEA approach was 3.07 ± 0.68 GtC ppm⁻¹ which is close to the observation-based β (3.22 ± 0.32 GtC ppm⁻¹) for the same period using historical CO₂ emissions as forcing. Previous studies demonstrated that the carbon-concentration feedback was strongly dependent on the growth rate of atmospheric CO₂ and hence on emission scenarios¹⁰. When extending the calculation to the future high emission scenario of the IPCC SRES A2 (close to the RCP8.5 pathway, the mean CO₂ growth rate was about 0.72% year⁻¹), we found that β<sup>BGC</sup> from C⁴MIP models was largely reduced by 19.5% to 2.47 ± 0.60 GtC ppm⁻¹ for the period of 1880–2100 (Fig. 3a, Supplementary Table 3). As the 1pccO₂ experiments of CMIP5 models were configured under a higher emission scenario (1% year⁻¹ of the CO₂ growth rate), we found that the estimated β<sup>BGC</sup> from nine CMIP5 models using the FEA approach was 1.71 ± 0.44 GtC ppm⁻¹, which is 45% smaller than the observation-based or C⁴MIP-based β for the historical period 1880–2017 (Fig. 3a), indicating that carbon-concentration feedbacks became smaller under higher growth rates of CO₂ in the climate-carbon cycle system.

The γCOU-BGC is theoretically the γ* in this study, as the γCOU-BGC was calculated from all feedback effects in COU simulations minus the direct β-feedback in BGC simulations. We found that the estimated γCOU-BGC from 11 C⁴MIP models were on magnitude increased from −27.52 ± 11.93 GtC K⁻¹ for 1880–2017 to −52.18 ± 26.54 GtC K⁻¹ for 1880–2100, when the βBGC was decreased for 1880–2100 (Fig. 3a, b, Supplementary Table 3). The estimated γCOU-BGC from nine CMIP5 models was −70.14 ± 32.43 GtC K⁻¹, which theoretically came from direct feedback γRAD (−65.08 ± 30.74 GtC K⁻¹) and the nonlinear
feedback contribution $f(\beta, \gamma)\Delta C_A$ ($-9.6 \pm 10.03 \text{ GtC K}^{-1}$, ~15 ± 17% of the $\gamma$-feedback) under the high emission scenario of the 1pctCO$_2$ experiments (Fig. 3b, Supplementary Table 4). Without CO$_2$ impact, the estimated nonlinear feedback parameter $f(\beta, \gamma)$ was $-11.22 \pm 11.72 \times 10^{-3} \text{ GtC ppm}^{-1} \text{ K}^{-1}$ (Supplementary Table 4). On the other hand, as the observation-based $\gamma_k$ for 1000–1850 and the model-based $\gamma^{\text{COU-BGC}}$ refer to different periods and timescales which exaggerates their differences owing to the timescale dependency of $\gamma$-feedback, it should be noted that there is limited comparability between them.

To further assess whether $\beta$ varies with timescale, we estimated $\beta_k^{\text{BGC}}$ and $\gamma_k^{\text{RAD}}$ by applying the FEA approach to the simulations from the nine CMIP5 models from inter-annual to centennial timescales (see “Methods” and Supplementary Texts 1, 2). The results supported that $\beta$-feedback was approximately constant across different timescales while $\gamma$-feedback increased significantly with timescale for most models (Supplementary Fig. 7a, b). Similar timescale dependences were also found for $\beta_k^{\text{BGC}}$ and $\gamma_k^{\text{RAD}}$ by applying the FEA approach to the simulations from C$^4$MIP models (Supplementary Fig. 8).

### Estimates of feedback gain (g) for 1880–2017 and high CO$_2$ emissions scenarios.
For further comparisons, we estimated the observational and model-based airborne fraction (AF) of cumulative CO$_2$ emissions ($AF = \Delta C_A/\Delta C_E$) and climate-carbon cycle feedback gain factor (g) across timescales ($g = 1 - 1/[AF(1 + \beta)]$, see “Methods”). We estimated AF$_k$ over different timescales from amplitudes of atmospheric CO$_2$ concentration and cumulative CO$_2$ emissions using Fourier analysis (Fig. 2a), then calculated the average and standard deviation of AF from AF$_k$ over timescales (Supplementary Fig. 7). We showed that the observation-based AF for the period 1880–2017 was nearly constant across timescales with an average of $0.40 \pm 0.05$ (Supplementary Fig. 7c), suggesting that about 60% of CO$_2$ emission was taken up by land and ocean. Previous studies found the relationship between AF and the two feedback parameters: $AF = 1/(1 + \beta + \alpha\gamma^{*})^{10}$. As for the period 1880–2017, we found a relatively smaller contribution of $\alpha\gamma^{*}$ (~0.05 GtC GtC$^{-1}$) from $\alpha$ ($0.005 \text{ K GtC}^{-1}$) and $\gamma^{*}$ (~10.9 GtC K$^{-1}$) compared to $\beta$ (~1.52 GtC GtC$^{-1}$) from the observational estimates, i.e., $\beta \gg \alpha\gamma^{*}$, suggesting that the carbon-concentration feedback ($\beta$) dominated the relatively stable cumulative airborne fraction over the industrial period.

Further analysis showed that the observation-based feedback gain (g) was very small ($0.01 \pm 0.05$) for 1880–2017 (Fig. 3c). We showed that C$^4$MIP-based AF for 1880–2017 was also nearly constant ($0.45 \pm 0.06$) across different timescales (Supplementary Fig. 7c), slightly larger than the observed estimate. However, the C$^4$MIP-based g for the same period 1880–2017 was $0.09 \pm 0.04$, larger than the observational value by about an order of magnitude (Fig. 3c). As a result, the observation-based and C$^4$MIP-based feedback amplification G ($G = 1/(1-g)$, see “Methods”) are $1.01 \pm 0.05$ and $1.10 \pm 0.04$, respectively, suggesting the modeled amplification effect is about 9 ± 7% larger. Under high emission scenarios, the C$^4$MIP-based g increased to $0.15 \pm 0.08$ with increased AF ($0.56 \pm 0.09$) for 1880–2100, and the CMIP5-based g increased to $0.13 \pm 0.08$ for the 1pctCO$_2$ (Fig. 3c), which are much higher than the observation-based g.

## Discussion

This study expanded the traditional climate-carbon cycle feedback framework$^{7,10,14,19}$ by including a nonlinear feedback term (see Eq. (1)). Using CMIP5 modeling experiments under a scenario of high CO$_2$ growth rate (1% yr$^{-1}$), we found that the contribution of the nonlinear term to land and ocean $\beta$-feedback is relatively small ($3 \pm 3\%$), while its contribution to land and

ocean $\gamma$-feedback is $15 \pm 23\%$. These estimates are noticeably smaller than those simulated from previous modeling studies on land and ocean (20% for $\beta$ and 45% for $\gamma$), or on ocean only (6% for $\beta$ and 60% for $\gamma$)$^{10,14}$ (see also Supplementary Text 3 for detailed discussion). However, the estimated nonlinear contributions vary significantly across CMIP5 models (0.2–9.6% for $\beta$ and 0.8–45% for $\gamma$, see Supplementary Table 4), which suggests a large uncertainty in the modeled nonlinear feedback among the advanced Earth system models.

This study also stated the relationship among the four commonly used quantities ($\beta$, $\gamma$, $\alpha$, and TCRE) which is consistent to some similar relationships as stated by previous carbon-climate feedback studies$^{10,21}$. Using historical temperature and CO$_2$ records over 1880–2017, we estimated that the carbon-concentration feedback parameter $\beta$ is $3.22 \pm 0.32 \text{ GtC ppm}^{-1}$ and the climate-carbon cycle feedback gain factor $g$ is $0.01 \pm 0.05$. These estimates were nearly constant across inter-annual to decadal timescales. On the other hand, we also found that the $\gamma$-feedback parameter increased with timescale from $-33 \pm 14 \text{ GtC K}^{-1}$ on a decadal scale to $-122 \pm 60 \text{ GtC K}^{-1}$ on a centennial scale based on reconstructions over 1000–1850. Furthermore, the estimated climate amplification from carbon-climate feedback based on observations in this study is much smaller than the previous estimates by IPCC reports based on model simulations under high emission scenarios$^{12,13}$. Our results based on observations have significant implications for understanding the strength of the climate-carbon cycle feedback and the allowable CO$_2$ emissions to mitigate future climate change. For example, the allowable emissions based on Earth system models for a 68% probability of limiting warming to 1.5 °C above the preindustrial level by ~2050 as stipulated in the Paris Agreement$^{42}$, are ~115 GtC (or 420 Gt CO$_2$)$^{43}$, using the much smaller feedback gain of $0.01 \pm 0.05$ than the C$^4$MIP-based estimates during the same period, we estimate that the allowable emissions would be 9 ± 7% more, or 125 ± 8 GtC.

In this study, the observation-based $\beta$ and AF were found to be nearly constant (with <10% change), which together determined the nearly constant feedback gain (g). This can be explained by the linear system in response to exponential increase of forcing (LinExp) theory$^{44,45}$, in which the carbon-climate system over the industrial period can be approximated as a linear system of the carbon cycle forced by exponentially growing CO$_2$ emissions with $y = 0.27e^{0.018t}$ for the 1850–2017 (Supplementary Fig. 9a, b), then all ratios of responses to forcings are constant$^{44}$. Using the simulations over the historical period (1901–2010) from 14 terrestrial ecosystem models forced by observational climate, land-use change, and atmospheric CO$_2$$^{46}$, we showed that when excluding the effects of climate change and land use change, both annual global GPP and cumulative land carbon sink increased exponentially in response to exponentially growing atmospheric CO$_2$ over 1901–2010 (Supplementary Fig. 9c, d). Therefore, the ratio of the exponential increase in carbon uptake and the exponential increase in CO$_2$ sustained a nearly constant value of $\beta$ over the historical period.

The timescale dependency of the $\gamma$-feedback can be used to reconcile the diverging estimates of $\gamma$ parameters in the previous studies. The estimated $\gamma$ at centennial timescales for the pre-industrial period ($-180 \text{ to } -40 \text{ GtC K}^{-1}$) was found to be much more negative than the value of $\gamma^*$ at multi-decadal timescales ($-14 \text{ to } -7 \text{ GtC K}^{-1}$) for the industrial period, implying that the climate sensitivity of carbon cycle depends on the base climate state. If we allow the present transient climate system to reach an equilibrium state, the estimated climate sensitivity of the carbon cycle would be more negative. Furthermore, the $\gamma$ (or $\gamma^*$)-feedback also depends on CO$_2$ emissions scenario. By comparing observation-based $\gamma^*$ and model-based $\gamma^{\text{COU-BGC}}$ from C$^4$MIP
models for the industrial period (1880–2017) with the model-based $γ^{COU-BGC}$ estimate for the future period with higher CO₂ emissions, we showed that under the conditions of the higher CO₂ growth and warmer climate states, the $γ$-feedback and the nonlinear feedback contribution to the $γ^*$ became more negative, while the $β$-feedback decreased (Fig. 3, Supplementary Table 3). The contribution of the nonlinear carbon-climate feedback to total carbon-climate feedback was estimated to be 15 ± 17% for the 1pctCO₂ high emission scenario as used in the CMIP5 study. For the ocean, the nonlinear feedback may become greater due to the decreased downward carbon transport to the deep ocean owing to reduced overturning in a warmer climate¹⁴. The $γ$ encapsulates both the direct effect of warming on plant and soil respiration⁴⁷,⁴⁸ and sea water solubility ($γ$ gets more negative), but also the indirect effect of warming on biological productivity and phenology⁴⁹ and on stratification/reduction of the overturning¹⁴,⁴¹ ($γ$ gets less negative). How the $γ$-feedback will change in a warmer climate with high CO₂ emissions depends on the contributions by these competing processes, and the uncertainties in the earth system models used for the analysis⁴⁵⁰. Overall, the $γ$-feedback at equilibrium under a high emissions scenario would become much more negative than the $γ$-feedback for the industrial or preindustrial periods, and likely play the dominant role in the positive net climate-carbon cycle feedback.

Uncertainties in the estimates of $β$, $γ$, and corresponding $g$ arise from the errors in observations, model uncertainties, methodology, and study periods chosen. First, we did not consider separately the contributions of non-CO₂ greenhouse gases (e.g., CH₄, N₂O, O₃) from those of CO₂ on climate variation, which may have led to some biases in the estimates of $β$ and $γ$. Second, the overestimation of C⁴MIP-based $g$ compared to the observation-based $g$ for the same period 1880–2017, might be due to several reasons in model descriptions, e.g., poor descriptions of terrestrial ecosystem processes (carbon pools, carbon-water coupling, soil respiration, vegetation phenology, and ecosystem climate adaptation, etc.), omission of land-use change, and biases in model initial climate/carbon pool base-states, etc. Third, errors could arise from the two approaches (the Fourier analysis approach and the FEA approach) for estimating observation-based $β$ and $γ^*$ or model-based $β^{BGC}$ and $γ^{COU-BGC}$ and corresponding $g$. Comparative analysis showed that the ensemble means of C⁴MIP-based $β^{BGC}$ and $γ^{COU-BGC}$ using the FEA approach are consistent with those of C⁴MIP-based $β$ and $γ^*$ using the Fourier analysis approach for the 1880–2017 period, indicating a relatively small difference of C⁴MIP-based $g$ between using the two approaches (Supplementary Fig. 11).

In addition, although the $β$-feedback was found to be stable over the industrial period with the approximately exponential growth of CO₂, applying the same $β$ value to the preindustrial period may result in bias in estimated $γ$ for the preindustrial period. Sensitivity analysis found that an overestimation (or underestimation) of 50% in $β$ would result in an underestimation (or overestimation) of about 30% in $γ_k$ at timescales from 10 to 1000 years for the 1000–1850 period (Supplementary Fig. 5). This change of 30% (~25 GtC K⁻¹) in $γ_{100yr}$ is still smaller than the uncertainty of $γ_{100yr}$ (a standard deviation of 41.90 GtC K⁻¹) that was mainly caused by the large divergences in the three ice-core CO₂ records and reconstructed temperatures. Furthermore, this study ignored the contribution of early land-use change to climate change ($ζ_{100yr} ≈ 0$) during 1000–1850, therefore may have overestimated the $γ_{100yr}$, e.g., for the Little Ice Age (1400–1700). The observed atmospheric CO₂ drop in the 1400–1700 period occurred because of the cooling-induced increase in terrestrial C uptake²⁷, but it can also be partly explained by enhanced forest restoration from reduced land-use change as a result of aban-

donment of agricultural land from the collapse of native population in the Americas during 1500–1650⁵¹. Despite these uncertainties, we demonstrated the nature of the timescale dependency of $γ$ on climate-carbon cycle feedback and important implications for observational constraints on earth system models for projecting future climate changes.

## Methods

**Anthropogenic CO₂ emission data.** We calculated the global annual total anthropogenic CO₂ emission flux ($F_E$) as the sum of annual emission fluxes from fossil fuel combustion and industrial processes ($F_{FF}$) and from land use and land cover changes ($F_{LUC}$) by human activity covering the period 1850–2017, both of which were obtained from the Global Carbon Project's annual global carbon budget report¹. The $F_{FF}$ (uncertainty of ±5% for a ±1σ (68%) confidence interval) was estimated by ref. ²² and the $F_{LUC}$ (uncertainty of ±0.7 GtC year⁻¹ representing a ±1σ confidence interval) was averaged from emission estimates based on two bookkeeping models by ref. ²³ and ref. ²⁴. We then calculated the cumulative annual anthropogenic CO₂ emission: $C_E(t) = \int_0^t F_E dt = \int_0^t (F_{FF} + F_{LUC}) dt$. Units were converted from GtC to ppm by dividing $m$ (=2.12 GtC ppm⁻¹) for comparison with atmospheric CO₂ concentration⁵². The conversion factor $m$ was also used for unit conversions of $β$, $γ$, $α$, and $η$ parameters.

**Instrumental and ice-core reconstructed atmospheric CO₂ records.** Global annual atmospheric CO₂ concentration over 1850–2017 was reconstructed from a combination of ice-core CO₂ records and instrumental CO₂ measurements since 1950s⁵³,⁵⁴. Information about the inter-annual variability in the CO₂ data covering 1850–1940s from ice-core CO₂ would be lost, as ice cores smooth atmospheric CO₂ records by fin diffusion, which may lead to some biases in the timescale analysis.

The highly resolved CO₂ covering the preindustrial period of 1000–1850 was compiled from ice-core records in Antarctica including three datasets at Law Dome³⁰,⁵⁵, WAIS Divide²⁸, and Dronning Maud Land³¹, respectively. As the three records have different resolutions, CO₂ values were resampled to the nearest calendar year and then smoothed using 30-year splines (Fig. 1a). However, large variations remain among the three records of CO₂ from ice-core (see Fig. 1a and Supplementary Fig. 3a).

**Instrumental and reconstructed temperature datasets.** All available instrumental temperature datasets including HadCRUT4⁵⁶ and Berkeley Earth⁵⁷ over 1850–2017, and GISTEMP⁵⁸ and NOAA GlobalTemp⁵⁹ over 1880–2017 were used. Global annual temperature anomalies were calculated by area-weighted averaging over both land (2 m air temperature) and ocean (sea surface temperature) minus the mean surface air temperature of land and ocean for the period from 1961 to 1990.

Five reconstructions of Northern Hemispheric temperature were obtained from four groups which were abbreviated as PAGES2k³², Frank2010¹⁶, Mann2009³⁴, MannEIV³³, and Moberg2005³⁵ in this study (Fig. 1b). All temperature reconstructions were smoothed with 5-year splines, and adjusted to removed-means with respect to 1961–1990. Northern Hemispheric-mean temperature had been shown to be highly representative for global-mean temperature variations¹⁶.

The ensemble estimates of 521 calibrated temperature reconstructions were obtained from the Frank2010¹⁶, in which all amplitude and variability on >30-year timescales are well preserved for uncertainty analysis. The ensemble estimates were recalibrated by nine available datasets of Northern Hemispheric-mean temperature reconstructions using a reconstructing technique based on a state-space time series and Kalman filter algorithm⁶⁰. The ensemble means (Fig. 1b) and individual temperature reconstructions (Supplementary Fig. 4a) by Frank2010 are used for a large ensemble of >1500 estimates (EnOBS) of the $η$ and $γ$ over the preindustrial period. Our estimate based on Fourier analysis shows that uncertainty in the $η$ on centennial timescales (Supplementary Fig. 4b) was consistent with the estimate by Frank2010 which was based on a lag-regression method (searching highest correlations) of CO₂ and temperature reconstructions which were smoothed with spline using a range of cutoff time (50, 75,..., 200 years)¹⁶.

**Theoretical analysis for the climate-carbon cycle feedback.** Following previous analysis⁷,¹⁰, the increase from anthropogenic CO₂ emissions ($ΔC_E$) is the sum of changes in the three stores,

$$ \Delta C_E = \Delta C_A + \Delta C_L + \Delta C_O \quad (4) $$

where $ΔC_A$, $ΔC_L$, and $ΔC_O$ are the changes in carbon storage on atmosphere, land, and ocean, respectively, over a time period ($Δt$) evolved from the reference climate state. For the analysis for industrial period, we calculated $ΔC_E$, $ΔC_A$, $ΔC_L$, and $ΔC_O$ (in a unit of GtC) during $Δt$ since 1880. Note that change in atmospheric CO₂ concentration in a unit of ppm can be converted to GtC by a factor of $m$ (~2.12 GtC ppm⁻¹)⁵². Within the climate-carbon cycle feedback system, following the linear carbon cycle feedback framework by refs. ⁷,¹⁰,¹⁹ and the nonlinearity of carbon cycle feedback discovered by ref. ¹⁴, we consider the changes in land and
$$y^{\text{COU-BGC}} = \frac{\Delta C_B^{\text{COU}} - \Delta C_B^{\text{BGC}}}{\Delta T_A^{\text{COU}} - \Delta T_A^{\text{BGC}}} \approx \frac{\Delta C_B^{\text{COU}} - \Delta C_B^{\text{BGC}}}{\Delta T_A^{\text{COU}}} \quad (23)$$

and the COU-RAD approach:

$$\beta^{\text{COU-RAD}} = \frac{\Delta C_B^{\text{RAD}} \Delta T_A^{\text{COU}} - \Delta C_B^{\text{COU}} \Delta T_A^{\text{RAD}}}{\Delta C_A^{\text{COU}} \Delta T_A^{\text{RAD}}} \approx \frac{\Delta C_B^{\text{COU}} - \Delta C_B^{\text{RAD}}}{\Delta C_A^{\text{COU}}} \quad (24)$$

$$y^{\text{RAD}} = \frac{\Delta C_B^{\text{RAD}}}{\Delta T_A^{\text{RAD}}} \quad (25)$$

where $\Delta C_A^{\text{COU}} = \Delta C_A^{\text{BGC}}$ for CMIP5 or C4MIP simulations. To calculate the nonlinear feedback parameter $f(\beta, \gamma)$, we also defined:

$$f(\beta, \gamma) = \frac{\Delta C_B^{\text{COU}} - (\Delta C_B^{\text{BGC}} + \Delta C_B^{\text{RAD}})}{\Delta C_A^{\text{COU}}(\Delta T_A^{\text{COU}} - \Delta T_A^{\text{BGC}})} \approx \frac{\Delta C_B^{\text{COU}} - (\Delta C_B^{\text{BGC}} + \Delta C_B^{\text{RAD}})}{\Delta C_A^{\text{COU}} \Delta T_A^{\text{COU}}} \quad (26)$$

We set the reference time at 1880 and then calculated the $\beta^{\text{BGC}}$, $\gamma^{\text{RAD}}$, $\gamma^{\text{COU-BGC}}$, and $f(\beta, \gamma)$ for the observation-overlapped period of 1880–2017 and the future emission scenario of 1880–2100 for the C4MIP models, and the 1pctCO₂ 140-year period for the CMIP5 models, respectively.

### A box model for diagnosing the climate-carbon cycle feedback parameters.

The coupled climate-carbon cycle system can be simplified as the combination of variations in temperature ($T_A = T_0 + \Delta T_A$) and CO₂ ($C_A = C_0 + \Delta C_A$) over time intervals ($\Delta t$). The change in $T_A$ is assumed to increase logarithmically with CO₂¹⁷,⁶⁶,

$$\Delta T_A = \frac{s}{\ln(2)} * \ln\left(\frac{C_A}{C_0}\right) + \varepsilon \quad (27)$$

where $s$ is the impact of CO₂ on the temperature that is suggested to be 1.5–4.5 K by the IPCC reports⁶⁷, and $\varepsilon$ is the residual term from climate internal variability. From Eq. (6),

$$\Delta C_A = \frac{\Delta C_E - \gamma^* \Delta T_A}{1 + \beta} \quad (28)$$

We then have the three parameters box model with input from CO₂ emissions,

$$C_A = C_0 + \frac{\Delta C_E - \gamma^*(T_A - T_0)}{1 + \beta} \quad (29)$$

$$T_A = T_0 + \frac{s}{\ln(2)} \ln\left(\frac{C_A}{C_0}\right) + \varepsilon \quad (30)$$

where $C_0$ and $T_0$ are initial values at the first year, here referenced to 1850. To validate the estimated $\beta$ and $\gamma^*$ for the industrial period, we applied this box model to predict temperature and CO₂ over 1850–2017 using annual cumulated CO₂ emissions, by setting $\beta = 3.22 \text{ GtC ppm}^{-1} = 1.52 \text{ GtC GtC}^{-1}$, $\gamma^* = -10.9 \text{ GtC K}^{-1}$, $s = 3 \pm 1.5 \text{ K}$ and $\varepsilon$ being the detrended $T_A$ anomaly time series from HadCRUT4 (Supplementary Fig. 2). Result shows that the increasing trend in predicted annual CO₂ is very close to the observation ($R^2 = 0.99$, RMSE = 3.5 ppm), while predicted temperature has a larger trend since 1980s ($R^2 = 0.96$, RMSE = 0.17 K), compared to observed records (Supplementary Fig. 2).

**Uncertainty.** We quantified the uncertainty by a $\pm 1\sigma$ (standard deviation) that represents a 68% confidence interval. The uncertainty in $C_E$ (Fig. 1) resulted from the uncertainty in $F_E$ ($\sigma_{F_E}$) that was calculated by $\sigma_{F_E} = \sqrt{\sigma_{FF}^2 + \sigma_{FLUC}^2}$. The uncertainty in $F_{FF}$ ($\sigma_{FFF}$) was $\pm 5\%$ of $F_{FF}^{22}$. This is consistent with a more detailed recent analysis of uncertainty for $F_{FF}^{68}$. The uncertainty in $F_{FLUC}$ ($\sigma_{FFLUC}$) was $\pm 0.7 \text{ GtC year}^{-1}$ at each year¹. To estimate the $\beta$ and its uncertainty ($\sigma_\beta$), we first estimated the $\beta_i$ and $\sigma_{\beta,i}$ ($i = 1, \dots, 4$) for each combination of 1 CO₂ × 1 CO₂ emission × 4 temperature datasets (HadCRUT4, GISTEMP, Berkeley Earth, and NOAA GlobalTemp) that included uncertainty on $C_E$ and the regression. We then estimated the mean value of $\beta$ by $\bar{\beta} = \frac{1}{4} \sum \beta_i$ and the $\sigma_\beta$ by $\sigma_\beta = \frac{1}{4} \sqrt{\sum \sigma_\beta^2}$ (Supplementary Table 1). Similarly, we estimated uncertainty in $\gamma$ ($\sigma_\gamma$) for the industrial period at a given timescale by $\sigma_\gamma = \frac{1}{4} \sqrt{\sum \sigma_\gamma^2}$. As the uncertainty in $\gamma$ ($\sigma_\gamma$) for the preindustrial period arises from uncertainties in ice core CO₂ records and temperature reconstructions, we estimated the $\sigma_\gamma$ over 1000–1850 or 1400–1700, or 1000–1300 periods at a given timescale by calculating the standard deviation of the ensemble of $\gamma$ from 3 ice-core CO₂ records × 5 temperature reconstructions (or the EnOBS): $\sigma_\gamma = \sqrt{\frac{1}{n-1} \sum_{i=1}^{n-1} (\gamma_i - \bar{\gamma})^2}$ (Supplementary Table 2, Fig. 3c, d). The estimated uncertainties in $\beta$ and $\gamma$ for C⁴MIP or CMIP5 models were calculated as the standard deviations of their ensembles of $\beta$ and $\gamma$ for each model (Supplementary Tables 3, 4). As $G = \frac{1}{1-g} \approx 1 + g$, uncertainty in $G$ ($\sigma_G$) was $\sigma_g$, and uncertainty in

relative change of $G$ from C⁴MIP models compared to OBS ($\frac{G_{\text{MOD}} - G_{\text{OBS}}}{G_{\text{OBS}}}$) was $\frac{G_{\text{MOD}}}{G_{\text{OBS}}}$ * $\sqrt{\left(\frac{\sigma_{\text{MOD}}}{G_{\text{MOD}}}\right)^2 + \left(\frac{\sigma_{\text{OBS}}}{G_{\text{OBS}}}\right)^2}$.

### Data availability

Resulting data that support the findings in this study are available in the Supplementary Tables, and more related datasets can be downloaded at https://github.com/xuanzezhang/climate-carbon-cycle-feedback. Global Carbon Project report 2018 dataset (including anthropogenic CO₂ emission): https://www.icos-cp.eu/GCP/2018. Observed atmospheric CO₂ is available at http://www.esrl.noaa.gov/gmd/ccgg/trends/global.html. Ice-core CO₂ and temperature reconstruction datasets are retrieved from their published papers listed in References. The HadCRUT4 is available at https://www.metoffice.gov.uk/hadobs/hadcrut4/. The GISTEMP is retrieved from https://data.giss.nasa.gov/gistemp/. NOAA GlobalTemp is retrieved from https://climatedataguide.ucar.edu/climate-data/global-surface-temperature-data-mlost-noaa-merged-land-ocean-surface-temperature/. Berkeley Earth temperature is available at http://berkeleyearth.org/data/. The coupled/uncoupled simulations of C⁴MIP and CMIP5 models are available at https://www.c4mip.net.

### Code availability

The NCL processing codes are available via GitHub at https://github.com/xuanzezhang/climate-carbon-cycle-feedback (https://doi.org/10.5281/zenodo.4575812) or from the corresponding author upon request.

Received: 8 February 2020; Accepted: 10 March 2021; Published online: 19 May 2021

### References

1. Le Quéré, C. et al. Global carbon budget 2017. Earth Syst. Sci. Data 10, 405–448 (2018).

2. Shevliakova, E. et al. Historical warming reduced due to enhanced land carbon uptake. Proc. Natl Acad. Sci. USA 110, 16730–16735 (2013).

3. Holden, P. B. et al. Climate-carbon cycle uncertainties and the Paris Agreement. Nat. Clim. Change 8, 609–613 (2018).

4. Friedlingstein, P. et al. Uncertainties in CMIP5 climate projections due to carbon cycle feedbacks. J. Clim. 27, 511–526 (2014).

5. Friedlingstein, P. Carbon cycle feedbacks and future climate change. Philos. Trans. A Math. Phys. Eng. Sci. 373, 20140421 (2015).

6. Cox, P. M. et al. Acceleration of global warming due to carbon-cycle feedbacks in a coupled climate model. Nature 408, 184–187 (2000).

7. Friedlingstein, P. et al. Climate-carbon cycle feedback analysis: results from the C4MIP model intercomparison. J. Clim. 19, 3337–3353 (2006).

8. Dufresne, J. L. et al. On the magnitude of positive feedback between future climate change and the carbon cycle. Geophys. Res. Lett. 29, 43–41–43–44 (2002).

9. Bolin, B. et al. Carbon Cycle Modelling (SCOPE 16) (John Wiley, New York, 1981).

10. Gregory, J. M., Jones, C. D., Cadule, P. & Friedlingstein, P. Quantifying carbon cycle feedbacks. J. Clim. 22, 5232–5250 (2009).

11. Arora, V. K. et al. Carbon-concentration and carbon-climate feedbacks in CMIP5 Earth system models. J. Clim. 26, 5289–5314 (2013).

12. Ciais, P. et al. Carbon and other biogeochemical cycles. in Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (eds. Stocker, T. F. et al.) (Cambridge University Press, 2013).

13. Denman, K. L. et al. Couplings between changes in the climate system and biogeochemistry. in Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change (eds. Solomon, S. et al.) (Cambridge University Press, 2007).

14. Schwinger, J. et al. Nonlinearity of ocean carbon cycle feedbacks in CMIP5 Earth system models. J. Clim. 27, 3869–3888 (2014).

15. Cox, P. M. et al. Sensitivity of tropical carbon to climate change constrained by carbon dioxide variability. Nature 494, 341–344 (2013).

16. Frank, D. C. et al. Ensemble reconstruction constraints on the global carbon cycle sensitivity to climate. Nature 463, 527–530 (2010).

17. Scheffer, M., Browkin, V. & Cox, P. M. Positive feedback between global warming and atmospheric CO₂ concentration inferred from past climate change. Geophys. Res. Lett. 33, L10702 (2006).

18. Cox, P. & Jones, C. Climate change. Illuminating the modern dance of climate and CO₂. Science 321, 1642–1644 (2008).
19. Friedlingstein, P., Dufresne, J. L., Cox, P. M. & Rayner, P. How positive is the feedback between climate change and the carbon cycle? *Tellus B* 55, 692–700 (2003).

20. Matthews, H. D., Gillett, N. P., Stott, P. A. & Zickfeld, K. The proportionality of global warming to cumulative carbon emissions. *Nature* 459, 829–832 (2009).

21. Jones, C. D. & Friedlingstein, P. Quantifying process-level uncertainty contributions to TCRE and carbon budgets for meeting Paris Agreement climate targets. *Environ. Res. Lett.* 15, 074019 (2020).

22. Boden, T. A., Marland, G. & Andres, R. J. *Global, Regional, and National Fossil-Fuel CO₂ Emissions*. Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, Tenn., USA (2017).

23. Houghton, R. A. & Nassikas, A. A. Global and regional fluxes of carbon from land use and land cover change 1850-2015. *Glob. Biogeochem. Cycles* 31, 456–472 (2017).

24. Hansis, E., Davis, S. J. & Pongratz, J. Relevance of methodological choices for accounting of land use change carbon fluxes. *Glob. Biogeochem. Cycles* 29, 1230–1246 (2015).

25. Myhre, G., Highwood, E. J., Shine, K. P. & Stordal, F. New estimates of radiative forcing due to well mixed greenhouse gases. *Geophys. Res. Lett.* 25, 2715–2718 (1998).

26. He, J. et al. Transient climate sensitivity depends on base climate ocean circulation. *J. Clim.* 30, 1493–1504 (2017).

27. Rubino, M. et al. Low atmospheric CO₂ levels during the Little Ice Age due to cooling-induced terrestrial uptake. *Nat. Geosci.* 9, 691–694 (2016).

28. Bauska, T. K. et al. Links between atmospheric carbon dioxide, the land carbon reservoir and climate over the past millennium. *Nat. Geosci.* 8, 383–387 (2015).

29. Pongratz, J., Reick, C., Raddatz, T. & Claussen, M. A reconstruction of global agricultural areas and land cover for the last millennium. *Glob. Biogeochem. Cycles* 22, L02709 (2008).

30. Rubino, M. et al. A revised 1000 year atmospheric δ¹³C-CO₂ record from Law Dome and South Pole, Antarctica. *J. Geophys. Res. -Atmos.* 118, 8482–8499 (2013).

31. Siegenthaler, U. et al. Supporting evidence from the EPICA Dronning Maud Land ice core for atmospheric CO₂ changes during the past millennium. *Tellus B* 57, 51–57 (2005).

32. PAGES2k Consortium. Continental-scale temperature variability during the past two millennia. *Nat. Geosci.* 6, 339–346 (2013).

33. Mann, M. E. et al. Proxy-based reconstructions of hemispheric and global surface temperature variations over the past two millennia. *Proc. Natl Acad. Sci. USA* 105, 13252–13257 (2008).

34. Mann, M. E. et al. Global signatures and dynamical origins of the Little Ice Age and Medieval Climate Anomaly. *Science* 326, 1256–1260 (2009).

35. Moberg, A. et al. Highly variable Northern Hemisphere temperatures reconstructed from low- and high-resolution proxy data. *Nature* 433, 613–617 (2005).

36. Koch, A., Brierley, C., Maslin, M. M. & Lewis, S. L. Earth system impacts of the European arrival and Great Dying in the Americas after 1492. *Quat. Sci. Rev.* 207, 13–36 (2019).

37. Willeit, M. et al. Time-scale and state dependence of the carbon-cycle feedback to climate. *Clim. Dyn.* 42, 1699–1713 (2014).

38. Terrer, C. et al. Nitrogen and phosphorus constrain the CO₂ fertilization of global plant biomass. *Nat. Clim. Change* 9, 684–689 (2019).

39. Liu, Y. et al. Field-experiment constraints on the enhancement of the terrestrial carbon sink by CO₂ fertilization. *Nat. Geosci.* 12, 809–814 (2019).

40. Arora, V. K. et al. Carbon-concentration and carbon-climate feedbacks in CMIP6 models and their comparison to CMIP5 models. *Biogeosciences* 17, 4173–4222 (2020).

41. Zickfeld, K. et al. Nonlinearity of carbon cycle feedbacks. *J. Clim.* 24, 4255–4275 (2011).

42. Paris Climate Agreement (UNFCCC, 2015).

43. IPCC. in *Global Warming of 1.5 °C* (eds Masson-Delmotte, V. et al.) (World Meteorological Organization, 2018).

44. Raupach, M. R. The exponential eigenmodes of the carbon-climate system, and their implications for ratios of responses to forcings. *Earth Syst. Dyn.* 4, 31–49 (2013).

45. Gloor, M., Sarmiento, J. L. & Gruber, N. What can be learned about carbon cycle climate feedbacks from the CO₂ airborne fraction? *Atmos. Chem. Phys.* 10, 7739–7751 (2010).

46. Huntzinger, D. N. et al. Uncertainty in the response of terrestrial carbon sink to environmental drivers undermines carbon-climate feedback predictions. *Sci. Rep.* 7, 4765 (2017).

47. Cheng, L. et al. Warming enhances old organic carbon decomposition through altering functional microbial communities. *ISME J.* 11, 1825–1835 (2017).

48. Schuur, E. A. G. et al. Climate change and the permafrost carbon feedback. *Nature* 520, 171–179 (2015).

49. Piao, S. et al. Plant phenology and global climate change: current progresses and challenges. *Glob. Change Biol.* 25, 1922–1940 (2019).

50. Randerson, J. T. et al. Multicentury changes in ocean and land contributions to the climate-carbon feedback. *Glob. Biogeochem. Cycles* 29, 744–759 (2015).

51. Stocker, B. D., Yu, Z. C., Massa, C. & Joos, F. Holocene peatland and ice-core data constraints on the timing and magnitude of CO₂ emissions from past land use. *Proc. Natl Acad. Sci. USA* 114, 1492–1497 (2017).

52. Ballantyne, A. P. et al. Increase in observed net carbon dioxide uptake by land and oceans during the past 50 years. *Nature* 488, 70–72 (2012).

53. Keeling, C. D. et al. Atmospheric carbon dioxide variations at Mauna Loa Observatory, Hawaii. *Tellus B* 28, 538–551 (1976).

54. Dlugokencky, E. & Tans, P. Trends in atmospheric carbon dioxide. http://www.esrl.noaa.gov/gmd/ccgg/trends/global.html (2018).

55. MacFarling Meure, C. et al. Law Dome CO₂, CH₄ and N₂O ice core records extended to 2000 years BP. *Geophys. Res. Lett.* 33, https://doi.org/10.1029/2006GL026152 (2006).

56. Morice, C. P., Kennedy, J. J., Rayner, N. A. & Jones, P. D. Quantifying uncertainties in global and regional temperature change using an ensemble of observational estimates: the HadCRUT4 data set. *J. Geophys. Res. -Atmos.* 117, D08101 (2012).

57. Rohde, R. et al. A new estimate of the average Earth surface land temperature spanning 1753 to 2011. *Geoinfor. Geostat.: Overv.* 01, 1–7 (2013).

58. Hansen, J., Ruedy, R., Sato, M. & Lo, K. Global surface temperature change. *Rev. Geophys.* 48, RG4004 (2010).

59. Vose, R. S. et al. NOAA's merged land-ocean surface temperature. *Anal. B. Am. Meteorol. Soc.* 93, 1677–1685 (2012).

60. Lee, T. C. K., Zwiers, F. W. & Tsao, M. Evaluation of proxy-based millennial reconstruction methods. *Clim. Dyn.* 31, 263–281 (2007).

61. Woodwell, G. M. et al. Biotic feedbacks in the warming of the earth. *Clim. Change* 40, 495–518 (1998).

62. Torrence, C. & Compo, G. P. A practical guide to wavelet. *Anal. B. Am. Meteorol. Soc.* 79, 61–78 (1998).

63. Zhang, X. et al. Dominant regions and drivers of the variability of the global land carbon sink across timescales. *Glob. Change Biol.* 24, 3954–3968 (2018).

64. Athreya, K. B. & Lahiri, S. N. *Measure Theory and Probability Theory* (Springer, 2006).

65. Stoica, P. & Moses, R. *Spectral Analysis of Signals* (Prentice Hall, 2005).

66. Budyko, M. I. *The Earth's Climate, Past and Future*, 1st edn (Elsevier, 1982).

67. IPCC. Climate change 2013: the physical science basis. *Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change*. Cambridge, United Kingdom and New York, NY, USA (2013).

68. Andres, R. J., Boden, T. A. & Higdon, D. A new evaluation of the uncertainty associated with CDIAC estimates of fossil fuel carbon dioxide emission. *Tellus B* 66, 23616 (2014).

## Acknowledgements

This study was funded by the National Key Research and Development Program of China (2017YFA0604603), the CAS Pioneer Talent Program, and the National Key Research and Development Program of China (2016YFA0602501). X.Z.Z. was also sponsored by the National Natural Science Foundation of China (42001019) and the Shanghai Sailing Program (19YF1413100). Y.P.W. was supported by the National Environmental Science Program (climate change and earth system science). P.C. acknowledges support from the European Research Council Synergy grant ERC-2013-SyG-610028 IMBALANCE-P. We gratefully thank Dr. Pierre Friedlingstein for providing constructive comments. We acknowledge the following data providers and model developers: the Global Carbon Project, the University of East Anglia Climatic Research Unit (CRU), the NASA Goddard Institute for Space Studies, the NOAA, and Berkeley Earth. We also gratefully thank Drs Chris Jones and Vivek Arora for providing the outputs of 11 C⁴MIP and nine CMIP5 models.

## Author contributions

X.Z.Z. and Y.P.W. designed the study. X.Z.Z. performed the analysis. X.Z.Z., Y.P.W., Y.Q.Z., and P.C. drafted the paper. P.J.R. proposed the conception of box model. K.H., Y.Q.L., S.L.P., Z.L.W., J.Y.X., W.Z., X.G.Z., and J.T. contributed to the interpretation of the results and to the text.

## Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-021-22392-w.

**Correspondence** and requests for materials should be addressed to X.Z., Y.-P.W. or Y.Z.

**Peer review information** *Nature Communications* thanks the anonymous reviewers for their contributions to the peer review of this work. Peer review reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021