Recent Amplified Global Gross Primary Productivity Due to Temperature Increase Is Offset by Reduced Productivity Due to Water Constraints

Satellite remote sensing observations show an increased greenness trend over land in recent decades. While greenness observations can indicate increased productivity, estimation of total annual productivity is highly dependent on vegetation response to climate and environmental conditions. Models have been struggling to determine how much carbon is taken up by plants as a result of increased atmospheric CO2 fertilization. Current remote sensing light use efficiency (LUE) models contain considerable uncertainty due to the lack of spatial and temporal variability in maximum LUE parameter and climate sensitivity defined for global plant functional types (PFTs). We used the optimum LUE (LUEopt) previously derived from the global FLUXNET network to improve estimation of global gross primary productivity (GPP) for the period 1982–2016. Our results indicate increasing GPP in northern latitudes owing to reduced cold temperature constraints on plant growth, thereby suggesting increasing negative carbon‐climate feedback in high latitudes. In the tropics, by contrast, our results indicate an emerging positive climate feedback, mainly due to increasing atmospheric vapor pressure deficit (VPD). Further pervasive VPD increase is likely to continue to reduce global GPP and amplify carbon emissions.


Introduction
Life on Earth is supported by plant photosynthesis through gross primary productivity (GPP), which represents the largest annual carbon flux linked directly linked to environmental conditions and atmospheric CO 2 concentrations (Beer et al., 2010). Human reliance on plant photosynthesis includes GPP allocation to food, fiber, and fuel production, as well as ecosystem services provided by offsetting atmospheric CO 2 emissions from fossil fuel consumption (Norby et al., 2010;Quéré et al., 2018;Schimel et al., 2015). For the past three decades, global satellite remote sensing has provided direct observations of the amount of photosynthetic leaf area (Myneni et al., 2002). These observations serve as primary inputs for satellite-driven diagnostic models of GPP such as light use efficiency (LUE) models (Running et al., 2004). Satellite observations from the Advanced Very High Resolution Radiometer (AVHRR) and the MODerate resolution Image Spectrometer (MODIS) sensors provide consistent global measurements of changes in photosynthetic leaf area starting in June 1981 (Zhu et al., 2013). During this period, the global mean atmospheric CO 2 concentration has increased by 20%, from 340 ppm in 1981 to 407 ppm in 2016 (Etheridge et al., 1996;Keeling et al., 2005). This increase has coincided with widespread increase in leaf area (Zhu et al., 2016) and changes in vegetation phenology, including earlier spring green-up (Cleland et al., 2007;Cong et al., 2013;Zhu et al., 2016). Conversely, anomalous changes in global productivity associated with climate extremes driven by the El Niño-Southern Oscillation (ENSO) have also been observed and modeled Nemani et al., 2003;Zhao & Running, 2010;Zhu et al., 2018). In addition to ENSO, other factors coinciding with water scarcity, high temperatures, and large fires  have significantly impacted the global carbon cycle over the past few decades. Some of these satellite-observed events, including large-scale wind throw, biotic events, pest outbreaks, and deforestation, have significantly impacted global vegetation cover Zscheischler et al., 2013Zscheischler et al., , 2014. Extreme events associated with cold temperature events and heavy rain are also known to impact the global carbon cycle (Zscheischler et al., 2014). However, the current generation of remote sensing driven LUE models has several key limitations that make it difficult to properly estimate long-term GPP trends.
The biome property lookup table approach is a well-known shortcoming of LUE models (Madani et al., 2014;Turner et al., 2002;Wang et al., 2010;Way et al., 2005). In this approach, photosynthetic rate is constrained by biome-specific, predefined thresholds to represent optimum climatic conditions for plant productivity. In addition, the maximum LUE rate, which defines potential GPP, is typically assumed to be temporally constant (e.g., Kimball et al., 2017;Running et al., 2004). Improving these basic GPP model limitations will reduce uncertainty in global GPP estimates and advance the understanding of the terrestrial biosphere response to environmental change and climate extremes.
The variability in CO 2 sources and sinks in natural environments including ocean and land ecosystems is driven by the variability in atmospheric CO 2 accumulation rate . However, estimation of the carbon sources and sinks in land ecosystems remains challenging, where the range of variability in estimated annual GPP and its interannual variability and trend is large among Earth system, LUE, and machine learning-based models (Anav et al., 2015). Even though all of the global models reviewed by Anav et al. (2015) show positive annual GPP trends over the last few decades, there are large discrepancies in the estimated magnitude of GPP trends and interannual variability. Previous studies noted that global ecosystem net primary productivity models that use a satellite data-driven LUE modeling approach show an increasing trend for the period 1982-1999 (Nemani et al., 2003), but this productivity trend diverged after 2000 due to climatic changes, including severe droughts (Zhao & Running, 2010).
Here, we provide a quantitative and mechanistic multidecadal assessment of global GPP trends and anomalies using an enhanced remote sensing LUE model. Our primary goal is to identify the most important factors driving long-term GPP change across key bioclimatic regions. We model global monthly GPP using the third-generation Global Inventory Modeling and Mapping Studies (GIMMS3g) fraction of photosynthetically active radiation (FPAR) record for the period 1982-2016 (Zhu et al., 2013) as a primary model input. By building upon our previous experience (Madani et al., 2014;Madani, Kimball, Jones, et al., 2017;Madani, Kimball, & Running, 2017), our enhanced LUE model provides temporally and spatially explicit dynamic optimum LUE (LUE opt ) information that supports improved estimates of long-term (1982-2016) GPP trends across the globe.

Geospatial Data
We acquired the global semimonthly GIMMS3g FPAR data (Zhu et al., 2013) for the 35-year (1982-2016) study period. GIMMS3g FPAR is created based on the relationship between the new improved GIMMS3g normalized difference vegetation index (NDVI) and best quality MODIS leaf area index (LAI) and FPAR products for the overlapping period 2000-2009 using a neural network algorithm (Zhu et al., 2013). We obtained meteorological and other geospatial information, including monthly minimum air temperature, dew point temperature, incoming shortwave solar radiation, and surface-to-root-zone soil moisture (SM) from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2; Gelaro et al., 2017). We aggregated FPAR semimonthly data to monthly scales by averaging the FPAR values over each month and resampled the MERRA-2 ½°latitude by 5/8°longitude spatial resolution meteorological data using nearest neighboring approach to match the FPAR 0.08°spatial resolution for modeling GPP at monthly time scale for the entire GIMMS3g record. Adopting the finer FPAR resolution for the GPP model simulations, rather than the coarser resolution of the meteorological data, allowed us to better capture the effect of land cover and land use change on GPP (Robinson et al., 2018).

Spatially Explicit LUE opt Data
Global, spatially distributed LUE opt values were derived from the flux tower-based estimates of LUE opt (Madani, Kimball, & Running, 2017). The tower eddy covariance CO 2 flux measurement sites presented in Madani, Kimball, and Running (2017) represent a broad range of global biomes (Table S1 in the supporting information) with at least 2 years of daily ecosystem CO 2 exchange measurement records at each site. In this approach, the upper 98-99.5% bin of FLUXNET tower daily gap-filled GPP values for each tower site is selected to represent the maximum daily GPP (GPP max ). It is assumed that in the upper bin of GPP, plant productivity is not restricted by climate constraint factors (Kergoat et al., 2008;Madani et al., 2014). The FPAR data collocated with FLUXNET tower site locations are temporally matched with the tower GPP records, and PAR is resampled to FPAR resolution using nearest neighboring approach. For each tower site, LUE is defined as follows: In Equation 1, APAR is the absorbed photosynthetically active radiation (PAR), which is derived from the product of FPAR defined from the GIMMS3g record and the daily PAR, estimated as half of the global incoming shortwave solar radiation derived from the MERRA-2 global reanalysis (Gelaro et al., 2017). For each of the tower sites the averaged daily LUE observations from Equation 1 are used to represent the LUE opt value of that specific site.
We extrapolated the tower-based LUE opt values to the global domain based on a generalized additive model (GAM) framework (Hastie & Tibshirani, 1990). The model used several explanatory variables for LUE opt including average annual long-term temperature from MERRA-2 to determine climate sensitivity, the satellite observed solar-induced chlorophyll fluorescence (SIF) from the Global Ozone Monitoring Experiment-2 (GOME-2) (Köhler et al., 2015) to represent biome heterogeneity in productivity within the land cover classifications defined by MODIS MOD12-type-2 (Friedl et al., 2010) classes, maximum and minimum annual FPAR to represent annual changes in land cover as well as the potential effect of the atmospheric CO 2 concentration growth rate on plant leaf area (Zhu et al., 2016). We used the GAM model to provide annual LUE opt information distributed over the global vegetated land areas from 1982-2016 at 8-km spatial resolution. (Refer to Table S2 for the parametric and smoothed coefficient functions of selected environmental predictors used to extrapolate tower estimated LUE opt ). Our model was trained on measurements from a subset of global tower sites from the La Thuile FLUXNET synthesis data set (Baldocchi, 2008) and was tested using independent tower sites from the 2015 FLUXNET record (Pastorello et al., 2020; Refer to Figure S1 for location of tower sites). The trained model was then used along with dynamic annual FPAR observations to generate corresponding spatially explicit LUE opt data from 1982-2016.

Modeling Global GPP
The LUE model used here is similar to the National Aeronautics and Space Administration (NASA) Soil Moisture Active Passive (SMAP) mission's level 4 carbon model algorithm (L4C) , which uses SM as a water supply constraint factor, enabling improved GPP accuracy in water limited regions Kimball et al., 2012;Stocker et al., 2019). Our model also accounts for changes in atmospheric vapor pressure deficit (VPD), which we modeled following Murray (1967): where T a is the average daily temperature in degrees Celsius and T d is the dew point obtained from MERRA2.
Our LUE model provides enhanced GPP estimates (hereinafter termed GPP Enh ) as follows: 10.1029/2020AV000180 where f VPD, fSM, and f T min represent dimensionless environmental constraint functions ranging from 0 (fully constrained) to 1 (no effect) that describe reductions from optimal GPP due to water and temperature stress: The Min and Max subscripts in Equations 4 and 5 represent the minimum and maximum defined thresholds for minimum daily air temperature (T min ) and VPD functions derived from the bioclimatic factors controlling productivity at global scales (Madani, Kimball, & Running, 2017), in addition to ecosystems phenological patterns indicated from flux tower observations. In this regard, we used global tower sites to acquire minimum and maximum thresholds for the T min and VPD bio-climatic variables. T min in the LUE model defines the length of plant activity. We defined T Mmin and T Mmax as 10 and 20 quantiles of the daily GPP climatology and recorded SIF value of the corresponding time for a given tower location. We used a similar technique to establish the VPD thresholds with the exception of using the upper 90 daily GPP quantiles to assess the negative impact of high VPD on stomatal conductance. We then used the observed SIF seasonality to generate spatial maps of environmental constraint factors and used the constraint factors only for regions where seasonality in productivity, confirmed by SIF observations, was shown to be controlled by the specific constraint factor (Madani, Kimball, & Running, 2017).
Daily SM for the global simulations was normalized as a daily proportion of the maximum and minimum reported local SM values from the long-term (1982-2016) record for each pixel. The resulting normalized SM values were aggregated to monthly time steps and used in the following nonlinear constraint function built upon a nonlinear relationship between SM and LUE (Stocker et al., 2018) to estimate GPP in Equation 2:

Validation and Analysis
We validated GPP Enh against flux tower GPP observations for the 2007 to 2014 period obtained from the 2015 FLUXNET record, where the tower validation sites were independent from the sites used for model training ( Figure S1). The tower sites used for model training and validation were selected on the basis of being representative of the major global biomes and having at least 2 years of CO 2 flux measurements. To assess different factors contributing to changes in GPP, we executed the GPP Enh model on a monthly basis with static APAR (APAR climatology) and variable climatic factors ( f VPD × fSM × f T min ) and once again with static climatic factors (climatology of f VPD × fSM × f T min ) and dynamic APAR data. We extended our analysis by detrending GPP and underlying factors controlling productivity, including annual FPAR, SM, PAR, T min , and VPD, and performed annual and multidecadal assessment of GPP anomalies at regional and global scales. We calculated the anomalies in time series data by removing the linear trend. In this regard, residuals e t , or the differences between the data values ( y) and the corresponding linearly fitted values over time x t , are defined as follows: For comparison, we used the GIMMS FPAR data and MERRA-2 meteorology to model GPP using fixed LUE max values used by the MODIS MOD17 operational (Collection 6) GPP product (Zhao & 10.1029/2020AV000180 Running, 2010). We also compared our GPP Enh key findings with the TRENDY ensemble mean GPP, atmospheric CO 2 inversion model results, and SIF observations.
We used the ensemble mean GPP of 10 dynamic global vegetation models (DGVMs ; Table S3) from the TRENDY-v7 project (Quéré et al., 2018;Sitch et al., 2015) for comparison with our GPP model results.
The selected models with spatial resolutions of 0.5°to 2°for the period 1982-2016 use climate, land use, and CO 2 forcing effects on ecosystem productivity. We also acquired net biome productivity (NBP) data from CO 2 inversion model results with 1°monthly spatiotemporal resolution from six inversion models, including CT2017, CTE2018, CarboScope s76_v4.2, CarboScope s85_v4.2, JAMSTEC, and CAMS (see Table S4 for references and details) to compare with our LUE model findings. In addition to modeled ecosystem productivity data, we used SIF from the scanning imaging absorption spectrometer for atmospheric chartography (SCIAMACHY) for 2003-2011 (Joiner et al., 2012) and GOME-2 (Köhler et al., 2015) for 2007-2016 as remote sensing indicators of ecosystem productivity. To mitigate artifacts in the GOME-2 SIF retrievals after mid-2012 due to sensor degradation , we corrected the drift in time series data by matching the mean of observations after mid-2012 to the mean values from 2007 to mid-2012. We generated the anomalies in annual GPP Enh , NBP, TRENDY GPP, and SIF by calculating the departure from long-term average and normalized the values using: where x = (x 1 , …x n ) and z i denotes i th normalized data. We compared these data over tropical and northern high latitudes, the two highly important regions for carbon cycle dynamics.

Results
The new GPP Enh model explains 80% of the variation in annual GPP across flux tower sites, with an RMSE of 331 g C m −2 year −1 . The variance explained declines to 75%, and RMSE increases to 506 g C m −2 year −1 for the model with constant LUE max (and otherwise the same meteorology and FPAR data; Figure S2; Refer to Figure S3 for comparison between GPP Enh with the conventional LUE model over independent test sites). The improvement in the GPP estimate was a result of the environmental explanatory variables that explained~56% of the spatial variation (p < 0.0001) among tower observed LUE opt values. However, the fixed LUE max parameters defined for each land cover type could only explain~36% of the variance in tower observed LUE opt .
The 35-year linear trends in GPP demonstrate that, in~50% of the vegetated land areas, GPP is increasing by up to 20 g C m −2 year −1 , whereas the GPP of tropical regions is declining at the same rate. Black dots represent pixels with statistically significant trends (p < 0.05).
Our estimated global average GPP over the last three decades is 130 ± 1.6 Pg year −1 . The lowest GPP is estimated for 1983 (126 Pg year −1 ), and the highest GPP is for 2011 (133 Pg year −1 ). Over the study period, annual GPP trends indicate that GPP in Amazonia and the Southeast Asian tropics decreased at rates of up to 20 g C m −2 year −1 (at grid scale), while GPP in the northern latitudes increased at the same rate ( Figure 1). To further assess the regional GPP trends, we performed a multidecadal GPP assessment for selected latitudinal zones.
In the northern high latitudes (>45°N) GPP began to increase by 0.07 Pg year −1 from 1982 onward, which represents about 0.4% increase in GPP per year relative to the 35-year mean (17.9 Pg C). In contrast, equatorial GPP (10°S to 10°N) has steadily declined since the 1980s, leading to a reduction of 0.5-1 Pg over 35 years compared to the long-term average (Figure 2).
We further analyzed the main factors affecting global annual GPP anomalies. VPD contributes the highest variability to GPP in tropical regions. In arid environments, water restrictions defined by SM and VPD are the primary limiting factors. In the northern latitudes, where the growing season is short, seasonal cold temperatures primarily limit productivity (Figure 3).
In the Amazon forest, western and central United States, southern Australia, and Africa, GPP limitation is more related to water constraints (VPD and SM) through the negative impact of high VPD and the positive 10.1029/2020AV000180 impact of SM availability. In tropical rainforest, GPP is controlled by the amount of incoming radiation and the negative impact of high VPD (Figure 3).
We analyzed the relationships between GPP anomalies at decadal time scales and underlying environmental factors affecting plant activity by detrending the long-term annual data. During the 1980s, global GPP was most strongly correlated with PAR and FPAR, indicating that the water and temperature constraints had a relatively smaller influence on interannual variability in GPP (Figure 4). Since in the LUE model the cold temperature constraint factor controls the length of the growing season, it exerts the strongest control in seasonally cold environments such as temperate forests and northern high-latitude ecosystems. On the other hand, the sensitivity of productivity to SM, which has a stronger influence on the productivity in arid environments, has not significantly changed in 2000s compared to previous decades. However, with increasing temperature, the cold temperature constraint effect declines, while correlations with VPD, which limits the productivity during the growing season, increased after the 1980s, indicating that global GPP is shifting from being temperature limited to VPD limited (Figures 4 and S4).
Because GPP reductions occur primarily in tropical zones (Figure 2), we performed an anomaly analysis for GPP in the tropics. Horizon plots ( Figure 5) show how annual GPP in each tropical region has changed relative to the average of annual GPP for 1982-2016. GPP significantly declined after the early 2000s in the Amazon and, to a lesser degree, in Africa, whereas GPP in the Asian rainforests began to decline almost a decade earlier than on other continents. GPP in the Amazon has been declining since the 2005 drought, so its mean annual GPP was~0.13 Pg C lower during the 2000s (compared to the 35-year average) and has continued to drop by up to −1.2 Pg C year −1 after the 2010 drought. The annual average GPP of the African tropics was 8.14 Pg C for the study period; it began to slightly decline after 2000 by about 0.06 Pg C and increased by about 0.03 Pg C after 2010. In the Asian tropics after the 1990s, average GPP indicates a decline of 0.3-2% per decade (0.03-0.17 Pg C).
To unravel the underlying mechanism driving tropical GPP change, we performed a VPD anomaly analysis that directly influenced our modeled GPP results for the tropics. Figure 5b demonstrates that VPD in the Amazon began to increase in the early 2000s. In the African tropics, VPD increased from the mid-1990s to mid-2000s, resulting in decreased GPP compared to the 1982-2016 average. In the Asian tropics, where interannual VPD variability is much lower than in Africa and the Amazon, PAR is a larger limiting factor than VPD (Figure 3).
The correlations between interannual GPP variability and environmental factors that constrain our LUE model in the tropics indicate that, in the Amazon, GPP variability is significantly (p < 0.05) correlated with variability in VPD (R 2 = 0.43) and PAR (R 2 = 0.47) but has no significant correlation with FPAR variability over the 35-year record. However, FPAR showed a low but statistically significant correlation (p < 0.05) with GPP after 2005 in the Amazon (R 2 = 0.16). Over the 35-year period, the African tropics showed a significant correlation with VPD (R 2 = 0.21) but not with PAR and FPAR. GPP in the Asian tropics showed a significant We further compared anomalies as a departure from long-term mean values of GPP from the TRENDY models, inversion model CO 2 fluxes, and SIF from the GOME-2 and SCIAMACHY ( Figure 6) for tropical and northern northern middle and high latitudes. Our results indicate that, unlike the TRENDY GPP, the GPP Enh model shows a recent variable response of northern ecosystem productivity to climatic changes. GOME-2 SIF also shows a variable annual signal despite focusing on a shorter period compared to GPP Enh , TRENDY, and net biome production (NBP) data from the inversion models. The NBP data show increasing net CO 2 uptake after 2000, even though there is more interannual variation compared to the TRENDY data.   ) with average detrended FPAR, PAR, SM, T min , and VPD values for the corresponding decade. Interannual GPP before 2000 was highly correlated with FPAR and PAR variations, but global GPP was significantly controlled by higher atmospheric VPD after 2000.

Discussion
Our results indicated increasing trends in annual GPP in middle to high latitudes. The GPP increase shown in the northern tundra and boreal ecosystems (>45°N) supports previous evidence of greening trends observed from long-term satellite records (Myers-Smith et al., 2015;Zhu et al., 2016). Our results also provided evidence of a link to warming and longer growing seasons consistent with recent climate change (Mao et al., 2016;Zhu et al., 2016). Figure 6. Regional anomalies in long-term ecosystem productivity metrics and estimates. Ecosystem productivity metrics include GOME-2 (2007, Joiner et al., 2012) and SCIAMACHY (2003-2011, Köhler et al., 2015 SIF; GPP models included the enhanced GPP (GPP Enh ; this study) and TRENDY GPP (ensemble mean of 10 ecosystem models; Quéré et al., 2018;Sitch et al., 2015), compared with net biome productivity (NBP; ensemble mean from six inversion models, see Table S4 for references) for (a) northern latitudes (> 45 o N) and (b) tropical zones (10°S to 10°N). Anomalies as departures from the mean are calculated at regional scales for each year and normalized for visualization and comparisons.

10.1029/2020AV000180
The rapidly changing arctic and boreal ecosystems are crucial components of the Earth system that store more than 30% of terrestrial carbon stocks (Apps et al., 1993;Pan et al., 2011). While boreal ecosystems have remained a persistent terrestrial carbon sink , recent models and observations predict that increasing air temperatures will reduce the carbon uptake capacity of these biomes over the next century (Liu et al., 2019;Natali et al., 2019). Longer growing seasons and earlier observed photosynthesis from climate warming (Assmann et al., 2019;Box et al., 2019;Parazoo et al., 2018) lead to increased rate and duration of evapotranspiration which can deplete SM and plant available water in the late growing season (Buermann et al., 2013(Buermann et al., , 2018Lian et al., 2019;Parida & Buermann, 2014;Yi et al., 2014;Zhang et al., 2020). Prevalence of warming and browning in the Arctic (Bhatt et al., 2013;Phoenix & Bjerke, 2016;Treharne et al., 2019) also increases the risk of fire occurrences (Hu et al., 2010). All of these factors can affect satellite-observed FPAR and SIF, while our model results also confirm recent high interannual variability in productivity in the northern high latitudes consistent with variability in temperature and water constraints.
Although the GPP increase in the northern high latitudes indicates a persistent, increasing negative carbon-climate feedback, our results suggest an emerging positive feedback to climate in the tropics. The negative GPP trend in the tropics suggests that the increased atmospheric water demand is not balanced by increased available water supply. The changes in rainfall patterns and recent increase in forest mortality in the Amazon forest (Brienen et al., 2015;Phillips et al., 2009;Wigneron et al., 2020) are clear examples of the severe impact of episodic drought and changes in patterns of water supply on these critical ecosystems. These changes in water supply and precipitation forcing in the Amazon influence VPD through land-atmosphere feedback and the trends in PAR.
In contrast to the declining trends seen in GPP Enh in the tropical zones, the TRENDY models show an increase in GPP. LUE models have the advantage over prognostic vegetation models of a direct FPAR observational constraint and can thus potentially reflect anthropogenic effects such as deforestation and human pressure on the tropical hydroclimate system (Khanna et al., 2017) and indirectly impact ecosystem productivity. Like other LUE models, our model is directly constrained by remote sensing observed vegetation indices that have been at the center of debates, especially over dense tropical forests (Bi et al., 2015;Huete et al., 2006;Morton et al., 2014;Saleska et al., 2016). However, our model revealed that there is no significant correlation between interannual variability of GPP in tropical South America with FPAR variability. Instead, we report a strong sensitivity of tropical GPP to VPD variability, which has also been shown for spaceborne SIF (Lee et al., 2013). We analyzed monthly VPD and GPP climatology observed at a CO 2 flux tower site in the Amazon and found a reduction in GPP when VPD increased beyond 800 Pa ( Figure S5). Our results are consistent with other reports of increasing VPD at a global scale after the mid-1990s (Yuan et al., 2019) and highlight the potential constraining impact of increasing water limitations on global ecosystem productivity. This is especially true in the tropics, where changes in water constraints can lead to variable responses in net carbon exchange . However, this VPD impact on productivity seems to be less emphasized in Earth system models (Smith et al., 2016), which show increasing vegetation activity in the tropical zones after 2000 (Figure 7).
In tropical zones, where TRENDY models show increased GPP after 1997, GPP Enh estimates show divergent results, including a reduction in annual GPP after 2004. NBP obtained from inverse models generally indicates enhanced carbon uptake in the tropical zones after 2004 with some variation. It should be noted that inversion models have difficulty modeling the distribution of carbon sources and sinks in the tropics given the intensity of tropical convection, which can affect the spatial distribution of CO 2 concentration (Malhi & Phillips, 2004). In addition, the divergence in productivity estimates between DGVMs and LUE models can be related to DGVM oversensitivity to trends in atmospheric CO 2 fertilization (Smith et al., 2016) including lack of nutrient limitations, as these models tend to have higher sensitivity to CO 2 increase in tropical ecosystems than temperate and boreal ecosystems (Hickler et al., 2008;Schimel et al., 2015). However, it has been argued that most LUE models underestimate the CO 2 fertilization effect, as they do not explicitly account for atmospheric CO 2 concentrations (De Kauwe et al., 2016). LUE models are parametrized using carbon flux towers that have been operational since the late 1990s (Baldocchi et al., 2001). The dynamic effect of CO 2 fertilization on traditional LUE models is only reflected in FPAR observations that show long-term sensitivity to CO 2 trends (Chen et al., 2019;Zhu et al., 2016). Even though we used dynamic maximum and minimum annual FPAR for LUE opt extrapolation to represent changes in the atmospheric CO 2 growth rate and the effects of land use change on photosynthetic efficiency, it is likely that the long-term trend in our LUE based GPP is underestimated.
As we addressed, our study was limited by the LUE opt estimated for the majority of the flux towers that were operational mostly during the 2000s. However, water and temperature constraints also play a significant role in controlling the growing season length and vegetation phenology. Our results indicate that when climate remained static, the APAR-only model, driven implicitly by leaf area and CO 2 fertilization, increased global GPP for the period 1982-2016 at a rate of 0.1 Pg C year −1 ( Figure S6). The addition of climate constraints reduces the global APAR-driven GPP trend by 10% to 0.09 Pg C year −1 . However, it is important to note that the long-term trends here are affected by nonstationarity in time series. For example, large ENSO events affect these trends, but our results indicate that climate warming and drying in the tropics are gradually reducing the GPP growth rate at global scale. This estimated annual GPP trend is significantly lower than the TRENDY estimated GPP trend of 0.57 Pg C year −1 , which optimistically follows the atmospheric CO 2 growth rate pattern ( Figure S7).
Our GPP approach of using spatially and temporally variable LUE opt shows significant improvements over using fixed predefined LUE max values per biome type ( Figures S2 and S3). The LUE opt model is based on the concept that ecosystem processes differ based on plant community compositing and that consideration of the geographic location and key life history traits of plants better accounts for the range of plant functional relationships with climate (Madani et al., 2014). Improving the LUE concept should also lead to better understanding of the response of plant productivity to climate change, despite the limitations associated with our LUE model approach as a whole.
Here, we focused only on the uncertainties related to extrapolated LUE ( Figure S8) that were caused by random errors. At the global scale, these errors correspond to less than 10% of the LUE values for PFTs ( Figure S9) and 6 Pg C standard deviation in annual GPP estimates. The resulting standard deviation around GPP estimates ( Figure S10) does not affect our key findings. Nonetheless, uncertainties are involved in each of the LUE model inputs (Zhao et al., 2005) including the MERRA-2 surface meteorological data. Like all reanalysis data, MERRA-2 estimates may be impacted adversely by discontinuities in the assimilated satellite observing system record that impact the modeled water and energy fluxes (Gelaro et al., 2017;Robertson et al., 2016). Moreover, the use of gauge-based precipitation forcing in MERRA-2 can likewise result in discontinuities, especially in poorly observed regions such as the Amazon  their Figure 8).
Even though CO 2 fertilization and nutrient effects are indirectly considered in the remote sensing-derived FPAR observations and the spatially explicit estimation of LUE opt model, future work should more directly account for these effects. Further improvement in the LUE models by including higher spatiotemporal resolution meteorological information capturing local variations in SM (due to topography) and incoming shortwave radiation (due to clouds, diffuse and direct fraction), better representation of disturbance events such wildfires, and full representation of plant water availability, such as the inclusion of surface-to-groundwater information and the assimilation of satellite data (Madani et al., 2020;Smith et al., 2019), may further improve the model correspondence with productivity benchmark observations derived from the satellite SIF and global carbon flux tower record. These improvements will enable more accurate assessments and attribution of long-term climate and CO 2 effects and improved benchmarking of DGVMs, giving us better insight into future productivity changes.

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.