CO₂ fertilization of terrestrial photosynthesis inferred from site to global scales

Three metrics are used to quantify and assess the response of GPP to CO₂.

We find a strong increasing trend of GPP measured across 632 site-years of observations in the EC network of 9.1 gC m⁻² year⁻² between 2001 and 2014 [interquartile range (IQR) ∈ [8.5, 11.5], P < 0.01, ${\text{β}}_{\mathit{app}}^{\text{ln}}$ = 1.24]. This EC-inferred trend is captured by our EEO framework (8.3 gC m⁻² year⁻², IQR ∈ [7.0, 9.4], P < 0.01, ${\text{β}}_{\mathit{app}}^{\text{ln}}$ = 1.12) (Figs. 1A and 2A). Both ${\text{β}}_{\mathit{app}}^{\text{ln}}$ values are relatively high compared with expectations from measured direct responses to CO₂ (i.e., ${\text{β}}_{\mathit{dir}}^{\text{ln}}$) but are similar to several published indirect apparent GPP proxies (8): carbonyl sulfide records (${\text{β}}_{\mathit{app}}^{\text{ln}}$ = 0.95) (39), ice-core measurements of atmospheric O₂ isotopes (${\text{β}}_{\mathit{app}}^{\text{ln}}$ = 1.3) (40), and satellite monitoring of water-use efficiency (${\text{β}}_{\mathit{app}}^{\text{ln}}$ = 1.1) (41). Since ${\text{β}}_{\mathit{app}}^{\text{ln}}$ is derived from natural environments and influenced by CO₂ and climate variability, it differs from the direct ${\text{β}}_{\mathit{dir}}^{\text{ln}}$ discussed later. The detected overall trends are robust to the uncertainty caused by filtering of low-quality data, and the uneven distribution of sites and site-years (SI Appendix, Fig. S2). Variations in the EC-inferred GPP due to partitioning methods and friction velocity filtering methods have a minimal effect on the EEO-inferred GPP (shaded area in Fig. 1A), where the former is used to calibrate the canopy upscaling of the latter. In general, the EEO-inferred GPP reproduces the interannual variability (IAV) of the EC-inferred GPP (Pearson’s r = 0.91) (Fig. 1A; SI Appendix, Fig. S1). The high level of agreement from annual (Fig. 1B) to monthly (SI Appendix, Fig. S3) scales supports the use of the EEO framework for attribution analysis.

To quantify the GPP trends contributed by different factors, we use the EEO framework to perform a univariate sensitivity analysis (SI Appendix, Fig. S4). The aggregated trend from individual factors (10.3 gC m⁻² year⁻², IQR ∈ [7.7, 11.0]) is similar to the EEO-inferred and EC-inferred trends (Fig. 2A). We find a strong direct contribution from c_a at the site level, which accounts for 44% (4.5 gC m⁻² year⁻², P < 0.001) of the aggregated GPP trend (Fig. 2A). We also diagnosed another overall CO₂-induced GPP trend using the partial differential approach (4.9 gC m⁻² year⁻², the first term on the right-hand side of Eq. 1), resulting in a small difference (0.4 gC m⁻² year⁻²) compared with the univariate analysis. We then convert the estimate from the univariate analysis, which corresponds to a direct response ratio (${\text{β}}_{\mathit{dir}}^{\text{ln}}$) of 0.61 for GPP to CO₂. This ${\text{β}}_{\mathit{dir}}^{\text{ln}}$ is lower than those light-saturated leaf-level analyses in FACE experiments (${\text{β}}_{\mathit{dir}}^{\text{ln}}$ = 0.79) (4, 8) and those using deuterium isotopomers (${\text{β}}_{\mathit{dir}}^{\text{ln}}$ = 1.0) (8, 42) but is comparable to a theoretical ${\text{β}}_{\mathit{dir}}^{\text{ln}}$ of 0.6 in Ref. 8 that considers canopy radiative transfer, suggesting a coexisting of light-saturated and light-limited photosynthesis over our study sites. CFE is stronger under light-saturated conditions but is present regardless of light conditions (27, 28). An increase in c_a eventually elevates leaf c_i through optimization (SI Appendix, Eq. S16), which in the FvCB scheme leads to higher light use efficiency and carbon assimilation rates (7, 32).

In addition to c_a, T_a plays the second most important role, accounting for 28% (2.9 gC m⁻² year⁻², P < 0.05) of the aggregated trend (Fig. 2A), mainly because of the effect of warming on photosynthetic capacity, i.e., the maximum rate of Rubisco activity, V_cmax, and the maximum rate of electron transport, J_max (Fig. 2C) (38). Although increasing temperature negatively affects GPP through, for example, vapor pressure deficit (VPD) and some biochemical reaction pathways (Fig. 2C), these negative effects are weaker than the positive effects of temperature (Fig. 2 A and B). SWC and specific humidity, as proxies for plant water supply and atmospheric water demand, together contribute 14% of the aggregated trend (0.91 and 0.50 gC m⁻² year⁻², P > 0.05) but are statistically insignificant. These two factors (i.e., SWC and q_a), however, are the largest contributors to GPP IAV, and their contributions (median values are 48 and 33 gC m⁻² year⁻¹, respectively; Fig. 2B) are an order of magnitude larger than those of the GPP trends (Fig. 2A). These results are consistent with previous data and ESM-driven studies on plant responses to water (43–47) and are further supported by observations that GPP anomalies are regulated by stomatal conductance and leaf c_i, as plants attempt to maximize their efficiency of carbon gain per unit water loss (30, 32, 48, 49).

Other factors, such as long-term changes in LAI and radiation, could also lead to changes in GPP across the EC network. However, we find a small and statistically insignificant trend in GPP due to increased LAI (0.88 gC m⁻² year⁻²; Fig. 2A), which is consistent with a previous study that solely considered the response of GPP to changes in LAI (50). This small LAI-induced trend in GPP is comparable to the LAI trend, both of which are on the order of ‰ year⁻¹ (51), although our framework diagnoses a large underlying sensitivity of GPP to LAI (mean $\frac{\partial \mathit{GPP}}{\partial \mathit{LAI}}$ = 119, intersite SD = 178, unit: gC m⁻² year⁻¹). In addition, light saturation in canopy radiative transfer (52) may lead to a low LAI-induced GPP trend as leaf area increases do not effectively translate to increases in the fraction of absorbed photosynthetically active radiation at high LAI (3, 20). For the changes in incoming shortwave radiation, their contribution to the aggregated GPP is also small (0.43 gC m⁻² year⁻²; Fig. 2A) since there is no large trend in radiation over the study sites at the annual scale. Therefore, we conclude that neither changes in LAI nor radiation significantly contribute to the observed upward trend in GPP. The effect of changes in surface pressure is also negligible (Fig. 2 A and B).

To examine global CFE patterns, we apply our framework to Moderate Resolution Imaging Spectroradiometer (MODIS) LAI and European Centre for Medium-Range Weather Forecasts Reanalysis Version 5 – Land data (hereafter, ERA5-land). Here, our EEO framework is calibrated by the ensemble mean of all eight satellite-derived GPP products. We then assess both the GPP trend directly estimated by the EEO framework and the underlying differential sensitivities, which themselves can be used to reconstruct a composite trend estimate (Eq. 1). The two approaches show almost identical apparent GPP trends of about 4.7% decade⁻¹ (P < 0.001) relative to their mean GPP throughout 2001 and 2016 (A0 and A1 in Fig. 3A). The resulting CO₂-induced GPP trend is about 4.1% decade⁻¹ (4.4 gC m⁻² year⁻², P < 0.001) (A2 in Fig. 3A; SI Appendix, Fig. S5). We note that the fractional contribution of c_a to the apparent GPP trend at the global-scale analysis may be uncertain, and we do not report it because of the uncertainties in apparent GPP trend caused by the uncertain trends in climate forcings. However, estimates of the CO₂-induced GPP trends (A2 in Fig. 3, calculated as the product of β_CO2 and Δc_a) are independent of the uncertainties in the climate forcing trends because of the nature of our partial differential approach: 1) The spatial variation of β_CO2 is rational because the climate forcings well capture their spatial variations at the global scale (53–55), and 2) Δc_a is derived from the National Oceanic and Atmospheric Administration (NOAA)’s observations with high confidence. As a result, after excluding the spatial heterogeneity introduced by LAI and canopy structure, we show a good agreement when comparing the leaf-level CO₂-induced GPP trend at the global-scale analysis with their EC site counterparts (r = 0.88; SI Appendix, Fig. S6A).

We compare our CO₂-induced GPP trends (A2 in Fig. 3) with simulations from 13 dynamic global vegetation models (DGVMs) and eight satellite-derived GPP products during their overlap period (2001 to 2016). Our estimates are at the upper end of their trend distribution, at least one-third stronger than their median (Fig. 3A). The lower GPP trends in the median of DGVMs (B1) and satellite-derived products (C1) is likely due to their smaller CFE in evergreen broadleaf forests (EBFs) (Fig. 3B). Our EEO framework suggests similar relative GPP trends caused by CO₂ among all C3 species (including EBFs), ranging from 4.3 to 5.1% decade⁻¹ (A2 in Fig. 3B). We further replace the ensemble mean of all satellite-derived GPP with each individual product to calibrate the EEO framework, finding a small spread in the diagnosed CO₂-induced relative GPP trends (SI Appendix, Table S2), despite the fact that these individual satellite GPP products used for calibration vary markedly (global GPP in 2001 with a minimum of MODIS Collection 5.5 [MOD-C55] at 104 Pg C year⁻¹ and a maximum of Pmodel-s0 at 131 Pg C year⁻¹). In EBFs, it is noteworthy that the satellite GPP products that only consider the indirect CFE on satellite reflectance, but not the direct CFE on gas exchange [i.e., FluxCom, MOD-C55, MOD-C6, and vegetation photosynthesis model (VPM)], do not report an increasing GPP trend (Fig. 3B). Their GPP trends are not a proxy of CFE because they model GPP solely based on climate forcings and leaf greenness or its surrogates. The climate forcings in reanalysis data may show negative impacts on GPP because of excessive temperature and drought stress in the tropics (see the EBF column in SI Appendix, Table S3), and the weak trends in leaf greenness cannot fully extrapolate the direct CFE effect because of the saturation of reflectance in EBFs (56–62). While the choice of climate forcing could introduce some uncertainties to our CFE estimates, we find similar diagnostic results with another forcing dataset [i.e., Climatic Research Unit and Japanese 55-year Reanalysis (CRU-JRA55) in SI Appendix, Fig. S7]. Furthermore, we calculate ${\text{β}}_{\mathit{app}}^{\text{ln}}$ and ${\text{β}}_{\mathit{dir}}^{\text{ln}}$ for the global-scale analysis (SI Appendix, Table S3), which results in the same conclusions as those using the relative GPP trend metric.

To further characterize the underlying CFE, we show the spatial distribution of β_CO2, which disentangles the confounding effects from other factors. The magnitude of β_CO2 is largest in hot and humid regions at the annual scale (Fig. 4A). Since the CO₂ trend is assumed to be constant worldwide, the spatial pattern of β_CO2 represents that of CO₂-induced GPP trends (A2 in SI Appendix, Fig. S5). At the ecosystem level, β_CO2 is highest in EBFs at about 5.8 gC m⁻² year⁻¹ per ppmv CO₂, followed by croplands and other forests from the temperate to boreal regions at about 2.2 gC m⁻² year⁻¹ per ppmv CO₂ (Fig. 4B). EBFs have the longest growing season, and their high temperature and radiation determine the high photosynthetic capacity (V_cmax and J_max) and thus β_CO2 (4, 28, 35). On the other hand, croplands cultivated in mild regions also have high photosynthetic capacity (63) and additional resource inputs from human management, including multiple cropping and irrigation (51), leading to a relatively high annual β_CO2. Compared with C3 plants, our results suggest a positive but much weaker β_CO2 in C4 plants (3, 4, 13), on the order of 0.6 gC m⁻² year⁻¹ per ppmv CO₂, because the C4 CO₂ pump inhibits the benefit of CFE (49).

According to plant photosynthesis optimization theories and gauged by the GPP magnitude estimated from an EC network or satellites, our EEO framework analytically constrains the sensitivity of GPP to CO₂ (i.e., β_CO2) from the gas exchange perspective. Our results suggest that the analytical solution of β_CO2 is insensitive to the GPP uncertainty used for calibration of the framework (Fig. 1A; SI Appendix, Table S2) and largely unaffected by interannual fluctuations of climatic conditions in recent years (Fig. 4B). This highlights the advantage of the EEO framework in that the partial derivatives decouple the confounding effects of climate on CFE. A recent study using regression methods showed a 40% decline in CFE from 1982 to 2015 (16). However, their key findings are subject to the large temporal uncertainty in the satellite data and limited by the causal ambiguity of their regression (64–66). In addition, our EEO framework assumes that plants coordinate nutrient investment to optimize photosynthetic capacity (V_cmax and J_max) with environmental constraints (34–37), rather than prescribing specific values for plant functional types. Photosynthetic capacity not only controls the slope of the GPP response to CO₂ but also strongly influences it by determining whether photosynthesis is light saturated or light limited (27, 28). We explore the influence of coordination at different timescales, comparing the results with the EC-inferred GPP (Fig. 1A; SI Appendix, Fig. S8), and adopt a decadal coordination to the environments in the peak LAI month for nutrient investment (67, 68) while adjusting the enzyme kinetics monthly with temperature (38). These assumptions result in a fixed reference J_max/V_cmax ratio for each site (or grid) over the study period while still allowing both to vary spatially to account for different growing climates. The decadal coordination allows for the presence of light-saturated and light-limited photosynthesis because the coordination is set at a much longer timescale than the stomatal and leaf c_i optimization (35). Under these conditions, the EEO framework diagnoses a large proportion of marginally light-saturated photosynthesis, which responds more strongly to changes in c_a than does light-limited photosynthesis (27, 28). Specifically, 63 to 69% of C3 canopy photosynthesis originates from light-saturated conditions for all seasons using the ERA5-land forcing. Together, these features explain the high β_CO2 predicted by our framework and also lead to a good agreement with the overall trend in EC-inferred GPP and supported by multiple independent studies (39–42).

We detect a large increase in GPP from a globally distributed EC network. Our EEO framework successfully captures the trends and IAVs in the EC-inferred GPP. The site-scale analysis suggests that CO₂ is the dominant contributor to the GPP trends, while SWC and specific humidity control the IAVs. Using this framework to scale the GPP globally, the estimated absolute CO₂-caused GPP trend is comparable to its EC-inferred counterpart and translates this CO₂ fertilization effect to a global increase in photosynthesis of 4.1% decade⁻¹ since the 2000s (or 5.1 PgC decade⁻² for global sum, 4.4 gC m⁻² year⁻² for global average). However, our EEO framework diagnoses a high CFE compared with the median values of the DGVMs and the satellite-derived GPP products, with the largest disagreement in tropical forests. These tropical forests are responsible for about one-third of global GPP, and thus, an accurate estimation of CFE in tropical forests is critical to global carbon cycle modeling. Finally, we urge expansion of the currently sparse EC observing network in tropical ecosystems to improve understanding of their physiological responses to CO₂ and climate change.

C.C., T.F.K., and W.J.R. are supported by the US Department of Energy Office of Science Biological and Environmental Research, Reducing Uncertainties in Biogeochemical Interactions through Synthesis and Computation Scientific Focus Area. T.F.K. acknowledges support from NASA Terrestrial Ecology Program Awards NNH17AE86I and 80NSSC21K1705. T.F.K. and I.C.P. acknowledge additional support from the LEMONTREE (Land Ecosystem Models based On New Theory, obseRvations and ExperimEnts) project, funded through the generosity of Eric and Wendy Schmidt by recommendation of the Schmidt Futures programme.