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

HAL
open science

# Global vegetation biomass production efficiency constrained by models and observations

Yue He, Shushi Peng, Yongwen Liu, Xiangyi Li, Kai Wang, Philippe Ciais, M.

Altaf Arain, Yuanyuan Fang, Joshua Fisher, Daniel S Goll, et al.

► To cite this version:

Yue He, Shushi Peng, Yongwen Liu, Xiangyi Li, Kai Wang, et al. Global vegetation biomass production efficiency constrained by models and observations. Global Change Biology, 2020, 26 (3), pp.1474-1484. 10.1111/gcb.14816. hal-02899467

HAL Id: hal-02899467

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

Submitted on 4 Aug 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Article Type: Primary Research Articles

2

3 **Global vegetation biomass production efficiency constrained by models and observations**

4 Yue He¹, Shushi Peng¹*, Yongwen Liu¹,²,³, Xiangyi Li¹, Kai Wang¹, Philippe Ciais⁴, M. Altaf

5 Arain⁵, Yuanyuan Fang⁶, Joshua B. Fisher⁷, Daniel Goll⁴, Daniel Hayes⁸, Deborah N.

6 Huntzinger⁹, Akihiko Ito¹⁰,¹¹, Atul Jain¹², Ivan A. Janssens¹³, Jiafu Mao¹⁴, Campioli Matteo¹³,

7 Anna M. Michalak⁶, Changhui Peng¹⁵,¹⁶, Josep Peñuelas¹⁷,¹⁸, Benjamin Poulter¹⁹, Dahe

8 Qin²⁰,²¹, Daniel M Ricciuto¹⁴, Kevin Schaefer²², Christopher R. Schwalm²³,²⁴, Xiaoying Shi¹⁵,

9 Hanqin Tian²⁵, Sara Vicca¹³, Yaxing Wei¹⁴, Ning Zeng²⁶, Qiuan Zhu¹⁵,¹⁶

10

11 ¹Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University, Beijing 100871, China.

12 ²Key Laboratory of Alpine Ecology, Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing 100101, China.

13 ³CAS Center for Excellence in Tibetan Earth Science, Chinese Academy of Sciences, Beijing 100101, China.

14 ⁴Laboratoire des Sciences du Climat et de l'Environnement, CEA CNRS UVSQ, 91191 Gif-sur-Yvette, France.

15 ⁵School of Geography and Earth Sciences and McMaster Centre for Climate Change, McMaster University, Hamilton, Ontario L8S 4K1, Canada.

16 ⁶Department of Global Ecology, Carnegie Institution for Science, Stanford, California 94305, USA.

17 ⁷Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California 91109, USA.

18 ⁸School of Forest Resources, University of Maine, Orono, Maine 04469, USA.
26 ⁹School of Earth Sciences and Environmental Sustainability, Northern Arizona University, Flagstaff, Arizona 86011, USA.

28 ¹⁰National Institute for Environmental Studies, Tsukuba 305-8506, Japan.

29 ¹¹Japan Agency for Marine-Earth Science and Technology, Yokohama 236-0001, Japan.

30 ¹²Department of Atmospheric Sciences, University of Illinois, Urbana, Illinois 61801, USA.

31 ¹³Centre of Excellence PLECO (Plant and Vegetation Ecology), Department of Biology, University of Antwerp, 2610 Wilrijk, Belgium.

33 ¹⁴Environmental Sciences Division and Climate Change Science Institute, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA.

35 ¹⁵Institute of Environment Sciences, Biology Science Department, University of Quebec at Montreal, Montreal, Quebec H3C 3P8, Canada.

37 ¹⁶State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau, College of Forestry, Northwest A & F University, Yangling 712100, China.

39 ¹⁷CSIC, Global Ecology Unit CREAF-CEAB-UAB, Cerdanyola del Vallés, Barcelona 08193, Catalonia, Spain.

41 ¹⁸CREAF, Cerdanyola del Vallés, Barcelona, 08193, Catalonia, Spain.

42 ¹⁹Montana State University, Institute on Ecosystems and the Department of Ecology, Bozeman, Montana 59717, USA.

44 ²⁰State Key Laboratory of Cryospheric Sciences, Northwest Institute of Eco-environment and Resources, Chinese Academy of Sciences, Lanzhou 730000, China

46 ²¹National Climate Center, China Meteorological Administration, Beijing 100081, China

47 ²²National Snow and Ice Data Center, Boulder, Colorado 80309, USA.

48 ²³Woods Hole Research Center, 149 Woods Hole Road, Falmouth, Massachusetts 02540, USA.

50 ²⁴Center for Ecosystem Science and Society, Northern Arizona University, Flagstaff, Arizona
51 86011, USA.

52 <sup>25</sup>International Center for Climate and Global Change Research and School of Forestry and

53 Wildlife Sciences, Auburn University, Auburn, Alabama 36849, USA.

54 <sup>26</sup>Department of Atmospheric and Oceanic Science, University of Maryland, College Park,

55 Maryland 20740, USA.

56

57

58 Running title: Emergent constraint on global vegetation BPE

59 Keywords: biomass production, BPE, emergent constraint, carbon sink, terrestrial

60 biosphere model

61

62

Manuscript for *Global Change Biology*

63

64 * Correspondence: Dr. Shushi Peng, Email: speng@pku.edu.cn

65

66

67
68 **Abstract**

Plants use only a fraction of their photosynthetically derived carbon for biomass production. The biomass production efficiency (BPE), defined as the ratio of biomass production to photosynthesis, and its variation across and within vegetation types is poorly understood, which hinders our capacity to accurately estimate carbon turnover times and carbon sinks. Here, we present a new global estimation of BPE obtained by combining field measurements from 113 sites with 14 carbon cycle models. Our best estimate of global BPE is 0.41 ± 0.05, excluding cropland. The largest BPE is found in boreal forests (0.48 ± 0.06) and the lowest in tropical forests (0.40 ± 0.04). Carbon cycle models overestimate BPE, although models with carbon-nitrogen interactions tend to be more realistic. Using observation-based estimates of global photosynthesis, we quantify the global biomass production of non-cropland ecosystems of 41 ± 6 Pg C yr⁻¹. This flux is less than net primary production as it does not contain carbon allocated to symbionts, used for exudates or volatile carbon compound emissions to the atmosphere. Our study reveals a positive bias of 24 ± 11% in the model estimated biomass production (10 of 14 models). When correcting models for this bias while leaving modeled carbon turnover times unchanged, we found that the global ecosystem carbon storage change during the last century is decreased by 67% (or 58 Pg C).

85
86

**Introduction**

87

Biomass production (BP) consists of photosynthetically derived carbon used for biomass growth, i.e., to buildup leaves, wood and roots (Vicca *et al.*, 2012). BP is smaller than net primary productivity (NPP), which is the sum of BP and carbon allocated to non-structural organic compounds that are used to maintain functions of rhizosphere and plant (e.g., mycorrhizae, root exudates, and volatile carbon compound emissions) (Chapin *et al.*, 2006; Vicca *et al.*, 2012; Campioli *et al.*, 2016). A large fraction of BP (about 8.2 Pg C yr⁻¹) is being appropriated by humans through the harvest of food, wood and fibers (Haberl *et al.*, 2007; Erb *et al.*, 2016). The remaining BP enters into a cascade of ecosystem carbon pools, and determines the land carbon balance (Friend *et al.*, 2014). Currently non-harvested BP forces a carbon sink offsetting one third of anthropogenic emissions (Le Quéré *et al.*, 2015). The fate of this carbon sink, however, is highly uncertain under climate change and anthropogenic activities (Ciais *et al.*, 2013; Peñuelas *et al.*, 2017). In particular, since BP is a fraction of photosynthetically assimilated carbon (gross primary production, GPP), both dynamics in GPP and the fraction of GPP allocated to BP could affect the outcome of BP.

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

Progress has been made to derive observation-based GPP estimates from eddy-covariance networks (Beer *et al.*, 2010), and from new remote sensing data (e.g., solar-induced chlorophyll fluorescence, Guanter *et al.*, 2014; near-infrared reflectance of vegetation, Badgley *et al.*, 2017). However, we still know little about the allocation of GPP to BP, hereafter defined as biomass production efficiency (BPE), i.e., the ratio of BP to GPP. This is because *in situ* estimates of BP are rarely at the same locations as GPP. Earlier studies usually suggest a constant BPE universally (McCree & Troughto, 1966; Gifford, 1995; Waring *et al.*, 1998). However, more recent syntheses of both field measurements and remote sensing products all indicate that BPE varies greatly across vegetation types, environmental

103

104

105

106

107

108

109

110
conditions, ecosystem management intensity and soil fertility gradients (Zhang et al., 2009, 2014; Vicca et al., 2012; Fernández-Martínez et al., 2014; Campioli et al., 2015). This spatial heterogeneity of BPE hinders the up-scaling of this quantity to a global scale.

114

Terrestrial carbon cycle models simulate the dynamics of CO₂ fluxes and carbon pools in response to variations in climate and atmospheric composition, and can be run with climate scenarios to project future carbon storage change (Cramer et al., 2001; Friedlingstein et al., 2006; Piao et al., 2008; Yao et al., 2018). Previous models often assumed a constant fraction of autotrophic respiration (Ra) to GPP as a simplifying concept (e.g., Landsberg & Waring 1997; Sands et al., 2000; Veroustraete et al., 2002; Nemani et al., 2009), whereas ecosystem models at present employ more detailed descriptions of the carbon cycle (e.g., Krinner et al., 2005; Peng et al., 2011; Tian et al., 2012; Zhu et al., 2014). These latter carbon cycle models have structural similarities in the equations of leaf scale photosynthesis (Rogers et al., 2017) and soil carbon dynamics (Luo et al., 2016), but they differ greatly regarding how GPP assimilates are allocated to different plant organs or respired for maintenance and growth purposes, the processes that determine BPE. In this respect, BPE in a model is not a parameter, but rather a diagnostic variable that can be calculated from model outputs and reflects plant processes driven by environmental conditions and ecological stoichiometry in the case of models integrating carbon-nutrient interactions. Note that in current models BP = NPP, because mycorrhizae and exudates are not considered.

131

In this study, we aim to quantify the global value of BPE by applying an emergent constraint approach (Cox et al., 2013; Wenzel et al., 2014; Zhao et al., 2016; Kwiatkowski et al., 2017) to model outputs from a recent dataset of BPE measured globally (Campioli et al., 2015). This approach can reduce the inter-model uncertainties by identifying a robust linear

111

112

113

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

132

133

134

135
inter-model relationship between one quantity that cannot be directly measured and the other quantity with available field observations. When combined with observations of the second quantity, this relationship can be used to constrain the first quantity (see Methods). The set of global carbon cycle models we used is known as the Multi-scale Synthesis and Terrestrial Model Intercomparison Project (MsTMIP) (Huntzinger et al., 2013; Wei et al., 2014; see model list in Table S1). We also extend this analysis to an evaluation of the model performance in simulating global biomass production, and ecosystem carbon storage change during 1901-2010, with a distinction between carbon-only and carbon-nitrogen interactions.

144

## 145 Materials and methods

### 146 BPE data sets

The newly-assembled biomass production efficiency (BPE) dataset developed by Campioli et al. (2015) contains 131 field observations, including forests, grasslands, croplands, wetlands and tundra. For the current study, 18 sites were excluded because carbon cycle models do not have an explicit consideration of wetlands and tundra at present. In total, 113 observation-based BPE values covering almost all major vegetation types were obtained, including 33 data points of boreal forests (BoF), 31 data points of temperate forests (TeF), 5 data points of tropical forests (TrF), 20 data points of grasslands, and 24 data points of croplands.

155

A binary management classification scheme was also adopted in our study (Campioli et al., 2015). This classification scheme contained both managed status (when the site was dominated by human activity, such as thinning, harvesting or planting for forests, or newly established and fertilized grasslands) and unmanaged status (when the site was dominated by natural processes with a low human activity impacts). The anthropogenic ecosystem, i.e.,
cropland, is inherently under managed conditions. This BPE dataset also contained
information of site fertility, which was divided into three levels (i.e., low, medium or high
nutrient availability), according to information in previous studies or provided by associated
authors (Campioli *et al.*, 2015). In addition, climate data (i.e., mean annual temperature and
mean annual precipitation) for each site was retrieved from the WorldClim database
(http://worldclim.org/version2) (Fick & Hijmans, 2017).

167

The uncertainty ranges of biomass production (BP) and gross primary production (GPP)
for each site were estimated based on method-specific reduction factor determined by expert
judgment, the number of measurement years and vegetation type (Luyssaert et al., 2007;
Campioli et al., 2015). Campioli et al. (2015) calculated the uncertainty range of BPE for
each site by the means of the error propagation theory. And we then used this error
propagation method to obtain the statistical uncertainty of averaged BPE across all sites (e.g.,
the red shaded area shown in Fig 1a represents the uncertainty range of site-level mean BPE).

175

176 GPP datasets

We used three different data-driven GPP datasets in our study, as they are often
considered as benchmark for model evaluation (Piao et al., 2013; Li et al., 2019). The first
one is an eddy covariance flux tower data-driven GPP product at a spatial resolution of 0.5°
during 1982-2011 (Jung et al., 2011). The second one is Moderate Resolution Imaging
Spectroradiometer (MODIS) GPP product (collection 055) at 1 km resolution over 2000-2012
(Zhao et al., 2005). The third one is derived from a process-based model, the Breathing Earth
System Simulator with a spatial resolution of 1 km over 2000-2015 (Jiang & Ryu, 2016).
Here, all these GPP datasets were regridded to a common 0.5° grid to match the carbon cycle
models used in our study. We separately calculated the mean annual global GPP during their
177
178
179
180
181
182
183
184
185
186

observation period, and the average of them was further used in the calculation of global BP
(excluding cropland).

187

188

189

Carbon cycle models

190

We used modeling results from the Multi-scale Synthesis and Terrestrial Model Intercomparison Project (MsTMIP) (Huntzinger et al., 2013; Wei et al., 2014) to calculate model-derived BPE. All model simulations participating in this project were performed with the same protocol for the period of 1901-2010 and at a 0.5° spatial resolution. Here we focused on model outputs over the past three decades (1981-2010) from the ‘SG3’ simulation, which was forced with time-varying climate, atmospheric CO2 concentrations and land use and land-cover change. A total of 14 models drawn from MsTMIP version 1 (http://dx.doi.org/10.3334/ORNLDAAC/1225) were used in this study. These 14 MsTMIP models are: CLASS-CTEM-N, CLM4, CLM4VIC, DLEM, GTEC, ISAM, LPJ-wsl, ORCHIDEE-LSCE (ORCHIDEE hereafter), SiB3, SiBCASA, TEM6, TRIPLEX-GHG, VEGAS2.1 and VISIT (see Table S1). These 14 MsTMIP models were further divided into two groups based on whether an explicit nitrogen cycling was considered in each model, namely models with carbon-nitrogen interactions (CLASS-CTEM-N, CLM4, CLM4VIC, DLEM, ISAM, TEM6, TRIPLEX-GHG) and models without carbon-nitrogen interactions (GTEC, LPJ-wsl, ORCHIDEE, SiB3, SiBCASA, VEGAS2.1, VISIT). The direct outputs of four variables, i.e., the net primary productivity (NPP), GPP, total living biomass and total soil carbon were used in this analysis. Because TRIPLEX-GHG did not report the NPP data in the MsTMIP variable list online, we calculated NPP indirectly as the difference between GPP and autotrophic respiration (Ra) for this model. To estimate the changes in modeled carbon storage, three models were excluded for reasons being either unavailability of carbon pools (SiB3), or nonphysical values of soil carbon stocks (CLASS-CTEM-N and CLM4VIC) due to

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210
issue such as non-linear interactions among hydrologic and carbon cycle parameters (Sargsyan et al., 2014; Huang et al., 2016). In MsTMIP, NPP equals BP as all carbon cycle models do not deal with non-structural NPP components (e.g. root exudation and root symbionts) explicitly at present. Thus the modeled NPP is hereafter called modeled BP in our study.

All of the global land cover maps used in MsTMIP models originated from the Synergetic Land Cover Product (SYNMAP) (Jung et al., 2006; Wei et al., 2014). To be compatible with the five vegetation types adopted by this study, the 47 SYNMAP classes (excluding shrubs, snow, barren and urban areas) used in SYNMAP were first reclassified into three plant functional types (PFTs) (i.e., forests, grasslands and croplands) based on a translation method defined by Jung et al. (2006). Here, croplands included all those SYNMAP vegetation classes with “crops” as a life form. According to the Köppen-Geiger climate classification (Peel et al., 2007), the forests located in boreal/tropical regions were then merged into boreal/tropical forests (BoF/TrF), while those over temperate areas were renamed as temperate forests (TeF). Using this reconstructed land cover map, represented PFT was defined as the PFT covering more than 50% of the fractional area in each grid cell. To obtain global natural vegetation BP and GPP, the grid cells with cropland as represented PFT were excluded in the calculation.

## Model constraining

We further adopted the recently widely used “emergent constraint” approach (e.g. Cox et al., 2013; Wenzel et al., 2014; Zhao et al., 2016; Kwiatkowski et al., 2017) to reduce the uncertainties in MsTMIP-simulated global mean BPE. This approach relies on a strong linear inter-model relationship between one quantity that cannot be directly measured and the other
quantity with available field observations. A more accurate estimate can be achieved by
constraining the first quantity through available field observations of the second quantity
(illustrated schematically in the Box 1 from Eyring et al. (2019)). Here, using a similar
method like Zhao et al. (2016), we first built a linear relationship between global and
site-level mean BPE values across the 14 MsTMIP models (as shown in Fig. 1a). This
model-derived relationship between global and site-level mean BPE, combined with field
observations, provided an emergent constraint on the simulated global mean BPE. Specifically,
the conditional probability density function (PDF) for the constrained global BPE was
calculated by integrating the PDF of the observations with the PDF of the model-derived
regression line. More detailed information of the “emergent constraint” approach can be
found in Appendix S1.

247

We first computed the values of site-level BPE from modeled BP and GPP at the same
location and for the same vegetation type than observed data. We only used BPE data from 83
unmanaged sites (except for cropland) because forest and grassland management was not
taken into account in the MsTMIP models. As cropland management is not well represented
in current models with simplified and incomplete description (Piao et al., 2018), the cropland
sites are not included in the emergent constraint for the global mean BPE, while they are
included for the vegetation type level BPE only for comparison. We also excluded 28 sites
because their vegetation type was a fraction smaller than 50% of the modeled vegetation in
the corresponding grid-cell (hereafter these data are called non-representative data). This
leaves 55 sites for our analysis, including 36 non-cropland sites. Hence, the constraining
approach was applied to all PFTs combined (Fig. 1, excluding cropland, n=36) and to each
individual PFT (Fig. 2, Figs S3 and S4).

248

249

250

251

252

253

254

255

256

257

258

259

260
261 Changes in global ecosystem carbon storage during 1901-2010

262 We used a total of 11 MsTMIP models (i.e., CLM4, DLEM, GTEC, ISAM, LPJ-wsl, 263 ORCHIDEE, SiBCASA, TEM6, TRIPLEX-GHG, VEGAS2.1, and VISIT) ‘SG3’ simulations 264 with GPP, BP, soil carbon storage, and biomass carbon storage available. Based on the 265 constrained global BPE, we further re-estimated the changes in global (excluding cropland) 266 ecosystem carbon storage during 1901-2010 for each model as following steps. 267

268 First, we calculated the turnover times of global ecosystem biomass carbon and soil 269 carbon in each year during 1902-2010 for each model using the following equation set.

$$
\begin{equation}
\begin{cases}
\tau_B(i) = \frac{c_B(i)}{\text{NPP}(i) - \Delta c_B(i)} \\
k_B(i) = \frac{1}{\tau_B(i)} \\
\tau_S(i) = \frac{c_S(i)}{c_B(i) \times k_B(i) - \Delta c_S(i)} \\
k_S(i) = \frac{1}{\tau_S(i)}
\end{cases}
\tag{1}
\end{equation}
$$

271 where,

272 C<sub>B</sub>(i) is the biomass carbon storage in year i,

273
274 ΔC<sub>B</sub>(i) is the change in biomass carbon storage from year i-1 to year i [ΔC<sub>B</sub>(i)=C<sub>B</sub>(i)-
275 C<sub>B</sub>(i-1)],

276 NPP(i) is the NPP in year i,

277 $\tau_B(i)$ is the turnover time of biomass carbon in year i,

278 $k_B(i)$ is the turnover rate of biomass carbon in year i,

279 C<sub>S</sub>(i) is the soil carbon storage in year i,

280 ΔC<sub>S</sub>(i) is the change in soil carbon storage from year i-1 to year i [ΔC<sub>S</sub>(i)= C<sub>S</sub>(i)-C<sub>S</sub>(i-1)],
281

282 $\tau_S(i)$ is the turnover time of soil carbon in year i,
$k_S(i)$ is the turnover rate of soil carbon in year i,

i is year during 1902-2010 (i = 2, 3, 4, ..., 110).

285

Second, assuming that simulated carbon pools in 1901 are in equilibrium (Fig. S6) and BPE is invariant during 1901-2010 (equals to, in numeral, the constrained BPE). The biomass ($C'_B$) and soil carbon ($C'_S$) storage in each model in 1901 were recalculated as $C'_B(1) = GPP(1) \times BPE_{Constrained} \times \tau_B(1)$ and $C'_S(1) = GPP(1) \times BPE_{Constrained} \times \tau_S(1)$, respectively. Here, $C'_B(1)$ and $C'_S(1)$ is the re-estimated biomass and soil carbon storage in each model in 1901, respectively. GPP(1) is the GPP in the first year (1901). $BPE_{Constrained}$ is the constrained global natural vegetation BPE (0.41). $\tau_B(1)$ is the turnover time of biomass carbon in the first year (1901), assuming that $\tau_B(1)$ equals $\tau_B(2)$ in equation (1). $\tau_S(1)$ is the turnover time of soil carbon in the first year (1901), assuming that $\tau_S(1)$ equals $\tau_S(2)$ in equation (1). Then, we reconstructed the global biomass carbon and soil carbon in each year during 1902-2010 for each model using the following equations.

286

287

288

289

290

291

292

293

294

295

296

297

$$ \begin{cases} C'_{B}(i) - C'_{B}(i-1) = GPP(i) \times BPE_{Constrained} - C'_{B}(i) \times k_{B}(i) \\ C'_{S}(i) - C'_{S}(i-1) = C'_{B}(i) \times k_{B}(i) - C'_{S}(i) \times k_{S}(i) \end{cases} \quad (2) $$

298

299

300 where,

301 $C'_B(i)$ is the re-estimated biomass carbon storage in year i,

302 GPP(i) is the GPP in year i,

303 $BPE_{Constrained}$ is the constrained global natural vegetation BPE (0.41),

304 $k_B(i)$ is the turnover rate of biomass carbon in year i,

305 $C'_S(i)$ is the re-estimated soil carbon storage in year i,

306 $k_S(i)$ is the turnover rate of soil carbon in year i,

307 i is year during 1902-2010 (i = 2, 3, 4, ..., 110).
308

309

310

311

312

313

314

315

# Results

316

## Emergent constraint on the global mean BPE

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

Fig. 1a shows the relationship between modeled global and site-level mean BPE (cropland is excluded throughout this section, see Methods). MsTMIP models simulate a large spread in global non-cropland BPE, with a range going from 0.36 for the TRIPLEX-GHG model to 0.62 for the SiB3 model (Fig. 1a). We found a strong linear relationship across different models between their site-level and global mean BPE ($R^2$ = 0.77, *P* < 0.05). This allows us to define an emergent constraint relationship (see Methods; Cox *et al.*, 2013; Wenzel *et al.*, 2014; Zhao *et al.*, 2016; Kwiatkowski *et al.*, 2017) whereby the selected field observations can be used to constrain the global modeled value of BPE. According to the emerging relationship shown in Fig. 1a and the average of selected field BPE observations (0.44), the constrained global best estimate of BPE is 0.41. It is lower than the original multi-model average (0.46). When compared to the unconstrained probability density function (PDF) of models (black line), the constrained PDF (red line) is more narrow, so that the uncertainty on BPE is reduced from ± 0.08 to ± 0.05 (1-sigma standard deviation, Fig. 1b).

Given that several studies showed nutrient scarcity to be associated with lower BPE for comparable GPP level (Vicca *et al.*, 2012; Fernández-Matínez *et al.*, 2015), we divided the
models into two groups, those with (n = 7) and without (n = 7) carbon-nitrogen interactions, to test possible controls of BPE by nutrient limitations. Result shows that the constrained global BPE and uncertainty range are almost the same between the two groups of models (0.40 ± 0.05 for models with carbon-nitrogen interactions and 0.41 ± 0.05 for the other models; Fig S1). Models with carbon-nitrogen interactions (0.47 ± 0.08) tend to better match site observations (0.44 ± 0.04) compared to carbon-only models (0.53 ± 0.09), because of their lower BPE (Fig. S1). At the global scale, the BPE from models with carbon-nitrogen interactions (0.43 ± 0.07) is also smaller than that from carbon-only models (0.49 ± 0.08) (Fig. S1). This is consistent with a decrease of BPE observed at forest sites under nitrogen limitation (Vicca *et al.*, 2012). Furthermore, the constrained global mean BPE for models with carbon-nitrogen interaction is robust when separately calculated using site-level BPE with different nutrient availability classes (categorized as low, medium and high nutrient level; see Methods) (Fig. S2).

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

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Emergent constraint on the global mean BPE. (a) Relationship between the modeled global (excluding cropland) and site-level mean BPE of unmanaged vegetation. The black line represents the best fit across different models (black dots), the vertical red dashed</figcaption></figure>
line and shaded area represent the mean value of BPE and its statistical uncertainty from field
observations, respectively (see Methods), and the horizontal red dashed line represents the
constrained global mean BPE. (b) Probability density function for global BPE. The black line
represents the probability distribution of the multi-model mean BPE, assuming that all BPE
values estimated from terrestrial carbon models conform to a Gaussian distribution, the red
line represents the probability distribution of global BPE constrained by field observations,
and the black and red letters represent the values of global BPE (mean ± standard deviation)
before and after the field constraints, respectively.

359

360

Emergent constraint on the mean BPE for different vegetation types

361
Applying a similar emerging constraint method, we also constrained BPE
362
for the different vegetation types, including forests, grasslands and croplands (see Methods).
363
To do so, we conventionally divided forest sites into boreal, temperate and tropical groups
364
(Luyssaert et al., 2007; Campioli et al., 2015). Strong linear relationships like the one in Fig.
365
1 were also found within each vegetation type (P < 0.05, Fig. 2) particularly for boreal and
366
tropical forests, with R² > 0.9. The constrained BPE is substantially lower than the original
367
multi-model average for each vegetation type, except cropland (Fig. 2). A noticeable decrease
368
of BPE (constrained value minus original model mean) is found in boreal forests (Fig. 2a).
369
The multi-model mean BPE for boreal forests is 0.55 ± 0.11, but the magnitude of BPE after
370
the constraint from observations is decreased to 0.48 ± 0.06. In contrast, crops have a higher
371
BPE when constrained by field observations (Fig. 2e). The higher BPE for crops from
372
observations (0.55 ± 0.06) than from models (0.50 ± 0.10) may be related to several traits and
373
processes that are not included in models, including crop varieties and more generally the
374
allocation processes specific to crop plant types, and management operations such as
375
irrigation and fertilization. The sample sizes of the site data (n) used to constrain BPE in some
vegetation types are small (n=5), but the constrained values of BPE are robust even if we
include non-representative data (Figs S3), or if we base the site level model estimates on all
pixels with similar climate as the site observations instead of only the pixel of the actual site
location (the difference in climate conditions within an interval of 1°C × 50 mm) (Figs S4).
Overall, among the natural vegetation types, the lowest value of BPE after applying the
constraint is found for tropical forests (0.40 ± 0.04) and the largest is for boreal forests (0.48 ±
0.06) (Table S2).
<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" /></figure>

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

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

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

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

<figure><img src="image_14.png" /></figure>
Figure 1 but the field observations are divided into five vegetation types: (a, b) Boreal forests (BoF). (c, d) Temperate forests (TeF). (e, f) Tropical forests (TrF). (g, h) Grassland. (i, j) Cropland. Only unmanaged observations (except for cropland) for which the fraction of their represented vegetation type is higher than 50% in area in the corresponding model pixel are included. Right-hand column shows the relationship between modeled biome-specific and site-level mean BPE across 14 MsTMIP models, for 5 vegetation types respectively. In each panel, n is the sample size of the field observations, the black line represents the best fit across different models (black dots), the vertical red dashed line and shaded area represent the mean value of BPE and its statistical uncertainty from field observations, respectively, and the horizontal red dashed line represents the constrained mean BPE. Left-hand column shows the probability distributions of mean BPE across 14 MsTMIP models, for 5 vegetation types respectively. In each panel, the black line represents the probability distribution of multi-model mean BPE, assuming that all BPE values estimated from terrestrial carbon models can be represented by a Gaussian distribution, the red line represents the probability distribution of global mean BPE constrained by field observations, and the black and red letters represent the values of BPE (mean ± standard deviation) before and after constraint, respectively.

## Discussion

### Possible reasons for overestimated BPE in models

As mentioned above, BPE can be used as a critical diagnostic parameter to characterize how carbon is partitioned among different pathways, providing us an opportunity for testing model performance. Our study suggested an overestimation of BPE by most models in natural ecosystems. That is to say, in current models, the fraction of carbon sequestrated by plants allocated to BP is overestimated, but less for respiration purpose. As current models have
similar structures in representing photosynthesis, they differ greatly in simulating autotrophic
respiration ($R_a$). $R_a$ is usually separated into two components: growth ($R_g$) and maintenance
respiration ($R_m$) (Ryan 1995; Dewar et al., 1999; Thornley & Cannell, 2000; Gifford, 2003;
Atkin et al., 2014). This concept has been widely utilized in carbon cycle models, where $R_m$ is
first predicted as a function of temperature and biomass, and $R_g$ is a fixed proportion of
productivity (Cox, 2001; Piao et al., 2010). The simulations and parameterizations of $R_m$ and
$R_g$ vary greatly among models, for instance, constant growth respiration coefficient in CLM is
0.3 (Lawrence et al., 2011), but in TEM6 is 0.2 (Hayes et al., 2011), potentially resulting in
large biases of BPE outputs.

420

On the other hand, whether a model has an explicit consideration of carbon-nitrogen
interaction appears to strongly influence modeled BPE values. Although carbon cycle models
overestimate plant BPE, models without carbon-nitrogen interactions show much more
model-data mismatch. These models, without a consideration of nitrogen limitation on carbon
cycle, may overestimate the CO₂ fertilization effect of vegetation growth and therefore
overestimate plant BPE (Vitousek & Howarth, 1991). As for models with carbon-nitrogen
interactions, although modeled BPE is smaller than that in carbon-only models and better
match with the benchmark (Fig. 3a), this does not imply that carbon-nitrogen models perform
better, because if modeled GPP and BP are both well simulated, $R_a$ in the model will have to
compensate for ignoring the production of non-structural compounds, and probably
overestimating $R_a$ in this case. On condition that the production of non-structural compounds
are incorporated in modeling processes, simulated $R_a$ will have to be (objectively or
inadvertently) tuned down. Thus, the processes related to BPE still need improvements in
carbon-nitrogen models. This BPE case study also warns that future model development
should avoid parameters over-tuning for “perfect” model performance when there are key
436

processes missing yet in models, to generate the right causality in carbon cycle processes and
vegetation dynamics.

437

438

439

Implications of overestimated BPE in models

440

As shown in Fig. 3a, global BP (excluding cropland) ranges from 35 Pg C yr⁻¹ in CLM4VIC to 88 Pg C yr⁻¹ in SiB3. Since most models overestimate BPE, they most likely also overestimate BP (in case of similar or overestimated global GPP in models). To test this hypothesis, we quantify global natural vegetation BP (excluding cropland) by multiplying the constrained estimate of BPE presented above by observation-based GPP estimates from three different products: the Moderate Resolution Imaging Spectroradiometer (MODIS) GPP product (Zhao *et al.*, 2005), a eddy covariance flux tower data-driven GPP product (Jung *et al.*, 2011) and a new GPP product derived by Breathing Earth System Simulator (Jiang & Ryu, 2016) — all averaged over non-cropland vegetation only. The result shows a global BP (excluding cropland) of about 41 ± 6 Pg C yr⁻¹ which is lower by 24 ± 11% than that originally simulated by the models (10 out of 14 models) (Fig. 3).

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

The overestimation of BPE and BP by carbon cycle models implies reduced carbon
storage change during the historical period in models, because carbon storage change is
controlled by BP changes and by the turnover rate of ecosystem carbon pools (Parton et al.,
1987; Potter et al., 1993; Luo et al., 2003; Xia et al., 2013). By correcting the MsTMIP model
outputs for their bias of BPE based on equations for substitute models (see Methods) and not
accounting for BPE change in response to climate and rising CO₂ (see Fig. S5), we found that
the lower BPE implies a 67% smaller increase of ecosystem carbon storage during 1901-2010
(from 86 ± 144 Pg C of increase in average ecosystem carbon storage in original model
simulations to 28 ± 115 Pg C of increase in average ecosystem carbon storage, *P* < 0.01, Fig.
3b), assuming unchanged carbon turnover rate. This overestimation of ecosystem carbon storage increase due to the bias of BPE is particularly significant for carbon-only models ($P < 0.01$, Fig. 3b). Compared to carbon-only models, models with carbon-nitrogen generally reproduce smaller ecosystem carbon storage changes with original model simulations ($P = 0.09$, Fig. 3b), because their BPE is more comparable with the constrained BPE.

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

<figure><img src="image_16.png" /><figcaption>**Figure 3.** Global BP and change in global ecosystem carbon storage during 1901-2010 (excluding cropland). (a) Estimates of GPP, BPE and BP expressed as the dots in the space,</figcaption></figure>
with each represents a model. The black star represents the reference values (i.e., the
observation-based estimate of GPP, data-constrained BPE and newly-estimated global BP).
Shaded gray areas are the uncertainty ranges of the referenced GPP and BP. The uncertainty
bars for BP and GPP for the models represent the standard deviations for 1981-2010. (b)
Change in global terrestrial carbon storage from 1901 to 2010 estimated by original MsTMIP
model simulations and by collecting the mean bias of BPE derived from the MsTMIP models
using constrained BPE. The black triangle represents the mean value of all models, and the
orange and blue triangles represent the mean values of C-N coupled and C-only models,
respectively. The uncertainty bars for the triangles are the standard deviations between models.
The solid and dashed bars indicate significant ($P < 0.05$) and marginally significant ($0.05 < P < 0.10$) differences, respectively, between the original and the reproduced results. Orange and
blue dots in both panels represent C-N coupled and C-only models, respectively.

In summary, our results show that carbon cycle models generally overestimate BPE, not
only suggesting that the capability of future terrestrial carbon sequestration is less than current
model predictions, but also implying that future reduction of greenhouse gas emission to meet
climate target may be underestimated. Furthermore, models with carbon-nitrogen interactions
simulate BPE closer to observed values compared to carbon-only models, although Ra in
these models are probably overestimated to compensate for the omission of non-structural
carbon fluxes. These results suggest that more efforts should be put into improving carbon
cycle modeling with processes treated either very simplistically or insufficiently, in order to
represent correct carbon cycle processes and vegetation dynamics for correct reasons.
Measurements of BPE across different vegetation types, for different management intensities
and different levels of soil-available nitrogen and phosphorus are also needed to refine the
global estimation of this key carbon cycle variable.
495

496

References

497

Atkin, O. K., Meir, P., & Turnbull, M. H. (2014). Improving representation of leaf respiration in large-scale predictive climate-vegetation models. *New Phytologist*, 202(3), 743-748.

498

499

Badgley, G., Field, C. B., & Berry, J. A. (2017). Canopy near-infrared reflectance and terrestrial photosynthesis. *Science Advances*, 3(3), e1602244.

500

501

Beer, C., Reichstein, M., Tomelleri, E., Ciais, P., Jung, M., Carvalhais, N., ... & Bondeau, A. (2010). Terrestrial gross carbon dioxide uptake: global distribution and covariation with climate. *Science*, 329(5993), 834-838.

502

503

504

Campioli, M., Vicca, S., Luyssaert, S., Bilcke, J., Ceschia, E., Chapin Iii, F. S., ... & Olefeldt, D. (2015). Biomass production efficiency controlled by management in temperate and boreal ecosystems. *Nature Geoscience*, 8(11), 843.

505

506

507

Chapin, F. S., Woodwell, G. M., Randerson, J. T., Rastetter, E. B., Lovett, G. M., Baldocchi, D. D., ... & Wirth, C. (2006). Reconciling carbon-cycle concepts, terminology, and methods. *Ecosystems*, 9(7), 1041-1050.
Ciais, P., Sabine, C., Bala, G., Bopp, L., Brovkin, V., Canadell, J., ... & Jones, C. (2013). Carbon and other biogeochemical cycles. *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 University Press Cambridge United Kingdom and New York NY USA, 465-570.

508

509

510

511

512

513

514

Cox, P. M. (2001). Description of the “TRIFFID” dynamic global vegetation model. *Hadley Centre technical note*, 24, 1-16, 2001.

515

516

Cox, P. M., Pearson, D., Booth, B. B., Friedlingstein, P., Huntingford, C., Jones, C. D., & Luke, C. M. (2013). Sensitivity of tropical carbon to climate change constrained by carbon dioxide variability. *Nature*, 494(7437), 341.

517

518

519

Cramer, W., Bondeau, A., Woodward, F. I., Prentice, I. C., Betts, R. A., Brovkin, V., ... &
520

Kucharik, C. (2001). Global response of terrestrial ecosystem structure and function to CO₂ and climate change: results from six dynamic global vegetation models. Global change biology, 7(4), 357-373.

521

522

523

De Kauwe, M. G., Medlyn, B. E., Zaehle, S., Walker, A. P., Dietze, M. C., Wang, Y. P., ... & Wårlind, D. (2014). 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 Phytologist, 203(3), 883-899.

524

525

526

527

Dewar, R. C., Medlyn, B. E., & Mcmurtrie, R. E. (1999). Acclimation of the respiration/photosynthesis ratio to temperature: insights from a model. Global Change Biology, 5(5), 615-622.

528

529

530

Erb, K. H., Fetzel, T., Plutzar, C., Kastner, T., Lauk, C., Mayer, A., ... & Haberl, H. (2016). Biomass turnover time in terrestrial ecosystems halved by land use. Nature Geoscience, 9(9), 674.

531

532

533

Eyring, V., Cox, P. M., Flato, G. M., Gleckler, P. J., Abramowitz, G., Caldwell, P., ... & Hurtt, G. C. (2019). Taking climate model evaluation to the next level. Nature Climate Change, 9, 102-110.

534

535

536

Fernández-Martínez, M., Vicca, S., Janssens, I. A., Sardans, J., Luyssaert, S., Campioli, M., ... & Papale, D. (2014). Nutrient availability as the key regulator of global forest carbon balance. Nature Climate Change, 4(6), 471.

537

538

539

Fick, S. E., & Hijmans, R. J. (2017). WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International journal of climatology, 37(12), 4302-4315.

540

541

Friedlingstein, P., Cox, P., Betts, R., Bopp, L., von Bloh, W., Brovkin, V., ... & Bala, G. (2006). Climate-carbon cycle feedback analysis: results from the C4MIP model intercomparison. Journal of climate, 19(14), 3337-3353.

542

543

544

Friend, A. D., Lucht, W., Rademacher, T. T., Keribin, R., Betts, R., Cadule, P., ... & Ito, A.
545

(2014). Carbon residence time dominates uncertainty in terrestrial vegetation responses to future climate and atmospheric CO₂. Proceedings of the National Academy of Sciences, 111(9), 3280-3285.

546

547

548

Gifford, R. M. (1995). Whole plant respiration and photosynthesis of wheat under increased CO₂ concentration and temperature: long-term vs. short-term distinctions for modelling. Global Change Biology, 1(6), 385-396.

549

550

551

Gifford, R. M. (2003). Plant respiration in productivity models: conceptualisation, representation and issues for global terrestrial carbon-cycle research. Functional Plant Biology, 30(2), 171-186.

552

553

554

Guanter, L., Zhang, Y., Jung, M., Joiner, J., Voigt, M., Berry, J. A., ... & Moran, M. S. (2014). Global and time-resolved monitoring of crop photosynthesis with chlorophyll fluorescence. Proceedings of the National Academy of Sciences, 111(14), E1327-E1333.

555

556

557

Haberl, H., Erb, K. H., Krausmann, F., Gaube, V., Bondeau, A., Plutzar, C., ... & Fischer-Kowalski, M. (2007). Quantifying and mapping the human appropriation of net primary production in earth's terrestrial ecosystems. Proceedings of the National Academy of Sciences, 104(31), 12942-12947.

558

559

560

561

Hayes, D. J., McGuire, A. D., Kicklighter, D. W., Gurney, K. R., Burnside, T. J., & Melillo, J. M. (2011). Is the northern high-latitude land-based CO₂ sink weakening? Global Biogeochemical Cycles, 25(3).

562

563

564

Huang, M., Ray, J., Hou, Z., Ren, H., Liu, Y., & Swiler, L. (2016). On the applicability of surrogate-based Markov chain Monte Carlo-Bayesian inversion to the Community Land Model: Case studies at flux tower sites. Journal of Geophysical Research: Atmospheres, 121(13), 7548-7563.

565

566

567

568

Huntzinger, D. N., Schwalm, C., Michalak, A. M., Schaefer, K., King, A. W., Wei, Y., ... & Berthier, G. (2013). The north american carbon program multi-scale synthesis and

569
570

terrestrial model intercomparison project-part 1: Overview and experimental
design. Geoscientific Model Development, 6(6), 2121-2133.

571

572

Jiang, C., & Ryu, Y. (2016). Multi-scale evaluation of global gross primary productivity and evapotranspiration products derived from Breathing Earth System Simulator (BESS). Remote Sensing of Environment, 186, 528-547.

573

574

575

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

576

577

578

Jung, M., Reichstein, M., Margolis, H. A., Cescatti, A., Richardson, A. D., Arain, M. A., ... & Gianelle, D. (2011). Global patterns of land-atmosphere fluxes of carbon dioxide, latent heat, and sensible heat derived from eddy covariance, satellite, and meteorological observations. Journal of Geophysical Research: Biogeosciences, 116(G3).

579

580

581

582

Krinner, G., Viovy, N., de Noblet - Ducoudré, N., Ogée, J., Polcher, J., Friedlingstein, P., ... & Prentice, I. C. (2005). A dynamic global vegetation model for studies of the coupled atmosphere - biosphere system. Global Biogeochemical Cycles, 19(1).

583

584

585

Kwiatkowski, L., Bopp, L., Aumont, O., Ciais, P., Cox, P. M., Laufkötter, C., ... & Séférian, R. (2017). Emergent constraints on projections of declining primary production in the tropical oceans. Nature Climate Change, 7(5), 355.

586

587

588

Landsberg, J. J., & Waring, R. H. (1997). A generalised model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning. Forest ecology and management, 95(3), 209-228.

589

590

591

Lawrence, D. M., Oleson, K. W., Flanner, M. G., Thornton, P. E., Swenson, S. C., Lawrence, P. J., ... & Bonan, G. B. (2011). Parameterization improvements and functional and structural advances in version 4 of the Community Land Model. Journal of Advances in Modeling Earth Systems, 3(1).

592

593

594
595 Le Quéré, C., Moriarty, R., Andrew, R. M., Canadell, J. G., Sitch, S., Korsbakken, J. I., ... & 596 Houghton, R. A. (2015). Global carbon budget 2015. Earth System Science Data, 7, 597 349-396.

598 Luo, Y., White, L. W., Canadell, J. G., DeLucia, E. H., Ellsworth, D. S., Finzi, A., ... & 599 Schlesinger, W. H. (2003). Sustainability of terrestrial carbon sequestration: A case study 600 in Duke Forest with inversion approach. Global biogeochemical cycles, 17(1).

601 Luo, Y., Ahlström, A., Allison, S. D., Batjes, N. H., Brovkin, V., Carvalhais, N., ... & 602 Georgiou, K. (2016). Toward more realistic projections of soil carbon dynamics by Earth 603 system models. Global Biogeochemical Cycles, 30(1), 40-56.

604 Luyssaert, S., Inglima, I., Jung, M., Richardson, A. D., Reichstein, M., Papale, D., ... & 605 Aragao, L. E. O. C. (2007). CO₂ balance of boreal, temperate, and tropical forests 606 derived from a global database. Global change biology, 13(12), 2509-2537.

607 McCree, K. J., & Troughton, J. H. (1966). Prediction of growth rate at different light levels 608 from measured photosynthesis and respiration rates. Plant Physiology, 41(4), 559-566.

609 Nemani, R., Hashimoto, H., Votava, P., Melton, F., Wang, W., Michaelis, A., ... & White, M. 610 (2009). Monitoring and forecasting ecosystem dynamics using the Terrestrial 611 Observation and Prediction System (TOPS). Remote Sensing of Environment, 113(7), 612 1497-1509.

613 Parton, W. J., Schimel, D. S., Cole, C. V., & Ojima, D. S. (1987). Analysis of factors 614 controlling soil organic matter levels in Great Plains Grasslands. Soil Science Society of 615 America Journal, 51(5), 1173-1179.

616 Peel, M. C., Finlayson, B. L., & McMahon, T. A. (2007). Updated world map of the 617 Köppen-Geiger climate classification. Hydrology and earth system sciences 618 discussions, 4(2), 439-473.

619 Peng, C., Liu, J., Dang, Q., Apps, M. J., & Jiang, H. (2002). TRIPLEX: a generic hybrid
model for predicting forest growth and carbon and nitrogen dynamics. Ecological modelling, 153(1-2), 109-130.

Peñuelas, J., Ciais, P., Canadell, J. G., Janssens, I. A., Fernández-Martínez, M., Carnicer, J., ... & Sardans, J. (2017). Shifting from a fertilization-dominated to a warming-dominated period. Nature ecology & evolution, 1(10), 1438.

Piao, S., Ciais, P., Friedlingstein, P., Peylin, P., Reichstein, M., Luyssaert, S., ... & Grelle, A. (2008). Net carbon dioxide losses of northern ecosystems in response to autumn warming. Nature, 451(7174), 49.

Piao, S., Luyssaert, S., Ciais, P., Janssens, I. A., Chen, A., Cao, C., ... & Wang, S. (2010). Forest annual carbon cost: A global-scale analysis of autotrophic respiration. Ecology, 91(3), 652-661.

Piao, S., Sitch, S., Ciais, P., Friedlingstein, P., Peylin, P., Wang, X., ... & Huntingford, C. (2013). Evaluation of terrestrial carbon cycle models for their response to climate variability and to CO₂ trends. Global change biology, 19(7), 2117-2132.

Piao, S., Huang, M., Liu, Z., Wang, X., Ciais, P., Canadell, J. G., ... & Le Quéré, C. (2018). Lower land-use emissions responsible for increased net land carbon sink during the slow warming period. Nature Geoscience, 11(10), 739.

Potter, C. S., Randerson, J. T., Field, C. B., Matson, P. A., Vitousek, P. M., Mooney, H. A., & Klooster, S. A. (1993). Terrestrial ecosystem production: a process model based on global satellite and surface data. Global Biogeochemical Cycles, 7(4), 811-841.

Rogers, A., Medlyn, B. E., Dukes, J. S., Bonan, G., Von Caemmerer, S., Dietze, M. C., ... & Prentice, I. C. (2017). A roadmap for improving the representation of photosynthesis in Earth system models. New Phytologist, 213(1), 22-42.

Ryan, M. G., Gower, S. T., Hubbard, R. M., Waring, R. H., Gholz, H. L., Cropper, W. P., & Running, S. W. (1995). Woody tissue maintenance respiration of four conifers in
645 contrasting climates. Oecologia, 101(2), 133-140.

646 Sands, P. J., Battaglia, M., & Mummery, D. (2000). Application of process-based models to
647 forest management: experience with PROMOD, a simple plantation productivity model.

648 Tree Physiology, 20(5-6), 383-392.

649 Sargsyan, K., Safta, C., Najm, H. N., Debusschere, B. J., Ricciuto, D., & Thornton, P. (2014).
650 Dimensionality reduction for complex models via Bayesian compressive

651 sensing. International Journal for Uncertainty Quantification, 4(1).

652 Tian, H., Chen, G., Zhang, C., Liu, M., Sun, G., Chappelka, A., ... & Chen, H. (2012).
653 Century-scale responses of ecosystem carbon storage and flux to multiple environmental

654 changes in the southern United States. Ecosystems, 15(4), 674-694.

655 Veroustraete, F., Sabbe, H., & Eerens, H. (2002). Estimation of carbon mass fluxes over
656 Europe using the C-Fix model and Euroflux data. Remote Sensing of Environment, 83(3),

657 376-399.

658 Vicca, S., Luyssaert, S., Penuelas, J., Campioli, M., Chapin III, F. S., Ciais, P., ... & Malhi, Y.
659 (2012). Fertile forests produce biomass more efficiently. Ecology letters, 15(6), 520-526.

660 Vitousek, P. M., & Howarth, R. W. (1991). Nitrogen limitation on land and in the sea: how
661 can it occur? Biogeochemistry, 13(2), 87-115.

662 Waring, R. H., Landsberg, J. J., & Williams, M. (1998). Net primary production of forests: a
663 constant fraction of gross primary production? Tree physiology, 18(2), 129-134.

664 Wei, Y., Liu, S., Huntzinger, D. N., Michalak, A. M., Viovy, N., Post, W. M., ... & Tian, H.
665 (2014). The North American carbon program multi-scale synthesis and terrestrial model

666 intercomparison project-Part 2: environmental driver data. Geoscientific Model
667 Development Discussions, 7, 2875-2893.

668 Wenzel, S., Cox, P. M., Eyring, V., & Friedlingstein, P. (2014). Emergent constraints on
669 climate-carbon cycle feedbacks in the CMIP5 Earth system models. Journal of
670 Geophysical Research: Biogeosciences, 119(5), 794-807.

671 Xia, J., Luo, Y., Wang, Y. P., & Hararuk, O. (2013). Traceable components of terrestrial carbon storage capacity in biogeochemical models. *Global Change Biology*, 19(7), 2104-2116.

674 Yao, Y., Piao, S., & Wang, T. (2018). Future biomass carbon sequestration capacity of Chinese forests. *Science Bulletin*, 63(17), 1108-1117.

676 Zhang, Y., Xu, M., Chen, H., & Adams, J. (2009). Global pattern of NPP to GPP ratio derived from MODIS data: effects of ecosystem type, geographical location and climate. *Global Ecology and Biogeography*, 18(3), 280-290.

679 Zhang, Y., Yu, G., Yang, J., Wimberly, M. C., Zhang, X., Tao, J., ... & Zhu, J. (2014). Climate-driven global changes in carbon use efficiency. *Global Ecology and Biogeography*, 23(2), 144-155.

682 Zhao, C., Piao, S., Wang, X., Huang, Y., Ciais, P., Elliott, J., ... & Liu, Y. (2017). Plausible rice yield losses under future climate warming. *Nature plants*, 3(1), 16202.

684 Zhao, M., Heinsch, F. A., Nemani, R. R., & Running, S. W. (2005). Improvements of the MODIS terrestrial gross and net primary production global data set. *Remote sensing of Environment*, 95(2), 164-176.

687 Zhu, Q., Liu, J., Peng, C., Chen, H., Fang, X., Jiang, H., ... & Zhou, X. (2014). Modelling methane emissions from natural wetlands by development and application of the TRIPLEX-GHG model. *Geoscientific Model Development*, 7(3), 981-999.

690

691 Acknowledgements

692 This study was supported by the National Natural Science Foundation of China (grant numbers 41722101, 41701089 and 41671079), and National Key Research and Development Program of China (2016YFC0500203 and 2016YFA0600202). P.C., I.J. and J.P.
695 acknowledges support from the European Research Council Synergy grant
696 ERC-2013-SyG-610028 IMBALANCE-P. We thank Multi-scale Synthesis and Terrestrial
697 Model Intercomparison Project (MsTMIP) providing simulation results.

698

699 **Conflict of Interest** The authors declare no competing financial interests.

700
701 **Supporting Information Captions**

702 **Tables**

703 **Table S1.** Information for the 14 terrestrial biosphere models from MsTMIP used in this study.
704

705 **Table S2.** Multi-model mean BPE and constrained BPE for non-cropland vegetation types.
706 Only observations from unmanaged sites were used in the constraint approach. 14 MsTMIP
707 models are classified based on whether they have an explicit consideration for
708 carbon-nitrogen (C-N) interactions.
709

710 **Figures**

711 **Figure S1.** Emergent constraint on the global mean BPE for two groups of the 14 models.
712 Same as Figure 1 in the main text but the models are classified based on whether they have an
713 explicit consideration for nitrogen cycling. (a, b) Models with carbon-nitrogen interactions
714 (i.e. TRIPLEX-GHG, TEM6, ISAM, DLEM, CLM4, CLM4VIC, and CLASS-CTEM-N). (c,
715 d) Models without carbon-nitrogen interactions (i.e. VISIT, VEGAS2.1, SiBCASA, SiB3,
716 ORCHIDEE, LPJ-wsl, and GTEC).

717 **Figure S2.** Emergent constraint on the global mean BPE for the models with carbon-nitrogen interactions. Same as Figure 1 in the main text but the field observations are divided into three classes based on site fertility: (a, b) Low nutrient availability. (c, d) Medium nutrient availability. (e, f) High nutrient availability.

721 **Figure S3.** Emergent constraint on the mean BPE for different vegetation types. Same as Figure 2 in the main text but non-representative data (observations for which the fraction of their vegetation type in the corresponding model pixel is higher than 50%) are also included.

724 **Figure S4.** Emergent constraint on the mean BPE for different vegetation types. Same as Figure 2 in the main text but only unmanaged observations (except for cropland) with similar
726

climate to the model pixels and for which the fraction of their vegetation type in the
727 corresponding pixel is higher than 50% are included.

728

**Figure S5.** Comparison between global mean BPE during 1901-1930 and 1981-2010 in the MsTMIP models. Orange and blue dots represent C-N coupled and C-only models, respectively. The black triangle represents the mean value of all models. The uncertainty bars are the standard deviations between models.

729

730

731

732

**Figure S6.** Comparison between changes in global ecosystem carbon storage during 1901-2010 from original model simulations, and based on reproduced initial carbon storage, assuming that the simulated carbon pools in 1901 are in equilibrium. Orange and blue dots represent C-N coupled and C-only models, respectively.

733

734

735