Increased atmospheric vapor pressure deficit reduces global vegetation growth.

Atmospheric vapor pressure deficit (VPD) is a critical variable in determining plant photosynthesis. Synthesis of four global climate datasets reveals a sharp increase of VPD after the late 1990s. In response, the vegetation greening trend indicated by a satellite-derived vegetation index (GIMMS3g), which was evident before the late 1990s, was subsequently stalled or reversed. Terrestrial gross primary production derived from two satellite-based models (revised EC-LUE and MODIS) exhibits persistent and widespread decreases after the late 1990s due to increased VPD, which offset the positive CO2 fertilization effect. Six Earth system models have consistently projected continuous increases of VPD throughout the current century. Our results highlight that the impacts of VPD on vegetation growth should be adequately considered to assess ecosystem responses to future climate conditions.

INTRODUCTION

Vapor pressure deficit (VPD), which describes the difference between the water vapor pressure at saturation and the actual water vapor pressure for a given temperature, is an important driver of atmospheric water demand for plants (). Rising air temperature increases saturated water vapor pressure at a rate of approximately 7%/K according to the Clauius-Clapeyron relationship, which will drive an increase in VPD if the actual atmospheric water vapor content does not increase by exactly the same amount as saturated vapor pressure (SVP). Numerous studies have indicated substantial changes of relative humidity (ratio of actual water vapor pressure to saturated water vapor pressure) not only in continental areas located far from oceanic humidity () but also in humid regions (). Although the long-term trend of globally averaged land surface relative humidity remains insignificant (, ), a sharp decrease has been observed since 2000 (, ), implying a sharp increase in land surface VPD. However, the causes of changing atmospheric water demand are still unclear ().

Changes of VPD are important for terrestrial ecosystem structure and function. Leaf and canopy photosynthetic rates decline when atmospheric VPD increases due to stomatal closure (). A recent study highlighted that increases in VPD rather than changes in precipitation substantially influenced vegetation productivity (). Increasing VPD notably affects vegetation growth (–), forest mortality (), and maize yields (). In addition, rising VPD greatly limits land evapotranspiration in many biomes by altering the behavior of plant stomata (). Given that the global precipitation is projected to remain steady (), the changing VPD and soil drying would likely constrain plant carbon uptake and water use in terrestrial ecosystems (). However, the large-scale constraints of VPD changes on vegetation growth have not yet been quantified. In this study, we determined the changes in VPD trends through observation-based global climate datasets, and then quantified the impacts of these VPD changes on vegetation growth and productivity, using satellite-based vegetation index [i.e., normalized difference vegetation index (NDVI)] and leaf area index (LAI), tree-ring width chronologies, and remotely sensed estimates of gross primary production (GPP).

RESULTS

This study used four observation-based globally gridded climate datasets—CRU (Climatic Research Unit), ERA-Interim, HadISDH, and MERRA (Modern-Era Retrospective analysis for Research and Applications) (table S1)—to analyze the long-term trend of VPD over vegetated land. Similar to previous analyses (, , ), anomalies in all four datasets showed that VPD trends were temporally and spatially heterogeneous over recent decades (Fig. 1). A piecewise linear regression method was used to quantify the change in trends and detect the potential turning point (TP) in each dataset. It was observed that VPD increased slightly before the late 1990s but increased more strongly afterward with 1.66 to 17 times larger trends according to the four datasets (fig. S1). The datasets showed that 53 to 64% of vegetated areas experienced increased VPD trends since the late 1990s (fig. S2). To illustrate the magnitude and spatial variability of VPD change, we calculated the global pattern of the percentage change of annual growing season mean VPD between two periods of 1982–1986 and 2011–2015 (fig. S3A). On average, the annual growing season mean VPD of 2011–2015 was 11.26% higher than that of 1982–1986, and the VPD increased larger than 5% in more than 53% area. In addition, the increases of global mean VPD over 12 months were positively correlated with the mean VPD values of 1982–1986 at more than 64.5% areas (fig. S3B), which implies that the higher VPD increases in the months with high VPD.

Fig. 1

Global mean vapor pressure deficit (VPD) anomalies of vegetated area over the growing season.

Anomalies are relative to the mean of 1982–2015 when data from all datasets are available. Vegetation areas were determined using the MODIS land cover product. Blue line and gray area illustrate the mean and SD of VPD simulated by six CMIP5 models under the RCP4.5 scenario.

Apart from HadISDH, datasets showed that the increased saturated water vapor pressure and decreased actual water vapor pressure jointly determined the increases of VPD after the TP. On average, the rate of increase in saturated water vapor was 1.43 to 1.64 times higher after the TP year than before, and the actual water vapor exhibited stalled or decreased trends (fig. S4). Increased air temperature explains the changes in saturated water vapor pressure (fig. S4). The HadISDH dataset indicates a decrease in saturated water vapor because of large spatial gaps in the dataset.

A change of oceanic evaporation is the most important mechanism for the observed decrease in actual water vapor pressure over the land (). Oceanic evaporation is the most important source of atmosphere water vapor, and approximately 85% of atmospheric water vapor is evaporated from oceans, with the remaining 15% coming from evaporation and transpiration over land (). Most of the moisture over land is transported from the oceans, which accounts for 35% of precipitation and 55% of evapotranspiration over land (). We analyzed long-term changes of oceanic evaporation based on a global oceanic evaporation dataset [Objectively Analyzed Air–Sea Fluxes (OAFlux)] (). The almost 60-year time series showed that the decadal change of global oceanic evaporation (Eocean) was marked by a distinct transition from an upward to a downward trend around 1998 (Fig. 2A). The global oceanic Eocean has decreased by approximately 2.08 mm year−1, from a peak of 1197 mm year−1 in 1998 to a low 1166 mm year−1 in 2015 (Fig. 2A), and 76% of the sea surface revealed a decreased Eocean after 1999 (Fig. 2B). Rhein et al. () reported stalled increases of sea surface temperature after the late 1990s based on multiple global datasets, which substantially limited oceanic evaporation (). Some studies using global climate models (GCMs) also highlighted that VPD trends over land were predominantly explained by dynamic mechanisms related to moisture supply from oceanic source regions (, ). Changes in the recycling of atmospheric moisture over land controlled by soil moisture in supply-limited regions may be an additional contribution to the observed increase of VPD. Koster et al. () showed that moisture variability contributed to total precipitation variance in mid-northern latitude regions such as the western United States. Drier soils evaporate less and thus lead to lower water vapor in the atmosphere (). Previous study reported a decreased trend in the global land evapotranspiration after the late 1990s limited by soil moisture supply ().

Fig. 2

Comparison of oceanic evaporation (Eocean) trends during the two periods of 1957–1998 and 1999–2015.

(A) Time series of globally averaged oceanic evaporation. (B) Spatial pattern on differences of oceanic evaporation trends between 1999–2015 and 1957–1998. Gray shaded area in (A) indicates ±1 SD. The inset in (B) shows the frequency distributions of the corresponding differences.

Figure 3 illustrates that the satellite-based NDVI substantially increased from 1982 to 1998 (y = 0.0014x − 1.86, R2 = 0.43, P < 0.05), while NDVI remained constant and then stalled after 1999 (y = −0.0004x + 1.23, R2 = 0.06, P = 0.65) (Fig. 3A). From 1982 to 1998, approximately 84% of the vegetation surface showed an increased NDVI trend (28.50% with a significant increase; Fig. 4A). In comparison, after 1999, the trends of NDVI over many regions reversed, and 59% of vegetation areas showed a pronounced NDVI browning (decreasing) trend (21.50% with a significant decrease; Figs. 3B and 4). Mean NDVI trends for 12 months after 1999 were lower than those from 1982 to 1998 over globally vegetated areas (Fig. 3C). Moreover, we analyzed long-term trends of LAI based on four global LAI datasets [Global Land Surface Satellite (GLASS), GLOMap, LAI3g, and Terrestrial Climate Data Record (TCDR); table S1] (). Despite the large variability of the estimated interannual LAI among the four products, all four LAI datasets exhibited a transition from increasing trends before the late 1990s to decreasing trends afterward (fig. S5). The LAI showed a decreasing trend since the late 1990s over vegetated areas of 64.72, 72.62, 62.73, and 80.11% for GLASS, GLOMap, LAI3g, and TCDR datasets, respectively (fig. S6). The differences of NDVI and LAI trends during these two periods are the opposite of VPD trends derived from four VPD datasets.

Fig. 3

Comparisons of NDVI trends over the globally vegetated areas from 1982 to 2015.

(A) Time series of NDVI. The numbers show the change rates of NDVI, and * indicates the significant changes at a significance level of P < 0.05. (B) Probability density function of NDVI trends during the two periods, with bars indicating the proportion of increased (gray) and decreased (black) responses. (C) Mean monthly NDVI trends between the two periods. Shaded area in (A) and error bars in (C) indicate ±1 SD.

Fig. 4

Comparison of NDVI trends over the globally vegetated areas between two periods of 1982–1998 and 1999–2015.

(A) NDVI trend of 1982–1998. (B) NDVI trend of 1999–2015. (C) Differences of NDVI trend between 1999–2015 and 1982–1998. The insets (I) show the relative frequency (%) distribution of significant decreases (Dec*; P < 0.05), decreases (Dec), increases (Inc), and significant increases (Inc*), and the insets (II) show the frequency distributions of the corresponding ranges.

Partial correlation analysis indicated significant correlations of detrended VPD with detrended NDVI and LAI when the impacts of air temperature, radiation, and atmospheric CO2 concentration were excluded (Fig. 5). Detrended NDVI over 62% of the vegetated areas shows a negative correlation with detrended VPD (about 14% with a significant negative correlation) (Fig. 5A). Similarly, four detrended satellite-based LAI correlated negatively with detrended VPD over 65 to 70% of vegetated areas (16 to 22% with a significant negative correlation) (Fig. 5, B to E). In addition, all five satellite-based datasets show highly consistent signs of correlation with VPD, and at least three datasets revealed consistently negative correlations with VPD over 72% of vegetated area (Fig. 5F). A machine learning method [i.e., random forest (RF)] was used to reconstruct NDVI based on atmospheric [CO2] concentration and five climate factors (air temperature, precipitation, radiation, wind speed, and VPD) over the last 34 years in each pixel (fig. S7) and then model experiments were applied to separate the impacts of VPD as well as of other variables (see Materials and Methods). Globally, the model experiments suggest that the atmospheric CO2 concentration, air temperature, and VPD are the most important contributors for the variability of NDVI (fig. S8A). Rising VPD was found to significantly decrease NDVI, indicated by the larger negative NDVI differences from 1999 to 2015, suggesting that substantial increases of VPD strongly limited NDVI (fig. S8B).

Fig. 5

Spatial patterns of correlations between VPD and satellite-based NDVI/LAI.

Partial correlations between detrended CRU VPD and detrended satellite-based NDVI/LAI were shown: GIMMS NDVI (A), GLASS LAI (B), GLOBMap LAI (C), LAI3g LAI (D), and TCDR LAI (E) during 1982–2015 (GLOBMap and LAI3g from 1982–2011). The insets in (A) to (E) show the relative frequency (%) distribution of significant negative correlations (Neg*; P < 0.05; dark green), negative correlations (Neg; light green), positive correlations (Pos; light red), and significant positive correlations (Pos*; P < 0.05; dark red). (F) Number of satellite-based NDVI/LAI datasets with the same sign of correlation: e.g., (5, –) indicates that all five satellite-based NDVI/LAI datasets showed negative correlations with VPD.

This study used two satellite-based models [revised eddy covariance–light use efficiency (EC-LUE) and Moderate Resolution Imaging Spectroradiometer (MODIS)] to investigate the impacts of VPD on long-term changes of global GPP (, ). EC-LUE and MODIS showed quite similar long-term trends of GPP, with a significantly increased trend from 1982 to the late 1990s, averaged at 0.73 Pg C year−1 (P < 0.05; from 1982 to 1998) and 0.26 Pg C year−1 (P < 0.05; from 1982 to 1997) over globally vegetated area, respectively (Fig. 6A). The GPP trends then stalled and decreased afterward (−0.016 Pg C year−1, P = 0.67 and −0.032 Pg C year−1, P = 0.44) (Fig. 6A). The GPP trends derived from the two models during the two periods are the opposite of VPD trends derived from the four VPD datasets.

Fig. 6

Long-term changes of global GPP and environmental regulations.

(A) Time series of global GPP estimates derived from EC-LUE and MODIS-GPP models. (B) GPP sensitivity to climate variables, NDVI/fPAR, and atmospheric CO2 concentration. (C) Contributions of climate variables, NDVI/fPAR, and atmospheric CO2 concentration to GPP changes over the two periods. Three climate variables are included: vapor pressure deficit (VPD), air temperature (Ta), and photosynthetically active radiation (PAR).

To quantify the impacts of VPD on GPP, we further explored GPP sensitivity to climate variables (i.e., air temperature, VPD, and radiation), atmospheric CO2 concentration, and satellite-based NDVI/fPAR (see Materials and Methods; Fig. 6B). Two satellite-based models showed the similar GPP sensitivity to VPD, whereby global GPP decreased by 13.82 ± 3.12 Pg C and 18.29 ± 3.65 Pg C with a VPD increase of 0.1 kPa (Fig. 6B), which is comparable to the GPP increase with a 100–parts per million (ppm) rise of atmospheric [CO2] (i.e., βCO2 = 19.01 ± 4.01 Pg C 100 ppm−1). On the basis of the estimated GPP sensitivity, we estimated the contributions of climate variables, CO2 fertilization, and vegetation index to global GPP over the two study periods (table S2). After the late 1990s, VPD increased by 0.0017 ± 0.0001 kPa year−1 according to the CRU dataset (fig. S1), which resulted in GPP decreases of 0.23 ± 0.09 Pg C year−1 and 0.31 ± 0.11 Pg C year−1 according to the EC-LUE and MODIS models, respectively (Fig. 6C and table S2). The VPD-induced GPP decreases partly counteract the CO2 fertilization effect (0.38 ± 0.08 Pg C year−1) after the late 1990s with the rising rate of atmospheric CO2 concentration by 2.02 ± 0.01 ppm year−1. From 1982 to the late 1990s, CO2 fertilization played a dominant role in the GPP increase (Fig. 6C). According to the EC-LUE model, GPP increases of 0.28 ± 0.15 Pg C year−1 occurred because of the rising atmospheric [CO2] (Fig. 6C and table S2).

We further investigated the impacts of VPD on LUE using measurements of global EC towers (see Materials and Methods; table S3). Because VPD correlated strongly with air temperature, we excluded the effects of air temperature to investigate the impacts of VPD by binning the observations for ranges of air temperature. We binned the LUE data for different restricted ranges of air temperature and found strong negative correlations between VPD and LUE for almost all air temperature ranges (fig. S9 and table S3). In addition, we investigated the impacts of VPD on vegetation growth using a global comprehensive dataset of tree-ring width measurements from 171 locations with temporal coverage from 1982 until at least 2005. Partial correlation analysis showed that the detrended VPD derived from the four datasets correlated with detrended tree-ring width at most sites (56 to 72%) when the impacts of air temperature, radiation, and atmospheric CO2 concentration were excluded (fig. S10, A to D). We compared the differences of mean tree-ring width values between before and after 1998 and observed smaller tree-ring widths after 1998 compared to those before 1998 at 64% of sites (25% sites with significant level) (fig. S10E).

DISCUSSION

Our results support increased VPD being part of the drivers of the widespread drought-related forest mortality over the past decades, which has been observed in multiple biomes and on all vegetated continents (, ). Increased VPD may trigger stomatal closure to avoid excess water loss due to the high evaporative demand of the air (), leading to a negative carbon balance that depletes carbohydrate reserves and results in tissue-level carbohydrate starvation (). In addition, reduced soil water supply coupled with high evaporative demand causes xylem conduits and the rhizosphere to cavitate (become air-filled), stopping the flow of water, desiccating plant tissues, and leading to plant death (). Previous studies reported that increased VPD explained 82% of the warm season drought stress in the southwestern United States, which correlated to changes of forest productivity and mortality (). In addition, enhanced VPD limits tree growth even before soil moisture begins to be limiting (, ).

We examined whether terrestrial ecosystem models can adequately capture the observed responses of vegetation growth to increased VPD after the late 1990s from 10 terrestrial ecosystem models. We found that the simulated GPP trends of most models did not match the GPP trends documented above (fig. S11). Only the CLASS (Canadian Land Surface Scheme) model showed a decreased GPP after the late 1990s in response to increased VPD, similar to satellite-based GPP estimates (fig. S11). The terrestrial ecosystem models showed lower GPP sensitivity to VPD than two satellite-based models (i.e., EC-LUE and MODIS) (Fig. 6B and table S2).

Our results imply that most terrestrial ecosystem models cannot capture vegetation responses to VPD. Thus, problems reproducing the observed long-term vegetation responses to climate variability may challenge their ability to predict the future evolution of the carbon cycle. Earth system models (ESMs) participating in the CMIP5 (Coupled Model Intercomparison Project Phase 5) (table S5) project a continuous increase of VPD until the end of this century (Fig. 1). The globally averaged VPD is 0.12 kPa higher in 2090–2100 than in 1980–1999 (Fig. 1). The ESMs used in Fig. 1 showed good performance when reproducing historical variations of VPD (table S6), providing confidence in the projected increases of VPD during future decades. The results of our analysis suggest that this projected increased VPD might have a substantially negative impact on vegetation, which must be examined carefully when evaluating future carbon cycle responses.

MATERIALS AND METHODS

Datasets

Four global climate datasets were used to investigate the long-term changes of atmospheric VPD, including CRU, ERA-Interim, HadISDH, and MERRA. Monthly gridded CRU and HadISDH datasets were based on climate observations from global meteorological stations (, ). ERA-Interim and MERRA datasets were reanalysis products based on Integrated Forecast System of European Centre for Medium-Range Weather Forecasts (ECMWF-IFS) () and the Goddard Earth Observing System Data Assimilation System Version 5 (GEOS-5) (), respectively. VPD was calculated on the basis of different variables of four datasets ()

CRU:

ERA-Interim:

HadISDH and MERRA:where SVP and AVP are saturated vapor pressure and actual vapor pressure (kPa), respectively. Td is the dew point temperature (°C). RH is the land relative humidity (%).where Ta is the land air temperature (°C). Z is the altitude (m). Pmst is the air pressure (hPa), and Pmsl is the air pressure at mean sea level (1013.25 hPa). In addition, the OAFlux dataset was used to examine the variability of oceanic evaporation (table S1) ().

We used the newest release of the advanced very high resolution radiometer (AVHRR) NDVI to indicate vegetation growth from 1982 to 2015. The AVHRR is a nonstationary NDVI version 3 dataset made available by NASA’s Global Inventory Modeling and Monitoring Study third-generation dataset (GIMMS3g) group (). GIMMS3g contains global NDVI observations at approximately 8-km spatial resolution and bimonthly temporal resolution, derived from AVHRR channels 1 and 2, corresponding to red (0.58 to 0.68 μm) and infrared (0.73 to 1.1 μm) wavelengths, respectively. Each 15-day data value is the result of maximum value compositing, a process that aims to minimize the influence of atmospheric contamination from aerosols and clouds. Moreover, this study analyzed long-term trends of LAI based on four global satellite LAI products (table S1): GLASS (version 4) (), GLOBMap (), LAI3g (), and the TCDR ().

We calculated the annual growing season mean NDVI and LAI by averaging monthly NDVI and LAI values with monthly mean temperatures above 0°C. We also calculated multiyear averaged monthly mean temperatures from the CRU dataset to ensure that the same growing season land mask was used over the entire period (1982–2015). The global mean NDVI and LAI values were calculated by the average of the annual growing season mean NDVI and LAI, excluding unvegetated regions. The MODIS land cover type product (MCD12Q1) was used to identify the vegetated regions.

We calculated the LUE (g C m−2 MJ−1) based on EC measurements from the FLUXNET2015 dataset (www.fluxdata.org) to examine the correlation between LUE and VPD (table S4)where GPP indicates the estimated GPP values from EC measurements, PAR is photosynthetically active radiation (MJ m−2), and fPAR is the fraction of PAR absorbed by the vegetation canopy calculated by GIMMS3g NDVI ().

The tree-ring width measurements around the world were used from the International Tree-Ring Data Bank (ITRDB) (). The wood samples were taken and processed following standard protocols and taking two radial cores per tree at 1.3 m. Tree-ring width measurements were detrended and standardized by the scientists who contributed the chronologies to the ITRDB. Each local chronology represents the average growth of several trees (typically more than 10) of the same species growing at the same site. The temporal span of the tree-ring data series selected began at 1982, lasting at least until 2005. Eventually, 171 sites were analyzed and each chronology of the sites is a representation of annual tree-ring width.

This study conducts the partial correlation analysis between VPD and tree-ring width by excluding the impacts of air temperature, radiation, and atmospheric CO2 concentration. Air temperature and PAR from MERRA dataset were used. For atmospheric CO2 concentration, this study used the GLOBALVIEW-CO2 product, which provides observations of atmospheric CO2 concentration at 7-day intervals over 313 global air-sampling sites (). If missing 7-day data accounted for >20% of all data for an entire year, then the value for that year was indicated as “missing.” For a site to be included in this study, it had to have at least 10 years of observations. Eventually, 77 sites were included equally in the calculation of global monthly mean CO2 concentration without any weighting of individual sites.

Satellite-based GPP model

We used two satellite-based GPP models to investigate the impacts of VPD on vegetation GPP. The first model is the MODIS-GPP model, and this study used long-term global MODIS GPP dataset driven by GIMMS fPAR data ().

The second model is the revised EC-LUE model (), derived by (i) integrating the impact of atmospheric CO2 concentration on GPP and (ii) adding the limit of VPD to GPP. The revised EC-LUE model simulates terrestrial ecosystem GPP aswhere PAR is the incident photosynthetically active radiation (MJ m−2) per time period (e.g., day); fPAR is the fraction of PAR absorbed by the vegetation canopy calculated by the GIMMS3g NDVI dataset; εmax is the maximum LUE; Cs, Ts, and Ws represent the downward-regulation scalars for the respective effects of atmospheric CO2 concentration ([CO2]), temperature (Ta), and atmospheric water demand (VPD) on LUE; and min denotes the minimum value of Ts and Ws.

The effect of atmospheric CO2 concentration on GPP was calculated according to Farquhar et al. () and Collatz et al. ()where θ is the CO2 compensation point in the absence of dark respiration (ppm) and Ci is the CO2 concentration in the intercellular air spaces of the leaf (ppm), which is the product of atmospheric CO2 concentration (Ca) and the ratio of leaf internal to ambient CO2 (χ). χ is estimated (–) bywhere Kc and Ko are the Michaelis-Menten coefficient of Rubisco for carboxylation and oxygenation, respectively, expressed in partial pressure units, and Po is the partial pressure of O2 (ppm). R is the molar gas constant (8.314 J mol−1 K−1), and η* is the viscosity of water as a function of air temperature ().

Ts and Ws were calculated using the following equationswhere Tmin, Tmax, and Topt are the minimum, maximum, and optimum air temperature (°C) for photosynthetic activity, which were set to 0°, 40°, and 20.33°C, respectively (). VPD0 is an empirical coefficient of the VPD constraint equation.

Parameters εmax, θ, and VPD0 were calibrated using estimated GPP at EC towers (table S7). The nonlinear regression procedure (Proc NLIN) in the Statistical Analysis System (SAS; SAS Institute Inc., Cary, NC, USA) was applied to estimate the three parameters in the revised EC-LUE model. The revised EC-LUE model was calibrated at 50 EC towers and validated at 41 different towers (table S3). The results showed good model performance of the revised EC-LUE model for simulating biweekly GPP variations (fig. S13). To estimate global GPP, EC-LUE and MODIS models used the MERRA dataset (i.e., air temperature, VPD, PAR). Because of the different model algorithm, GIMMS3g NDVI and fPAR products were used to indicate vegetation conditions for EC-LUE and MODIS, respectively.

We performed two types of experimental simulation to evaluate the relative contribution of three main driving factors: CO2 fertilization, climate change, and satellite-based NDVI/fPAR changes. The first simulation experiment (SALL) was a normal model run, and all drivers were set to change over time to examine the responses of GPP to all environmental changes, including climate, atmospheric [CO2], and NDVI/fPAR. The second type of simulation experiments (SCLI0, SNDVI0, and SCO20) allowed two driving factors to change with time while holding the third constant at an initial baseline level. For example, the SCLI0 simulation experiment allowed NDVI and atmospheric [CO2] to change with time, while climate variables were held constant at 1982 values. SNDVI0 and SCO20 simulation experiments kept NDVI and atmospheric [CO2] constant at 1982 values and varied the other two variables.

We considered the differences between simulation results of the first type (SALL) and second type (SCO20 and SNDVI0) of experiments to estimate the sensitivity of GPP to atmospheric [CO2] (βCO2) and NDVI/fPAR (βNDVI). βCO2 and βNDVI were calculated on the basis of the following equationswhere ΔGPPi, ΔCO2i, and ΔNDVIi represent the differences of GPP simulations, atmospheric [CO2], and NDVI between two model experiments from 1982 to 2015, and ε is the residual error term.

A multiple regression approach was used to estimate GPP sensitivities to three climate variables: air temperature (βTa), VPD (βVPD), and PAR (βPAR)where ΔTai, ΔVPDi and ΔPARi represent the differences of air temperature, VPD, and photosynthetically active radiation time series between two model experiments (SALL and SCLI0), respectively. The regression coefficient β was estimated using maximum likelihood analysis.

The EC-LUE model suggested a CO2 sensitivity (βCO2) of 19.01 ± 4.01 Pg C 100 ppm−1 (Fig. 6B and table S2), which indicates a 15.7% increase of GPP with a rise of atmospheric [CO2] of 100 ppm. Our estimate is close to CO2 sensitivity derived from ecosystem models (βCO2 = 21.92 ± 4.55 Pg C 100 ppm−1; Fig. 6B) and is comparable to the observed response of NPP (net primary production) to the increased CO2 at the FACE experiment locations (13% per 100 ppm) and estimates of other ecosytem models (5 to 20% per 100 ppm) ().

A machine learning method (i.e., RF) was used to model the effects of VPD on NDVI. RF combines tree predictors such that each tree depends on the values of a random vector that is sampled independently, with the same distribution for all trees in the forest. We constructed RF models for simulating annual growing season mean NDVI at each pixel driven by air temperature, precipitation, radiation, wind speed, atmospheric [CO2] concentration, and VPD. The training data were the GIMMS3g NDVI dataset from 1982 to 2015. The R package “randomForest” used in the study was modified by A. Liaw and M. Wiener from the original Fortran by L. Breiman and A. Cutler (https://cran.r-project.org/web/packages/randomForest/).

The RF model was driven by all variables (climate and atmospheric [CO2]) changing over time (RFALL), and two factorial simulations of NDVI (RFCO20 and RFCLI0) were produced by holding one driving factor (climate or atmospheric [CO2]) constant at its initial level (first year of data) while allowing the other driving to change with time. The RFCLI0 simulation experiment allowed atmospheric [CO2] other than climate variables to vary since 1982. RFCO20 simulation experiments kept atmospheric [CO2] constant at 1982 values and varied the climate variables. At each pixel, we selected 33 years of NDVI observations out of the total 34 years (1982–2015) to develop the RF model, and the remaining 1 year of NDVI observations was used for cross-validation. The model was run 34 times to ensure that the data of each year can be selected to do model validation. The simulated NDVI of three model experiments (i.e., RFALL, RFCO20, and RFCLI0) are mean values of all 34 times simulations. The simulated NDVI only from the validation year constitutes the RFVLI dataset, which was used to examine the performance of random forest for reproducing NDVI. The simulated NDVI of RFVLI matched the GIMMS3g NDVI very well (fig. S7), and the correlation coefficient (R2) is larger than 0.90 at the 88% vegetated areas globally. The tropical forest showed the relative low R2. The relative predictive errors range from −1.2 to 1.04% globally and imply that the RF model can accurately simulate interannual variability and magnitude of NDVI.

On the basis of the three model experiments, we used the same method above shown at Eqs. 20 and 22 to estimate the sensitivity of NDVI to atmospheric [CO2] (δCO2) and five climate variables: air temperature (δTa), VPD (δVPD), PAR (δPAR), precipitation (δPrec), and wind speed (δWS)where Δ represents the differences of NDVI simulations, atmospheric [CO2], and climate variables between two model experiments from 1982 to 2015, and ε is the residual error term. The regression coefficient δ was estimated using maximum likelihood analysis. We quantified the contributions of atmospheric [CO2] and five climate variables to NDVI changes during the two periods (1982–1998 and 1999–2015) by multiplying the magnitude of their changes and sensitivity of NDVI (δ).

Terrestrial carbon cycle models

This study used a set of 10 terrestrial carbon cycle models included in the TRENDY project (version 5), which aims to further investigate the spatial trends in terrestrial ecosystem carbon cycles (): CABLE (Community Atmosphere Biosphere Land Exchange) (), CLASS (), CLM (Community Land Model) (), ISAM (Integrated Science Assessment Model) (), JSBACH (Jena Scheme for Biosphere-Atmosphere Coupling in Hamburg) (), JULES (Joint UK Land Environment Simulator) (), LPJ-GUESS (Lund-Potsdam-Jena General Ecosystem Simulator) (), LPX (Land surface Processes and eXchanges) (), ORCHIDEE (Organizing Carbon and Hydrology in Dynamic Ecosystems) (), and VISIT (Vegetation Integrated Simulator for Trace gases) (). Three TRENDY model experiments were used to evaluate the relative contribution of atmospheric CO2 concentration and climate change to GPP: (S0) no forcing change, (S1) varying CO2 only, and (S2) varying CO2 and climate. The model differences of S1 and S0 and Eq. 20 were used to estimate GPP sensitivities to atmospheric CO2 concentration, and the differences of S2 and S1 and Eq. 22 were used to estimate GPP sensitivities to three climate variables: air temperature (βTa), VPD (βVPD), and PAR (βPAR).

All models were forced with reconstructed historical climate fields and atmospheric CO2 concentrations. All models used the same forcing files, of which historical climate fields were obtained from the CRU-NCEP v4 dataset (http://dods.extra.cea.fr/data/p529viov/cruncep/), and global atmospheric CO2 concentration was obtained from a combination of ice core records and atmospheric observations ().

Analysis

A piecewise linear regression approach was used to detect changes in the trends of various variables ()where y is VPD, NDVI, LAI, GPP, and oceanic evaporation (Eocean); t is the year; α is the estimated TP of the time series, defining the timing of a trend change; β0, β1, and β2 are the regression coefficients; and ε is the residual of the fit. The investigated variable linear trend is β1 before the TP and β1 + β2 after the TP. Least squares linear regression was used to estimate α and other coefficients, where a P value of ≤0.05 was considered significant. The 5-year running means were used to quantify the TP.

In addition, the partial correlation method was used to analyze the correlation between detrended VPD and detrended NDVI/LAI of five satellite-based datasets, excluding the impacts of air temperature, radiation, and atmospheric [CO2] concentration. Because of the obvious transitions of VPD, NDVI, and LAI from 1982 to 2015, we used piecewise linear regression approach to determine the TP and then removed the trends before and after the TP before the correlation analysis.