Globally Increasing Atmospheric Aridity Over the 21st Century

Abstract

Vapor pressure deficit (VPD) is of great importance to control the land-atmosphere exchange of water and CO₂. Here we use in situ observations to assess the performance of monthly VPD calculated from state-of-the-art data sets including CRU, ERA5, and Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2). We investigate trends in VPD at global scale and for different climatic zones for 1981–2020 and future trends (2021–2100) from Coupled Model Inter-comparison Project phase 6 (CMIP6) outputs. The results show that monthly VPD estimated from CRU, ERA5, and MERRA2 correlated well against in situ estimates from 15,531 World Meteorological Organization stations, with R² ranging between 0.92 and 0.96. Moreover, robust correlations were also found across in situ stations and when analyzing different months separately. During 1981–2020, VPD increased in all climatic zones, with the strongest increase in the arid zone, followed by tropical, temperate, cold and polar zones. CMIP6 simulations show a continuously increasing trend in VPD (0.028 hPa year⁻¹), with the largest increase in the arid zone (0.063 hPa year⁻¹). The magnitudes of trends are found to increase following the magnitude of CO₂ increases in the future emission scenarios. We highlight that atmospheric aridification will continue under global warming, which may pose an increasing threat to terrestrial ecosystems and particularly dryland agricultural systems.

Introduction

Vapor pressure deficit (VPD) is a critical determinant of the atmospheric demand for water vapor (Anderson, 1936). VPD is defined as the difference between saturated vapor pressure (SVP), which is the amount of water vapor that air can hold when it is saturated, and actual vapor pressure (AVP), which is the amount of water vapor that air actually hold at a given condition. According to the Clausius-Clapeyron relation, increasing air temperature generally leads to increases in SVP (Bohren & Albrecht, 2000; Breshears et al., 2013). Changes in AVP can be caused by many factors controlling air temperature and moisture in the air. Decreases in humidity and/or stomatal conductance and/or evapotranspiration and/or soil moisture possibly all lead to a lower increase in AVP than SVP, consequently resulting in non-linear increases in VPD under rising air temperature (Ficklin & Novick, 2017). VPD directly influences plant physiology through regulation of stomatal conductance (Buckley, 2016; Grossiord et al., 2020), vegetation growth (Yuan et al., 2019), productivity (Novick et al., 2016) and ecosystem transpiration (Liu et al., 2020). An increasing trend in VPD has been observed across the terrestrial biosphere (Grossiord et al., 2020), which has been linked to changes in soil moisture (Zhou et al., 2019), the seasonality of leaf phenology (Chen et al., 2021), and the exchange of water and CO₂ between the biosphere and the atmosphere (Novick et al., 2016). Moreover, an increasing VPD can alter the water supply-demand of vegetation with potentially adverse impact on agricultural systems and food security (Lobell et al., 2013). A further increase in VPD is expected in the future as growing global temperatures lead to a larger increase in SVP than AVP (Westra et al., 2013), reducing stomatal conductance and ultimately affecting the carbon and water exchange of terrestrial ecosystems (Ficklin & Novick, 2017).

Our knowledge on VPD and its importance for the terrestrial biosphere has greatly improved over recent decades. However, our understanding of VPD and its role in terrestrial vegetation ecosystems depends on the data quality. The CRU TS (Harris et al., 2020) climate data are developed using a global network of weather station observations and the ERA5 (Hersbach et al., 2020) and Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2) (Gelaro et al., 2017) data sets are the current state-of-the-art global climate reanalysis data (Gelaro et al., 2017; Muñoz Sabater, 2019). These products have been widely used to study climate change (Van Houtan et al., 2021), including atmospheric aridity (Ambika & Mishra, 2020). However, no studies have yet evaluated the consistency of the VPD products derived from these data sets, which might be subject to some differences originating from the product-specific handling of the aggregation, assimilation and reanalysis processes that might not scale linearly between the variables used to calculate VPD. These differences could further introduce bias in the spatio-temporal trends in VPD, which could lead to misinterpretations regarding the underlying mechanisms of VPD controlling vegetation functioning.

With the background of an intensified global warming over the 21st century (Wang et al., 2017), a projection of the change in VPD is important to study its spatio-temporal response to an increasing temperature. This is because increases in future VPD can further enhance air-drying that profoundly affects water cycle and threatens ecosystem functioning (carbon uptake) (Yuan et al., 2019). However, spatially explicit predictions of future changes in VPD have not been reported, leaving a gap in our current understanding of the implications for the future functioning of global vegetation ecosystems. The Coupled Model Inter-comparison Project Phase 6 (CMIP6) provides the relevant variables for calculating VPD under different future climate scenarios (O'Neill et al., 2016), enabling the assessment of future changes in VPD. In this study, we first evaluate the performance of VPD derived from the CRU, ERA5 and MERRA2 datasets using in situ data from the global surface summary of day product (GSOD) (Table 1). We then analyze the spatio-temporal variations of VPD at a global scale and at the level of climatic zones for the period 1981–2020 using the ensemble mean of VPD from CRU, ERA5 and MERRA2. Future trends of VPD for the period 2021–2100 are finally analyzed using the outputs of CMIP6 under different climate scenarios, representing different climate change mitigation pathways.

Table 1 Data Sets Used in This Study

Datasets	Property	Temporal coverage	Temporal resolution	Spatial resolution	Variables (units)
GSOD	In situ	1929–2021	Daily	–	T (°F), Td (°F)
ERA5-land	Reanalysis	1981–2020	Monthly	0.1° × 0.1°	T (K), Td (K)
MERRA2	Reanalysis	1980–2020	Monthly	0.5° × 0.625°	T (°C),Td (°C)
CRU TS 4.05	Interpolation	1901–2020	Monthly	0.5° × 0.5°	T (°C), AVP (hPa)
CMIP6	Models	2015–2100	Monthly	Table S1 in Supporting Information S1	T (°C), RH (%)

Note. T = air temperature; Td = dew point temperature, AVP = actual vapor pressure, RH = relative humidity.

Materials and Methods In Situ Data

The in situ data used in this study was the GSOD (Menne et al., 2012), which includes 18 surface meteorological variables derived from synoptic/hourly observations from more than 9000 stations worldwide (Figure S1a in Supporting Information S1). This data summary is based on data exchanged under the World Meteorological Organization World Weather Watch Program. The data are produced by the National Centers for Environmental Information and available from 1929 to present, but the spatio-temporal coverage is most complete after 1973 (Menne et al., 2012). The daily GSOD based climate variables are provided at a daily basis calculated as the mean of hourly records per day. The number of hourly records can vary between different days and stations. We averaged the daily dew point temperature (Td) and air temperature (T) to a monthly scale, which were used to calculate monthly VPD. To reduce the uncertainties caused by insufficient records, we excluded VPD estimates for a given month if the daily T and Td were not available for at least 25 days. The monthly VPD from 1981 to 2020, corresponding to the period for which all gridded datasets ERA5, MERRA2 and CRU are available, were extracted for our analysis.

CRU and Reanalysis Data

CRU. The Climatic Research Unit gridded Time Series 4 (CRU TS 4.05) is produced from an extensive network of weather station observations gridded at 0.5° spatial resolution (Harris et al., 2020). The monthly land air temperature (T) and AVP of the newest CRU data set (CRU TS 4.05) from 1981 to 2020 were used to calculate VPD.

ERA5. The monthly ERA5 data set was developed by the Copernicus Climate Change Service at the European Centre for Medium-Range Weather Forecasts. The data are the fifth generation reanalysis of the global climate (ERA5) (Muñoz-Sabater et al., 2021). Land air temperature (T) and dew point temperature (Td) of ERA5 were used in this study to calculate VPD.

MERRA2. Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2) reanalysis data were generated by NASA's Global Modeling and Assimilation Office (Gelaro et al., 2017). MERRA2 data are available from 1980 to present with a spatial resolution of 0.5° latitude by 0.625° longitude. We extracted monthly T and Td from MERRA2 to calculate VPD.

For the CRU temperature data, monthly temperature and dew point temperature were calculated from averaging daily temperatures, which in turn were calculated from mean hourly temperature during a day.

All of the above climate datasets (CRU, MERRA2, and ERA5) variables (T, Td, and AVP) were resampled to a spatial resolution of 1° using the weighted averaging approach for further analysis, except for the product evaluation against in situ observations, where the original spatial resolution was retained.

CMIP6 Data

We used 14 Earth system models (Table S1 in Supporting Information S1) included in CMIP6 to estimate VPD under different future climate scenarios (Eyring et al., 2016; O'Neill et al., 2016). CMIP6 represents a substantial refinement over CMIP5, including more modeling groups and different socioeconomic factors in the emission scenarios, in additional to greenhouse gases, aerosols and other climateforcing used in CMIP5. The scenarios used in CMIP6 combine Shared Socioeconomic Pathways (SSP) and target radiative forcing levels at the end of the 21st century (Gidden et al., 2019). SSP126, SSP245, SSP370 and SSP585 indicate the socioeconomic development together with CO₂ emissions from the lowest to the highest level (Gidden et al., 2019) representing a continuum from more aggressive mitigation strategies targeting a sustainable path of lower CO₂ emissions (SSP1) to a fossil-fueled development with higher emissions (SSP5). Air temperature and relative humidity were selected to analyze trends in VPD for 2021–2100. All datasets from the 14 models were resampled to a spatial resolution of 1° using cubic convolution. Subsequently, the projected yearly VPD were produced by the mean annual VPD of the 14 models for further analysis.

Ancillary Data

Digital elevation model (DEM) data. The DEM used in this study is generated from the Panchromatic Remote-sensing Instrument for Stereo Mapping, one of sensors carried by Advanced Land Observing Satellite (ALOS). The DEM was produced in 2016 and updated in 2020, with a 30 m spatial resolution (Tadono et al., 2016; Takaku et al., 2020) (Figure S1b in Supporting Information S1).

Climate zones. The map of present climatic zones used in this study is an improved Köppen-Geiger climate classification map provided in a 1/12° spatial resolution (Figure S1c in Supporting Information S1). Multiple independent data sources have been used to maximize the accuracy of the classification of climatic zones (Beck et al., 2018).

Moderate Resolution Imaging Spectroradiometer (MODIS) land cover. The MODIS Land Cover Climate Modeling Grid (MCD12C1) Version 6 data product for 2019 with a spatial resolution of 0.05° was used to delineate major cropland areas (Figure S1d in Supporting Information S1).

We resampled these three ancillary datasets to a spatial resolution of 1° using majority interpolation to match the spatial resolution of all other gridded datasets.

Calculation of VPD

We calculated VPD using the following equations (Yuan et al., 2019): [Image Omitted. See PDF] [Image Omitted. See PDF]where, SVP, AVP, and RH are the SVP (hPa), AVP (hPa), and relative humility (%), respectively. SVP and AVP can be derived from the following equations: [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]where T, Td, Z, Pmst, and Pmsl are the air temperature (°C), dew point temperature (°C), altitude (m), air pressure (hPa), and the air pressure at the mean sea level (1013.25 hPa), respectively.

Normalization of VPD

To study the relative increase of VPD in different climatic zones, the yearly VPD was normalized using the following equation: [Image Omitted. See PDF]where X, X_i, mean, and std are the yearly normalized VPD, yearly original VPD, mean of yearly VPD, and standard deviation of yearly VPD.

Analyses

We extracted monthly VPD values from gridded data (ERA5, MERRA2 and CRU) for pixels covering the locations of the in situ stations. In total, 15,531 pairs of observations were derived for the VPD validation. We first assessed the global agreement between 15531 pairs of VPD using the slope of the linear regression and the coefficient of determination (R²) between in situ VPD observations and gridded data. We assessed monthly VPD at in situ station level for the entire period and as a function of different months to check for possibly differences in agreement between seasons. It should be noted that the original spatial resolution of the gridded VPD data were retained, when validated against in situ observations.

To detect linear trends in VPD at annual scale, we first averaged monthly VPD to annual scale for all the three gridded VPD products (ERA5, MERRA2 and CRU), and subsequently the three gridded annual VPD products were averaged to derived an ensemble average VPD. Linear regression analysis of ensemble average annual VPD against year was thus conducted at the per-pixel level to assess historical (1981–2020) and projected (2021–2100) trends in VPD. A statistical F-test was used to determine pixels with significant increasing or decreasing trend. The statistical significance level of p < 0.05 was applied to retain clear and coherent spatial clusters/patterns of trends, which showed similar patterns when relaxing the p level from 0.01 to 0.1 (Figure S8 in Supporting Information S1). We used the improved Köppen-Geiger Climate zones to study changes in VPD at biome level.

Results Evaluation of VPD

The performance of VPD derived from CRU, ERA5 and MERRA2 was evaluated against in situ observations. The results show that VPD derived from these three products performs very well against in situ observations of VPD, with R² ranging from 0.92 to 0.96 (Figures 1a–1c). ERA5 VPD data, which have the highest spatial resolution amongst the three products, had the best performance, which might be related to smaller miss-match between the spatial footprint of a given grid-cell and the associated point observation of the in situ measurement. All products slightly underestimated the magnitude of VPD as compared to in situ observations, reflected by the slope <1.0 and also confirmed by an averaged VPD of 6.31 hPa from all monthly in situ observations, which was larger than from CRU (5.75 hPa), ERA5 (5.91 hPa) and MERRA2 (5.81 hPa) (Figure 1d).

Figure 1. Evaluation of vapor pressure deficit (VPD) from three climate data set products using in situ observations at monthly scale. Panels (a–c) are the correlations between VPD derived from three ERA5, Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2), CRU and VPD derived from all in situ stations (n = 15531). Panel (d) shows the mean VPD from in situ data and from the corresponding pixel value of three climate datasets, respectively. The green dotted lines in the boxes indicate the mean values.

The performance of VPD derived from CRU, ERA5 and MERRA2 across space were assessed and generally showed a high consistency between the three products and in situ observations (Figure 2). However, lower correlations were observed in tropical areas, which was especially pronounced for the CRU data set (Figures 2a–2c; Figures S2 and S3 in Supporting Information S1). Changes in the global average correlation coefficients between different months of the year were observed, with higher correlations in the summer associated with higher temperature, and lower correlations in the winter associated with lower temperature for all products (Figure 2d and Figure S2 in Supporting Information S1). The patterns of changes in the global correlations agreed well with correlation changes in the Northern Hemisphere, whereas the inverse pattern was observed in the Southern Hemisphere (Figures 2e–2f). The correlation coefficients between ERA5 and in situ observations are the highest for all months, followed by MERRA2 and CRU.

Figure 2. Spatial and temporal validation of vapor pressure deficit (VPD) from three climate data set products against in situ observations. Panels (a–c) are the spatial patterns of the correlation coefficient of monthly VPD between ERA5, Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2) and CRU and in situ stations. Panel (d–f) shows the changes in the correlation coefficients as a function of the month of year across the globe and for the Northern and Southern hemispheres, respectively. The solid lines and shading indicate the correlations and their standard errors, respectively.

Historical Changes in VPD

Large spatial variability was observed using the ensemble of averaged annual VPD for 1981–2020 (Figure 3a). The Sahara desert, the Arabian peninsula, and the Australian desert show the highest values of VPD. For different climate zones, arid areas have the highest VPD (15.77 ± 7.94 hPa), followed by topical (8.81 ± 4.06 hPa), temperate (6.82 ± 3.20 hPa), and cold (2.44 ± 1.51 hPa) zones. The polar zone has the lowest VPD (1.05 ± 1.27 hPa) (Figure 3b). The global variability in VPD is driven primarily by the arid climate zone.

Figure 3. Spatial distribution of mean Vapor pressure deficit (VPD) (ERA5, Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2), and CRU) averaged over 1981–2020. (a) Global spatial pattern of mean VPD. (b) Mean VPD globally and for different climatic zones. The green dotted lines in the boxes indicate the mean values.

Normalized VPD changes from 1981 to 2020 show a similar inter-annual variability at global scale and for different climatic zones (Figure 4a). An increase in VPD is observed globally and for each climate zone, with the highest increase in VPD observed in the arid zone, closely followed by the tropical zone, and the lowest increase in the polar zone (Figure 4b; Figure S4 in Supporting Information S1). Moreover, inter-annual variability of normalized VPD was found to be similar to SVP (Figure 4a; Figure S5 in Supporting Information S1). Trends in VPD (Figure 4c; Figure S6 in Supporting Information S1) show a more pronounced positive relation to trends in temperature (Figure S7 in Supporting Information S1) as compared to the corresponding relation of VPD trends to trends in AVP (Figure S7 in Supporting Information S1). Spatially, by far most areas worldwide show an increasing VPD from 1981 to 2020 (Figure 4c). However, few areas stand out by showing a clear negative trend in VPD, for example, Himalayas and Central America, and areas such as central Canada and Greenland showing a trend of slightly decreasing VPD (Figure 4c).

Figure 4. Temporal changes in mean vapor pressure deficit (VPD) (ERA5, Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2), and CRU) from 1981 to 2020. (a) Temporal change in global VPD as well as for individual climatic zones. (b) Averaged trends in VPD globally and for different climatic zones. The green dotted lines in the boxes indicate the mean values. (c) Spatial distribution of trends in VPD. White areas indicate non-significant trends (p > 0.05).

Projected Changes in VPD

We first tested the agreement between VPD from the ensemble average (ERA5, MERRA2, and CRU) and VPD based on historical multi-model mean CMIP6 data from 1981 to 2014. A good agreement between products was observed (Figure S9 in Supporting Information S1), which reinforces our confidence in the following analysis of projected changes in VPD based on the outputs from CMIP6. Overall, an increase in VPD is observed at the global scale as well as for different climate zones for all different climate scenarios (Figure 5; Figure S10). However, the SSP126 scenario reaches a plateau around 2060, from whereon VPD values remain rather constant. Vapor pressure deficit increases are more pronounced in the other scenarios, with the strongest increase under SSP585 for all climatic zones. Specifically, VPD increases globally by more than 3 hPa under SSP585, but increases by less than 1 hPa under SSP126 from 2021 to 2100. The arid zone shows the highest increase in VPD (6.44 hPa under SSP585), followed by tropical areas (5.24 hPa under SSP585), both of which show considerably larger projected increase in VPD than the other zones: temperate 3.44 hPa under SSP585, cold 1.80 hPa under SSP585, and polar 0.50 hPa under SSP585 from 2021 to 2100.

Figure 5. Projected inter-annual variability of vapor pressure deficit (VPD) 2021–2100 at global land and for different climatic zones (tropical, arid, temperate, cold, and polar) and different climate scenarios (SSP126, SSP, 245, SSP370 and SSP 585). The shaded areas indicate the standard deviation of VPD in each year.

Spatial patterns show significantly increasing trends in VPD over all climate scenarios (Figures 6a–6e), especially under SSP585 (0.028 ± 0.028 hPa year⁻¹) (Figure 6d). The areas with the largest magnitude of increasing VPD trends are located in tropic (0.051 ± 0.028 hPa year⁻¹ under the SSP585 scenario) and arid areas (0.063 ± 0.034 hPa year⁻¹ under the SSP585 scenario). This is in particular the case for the United States, South America, India, and Australia (Figures 5a–5d). Contrastingly, the temperate (0.040 ± 0.021 hPa year⁻¹), cold (0.019 ± 0.011 hPa year⁻¹), and polar areas (0.006 ± 0.004 hPa year⁻¹) show only slightly increasing trends under the SSP585 scenario (Figure 6e). Furthermore, we compared the difference in VPD values between the projected (from 2061 to 2100) and historical estimates (from 1981 to 2020) (Figure 7). VPD in the future period is considerably higher than in the recent past in most climatic zones, but particularly for the arid and tropical zones (Figure 7e), with the largest difference observed in the scenarios of SSP585 (3.22 ± 3.57 hPa year⁻¹ under SSP585 scenario). However, projected VPD values that are lower than the historical estimates are observed in many areas of high elevation, such as the Tibetan Plateau (Asia), Andes (South America), Rocky Mountains (North America), Ethiopian Highlands (East Africa), and Greenland (Figures 7a–7d).

Figure 6. Projected trend in vapor pressure deficit (VPD) from 2021 to 2100. (a–d) Show spatial patterns of the projected trend (p [less than] 0.05) in VPD under SSP126, SSP245, SSP370, SSP585, respectively. White areas denote non-significant trends (p > 0.05). (e) Projected mean trend in VPD (p [less than] 0.05) globally and for different climatic zones under different scenarios. The error bars indicate the 75% (upper) and 25% (lower) quartile of trend in VPD.

Figure 7. Difference between projected vapor pressure deficit (VPD) (mean VPD from 2061 to 2100) and historic VPD (mean VPD from 1981 to 2020). (a–d) Show spatial patterns of the difference between projected VPD and historic VPD under SSP126, SSP245, SSP370, SSP585, respectively. (e) Mean difference in VPD globally and for different climatic zones under different scenarios. The error bars indicate the 75% (upper) and 25% (lower) quartile of the difference in VPD.

Discussion Comparison of VPD Products From Different Datasets

Our findings suggest that VPD extracted from all three datasets (i.e., ERA5, MERRA2 and CRU) generally shows a good agreement with a global network of in situ observations, especially the ERA5 data (R² = 0.95, slope = 0.94) (Figure 1). This is due to the high accuracies of the relevant variables temperature, dew point temperature, and relative humidity, which were used to derive VPD (Gupta et al., 2020; Luo et al., 2020; Tarek et al., 2020). For example, Luo et al. (2020) showed that both relative humidity and air temperature from MERRA2 have good consistencies with the data from 553 radiosonde measurements. Tarek et al. (2020) suggested that the performance of ERA5-based hydrological modeling is equivalent to in situ observations, and more accurate than previous versions of ERA-Interim. However, uncertainties with respect to different number of hourly records of T and Td in ERA5 and MERRA2 to calculate VPD can introduce some uncertainty in the data, even being generated using an advanced assimilation system, with multiple data sources as inputs including in situ and remote sensing observations (Iqbal et al., 2022). By contrast, scaling up the point observations to a gridded data set dependent on the interpolation method (angular-distance weighting) can lead to a relatively lower correlation between in situ observations and CRU based VPD. While the R² value of the CRU validation is in the same range as ERA5 and MERRA2 (Figure 2), it can be seen that the slope/offset is lower/higher in the case of the CRU comparison. When conducting point to pixel comparisons for coarse resolution pixels/gridded data a common phenomenon is that lowest/highest observations recorded in a given point will be averaged out in the pixel/grid representation, as it is less likely that these low/high values will be representative for the entire pixel/grid (Fensholt et al., 2006; Rammig et al., 2018). The negative offset and slope values in ERA5, MERRA2 and CRU regressions indicate that the underestimation applies to the entire range of data points, and it is possible that weather stations represents a slightly biased sample of the land surface (located in clearings with mowed or removed vegetation or pavement, such as at airports and in urban areas with high impervious surface cover). Moreover, relatively lower agreement between in situ data and the gridded/reanalysis datasets observed in tropical areas can be attributed to the sparse distributions of in situ stations. Overall, our results highlight the use of VPD derived from MERRA2 and ERA5 reanalysis data for different domains, such as terrestrial ecology studies focusing on vegetation functional traits in relation to water and carbon fluxes.

Historical and Projected Changes in VPD

A considerable increase of VPD was observed in most areas around the world during the last four decades (Figure 4). The increase in VPD was projected to continue from 2021 to 2100 (Figure 5). This can be primarily attributed to increases in temperature as higher temperatures lead to an increase in SVP according to Clausius-Clapeyron relation (Bohren & Albrecht, 2000; Breshears et al., 2013). However, increases in AVP can be modulated by reduced soil moisture or evapotranspiration (Deng et al., 2020; Jung et al., 2010), which would possibly lead to a higher increase in SVP than AVP. However, a decreasing trend in historic VPD was observed in some areas, such as the Himalaya with persistent ice and snow cover (Figure 4). This is likely caused by warming induced melting of ice and snow which is expected to accelerate evaporation, thereby resulting in increases in AVP. There is an obvious slowdown in VPD in the scenarios of SSP126 and SSP245 from the 2050s and onwards (Figure 5), which is mainly because the simulated warming trend slows down during this period (Fan et al., 2020). In addition, diverging trends in daytime and nighttime VPD are expected due to differences in the temperature during the daytime (higher) and nighttime (lower). The associated differences in evapotranspiration and stomatal conductance can also contribute to such diverging trends. Future studies are thus called for to better understand the implications of different trends in daytime and nighttime VPD.

Implications for Terrestrial Ecosystems

The increasing VPD in the recent past and future indicates an ongoing increasing atmospheric drought globally, which is supposedly co-occurring with water related droughts. This threatens vegetation including agricultural systems, in particular in the drier parts of the word and/or areas with smallholder rainfed agriculture in developing countries, where irrigation is not a viable alternative (D B Lobell et al., 2014; Tack et al., 2015; Zhou et al., 2019). For example, a profoundly increasing VPD was observed in croplands around the world (Figure 8). Moreover, rising VPD will also accelerate tropical tree mortality (Bauman et al., 2022). A high VPD generally leads to a decline in ecosystem productivity (Goodrich et al., 2015) through increasing water stress, independent of rising temperatures (Tack et al., 2015). This is because the stomatal conductance (regulating the loss of water) is reduced at the cost of a decrease in CO₂ exchange (productivity or yield) (Grossiord et al., 2020). However, this negative effect can to some extent be alleviated by elevated CO₂, depending on precipitation gradient (Hsiao et al., 2019). On the other hand, increasing VPD could amplify soil moisture dryness by causing an increase in evapotranspiration, resulting in severe droughts (Farahmand et al., 2021) and even vegetation mortality (McDowell et al., 2018).

Figure 8. Projected trend in vapor pressure deficit (VPD) in clusters of major cropland areas (Figure S1d in Supporting Information S1). (a–d) Show spatial patterns of the projected trend (p [less than] 0.05) in VPD under SSP126, SSP245, SSP370, SSP585, respectively. White areas denote croplands non-significant trends (p [less than] 0.05). (e) Projected mean trend in VPD (p [less than] 0.05) for different cropland areas under different scenarios. Error bars indicate the 75% (upper) and 25% (lower) quartile of trend in VPD.

In summary, our study documented that the reanalysis (i.e., ERA5 and MARRE2) and gridded (i.e., CRU), as well as CMIP6 simulated VPD data provide reliable estimates of historic and future changes in VPD. Increasing VPD is expected in the future and poses a global threat to vegetation growth as we know it today, yet varying with local conditions in relation to environmental and socio-economic resources. More research is needed to gain knowledge on the impact of increasing VPD on carbon fluxes and water cycle at the regional to global scale. One pathway toward such improved understanding could be to focus on improved implementation on the adverse impacts of increasing VPD on vegetation production, crop yields, or stomatal conductance. This would allow for a better understanding of global vegetation changes in the context of climate change and intensified human land management.

Acknowledgments

ZF is funded by the China Scholarship Council (CSC) (Grant 201906410082). WZ and MB are supported by ERC project TOFDRY (Grant 947757). WZ also acknowledges funding from the National Natural Science Foundation of China (Grant 42001349). RF acknowledge support by the Villum Foundation through the project “Deep Learning and Remote Sensing for Unlocking Global Ecosystem Resource Dynamics” (DeReEco) (Project Number 34306).

Conflict of Interest

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

Data Availability Statement

Global surface summary of day product is available from: https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.ncdc:C00516; ERA5 is available from: https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land-monthly-means?tab=overview; CRU is available from: https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.05/; MERRA2 is available from: https://disc.gsfc.nasa.gov/datasets?project=MERRA-2; DEM is available from: https://www.eorc.jaxa.jp/ALOS/en/dataset/aw3d30/aw3d30_e.htm; Köppen-Geiger climate classification is available from: http://www.gloh2o.org/koppen/; MODIS land cover is available from: https://ladsweb.modaps.eosdis.nasa.gov/search/order; CMIP6 data set is available from: https://esgf-node.llnl.gov/search/cmip6/.