1 INTRODUCTION

Four independent lines of evidence indicate intensifying terrestrial biosphere activity (Campbell et al., 2017; Forzieri, Alkama, Miralles, & Cescatti, 2017; Graven et al., 2013; Le Quéré et al., 2018), namely the increasing positive trends in (a) the global net land carbon sink (NLS; Jinfeng et al., 2017; Le Quéré et al., 2018; Li et al., 2018); (b) the amplitude of the seasonal cycle (ASC) of atmospheric CO₂ (c_a) in the Northern Hemisphere (Forkel et al., 2016; Graven et al., 2013; Thomas et al., 2016); (c) satellite-observed leaf area (Forzieri et al., 2017; Zhu et al., 2016); and (d) global gross primary production (GPP), as inferred from long-term atmospheric carbonyl sulfide concentration (COS) records (Campbell et al., 2017). Attribution and future projection of these trends relies on terrestrial biosphere models (TBMs) that encapsulate mechanistic understanding of the terrestrial carbon cycle and are used to inform the global carbon budget (Le Quéré et al., 2018; Sitch et al., 2015).

Our focus in this work is on the historical increase in global GPP. The estimate presented by Campbell et al. (2017) [(31 ± 5)% increase (2σ, 1900–2010)] is based on the long-term atmospheric carbonyl sulfide (COS) records, derived from ice-core, firn and ambient air samples. Campbell et al. (2017) interpreted these records using a model that simulates the changes in atmospheric COS concentration according to the changes in its sources and sinks—including a large sink that is related to GPP—and suggested that their results provide a global-scale benchmark for historical carbon cycle simulations. The accuracy of COS-based GPP estimates depends in part on accounting for different environmental responses of COS and CO₂ uptake (Whelan et al., 2018). Despite recent progress in accounting for these different responses (Kooijmans et al., 2019), the estimate by Campbell et al. (2017) remains the only COS-based estimate of the historical change in global GPP and has yet to be utilized as a benchmark for global TBMs. Previous estimates of the global GPP increase over this timescale are restricted to process model simulations that suggest a much lower increase of 16 ± 4% (Keenan et al., 2016). While there is consensus that the simulated increase is largely driven by CO₂ fertilization (Baldocchi, Ryu, & Keenan, 2016; Keenan et al., 2016; Schimel, Stephens, & Fisher, 2015), the general underprediction of the global GPP increase means that a large fraction (>one-third) of the 31% ± 5% growth is unaccounted for by current TBMs, impacting their future projections of the dynamics of the NLS. Furthermore, it remains to be explored whether the large COS-based increase in global GPP is compatible with the other lines of evidence for increased terrestrial biospheric activity, particularly trends in the NLS and ASC.

The NLS may be estimated as the sum of atmosphere and ocean sinks minus the emissions from cement productions and combustion of fossil fuels. The NLS is positive (2.3 ± 0.4 PgC/year; 2007–2016) and has been increasing at a rate of 0.061 ± 0.02 (1 SE) PgC/year² (1980–2016; Le Quéré et al., 2018), implying that, globally, GPP is in excess of ecosystem carbon loss by fire and respiration, and that GPP is increasing faster than ecosystem carbon loss. As such, the magnitude and trend in the NLS provide indirect constraints on the magnitude of global GPP increase. Interannual variations in the NLS—driven largely by anomalies in plant productivity (Ahlström et al., 2015)—also provide a constraint on the sensitivity of simulated productivity to temperature and precipitation anomalies (Piao et al., 2013).

The ASC in the Northern Hemisphere (averaged over latitudes >45°N, 500 mb), as calculated from aircraft data in two separate 3- to 4-year periods (IGY data; 1958–1961; Keeling, Harris, & Wilkins, 1968) and (HIPPO and NOAA data; 2009–2011; Wofsy, 2011), has been reported to show a large increase of (57 ± 7)% (Graven et al., 2013). Atmospheric transport simulations reveal that these ASC observations are sensitive to net uptake on land north of 40°N, with negligible contributions from lower latitudes or from the ocean (Graven et al., 2013). The increase in ASC in studies examining the multi-model TBM simulations combined with atmospheric transport modeling show that the simulated increase in seasonal amplitude of c_a is generally underpredicted both at high altitude (Graven et al., 2013) and at remote surface stations (Thomas et al., 2016). It has been speculated that this is because TBMs underpredict CO₂ and/or warming effects on GPP (Forkel et al., 2016; Thomas et al., 2016), but changes in the seasonal cycle of ecosystem respiration (e.g., Liptak, Keppel-Aleks, & Lindsay, 2017), including winter respiration (Parazoo et al., 2016), must also be implicated.

The magnitude of the CO₂ fertilization effect on photosynthesis depends on biochemical demand for CO₂ and supply of CO₂ from the atmosphere to the chloroplast via the stomata. Global TBMs typically simulate biochemical demand by C₃ plants (accounting for ~77% of photosynthesis globally; Still, Berry, Collatz, & Defries, 2003) using the photosynthesis model of Farquhar and von Caemmerer (1981), which equates gross assimilation (A) with the lesser of Rubisco and electron transport-limited assimilation rates (A_c and A_j, respectively). The sensitivities of A_c and A_j to c_a are well known (Lloyd & Farquhar, 1996; Long, 1991) and quite different, meaning that the relative contributions of A_c (high sensitivity) and A_j (low sensitivity)-limited rates to A_n are an important determinant of the leaf-level CO₂ fertilization effect. At the leaf level, the dimensionless elasticity of gross carbon assimilation by photosynthesis (A) to c_a characterizes the instantaneous sensitivity of A to c_a. Leaf gas-exchange observations (Farquhar & von Caemmerer, 1981; Maire et al., 2012) support the coordination theory (Chen, Reynolds, Harley, & Tenhunen, 1993; Farquhar & von Caemmerer, 1981; Wang et al., 2017) that plants optimize the productivity in their environment through relative nitrogen investment in electron transport and Rubisco-limited steps in the photosynthesis chain, such that they are co-limiting. Such co-limitation occurs for average growing conditions spanning wide-ranging c_a, vapor pressure deficit, light, and temperature (Maire et al., 2012), with the result that the long-term contributions of A_c and A_j to A are approximately equal. The results of Maire et al. (2012; reproduced in Figure 1a) represent the predictions of A_c and A_j for mean growth conditions of the previous month, based on the experimental values of photosynthetic leaf traits. The response to increasing c_a appears to be most marked, with a reduction in Rubisco content of leaves, in plants growing in nutrient-limited environments (Wong, 1979). This coordination of photosynthesis, and its associated constraint on β, is not generally represented in TBMs, although approaches are emerging that enable such representation (Ali et al., 2016; Haverd et al., 2018; Wang et al., 2017). The usefulness of this constraint on CO₂ fertilization as a driver of global GPP variation has so far been largely neglected.

Here we reconcile the leaf-level constraint on the CO₂ fertilization effect imposed by coordination of photosynthesis with the constraints implied by trends in biospheric activity reflected in globally observed increases in ASC, NLS, and GPP inferred from atmospheric COS. This is achieved using mechanistic principles encapsulated by a TBM framework (Pugh et al., 2016) that links biophysics, ecophysiology, biogeochemistry (including nitrogen cycle), and vegetation structural dynamics via dependencies on climate, c_a, disturbance and land-use drivers, and evolving system state. We quantify the contribution of the biospheric CO₂ response to the trend in global GPP and explore the potential implications for the future evolution of the NLS.

2 MATERIALS AND METHODS

2.4 CABLE: Historic simulations

Global simulations were performed at 0.5 × 0.5° spatial resolution, with time steps of 3 hr (biophysics); 1 day (biogeochemistry); and 1 year (woody demography, disturbance, land use change). The nitrogen cycle was enabled, but not the phosphorus cycle. Recently developed parameterizations for drought response of stomatal conductance and effects of leaf litter on soil evaporation were enabled (Haverd et al., 2016). The relationship between V_c,max,0 and leaf nutrient status was prescribed using the meta-analysis of globally distributed leaf gas-exchange data by Walker et al. (2014). Other photosynthetic parameters were prescribed to be consistent with this analysis, in which any effect of mesophyll conductance (g_m) is implicit in the inferred, effective value of V_cmax,0. This effective value will be lower than the true Rubisco capacity owing to the resistance to gas transport between the stomatal chamber and the chloroplast stroma where carboxylation takes place. Following the recent recommendation of Bahar, Hayes, Scafaro, Atkin, and Evans (2018): “while incorporating g_m into models to estimate GPP is theoretically appealing, in practice it is not yet possible as it requires knowledge of Rubisco kinetic parameters for a range of evergreen tree species, an estimate of g_m and its temperature response for representative species,” we chose not to account for g_m explicitly.

Simulations were driven by (a) daily CRU-NCEP V7 (1901–2016; Viovy, 2009), down-scaled to 3-hourly resolution using a weather generator (Haverd, Raupach, et al., 2013); (b) CO₂ (1-year) resolution (Dlugokencky & Tans, 2017); (c) gridded nitrogen deposition (10-year resolution; Lamarque et al., 2011); and (d) gridded gross land use transitions and harvest (1500–2015) and initial land use states (1,500) from the LUH2 harmonized land use data set (Hurtt et al., 2011; Hurtt, Chini, Frolking, & Sahajpal, 2016).

The “Excl. Coord.” configuration of the model that excluded the coordination of C₃ photosynthesis differed from the full model in three ways: (a) dynamic optimization of the J_max,0/V_cmax,0 ratio was disabled; (b) the J_max,0/V_cmax,0 ratio was set to a low value of 1.3 to ensure that A_J was always less than A_c; (c) V_cmax,0 was scaled up by a factor of 1.25 to ensure that the simulated present day GPP from this configuration and the full model configuration are approximately equal (Figure S3a).

The partitioning of CABLE’s historic global GPP between climate and CO₂ effects (Figure 2a; Figure S5) was achieved by subtracting the results of a “CO₂ only experiment” from an “all driver experiment.” In the CO₂-only experiment, CO₂ and land use are varied and recycled meteorology (1901–1930) is prescribed.

In contrast, the CO₂-driven increase in GPP is derived using the “S1” scenario from the TRENDY protocol, in which fixed 1,860 land cover is assumed and only CO₂ is varied. This additional simulation allows for direct comparison of the simulated CO₂ fertilization effect with the TRENDYv6 ensemble.

We performed further single and double-factorial experiments to isolate the effects of single drivers and driver combinations: c_a only; temperature only; c_a and temperature together; precipitation; land use change.

3 RESULTS

Using the CABLE TBM, enabled to represent the coordination of photosynthesis (Haverd et al., 2018; Figure 1a), we successfully capture four major independent global-scale benchmarks of trends in biospheric activity. They are as follows: (a) multi-decadal increases in global GPP (Figure 1b), (b) amplitude of the seasonal CO₂ cycle in the Northern Hemisphere (Figure 1c), (c) global NLS (Figure 1d) including its mean, trend, and interannual variations that are considered a proxy for tropical temperature sensitivity (Figure S1), and (d) leaf area index (Figure S2). Using model experiments to attribute the effects of changing climate and c_a (Section 2), we find that the total increase in GPP is predominantly driven by CO₂ (Figure 2a). The remaining climate-driven component (14% of the full increase, Figure 2a) is largely attributable to warming and (temperature × c_a) interactions in Northern Hemisphere extratropical ecosystems (Figure 3c; Figure S5).

We estimate a historic global CO₂ fertilization effect on photosynthesis of 30% (1900–2010; 296–389 ppm c_a) that is significantly higher than current TBM (TRENDYv6) estimates (17 ± 4%; 1σ), Figure 2b. We project a GPP increase of (47 ± 2)% for a doubling of c_a (300–600 ppm, constant climate and land use, Section 2). Our corresponding extratropical Northern Hemisphere estimate (>30°N) of (42 ± 2)% (1 SE) is at the high end of the previously reported estimate (32 ± 9)% (2σ) obtained by constraining an ensemble of climate-carbon cycle model projections with observed trends in the ASC of c_a (Wenzel et al., 2016).

The simulation of coordination of photosynthesis contributes strongly to the high-CO₂ fertilization effect in CABLE (Figure 2b). Both CABLE configurations agree well with present-day GPP estimates based on upscaled eddy flux data (Jung et al., 2011) in terms of magnitude and latitudinal distribution (Figure S3b), highlighting that the contrasting CO₂ elasticities (Figure S3a) lead to different magnitudes and latitudinal distributions of the NLS (Figure S3a,c).

We partition the CO₂ fertilization effect on GPP between leaf-level and greening effects, using the dimensionless β diagnostic, which expresses the instantaneous sensitivity of photosynthesis to CO₂. This allows the contribution of greening (i.e., increased light harvesting capacity due to the investment of additional plant carbon in a higher leaf area) to be calculated as the difference between overall GPP increase and the leaf-level contribution to that increase. We find that 70% of the CO₂-driven increase in global GPP is attributable to the leaf-level effect, that is the CO₂-induced increase in GPP that occurs at constant leaf area. The remaining 30% is attributable to the CO₂-greening effect (Figure 2a), that is the positive feedback resulting from the investment by vegetation of a proportion of increased productivity in enhanced leaf area. The leaf-level effect from rising c_a includes both the direct biochemical effect (effect of c_a on Rubisco carboxylase as a substrate and inhibition of ribulose bisphosphate oxygenation and hence photorespiration; Long, 1991) and the indirect effect of enhanced water-use efficiency resulting from the reductions in stomatal conductance that adjusts to maximize carbon gain with respect to water loss (Franks et al., 2013).

Two-thirds of the global CO₂-greening effect is attributable to semiarid ecosystems. Over the (1982–2010) period, the 16% increase in semiarid GPP may be partitioned between leaf-level (6%) and CO₂-greening (10%) components (C), in agreement with the 11% increase in foliage cover reported for the same period in warm semiarid regions where relative change in foliage cover may be equated with the CO₂-greening effect on GPP (Donohue, Roderick, Mcvicar, & Farquhar, 2013). In contrast, the tropical forest GPP (1982–2010) increases proportionately with c_a (both 12%), almost entirely because of the leaf-level effect (Figure 4a), in agreement with the increase (12%) inferred from tropical forest catchment water balance (Yang, Donohue, Mcvicar, Roderick, & Beck, 2016). Globally, the leaf-level CO₂ effect dominates total GPP increase (Figure 2a), meaning that vegetation greenness trends are not generally a direct proxy for GPP trends (see also De Kauwe, Keenan, Medlyn, Prentice, & Terrer, 2016).

To explore the implications of a higher than expected CO₂ fertilization effect for the future fate of the land carbon sink, we employed CABLE, constrained by observed global trends (Figure 1b–d) to project the net land sink (2006–2099) under a low emissions, moderate land use change scenario (RCP2.6). In this scenario, c_a peaks at 440 ppm in 2050, global mean air temperature stabilizes at 2.4° above the 1900 level by 2040, carbon-climate feedbacks are expected to be small (Ciais, Sabine, & Bala, 2013), and emissions from carbon-stored permafrost are negligible (McGuire et al., 2018). Given these assumptions, we find that the terrestrial biosphere would continue to sequester carbon well into the second half of the century, leading to a cumulative net land sink (CO₂, climate, and land use change driven) of 149 PgC by the end of the century, significantly larger than the mean of the CMIP5 carbon-climate model ensemble under the RCP2.6 scenario (65 PgC; Ciais et al., 2013), but well within the among-model spread of the same ensemble (−50 to +195 PgC; Figure 5c). Adjusting to include only biophysical effects (CO₂ and climate-driven sink), we estimate a cumulative biophysical sink of 174 PgC (spatial distribution; Figure 5a,b), equivalent to 17 years of anthropogenic CO₂ emissions at current rates. This sink is large compared to the potential carbon uptake attainable by the recovery of existing secondary forests (117 Pg; Houghton & Nassikas, 2018), and 57 PgC larger than if a lower CO₂ fertilization effect (simulated by CABLE disregarding coordination) is assumed (Figure 4c).

4 DISCUSSION

The magnitude of the simulated CO₂ effect on electron transport-limited photosynthesis (A_j), the increase above this effect due to coordination, and the simulated reduction in β with c_a emerging from our study are supported by other evidence. First, the historical CO₂ fertilization effect on global GPP that is predicted by CABLE Excl. Coord. (A limited by A_j) is comparable to that predicted by the diagnostic photosynthesis algorithm of Wang et al. (2017) that is unconstrained by coordination theory, and explicitly equates A with A_j (Figure 2b). Furthermore, the same algorithm predicts a decline of β of (37%→19%) for a doubling of c_a (400→800 ppm) at 25°C (Keenan et al., 2016). Using CABLE Excl. Coord., we find similar corresponding sensitivities (35%→17%), but much higher values of β (56%→34%) for full CABLE, in which coordination is allowed to occur (see also Figure S3a for latitudinal profiles of β at different c_a). Second, by applying the full diagnostic photosynthesis algorithm of Wang et al. (2017) that is constrained by their alternate implementation of coordination theory, we confirm a large increase in c_a sensitivity of historic global GPP when coordination is simulated, compared with the simulations in which A is limited to A_j (Figure 2b). Third, as noted above, β decreases as c_a increases, and this is due in part to a decline in maximum carboxylation capacity (V_cmax) as relative investment of leaf nitrogen in electron transport increases such that coordination of photosynthesis is maintained. The simulated reduction in V_cmax agrees well with CO₂-acclimation effects on V_cmax inferred from free-air CO₂ enrichment (FACE) studies (~10% reduction for an increase in c_a from 366 ppm to 567 ppm; Ainsworth & Rogers, 2007). We simulate a decrease of V_cmax/J_max of around 20% (400–600 ppm c_a), whereas a synthesis of FACE results reports a decrease of only 5% (Ainsworth & Rogers, 2007). The robustness of this average estimate is however questionable, since it relies on a fixed linear relationship between J_max and V_cmax (Ainsworth & Rogers, 2007). In fact, the proportionality between these quantities probably depends on nutrient status; as Wong (1979) first showed, cotton plants grown under low N (0.6 mmol/L nutrient solution) show a larger reduction of initial slope (V_cmax) than do plants grown with plentiful N (24 mmol/L), but even in the latter case V_cmax/J_max appears to decrease with elevated c_a.

Fourth, our simulated CO₂ responses of GPP and net primary production (NPP; plant growth) in temperate forests globally agree well with the syntheses of data from FACE studies in temperate forests (Ainsworth & Rogers, 2007; Norby et al., 2005; Figure 6). EucFACE provides an exceptional example: here there is a negligible response of aboveground plant growth, probably because of strong phosphorus limitation (Ellsworth et al., 2017). Ecosystem models broadly reproduce the response of ecosystem NPP to elevated CO₂ that is determined experimentally at FACE sites (Hickler et al., 2008; Zaehle et al. 2014). Differences in CO₂ response between CABLE and other ecosystem models are mainly apparent in latitudes south of 15°N (tropics and savanna regions; see Figure S9, showing latitudinal contributions to each TRENDY model's global change in GPP since 1901). Hickler et al. (2008) likewise noted higher simulated responses of NPP to CO₂ in the tropics than seen in temperate forest FACE sites.

Three other models in the TRENDYv6 ensemble simulate a centennial increase (1901–2010) of global GPP of ~30%, similar to that of CABLE and in agreement with the COS observation-based estimate. These provide an opportunity to explore the alternative hypotheses to a higher-than-expected CO₂ fertilization effect as explanations for the high increase in global GPP: (a) LPJ-GUESS simulates very high land use change-induced increase in Northern Hemisphere extratropical GPP (30–60°N), (b) JSBACH and JULES simulate high-land use change and climate-induced increases in GPP in the (−12 to −30°S) latitude band (Figure S9C-D). Land use change is unlikely to induce the increases in GPP in this latitude band because agricultural land area over the last century has been dominated by clearing of tropical forest—a land use conversion that is expected to contribute negatively to trends in NPP (Nyawira, Nabel, Don, Brovkin, & Pongratz, 2016, table 4), and hence GPP. Climate-induced increases in GPP simulated by JULES and JSBACH also appear to be unrealistic, based on apparent temperature and precipitation sensitivities derived from interannual variations in the global NLS (Figure S1b,c), which are known to be dominated by interannual variations in plant productivity in tropical forest and savanna regions (Ahlström et al., 2015). In particular, the apparent temperature sensitivity of JSBACH is positive whereas the atmospheric CO₂ anomalies point to a negative relationship with the temperature (Figure S1b), while JULES exhibits an exaggerated sensitivity to the precipitation (Figure S1c). Thus, the high GPP increases in these models appear to stem partly from simulated climate sensitivities that are too high. Furthermore, the high-relative changes are confounded with significant overestimates of GPP in these same latitude bands, compared with the observation-based FLUXNET estimate (Figure S9a). This suggests that the contributions of land use change and climate alone in other models simulating a large GPP increase are likely to be too high. Our hypothesis of a large CO₂ contribution, particularly in tropical forests and semiarid regions, is more compatible with the spatial distribution of GPP suggested by the FLUXNET data, and with the hydrological constraints (Yang et al., 2016) on GPP trends in tropical forests, as discussed above. Nonetheless, the role of agricultural expansion and intensification on global greening and GPP remains important and unresolved (Mueller et al., 2014; Zeng et al., 2014).

Our simulated latitudinal and regional distributions of trends in GPP and the NLS (Figure 3a,b; Figures S5–S8) are consistent with the observed increase in the ASC of CO₂ in the Northern Hemisphere (56 ± 10% north of 45°N, 1960–2010; Graven et al., 2013), that is underestimated by current TBMs (Figure 1b), as previously reported (Graven et al., 2013; Thomas et al., 2016). Alternate explanations for the observed trend are increasing light-use efficiency because of CO₂ fertilization (Bastos et al., 2019; Piao et al., 2017; Thomas et al., 2016; Wenzel et al., 2016) or high-latitude warming effects on biome distribution and plant productivity (Forkel et al., 2016). Warming effects emerged as important in simulations using the LPJmL TBM which accurately simulates the large observed ASC trend, (Forkel et al., 2016), but simulates a global GPP trend ( PgC/year²; 1970–2011; Forkel et al., 2016) that is less than half the corresponding CABLE-simulated trend (0.53 ± 0.02 PgC/year²; 1970–2011). The CABLE trend is compatible with the COS-based estimate (Campbell et al., 2017), whereas the LPJmL estimate appears too low. LPJmL is thus an example of a model that satisfies the regional ASC constraint, but not the global constraint on the trend in GPP. The high-CO₂ fertilization effect in CABLE, which is dominated by a leaf-level effect in the Northern Hemisphere (Figure 4b,d), supports the hypothesis of increasing light-use efficiency as a key driver of the amplitude of seasonal cycle increase. However, we note that this increase is highly dependent on lags between photosynthesis and ecosystem respiration and cannot be readily disentangled into warming versus c_a-driven trends in GPP.

The uncertain contribution of agricultural intensification to the increase in seasonal amplitude of CO₂ in the Northern Hemisphere has been the subject of much debate (Gray et al., 2014; MacBean & Peylin, 2014; Zeng et al., 2014), but the following four points argue against agricultural intensification being a major contributor to global GPP increase: (a) the seasonal amplitude of CO₂ in the Northern Hemisphere is sensitive to CO₂ uptake and release north of 35°N (Graven et al., 2013) and not the whole globe, and land north of 35°N contributes less than 25% to global GPP (Figure S9). (b) Globally, it has been estimated that the combined effects of land use change and land management over the last century—including expansion of agricultural lands, and intensification by irrigation and fertilization—have had almost no impact on global NPP (Bondeau et al., 2007, figure 15f), although the effects on the seasonal cycle are nontrivial (Calle, Poulter, & Patra, 2019). (c) Expansion of agricultural land area over the last century has been dominated by clearing of tropical forest—a land use conversion that is expected to contribute negatively to trends in NPP (Nyawira et al., 2016). (d) Recent studies using multiple models and data sets have concluded that the effect of land use change and land management on the increase in seasonal amplitude of CO₂ in the Northern Hemisphere is very small compared with the effects of CO₂ fertilization and warming-induced lengthening of growing season (Bastos et al., 2019; Piao et al., 2017). Bastos et al. (2019) identify Eurasia as the region contributing most to increasing seasonal amplitude of CO₂—a region that is dominated by natural ecosystems and has experienced very little land use change (Verburg et al., 2015) over the last half-century. A limitation of the study by Piao et al. (2017) is that it relies on an ensemble of models, many of which do not account for agricultural intensification. Bastos et al. (2019) address this limitation in their more recent study, by complementing the ecosystem model results with trend driver attribution using statistical models, confirming that the observed trends in seasonal amplitude of CO₂ are not significantly related to agricultural intensification.

By reconciling multiple global-scale observational constraints, we identified a CO₂ fertilization effect on historical global GPP that is significantly higher than current estimates. Independent regional studies using amplitude of seasonal cycle data (Northern Hemisphere extra-tropics; Wenzel et al., 2016) and catchment water balance (tropical forests; Yang et al., 2016) have also inferred larger CO₂ fertilization effects than predicted by TBM ensembles. The causes of such model-data discrepancies are poorly known, but biases associated with the representation of nutrient limitations on GPP have been invoked as one possible cause (Wenzel et al., 2016; Yang et al., 2016). Our results, that account for nitrogen-cycle effects on ecosystem productivity, suggest that underprediction of GPP trends and CO₂ responses is associated with a failure by current TBMs to account for plant coordination of photosynthesis. This finding is important for the future role of land carbon sinks, suggesting an underestimate by current models of potential CO₂ removal under low-emission scenarios consistent with the Paris Agreement targets.