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 CO 2 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.

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 ( 9 ). A recent study highlighted that increases in VPD rather than changes in precipitation substantially influenced vegetation productivity ( 10 ). Increasing VPD notably affects vegetation growth ( 11 – 13 ), forest mortality ( 14 ), and maize yields ( 15 ). In addition, rising VPD greatly limits land evapotranspiration in many biomes by altering the behavior of plant stomata ( 9 ). Given that the global precipitation is projected to remain steady ( 16 ), the changing VPD and soil drying would likely constrain plant carbon uptake and water use in terrestrial ecosystems ( 17 ). 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).

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 ( 1 ). 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 ( 2 ) but also in humid regions ( 3 ). Although the long-term trend of globally averaged land surface relative humidity remains insignificant ( 4 , 5 ), a sharp decrease has been observed since 2000 ( 6 , 7 ), implying a sharp increase in land surface VPD. However, the causes of changing atmospheric water demand are still unclear ( 8 ).

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 (4, 5, 7), 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 (18). 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 (19). Most of the moisture over land is transported from the oceans, which accounts for 35% of precipitation and 55% of evapotranspiration over land (19). We analyzed long-term changes of oceanic evaporation based on a global oceanic evaporation dataset [Objectively Analyzed Air–Sea Fluxes (OAFlux)] (20). The almost 60-year time series showed that the decadal change of global oceanic evaporation (E ocean ) was marked by a distinct transition from an upward to a downward trend around 1998 (Fig. 2A). The global oceanic E ocean 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 E ocean after 1999 (Fig. 2B). Rhein et al. (16) reported stalled increases of sea surface temperature after the late 1990s based on multiple global datasets, which substantially limited oceanic evaporation (20). 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 (8, 21). 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. (22) 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 (23). Previous study reported a decreased trend in the global land evapotranspiration after the late 1990s limited by soil moisture supply (24).

Fig. 2 Comparison of oceanic evaporation (E ocean ) 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] (25). 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 CO 2 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 [CO 2 ] 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 CO 2 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 (26, 27). 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 CO 2 concentration. (C) Contributions of climate variables, NDVI/fPAR, and atmospheric CO 2 concentration to GPP changes over the two periods. Three climate variables are included: vapor pressure deficit (VPD), air temperature (T a ), 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 CO 2 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 [CO 2 ] (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, CO 2 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 CO 2 fertilization effect (0.38 ± 0.08 Pg C year−1) after the late 1990s with the rising rate of atmospheric CO 2 concentration by 2.02 ± 0.01 ppm year−1. From 1982 to the late 1990s, CO 2 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 [CO 2 ] (Fig. 6C and table S2).