Higher than expected CO2 fertilization inferred from leaf to global observations

Abstract Several lines of evidence point to an increase in the activity of the terrestrial biosphere over recent decades, impacting the global net land carbon sink (NLS) and its control on the growth of atmospheric carbon dioxide (ca). Global terrestrial gross primary production (GPP)—the rate of carbon fixation by photosynthesis—is estimated to have risen by (31 ± 5)% since 1900, but the relative contributions of different putative drivers to this increase are not well known. Here we identify the rising atmospheric CO2 concentration as the dominant driver. We reconcile leaf‐level and global atmospheric constraints on trends in modeled biospheric activity to reveal a global CO2 fertilization effect on photosynthesis of 30% since 1900, or 47% for a doubling of ca above the pre‐industrial level. Our historic value is nearly twice as high as current estimates (17 ± 4)% that do not use the full range of available constraints. Consequently, under a future low‐emission scenario, we project a land carbon sink (174 PgC, 2006–2099) that is 57 PgC larger than if a lower CO2 fertilization effect comparable with current estimates is assumed. These findings suggest a larger beneficial role of the land carbon sink in modulating future excess anthropogenic CO2 consistent with the target of the Paris Agreement to stay below 2°C warming, and underscore the importance of preserving terrestrial carbon sinks.

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 sinksincluding 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 2 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% . While there is consensus that the simulated increase is largely driven by CO 2 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; and has been increasing at a rate of 0.061 ± 0.02 (1 SE) PgC/year 2 (1980Le 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-1961Keeling, Harris, & Wilkins, 1968) and (HIPPO and NOAA data;2009-2011Wofsy, 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 2 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 2 fertilization effect on photosynthesis depends on biochemical demand for CO 2 and supply of CO 2 from the atmosphere to the chloroplast via the stomata. Global TBMs typically simulate biochemical demand by C 3 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 2 fertilization effect. At the leaf level, the dimensionless elasticity of gross carbon assimilation by photosynthesis (A) to c a = A c a c a 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 2 fertilization as a driver of global GPP variation has so far been largely neglected.
Here we reconcile the leaf-level constraint on the CO 2 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 2 response to the trend in global GPP and explore the potential implications for the future evolution of the NLS.

| Community Atmosphere-Biosphere Land Exchange
The Community Atmosphere-Biosphere Land Exchange model (CABLE) is a biogeochemical land surface model (including N-cycle) that represents vegetation structural dynamics. CABLE consists of a biophysics core (Haverd, Cuntz, Nieradzik, & Harman, 2016;Kowalczyk et al., 2013;Wang et al., 2011), the CASA-CNP biogeochemistry module (Wang, Law, & Pak, 2010), the POP module for woody demography and disturbance-mediated landscape F I G U R E 1 Constraints on trends in global-scale terrestrial biospheric activity. (a) Coordination of photosynthesis. Black: CABLE predictions of A c versus A j contributions to C 3 photosynthesis  at 1,000 randomly distributed grid cells across the global ice-free land surface. Red lines reproduced from Maire et al. (2012): A c versus A j under average plant growth conditions predicted from observed leaf-level traits (solid); predictive error (dashed). (b) Increase in global gross primary production relative to 1900, simulated by CABLE (black) and TRENDYv6 (Le Quéré et al., 2018;Sitch et al., 2015 grey), and inferred from observations of [COS] (4) (red); Increase in c a (purple; Dlugokencky & Tans, 2017). (c) Increase in amplitude of the seasonal cycle (ASC) of c a (1958-1961 to 2009-2011, >45°N, 500 mb) versus ASC of net land carbon sink (NLS) north of 30°N: aircraft observations (Keeling et al., 1968;Wofsy, 2011) and corresponding TM3 (Heimann & Körner, 2005) atmospheric inversion reproduced from Graven et al. (2013;red); CABLE-TM3 (black); TRENDYv6-TM3 (grey); average for each model (large circles); individual ensemble members (small circles); ensemble range (bars; Section 2). (d) Time-series of annual global NLS: CABLE (black); TRENDYv6 (grey); Global Carbon Project (GCP) mass balance estimate (note 0.4 PgC/year uncertainty), as the sum of atmosphere and ocean sinks, minus fossil fuel emissions (red). Linear trends : 0.061 ± 0.02 (1 SE) PgC/year 2 (GCP) and 0.067 ± 0.01 PgC/year 2 (CABLE). See also Figure S1 for metrics related to the NLS heterogeneity (Haverd, Smith, et al., 2013), and a module for land use and land management (POPLUC; Haverd et al., 2018). The land use and land cover change module is driven by gross land use transitions and wood harvest area. CABLE is one of an ensemble of ecosystem and land surface models contributing to the Global Carbon Project's annual update of the global carbon budget (Le Quéré et al., 2018) via the TRENDY ensemble.
This article uses a version of CABLE that was recently extended to account for (a) tree demography and disturbance-mediated landscape heterogeneity; (b) land use and land cover change; and (c) coordinated leaf nitrogen investment in electron transport and Rubisco-limited steps in the photosynthesis chain, such that they adjust seasonally and across biomes to be co-limiting (Haverd et al., 2018). This version also includes a new diagnostic elasticity = A c s c s A that allows the GPP increase attributable to the leaf-level component of the CO 2 fertilization effect to be isolated as A c s Δc s , where c s is the CO 2 concentration at the leaf surface.

| CABLE: Coordination of C 3 photosynthesis
CABLE represents the coordination of C 3 photosynthesis, adjusting electron transport and Rubisco-limited rates (Farquhar, Caemmerer, & Berry, 1980) to be co-limiting (Haverd et al., 2018). This is achieved by dynamically optimizing the ratio b JV of the maximum rate of electron transport (J max ) to the Rubisco-catalyzed maximum rate of carboxylation (V cmax ) at standard temperature. The optimization aims to maximize the rate of photosynthesis per unit leaf nitrogen.
The approach to optimization of b JV is based on the following three assumptions: 1. Leaf nitrogen resources may be dynamically redistributed at a 5-day timescale at no cost, that is b JV is optimized, such that net photosynthesis (given total available leaf nitrogen) accumulated over the last 5 days (approximately the timescale for turnover of photosynthetic nitrogen; Thornley, 2004) would have been maximized.

Leaf nitrogen resources available for partitioning between
Rubisco and electron transport capacity are proportional to effective nitrogen content (N eff ), defined as the sum of prior estimates of V cmax,0 and J cmax,0 , weighted by relative cost factor c cost,JV .
(There are no explicit nitrogen pools corresponding to Rubisco and electron transport capacities because there are no explicit nitrogen pools in the photosynthesis model of Farquhar et al. 1980.) 3. The prior values of V cmax,0 (related to time-varying prognostic leaf nitrogen content) and b JV are prescribed according to the synthesis of globally distributed leaf gas-exchange measurements by Walker et al. (2014).
The method for implementing these assumptions in CABLE is: 1. Maintain a 5-day history of sub-diurnal leaf-level meteorology (absorbed PAR; leaf-air vapor pressure difference; leaf temperature, CO 2 concentration at the leaf surface) for sunlit and shaded leaves, such that total photosynthesis over the last 5 days can be reconstructed.
2. Construct a function that calculates leaf nitrogen cost per unit net photosynthesis. Inputs to this function are as follows: (a) current estimate of b JV , (b) effective leaf nitrogen content (related to prior value of V cmax,0 ), (c) 5-day history of sub-diurnal leaf-level meteorology.
3. At the end of each day, solve for b JV that maximizes the total photosynthesis over the last 5 days per unit effective leaf nitrogen, and update V c,max,0 and J max,0 values.
The emerging contributions of electron transport and Rubisco-limited rates contribute approximately equally to total net photosynthesis.

| Ensemble of current TBMs
We use simulations from the TRENDYv6 ensemble of TBMs to represent current TBMs: these models (Table S1)

| 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 . 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.
The "Excl. Coord." configuration of the model that excluded the coordination of C 3 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 2 effects (Figure 2a; Figure S5) was achieved by subtracting the results of a "CO 2 only experiment" from an "all driver experiment." In the CO 2 -only experiment, CO 2 and land use are varied and recycled meteorology  is prescribed.
In contrast, the CO 2 -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 2 is varied. This additional simulation allows for direct comparison of the simulated CO 2 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.

| CABLE: Future projections
The Representative Concentration Pathways (RCPs) were designed to inform the mitigation and adaptation policies, and describe dif- To assess the CO 2 fertilization effect on GPP for a doubling of c a , we also performed a projection (2006-2099) using CO 2 forcing from RCP8.5, combined with fixed IPSL-CM5A-LR climate (repeated 1990-2010) and land cover (2015). The simulated trend in GPP over this period corresponds to projected CO 2 fertilization effects of (47 ± 2)% for a doubling of c a (300-600 ppm, constant climate and land use) for the globe and (42 ± 2)%, for the Northern Hemisphere extra-tropics >30°N.

| Amplitude of CO 2 seasonal cycle in the Northern Hemisphere
The increase in amplitude of seasonal cycle in c a was simulated by using air-land fluxes from the CABLE and TRENDYv6 TBMs as F I G U R E 2 CO 2 fertilization effect on global gross primary production (GPP). (a) Full CABLE simulation with dynamic land use (solid black lines) attributed to: (i) climate effect and c a × climate interaction effect (orange); (ii) CO 2 greening effect (light green); (iii) CO 2 leaf-level effect (dark green). (b) CO 2 -driven increase in global GPP relative to 1900 (fixed 1,860 land use), as simulated by: (i) full CABLE (thick black line); (ii) CABLE without coordination of photosynthesis, that is assimilation is always limited by electron transport (thin black line); (iii) the TRENDYv6 ensemble (grey lines) and (iv) full diagnostic photosynthesis algorithm of Wang et al. (2017), applied using prescribed vegetation cover from CABLE (thick purple dashed line); (v) diagnostic photosynthesis algorithm of Wang et al. (2017) without coordination (i.e., assimilation is always limited by electron transport), applied using prescribed vegetation cover from CABLE (thin purple-dashed line) boundary conditions for "fine grid" version of the global atmospheric tracer transport model TM3 (Heimann & Körner, 2005), with an approximate horizontal resolution of 3.8 × 5° and 19 vertical levels. It was driven by winds from NCEP reanalysis (Kalnay et al., 1996). Influences from outside of the terrestrial biosphere on seasonal CO 2 amplitude in the NH, including fossil fuel emissions and air-ocean fluxes have been shown to be negligible (Graven et al., 2013), and were omitted in this analysis. The simulated monthly time series corresponding to each 3-to 4-year aircraft campaigns was interpolated to the 500 mb isobar and linearly detrended. The seasonal CO 2 amplitude was quantified by taking the difference of the maximum and minimum monthly detrended concentrations.
The NLS in the Northern Hemisphere is the net sum of two similarly large flux components (GPP and ecosytem carbon loss-largely ecosystem respiration, but also including fire loss-to the atmosphere: which both have summer seasonal peaks. Poor predictions of the ASC magnitude may result from small errors in the relative phase and/or the amplitude of either flux, which in turn may influence the simulated ASC trends. To address this source of uncertainty, we made the assumption that the seasonal amplitude in GPP is well-simulated (as confirmed by comparison with latitudinally averaged upscaled FluxNET data (Jung et al., 2011)), but that the relative phase of R ′ eco may vary by ±14 days and the amplitude by ±20% compared with the simulated values. Such errors in the phase and amplitude of simulated R ′ eco may arise from a number of factors such as: (a) approximating soil carbon stocks as pools with discrete turnover times and associated effective soil depths (Koven et al., 2013) Wenzel, Cox, Eyring, and Friedlingstein (2016). Observed monthly mean in situ atmospheric CO 2 concentrations at BRW and KMK were obtained from the Scripps Institute of Oceanography records (http:// scrip psco2.ucsd.edu/data/atmos pheric_co2; Keeling et al., 2001). The ensemble members corresponding to each TBM were selected such that perturbed ASC (2010) at both stations matched the actual ASC to within ±3 ppm. These selected subsets were then used as input to TM3 to simulate the relative change in ASC, with mean and uncertainties for each model expressed as the mean and range of ensemble members (Figure 1c).  Figure S5).

| RE SULTS
We estimate a historic global CO 2 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 F I G U R E 3 Spatial and driver attribution of gross primary production (GPP) trend . (a) Spatial distribution. (b) Latitudinal distribution of driver components. (c) Driver attribution of global GPP trend. Driver attribution is achieved by factorial experiments in which only one or a pair of drivers is varied while all others are fixed (Section 2). Tot = total effect; C = c a only; T = temperature; C × T = c atemperature interaction; P = precipitation; LUC = land use change. Only major drivers are isolated. Additional drivers (nitrogen deposition; specific humidity; wind speed; incoming short and longwave radiation) and their interactions are not isolated but may contribute to the total effect. Error bars represent standard errors from linear regression fits to each driver component (Section 2) We partition the CO 2 fertilization effect on GPP between leaflevel and greening effects, using the dimensionless β diagnostic, which expresses the instantaneous sensitivity of photosynthesis to CO 2 . 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 2 -driven increase in global GPP is attributable to the leaf-level effect, that is the CO 2 -induced increase in GPP that occurs at constant leaf area. The remaining 30% is attributable to the CO 2 -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 2 -greening effect is attributable to semiarid ecosystems. Over the  period, the 16% increase in semiarid GPP may be partitioned between leaf-level (6%) and CO 2 -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 2 -greening effect on GPP (Donohue, Roderick, Mcvicar, & Farquhar, 2013). In contrast, the tropical forest GPP (1982GPP ( -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 leaflevel CO 2 effect dominates total GPP increase (Figure 2a To explore the implications of a higher than expected CO 2 fertilization effect for the future fate of the land carbon sink, we em- Given these assumptions, we find that the terrestrial biosphere F I G U R E 4 CO 2 fertilization effect on gross primary production (GPP) by region. Regions correspond to land cover classification ( Figure S3). Solid black line: full CABLE simulation with all drivers varying ("all driver expt"). Dotted black line: leaf-level CO 2 effect as a component of the all driver expt. Solid green line: CABLE simulation in which CO 2 is varying, actual land cover assumed and recycled meteorology  is prescribed ("CO 2 only expt"). Dotted green line: leaf-level CO 2 effect as a component of the CO 2 only experiment. The difference between the green-dotted and solid lines represents the CO 2 greening effect. The difference between the solid green and black line represents the climate effect and c a × climate interaction effect on GPP

| D ISCUSS I ON
The magnitude of the simulated CO 2 effect on electron transportlimited 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 2 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 . 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 2 -acclimation effects on V cmax inferred from free-air CO 2 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 2 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 2 that is determined experimentally at FACE sites (Hickler et al., 2008;Zaehle et al. 2014). Differences in CO 2 response between CABLE and other ecosystem models are mainly apparent in latitudes south of 15°N (tropics and savanna regions; see Figure S9,  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 2 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 (56 ± 10% north of 45°N, 1960-2010Graven 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 2 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 2 ;1970-2011Forkel et al., 2016) that is less than half the corresponding CABLE-simulated trend (0.53 ± 0.02 PgC/year 2 ; 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 2 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 2 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 2 in the Northern Hemisphere is sensitive to CO 2 uptake and release north of 35°N F I G U R E 6 CO 2 fertilization effect on gross primary production (GPP; a) and net primary production (NPP; b). The solid lines represent the median predictions of CABLE (fixed land use and climate, rising CO 2 following the trajectory of RCP8.5) for each of temperate, tropical and boreal forests, referenced to a 400 ppm baseline. Shading represents 25th and 75th percentiles. Data points (triangles) are taken from meta-analyses by Ainsworth and Rogers (2007; (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 fertilizationhave 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 2 in the Northern Hemisphere is very small compared with the effects of CO 2 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 2 -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 2 are not significantly related to agricultural intensification.
By reconciling multiple global-scale observational constraints, we identified a CO 2 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 2 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 2 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 2 removal under low-emission scenarios consistent with the Paris Agreement targets.