CO2 fertilization contributed more than half of the observed forest biomass increase in northern extra‐tropical land

The existence of a large‐biomass carbon (C) sink in Northern Hemisphere extra‐tropical ecosystems (NHee) is well‐established, but the relative contribution of different potential drivers remains highly uncertain. Here we isolated the historical role of carbon dioxide (CO2) fertilization by integrating estimates from 24 CO2‐enrichment experiments, an ensemble of 10 dynamic global vegetation models (DGVMs) and two observation‐based biomass datasets. Application of the emergent constraint technique revealed that DGVMs underestimated the historical response of plant biomass to increasing [CO2] in forests ( βForestMod ) but overestimated the response in grasslands ( βGrassMod ) since the 1850s. Combining the constrained βForestMod (0.86 ± 0.28 kg C m−2 [100 ppm]−1) with observed forest biomass changes derived from inventories and satellites, we identified that CO2 fertilization alone accounted for more than half (54 ± 18% and 64 ± 21%, respectively) of the increase in biomass C storage since the 1990s. Our results indicate that CO2 fertilization dominated the forest biomass C sink over the past decades, and provide an essential step toward better understanding the key role of forests in land‐based policies for mitigating climate change.


| INTRODUC TI ON
Terrestrial ecosystems have absorbed approximately 32% of the anthropogenic emissions of carbon dioxide (CO 2 ) since preindustrial times, broadly taken as the years 1850 to the present day (Friedlingstein et al., 2020).Such absorption has appreciably slowed the rate of global warming (Ballantyne et al., 2012;Schimel et al., 2015;Shevliakova et al., 2013).The substantial increase in the amount of biomass carbon (C) stored in Northern Hemisphere extratropical ecosystems (NHee) is thought to be an important part of the overall enhanced global terrestrial C sink (Liu et al., 2015;Pan et al., 2011;Tagesson et al., 2020;Xu et al., 2021).Such an increase in biomass is simultaneously regulated by multiple covarying factors, such as the CO 2 fertilization effect (CFE), climate change, nitrogen (N) deposition and land-use change, complicating the quantification of their individual effects (Fernández-Martínez et al., 2019;Peylin et al., 2013;Reich et al., 2020;Sitch et al., 2015).Previous studies have generally suggested that CO 2 fertilization drove the primary gain in stores of biomass C (Huntzinger et al., 2017;Piao et al., 2013;Schimel et al., 2015;Walker et al., 2021;Wenzel et al., 2016), but the direct impacts of climate change could ultimately suppress such gains in the future (Jiao et al., 2021;Peñuelas et al., 2017;Piao et al., 2020;Yuan et al., 2019), while limited N availability may eventually lower the ability to "draw down" atmospheric CO 2 concentrations ([CO 2 ]; Norby et al., 2010;Terrer et al., 2019;Zaehle et al., 2015).Isolating and accurately quantifying the magnitude of CO 2 fertilization and its contribution is accordingly essential, with implications for the formulation and implementation of land-based mitigation strategies that may maximize the offset of CO 2 emissions.
Considerable efforts have been made in recent decades to improve our knowledge of the CFE on the terrestrial C cycle, in particular by conducting CO 2 -enrichment experiments, many, but not all, of which report that elevated [CO 2 ] has stimulated plant growth and biomass C sequestration (McCarthy et al., 2010;Terrer et al., 2019;Walker et al., 2021).Site-scale measurements can provide critical local information on CFE, but they have limited ability to represent geographically extensive changes or long-term results due to their relative scarcity and short duration (~5 years on average; Terrer et al., 2021).As an alternative, dynamic global vegetation models (DGVMs) can provide complete geographical coverage and be configured to isolate and quantify long-term CFE around the globe, but require accurate parameterization.DGVMs collectively agree that CFE is the dominant cause of the increase in biomass, but large inter-DGVMs exist, causing uncertainty in their projections (Haverd et al., 2020;Schimel et al., 2015).For example, an ensemble of DGVMs from the Trends in Net Land Atmosphere Carbon Exchanges project version 9 (TRENDYv9) showed that the simulated increase in CFE-induced plant biomass in NHee during 1850-2019 varied from 11.25 Pg C in the LPX-Bern model to 54.15 Pg C in the CABLE-POP model.Observations from site-based CO 2 -enrichment experiments and the outputs from DGVMs are not directly comparable, so a newly emerged mathematical approach, the so-called emergent constraint (Cox et al., 2013;Eyring et al., 2019;He et al., 2020;Lian et al., 2018; see Section 2), can be used to reconcile the spatialtemporal mismatch between them, for better understanding the historical effect of CO 2 fertilization on plant biomass.
Our overall aim is to quantify the specific role of CO 2 fertilization, based on both models and observations, in driving the biomass C sequestration in NHee ecosystems.We first used emergent constraint approach to derive the response of plant biomass to increasing [CO 2 ] (parameter ) in NHee from 1850 to 2019, resolving into all total vegetation, or forest/grassland only (see Section 2).Then, for the shorter period of 1990-2019, we presented much more specific information, with the fertilization effect disaggregated into further biome classifications and for six key different countries/regions in NHee.Using the inferred constraints in combination with observed forest biomass changes reported previously (Pan et al., 2011;Xu et al., 2021), we thus estimated the contribution of CFE for whole NHee and each of the six major countries/regions since the 1990s.
Our database is comprised of 41 field observations collected at 24 CO 2 -enrichment experiments in forest and grassland systems (Table S1; Figure S1), together with the outputs of the 10 DGVMs from the TRENDY ensemble.Since CO 2 -enrichment experiments and DGVMs have very different settings for [CO 2 ] increases (∆CO 2 ), there are two important transformations to be addressed to align the models with the data before emergent constraint (Figure 1).First, to simulate the step increase in [CO 2 ] common in CO 2 -enrichment experiments, we calibrated the DGVM outputs, originally simulating a gradual increase in [CO 2 ], over the grid cells corresponding to the experimental sites using a two-box model (Liu et al., 2019).Second, to match the magnitude of ∆CO 2 at historical levels (e.g., from 286.46 to 409.39 ppm during 1850-2019), we converted the experimental CFE, originally conducted at future [CO 2 ] levels (e.g., from ~380 to ~610 ppm) using a theoretical model of photosynthesis (P-model; see Section 2).

| Study region
Our study region was restricted to NHee, defined as the areas between 23° N and 90° N including a southern part of China below 23° N for integrity (Figure S1).We used the Collection-5 MODIS land-cover product (MCD12C1) with the classification scheme of the International Geosphere-Biosphere Programme at a resolution of 0.05° for 2001.This NHee vegetation map was first aggregated to a resolution of 0.5° by dominant types of land cover, and the 17 land-cover types were then reclassified into four classes by biome: studies from Song et al. (2019) and Terrer et al. (2019).We screened all these data using the following criteria: (1) experiments containing fertilized/high water/high temperature/drought treatments were all excluded because of the interactions; (2) measurements from the same experiment but different years were not independent, so only observations in the last year were selected; (3) different species in the same experiment were considered independent; and (4) experiments were only included when both above-and below-ground biomass per unit area were provided.A factor of 0.5 was used to convert dry-matter content to carbon content.In total, we compiled 41 plant biomass data (under the elevated [CO 2 ] and control treatments) from 24 CO 2 -enrichment experiments (Figure S1; Table S1).
For each site, the response of plant biomass to increasing [CO 2 ] (

| Simulated CFE on plant biomass from an ensemble of DGVMs
We used the outputs of biomass C stock (variable 'cVeg') during 1850-2019 from an ensemble of 10 DGVMs in the TRENDYv9 project (CABLE-POP, DLEM, IBIS, ISAM, ISBA-CTRIP, JSBACH, JULES-ES-1.0,LPJ-GUESS, LPX-Bern, and ORCHIDEEv3).All DGVMs were coordinated to perform consistent factorial simulations based on the TRENDY intercomparison protocol using the same climate forcings, increasing atmospheric [CO 2 ], and products of changes in land use and land cover (Ballantyne et al., 2012;Sitch et al., 2015).All modeled outputs were resampled to a common spatial resolution of 0.5° × 0.5° using the nearest neighbor method.
We first extracted the CFE-induced change in biomass over time   where ΔBiomass(i) is the CFE-induced change in biomass C density (S1-S0) in year i relative to the first year, and NPP(i) is the CFE-induced mean annual NPP (S1-S0) in year i.
We then calculated the analytical solution to Equation (3) as: Second, the original simulated ΔBiomass and NPP for 1850-2019 were used to fit the parameter μ for each site and each model based on Equation (4).As a result, site-scale ΔBiomass (ΔBiomass site ) emulated by the two-box model generally agreed well with the original simulated ΔBiomass site in most of the DGVMs, indicating that the fitted parameter μ could be used to represent the turnover rate of biomass in this study (Figure S9).
Third, we assumed that [CO 2 ] increased abruptly from CO First 2 to CO Last 2 and then remained stable during the simulation period for each experimental site and each model, that is, simulated NPP would also increase abruptly in response to increasing [CO 2 ], which was assumed to be constant during the experimental period (i in Equation 4) and calculated as the difference between the CFE-induced mean annual NPP under CO Last 2 and CO First 2 .ΔBiomass site for each model could thus be reproduced using Equation ( 4), applying the constant parameter μ obtained above.At last, the simulated site-scale CFE on plant biomass C density ( Mod Site ) was equal to ΔBiomass site divided by the difference between CO Last 2 and CO First 2 (Equation 2).

| Conversion of observed site-scale CFE on plant biomass at future [CO 2 ] to CFE at historical [CO 2 ]
The difference in atmospheric [CO 2 ] between the control (e.g., about where ΔNPP Hist Site and ΔNPP Obs Site are changes in CFE-induced NPP at historical and future [CO 2 ], respectively.C-use efficiency (CUE, the ratio of NPP to gross primary productivity (GPP)) was assumed to be constant across historical and future [CO 2 ] in the models (Figure S11), so the site-scale CFE on plant biomass C density at historical [CO 2 ] ( Hist Site ) could be expressed as: (2) where ΔGPP Hist Site and ΔGPP Obs Site are changes in CFE-induced GPP at historical and future [CO 2 ], respectively, which could be simulated using an optimality-based model of light-use efficiency (P-model; Wang et al., 2017).
The validity of the P-model for predicting GPP has been successfully evaluated based on extensive measurements of ecosystem flux and field studies (Wang et al., 2017).The rpmodel R package provides us with a simple but effective way to estimate CFE-induced GPP change under historical and future [CO 2 ] on the basis of several environmental variables (Table S2).A summary of how the P-model derives CO 2 fertilization effects is provided in Supporting Information.
The underlying optimality assumption of the P-model implies that it can represent the effects of climate (e.g., temperature, water vapor pressure deficit, atmospheric pressure) based on the first principles rather than imposed empirical functions.The formulation of the Pmodel thus allows us to account for the CO 2 fertilization effects on the LUE, which were mainly ignored in previous LUE frameworks (De Kauwe et al., 2016).We simulated GPP for each site at four [CO 2 ] , and CO C 2 ) based on the P-model corresponding to the sampling years in the CO 2 -enrichment experiments (Table S1).

ΔGPP Obs
Site and ΔGPP Hist Site in Equation ( 8) were thus equal to the difference between GPP T and GPP C and between GPP Last and GPP First , respectively.Based on this, the observed site-scale CFE on plant biomass at future [CO 2 ] levels ( Obs Site ) can be extrapolated to the CFE at historical [CO 2 ] levels ( Hist Site ).We note that there is a substantial geographical spread in Hist Site , with a small proportion of observations yielding negative values.
To conduct further emergent constraint analysis, it is essential to obtain the cross-site-average Hist Site and its corresponding uncertainty.This can be accomplished through a bootstrap resampling approach.
First, we randomly selected one Hist Site from each site and then averaged all selected Hist Site across sites.Second, we repeated the first step 1000 times to generate a sample of 1000 Hist Site .Third, we calculated the mean ( Hist Site ) and standard deviation Hist Site of all the 1000 Hist Site samples, assuming a Gaussian distribution for all observations.

| Emergent constraints of CFE on plant biomass
The concept of emergent constraint is founded on the premise that the ensemble of models may be applicable for establishing a heuristic relationship between a variable that can be directly measured/observed and another (i.e., the variable of interest) that cannot, even though large uncertainties remain between models (Cox et al., 2013).To constrain the response of plant biomass C density to increasing [CO 2 ] in NHee, we first used least-squares linear regression to assess the relationship between Mod Site and Mod NHee across the ensemble of DGVMs, defined as: where a is the slope and b is the intercept of the regression line, respectively.
The least-squares error of this regression model was calculated as: where M is the number of DGVMs, (M − 2) is the degrees of freedom of the linear regression, and ̃ Mod NHee represents the predicted Mod NHee based on Equation (9).
We then re-estimated The PDF for emergent-constrained ̃ Mod NHee could thus be estimated as: where P Hist Site is the PDF of the variable Hist Site .We also applied the same method of emergent constraint to explore the model performance in different biomes and different regions/countries separately.Correlating simulated Mod Reg and Mod

Site
for each of the six major countries/regions, however, was difficult, because not all the regions have the corresponding CO 2 experimental sites (Figure S1).We found that simulated Mod Reg in the United States increased linearly across the DGVMs, with both regional site- Site for all temperate-forest sites (Figure S6) and forest sites in the United States only (Figure S12).That is to say, Mod Site of boreal or temperate forests could be used to replace regional Mod Site to correlate with regional CFE-induced increase, as an alternative.
We thus used observed site-average Hist Site of boreal and temperate forests to constrain the simulated Mod Reg in boreal (Figures S3-S5) and temperate (Figures S6-S8) countries/regions, respectively.

| Analysis of robustness
The robustness of the constrained CFE on plant biomass was confirmed by reducing the number of experimental sites.We first selected (n − m) experimental sites from a total of n sites to constrain the simulated response in the DGVMs, where m is the number of ( 9) experimental sites (m = 1, 2, 3) removed before the calculation.We then repeated this step until all combinations were selected.The combinatorial number was Finally, we calculated the mean and standard deviation for all the above combinations to represent the average constrained CFE and its uncertainty.

| Calculation of the observed changes in forest biomass since the 1990s
The observed changes in forest biomass over the last three decades in NHee were collected from an inventory-based study (Pan et al., 2011) and a satellite-based product (Xu et al., 2021) (Xu et al., 2021).This data set was produced at a spatial resolution of 10 km by compiling a large number of ground inventories, airborne laser scanning, and satellite lidar data to train a self-improving machine-learning model using systematic time series observations from microwave and optical satellite imagery.The conversion of above-to below-ground live biomass was based on vegetation specific root: shoot ratios (see table S5 in Xu et al., 2021).Note that detecting biomass gain in intact forests with high biomass can be challenging due to the "slow in-fast out" characteristic, as well as the potential saturation of satellite signals in dense forest (Harris et al., 2016;Smith et al., 2020).To address this, Xu et al. (2021) established an approach to adjust the long-term estimates of C changes in large-biomass forests.We follow the same approach.We first separate the biomass changes in large-biomass forests (defined as areas with AGB > 100 Mg ha −1 ) and low-biomass forests (defined as areas with AGB ≤ 100 Mg ha −1 ).Specifically, the biomass changes in large-biomass forests were estimated by multiplying the biome-specific growth rates (see table S4 in Xu et al., 2021) by the areas of the regions, whereas the biomass changes in low-biomass forests were estimated using linear regression over time.The sum of the biomass changes in large-biomass forests and low-biomass forests was considered as the ΔStorage Obs .
Furthermore, the change in the density of forest biomass C was estimated as ΔDensity Obs = ΔStorageObs Area , where Area is the total forest area in 2001 in NHee based on MODIS land-cover data, without considering the changes in forest area over time.

| Calculation of the CFE-induced changes in forest biomass since the 1990s
We calculated the CFE-induced changes in the density of biomass C (ΔDensity CFE , g C m −2 year −1 ) and changes in the storage of biomass C (ΔStorage CFE , Tg C year −1 ) in NHee as follows: where

| Constrained CFE on plant biomass in extra-tropical Northern Hemisphere
The simulated response of plant biomass to time-evolving smoothly increasing [CO 2 ] in the entire NHee ( Mod

NHee
) was tightly correlated with the simulated site-average response of plant biomass to shortterm step increases in [CO 2 ] ( Mod Site ) across the DGVMs (r = .73,p = .0175;Figure 2b).This strong linear relationship confirmed the existence of our proposed emergent constraint, and models with a large site-scale Mod Site also produced a large NHee-scale Mod NHee , as expected.Combined with the observed site-average response measurement, with a value at historical [CO 2 ] levels ( Hist Site ) of 0.29 ± 0.02 kg C m −2 [100 ppm] −1 (mean ± standard deviation), the model-derived heuristic relationship provided a constrained Mod NHee of 0.49 ± 0.12 kg C m −2 [100 ppm] −1 for period 1850-2019.This more refined emergent constraint-based Mod NHee was slightly lower than the original multi-model average (0.53 ± 0.15 kg C m −2 [100 ppm] −1 ).Furthermore, the relationship shown in Figure 2b generated   NHee by 20% after constraint (Figure 2c).Additional sensitivity analyses are performed that confirm the robustness of the constrained Mod NHee when using less observational CO 2 -enrichment sites in the emergent constraint, separately (Figure S13).

| Contribution of CFE to forest biomass change
The growing amount of inventory data and the development of satellite remote sensing enable a rigorous characterization of the growth trajectory of forest biomass over the past three decades (Harris et al., 2021;Liu et al., 2015;Pan et al., 2011;Xu et al., 2021).Such data provided us with an excellent opportunity to examine the contribution of CFE to the observed overall increases in NHee forest biomass (see Section 2).(Pan et al., 2011) and satellites (Xu et al., 2021), our emergent constraint analysis of Figure 2 was revisited, but now for a shorter period of years 1990-2019.Using the same data from CO 2 -enrichment sites, in tandem with the emergent constraint, we obtained a lower response of plant biomass to rising [CO 2 ] in NHeescale total forests ( Mod Forest , 0.86 ± 0.28 kg C m −2 [100 ppm] −1 , Figure 3) compared to that for 1850-2019.As before, the robustness of the constrained Mod Forest was confirmed by the analyses using fewer observational CO 2 -enrichment sites for the emergent constraint (Figure 3; Figure S2).We then multiplied this constrained Mod Forest with a historical ∆CO 2 of 56.19 ppm (from 353.20 in 1990 to 409.39 ppm in 2019), the result indicated that CFE has increased the forest biomass C storage by 248.04 ± 80.87 Tg C year −1 (ΔStorage CFE ; Figure S14).This CFEdriven change was then compared against the observed total changes (ΔStorage Obs ) from our two datasets for biomass.We therefore identified that CO 2 fertilization alone explained 54 ± 18% of the observed increase (461.11Tg C year −1 ) reported from forest inventories over the past three decades (Table S3).Using an independent satellite-derived biomass product, the dominant role of CFE was further confirmed, suggesting a contribution of 64 ± 21% for the same period (Table S3).
The CFE may be different between temperate and boreal forests due to their contrasting environmental limitations (Terrer et al., 2019).Splitting our analysis into these biomes, we found a major difference, where the constrained Mod Temp of 1.04 ± 0.31 kg C m −2 [100 ppm] −1 for temperate forests was three times the size of the constrained value for boreal forests, Mod Boreal of 0.36 ± 0.26 kg C m −2 [100 ppm] −1 (Figure 3; Figure S2).Combining these two values with the estimates of satellite-derived biomass increase (countrylevel inventory data cannot be recalculated to match the regional extent of each biome), the contribution of CFE was found to have made a larger contribution to temperate forests (71 ± 21%) than boreal forests (34 ± 25%).
F I G U R E 3 Constrained CO 2 fertilization effect (CFE) on temperate and boreal forest biomass in Northern Hemisphere extra-tropical ecosystems during 1990-2019.The first number in each circle indicates the multi-model average of originally simulated CFE (Unconstrained), constrained CFE using all observational sites (Constrained) and constrained CFE using one (Robustness 1 ), two (Robustness 2 ), or three (Robustness 3 ) fewer observational sites for total forest, boreal forest, and temperate forest, respectively.The standard deviation of CFE is shown in brackets.

| Regional differences in the contribution of CFE to forest biomass change
Carbon loss and gain in forest biomass may be strongly influenced by direct human disturbance and other environmental changes, which may alter the relative contribution of CFE at the regional scale.To further identify the role of CFE for the six major countries/regions in NHee, that is, Canada, northern Europe, Russia, the United States, temperate Europe, and China, we used the same fractional approach as outlined above (i.e., the calculation of ΔStorage CFE /ΔStorage Obs ; Materials and Methods).Our analyses showed that CFE contributed the most to the inventory-based forest biomass increase in China (56 ± 19%), followed by northern Europe (39 ± 20%), the United States (29 ± 9%), and Russia (28 ± 24%), with a lesser contribution for about 13 ± 5% in temperate Europe (Figure 4; Table S3).Of particular note is that in Canada, a negative contribution of CFE was detected (−117 ± 71%; Figure 4; Table S3).Additional to our analysis based on forest inventory data (black regular triangles in Figure 4), we also examined results from an independent satellite-based biomass dataset (red inverted triangles in Figure 4), which inferred similar or larger proportions of CFE among most countries/regions.These satellitederived values have some similarities with the estimates from forest inventories, although there are particular differences for temperate Europe.

| DISCUSS ION
By constraining the DGVM simulations against CO 2 -enrichment experiments, we have shown that the historical response of forest biomass to increasing [CO 2 ] was underestimated by DGVMs (Figure 2d,e).Specifically, this underestimation may be associated with a simulated N-cycle in DGVMs that overly restricts predicted CO 2 fertilization on photosynthesis, as models with C-N interactions generally produced a strong and progressive N limitation on net primary production for increasing [CO 2 ].These models, specifically, may underestimate plant N uptake compared to observations (Zaehle et al., 2014), probably because of the missing simulation of high N uptake through ectomycorrhizal tree-fungal associations (Terrer et al., 2016(Terrer et al., , 2018)).Furthermore, the effect of CO 2 fertilization on plant biomass is dependent on a fine balance between changes in photosynthesis and changes in turnover, where the latter is mainly dependent on the C allocation to different plant components, simulated tissue lifespan as well as whole-plant mortality rates (Koven et al., 2015).Current models focus more on photosynthesis, with a common lack of representation of the mechanisms behind turnover-related processes (Friend et al., 2014;He et al., 2021;Koven et al., 2015;Sullivan et al., 2022).For example, most models estimate an overall higher baseline woody turnover rate, as well as divergent turnover response to rising [CO 2 ], when compared to the observations at Oak Ridge and Duke free-air CO 2 enrichment sites.This deficiency in process representation and parameterization further leads to large uncertainties in biomass carbon accumulation (De Kauwe et al., 2014).
In addition, our results show that the historical response of grassland biomass to increasing [CO 2 ] in grassland was generally overestimated by DGVMs (Figure 2f,g).We considered four possible explanations that may account for this major difference.First, grasslands are mainly distributed in temperate regions (Figure S1), where seasonal variations in precipitation could strongly regulate terrestrial plant growth (Hovenden et al., 2014;Obermeier et al., 2017).
A previous study suggested that annual CFE on biomass production was largely offset by the opposite response of CFE between the increase in spring and non-spring precipitation (Hovenden et al., 2019), but such mechanism is not considered in current models.Second, most land-surface models have only a limited description of plant functional types, including models for grasses.Recent studies have found that changes in the composition of grassland species due to higher [CO 2 ] may lead to the loss of C, generating a lower (and observed) impact of CFE on grassland biomass (Mueller et al., 2016;Song et al., 2019;Zhu et al., 2020).Third, a global synthesis of CO 2 experiments suggested that grasslands tend to accumulate more C in soils, rather than in biomass, when exposed to elevated [CO 2 ] (Terrer et al., 2021).However, this tradeoff, as well as the complex interactions of root−microbe−mineral in grass rhizosphere, are not currently reproduced in models.Finally, our inferred low CFE-induced biomass increases in grassland may be also associated with the low efficiency in N uptake in arbuscular mycorrhizal plants, including grassland species (Terrer et al., 2018), and that once again is not represented in DGVMs.All these possible omissions in models may account for the overestimation of simulated CFE, indicating that the biomass C sink capacity of grassland is far lower than previously expected.
In terms of the contribution of CFE to plant biomass since the 1990s, our findings suggest a larger contribution of CFE in temperate forests compared to boreal forests.This difference could be explained by the stronger limitations on plant growth due to N constraints (Du et al., 2020;Terrer et al., 2019;Vallicrosa et al., 2021) and low temperatures (Keenan & Riley, 2018;Zhu et al., 2016) in boreal ecosystems.However, we also observed a negative contribution of CFE in Canada.One possible reason for this negative contribution is that the storage-based metric (ΔStorage CFE /ΔStorage Obs ) used here is more reliable in regions with substantial increases in forest area (Figure S15), while a density-based metric (ΔDensity CFE / ΔDensity Obs ) may be better suited for regions undergoing forest loss, that is, Canada.Nevertheless, even when using the densitybased metric, the contribution of CFE in Canada was still negative (−97 ± 59%; Figure 4; Table S3).Therefore, we suggested that the negative impacts on forest biomass were primarily due to the significant influence of disturbances in Canadian forests, and in particular pest attacks and wildfires (Kasischke & Turetsky, 2006;Wang et al., 2015).These disturbances were so severe that they exceed any potential positive effects of CFE, as well as any other factors, resulting in a net decrease in forest biomass.
Considering the differences in CFE contribution between satelliteand inventory-based estimates (Table S3), there are several possible explanations.For example, some lack of agreement may be attributable to the differences in the observational periods.The satellite data were integrated over years 2000-2019, and therefore potentially captured more loss of forest biomass C from post-2007 disturbances in boreal and temperate regions (Forzieri et al., 2021;Seidl et al., 2014), which was not captured by the forest inventory data based on the years 1990-2007.In addition, previous studies have suggested that optical and microwave satellite data with short wavelengths may have inadequate sensitivity to detect gradual forest gain, compared with their good ability to observe more abrupt forest loss signals, such as deforestation (Bartels et al., 2016;Li et al., 2017;Tian et al., 2015).
This deficiency is more evident in dense forests, where the data may underestimate the increase in observed woody biomass due to CFEinduced increases in light-use efficiency (Smith et al., 2020), and thus result in a greater fraction of CFE contribution.
An additional theoretical possibility of finding low CFE contributions is for regions characterized by large-scale forest areal expansion.A large C accumulation for trees may be expected due to rapid growth through direct afforestation/reforestation (Chen et al., 2019), rather than the fertilization effect, in areas of planting and regrowth.Contrary to expectations, however, the country with the largest gain in forest area (Pan et al., 2011), that is, China, also had the highest percentage contribution of CFE compared to that of other regions (Figure 4; Table S3).We do, though, suggest some caution regarding the robustness of our predicted contributions at the regional scale due to a dearth of sufficient CO 2 experiments, including both above-and below-ground information, to adequately constrain regional CFE.Further validation is also necessary, potentially through enhanced spatial sampling with longer-duration CO 2enrichment experiments (Jones et al., 2014).
Furthermore, recent studies have indicated an increasing constraint of climate on vegetation growth, which may ultimately counteract the potential growth benefits of CFE (Green et al., 2019;Jiao et al., 2021;Peñuelas et al., 2017;Yuan et al., 2019).Our analysis revealed a small decrease in constrained CFE for forests from 1850-2019 (0.94 ± 0.27 kg C m −2 [100 ppm] −1 ) to 1990-2019 (0.86 ± 0.28 kg C m −2 [100 ppm] −1 ), caution should be exercised in drawing conclusions about a decline in CFE solely based on these findings.We note that same site-scale CO 2 -enrichment observations at future Additionally, the application of the P-model, as needed to convert site-scale response of GPP to rising [CO 2 ] to expected changes between historical [CO 2 ] levels and future [CO 2 ] levels, contains implicitly a saturation effect over time (Figure 1; Liu et al., 2019), which may eventually result in a predicted slowdown of CFE between the two time periods.Whether the positive effects of CFE on biomass carbon enhancement may be expected to continue to outweigh the negative effects of climate-induced biomass losses in the immediate future remains an important issue to be addressed.

| CON CLUS ION
In summary, our study presents a simple, highly intuitive, and effective approach to determine the magnitude of the large-scale influence of CO 2 -induced fertilization on plant biomass.Of particular interest is that we have formally isolated the contribution of CFE from other non-CFE drivers.We have achieved this by exploiting multiple observational constraints, dynamical global vegetation models, and an emergent constraint technique to fuse together such observational data and models.Our study indicated that CO 2 fertilization is the dominant driver of the observed forest biomass increase over the last decades across the NHee.Specifically, inventory-and satellite-based evidence suggested the fertilization contribution was 54 ± 18% and 64 ± 21%, respectively.Eventually, the positive effect of CO 2 fertilization may slow down and saturate in the process of achieving carbon neutrality (Franks et al., 2013;Peñuelas et al., 2017;Piao et al., 2022;Winkler et al., 2021), which will further lower the terrestrial ecosystem C sink capacity.Consequently, the presentation of CFE as a solution for long-term climate mitigation merits caution, potentially requiring additional CO 2 -offset options (e.g., renewable power utilization, deliberate afforestation/reforestation), to achieve carbon neutrality.

AUTH O R CO NTR I B UTI O N S
Shilong Piao and Yongwen Liu designed the research; Yue He and Lingjie Lei performed the analysis and wrote the draft of the manuscript; and all authors contributed to the interpretation of the results and to the text.

(
written here as Biomass) as the difference in simulated cVeg between the S1 (time-varying [CO 2 ] only) and S0 (constant [CO 2 Schematic overview of the methodology used for the transformation of [CO 2 ]-related settings between model and data.Green arrows show the methodological steps to convert observed site-scale CO 2 fertilization effect (CFE) on plant biomass at future [CO 2 ] levels ( Obs Site ) to CFE at historical [CO 2 ] levels ( Hist Site ) for the CO 2 experiments.Brown arrows represent the methodological steps to convert simulated site-scale CFE on plant biomass with a gradual increase in [CO 2 ] to a step increase in [CO 2 ] ( Mod Site ), and this value together with the simulated Northern Hemisphere extra-tropical ecosystems (NHee)-scale CFE ( Mod NHee ) were then used to build a linear relationship across the dynamic global vegetation models for further emergent constraint.All [CO 2 ] values shown are rounded down to integral numbers for clarity.More detailed information can be found in Section 2. for each model.The simulated response of plant biomass to increasing [CO 2 ] in NHee ( Mod NHee ) for each model was then calculated as: where Biomass Last and Biomass First represent the area-weighted average of CFE-induced change in biomass C density (S1-S0) in NHee for the last year and the first year of the simulated period, respectively, [CO 2 ] for the last year and the first year of the simulated period, respectively.The first year was 1850 when estimating

First, we obtained
the simulated CFE-induced change in net primary productivity (NPP) and biomass C density (S1-S0) for 1850-2019 for each experimental site and each model by averaging the simulated values within a 4.5° × 4.5° window around the corresponding site.Our choice of the window size was adopted from Liuet al. (2019).We assumed that the turnover rate of the biomass C pool (μ) was constant, so the change in biomass C density resulted from the change in NPP could be represented by: 380 ppm) and elevated [CO 2 ] treatments (e.g. about 610 ppm) in the CO 2 -enrichment experiments was much larger than the difference in the DGVMs (e.g., from 286.46 ppm in 1850 to 409.39 ppm in 2019), preventing the direct comparison between the DGVMs and the CO 2enrichment experiments (Figure 1).To transform the observed CFE on plant biomass at future [CO 2 ] levels in the experiments to that at historical [CO 2 ] levels, we needed to extrapolate the change in plant biomass from CO C 2 to CO T 2 to the change in plant biomass from CO First 2 to CO Last 2 .Assuming that the turnover rates of the biomass C pool were unchanged with an increase in [CO 2 ] in the DGVMs (Figure S10), changes in biomass caused by a step increase in [CO 2 ] at historical (ΔBiomass Hist Site ) and future (ΔBiomass Obs Site ) [CO 2 ] levels could be calculated using Equation (4) as: thus, , respectively.Since forest inventory data have documented changes in forest area for different countries/regions, both the observed change in the density of biomass C (ΔDensity Obs , g C m −2 year −1 ) and change in the storage of biomass C (ΔStorage Obs , Tg C year −1 ) for whole NHee and each of the six major countries/regions in NHee were analyzed.First, we extracted the ΔDensity Obs and ΔStorage Obs during 1990-2007 from inventories as follows, since the area of forest cover has changed over time: where Area 1990 and Area 2007 denote the total forest area for 1990 and 2007, respectively, and Storage 1990 and Storage 2007 denote the total storage of biomass C for 1990 and 2007, respectively, all of which could be obtained from tables S2 and S3 in Pan et al. (2011).Δyear is the length of the study period which is set to 18. Second, we estimated ΔDensity Obs and ΔStorage Obs using a newly released satellite-based dataset of global C stocks in woody vegetation during 2000-2019

F
Constrained CO 2 fertilization effect (CFE) on plant biomass in Northern Hemisphere extra-tropical ecosystems (NHee) during 1850-2019.(a) CFE-induced change in plant biomass during 1850-2019 in the dynamic global vegetation models (DGVMs).Note that changes in plant biomass were all shown with a 30-year moving window for demonstration purposes.(b) Relationship between the simulated response of plant biomass to increasing [CO 2 ] in NHee ( Mod NHee ) and the simulated site-average response of plant biomass to increasing [CO 2 ] ( Mod Site ) across the DGVMs.The vertical gray area represents the observed response of plant biomass to increasing [CO 2 ] from CO 2 -enrichment experiments ( Obs Site , mean ± standard deviation).The horizontal pink lines represent the constraint estimate of Mod NHee (solid line) ± standard deviation (dashed lines).(c) Probability density function of Mod NHee before (black line) and after (pink line) constraint.(d-g) Same as (b) and (c) except the field observations of Obs Site are divided into forest and grassland to constrain CFE on plant biomass in (d, e) forest ( , 15, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16806by JOSEP Penuelas -Csic Organización Central Om (Oficialia Mayor) (Urici) , Wiley Online Library on [05/07/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License DGVMs, lowering the range of uncertainty of Mod Unlike the CO 2 -enrichment experiments, which are designed to capture in isolation the impacts of elevated [CO 2 ], the satellite retrievals and inventories capture the trends in plant biomass due to all forcings combined (notably climate change, CFE, N deposition, and land-use change).To match the temporal extent of the biomass observations from inventories

F
Contribution of CFE to changes in forest biomass for Northern Hemisphere extra-tropical ecosystems (NHee) and each of the six major countries/regions since the 1990s.Panels (a) to (g) represent the regions of NHee, Canada, northern Europe, Russia, the United States (USA), temperate Europe, and China, respectively.For each panel, the left side shows the CFE-induced change in the density of biomass carbon (ΔDensity CFE ), and the right side shows the CFE-induced change in the storage of biomass carbon (ΔStorage CFE ).Bars represent the changes in biomass caused by CFE before (gray bars) and after (blue bars) the emergent constraint.The error bars represent the standard deviations of the CFE-induced changes in biomass.Observed changes in biomass carbon density (ΔDensity Obs ) and storage (ΔStorage Obs ) over the last three decades are derived from Pan et al. (2011; black regular triangles) and Xu et al. (2021; red inverted triangles) respectively.The green areas represent forest areas, and the gray areas represent non-forest areas.Map lines delineate study areas and do not necessarily depict accepted national boundaries.

| Conversion of simulated site-scale CFE on plant biomass with a gradual increase in [CO 2 ] to a step increase in [CO 2 ]
Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16806by JOSEP Penuelas -Csic Organización Central Om (Oficialia Mayor) (Urici) , Wiley Online Library on [05/07/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ModForest is the response of forest biomass C density to increasing [CO 2 ] in NHee for 1990-2019, Mod Forest,S is the response of forest biomass C storage to increasing [CO 2 ] in NHee for 1990-2019, ΔCO 2 is the difference in the global atmosphere [CO 2 ] between 2019 and 1990, and Δyear is the length of the study period which is set to 30.