# Field-experiment constraints on the enhancement of the terrestrial carbon sink by CO2 fertilization

Yongwen Liu, Shilong Piao, Thomas Gasser, Philippe Ciais, Hui Yang, Han

Wang, Trevor Keenan, Mengtian Huang, Shiqiang Wan, Jian Song, et al.

► To cite this version:

Yongwen Liu, Shilong Piao, Thomas Gasser, Philippe Ciais, Hui Yang, et al. Field-experiment constraints on the enhancement of the terrestrial carbon sink by CO2 fertilization. Nature Geoscience, 2019, 12 (10), pp.809-814. 10.1038/s41561-019-0436-1. hal-02895151

HAL Id: hal-02895151

https://hal.science/hal-02895151v1

Submitted on 2 Oct 2025

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
**Title**

Field-experiment constraints on the enhancement of the terrestrial carbon sink by CO2 fertilization

**Permalink**

https://escholarship.org/uc/item/1rx3n2tk

**Journal**

Nature Geoscience, 12(10)

**ISSN**

1752-0894

**Authors**

Liu, Yongwen
Piao, Shilong
Gasser, Thomas
et al.

**Publication Date**

2019-10-01

**DOI**

10.1038/s41561-019-0436-1

Peer reviewed
Field-experiment constraints on the enhancement of the
terrestrial carbon sink by CO₂ fertilization

Yongwen Liu¹,², Shilong Piao¹,²,³*, Thomas Gasser⁴, Philippe Ciais⁵, Hui Yang³, Han Wang⁶, Trevor F. Keenan⁷,⁸, Mengtian Huang³, Shiqiang Wan⁹, Jian Song⁹, Kai Wang³, Ivan A. Janssens¹⁰, Josep Peñuelas¹¹,¹², Chris Huntingford¹³, Xuhui Wang³, Muhammad Altaf Arain¹⁴, Yuanyuan Fang¹⁵, Joshua B. Fisher¹⁶, Maoyi Huang¹⁷, Deborah N. Huntzinger¹⁸, Akihiko Ito¹⁹, Atul K. Jain²⁰, Jiafu Mao²¹, Anna M. Michalak¹⁵, Changhui Peng²²,²³, Benjamin Poulter²⁴, Christopher Schwalm²⁵,²⁶, Xiaoying Shi²¹, Hanqin Tian²⁷, Yaxing Wei²¹, Ning Zeng²⁸, Qiuan Zhu²⁹ and Tao Wang¹,²

¹ Key Laboratory of Alpine Ecology, Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing, China. ² CAS Center for Excellence in Tibetan Plateau Earth Sciences, Chinese Academy of Sciences, Beijing, China. ³ Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University, Beijing, China. ⁴International Institute for Applied Systems Analysis, Laxenburg, Austria. ⁵ Laboratoire des Sciences du Climat et de l'Environnement, CEA-CNRS-UVSQ, Gif-sur-Yvette, France. ⁶Department of Earth System Science, Tsinghua University, Beijing, China. ⁷ Climate and Ecosystem Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA. ⁸Department of Environmental Science, Policy and Management, UC Berkeley, Berkeley, CA, USA. ⁹College of Life Sciences, Hebei University, Baoding, Hebei, China. ¹⁰Department of Biology, University of Antwerp, Wilrijk, Belgium. ¹¹CREAF, Cerdanyola del Valles, Barcelona, Catalonia, Spain. ¹²CSIC, Global Ecology Unit CREAF- CSIC-UAB, Bellaterra, Barcelona, Catalonia, Spain. ¹³Centre for Ecology and Hydrology, Wallingford, UK. ¹⁴School of Geography and Earth Sciences and McMaster Centre for Climate Change, McMaster University, Hamilton, Ontario, Canada. ¹⁵Department of Global Ecology, Carnegie Institution for Science, Stanford, CA, USA. ¹⁶Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA. ¹⁷Atmospheric Sciences and Global Change Division, Pacific Northwest National Laboratory, Richland, WA, USA. ¹⁸School of Earth Sciences and Environmental Sustainability, Northern Arizona University, Flagstaff, AZ, USA. ¹⁹National Institute for Environmental Studies, Tsukuba, Japan. ²⁰Department of Atmospheric Sciences, University of Illinois at Urbana-Champaign, Urbana, IL, USA. ²¹Environmental Sciences Division and Climate Change Science Institute, Oak Ridge National Laboratory, Oak Ridge, TN, USA. ²²Department of Biology Sciences, Institute of Environment Sciences, University of Quebec at Montreal, Montreal, Quebec, Canada.

²³Laboratory for Ecological Forecasting and Global Change, College of Forestry, Northwest A&F University, Yangling, Shaanxi, China. ²⁴NASA GSFC, Biospheric Sciences Lab., Greenbelt, MD, USA. ²⁵Woods Hole Research Center, Falmouth, MA, USA. ²⁶Center for Ecosystem Science and Society, Northern Arizona University, Flagstaff, AZ, USA. ²⁷International Center for Climate and Global Change Research and School of Forestry and Wildlife Sciences, Auburn University, Auburn, AL, USA. ²⁸Department of Atmospheric
and Oceanic Science, University of Maryland, College Park, MD, USA. ²⁹State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau, Northwest A&F University, Yangling, China.

*e-mail: slpiao@pku.edu.cn

## Abstract

Clarifying how increased atmospheric CO₂ concentration (eCO₂) contributes to accelerated land carbon sequestration remains important since this process is the largest negative feedback in the coupled carbon-climate system. Here, we constrain the sensitivity of the terrestrial carbon sink to eCO₂ over the temperate Northern Hemisphere for the past five decades, using 12 terrestrial ecosystem models and data from seven CO₂ enrichment experiments. This constraint uses the heuristic finding that the northern temperate carbon sink sensitivity to eCO₂ is linearly related to the site-scale sensitivity across the models. The emerging data-constrained eCO₂ sensitivity is 0.64 ± 0.28 PgC yr⁻¹ per hundred ppm of eCO₂. Extrapolating worldwide, this northern temperate sensitivity projects the global terrestrial carbon sink to increase by 3.5 ± 1.9 PgC yr⁻¹ for an increase in CO₂ of 100 ppm. This value suggests that CO₂ fertilization alone explains most of the observed increase in global land carbon sink since the 1960s. More CO₂ enrichment experiments, particularly in boreal, arctic and tropical ecosystems, are required to explain further the responsible processes.

## Introduction

Human activities have profoundly altered the global carbon cycle¹. About 600 PgC was emitted as CO₂ to the atmosphere during 1750-2015 due to the burning of fossil fuel, cement production and land-use change². It is estimated that terrestrial ecosystems responded to this perturbation by absorbing about 32% of the cumulative anthropogenic emissions² and therefore played a key role in mitigating climate change³,⁴,⁵. The terrestrial CO₂ sink (excluding emissions from land-use change) estimate doubled between the 1960s (period 1960-1969) and the current decade (2007-2016)² (Fig. 1a). Terrestrial ecosystem model simulations have indicated that increasing CO₂ (eCO₂) is the main driver of the increase in the terrestrial carbon sink⁵,⁶,⁷,⁸. However, the magnitude of this eCO₂ fertilization differs strongly between models (Fig. 1b). Reducing this uncertainty is critical for refining understanding of the role of land in the future evolution of the coupled terrestrial carbon cycle/climate system, one of the Grand Challenges acknowledged by the World Climate Research Programme⁹.
<figure><img src="image_2.png" /><figcaption>**Fig. 1** | Effect of eCO₂ on the change in the global residual terrestrial sink during the past five decades. a, Changes in the residual terrestrial sink (excluding emissions from land-use changes) from global budget² and the effect of eCO₂ on the residual terrestrial sink changes from 12 MsTMIP models (shown as multi-model mean ± s.d., shaded green area) relative to 1959. The dashed lines indicate unconstrained linear regressions. The P-values indicate the statistical significance of the linear regressions. b, The eCO₂-induced change in global residual terrestrial sink during the 1960s (1960-1969) and the 2000s (2000-2009) in the MsTMIP models. The grey rectangle spans the first to the third quartiles. The red and green segments inside the rectangle indicate the mean and median, respectively. The whiskers above and below the rectangle indicate the maximum and minimum, respectively. The triangles indicate the models with dynamic nitrogen cycles.</figcaption></figure>

Free-air CO₂-enrichment (FACE) experiments have been conducted at the ecosystem scale since the 1980s to measure the response of carbon fluxes and stocks to eCO₂ (ref. 10). These comprehensive experiments and data sets provide insights into local nutrient limitations and physiological mechanisms controlling the effect of eCO₂ on net primary productivity (NPP) and ecosystem carbon storage¹⁰,¹¹. However, the implications of their findings have not been scaled up to estimate large-scale effects of eCO₂ on carbon sinks. Overall, meta-analyses of FACE experiment results find that eCO₂ generally stimulates plant growth and increases ecosystem carbon storage¹²,¹³. The FACE experiments have been mainly in temperate ecosystems and generally lasted for only a few years¹⁰. The effects of eCO₂ differ with environmental conditions and local biotic factors¹⁴. They also change over the duration of experiments, in relation to the ability of plants to take up more nutrients¹⁵, for example by increased below-ground allocation to roots¹⁶ and mycorrhizae symbiosis¹⁷. New approaches are needed to link these effects, which are implicit in the FACE experiments, to large-scale model simulations of the CO₂ fertilization effect on the terrestrial carbon sink.

We quantify the sensitivity of terrestrial CO₂ sink to eCO₂ (parameter *B*) in the Northern Hemisphere temperate zone (23–50° N) by combining the results of terrestrial ecosystem models with FACE experiments using an emergent-constraint approach¹⁸ (Methods). Emergent-constraint approaches have been successfully applied to constrain the effects of climatic warming on tropical carbon storage¹⁸, the sensitivity of gross primary productivity (GPP) to eCO₂ (ref. 19), the projections of future global carbon sinks²⁰, the sensitivity of global rice yield to long-term climate change²¹ and the ratio of plant transpiration to land evapotranspiration²². We use the output of 12
terrestrial models from the Multi-scale Synthesis and Terrestrial Model Intercomparison Project (MsTMIP)²³ (Supplementary Table 1) and the observed site-scale sensitivity of the net ecosystem production (NEP) to eCO₂ (Methods) from seven FACE experiments. All seven FACE experiments have a duration longer than one growing season (Supplementary Tables 2 and 3). They additionally include both nitrogen (N)-limited and N-rich sites and both arbuscular-mycorrhizae- and ectomycorrhizae-dominated ecosystems. The emergent constraint used here is a heuristic relationship between modelled *B* at different spatial scales. That is, it links the sensitivity *B*<sub>NH</sub> for all northern temperate ecosystems and the modelled mean sensitivity from the experiment site locations, *B*<sub>site</sub> (Methods). This relationship can then be used with FACE data-based knowledge of *B*<sub>site</sub> to constrain the value of *B*<sub>NH</sub>.

## Comparison of observed and modelled sensitivities

The FACE experimental configuration is the introduction of a step increase of atmospheric CO₂. After this, NPP generally increases, but this stimulation often slows down after a few years, which many attribute to nutrient limitations (for example, Norby et al.²⁴). Following the initial step increase of NPP, heterotrophic respiration (HR) also increases. This lag between HR and NPP leads to a transient sink (NEP), which integrates the direct effects of eCO₂ on NPP and delayed effects on respiration. The latter is a function of carbon allocation²⁵ and the turnover times of excess carbon in biomass, litter and soil organic matter. By contrast, in the historical simulations of the global terrestrial ecosystem models, the forcing of atmospheric CO₂ on different terrestrial components can be isolated through separate factorial simulations. However, in those simulations, CO₂ does not increase abruptly but follows the historical atmospheric CO₂ trend (roughly corresponding to an increase of 1.4 ppm yr⁻¹ for 1959–2010). In response to such a gradual increase of CO₂, NPP simulated by models increases slowly, and HR follows this increase with a lag. The resulting historical response of NEP to eCO₂ in models is thus not directly comparable to the short-term response to a step change in the FACE experiments. To link the two responses (FACE-based short-term step increase of CO₂ versus land model-based gradual rise of CO₂), we develop substitute pulse-response models (Methods). These can infer the step-response of NEP to eCO₂ in each MsTMIP model to mimic the conditions of the FACE experiment sites (Supplementary Fig. 1). The parameters of substitute models are the NPP sensitivity to eCO₂ and the turnover rates of excess carbon in total biomass and soil carbon pools. As an initial test, the substitute models generally emulate successfully the eCO₂-induced historical total biomass and soil carbon evolution from the original complex process-oriented MsTMIP models (generally R² > 0.95, P < 0.0001, Supplementary Figs. 2–13). Nevertheless, the substitute models do not work well to reproduce the biomass carbon pool evolution at grassland sites for 5 of 12 MsTMIP models (5 models at Duolun and 3 models at prairie heating and CO₂ enrichment (PHACE)), particularly when NPP monotonically increases but biomass carbon pool suddenly declines (Supplementary Figs.
2, 5, 7-9). This, therefore, may further result in uncertainty in the estimates
of both Northern Hemisphere temperate and global terrestrial B.

A second issue needs to be addressed to enable comparison of MsTMIP
estimates with FACE data. Atmospheric CO₂ was increased abruptly at the
FACE sites from recent amounts (~381 ppm) to about twice pre-industrial
CO₂ levels. That is a step increase (ΔCO₂) of about 200 ppm (Supplementary
Table 2). However, ΔCO₂ at FACE sites is significantly larger than ΔCO₂ in
models during the recent historical period (that is, about 60 ppm between
1959-1968 and 2001-2010). Because the response of photosynthesis to
eCO₂ is non-linear and convex²⁶, we extrapolate the observed FACE-based
NEP sensitivity to ΔCO₂ = 200 ppm to ΔCO₂ = 60 ppm on a background of 320
ppm as observed during the past 50 years. This is performed using (1) a
theoretical model of photosynthesis (Pmodel), which is extensively supported
by eddy-covariance measurements of CO₂ exchange and FACE
observations²⁷, and (2) a logarithmic function of photosynthetic response to
eCO₂ (refs. 3·28). These algorithms provide us with observation-based

estimates of $B_{\text{Site}}$ consistent with historical CO$_2$ levels, that is,
(Methods). Nevertheless, such scaling is based on theoretical response
curves of photosynthesis to rising CO$_2$ and thus subject to an unknown level
of uncertainty due to incomplete knowledge on acclimation and the shape of
the true response of ecosystem photosynthesis to eCO$_2$.

Northern Hemisphere temperate terrestrial *B* constrained by FACE
observations

The sensitivity of the Northern Hemisphere temperate ecosystems carbon

sink to yearly changing CO₂, $B_{\text{NH}}^{\text{Mod}}$ corresponding to the period 1959-2010, is found to be linearly correlated across the 12 terrestrial ecosystem models, with the modelled short-term response to a CO₂ step-increase (of 60 ppm) at each FACE location ($R^2$ = 0.86, *P* < 0.0001; Fig. 2a). The latter is diagnosed from the emulator of each model, as outlined above. Critically, this emergent constraint result indicates that a correspondence can be found between expected response to slowing changing CO₂ levels and a step change. That is, models with a large step-response to an abrupt increase in CO₂ at site-

scale generally also produce a large $B_{\text{NH}}^{\text{Mod}}$. The linear emergent relationship
across the models in Fig. 2a therefore provides a mechanism to constrain the
impact of eCO₂ on Northern Hemisphere temperate ecosystems' NEP using

the sensitivities inferred from FACE experiments, $B_{\text{Site}}^{\text{Hist}}$ = 75 ± 28 g C m⁻² yr⁻¹
[100 ppm]⁻¹ (mean ± s.d., indicated by the light grey area in Fig. 2a). This

$B_{\text{Site}}^{\text{Hist}}$

quantity is extrapolation of FACE observations to smaller jumps in CO₂
(ΔCO₂ = 60 ppm) added to a background value of 320 ppm on the basis of the Pmodel²⁷, again as outlined above.
<figure><img src="image_3.png" /></figure>

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

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

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

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

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

<figure><img src="image_9.png" /><figcaption>**Fig. 2 |** Emergent constraints on the sensitivity of the Northern Hemisphere temperate terrestrial carbon sink to eCO₂. **a.** Relationship between the Northern Hemisphere temperate zone and the site-scale responses of the terrestrial carbon sink to eCO₂ (*B*<sub>NH</sub> and *B*<sub>Site</sub>) for 1959–2010 across the MsTMIP models. *B*<sub>Site</sub> is based on MsTMIP model simulation to be appropriate for a jump from 320 ppm to 380 ppm and extrapolation of FACE experiments also to a jump from 320 ppm to 380 ppm. The light grey area shows the observed *B*<sub>Site</sub> in the FACE experiments (*B*<sub>Site</sub><sup>Hist</sup>, shown as mean ± s.d.). The light red area shows FACE data-constrained *B*<sub>NH</sub> (*B*<sub>NH</sub><sup>Constrained</sup>, shown as mean ± s.d.). **b.** Probability density function for *B*<sub>NH</sub>. The histogram indicates the probability density function of unconstrained *B*<sub>NH</sub> in the MsTMIP models (*B*<sub>NH</sub><sup>Mod</sup>), assuming they are Gaussian distributed. The red line indicates the probability density function of *B*<sub>NH</sub><sup>Constrained</sup>. **c,d.** Same as **a** and **b** but using the carbon-only models. **e,f.** Same as **a** and **b** but using only the C-N models. **g.** Relationship between the Northern Hemisphere temperate zone and global responses of the terrestrial carbon sink to eCO₂ (*B*<sub>NH</sub> and *B*<sub>Globe</sub>) for 1959–2010 across the MsTMIP models. In **a,c,e** and **g**, the solid lines indicate unconstrained linear regressions and the triangles indicate the models with dynamic nitrogen cycles.</figcaption></figure>
The relationship shown in Fig. 2a produces a probability density function of $B_{\text{NH}}^{\text{Mod}}$ constrained by the eCO$_2$ experiment-based probability density

function of $B_{\text{Site}}^{\text{Hist}}$ (Fig. 2b). The unconstrained $B_{\text{NH}}^{\text{Mod}}$ of $0.58 \pm 0.33 \text{ PgC yr}^{-1}$ [$100 \text{ ppm}]^{-1$ (mean $\pm$ s.d.) is found to be lower than the observation-constrained value of $0.64 \pm 0.28 \text{ PgC yr}^{-1}$ [$100 \text{ ppm}]^{-1$ (indicated by the red line in Fig. 2b and the light red area in Fig. 2a). Multiplying the constrained $B_{\text{NH}}$ by the historical increase in atmospheric CO$_2$ (57 ppm) provides an estimate of the contribution of historical CO$_2$ fertilization to the Northern Hemisphere terrestrial carbon sink change as $0.36 \pm 0.16 \text{ PgC yr}^{-1}$ between the 1960s (1960–1969) and the 2000s (2000–2009). Despite the small quantity of FACE data, our estimate of constrained $B_{\text{NH}}$ is robust when calculated using six, five and four eCO$_2$ experiment sites separately (Supplementary Fig. 14). The constrained $B_{\text{NH}}$ is also robust using either the Pmodel (Methods equation [5])$^{27}$ at FACE sites or the alternative logarithmic function for the non-linear response of photosynthesis to rising CO$_2$ (Methods equation [14])$^{3,28}$ with all other conditions unchanged (Supplementary Fig. 15 and Methods).

Six of the twelve MsTMIP models incorporated processes controlling carbon-nitrogen (C-N) interactions (Supplementary Table 1), so we constrained $B_{\text{NH}}$ values from carbon-only models and from C-N models separately. The Northern Hemisphere temperate ecosystems and site-level B are correlated for both the carbon-only models ($R^2 = 0.90, P = 0.004$; Fig. 2c) and the C-N models ($R^2 = 0.82, P = 0.013$; Fig. 2e). The constrained value of $B_{\text{NH}}$ is $0.67 \pm 0.24 \text{ PgC yr}^{-1}$ [$100 \text{ ppm}]^{-1}$ for the carbon-only models (Fig. 2c,d) and $0.71 \pm 0.39 \text{ PgC yr}^{-1}$ [$100 \text{ ppm}]^{-1}$ for the C-N models (Fig. 2e,f); thus, we find no notable difference of constrained sensitivity depending on whether N cycling is included in models. Nevertheless, the slope of the linear relationship between $B_{\text{NH}}$ and $B_{\text{Site}}$ across the C-N models (slope = 0.013; Fig. 2e) is higher than that across the carbon-only models (slope = 0.008; Fig. 2c). This indicates that the emergent relationship across models can change when most models consider a new process, such as here for the N cycle. The relationship may be refined further when, for example, models routinely model the phosphorus cycle$^{29}$. In addition, the temperate sites in our study mainly reflect forest ecosystems, with only two grassland sites and no site in shrublands or other ecosystems. Therefore, the relationship between Northern Hemisphere temperate ecosystems and site-level B across models, and the data-constrained $B_{\text{NH}}$ value, may adjust when new FACE sites become available and are included.

Compared to the constrained $B_{\text{NH}}$ from the 12 models ($0.64 \pm 0.28 \text{ PgC yr}^{-1}$ [$100 \text{ ppm}]^{-1}$), $B_{\text{NH}}$ is slightly overestimated by carbon-only models ($0.74 \pm 0.33 \text{ PgC yr}^{-1}$ [$100 \text{ ppm}]^{-1}$), but is generally underestimated by C-N models ($0.41 \pm 0.26 \text{ PgC yr}^{-1}$ [$100 \text{ ppm}]^{-1}$) (Supplementary Table 4). Specifically, for the C-N models, the observation-constrained $B_{\text{NH}}$ is comparable with three of
the six C-N models (ISAM, TEM6 and TRIPLEX-GHG) but is underestimated by
the others (CLM4, CLMVIC and DLEM) (Fig. 2a). This is consistent with the
site-scale comparison between the observed and the modelled NEP

sensitivity to eCO₂ (B<sup>Hist</sup><sub>Site</sub> and B<sup>Mod</sup><sub>Site</sub>, Fig. 2a, Supplementary Fig. 16). In FACE experiments, N limitation regulated the eCO₂-induced increase in NPP at the Oak Ridge forest²⁴ and at the moderately fertile Duke forest³⁰, and regulated the decade-long eCO₂-stimulation on plant biomass in a temperate grassland³¹. A large influence of mycorrhizal association across sites regulated the biomass increase in response to eCO₂ through mycorrhizal-N uptake¹⁷. The observed plant N uptake increase by eCO₂ is generally underestimated by C-N models at Duke-FACE and ORNL-FACE sites¹¹, but the eCO₂ effects on plant N uptake remain fully characterized in the other experiments used in this study. The evaluation of N limitation on B<sub>NH</sub> in the C-N models in this study must be considered with caution because we do not know whether our set of eCO₂ experiments is representative of the real-world extent of N limitation in northern temperate ecosystems due to the limited number of eCO₂ experiments. Nevertheless, there is the suggestion that some models may have an overly strong N response.

Global terrestrial *B* constrained by FACE observations

We can also expand our findings to the global scale. Across the terrestrial
ecosystem models, B<sub>NH</sub> was linearly correlated with the sensitivity of the
global terrestrial carbon sink to eCO<sub>2</sub>, B<sub>Globe</sub> (R<sup>2</sup> = 0.79, P = 0.0001; Fig. 2g).
This means that models with a large B<sub>NH</sub> also simulate a large B<sub>Globe</sub>. Using the
linear relationship between temperate B<sub>Site</sub> and B<sub>Globe</sub> across the terrestrial
ecosystem models, we estimated a constrained value for B<sub>Globe</sub> of 3.5 ± 1.9
PgC yr<sup>-1</sup> [100 ppm]<sup>-1</sup> (Supplementary Fig. 17a,b and Methods), similar to
carbon-only and C-N models (Supplementary Fig. 17c-f). The simple yet
powerful emerging relationship in the MsTMIP model ensemble between
temperate and global ecosystem responses to eCO<sub>2</sub> does have caveats. It
may be due to terrestrial ecosystem models having similar structures and
thus possibly common biases spatially coherent over the globe. For example,
they do not include specific response processes that could limit the NEP
response to eCO<sub>2</sub> in tropical ecosystems, such as phosphorus limitation<sup>29</sup>,
and do not distinguish between the responses of arbuscular-mycorrhizae-
and ectomycorrhizae-dominated ecosystems, which was recently found to be
a first-order explanation of increased biomass under eCO<sub>2</sub> (ref. 17). The
biases in modelling the eCO<sub>2</sub> effect across spatial scales (from site to globe)
cannot be verified by FACE experiments only, due to the lack of FACE-based
B data in boreal and tropical ecosystems. The planned extension of FACE
experiments in tropical (for example, AmazonFACE) and boreal (for example,
SwedFACE) ecosystems<sup>32</sup>, particularly designed with a gradual rise of CO<sub>2</sub>
(ref. 33), will improve the understanding of the mechanisms underlying the
terrestrial carbon sink responses to eCO<sub>2</sub> away from mid-latitudes. Such
additional FACE experiments would confirm or disconfirm our model-based extrapolation to estimate the global fertilization effect.

Combining the constrained $B_{\text{Globe}}$ from the 12 models with the historical increase in atmospheric CO$_2$ (57 ppm) estimates the contribution of CO$_2$ fertilization to the global terrestrial carbon sink change of $2.01 \pm 1.06 \text{ PgC yr}^{-1}$ between the 1960s and the 2000s (Supplementary Fig. 18). This constrained value is comparable with the multi-model-mean of the carbon-only models but is underestimated by three of six C-N models (Supplementary Fig. 18). Subtracting the estimated eCO$_2$ effect and multi-model-mean of N deposition effect ($0.23 \pm 0.14 \text{ PgC yr}^{-1}$) (Supplementary Fig. 19) from the observed increase in residual land sink between the 1960s and the 2000s ($1.2 \text{ PgC yr}^{-1}$)$^2$, the effect of long-term warming ($\sim 0.9 \,^\circ\text{C}$) on the residual land sink change between the 1960s and the 2000s is roughly estimated as $-1.04 \text{ PgC yr}^{-1}$. This indicates a temperature sensitivity of global terrestrial carbon sink of $\Gamma_{\text{Long}} = -1.16 \text{ PgC yr}^{-1} \text{ K}^{-1}$. This value falls within the range of the Global Carbon Budget data-constrained $\Gamma_{\text{Long}}$ ($-1.57 \pm 0.51 \text{ PgC yr}^{-1} \text{ K}^{-1}$) on the basis of the linear relationship between long-term and interannual temperature sensitivities of the land net carbon sink across MsTMIP models (Supplementary Fig. 20 and Methods). Hence, the findings from the 'top-down' closure of the global carbon cycle$^2$ are comparable with our FACE-based 'bottom-up' emergent constraint value for $B_{\text{Globe}}$ of $3.5 \pm 1.9 \text{ PgC yr}^{-1} [100 \text{ ppm}]^{-1}$.

In summary, our results provide a constrained estimate of the sensitivity of the terrestrial carbon sink to eCO$_2$ for a transient climate responding to anthropogenic forcing. This is inferred by combining pulse-response functions implicit in terrestrial carbon-cycle models with FACE experiment observations, which can then be scaled to the CO$_2$ increase observed during the past five decades. We find a robust emergent constraint between these model estimates, mimicking FACE conditions via our inferred pulse-response of models, and their estimates of evolving temperate ecosystems NEP response to rising atmospheric CO$_2$ concentrations due to fossil-fuel burning. This allows us to extrapolate the direct FACE measurements to estimate an Northern Hemisphere temperate ecosystems NEP fertilization response of $0.64 \pm 0.28 \text{ PgC yr}^{-1} [100 \text{ ppm}]^{-1}$ for the evolving conditions over recent decades. The consistency between model temperate regions and elsewhere allows us to make a unique data-based constraint on global vegetation direct physiological NEP response to rising CO$_2$ concentrations of $3.5 \pm 1.9 \text{ PgC yr}^{-1} [100 \text{ ppm}]^{-1}$. Further verification is required, potentially through longer-term FACE experiments. These FACE experiments would enable a more comprehensive investigation of slowly evolving processes such as vegetation and soil carbon turnover rates, including their relationship to nutrients that may change for long-term eCO$_2$ exposure. To clarify the mechanisms underlying elevated CO$_2$ effects on terrestrial NEP, continued collaboration between experimentalists and modellers, for example, the FACE Model-Data Synthesis project$^{32,34}$, remains necessary.
# Methods

Observed site-scale sensitivity of NEP to eCO₂ (B<sub>Site</sub>) from CO₂ enrichment experiments

We collected the data at the eCO₂ experiments for the effect of eCO₂ on NEP from the publications (Supplementary Table 2). Experimental sites with observational periods shorter than one growing season in each measurement year were excluded. A total of seven eCO₂ experimental sites were used in

our study (Supplementary Table 2). The B<sub>Site</sub> observed at each site (B<sup>Obs</sup><sub>Site</sub>) was calculated as:

$$B_{\text{Site}}^{\text{Obs}} = \frac{\Delta \text{NEP}}{\Delta \text{CO}_2} \qquad (1)$$

where ΔNEP was calculated as the difference in NEP between the eCO₂ and control treatments, and ΔCO₂ is the atmospheric CO₂ concentrations in the eCO₂ treatment minus the control ambient concentration during the experiment; ΔCO₂ differs between experiments. For the FACE experiments at Oak Ridge National Laboratory (ORNL-FACE), Duke-FACE, Rhinelander (Aspen-FACE) and Merritt Island, ΔNEP caused by eCO₂ was calculated as the eCO₂-induced change in ecosystem carbon storage divided by the experimental period.

MsTMIP simulations used to analyse eCO₂ impact on terrestrial carbon sink

We used net ecosystem CO₂ exchange (NEE, variable 'NEE' in the MsTMIP product) from the simulations of 12 terrestrial carbon-cycle models from the MsTMIP (Supplementary Table 1). The outputs of all models were downloaded from https://doi.org/10.3334/ORNLDAAC/1225 (ref. 35). Driver data of MsTMIP were downloaded from https://doi.org/10.3334/ORNLDAAC/1220 (ref. 36). The outputs and driver data cover all land surface areas excluding Antarctica with a 0.5° × 0.5° spatial resolution. The models simulated four sensitivities to test the influence of climate, land-use and land-cover change (LULCC), CO₂ concentrations and nitrogen (N) deposition on the terrestrial carbon cycle: (1) the models were forced by time-varying climate (scenario SG1), (2) the models were forced by time-varying climate and LULCC (scenario SG2), (3) the models were forced by time-varying climate, LULCC and CO₂ concentration (scenario SG3), and (4) the models with C-N interactions were forced by time-varying climate, LULCC, CO₂ concentration and N deposition (scenario BG1). The effect of rising CO₂ was isolated as the differences between SG3 and SG2. To analyse eCO₂ impact on terrestrial carbon sink, we used a total of 12 models (CLM4, CLM4VIC, DLEM, GTEC, ISAM, LPJ-wsl, ORCHIDEE-LSCE, SiBCASA, TEM6, TRIPLEX-GHG, VEGAS2.1 and VISIT) with both SG2 and SG3 simulations, complying with the following conditions: they
provided variables for NEE (variable 'NEE'), NPP (variable 'NPP'), total biomass carbon pool (variable 'TotLivBiom', except 'CarbPools' for TRIPLEX-GHG SG2 simulations) and soil carbon pool (variable 'TotSoilCarb', except 'CarbPools' for CLM4VIC and VISIT), and we could apply the two-box model (see below) for calibrating the outputs. Here, soil carbon pool includes soil carbon and litter carbon.

Calculation of the changes in the global residual terrestrial carbon sink relative to 1959 and its drivers

The residual terrestrial sink ($S_{LAND}$) in MsTMIP was calculated as the integrated effect of time-varying climate, CO₂ concentrations and N deposition on global terrestrial carbon sink (Fig. 1a). The effect of climate change on $S_{LAND}$ was given by the SG1 scenario. The effects of the other forcing factors on $S_{LAND}$ were calculated using differences between the simulations: differences between SG3 and SG2 for the effect of rising CO₂ (Fig. 1b and Supplementary Fig. 18) and differences between BG1 and SG3 for the effect of variations in N deposition (Supplementary Fig. 19). The $S_{LAND}$ from Global Carbon Budget (GCB) 2017 (refs. 2·37) was estimated by subtracting the growth rate in atmospheric CO₂ concentration ($G_{ATM}$) and ocean CO₂ sink ($S_{OCEAN}$) from global emissions from fossil fuels and industry ($E_{FF}$) and land-use change ($E_{LUC}$): $S_{LAND} = E_{FF} + E_{LUC} - (G_{ATM} + S_{OCEAN})$ (Fig. 1a).

Sensitivity of Northern Hemisphere temperate (23–50° N) terrestrial carbon sink to eCO₂ ($B_{NH}$) in MsTMIP

The sensitivity of the Northern Hemisphere temperate terrestrial carbon sink

to eCO₂ ($B_{NH}^{Mod}$) for 1959–2010 was calculated as $B_{NH}^{Mod} = \Delta NEP / \Delta CO₂$, where $\Delta NEP$ is the difference in the CO₂-caused Northern Hemisphere temperate terrestrial NEP (variable 'NEE' in the MsTMIP product, SG3-SG2) between 2001–2010 and 1959–1968, and $\Delta CO₂$ is the difference in the global atmospheric CO₂ concentration between 2001–2010 and 1959–1968.

Sensitivity of the modelled site-scale NEP to eCO₂ ($B_{Site}^{Mod}$) for the seven eCO₂ experimental sites in MsTMIP

We established a two-box model for each model and each site to obtain the effects of short-term pulses of eCO₂ on NEP from the MsTMIP models. We then used these two-box models to replicate the eCO₂ experiments with the form of changing CO₂ identical to that in the eCO₂ experiments.

First, we extracted the NPP, total biomass carbon and soil carbon variables of each MsTMIP model for each eCO₂ experimental site by averaging the values in the grid cells with the same dominant vegetation type as in the field experiment within a 4.5° × 4.5° window around the eCO₂ experimental site. We use a 4.5° × 4.5° window in this study partly because it is the minimum spatial scale for which the same vegetation type that, corresponding to the
FACE sites, is represented in all the MsTMIP models. The LULCC data set driving the MsTMIP simulations contained 47 classes of synergetic land-cover product (SYNMAP)³⁸,³⁹. We reclassified this data set ('biome_frac' variable) into the major life forms (SIMPLE legend) categories defined by Jung et al.³⁸ to obtain the spatial distributions of the forest and grassland fractions ('Trees' and 'Grasses' in the SIMPLE legend). Grid cells with mean fractions of trees or grasses during 1959-2010 not less than 50% were tagged 'forest' or 'grassland', respectively. Moreover, we analysed the relationship between the climate conditions in the 4.5° × 4.5° window and central 0.5° × 0.5° grid across the seven eCO₂ experiment sites (Supplementary Fig. 21). Results show that, as expected, there are large variations in climate conditions (temperature and precipitation) among the seven sites, but for each individual site the climate conditions (temperature and precipitation) in the 4.5° × 4.5° window are very close to those at the central 0.5° × 0.5° grid (Supplementary Fig. 21). The 4.5° × 4.5° window also generally captures the spatial variability of soil texture (clay, sand and silt fraction) across the seven eCO₂ experiment sites (Supplementary Fig. 22).

Second, we assumed that the total biomass-carbon pool ($C_B$) had a constant rate of decay ($\mu$) and that the soil-carbon pool ($C_S$) also had a constant rate of decay ($\rho$). The decaying carbon from the total biomass pool (B2S) enters the soil pool, and the decaying carbon from the soil pool is emitted to the atmosphere as HR. This simple system was represented by a linear two-box model:

$$ 
\begin{cases} 
\frac{d\Delta C_B}{dt} = \Delta NPP(t) - \mu \times \Delta C_B(t) \\ 
\frac{d\Delta C_S}{dt} = \mu \times \Delta C_B(t) - \rho \times \Delta C_S(t) 
\end{cases} 
\quad (2) 
$$

The analytical solution was:

$$ 
\left\{ 
\begin{aligned} 
\Delta C_B(t) &= \int_0^t [\exp(-\mu t') \times \Delta NPP(t-t')] dt' \\ 
\Delta B2S(t) &= \mu \times \Delta C_B(t) \\ 
\Delta C_S(t) &= \int_0^t [\exp(-\rho t') \times \Delta B2S(t-t')] dt' \\ 
\Delta HR(t) &= \rho \times \Delta C_S(t) 
\end{aligned} 
\right. 
\quad (3) 
$$

Where $\Delta NPP(t)$ is the CO₂-induced NPP (SG3 - SG2) in year $t$, $\Delta C_B(t)$ is the change in CO₂-induced total biomass-carbon storage (SG3 - SG2) in year $t$ relative to the first year, $\Delta C_S(t)$ is the change in CO₂-induced soil-carbon storage (SG3 - SG2) in year $t$ relative to the first year, $\mu$ is the constant decay rate of $\Delta C_B$ ($0 \le \mu \le 1$), $\Delta B2S(t)$ is the carbon flux from $\Delta C_B$ to $\Delta C_S$ in year $t$, $\rho$ is the constant decay rate of $\Delta C_S$ ($0 \le \rho \le 1$), $\Delta HR(t)$ is the CO₂-induced heterotrophic respiration in year $t$.
The parameters μ and ρ for each eCO₂ experimental site in each model were fitted with equation (3) using NPP, total biomass-carbon pool and soil-carbon pool for 1901–2010 from the MsTMIP outputs using Matlab (R2018a) software. The codes are shown in Supplementary Information. The two-box models generally emulated successfully the eCO₂-induced historical biomass and soil carbon evolution at the seven eCO₂ sites from the original complex process-oriented MsTMIP models (Supplementary Figs. 2–13). Nevertheless, the substitute models do not work well at reproducing the biomass-carbon pool evolution at grassland sites for 5 of 12 MsTMIP models (CLM4, GTEC, LPJ-wsl, ORCHIDEE-LSCE and SiBCASA at Duolun, and CLM4, LPJ-wsl and SiBCASA models at PHACE). This model underperformance is particularly when net primary productivity monotonically increases but biomass-carbon pool suddenly declines (Supplementary Figs. 2, 5, 7–9).

Third, we reproduced the eCO₂ experiment using the two-box model for each eCO₂ experimental site in each MsTMIP model, assuming that CO₂ abruptly increased by 60 ppm (from 320 ppm to 380 ppm) and was then held fixed. This analysis was performed with the following steps. (1) Calculate eCO₂-induced NPP in the reproduced eCO₂ experiment (ΔNPP<sup>site_FACE_M</sup>). We assumed that ΔNPP<sup>site_FACE_M</sup> was a constant and the same as the long-term response of NPP to eCO₂ in the MsTMIP simulations. Therefore,

$$ \Delta \text{NPP}^{\text{site\_FACE\_M}} = \frac{\Delta \text{NPP}}{\Delta \text{CO}_2} \times 60, $$

, where ΔNPP is the modelled difference in CO₂-induced NPP in the same months as the corresponding eCO₂ experiment between 2001–2010 and 1959–1968 and ΔCO₂ is the difference in the atmospheric CO₂ concentration between 2001–2010 and 1959–1968. (2) Calculate eCO₂-induced HR in the reproduced eCO₂ experiment (ΔHR<sup>site_FACE_M</sup>). We calculated ΔHR<sup>site_FACE_M</sup> using equation (3) for each MsTMIP model on the basis of its fitted parameters (μ and ρ) and ΔNPP<sup>site_FACE_M</sup>. (3) Calculate eCO₂-induced NEP in the reproduced FACE experiment (ΔNEP<sup>site_FACE_M</sup>) as

$$ \Delta \text{NEP}^{\text{site\_FACE\_M}} = \Delta \text{NPP}^{\text{site\_FACE\_M}} - \Delta \text{HR}^{\text{site\_FACE\_M}}. \quad (4) \text{ Calculate} $$

the sensitivity of NEP to eCO₂ in the reproduced eCO₂ experiment (B<sup>Mod</sup><sub>Site</sub>) as

$$ B_{\text{Site}}^{\text{Mod}} = \frac{\Delta \text{NEP}^{\text{site\_FACE\_M}}}{60} \quad (\text{Supplementary Fig. 1}). \quad (5) \text{ When comparing} $$

with FACE observations, the duration for the calculation of B<sup>Mod</sup><sub>Site</sub> is site specific and is the same as the duration in the real FACE experiment at each site. For example, the duration of Duke-FACE is nine years, from 1997 to 2005 (Supplementary Table 2). When comparing with FACE observations, we

used the mean of B<sup>Mod</sup><sub>Site</sub> in the nine years from the first to the ninth experimental years in the two-box-model-based FACE experiment at Duke-
FACE site (Supplementary Fig. 16). The simulated $B_{Site}$ in Fig. 2a,c,e is the mean $B_{Site}^{Mod}$ for the seven eCO₂ sites for each model.

Transform observed site-scale NEP sensitivity to elevated atmospheric CO₂ in eCO₂ experiments ($B_{Site}^{Obs}$) to the NEP sensitivity at a lower CO₂ concentration level ($B_{Site}^{Hist}$)

In the eCO₂ experiments, atmospheric CO₂ concentration ($c_a$) was increased from ambient (average ~381 ppm) to a very high concentration roughly representative of two times pre-industrial values (average ~588 ppm, indicated as ≈ 2 × CO₂). This change in $c_a$ is much larger than that used in MS-TMIP model simulations (from ~320 ppm during 1959-1968 to ~380 ppm during 2001-2010). To convert NEP responses at ≈ 2 × CO₂ $c_a$ in eCO₂ experiments to responses at historical $c_a$ values, for each eCO₂ experiment, we extrapolated the observed NEP sensitivity with CO₂ increasing from control $c_a$ ($c_{a2}$) to treatment $c_a$ ($c_{a2p}$) to the NEP sensitivity with CO₂ increasing from 320 ppm ($c_{a1}$) to 380 ppm ($c_{a1p}$).

Assuming that carbon use efficiency (CUE = NPP/GPP) and turnover rates of carbon pools in short-term eCO₂ experiments are unchanged when CO₂ varies from $c_{a1}$ to $c_{a1p}$ compared to from $c_{a2}$ to $c_{a2p}$ (Supplementary Figs. 23 and 24), eCO₂-induced NEP changes proportionally with eCO₂-induced GPP according to Supplementary equation (5), equivalent to the following equation:

$$ \frac{B_{NEP1}}{B_{NEP2}} = \frac{B_{GPP1}}{B_{GPP2}} \qquad (4) $$

where $B_{NEP1}$ and $B_{GPP1}$ are the CO₂ fertilization effects on NEP and GPP over $c_a$ ranging from $c_{a1}$ to $c_{a1p}$ and $B_{NEP2}$ and $B_{GPP2}$ are the CO₂ fertilization effects on NEP and GPP over $c_a$ ranging from $c_{a2}$ to $c_{a2p}$.

The response of GPP to rising $c_a$ can be estimated with a first-principles-based universal model of photosynthesis (Pmodel), which is extensively supported by flux observations and FACE experiment²⁷. The form of Pmodel is a 'light-use efficiency' (LUE) model that predicted GPP is proportional to absorbed photosynthetically active radiation (PAR), where the LUE is predicted from first principles²⁷. The standard biochemical model of photosynthesis proposed by Farquhar et al.⁴⁰ predicts the instantaneous rates of photosynthesis limited by Rubisco ($A_c$) and electron-transport ($A_j$), respectively. The coordination or colimitation hypothesis states that the maximum capacity of carboxylation acclimates to the prevailing environmental conditions at weekly or longer time scales, allowing the two
photosynthetic processes of A_c and A_j to be coordinated with each other (A_j = A_c) under typical daytime conditions<sup>27,41,42</sup>. The LUE model descriptions are in the methods of ref. 27 and more detailed information can be found in section 5 of "A light-use efficiency model for GPP" in the supplementary information. The calculation of GPP is shown as following equations (5) and (6) (equations [2] and [3] in ref. 27).

$$GPP = \varphi_0 I_{\text{abs}} m \sqrt{[1 - (c^*/m)^{2/3}]} \quad (5)$$

where

$$m = (c_a - \Gamma^*) / \left\{ c_a + 2\Gamma^* + 3\Gamma^* \sqrt{[1.6\eta^* D_0 \beta^{-1} (K + \Gamma^*)^{-1}]} \right\} \quad (6)$$

Here, c_a is atmospheric CO₂ concentration (ppm), φ₀ is the intrinsic quantum yield (g C mol⁻¹), m is a component of light-use efficiency, I_abs is the absorbed photosynthetic photon flux density (mol m⁻² s⁻¹), Γ^* is the photorespiratory compensation point (Pa), K is the effective Michaelis-Menten coefficient of Rubisco (Pa), η^* is the viscosity of water relative to its value at 25 °C, D₀ is vapour pressure deficit (Pa), and c^* and β represent the cost factor of the maintaining electron transport capacity (~0.41) and the ratio of carboxylation to transpiration cost factors (~240). Both are estimated from independent observational data.

According to this model, the ratio between GPP sensitivities to elevated CO₂ from c_a1 to c_a1p and from c_a2 to c_a2p (Supplementary Fig. 25) is estimated as:

$$\frac{B_{\text{GPP1}}}{B_{\text{GPP2}}} = \frac{m_{1\text{p}}\sqrt{[1 - (c^*/m_{1\text{p}})^{2/3}]} - m_1\sqrt{[1 - (c^*/m_1)^{2/3}]} \times \frac{(c_{\text{a2p}} - c_{\text{a2}})}{(c_{\text{a1p}} - c_{\text{a1}})}}{m_{2\text{p}}\sqrt{[1 - (c^*/m_{2\text{p}})^{2/3}]} - m_2\sqrt{[1 - (c^*/m_2)^{2/3}]} \times \frac{(c_{\text{a2p}} - c_{\text{a2}})}{(c_{\text{a1p}} - c_{\text{a1}})}} \quad (7)$$
Here, $m_1$, $m_{1p}$, $m_2$ and $m_{2p}$ are m calculated at $c_{a1}$, $c_{a1p}$, $c_{a2}$ and $c_{a2p}$, respectively, using equation (6). We estimated m from growing-season mean temperature (tmp), diurnal temperature range (dtr), relative humidity (reh) and elevation (elv) extracted from CRU climatology (CRU CL v. 2.0) at 10-minute resolution$^{43}$.

The observed NEP sensitivity measured at ≈ 2 × CO₂ in eCO₂ experiments ($B_{Site}^{Obs}$) was extrapolated to the sensitivity at CO₂ ranging from 320 to 380 ppm in the simulations ($B_{Site}^{Hist}$):

$$B_{Site}^{Hist} = B_{Site}^{Obs} \times \frac{B_{GPP1}}{B_{GPP2}} \qquad (8)$$

Probability density function of $B_{Site}^{Hist}$ from eCO₂ experiments

The average $B_{Site}^{Hist}$ from the eCO₂ experiments and its uncertainty were estimated using a bootstrap approach with the following steps. (1) One $B_{Site}^{Hist}$ was randomly drawn for each eCO₂ experiment site, and the selected $B_{Site}^{Hist}$ from all eCO₂ experiment sites were then averaged to obtain $B_{Site}^{Hist}$. (2) Step (1) was performed 1,000 times to obtain 1,000 $B_{Site}^{Hist}$. (3) The mean and standard deviation of the 1,000 $B_{Site}^{Hist}$ from step (2) were calculated to obtain mean $B_{Site}^{Hist}$ ($\overline{B_{Site}^{Hist}}$) and its uncertainty $[\sigma(B_{Site}^{Hist})]$ from the eCO₂ experiments. This uncertainty represented the site-scale uncertainty due to experimental procedures, such as uneven measurements and sample time, rather than to spatial heterogeneity across sites. The uncertainty could

therefore be comparable to the mean of simulated $B_{Site}$ ($B_{Site}^{Mod}$) for the seven eCO₂ experiment sites for each model (Fig. 2a,c,e). We assumed that all observations could be represented by a Gaussian distribution, with mean and standard deviation obtained using this bootstrap approach.

Least-squares linear regression between $B_{NH}^{Mod}$ and $B_{Site}^{Mod}$ across the MS-TMIP models
This study used an emergent-constraint approach proposed by Cox et al.¹⁸.

We identified a linear least-squares regression between $B_{NH}^{\text{Mod}}$ and $B_{Site}^{\text{Mod}}$ across the MsTMIP models:

$$B_{NH}^{\text{Mod}}(i) = a \times B_{Site}^{\text{Mod}}(i) + b \quad (9)$$

where $B_{NH}^{\text{Mod}}(i)$ is $B_{NH}^{\text{Mod}}$ of model i, and $B_{Site}^{\text{Mod}}(i)$ is $B_{Site}^{\text{Mod}}$ of model i.

The least-squares error of the regression model was calculated as:

$$s^2 = \frac{1}{N-2} \sum_{i=1}^{N} \left[ B_{NH}^{\text{Mod}}(i) - \tilde{B}_{NH}^{\text{Mod}}(i) \right]^2 \quad (10)$$

where N is the number of models used in the regression model, $B_{NH}^{\text{Mod}}(i)$ is $B_{NH}^{\text{Mod}}$ of model i, and $\tilde{B}_{NH}^{\text{Mod}}(i)$ is the predicted $B_{NH}^{\text{Mod}}$ using the regression model at $B_{Site}^{\text{Mod}}$ of model i.

The 'prediction error' of $\tilde{B}_{NH}^{\text{Mod}}$ based on the regression model for a given x was:

$$\sigma(\tilde{B}_{NH}^{\text{Mod}}|x) = s\sqrt{1 + \frac{1}{N} + \frac{(x - \bar{x})^2}{N\sigma_x^2}} \quad (11)$$

where $\bar{x}$ is the mean of $x_1, x_2 ... x_N$ and $\sigma_x$ is the standard deviation of $x_1, x_2 ... x_N$.

The probability density function of $B_{NH}^{\text{Mod}}$ for a given x was:
$$P_{\text{reg}}(B_{\text{NH}}^{\text{Mod}}|x) = \frac{1}{\sqrt{2\pi\sigma(\tilde{B}_{\text{NH}}^{\text{Mod}}|x)^2}} e^{-\frac{[\tilde{B}_{\text{NH}}^{\text{Mod}} - (\tilde{B}_{\text{NH}}^{\text{Mod}}|x)]^2}{2\sigma(\tilde{B}_{\text{NH}}^{\text{Mod}}|x)^2}} \quad (12)$$

Constraint on $B_{NH}$

The probability density function for $B_{NH}$ derived from observations ($B_{NH}^{\text{Constrained}}$), the red lines in Fig. 2b,d,f, was calculated as:

$$P(B_{\text{NH}}^{\text{Constrained}}) = \int_{-\infty}^{+\infty} P_{\text{reg}}(B_{\text{NH}}^{\text{Mod}}|x) \times \frac{1}{\sqrt{2\pi\sigma(B_{\text{Site}}^{\text{Hist}})^2}} e^{-\frac{(x-B_{\text{Site}}^{\text{Hist}})^2}{2\sigma(B_{\text{Site}}^{\text{Hist}})^2}} dx \quad (13)$$

Analysis of the robustness of $B_{NH}^{\text{Constrained}}$ when using fewer than seven eCO₂ experiment sites

We performed the following steps. (1) We randomly selected N eCO₂ experiment sites to calculate the constrained $B_{NH}$ using all 12 models. (2) Step (1) was performed 1,000 times to obtain 1,000 constrained $B_{NH}$. (3) All constrained $B_{NH}$ in step (2) were averaged to obtain the mean constrained $B_{NH}$. We took N = 6, 5 and 4 (Supplementary Fig. 14).

Analysis of the robustness of $B_{NH}^{\text{Constrained}}$ when using a logarithmic function

The responses of GPP to rising $c_a$ can be estimated with a logarithmic function (equation 15 in ref. 3) that was first introduced by Bacastow et al.²⁸:

$$P_t = P_0 \left[ 1 + B_{\text{log}} \ln \left( \frac{C_t}{C_0} \right) \right] \quad (14)$$

Here, $P_t$ and $P_0$ are ecosystem carbon fluxes (GPP, NPP or NEP) for $c_t$ and $c_0$ CO₂ concentration, $B_{\text{log}}$ is a constant factor. Thus, $B_{\text{GPP1}}/B_{\text{GPP2}}$ in equation (8) is expressed in this case as:
$$ \frac{B_{GPP1}}{B_{GPP2}} = \frac{\{GPP_0[1+B_{\text{log}}\ln(c_{a1p}/c_{a0})] - GPP_0[1+B_{\text{log}}\ln(c_{a1}/c_{a0})]\}/(c_{a1p}-c_{a1})}{\{GPP_0[1+B_{\text{log}}\ln(c_{a2p}/c_{a0})] - GPP_0[1+B_{\text{log}}\ln(c_{a2}/c_{a0})]\}/(c_{a2p}-c_{a2})} \\ = \frac{[\ln(c_{a1p}) - \ln(c_{a1})]/(c_{a1p}-c_{a1})}{[\ln(c_{a2p}) - \ln(c_{a2})]/(c_{a2p}-c_{a2})} \qquad (15) $$

We analysed the emergent constraints on the sensitivity of the Northern Hemisphere terrestrial carbon sink to elevated atmospheric CO₂ concentration using equation (15), rather than equation (7), with the other conditions unchanged (Supplementary Fig. 15).

Constraining the sensitivity of global terrestrial carbon sink to eCO₂

Using all 12 MsTMIP models and the observations at the seven eCO₂ experiment sites, we analysed the constrained sensitivity of global terrestrial

carbon sink to eCO₂ during 1959–2010 (B<sup>Constrained</sup><sub>Globe</sub>). The approach was the

same as for the calculation of B<sup>Constrained</sup><sub>NH</sub> but using the sensitivity of global terrestrial carbon sink to eCO₂ (B<sub>Globe</sub>) (Supplementary Fig. 17).

Constraining the sensitivity of global terrestrial carbon sink to long-term climate change (Γ<sub>Long</sub>)

We identified an emergent relationship between Γ<sub>Long</sub> and the sensitivity of global terrestrial carbon sink to interannual variability in climate (Γ<sub>IAV</sub>) across 13 MsTMIP models (BIOME-BGC, CLM4, CLM4VIC, DLEM, GTEC, ISAM, LPJ-wsl, ORCHIDEE-LSCE, SiB3, SiBCASA, TEM6, VEGAS2.1 and VISIT) with SG1 and SG3/BG1 (SG3 for carbon-only models and BG1 for C-N models) outputs (Supplementary Fig. 20), similar to Huntzinger et al.⁶ but focusing on the Γ<sub>Long</sub> during 1959–2010. CLASS-CTEM-N was not used following Ito et al.⁴⁴ because its net CO₂ exchange was always negative. TRIPLEX-GHG was not used because SG1 simulation is not available from https://doi.org/10.3334/ORNLDAAC/1225 (ref. 35).

Γ<sub>Long</sub> and Γ<sub>IAV</sub> were calculated using the approach from Huntzinger et al.⁶. Γ<sub>Long</sub> was calculated by regressing the land annual NEE given by the SG1 scenario against temperature and precipitation during 1959–2010 from the MsTMIP driver data set. The regression coefficient on temperature was used to determine Γ<sub>Long</sub>. Γ<sub>IAV</sub> was calculated by regressing the detrended land annual NEE given by the SG3 (for carbon-only models) or BG1 (for C-N models) scenario against detrended temperature during 1959–2010 from the MsTMIP driver data set. The detrended NEE and the detrended temperature were calculated by subtracting 11-year running means. The 11-year running means were calculated using ‘runmean’ function (https://www.mathworks.com/matlabcentral/fileexchange/10113-runmean, last accessed June 25, 2018) with its parameters ‘M’ and ‘MODESTR’ set as
‘5’ and ‘edge’, respectively. The post-volcano years (1963, 1964, 1982, 1983, 1991 and 1992) were removed before calculating both Γ<sub>Long</sub> and Γ<sub>IAV</sub>.

The Γ<sub>IAV</sub> from GCB 2017 (refs. 2·37) was used as an observational constraint
on the Γ<sub>Long</sub> (Supplementary Fig. 20). The GCB Γ<sub>IAV</sub> was calculated by
regressing the detrended GCB net land sink (NLS) against detrended
temperature from the MsTMIP driver data set, using the same approach as
the calculation of Γ<sub>IAV</sub> in MsTMIP models. The GCB NLS was estimated by
subtracting the growth rate in atmospheric CO<sub>2</sub> concentration and ocean CO<sub>2</sub>sink from global emissions from fossil fuels and industry: NLS = E<sub>FF</sub> - G<sub>ATM</sub> - S<sub>OCEAN</sub>. The GCB Γ<sub>IAV</sub> is estimated as -2.36 ± 0.77 PgC yr<sup>-1</sup> K<sup>-1</sup> (mean ± 1 standard error). Using this as an observational constraint, the constrained Γ<sub>Long</sub> is estimated as -1.57 ± 0.51 PgC yr<sup>-1</sup> K<sup>-1</sup> on the basis of the linear relationship between Γ<sub>Long</sub> and Γ<sub>IAV</sub> across MsTMIP models (Supplementary Fig. 20).

Data availability

Driver data of MsTMIP models are available from
https://doi.org/10.3334/ORNLDAAC/1220. The outputs of MsTMIP models are
available from https://doi.org/10.3334/ORNLDAAC/1225. Global Carbon
Budget 2017 data (Global Carbon Budget 2017v1.3.xlsx) are available from
https://doi.org/10.18160/GCP-2017. CRU climatology data (CRU CL v. 2.0) are
available from https:// rudata.uea.ac.uk/cru/data/hrg/.

Code availability

The code for fitting the two-box model, as given in equation (3), is shown in
Supplementary Information.

References

1.

Field, C. B. & Raupach, M. R. *The Global Carbon Cycle: Integrating Humans, Climate, and the Natural World* (Island Press, 2004).

2.

Le Quéré, C. et al. Global Carbon Budget 2017. *Earth Syst. Sci. Data* 10, 405–448 (2018).

3.

Friedlingstein, P. et al. On the contribution of CO₂ fertilization to the missing biospheric sink. *Glob. Biogeochem. Cycles* 9, 541–556 (1995).

4.

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

5.
Schimel, D., Stephens, B. B. & Fisher, J. B. Effect of increasing CO₂ on the terrestrial carbon cycle. *Proc. Natl. Acad. Sci. USA* 112, 436-441 (2015).

6.

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).

7.

Keenan, T. F. et al. Recent pause in the growth rate of atmospheric CO₂ due to enhanced terrestrial carbon uptake. *Nat. Commun.* 7, 13428 (2016).

8.

Sitch, S. et al. Recent trends and drivers of regional sources and sinks of carbon dioxide. *Biogeosciences* 12, 653-679 (2015).

9.

Ilyina, T. & Friedlingstein, P. WCRP Grand Challenge: Carbon feedbacks in the climate system. White paper (37th Session of the WCRP Joint Scientific Committee, 2016); https://www.wcrp-climate.org/component/content/article/872-carbon-feedbacks-overview?catid=139&Itemid=539

10.

Norby, R. J. & Zak, D. R. Ecological lessons from free-air CO₂ enrichment (FACE) experiments. *Annu. Rev. Ecol. Evol. Syst.* 42, 181-203 (2011).

11.

Zaehle, S. et al. Evaluation of 11 terrestrial carbon-nitrogen cycle models against observations from two temperate free-air CO₂ enrichment studies. *New Phytol.* 202, 803-822 (2014).

12.

Jastrow, J. D. et al. Elevated atmospheric carbon dioxide increases soil carbon. *Glob. Change Biol.* 11, 2057-2064 (2005).

13.

De Graaff, M. A., van Groenigen, K. J., Six, J., Hungate, B. & van Kessel, C. Interactions between plant growth and soil nutrient cycling under elevated CO₂: a meta-analysis. *Glob. Change Biol.* 12, 2077-2091 (2006).

14.

Hickler, T. et al. CO₂ fertilization in temperate FACE experiments not representative of boreal and tropical forests. *Glob. Change Biol.* 14, 1531-1542 (2008).

15.
Finzi, A. C. et al. Increases in nitrogen uptake rather than nitrogen-use efficiency support higher rates of temperate forest productivity under elevated CO₂. *Proc. Natl. Acad. Sci. USA* 104, 14014-14019 (2007).

16.

Drake, J. E. et al. Increases in the flux of carbon belowground stimulate nitrogen uptake and sustain the long-term enhancement of forest productivity under elevated CO₂. *Ecol. Lett.* 14, 349-357 (2011).

17.

Terrer, C., Vicca, S., Hungate, B. A., Phillips, R. P. & Prentice, I. C. Mycorrhizal association as a primary control of the CO₂ fertilization effect. *Science* 353, 72-74 (2016).

18.

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

19.

Wenzel, S., Cox, P. M., Eyring, V. & Friedlingstein, P. Projected land photosynthesis constrained by changes in the seasonal cycle of atmospheric CO₂. *Nature* 538, 499-501 (2016).

20.

Mystakidis, S., Davin, E. L., Gruber, N. & Seneviratne, S. I. Constraining future terrestrial carbon cycle projections using observation-based water and carbon flux estimates. *Glob. Change Biol.* 22, 2198-2215 (2016).

21.

Zhao, C. et al. Plausible rice yield losses under future climate warming. *Nat. Plants* 3, 16202 (2016).

22.

Lian, X. et al. Partitioning global land evapotranspiration using CMIP5 models constrained by observations. *Nat. Clim. Change* 8, 640-646 (2018).

23.

Huntzinger, D. N. et al. The North American Carbon Program Multi-scale Synthesis and Terrestrial Model Intercomparison Project - Part 1: overview and experimental design. *Geosci. Model Dev.* 6, 2121-2133 (2013).

24.

Norby, R. J., Warren, J. M., Iversen, C. M., Medlyn, B. E. & McMurtrie, R. E. CO₂ enhancement of forest productivity constrained by limited nitrogen availability. *Proc. Natl. Acad. Sci. USA* 107, 19368-19373 (2010).

25.
De Kauwe, M. G. et al. Where does the carbon go? A model-data intercomparison of vegetation carbon allocation and turnover processes at two temperate forest free-air CO₂ enrichment sites. *New Phytol.* 203, 883-899 (2014).

26.

Franks, P. J. et al. Sensitivity of plants to changing atmospheric CO₂ concentration: from the geological past to the next century. *New Phytol.* 197, 1077-1094 (2013).

27.

Wang, H. et al. Towards a universal model for carbon dioxide uptake by plants. *Nat. Plants* 3, 734-741 (2017).

28.

Bacastow, R. & Keeling, C. D. in *Carbon and the Biosphere* (eds Woodwell, G. M. & Pecan, E. V.) 86-135 (US Dep. of Comm., 1973).

29.

Peñuelas, J. et al. Human-induced nitrogen-phosphorus imbalances alter natural and managed ecosystems across the globe. *Nat. Commun.* 4, 2934 (2013).

30.

McCarthy, H. R. et al. Re-assessment of plant carbon dynamics at the Duke free-air CO₂ enrichment site: interactions of atmospheric [CO₂] with nitrogen and water availability over stand development. *New Phytol.* 185, 514-528 (2010).

31.

Reich, P. B. & Hobbie, S. E. Decade-long soil nitrogen constraint on the CO₂ fertilization of plant biomass. *Nat. Clim. Change* 3, 278-282 (2013).

32.

Norby, R. J. et al. Model-data synthesis for the next generation of forest free-air CO₂ enrichment (FACE) experiments. *New Phytol.* 209, 17-28 (2016).

33.

Drake, J. E. et al. Short-term carbon cycling responses of a mature eucalypt woodland to gradual stepwise enrichment of atmospheric CO₂ concentration. *Glob. Change Biol.* 22, 380-390 (2016).

34.

Medlyn, B. E. et al. Using ecosystem experiments to improve vegetation models. *Nat. Clim. Change* 5, 528-534 (2015).

35.
Huntzinger, D. N. et al. NACP MsTMIP: Global 0.5-deg Terrestrial Biosphere Model Outputs in Standard Format. v.1.0 (ORNL DAAC, 2016); https://doi.org/10.3334/ORNLDAAC/1225

36.

Wei, Y. et al. NACP MsTMIP: Global and North American Driver Data for Multi-model Intercomparison (ORNL DAAC, 2014); https://doi.org/10.3334/ORNLDAAC/1220

37.

Supplemental data of Global Carbon Budget 2017 v.1.0 (Global Carbon Project, 2017); https://doi.org/10.18160/gcp-2017

38.

Jung, M., Henkel, K., Herold, M. & Churkina, G. Exploiting synergies of global land cover products for carbon cycle modeling. *Remote Sens. Environ.* 101, 534-553 (2006).

39.

Wei, Y. et al. The North American Carbon Program Multi-scale Synthesis and Terrestrial Model Intercomparison Project: part 2: environmental driver data. *Geosci. Model Dev.* 7, 2875-2893 (2014).

40.

Farquhar, G. D., Caemmerer, S. V. & Berry, J. A. A biochemical model of photosynthetic CO₂ assimilation in leaves of C₃ species. *Planta* 149, 78-90 (1980).

41.

Maire, V. et al. The coordination of leaf photosynthesis links C and N fluxes in C₃ plant species. *PLoS ONE* 7, e38345 (2012).

42.

Togashi, H. F. et al. Thermal acclimation of leaf photosynthetic traits in an evergreen woodland, consistent with the coordination hypothesis. *Biogeosciences* 15, 3461-3474 (2018).

43.

Mark, N., David, L., Mike, H. & Ian, M. A high-resolution data set of surface climate over global land areas. *Clim. Res.* 21, 1-25 (2002).

44.

Ito, A. et al. Decadal trends in the seasonal-cycle amplitude of terrestrial CO₂ exchange resulting from the ensemble of terrestrial biosphere models. *Tellus B Chem. Phys. Meteor.* 68, 28968 (2016).

Acknowledgements
This study was supported by the Strategic Priority Research Program (A) of
the Chinese Academy of Sciences (Grant XDA20050101), the National Key
R&D Program of China (2017YFA0604702), and the National Natural Science
Foundation of China (41861134036, 41701089). The full list of
acknowledgements is provided in the Supplementary Information.