Recent trends in gross primary production and their drivers: analysis and modelling at flux-site and global scales

Gross primary production (GPP) by terrestrial ecosystems is the largest flux in the global carbon cycle, and its continuing increase in response to environmental changes is key to land ecosystems’ capacity to offset anthropogenic CO2 emissions. However, the CO2- and climate-sensitivities of GPP vary among models. We applied the ‘P model’—a parameter-sparse and extensively tested light use efficiency (LUE) model, driven by CO2, climate and remotely sensed greenness data—at 29 sites with multi-year eddy-covariance flux measurements. Observed (both positive and negative) GPP trends at these sites were predicted, albeit with some bias. Increasing LUE (due to rising atmospheric CO2 concentration) and green vegetation cover were the primary controls of modelled GPP trends across sites. Global GPP simulated by the same model increased by 0.46 ± 0.09 Pg C yr–2 during 1982–2016. This increase falls in the mid-range rate of simulated increase by the TRENDY v8 ensemble of state-of-the-art ecosystem models. The modelled LUE increase during 1900–2013 was 15%, similar to a published estimate based on deuterium isotopomers. Rising CO2 was the largest contributor to the modelled GPP increase. Greening, which may in part be caused by rising CO2, ranked second but dominated the modelled GPP change over large areas, including semi-arid vegetation on all continents. Warming caused a small net reduction in modelled global GPP, but dominated the modelled GPP increase in high northern latitudes. These findings strengthen the evidence that rising LUE due to rising CO2 level and increased green vegetation cover (fAPAR) are the main causes of increasing GPP, and thereby, the terrestrial carbon sink.


Introduction
Terrestrial gross primary production (GPP), the total ecosystem carbon uptake by photosynthesis, is the largest flux in the global carbon cycle (Friedlingstein et al 2019). Increasing GPP, primarily due to enhanced photosynthetic CO 2 uptake, is proposed to have driven the uptake of about one-third of anthropogenic CO 2 emissions on average during recent decades (Keenan et al 2014). But there is no consensus on the importance of different controls of spatial and temporal variations in GPP. Ecosystem models continue to differ in the relative effects of different environmental controls; better understanding of these controls is essential if carbon cycle dynamics are to be reliably projected into the future (Ciais et al 2013).
Photosynthesis in C 3 plant leaves increases with increasing CO 2 , the so-called CO 2 fertilization effect (Hungate et al 2003, Schimel et al 2015. However, climatic factors including solar radiation, temperature and soil moisture availability, and soil nutrient availability, also contribute to determining GPP (Terrer et al 2016, Xin et al 2016, von Buttlar et al 2018, Yuan et al 2019, Bartsch et al 2020. Since direct measurements of GPP at global scale do not exist (Ma et al 2015), global GPP can only be inferred indirectly, either using models, or based on independent proxies ( Models used to estimate spatial and temporal patterns of GPP fall into two broad categories. Data-driven models based on eddy-covariance measurements at flux tower sites rely on fitted relationships between inferred GPP and explanatory variables at site level, followed by upscaling with the help of remote-sensing measurements (Jung et al 2009, Tramontana et al 2016. Data-driven models can also be based on remote sensing of optical parameters linked to vegetation activity (Damm et al 2010, Zhao and Running 2010, Zhang et al 2016. The limitations of data-driven models are that they rely on data availability-this is a problem due to the highly clumped distribution of flux towers across the world, and their relatively short period of operation to date-and on the correct specification of explanatory variables (Beer et al 2010). In particular, the data-driven models illustrated here do not account for the effect of atmospheric CO 2 concentration on the light use efficiency (LUE) of GPP (De Kauwe et al 2016). Process-based models simulate the exchanges between different ecosystem components in a more explicitly mechanistic way (Friedlingstein et al 2006, Sitch et al 2008, with GPP represented by the Farquhar et al (1980) photosynthesis model or some variant of it. Limitations of process-based models include their complexity, number of (often poorly known) parameters, and different embedded hypotheses-not always clearly stated-that lead to differing outcomes (Anav et al 2013(Anav et al , 2015Prentice et al 2015).
Here we adopt a third approach, which aims to avoid some of the more troublesome limitations of both data-driven and complex process-based models. The 'P model' (Stocker et al 2019) is a new, firstprinciples model for GPP. It uses remotely-sensed greenness data as input (in common with datadriven models), but relies on a universal, processbased formulation of LUE. This formulation is based on the Farquhar et al (1980) photosynthesis model, combined with additional, individually tested ecoevolutionary optimality hypotheses (Prentice et al 2014;Wang et al 2017a, Smith et al 2019, Franklin et al 2020 that represent the acclimation and/or adaptation of photosynthetic capacities and stomatal behaviour to the growth environment. The P model achieves at least as good performance as other LUE models (e.g. Yuan et al 2007, Hilker et al 2008 in simulating GPP, as assessed by comparison with eddy-covariance measurements-and does so with far fewer parameters than other models (Wang et al 2017a, Stocker et al 2019. Our objectives were (1) to evaluate the P model's ability to reproduce historical trends in GPP at different flux-measurement sites; (2) to attribute these trends to environmental drivers; (3) to simulate the historical time-course of global GPP (for comparison with complex process-based model simulations, and with estimates based on independent proxies); and finally (4) to use the P model to assess the relative contributions of different drivers to the recent historical increase in global GPP.

Model description
The P model calculates GPP for C 3 photosynthesis as follows (see Wang et al 2017a, Stocker et al 2020: where I abs is the light absorbed by green plant tissues, expressed as the product of photosynthetic photon flux density (PPFD, µmol m -2 s -1 ) and fraction of absorbed photosynthetically active radiation (fAPAR, equivalent to the percentage of green vegetation cover), and LUE (= GPP/I abs ) is the product of three dimensionless quantities: the intrinsic quantum yield of C 3 photosynthesis (φ 0C3 ), a term (m) representing the CO 2 limitation of electron-transport limited photosynthesis, and the square-root term, which reduces GPP due to the cost of maintaining electron-transport capacity, expressed by the constant cost factor c * . The model embodies the coordination hypothesis: that is, carboxylation capacity is assumed to adjust to match the electron-transport limited rate of photosynthesis, which in turn is constrained by the absorbed PPFD. The CO 2 limitation term is given by: where c a is the ambient partial pressure of CO 2 (Pa), Γ * is the photorespiratory compensation point (Pa), D 0 is the leaf-to-air vapour pressure deficit (vpd, Pa), K is the effective Michaelis-Menten coefficient of Rubisco (Pa), η * is the ratio of viscosity of water to its value at 25 • C, and β is the ratio of carboxylation and transpiration cost factors at 25 • C. This expression combines the standard expression for electron-transport limited photosynthesis (Farquhar et al 1980) with the 'least-cost' formulation of the optimal ratio of leaf-internal to ambient CO 2 (Prentice et al 2014). A simple modification, presumed adequate within the c a range of study, is adopted for C 4 plants: for these we assume m = 1 (accounting for their CO 2 concentrating mechanism) and allow a different (lower) intrinsic quantum yield (φ 0C4 ): The temperature-dependence of φ 0 (Singsaas et al 2001) was modelled as follows: where T c is leaf temperature ( • C). Equation (4) is from Bernacchi et al (2003); equation (5) was fitted to measurements in Kubien et al (2003). φ 0C3 was set to zero for T c ≥ 0 • C and φ 0C4 was set to zero whenever equation (5) yielded negative values.
Rubisco kinetic parameters (K, Γ * ) and their temperature dependencies were obtained from the in vivo measurements on tobacco leaves by Bernacchi et al (2003). K and Γ * also depend on atmospheric pressure (Farquhar et al 1980, Wang et al 2017b. Values of β and c * were estimated from independent data by Wang et al (2017a): β from a global compilation of leaf stable carbon isotope measurements, and c * from published measurements of carboxylation and electron-transport capacities. For simplicity, leaf temperatures were equated to air temperatures.
Additional effects of low soil moisture were included via an empirical function that reduces light use efficiency (fLUE) from optimal values, following Stocker et al (2018). This function, β(θ), is: where θ is plant-available soil water as a fraction of available water-holding capacity, and θ * = 0.9 is the threshold value above which the function is held constant at 1. The sensitivity (q) of β(θ) to soil moisture is related to a moisture index, an estimate of the climatological ratio of actual to potential evapotranspiration (α) calculated by the SPLASH (Davis et al 2017) model: with a = 0.179 and b = 0.450 for grasslands, and a = 0.101 and b = 0.0063 otherwise, based on the global analysis of GPP data by Stocker et al (2018). The eddy-covariance flux tower sites were chosen based on their relatively long and continuous records, homogeneous vegetation cover in their surroundings. GPP for these locations was simulated for the period 1999-2014 using gap-filled CO 2 , temperature, shortwave radiation, vpd, pressure and precipitation data downloaded from the FLUXNET2015 (Tier 1 data) dataset (https://fluxnet.fluxdata.org/data/fluxnet2015dataset/). To avoid using poor-quality data, we only selected daily and monthly data flagged as high quality (quality flag ≤0.75) and gap-filled using linear interpolation (Zhang et al 2010, Verma et al 2015). The fAPAR data were downloaded from the Copernicus Global Land Service (https:// land.copernicus.eu/global/products/fapar) at a resolution of 0.5 × 0.5 degree. The mean value of 9 × 9 grid pixels surrounding the flux site was assigned as the site fAPAR value. The P model was run at a monthly timestep; SPLASH runs at a daily timestep.
Globally, the P model was run for the period 1982-2016 driven by monthly CO 2 data from Mauna Loa Laboratory (www. esrl.noaa.gov/gmd/ccgg/trends/data.html), vapour pressure and temperature from CRU TS 4.01 (Harris et al 2014), shortwave radiation from WFDEI (Weedon et al 2014), and greenness data from the GIMMS fAPAR 3 g dataset (Zhu et al 2013). Daily shortwave radiation, precipitation and air temperature were downloaded from WFDEI to run the SPLASH model.

In situ GPP
Monthly GPP data using the nighttime partitioning method (Reichstein et al 2005) and constant USTAR thresholds (GPP_NT_CUT_MEAN) were extracted from FLUXNET2015 (Pastorello et al 2020) and aggregated to annual totals. Twenty-nine sites were used (table 1) with record lengths ranging from eight to 16 years.

Model description 2.3.1. Process-based global GPP models
We compared our global results with simulations from 15 dynamic global vegetation models (DGVMs) in the TRENDY version 8 ensemble (table 2) . Three of these experiments (S0, S1, S2) were used. S0 is a baseline simulation with no change in forcing; S1 specifies only observed CO 2 changes; S2 specifies observed changes of CO 2 and climate. Monthly simulated GPP values over the period 1982-2016 were extracted and aggregated to annual global totals.   (2017) used long-term atmospheric COS trends derived from measurements in ice cores, firn and ambient air to infer changes in GPP.

Trends at flux sites
We calculated the temporal trend at each site, both in data and in simulations, during the studied period using the Theil-Sen slope estimator (mblm package in R) (Komsta and Komsta 2013). By choosing the median of the slopes through pairs of points, this regression method minimizes the influence of outliers and is robust against heteroscedasticity, especially for small samples.
As each site (both simulations and observations) has a different GPP trend and different uncertainties in the trend, we then performed 'linear regression with errors in X and Y' (the York_fit function in Matlab) (Wiens 2020) following York's method . This method performs a linear regression between two sets of data taking account of uncertainties that vary across the data points. We used this method to test the model's ability to represent the different trend slopes, across the whole set of sites.

Global trend comparisons
The spatial pattern of GPP trends over 1982-2016 was computed for each grid cell using the Theil-Sen slope estimator. To compare the global GPP trend derived from the P model with those from TRENDY, FLUXCOM and MTE, we first aggregated GPP from each grid cell weighted by grid-cell area to obtain global GPP, and calculated the trend using the Theil-Sen slope estimator. As FLUXCOM and MTE cover a shorter period, we calculated their trends separately over 2001-2015 and 1982-2011 so that combined they covered the study periods of TRENDY and the P model.

Contributions of different forcings
At site level, we performed experimental simulations to assess the relative contribution of different drivers: CO 2 , temperature, vpd, soil moisture, PPFD and fAPAR. For each experimental simulation (S driver ), only the driver to be assessed was allowed to change; other forcings were held at their baseline levels. The GPP trend estimated from that simulation was considered as the contribution of that driver to changes in GPP. A full model run (S all ), in which all drivers changed over time, was conducted to simulate for the overall response of GPP to environmental changes. We then calculated the average contribution of each driver over all study sites to obtain the overall contribution of the different forcings.
Simulations were also carried out at the global scale, and the contributions of CO 2 and climate factors were estimated to allow comparison with the DGVMs. The P model simulation in which  all variables were held constant is denoted S -. We calculated differences of GPP between pairs of simulations, (S CO2 −S -) and (S all −S CO2 ), in order to isolate the effects of CO 2 and climate respectively. The contributions of CO 2 and climate in the DGVMs were estimated analogously as trends in the differences (S 1 −S 0 ) and (S 2 −S 1 ). We further quantified the contribution of each driver to GPP by running simulations in which only one driver was allowed to change, and calculated the differences (S driver −S nought ). These differences were also used to quantify the sensitivity of GPP to different drivers, as follows: where ∆ driver represents the change of each driver between 1982 and 2016, and ε is the residual error term. A multiple regression was carried out at each grid cell in order to define the spatial pattern of the 'most influential' forcing. We fitted the model in each grid cell using all-subsets regression, and chose the model with the lowest AIC value to define the environmental factors that accounted for the change of GPP for that grid (the significant predictor). Then using Anova, the amount of variance explained by each significant predictor can be estimated. The predictor that explained greatest variance was selected as most influential forcing. The area occupied by each most influential forcing was calculated to produce the percentage of land area dominated by these forcings. The trend of LUE over the period 1982-2016 derived from TRENDY models and P model was calculated. LUE derived from the P model was computed by equation (11), while LUE values derived from TRENDY DGVMs were estimated using model-specific simulated LAI (converted to fAPAR using Beer's law with extinction coefficint k = 0.5) and solar radiation data. Global annual total LUE was aggregated from each grid cell weighted by grid-cell area, and the trend of LUE was computed using Theil-Sen slope estimator.

Comparison with independent estimates
We calculated the change in LUE over 1900-2013 to allow comparison with LUE trends inferred from the change in photorespiration/photosynthesis ratio (Ehlers et al 2015). This ratio of photorespiration  to photosynthesis (α * ) is given by the following expression: Assuming that: and expressing m as:

In situ GPP
Despite large uncertainties in the estimated trends of both observed and modelled GPP at individual sites, the P model succeeded in capturing the pattern of variation in GPP trends across sites. These observed trends were positive in most cases, but negative in others (figure 1). The regression slope is somewhat larger than 1 while the intercept is indistinguishable from 0, indicating that the model represents the trends at the different sites with no offset and only a modest bias.

GPP trend modelling at global scale
Mean annual simulated GPP in the TRENDY model ensemble ranged from 115 Pg C yr -1 to 190 Pg C yr -1 and the GPP trend from 0.1 Pg C yr -2 to 0.8 Pg C yr -2 (figure 2). The P model simulated a GPP value towards the higher end of the TRENDY range and a trend of 0.46 Pg C yr -2 , similar to the average of the TRENDY models. MTE showed trends close to zero. FLUXCOM, however, presented a negative trend, possibly due to the relatively short period covered compared with other three datasets. Most regions of the globe showed positive modelled GPP trends in the TRENDY ensemble (figure S2 (available online at https://stacks. iop.org/ERL/15/124050/mmedia)). Large areas with increased modelled GPP over northern high latitude indicate the positive impact of lengthening growing seasons (Piao et al 2007). Nevertheless, some regions showed negative modelled GPP trends. All models showed decreased GPP in NE Asia and SW North and South America, probably due to decreased precipitation or intensified dry seasons (Zhao and Running 2010). The spatial pattern of GPP trend produced by the P model is generally similar to that generated by other DGVMs. However, the P model does not show the negative trend over Australia indicated by some other models.

Modelled contributions of environmental drivers
Among flux sites, CO 2 was the only explanatory variable that was consistently and positively associated with the observed trend in GPP (figure 3). Variations (positive or negative) in fAPAR and PPFD also contributed to the observed trends. The GPP sensitivity for CO 2 (i.e. the ratio of the modelled GPP trend to the observed trend in CO 2 ) was estimated as 18.6 ± 13.3 gC m -2 yr -1 ppm -1 , significant but larger than the statistically derived estimate of 4.5 ± 0.75 g C m -2 yr -1 ppm -1 given by Fernández-Martínez et al (2017). The overall sensitivity of GPP to temperature was significant (-193.3 ± 25.8 gC m -2 yr -1 K -1 ), similar to −137.4 ± 415.4 g C m -2 yr -1 K -1 in the analysis by Fernández-Martínez et al (2017).
According to the 25% reduction in the ratio of photorespiration to photosynthesis decribed in Ehlers et al (2015), the estimated change of LUE by P model was 15% during 1900-2013, roughly consistent with 12% derived from Ehlers et al (2015). Campbell et al (2017) estimated an increase of GPP by 31% ± 5% over the 20th century based on COS. We cannot make a direct comparison with this estimate, however, because of the lack of remotely-sensed fAPAR data to drive the P model during most of the 20th century. LUE simulated by the P model over the period 1982-2016 fell in the high end of the TRENDY LUE range (figure 4), which varied between around 0.1 and 0.35, and was generally higher than average LUE derived from TRENDY DGVMs. Nervertheless, the P model presented a similar consistent positive trend and increase of LUE compared with TRENDY models.
The impacts of CO 2 and climate on global GPP in the TRENDY ensemble varied among models (figure 5). Contributions of CO 2 ranged from 0.22 Pg C yr -2 to 0.52 Pg C yr -2 , while contributions of climate varied from −0.05 Pg C yr -2 to0.19 Pg C yr -2 . The P model indicated a contribution of 0.33 Pg C yr -2 from CO 2 and 0.14 Pg C yr -2 from climate, lower than the TRENDY model average contribution of 0.39 Pg C yr -2 from CO 2 and higher than the average of 0.09 Pg C yr -2 from climate, but well within the range of the TRENDY models. Global GPP increased by 17.8 ± 0.8 Pg C (P model) and 21.7 ± 0.7 Pg C (TRENDY models) with a 100 ppm rise of atmospheric CO 2 .
A breakdown of the contribution of each environmental variable to modelled global GPP, and to the total modelled GPP of the northern and southern extratropics and tropics, is shown in figure 6 and tables S2-S3. CO 2 was the dominant contributor to the modelled GPP increase except in northern extratropical latitudes, where fAPAR surpassed CO 2 . With every degree of increase in temperature, global GPP was estimated to decrease by 1.5 ± 0.5 Pg C. The modelled impact of drought (vpd and soil moisture stress) was small compared to other factors. Nonethless, the modelled vpd-induced decrease in GPP offset the positive impact of other factors by −8.25 ± 1.9 Pg C kPa -1 . Figure 7 illustrates the spatial distribution of the most influential forcing over 1982-2016. fAPAR and CO 2 were the two most important controls of modelled GPP. Increasing CO 2 was the most important contributor to increasing modelled GPP over most temperate and tropical forest regions, but increasing fAPAR was most important over large areas of semi-arid vegetation in Asia, the Americas, Africa and Australia. In part of the tropical forest area, increasing PPFD was the dominant control. A positive effect of increasing soil moisture was modelled for scattered regions. Rising temperatures were a dominant cause of increasing modelled GPP in northern high latitudes; only very restricted areas of arid vegetation in South America and Australia showed a dominant negative effect of temperature. Only a few scattered areas showed a dominant effect (positive or negative) of vpd. Taken the area of each grid into account, fAPAR dominated highest land area (12.2%), followed by CO 2 (7.5%). The positive impact of soil moisture and temperature ranked third and fourth respectively, covering 2.2% and 1.6% of land area.

Discussion
The P model as used here is process-based (with an additional empirical component to represent soil moisture effects), but unlike the comprehensive DGVMs included in TRENDY and other ensembles, the model (a) requires only a very small set of parameters, all of which are estimated from independent data, and (b) is anchored in remotely-sensed observations of green vegetation cover-which is a key control on GPP, and yet one of the more uncertain quantities simulated by DGVMs (e.g. Kelley et al 2013). Moreover, we have shown that the P model can produce nearly unbiased simulations of different trends in GPP across flux measurement sites. The duration of this comparison is short, and the set of suitable sites is small and not globally representative. Nonetheless, the success of the comparison is encouraging, and provides the opportunity to examine the contributions of different environmental factors to the trends. Although the set of predictors used in this analysis is different from that used by Fernández-Martínez et al (2017) in their statistical analysis of GPP trends, there are commonalities in the results. In particular, both our results and theirs have proved that despite different trends (and dominant forcing factors) at each site, CO 2 is the one variable that consistently contributes positively to GPP across sites. In contrast, the impact of temperature varies greatly among sites. Increasing temperatures in temperate climates are often associated with lengthened growing seasons, which promote plant yield; but this positive effect can be offset by intensified water stress.
Our simulation of global GPP falls within the TRENDY range of modelled values for both GPP and its historical trend. Both the TRENDY models and the P model however differ from the data-driven FLUX-COM and MTE models, in showing a substantial GPP increase. The lack of such an increase in these datadriven models is, presumably, because they include no mechanism for LUE to increase in response to rising CO 2 (Anav et al 2015). Our simulated increase in LUE during 2000-2013 is consistent with the independent estimate based on deuterium isotopomers by Ehlers et al (2015). Furthermore, our 0.88% increase of LUE during 1982-2016 derived from P model simulations is similar to the 0.93% increase derived from TRENDY DGVMs.
Direct comparison with the COS-based estimate of a 20th century GPP increase of 26%-36% by Campbell et al (2017) is not possible, but this estimate would be consistent with our modelled LUE increase if increasing fAPAR made a positive contribution to GPP of similar magnitude to the contribution of increasing CO 2 . Our analysis suggests that CO 2 has exerted a stronger impact than climate change on global GPP during the 20th and 21st centuries. This finding is consistent with a range of~17%-30% in other studies (Keenan et al 2016, Haverd et al 2020. GPP increase driven by increased LUE due to rising CO 2 accounted for 83% of total modelled global GPP increase. However, regionally increased greenness and a changing climate are also important factors. The dominance of a positive temperature effect in northern high latitudes is consistent with the strong limiting effect of short growing seasons (Piao et al 2009), while globally, the P model indicates a net negative effect of warming on GPP. With over a quarter of the global land area sensitive to soil water availability, the projected temperature increase could also lead to a further soil water shortage, which could aggravate vpd-induced decreases in GPP (He et al 2016). Wang et al (2018) suggested that the weakening control of carbon uptake in spring and the emerging negative temperature response in summer during recent decades has led to a diminished positive impact of increasing temperature over high latitudes; once drought-induced water stress prevails, it is suggested that there will be a turning point towards a declining global carbon sink (Williams et al 2014). However, future projections of carbon cycle changes remain highly uncertain so long as there are large differences in the projections by different models. The P model as applied here represents a first step towards more robust quantification of the carbon cycle sensitivity to global environmental change. Required future developments include developing a wellfounded ability to model changes in green vegetation cover-a key source of uncertainty in contemporary DGVMs.

Acknowledgments
We thank TRENDY v8 modellers, namely, V Haverd and E Kato (VISIT) for contributing model outputs. We give special thanks to S Sitch and P Friedlingstein for organizing TRENDY, and advice on this work. This work used eddy covariance data acquired and shared by the FLUXNET community, including these networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux-TERN, TCOS-Siberia, and USCCC. The FLUXNET eddy covariance data processing and harmonization was carried out by the ICOS Ecosystem Thematic Center, AmeriFlux Management Project and Fluxdata project of FLUXNET, with the support of CDIAC, and the OzFlux, ChinaFlux and Asi-aFlux offices. This research was supported by the scholarship granted by the China Scholarship Council (CSC, 201908060024). This research is a contribution to the Imperial College initiative on Grand Challenges in Ecosystems and the Environment. It has received funding from the European Research Council (ERC) under the European Union's Horizon