SCIENCE ADVANCES | RESEARCH ARTICLE

ECOLOGY Copyright © 2019

The Authors, some

Increased atmospheric vapor pressure deficit reduces rights reserved;

exclusive licensee
global vegetation growth American Association
for the Advancement
of Science. No claim to
Wenping Yuan1,2*, Yi Zheng1, Shilong Piao3, Philippe Ciais4, Danica Lombardozzi5, original U.S. Government
Yingping Wang6,7, Youngryel Ryu8, Guixing Chen1,2, Wenjie Dong1,2, Zhongming Hu9, Atul K. Jain10, Works. Distributed
Chongya Jiang11, Etsushi Kato12, Shihua Li1, Sebastian Lienert13, Shuguang Liu14, under a Creative
Julia E.M.S. Nabel15, Zhangcai Qin1,2, Timothy Quine16, Stephen Sitch16, William K. Smith17, Commons Attribution
Fan Wang1,2, Chaoyang Wu18, Zhiqiang Xiao19, Song Yang1,2 NonCommercial
License 4.0 (CC BY-NC).

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

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

VPD, which offset the positive CO2 fertilization effect. Six Earth system models have consistently projected con-
tinuous increases of VPD throughout the current century. Our results highlight that the impacts of VPD on veg-
etation growth should be adequately considered to assess ecosystem responses to future climate conditions.

INTRODUCTION surface relative humidity remains insignificant (4, 5), a sharp decrease
Vapor pressure deficit (VPD), which describes the difference between has been observed since 2000 (6, 7), implying a sharp increase in land
the water vapor pressure at saturation and the actual water vapor pres- surface VPD. However, the causes of changing atmospheric water de-
sure for a given temperature, is an important driver of atmospheric mand are still unclear (8).
water demand for plants (1). Rising air temperature increases saturated Changes of VPD are important for terrestrial ecosystem structure
water vapor pressure at a rate of approximately 7%/K according to the and function. Leaf and canopy photosynthetic rates decline when at-
Clauius-Clapeyron relationship, which will drive an increase in VPD if mospheric VPD increases due to stomatal closure (9). A recent study
the actual atmospheric water vapor content does not increase by exactly highlighted that increases in VPD rather than changes in precipitation
the same amount as saturated vapor pressure (SVP). Numerous studies substantially influenced vegetation productivity (10). Increasing VPD
have indicated substantial changes of relative humidity (ratio of actual notably affects vegetation growth (11–13), forest mortality (14),
water vapor pressure to saturated water vapor pressure) not only in and maize yields (15). In addition, rising VPD greatly limits land
continental areas located far from oceanic humidity (2) but also in hu- evapotranspiration in many biomes by altering the behavior of plant
mid regions (3). Although the long-term trend of globally averaged land stomata (9). Given that the global precipitation is projected to remain
steady (16), the changing VPD and soil drying would likely constrain
1
School of Atmospheric Sciences, Guangdong Province Key Laboratory for Climate
plant carbon uptake and water use in terrestrial ecosystems (17). How-
Change and Natural Disaster Studies, Zhuhai Key Laboratory of Dynamics Urban ever, the large-scale constraints of VPD changes on vegetation growth
Climate and Ecology, Sun Yat-sen University, Zhuhai, Guangdong 510245, China.
2
have not yet been quantified. In this study, we determined the changes
Southern Marine Science and Engineering Guangdong Laboratory, Zhuhai in VPD trends through observation-based global climate datasets, and
519082, China. 3Sino-French Institute for Earth System Science, College of Urban
and Environmental Sciences, Peking University, Beijing 100871, China. 4Laboratoire then quantified the impacts of these VPD changes on vegetation
des Sciences du Climat et de l’Environnement, CEA CNRS UVSQ, Gif-sur-Yvette growth and productivity, using satellite-based vegetation index [i.e.,
91191, France. 5Terrestrial Sciences Section, Climate and Global Dynamics, National normalized difference vegetation index (NDVI)] and leaf area index
Center for Atmospheric Research, Boulder, CO 80305, USA. 6CSIRO, Oceans and At-
mosphere, Private Bag 1, Aspendale, Victoria 3195, Australia. 7South China Botanical (LAI), tree-ring width chronologies, and remotely sensed estimates
Garden, Chinese Academy of Sciences, Guangzhou 510650, China. 8Department of of gross primary production (GPP).
Landscape Architecture and Rural Systems Engineering, Seoul National University,
Seoul, Republic of Korea. 9School of Geography, South China Normal University,
Guangzhou 510631, China. 10Department of Atmospheric Sciences, University of Illinois
at Urbana-Champaign, Urbana, IL 61801, USA. 11College of Agricultural, Consumer & RESULTS
Environmental Sciences, University of Illinois at Urbana-Champaign, Urbana, IL 61801, This study used four observation-based globally gridded climate
USA. 12Global Environment Program, Research & Development Division, The Institute datasets—CRU (Climatic Research Unit), ERA-Interim, HadISDH, and
of Applied Energy (IAE), Shimbashi SY Bldg., 1-14-2 Nishi-Shimbashi Minato, Tokyo
105-0003, Japan. 13Climate and Environmental Physics, Physics Institute and Oeschger MERRA (Modern-Era Retrospective analysis for Research and Applica-
Centre for Climate Change Research, University of Bern, Bern, Switzerland. 14National tions) (table S1)—to analyze the long-term trend of VPD over vegetated
Engineering Laboratory for Applied Technology of Forestry & Ecology in South China and land. Similar to previous analyses (4, 5, 7), anomalies in all four datasets
College of Biological Science and Technology, Central South University of Forestry and
Technology, Changsha 410004, China. 15Max Planck Institute for Meteorology, 20146
showed that VPD trends were temporally and spatially heteroge-
Hamburg, Germany. 16Department of Geography, College of Life and Environmental neous over recent decades (Fig. 1). A piecewise linear regression
Sciences, University of Exeter, EX4 4RJ Exeter, UK. 17School of Natural Resources and method was used to quantify the change in trends and detect the
the Environment, University of Arizona, Tucson, AZ 85721, USA. 18Institute of Geographic potential turning point (TP) in each dataset. It was observed that
Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101,
China. 19Faculty of Geography, Beijing Normal University, Beijing 100875, China. VPD increased slightly before the late 1990s but increased more
*Corresponding author. Email: yuanwp3@mail.sysu.edu.cn strongly afterward with 1.66 to 17 times larger trends according to

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 1 of 12

SCIENCE ADVANCES | RESEARCH ARTICLE

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

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.

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.

the four datasets (fig. S1). The datasets showed that 53 to 64% of vege- growing season mean VPD between two periods of 1982–1986 and
tated areas experienced increased VPD trends since the late 1990s (fig. 2011–2015 (fig. S3A). On average, the annual growing season mean
S2). To illustrate the magnitude and spatial variability of VPD change, VPD of 2011–2015 was 11.26% higher than that of 1982–1986, and
we calculated the global pattern of the percentage change of annual the VPD increased larger than 5% in more than 53% area. In addition,

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 2 of 12

SCIENCE ADVANCES | RESEARCH ARTICLE

the increases of global mean VPD over 12 months were positively a peak of 1197 mm year−1 in 1998 to a low 1166 mm year−1 in 2015
correlated with the mean VPD values of 1982–1986 at more than (Fig. 2A), and 76% of the sea surface revealed a decreased Eocean after
64.5% areas (fig. S3B), which implies that the higher VPD increases 1999 (Fig. 2B). Rhein et al. (16) reported stalled increases of sea sur-
in the months with high VPD. face temperature after the late 1990s based on multiple global data-
Apart from HadISDH, datasets showed that the increased saturated sets, which substantially limited oceanic evaporation (20). Some
water vapor pressure and decreased actual water vapor pressure jointly studies using global climate models (GCMs) also highlighted that
determined the increases of VPD after the TP. On average, the rate of VPD trends over land were predominantly explained by dynamic
increase in saturated water vapor was 1.43 to 1.64 times higher after the mechanisms related to moisture supply from oceanic source regions
TP year than before, and the actual water vapor exhibited stalled or de- (8, 21). Changes in the recycling of atmospheric moisture over land
creased trends (fig. S4). Increased air temperature explains the changes controlled by soil moisture in supply-limited regions may be an ad-
in saturated water vapor pressure (fig. S4). The HadISDH dataset indi- ditional contribution to the observed increase of VPD. Koster et al.
cates a decrease in saturated water vapor because of large spatial gaps in (22) showed that moisture variability contributed to total precipita-
the dataset. tion variance in mid-northern latitude regions such as the western United
A change of oceanic evaporation is the most important mechanism States. Drier soils evaporate less and thus lead to lower water vapor in
for the observed decrease in actual water vapor pressure over the land the atmosphere (23). Previous study reported a decreased trend in the
(18). Oceanic evaporation is the most important source of atmosphere global land evapotranspiration after the late 1990s limited by soil mois-
water vapor, and approximately 85% of atmospheric water vapor is ture supply (24).

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

evaporated from oceans, with the remaining 15% coming from evapo- Figure 3 illustrates that the satellite-based NDVI substantially
ration and transpiration over land (19). Most of the moisture over land increased from 1982 to 1998 (y = 0.0014x − 1.86, R 2 = 0.43,
is transported from the oceans, which accounts for 35% of precipitation P < 0.05), while NDVI remained constant and then stalled after 1999
and 55% of evapotranspiration over land (19). We analyzed long-term (y = −0.0004x + 1.23, R2 = 0.06, P = 0.65) (Fig. 3A). From 1982 to 1998,
changes of oceanic evaporation based on a global oceanic evaporation approximately 84% of the vegetation surface showed an increased
dataset [Objectively Analyzed Air–Sea Fluxes (OAFlux)] (20). The NDVI trend (28.50% with a significant increase; Fig. 4A). In compari-
almost 60-year time series showed that the decadal change of global son, after 1999, the trends of NDVI over many regions reversed, and
oceanic evaporation (Eocean) was marked by a distinct transition from 59% of vegetation areas showed a pronounced NDVI browning
an upward to a downward trend around 1998 (Fig. 2A). The global (decreasing) trend (21.50% with a significant decrease; Figs. 3B and
oceanic Eocean has decreased by approximately 2.08 mm year−1, from 4). Mean NDVI trends for 12 months after 1999 were lower than those

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.

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 3 of 12

SCIENCE ADVANCES | RESEARCH ARTICLE

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

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.

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 4 of 12

SCIENCE ADVANCES | RESEARCH ARTICLE

from 1982 to 1998 over globally vegetated areas (Fig. 3C). Moreover, we with a significant negative correlation) (Fig. 5A). Similarly, four de-
analyzed long-term trends of LAI based on four global LAI datasets trended satellite-based LAI correlated negatively with detrended
[Global Land Surface Satellite (GLASS), GLOMap, LAI3g, and Terres- VPD over 65 to 70% of vegetated areas (16 to 22% with a significant
trial Climate Data Record (TCDR); table S1] (25). Despite the large var- negative correlation) (Fig. 5, B to E). In addition, all five satellite-
iability of the estimated interannual LAI among the four products, all based datasets show highly consistent signs of correlation with
four LAI datasets exhibited a transition from increasing trends before VPD, and at least three datasets revealed consistently negative cor-
the late 1990s to decreasing trends afterward (fig. S5). The LAI showed a relations with VPD over 72% of vegetated area (Fig. 5F). A machine
decreasing trend since the late 1990s over vegetated areas of 64.72, learning method [i.e., random forest (RF)] was used to reconstruct
72.62, 62.73, and 80.11% for GLASS, GLOMap, LAI3g, and TCDR da- NDVI based on atmospheric [CO2] concentration and five climate
tasets, respectively (fig. S6). The differences of NDVI and LAI trends factors (air temperature, precipitation, radiation, wind speed, and
during these two periods are the opposite of VPD trends derived from VPD) over the last 34 years in each pixel (fig. S7) and then model
four VPD datasets. experiments were applied to separate the impacts of VPD as well
Partial correlation analysis indicated significant correlations of as of other variables (see Materials and Methods). Globally, the model
detrended VPD with detrended NDVI and LAI when the impacts experiments suggest that the atmospheric CO2 concentration, air tem-
of air temperature, radiation, and atmospheric CO2 concentration perature, and VPD are the most important contributors for the varia-
were excluded (Fig. 5). Detrended NDVI over 62% of the vegetated bility of NDVI (fig. S8A). Rising VPD was found to significantly
areas shows a negative correlation with detrended VPD (about 14% decrease NDVI, indicated by the larger negative NDVI differences from

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

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.

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 5 of 12

SCIENCE ADVANCES | RESEARCH ARTICLE

1999 to 2015, suggesting that substantial increases of VPD strongly bCO2 = 19.01 ± 4.01 Pg C 100 ppm−1). On the basis of the estimated
limited NDVI (fig. S8B). GPP sensitivity, we estimated the contributions of climate variables,
This study used two satellite-based models [revised eddy CO2 fertilization, and vegetation index to global GPP over the two
covariance–light use efficiency (EC-LUE) and Moderate Resolution study periods (table S2). After the late 1990s, VPD increased by
Imaging Spectroradiometer (MODIS)] to investigate the impacts of 0.0017 ± 0.0001 kPa year−1 according to the CRU dataset (fig. S1),
VPD on long-term changes of global GPP (26, 27). EC-LUE and which resulted in GPP decreases of 0.23 ± 0.09 Pg C year−1 and
MODIS showed quite similar long-term trends of GPP, with a sig- 0.31 ± 0.11 Pg C year−1 according to the EC-LUE and MODIS models,
nificantly increased trend from 1982 to the late 1990s, averaged at respectively (Fig. 6C and table S2). The VPD-induced GPP decreases
0.73 Pg C year−1 (P < 0.05; from 1982 to 1998) and 0.26 Pg C year−1 partly counteract the CO2 fertilization effect (0.38 ± 0.08 Pg C year−1)
(P < 0.05; from 1982 to 1997) over globally vegetated area, respectively after the late 1990s with the rising rate of atmospheric CO2 concentra-
(Fig. 6A). The GPP trends then stalled and decreased afterward tion by 2.02 ± 0.01 ppm year−1. From 1982 to the late 1990s, CO2 fer-
(−0.016 Pg C year−1, P = 0.67 and −0.032 Pg C year−1, P = 0.44) tilization played a dominant role in the GPP increase (Fig. 6C).
(Fig. 6A). The GPP trends derived from the two models during the According to the EC-LUE model, GPP increases of 0.28 ± 0.15 Pg C
two periods are the opposite of VPD trends derived from the four year−1 occurred because of the rising atmospheric [CO2] (Fig. 6C and
VPD datasets. table S2).
To quantify the impacts of VPD on GPP, we further explored GPP We further investigated the impacts of VPD on LUE using measure-
sensitivity to climate variables (i.e., air temperature, VPD, and radia- ments of global EC towers (see Materials and Methods; table S3). Be-

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

tion), atmospheric CO2 concentration, and satellite-based NDVI/fPAR cause VPD correlated strongly with air temperature, we excluded the
(see Materials and Methods; Fig. 6B). Two satellite-based models effects of air temperature to investigate the impacts of VPD by binning
showed the similar GPP sensitivity to VPD, whereby global GPP de- the observations for ranges of air temperature. We binned the LUE data
creased by 13.82 ± 3.12 Pg C and 18.29 ± 3.65 Pg C with a VPD in- for different restricted ranges of air temperature and found strong neg-
crease of 0.1 kPa (Fig. 6B), which is comparable to the GPP increase ative correlations between VPD and LUE for almost all air temperature
with a 100–parts per million (ppm) rise of atmospheric [CO2] (i.e., ranges (fig. S9 and table S3). In addition, we investigated the impacts of

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).

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 6 of 12

SCIENCE ADVANCES | RESEARCH ARTICLE

VPD on vegetation growth using a global comprehensive dataset of on climate observations from global meteorological stations (31, 32).
tree-ring width measurements from 171 locations with temporal cover- ERA-Interim and MERRA datasets were reanalysis products based
age from 1982 until at least 2005. Partial correlation analysis showed on Integrated Forecast System of European Centre for Medium-Range
that the detrended VPD derived from the four datasets correlated with Weather Forecasts (ECMWF-IFS) (33) and the Goddard Earth
detrended tree-ring width at most sites (56 to 72%) when the impacts of Observing System Data Assimilation System Version 5 (GEOS-5)
air temperature, radiation, and atmospheric CO2 concentration were (34), respectively. VPD was calculated on the basis of different variables
excluded (fig. S10, A to D). We compared the differences of mean of four datasets (35)
tree-ring width values between before and after 1998 and observed CRU:
smaller tree-ring widths after 1998 compared to those before 1998 at
64% of sites (25% sites with significant level) (fig. S10E). VPD ¼ SVP  AVP ð1Þ

ERA-Interim:
DISCUSSION
Our results support increased VPD being part of the drivers of the 17:67Td

widespread drought-related forest mortality over the past decades, AVP ¼ 6:112  fw  eTd þ243:5 ð2Þ
which has been observed in multiple biomes and on all vegetated con-
tinents (28, 29). Increased VPD may trigger stomatal closure to avoid VPD ¼ SVP  AVP ð3Þ

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

excess water loss due to the high evaporative demand of the air (12),
leading to a negative carbon balance that depletes carbohydrate reserves HadISDH and MERRA:
and results in tissue-level carbohydrate starvation (28). In addition, re-
duced soil water supply coupled with high evaporative demand causes RH
AVP ¼  SVP ð4Þ
xylem conduits and the rhizosphere to cavitate (become air-filled), 100
stopping the flow of water, desiccating plant tissues, and leading to plant
death (28). Previous studies reported that increased VPD explained 82% VPD ¼ SVP  AVP ð5Þ
of the warm season drought stress in the southwestern United States,
which correlated to changes of forest productivity and mortality (14). where SVP and AVP are saturated vapor pressure and actual vapor
In addition, enhanced VPD limits tree growth even before soil moisture pressure (kPa), respectively. Td is the dew point temperature (°C).
begins to be limiting (17, 30). RH is the land relative humidity (%).
We examined whether terrestrial ecosystem models can adequately
capture the observed responses of vegetation growth to increased VPD 17:67Ta

after the late 1990s from 10 terrestrial ecosystem models. We found that SVP ¼ 6:112  fw  eTa þ243:5 ð6Þ
the simulated GPP trends of most models did not match the GPP trends
documented above (fig. S11). Only the CLASS (Canadian Land Surface
fw ¼ 1 þ 7  104 þ 3:46  106 Pmst ð7Þ
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
5:625
ðTa þ 273:16Þ
VPD than two satellite-based models (i.e., EC-LUE and MODIS) Pmst ¼ Pmsl ð8Þ
(Fig. 6B and table S2). ðTa þ 273:16Þ þ 0:0065  Z
Our results imply that most terrestrial ecosystem models cannot
capture vegetation responses to VPD. Thus, problems reproducing where Ta is the land air temperature (°C). Z is the altitude (m). Pmst is
the observed long-term vegetation responses to climate variability the air pressure (hPa), and Pmsl is the air pressure at mean sea level
may challenge their ability to predict the future evolution of the carbon (1013.25 hPa). In addition, the OAFlux dataset was used to examine
cycle. Earth system models (ESMs) participating in the CMIP5 the variability of oceanic evaporation (table S1) (20).
(Coupled Model Intercomparison Project Phase 5) (table S5) project We used the newest release of the advanced very high resolution
a continuous increase of VPD until the end of this century (Fig. 1). radiometer (AVHRR) NDVI to indicate vegetation growth from 1982
The globally averaged VPD is 0.12 kPa higher in 2090–2100 than in to 2015. The AVHRR is a nonstationary NDVI version 3 dataset made
1980–1999 (Fig. 1). The ESMs used in Fig. 1 showed good performance available by NASA’s Global Inventory Modeling and Monitoring Study
when reproducing historical variations of VPD (table S6), providing third-generation dataset (GIMMS3g) group (36). GIMMS3g contains
confidence in the projected increases of VPD during future decades. global NDVI observations at approximately 8-km spatial resolution
The results of our analysis suggest that this projected increased VPD and bimonthly temporal resolution, derived from AVHRR channels
might have a substantially negative impact on vegetation, which must 1 and 2, corresponding to red (0.58 to 0.68 mm) and infrared (0.73 to
be examined carefully when evaluating future carbon cycle responses. 1.1 mm) wavelengths, respectively. Each 15-day data value is the result of
maximum value compositing, a process that aims to minimize the in-
fluence of atmospheric contamination from aerosols and clouds. More-
MATERIALS AND METHODS over, this study analyzed long-term trends of LAI based on four global
Datasets satellite LAI products (table S1): GLASS (version 4) (37), GLOBMap
Four global climate datasets were used to investigate the long-term (38), LAI3g (39), and the TCDR (40).
changes of atmospheric VPD, including CRU, ERA-Interim, HadISDH, We calculated the annual growing season mean NDVI and LAI by
and MERRA. Monthly gridded CRU and HadISDH datasets were based averaging monthly NDVI and LAI values with monthly mean temper-

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 7 of 12

SCIENCE ADVANCES | RESEARCH ARTICLE

atures above 0°C. We also calculated multiyear averaged monthly mean scalars for the respective effects of atmospheric CO2 concentration
temperatures from the CRU dataset to ensure that the same growing ([CO2]), temperature (Ta), and atmospheric water demand (VPD) on
season land mask was used over the entire period (1982–2015). The LUE; and min denotes the minimum value of Ts and Ws.
global mean NDVI and LAI values were calculated by the average of The effect of atmospheric CO 2 concentration on GPP was
the annual growing season mean NDVI and LAI, excluding unvege- calculated according to Farquhar et al. (44) and Collatz et al. (45)
tated regions. The MODIS land cover type product (MCD12Q1) was
used to identify the vegetated regions. Ci  q
We calculated the LUE (g C m−2 MJ−1) based on EC measure- Cs ¼ ð11Þ
ments from the FLUXNET2015 dataset (www.fluxdata.org) to ex- Ci þ 2q
amine the correlation between LUE and VPD (table S4)
Ci ¼ Ca  c ð12Þ
GPP
LUE ¼ ð9Þ
fPAR  PAR where q is the CO2 compensation point in the absence of dark respira-
tion (ppm) and Ci is the CO2 concentration in the intercellular air
where GPP indicates the estimated GPP values from EC measurements, spaces of the leaf (ppm), which is the product of atmospheric CO2 con-
PAR is photosynthetically active radiation (MJ m−2), and fPAR is the centration (Ca) and the ratio of leaf internal to ambient CO2 (c). c is
fraction of PAR absorbed by the vegetation canopy calculated by estimated (46–48) by
GIMMS3g NDVI (41).

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

The tree-ring width measurements around the world were used g
c¼ pffiffiffiffiffiffiffiffiffiffi ð13Þ
from the International Tree-Ring Data Bank (ITRDB) (42). The wood gþ VPD
samples were taken and processed following standard protocols and
taking two radial cores per tree at 1.3 m. Tree-ring width measure- rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
356:51K
ments were detrended and standardized by the scientists who g¼ ð14Þ
contributed the chronologies to the ITRDB. Each local chronology rep- 1:6h
resents the average growth of several trees (typically more than 10) of
the same species growing at the same site. The temporal span of the

P0
tree-ring data series selected began at 1982, lasting at least until 2005. K ¼ Kc 1 þ ð15Þ
Eventually, 171 sites were analyzed and each chronology of the sites is a K0
representation of annual tree-ring width.
This study conducts the partial correlation analysis between VPD 79:43ðTa 298:15Þ

and tree-ring width by excluding the impacts of air temperature, ra- Kc ¼ 39:97  e 298:15RTa ð16Þ
diation, and atmospheric CO2 concentration. Air temperature and
PAR from MERRA dataset were used. For atmospheric CO2 concen- 36:38ðTa 298:15Þ

tration, this study used the GLOBALVIEW-CO2 product, which pro- Ko ¼ 27480  e 298:15RTa ð17Þ
vides observations of atmospheric CO2 concentration at 7-day
intervals over 313 global air-sampling sites (43). If missing 7-day where Kc and Ko are the Michaelis-Menten coefficient of Rubisco for
data accounted for >20% of all data for an entire year, then the carboxylation and oxygenation, respectively, expressed in partial pres-
value for that year was indicated as “missing.” For a site to be in- sure units, and Po is the partial pressure of O2 (ppm). R is the molar gas
cluded in this study, it had to have at least 10 years of observations. constant (8.314 J mol−1 K−1), and h* is the viscosity of water as a func-
Eventually, 77 sites were included equally in the calculation of tion of air temperature (49).
global monthly mean CO2 concentration without any weighting Ts and Ws were calculated using the following equations
of individual sites.
ðTa  Tmin Þ  ðTa  Tmax Þ
Satellite-based GPP model Ts ¼
ðTa  Tmin Þ  ðTa  Tmax Þ  ðTa  Topt Þ  ðTa  Topt Þ
We used two satellite-based GPP models to investigate the impacts of
VPD on vegetation GPP. The first model is the MODIS-GPP model, ð18Þ
and this study used long-term global MODIS GPP dataset driven by
GIMMS fPAR data (27). VPD0
The second model is the revised EC-LUE model (26), derived by Ws ¼ ð19Þ
(i) integrating the impact of atmospheric CO2 concentration on GPP VPD þ VPD0
and (ii) adding the limit of VPD to GPP. The revised EC-LUE model
simulates terrestrial ecosystem GPP as where Tmin, Tmax, and Topt are the minimum, maximum, and optimum
air temperature (°C) for photosynthetic activity, which were set to 0°,
GPP ¼ PAR  fPAR  emax  Cs  minðTs ; Ws Þ ð10Þ 40°, and 20.33°C, respectively (50). VPD0 is an empirical coefficient
of the VPD constraint equation.
where PAR is the incident photosynthetically active radiation (MJ m−2) Parameters emax, q, and VPD0 were calibrated using estimated
per time period (e.g., day); fPAR is the fraction of PAR absorbed by the GPP at EC towers (table S7). The nonlinear regression procedure
vegetation canopy calculated by the GIMMS3g NDVI dataset; emax is (Proc NLIN) in the Statistical Analysis System (SAS; SAS Institute
the maximum LUE; Cs, Ts, and Ws represent the downward-regulation Inc., Cary, NC, USA) was applied to estimate the three parameters

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 8 of 12

SCIENCE ADVANCES | RESEARCH ARTICLE

in the revised EC-LUE model. The revised EC-LUE model was cali- pends on the values of a random vector that is sampled independently,
brated at 50 EC towers and validated at 41 different towers (table S3). with the same distribution for all trees in the forest. We constructed RF
The results showed good model performance of the revised EC-LUE models for simulating annual growing season mean NDVI at each pixel
model for simulating biweekly GPP variations (fig. S13). To estimate driven by air temperature, precipitation, radiation, wind speed, atmo-
global GPP, EC-LUE and MODIS models used the MERRA dataset spheric [CO2] concentration, and VPD. The training data were the
(i.e., air temperature, VPD, PAR). Because of the different model GIMMS3g NDVI dataset from 1982 to 2015. The R package “random-
algorithm, GIMMS3g NDVI and fPAR products were used to indicate Forest” used in the study was modified by A. Liaw and M. Wiener from
vegetation conditions for EC-LUE and MODIS, respectively. the original Fortran by L. Breiman and A. Cutler (https://cran.r-project.
We performed two types of experimental simulation to evaluate org/web/packages/randomForest/).
the relative contribution of three main driving factors: CO2 fertilization, The RF model was driven by all variables (climate and atmospheric
climate change, and satellite-based NDVI/fPAR changes. The first [CO2]) changing over time (RFALL), and two factorial simulations of
simulation experiment (SALL) was a normal model run, and all drivers NDVI (RFCO20 and RFCLI0) were produced by holding one driving
were set to change over time to examine the responses of GPP to all factor (climate or atmospheric [CO2]) constant at its initial level (first
environmental changes, including climate, atmospheric [CO2], and year of data) while allowing the other driving to change with time. The
NDVI/fPAR. The second type of simulation experiments (SCLI0, SNDVI0, RFCLI0 simulation experiment allowed atmospheric [CO2] other than
and SCO20) allowed two driving factors to change with time while climate variables to vary since 1982. RFCO20 simulation experiments
holding the third constant at an initial baseline level. For example, the kept atmospheric [CO2] constant at 1982 values and varied the climate

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

SCLI0 simulation experiment allowed NDVI and atmospheric [CO2] to variables. At each pixel, we selected 33 years of NDVI observations out
change with time, while climate variables were held constant at 1982 of the total 34 years (1982–2015) to develop the RF model, and the re-
values. SNDVI0 and SCO20 simulation experiments kept NDVI and atmo- maining 1 year of NDVI observations was used for cross-validation.
spheric [CO2] constant at 1982 values and varied the other two variables. The model was run 34 times to ensure that the data of each year can
We considered the differences between simulation results of the first be selected to do model validation. The simulated NDVI of three model
type (SALL) and second type (SCO20 and SNDVI0) of experiments to experiments (i.e., RFALL, RFCO20, and RFCLI0) are mean values of all
estimate the sensitivity of GPP to atmospheric [CO2] (bCO2) and 34 times simulations. The simulated NDVI only from the validation
NDVI/fPAR (bNDVI). bCO2 and bNDVI were calculated on the basis of year constitutes the RFVLI dataset, which was used to examine the
the following equations 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
DGPPðSALL SCO20 Þi ¼ bCO2  DCO2ðSALL SCO20 Þi þ e ð20Þ
areas globally. The tropical forest showed the relative low R2. The rela-
tive predictive errors range from −1.2 to 1.04% globally and imply that
the RF model can accurately simulate interannual variability and mag-
DGPPðSALL SNDVI0 Þi ¼ bNDVI  DNDVIðSALL SNDVI0 Þi þ e ð21Þ
nitude of NDVI.
On the basis of the three model experiments, we used the same
where DGPPi, DCO2i, and DNDVIi represent the differences of GPP method above shown at Eqs. 20 and 22 to estimate the sensitivity of
simulations, atmospheric [CO2], and NDVI between two model NDVI to atmospheric [CO2] (dCO2) and five climate variables: air tem-
experiments from 1982 to 2015, and e is the residual error term. perature (dTa), VPD (dVPD), PAR (dPAR), precipitation (dPrec), and wind
A multiple regression approach was used to estimate GPP sensitiv- speed (dWS)
ities to three climate variables: air temperature (bTa), VPD (bVPD), and
PAR (bPAR) DNDVIðRFALL RFCO20 Þi ¼ dCO2  DCO2ðRFALL RFCO20 Þi þ e ð23Þ

DGPPðSALL SCLI0 Þi ¼ bTa  DTaðSALL SCLI0 Þi þ bVPD DNDVIðRFALL RFCLI0 Þi ¼ dTa  DTaðRFALL RFCLI0 Þi þ dVPD
 DVPDðSALL SCLI0 Þi þ bPAR  DVPDðRFALL RFCLI0 Þi þ dPAR
 DPARðSALL SCLI0 Þi þ e ð22Þ  DPARðRFALL RFCLI0 Þi þ dPrec
 DPrecðRFALL RFCLI0 Þi þ dWS
where DTai, DVPDi and DPARi represent the differences of air tempera-
ture, VPD, and photosynthetically active radiation time series between  DWSðRFALL RFCLI0 Þi þ e ð24Þ
two model experiments (SALL and SCLI0), respectively. The regression
coefficient b was estimated using maximum likelihood analysis. where D represents the differences of NDVI simulations, atmospheric
The EC-LUE model suggested a CO 2 sensitivity (bCO2 ) of [CO2], and climate variables between two model experiments from
19.01 ± 4.01 Pg C 100 ppm−1 (Fig. 6B and table S2), which indicates 1982 to 2015, and e is the residual error term. The regression coefficient
a 15.7% increase of GPP with a rise of atmospheric [CO2] of 100 ppm. d was estimated using maximum likelihood analysis. We quantified the
Our estimate is close to CO2 sensitivity derived from ecosystem contributions of atmospheric [CO2] and five climate variables to NDVI
models (bCO2 = 21.92 ± 4.55 Pg C 100 ppm−1; Fig. 6B) and is compa- changes during the two periods (1982–1998 and 1999–2015) by multi-
rable to the observed response of NPP (net primary production) to the plying the magnitude of their changes and sensitivity of NDVI (d).
increased CO2 at the FACE experiment locations (13% per 100 ppm)
and estimates of other ecosytem models (5 to 20% per 100 ppm) (51). Terrestrial carbon cycle models
A machine learning method (i.e., RF) was used to model the effects This study used a set of 10 terrestrial carbon cycle models included in
of VPD on NDVI. RF combines tree predictors such that each tree de- the TRENDY project (version 5), which aims to further investigate the

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 9 of 12

SCIENCE ADVANCES | RESEARCH ARTICLE

spatial trends in terrestrial ecosystem carbon cycles (52): CABLE (Com- Fig. S7. Model validation of random forest models for simulating NDVI.
Fig. S8. Environmental regulations on long-term changes of global NDVI.
munity Atmosphere Biosphere Land Exchange) (53), CLASS (54), CLM
Fig. S9. Correlations of LUE and VPD at the different temperature ranges taking DE-Tha site as
(Community Land Model) (55), ISAM (Integrated Science Assessment an example.
Model) (56), JSBACH (Jena Scheme for Biosphere-Atmosphere Cou- Fig. S10. Correlations between VPD and tree-ring width.
pling in Hamburg) (57), JULES (Joint UK Land Environment Simula- Fig. S11. Comparison on changes of global mean GPP trend simulated by ecosystem models.
tor) (58), LPJ-GUESS (Lund-Potsdam-Jena General Ecosystem Fig. S12. Projected future changes in VPD.
Fig. S13. Validation of EC-LUE model.
Simulator) (59), LPX (Land surface Processes and eXchanges) (60), Table S1. Climate and satellite datasets used in this study.
ORCHIDEE (Organizing Carbon and Hydrology in Dynamic Ecosys- Table S2. Responses of GPP simulated by EC-LUE, MODIS, and TRENDY models to climate
tems) (61), and VISIT (Vegetation Integrated Simulator for Trace gases) variables, satellite-based NDVI and fPAR, and atmospheric CO2 concentration.
(62). Three TRENDY model experiments were used to evaluate the rel- Table S3. Name, location, and durations of the study EC sites used for revised EC-LUE model
calibration and validation.
ative contribution of atmospheric CO2 concentration and climate
Table S4. Correlations between VPD and LUE at different temperature ranges.
change to GPP: (S0) no forcing change, (S1) varying CO2 only, and Table S5. CMIP5 models used to estimate VPD from 1850 to 2100.
(S2) varying CO2 and climate. The model differences of S1 and S0 Table S6. Correlation matrixes for global VPD simulated by the six CMIP5 ESMs and four
and Eq. 20 were used to estimate GPP sensitivities to atmospheric historical datasets (CRU, ERA-Interim, HadISDH, and MERRA).
CO2 concentration, and the differences of S2 and S1 and Eq. 22 were Table S7. Model parameters of EC-LUE for different vegetation types.

used to estimate GPP sensitivities to three climate variables: air tem-

perature (bTa), VPD (bVPD), and PAR (bPAR).

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

All models were forced with reconstructed historical climate fields REFERENCES AND NOTES
and atmospheric CO2 concentrations. All models used the same forcing 1. H. M. Rawson, J. E. Begg, R. G. Woodward, The effect of atmospheric humidity on
files, of which historical climate fields were obtained from the CRU- photosynthesis, transpiration and water use efficiency of leaves of several plant species.
Planta 134, 5–10 (1977).
NCEP v4 dataset (http://dods.extra.cea.fr/data/p529viov/cruncep/), 2. D. W. Pierce, A. L. Westerling, J. Oyler, Future humidity trends over the western United
and global atmospheric CO2 concentration was obtained from a com- States in the CMIP5 global climate models and variable infiltration capacity hydrological
bination of ice core records and atmospheric observations (63). modeling system. Hydrol. Earth Syst. Sci. 17, 1833–1850 (2013).
3. W. A. van Wijngaarden, L. A. Vincent, Trends in relative humidity in Canada from
1953–2003. Bull. Am. Meteorol. Soc., 4633–4636 (2004).
Analysis
4. A. G. Dai, Recent climatology, variability, and trends in global surface humidity. J. Clim. 19,
A piecewise linear regression approach was used to detect changes in 3589–3606 (2006).
the trends of various variables (64) 5. K. M. Willett, P. D. Jones, N. P. Gillett, P. W. Thorne, Recent changes in surface humidity:
Development of the HadCRUH dataset. J. Clim. 21, 5364–5383 (2008).
 6. A. J. Simmons, K. M. Willett, P. D. Jones, P. W. Thorne, D. P. Dee, Low-frequency variations
b0 þ b1 t þ e; t≤a in surface atmospheric humidity, temperature, and precipitation: Inferences from
y¼ ð25Þ
b0 þ b1 t þ b2 ðt  aÞ þ e; t > a reanalyses and monthly gridded observational data sets. J. Geophys. Res. 115, D01110
(2010).
7. K. M. Willett, R. J. H. Dunn, P. W. Thorne, S. Bell, M. de Podesta, D. E. Parker, P. D. Jones,
where y is VPD, NDVI, LAI, GPP, and oceanic evaporation (Eocean); t is C. N. Williams Jr., HadISDH land surface multi-variable humidity and temperature record
for climate monitoring. Clim. Past 10, 1983–2006 (2014).
the year; a is the estimated TP of the time series, defining the timing of a
8. S. M. Vicente-Serrano, R. Nieto, L. Gimeno, C. Azorin-Molina, A. Drumond, A. El Kenawy,
trend change; b0, b1, and b2 are the regression coefficients; and e is the F. Dominguez-Castro, M. Tomas-Burguera, M. Peña-Gallardo, Recent changes of relative
residual of the fit. The investigated variable linear trend is b1 before the humidity: Regional connections with land and ocean processes. Earth Syst. Dynam. 9,
TP and b1 + b2 after the TP. Least squares linear regression was used to 915–937 (2018).
estimate a and other coefficients, where a P value of ≤0.05 was 9. A. L. Fletcher, T. R. Sinclair, L. H. Allen Jr., Transpiration responses to vapor pressure deficit
in well watered ‘slow-wilting’ and commercial soybean. Environ. Exp. Bot. 61, 145–151
considered significant. The 5-year running means were used to quantify (2007).
the TP. 10. A. G. Konings, A. P. Williams, P. Gentine, Sensitivity of grassland productivity to aridity
In addition, the partial correlation method was used to analyze the controlled by stomatal and xylem regulation. Nat. Geosci. 10, 284–288 (2017).
correlation between detrended VPD and detrended NDVI/LAI of five 11. C. M. Restaino, D. L. Peterson, J. Littell, Increased water deficit decreases Douglas fir
growth throughout western US forests. Proc. Natl. Acad. Sci. U.S.A. 113, 9557–9562
satellite-based datasets, excluding the impacts of air temperature, radia-
(2016).
tion, and atmospheric [CO2] concentration. Because of the obvious 12. J. Carnicer, A. Barbeta, D. Sperlich, M. Coll, J. Peñuelas, Contrasting trait syndromes in
transitions of VPD, NDVI, and LAI from 1982 to 2015, we used piece- angiosperms and conifers are associated with different responses of tree growth to
wise linear regression approach to determine the TP and then removed temperature on a large scale. Front. Plant Sci. 4, 409 (2013).
the trends before and after the TP before the correlation analysis. 13. J. Ding, T. Yang, Y. Zhao, D. Liu, X. Wang, Y. Yao, S. Peng, T. Wang, S. Piao, Increasingly
important role of atmospheric aridity on Tibetan alpine grasslands. Geophys. Res. Lett. 45,
2852–2859 (2018).
14. A. P. Williams, C. D. Allen, A. K. Macalady, D. Griffin, C. A. Woodhouse, D. M. Meko,
SUPPLEMENTARY MATERIALS T. W. Swetnam, S. A. Rauscher, R. Seager, H. D. Grissino-Mayer, J. S. Dean, E. R. Cook,
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/ C. Gangodagamage, M. Cai, N. G. McDowell, Temperature as a potent driver of regional
content/full/5/8/eaax1396/DC1 forest drought stress and tree mortality. Nat. Clim. Chang. 3, 292–297 (2013).
Fig. S1. Five-year moving average of VPD. 15. D. B. Lobell, M. J. Roberts, W. Schlenker, N. Braun, B. B. Little, R. M. Rejesus, G. L. Hammer,
Fig. S2. Spatial distributions of the difference of VPD trends (kPa year−1) before and after Greater sensitivity to drought accompanies maize yield increase in the U.S. Midwest.
TP years. Science 344, 516–519 (2014).
Fig. S3. Spatial pattern of VPD changes between 1982–1986 and 2011–2015 derived from CRU dataset. 16. M. Rhein, S. R. Rintoul, S. Aoki, E. Campos, D. Chambers, R. A. Feely, S. Gulev,
Fig. S4. Interannual variability of SVP (red dots and lines), AVP (green dots and lines), and air G. C. Johnson, S. A. Josey, A. Kostianoy, C. Mauritzen, D. Roemmich, L. D. Talley, F. Wang,
temperature (Ta; blue dots and lines) derived from four datasets. Observations: Ocean, in Climate Change 2013: The Physical Science Basis: Working Group I
Fig. S5. Global mean LAI and linear trends during 1982–2015. Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate
Fig. S6. Differences of LAI trends over the globally vegetated areas between before and after Change, T. F. Stocker, D. Qin, G.-K. Plattner, M. Tignor, S. K. Allen, J. Boschung, A. Nauels,
TP years. Y. Xia, V. Be, P. M. Midgley, Eds. (Cambridge Univ. Press, 2013) ch. 3.

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 10 of 12

SCIENCE ADVANCES | RESEARCH ARTICLE

17. K. A. Novick, D. L. Ficklin, P. C. Stoy, C. A. Williams, G. Bohrer, A. C. Oishi, S. A. Papuga, 39. Z. C. Zhu, J. Bi, Y. Pan, S. Ganguly, A. Anav, L. Xu, A. Samanta, S. Piao, R. R. Nemani,
P. D. Blanken, A. Noormets, B. N. Sulman, R. L. Scott, L. Wang, R. P. Phillips, The increasing R. B. Myneni, Global data sets of vegetation leaf area index (LAI)3g and fraction of
importance of atmospheric demand for ecosystem water and carbon fluxes. Nat. Clim. Chang. photosynthetically active radiation (FPAR)3g derived from global inventory modeling
6, 1023–1027 (2016). and mapping studies (GIMMS) normalized difference vegetation index (NDVI3g) for the
18. Q. Fu, S. Feng, Responses of terrestrial aridity to global warming. J. Geophys. Res. Atmos. period 1981 to 2011. Remote Sens. 5, 927–948 (2013).
119, 7863–7875 (2014). 40. M. Claverie, J. Matthews, E. F. Vermote, C. O. Justice, A 30+ year AVHRR LAI and FAPAR
19. K. E. Trenberth, L. Smith, T. Qian, A. Dai, J. Fasullo, Estimates of the global water budget climate data record: Algorithm description and validation. Remote Sens. 8, 263 (2016).
and its annual cycle using observational and model data. J. Hydrometeorol. 8, 758–769 41. D. A. Sims, A. F. Rahman, V. D. Cordova, D. D. Baldocchi, L. B. Flanagan, A. H. Goldstein,
(2007). D. Y. Hollinger, L. Misson, R. K. Monson, H. P. Schmid, S. C. Wofsy, L. Xu, Midday values
20. L. S. Yu, R. A. Weller, Objectively analyzed air-sea heat fluxes for the global ice-free oceans of gross CO2 flux and light use efficiency during satellite over passes can be used to
(1981–2005). Bull. Am. Meteorol. Soc. 88, 527–539 (2007). directly estimate eight-day mean flux. Agric. For. Meteorol. 131, 1–12 (2005).
21. M. P. Byrne, P. A. O’Gorman, Understanding decreases in land relative humidity with 42. H. D. Grissino-Mayer, H. C. Fritts, The international tree-ring data bank: An enhanced
global warming: Conceptual model and GCM simulations. J. Clim. 29, 9045–9061 global database serving the global scientific community. The Holocene 7, 235–238 (1997).
(2016). 43. Cooperative Global Atmospheric Data Integration Project, Multi-laboratory compilation
22. R. D. Koster, M. J. Suarez, M. Heiser, Variance and predictability of precipitation at of synchronized and gap-filled atmospheric carbon dioxide records for the period
seasonal-to-interannual timescales. J. Hydrometeorol. 1, 26–46 (2000). 1979–2012 (obspack_co2_1_GLOBALVIEW-CO2_2013_v1.0.4_2013-12-23) (NOAA Global
23. R. D. Koster, S. D. Schubert, M. J. Suarez, Analyzing the concurrence of meteorological Monitoring Division, 2013); https://doi.org/10.3334/OBSPACK/1002.
droughts and warm periods, with implications for the determination of evaporative 44. G. D. Farquhar, S. von Caemmerer, J. A. Berry, A biochemical model of photosynthetic
regime. J. Clim. 22, 3331–3341 (2009). CO2 assimilation in leaves of C3 species. Planta 149, 78–90 (1980).
24. M. Jung, M. Reichstein, P. Ciais, S. I. Seneviratne, J. Sheffield, M. L. Goulden, G. Bonan, 45. G. J. Collatz, J. T. Ball, C. Grivet, J. A. Berry, Physiological and environmental regulation of
A. Cescatti, J. Chen, R. de Jeu, A. J. Dolman, W. Eugster, D. Gerten, D. Gianelle, N. Gobron, stomatal conductance, photosynthesis, and transpiration: A model that includes a

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

J. Heinke, J. Kimball, B. E. Law, L. Montagnani, Q. Mu, B. Mueller, K. Oleson, D. Papale, laminar boundary layer. Agric. For. Meteorol. 54, 107–136 (1991).
A. D. Richardson, O. Roupsard, S. Running, E. Tomelleri, N. Viovy, U. Weber, C. Williams, 46. B. E. Medlyn, R. A. Duursma, D. Eamus, D. S. Ellsworth, I. C. Prentice, C. V. M. Barton,
E. Wood, S. Zaehle, K. Zhang, Recent decline in the global land evapotranspiration trend K. Y. Crous, P. De Angelis, M. Freeman, L. Wingate, Reconciling the optimal and
due to limited moisture supply. Nature 467, 951–954 (2010). empirical approaches to modelling stomatal conductance. Glob. Chang. Biol. 17, 2134–2144
25. C. Jiang, Y. Ryu, H. Fang, R. Myneni, M. Claverie, Z. Zhu, Inconsistencies of interannual (2011).
variability and trends in long-term satellite leaf area index products. Glob. Chang. Biol. 23, 47. I. C. Prentice, N. Dong, S. M. Gleason, V. Maire, I. J. Wright, Balancing the costs of carbon
4133–4146 (2017). gain and water transport: Testing a new theoretical framework for plant functional
26. W. Yuan, S. Liu, G. Yu, J.-M. Bonnefond, J. Chen, K. Davis, A. R. Desai, A. H. Goldstein, ecology. Ecol. Lett. 17, 82–91 (2014).
D. Gianelle, F. Rossi, A. E. Suyker, S. B. Verma, Global estimates of evapotranspiration 48. T. F. Keenan, I. C. Prentice, J. G. Canadell, C. A. Williams, H. Wang, M. Raupach, G. J. Collatz,
and gross primary production based on MODIS and global meteorology data. Recent pause in the growth rate of atmospheric CO2 due to enhanced terrestrial carbon
Remote Sens. Environ. 114, 1416–1431 (2010). uptake. Nat. Commun. 7, 13428 (2016).
27. W. K. Smith, S. C. Reed, C. C. Cleveland, A. P. Ballantyne, W. R. L. Anderegg, W. R. Wieder, 49. L. Korson, W. Drost-Hansen, F. J. Millero, Viscosity of water at various temperatures.
Y. Y. Liu, S. W. Running, Large divergence of satellite and Earth system model estimates of J. Phys. Chem. 73, 34–39 (1969).
global terrestrial CO2 fertilization. Nat. Clim. Chang. 6, 306–310 (2016). 50. W. Yuan, S. Liu, G. Zhou, G. Zhou, L. L. Tieszen, D. Baldocchi, C. Bernhofer, H. Gholz,
28. N. McDowell, W. T. Pockman, C. D. Allen, D. D. Breshears, N. Cobb, T. Kolb, J. Plaut, A. H. Goldstein, M. L. Goulden, D. Y. Hollinger, Y. Hu, B. E. Law, P. C. Stoy, T. Vesala,
J. Sperry, A. West, D. G. Williams, E. A. Yepez, Mechanisms of plant survival and mortality S. C. Wofsy; AmeriFlux collaborators, Deriving a light use efficiency model from eddy
during drought: Why do some plants survive while others succumb to drought? covariance flux data for predicting daily gross primary production across biomes.
New Phytol. 178, 719–739 (2008). Agric. For. Meteorol. 143, 189–207 (2007).
29. C. D. Allen, A. K. Macalady, H. Chenchouni, D. Bachelet, N. McDowell, M. Vennetier, 51. S. Piao, S. Sitch, P. Ciais, P. Friedlingstein, P. Peylin, X. Wang, A. Ahlström, A. Anav,
T. Kitzberger, A. Rigling, D. D. Breshears, E. H. T. Hogg, P. Gonzalez, R. Fensham, Z. Zhang, J. G. Canadell, N. Cong, C. Huntingford, M. Jung, S. Levis, P. E. Levy, J. Li, X. Lin,
J. Castro, N. Demidova, J.-H. Lim, G. Allard, S. W. Running, A. Semerci, N. Cobb, M. R. Lomas, M. Lu, Y. Luo, Y. Ma, R. B. Myneni, B. Poulter, Z. Z. Sun, T. Wang, N. Viovy,
A global overview of drought and heat-induced tree mortality reveals emerging climate S. Zaehle, N. Zeng, Evaluation of terrestrial carbon cycle models for their response to
change risks for forests. For. Ecol. Manag. 259, 660–684 (2010). climate variability and to CO2 trends. Glob. Chang. Biol. 19, 2117–2132 (2013).
30. P. S. de Cárcer, Y. Vitasse, J. Peñuelas, V. E. J. Jassey, A. Buttler, C. Signarbieux, 52. S. Sitch, P. Friedlingstein, N. Gruber, S. D. Jones, G. Murray-Tortarolo, A. Ahlström,
Vapor–pressure deficit and extreme climatic variables limit tree growth. Glob. Chang. Biol. S. C. Doney, H. Graven, C. Heinze, C. Huntingford, S. Levis, P. E. Levy, M. Lomas, B. Poulter,
24, 1108–1122 (2018). N. Viovy, S. Zaehle, N. Zeng, A. Arneth, G. Bonan, L. Bopp, J. G. Canadell, F. Chevallier,
31. I. Harris, P. D. Jones, T. J. Osborn, D. H. Lister, Updated high-resolution grids of monthly P. Ciais, R. Ellis, M. Gloor, P. Peylin, S. L. Piao, C. Le Quéré, B. Smith, Z. Zhu, R. Myneni,
climatic observations—The CRU TS3.10 Dataset. Int. J. Climatol. 34, 623–642 (2014). Recent trends and drivers of regional sources and sinks of carbon dioxide. Biogeosciences
32. K. M. Willett, C. N. Williams Jr., R. J. H. Dunn, P. W. Thorne, S. Bell, M. de Podesta, 12, 653–679 (2015).
P. D. Jones, D. E. Parker, HadISDH: An updateable land surface specific humidity product 53. Y. P. Wang, R. M. Law, B. Pak, A global model of carbon, nitrogen and phosphorus cycles
for climate monitoring. Clim. Past 9, 657–677 (2013). for the terrestrial biosphere. Biogeosciences 7, 2261–2282 (2010).
33. D. P. Dee, S. M. Uppala, A. J. Simmons, P. Berrisford, P. Poli, S. Kobayashi, U. Andrae, 54. M. A. Arain, F. Yuan, T. A. Black, Soil-plant nitrogen cycling modulated carbon exchanges
M. A. Balmaseda, G. Balsamo, P. Bauer, P. Bechtold, A. C. M. Beljaars, L. van de Berg, in a western temperate conifer forest in Canada. Agric. For. Meteorol. 140, 171–192
J. Bidlot, N. Bormann, C. Delsol, R. Dragani, M. Fuentes, A. J. Geer, L. Haimberger, (2006).
S. B. Healy, H. Hersbach, E. V. Hólm, L. Isaksen, P. Kållberg, M. Köhler, M. Matricardi, 55. D. M. Lawrence, K. W. Oleson, M. G. Flanner, P. E. Thornton, S. C. Swenson, P. J. Lawrence,
A. P. McNally, B. M. Monge-Sanz, J.-J. Morcrette, B.-K. Park, C. Peubey, P. de Rosnay, X. Zeng, Z.-L. Yang, S. Levis, K. Sakaguchi, G. B. Bonan, A. G. Slater, Parameterization
C. Tavolato, J.-N. Thépaut, F. Vitart, The ERA-Interim reanalysis: Configuration improvements and functional and structural advances in version 4 of the community
and performance of the data assimilation system. Q. J. R. Meteorol. Soc. 137, 553–597 land model. J. Adv. Model. Earth Syst. 3, M03001 (2011).
(2011). 56. A. K. Jain, X. J. Yang, Modeling the effects of two different land cover change data sets on
34. R. H. Reichle, The MERRA-Land Data Product. GMAO Office Note No. 3 (Version 1.2) the carbon stocks of plants and soils in concert with CO2 and climate change.
(2012), vol. 38; https://gmao.gsfc.nasa.gov/reanalysis/MERRA/. Glob. Biogeochem. Cycles 19, GB2015 (2005).
35. A. L. Buck, New equations for computing vapor pressure and enhancement factor. 57. T. J. Raddatz, C. H. Reick, W. Knorr, J. Kattge, E. Roeckner, R. Schnur, K.-G. Schnitzler,
J. Appl. Meteorol. 20, 1527–1532 (1981). P. Wetzel, J. Jungclaus, Will the tropical land biosphere dominate the climate–carbon
36. J. E. Pinzon, C. J. Tucker, A non-stationary 1981–2012 AVHRR NDVI3g time series. Remote cycle feedback during the twenty-first century? Clim. Dyn. 29, 565–574 (2007).
Sens. 6, 6929–6960 (2014). 58. M. J. Best, M. Pryor, D. B. Clark, G. G. Rooney, R. L. H. Essery, C. B. Ménard, J. M. Edwards,
37. Z. Q. Xiao, S. L. Liang, J. D. Wang, Y. Xiang, X. Zhao, J. L. Song, Long time-series Global M. A. Hendry, A. Porson, N. Gedney, L. M. Mercado, S. Sitch, E. Blyth, O. Boucher, P. M. Cox,
Land Surface Satellite (GLASS) leaf area index product derived from MODIS and AVHRR C. S. B. Grimmond, R. J. Harding, The joint UK land environment simulator (JULES), model
data. IEEE Trans. Geosci. Remote Sens. 54, 5301–5318 (2016). description—Part 1: Energy and water fluxes. Geosci. Model Dev. 4, 677–699 (2011).
38. Y. Liu, R. Liu, J. M. Chen, Retrospective retrieval of long-term consistent global leaf area 59. B. Smith, I. C. Prentice, M. T. Sykes, Representation of vegetation dynamics in the
index (1981–2011) from combined AVHRR and MODIS data. J. Geophys. Res. Biogeosci. modelling of terrestrial ecosystems: Comparing two contrasting approaches within
117, G04003 (2012). European climate space. Glob. Ecol. Biogeogr. 10, 621–637 (2001).

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 11 of 12

SCIENCE ADVANCES | RESEARCH ARTICLE

60. K. M. Keller, S. Lienert, A. Bozbiyik, T. F. Stocker, O. V. Churakova, D. C. Frank, S. Klesse, S.Y. contributed to the interpretation of the results and to the writing of the paper. Competing
C. D. Koven, M. Leuenberger, W. J. Riley, M. Saurer, R. Siegwolf, R. B. Weigt, F. Joos, 20th interests: The authors declare that they have no competing interests. Data and materials
century changes in carbon isotopes and water-use efficiency: Tree-ring-based evaluation availability: This work used EC data acquired and shared by the FLUXNET community, including
of the CLM4.5 and LPX-Bern models. Biogeosciences 14, 2641–2673 (2017). these networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont,
61. G. Krinner, N. Viovy, N. de Noblet-Ducoudré, J. Ogée, J. Polcher, P. Friedlingstein, P. Ciais, ChinaFlux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux-TERN, TCOS-Siberia,
S. Sitch, I. C. Prentice, A dynamic global vegetation model for studies of the coupled and USCCC. The ERA-Interim reanalysis data are provided by ECMWF and processed by LSCE.
atmosphere-biosphere system. Glob. Biogeochem. Cycles 19, GB1015 (2005). The FLUXNET EC data processing and harmonization was carried out by the European Fluxes
62. A. Ito, The regional carbon budget of East Asia simulated with a terrestrial Database Cluster, AmeriFlux Management Project, and Fluxdata project of FLUXNET, with the
ecosystem model and validated using AsiaFlux data. Agric. For. Meteorol. 148, support of CDIAC and ICOS Ecosystem Thematic Center, and the OzFlux, ChinaFlux and
738–747 (2008). AsiaFlux offices. All data needed to evaluate the conclusions in the paper are present in the
63. C. D. Keeling, T. P. Whorf, Atmospheric CO2 records from sites in the SIO sampling paper and/or the Supplementary Materials. Additional data related to this paper may be
network, in Trends: A Compendium of Data on Global Change (Carbon Dioxide Information requested from the authors.
Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, 2005).
64. J. D. Toms, M. L. Lesperance, Piecewise regression: A tool for identifying ecological Submitted 25 February 2019
thresholds. Ecology 84, 2034–2041 (2003). Accepted 10 July 2019
Published 14 August 2019
Acknowledgments 10.1126/sciadv.aax1396
Funding: This study was supported by the National Basic Research Program of China
(2016YFA0602701), National Youth Top-notch Talent Support Program (2015-48), and Changjiang Citation: W. Yuan, Y. Zheng, S. Piao, P. Ciais, D. Lombardozzi, Y. Wang, Y. Ryu, G. Chen,
Young Scholars Programme of China (Q2016161). Author contributions: W.Y. designed the W. Dong, Z. Hu, A. K. Jain, C. Jiang, E. Kato, S. Li, S. Lienert, S. Liu, J. E.M.S. Nabel, Z. Qin,
research. Y.Z., S.P., G.C., and S. Li performed the analysis. W.Y., P.C., D.L., Y.W., and Y.R. drafted T. Quine, S. Sitch, W. K. Smith, F. Wang, C. Wu, Z. Xiao, S. Yang, Increased atmospheric

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

the paper. W.D., Z.H., A.K.J., C.J., E.K., S. Lie, S. Liu, J.E.M.S.N., Z.Q., T.Q., S.S., W.K.S., F.W., C.W., Z.X., and vapor pressure deficit reduces global vegetation growth. Sci. Adv. 5, eaax1396 (2019).

Yuan et al., Sci. Adv. 2019; 5 : eaax1396 14 August 2019 12 of 12

Increased atmospheric vapor pressure deficit reduces global vegetation growth
Wenping Yuan, Yi Zheng, Shilong Piao, Philippe Ciais, Danica Lombardozzi, Yingping Wang, Youngryel Ryu, Guixing Chen,
Wenjie Dong, Zhongming Hu, Atul K. Jain, Chongya Jiang, Etsushi Kato, Shihua Li, Sebastian Lienert, Shuguang Liu, Julia
E.M.S. Nabel, Zhangcai Qin, Timothy Quine, Stephen Sitch, William K. Smith, Fan Wang, Chaoyang Wu, Zhiqiang Xiao and
Song Yang

Sci Adv 5 (8), eaax1396.

DOI: 10.1126/sciadv.aax1396

Downloaded from http://advances.sciencemag.org/ on November 13, 2019

ARTICLE TOOLS http://advances.sciencemag.org/content/5/8/eaax1396

SUPPLEMENTARY http://advances.sciencemag.org/content/suppl/2019/08/12/5.8.eaax1396.DC1
MATERIALS

REFERENCES This article cites 59 articles, 2 of which you can access for free
http://advances.sciencemag.org/content/5/8/eaax1396#BIBL

PERMISSIONS http://www.sciencemag.org/help/reprints-and-permissions

Use of this article is subject to the Terms of Service

Science Advances (ISSN 2375-2548) is published by the American Association for the Advancement of Science, 1200 New
York Avenue NW, Washington, DC 20005. The title Science Advances is a registered trademark of AAAS.
Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of
Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial
License 4.0 (CC BY-NC).