Stomatal response to decreased relative humidity constrains the acceleration of terrestrial evapotranspiration

Terrestrial evapotranspiration (ET) is thermodynamically expected to increase with increasing atmospheric temperature; however, the actual constraints on the intensification of ET remain uncertain due to a lack of direct observations. Based on the FLUXNET2015 Dataset, we found that relative humidity (RH) is a more important driver of ET than temperature. While actual ET decrease at reduced RH, potential ET increases, consistently with the complementary relationship (CR) framework stating that the fraction of energy not used for actual ET is dissipated as increased sensible heat flux that in turn increases potential ET. In this study, we proposed an improved CR formulation requiring no parameter calibration and assessed its reliability in estimating ET both at site-level with the FLUXNET2015 Dataset and at basin-level. Using the ERA-Interim meteorological dataset for 1979–2017 to calculate ET, we found that the global terrestrial ET showed an increasing trend until 1998, while the trend started to decline afterwards. Such decline was largely associated with a reduced RH, inducing water stress conditions that triggered stomatal closure to conserve water. For the first time, this study quantified the global-scale implications of changes in RH on terrestrial ET, indicating that the temperature-driven acceleration of the terrestrial water cycle will be likely constrained by terrestrial vegetation feedbacks.


Introduction
On a global scale, approximately two-thirds of the annual precipitation over land originates from terrestrial evapotranspiration (ET) (Oki and Kanae 2006). More than half of the total solar energy absorbed by the land surface is used as latent heat to sustain ET (Trenberth et al 2009). Detecting and attributing variations in terrestrial ET is therefore crucial for improving our understanding of landatmosphere interactions and of the terrestrial water cycle, which is of key importance for water resource management, human activities, wellbeing and the associated carbon cycle (Wang and Dickinson 2012, Burke et al 2018. In a changing climate, rising temperatures increase saturated water vapour pressure at a rate of approximately 6-7%K −1 according to the Clausius-Clapeyron equation. There is evidence of increasing water vapour content at both the surface and upper tropospheric levels (Trenberth et al 2005); thus, numerous previous studies have suggested surface relative humidity (RH: ratio of actual water vapour pressure to saturated water vapour pressure) to be constant (Dai 2006, Willett et al 2008. However, recent studies have indicated that the actual atmospheric water vapour content does not always increase at exactly the same rate as saturated vapour pressure. A sharp decrease in RH has been observed since 2000 from both gridded observational and reanalysis data (Simmons et al 2010, Willett et al 2014. Increasing evidence suggests that the decrease in terrestrial RH are due to the greater warming over land than the ocean (Byrne and O'Gorman 2016) and are associated with a concurrent decline in the soil moisture over vast regions in the terrestrial biosphere (Byrne and O'Gorman 2016, Zhou et al 2019). However, details about the patterns of reduced RH still remain unclear (Vicente-Serrano et al 2018).
Reduced RH implies increased vapour pressure deficit (VPD), which is an important proxy of atmospheric water demand by plants (Mcadam andBrodribb 2015, Yuan et al 2019). Intuitively, we expect that terrestrial evapotranspiration (ET) would increase with increased water demand; however, it has been well established that plants close their stomata to prevent excessive water loss in response to increasing atmospheric VPD (Mcadam and Brodribb 2015). Recent studies have highlighted that a rising VPD greatly limits terrestrial ET in many biomes through stomata regulation (Novick et al 2016, Rigden andSalvucci 2017). Using the ET estimated from the relative humidity at equilibrium (ETRHEQ) method, Rigden and Salvucci (2017) detected a sharp decline in ET from 1998 to 2014 in the continental United States, and they attributed this decrease mostly to the decline in stomatal conductance. However, large amounts of observational data, such as half-hourly vertical profiles of RH, are required to estimate the ETRHEQ, which limits its application.
Reduced RH generally decreases actual ET while increasing potential ET. This relationship is consistent with the complementary relationship (CR) framework between actual and potential ET stating that the fraction of energy not used for evapotranspiration is dissipated as increased sensible heat flux that in turn increases potential ET (Bouchet 1963). Because the CR framework of ET requires no detailed knowledge of the complex processes and interactions between soil, vegetation, and the near-surface atmosphere, CR-based approaches are easily deployable to directly estimate ET (Kahler and Brutsaert 2006, Szilagyi and Jozsa 2008, Huntington et al 2011, Brutsaert 2015. However, the application of CR-based methods was limited by the need of a local parameters calibration, while in this study we developed an improved CR formulation not requiring any calibration. Given that the global-scale constraints of changes in RH on terrestrial ET have not yet been quantified, the objectives of this study are (1) to develop a new method to estimate global terrestrial ET and verify its accuracy; (2) to detect and attribute the influences of changes in RH on terrestrial ET.

Observation data-based evaluation of the influence of RH on ET
Based on the observed ET, net radiation, surface air temperature, wind speed, and RH, data-driven approaches were established to evaluate the influence of RH on ET. A linear relationship between RH and ET was initially analyzed. To exclude the effects of net radiation, surface air temperature and wind speed, the partial correlation between observed RH and ET was calculated with the R package 'ppcor' (Kim 2015). Secondly, the relative importance of factors affecting ET was assessed to evaluate the contribution of RH to the variation in ET. With ET as predicted variable and net radiation, surface air temperature, wind speed, and RH as predictor variables, a multiple regression was built. The change in the coefficient of determination (R 2 ) when a predictor variable was removed from the multiple regression was calculated. The change in R 2 is also called partial R 2 , and it represents the amount of unique variance that each predictor variable explains. Thereby, the partial R 2 in the multiple regression was used to illustrate the relative importance of factors affecting ET. Besides, the relative importance of factors affecting ET was further evaluated based on the permutation-based variable importance in a Random Forest method (Breiman 2001), as detailed in (Strobl et al 2007, Grömping 2009). The Random Forest method was implemented through the R package 'ranger' (Wright and Ziegler 2017). Half of the observed data was randomly selected as calibration data, while the other as verification data. The performances of the multiple regression method and the Random Forest method in estimating ET were evaluated using the metrics of R 2 and root mean square error (RMSE).

Estimation of actual ET by a theoretical derived method
We developed an improved complementary relationship (CR) formulation to estimate actual ET. The key idea behind the CR formulation is that as a surface dries, the fraction of energy not used for actual ET (ET a ) is dissipated as sensible heat flux that increases temperature and, hence, potential ET (ET p ), obtaining a complementary relationship between ET a and ET p (Bouchet 1963, Aminzadeh et al 2016. Implicit to the CR formulation, equilibrium ET (ET w ) is defined as the evaporation from an extensive wellwatered surface where energy is the limiting factor. With limited water availability, the reduction in energy used for ET a thus becomes available as sensible heat flux, thereby increasing ET p . The CR formulation can be expressed as (Kahler and Brutsaert 2006): Therefore, ET a can be estimated based on equation (1) as: where b is a coefficient that depicts the proportion of sensible heat flux that increases ET p and is calculated as follows (Granger 1989, Aminzadeh et al 2016: where e * s is saturated vapour pressure at surface temperature T s and e * w is saturated vapour pressure at a hypothetical wet surface temperature T w , γ = 0.665 × 10 −3 P (Allen et al 1998) is the psychrometric constant (kPa • C −1 ) and P is atmospheric pressure (kPa). Notedly, the surface air temperature (T s and T w ) were used in the absence of available surface temperature data. The saturated vapour pressure (unit of kPa) at temperature T is estimated following Allen et al (1998) as: Assuming that the net radiation remains constant as the surface wetness changes, the change in sensible heat flux (H) is opposite to the change in latent heat flux (LE) based on the energy balance, so that (Granger 1989, Aminzadeh et al 2016): where e s is the actual vapor pressure. T w is the only unknown in equation (5), so it can be solved and further substituted into equation (3) to estimate the parameter b.
In this study, ET p (mm d −1 ) is estimated using the approach proposed by Penman (1948): is the slope of the saturation vapor pressure curve at air temperature T s (kPa • C −1 ) (Allen et al 1998); R n and G are net radiation and soil heat flux (MJ m −2 d −1 ); RH is relative humidity; f (U 2 ) (mm d −1 kPa −1 ) is an empirical wind function to characterize the aerodynamic conditions of the atmosphere (Penman 1948), and it is defined as: where U 2 is the wind speed (m s −1 ) at a height of 2 m. Based on the concept of equilibrium evapotranspiration defined by Bouchet (1963), ET w (mm d −1 ) is mostly a function of available energy and can be approximated with the Bowen ratio of a wet environment (β wet ) based on the energy-budget equation: Similar to equation (5), the β wet for the hypothetical wet environment from temperature T w to T s is: It is worth to note that there is no necessity to calibrate any parameter in this improved CR formulation, whose accuracy was further verified in this study (see section 1 in Supplementary Information for details).

Data sets
The FLUXNET2015

Results and discussion
3.1. Contribution of RH to ET based on observation data 127 eddy covariance sites (table S2) selected from the FLUXNET2015 dataset were used to analyse the contribution of RH to ET. Partial correlation analysis (Kim 2015) indicated that 104 of 127 sites showed a positive correlation between RH and ET when the impacts of net radiation, air temperature, and wind speed were excluded ( figure 1(a)), and the positive correlation dominated in all 12 vegetation cover types or biomes except permanent wetlands where actual ET is close to potential ET ( figure S5). Based on the weekly data from all of the 127 eddy covariance sites, the multiple regression and Random Forest methods were used to evaluate the relative importance of net radiation, temperature, RH, and wind speed in controlling ET. As expected, net radiation resulted in the most important driver of ET representing the energy source for ET itself (figures 1(b) and (c)). However, both the multiple regression and Random Forest methods indicated that RH is a more important driver than temperature (figures 1(b) and (c)) highlighting its importance especially since a sharp decrease in RH has been widely observed at global scale in recent years (Simmons et al 2010, Willett et al 2014. The positive partial correlation ( figure 1(a)) and positive slope of RH vs ET in multiple regression ( figure 1(b)) are consistent with the evidence that plants close their stomata to prevent excessive water loss in response to increasing atmospheric VPD (Mcadam and Brodribb 2015). Because the relationship between RH and ET may be nonlinear, a scatter plot of dimensionless ET and RH was used to illustrate their relationship (see section 6 in the supplementary information for details on the method used to obtain dimensionless ET). The scatter plot of RH in relation to the normalized measured actual ET and estimated potential ET (refer to equation (6) for details on the estimation method) illustrates that the normalized actual ET generally increased with increasing RH, while the normalized potential ET generally decreased with increasing RH (figure 1(d)). This relationship is consistent with the framework of Bouchet (1963), who identified the complementary relationship between actual and potential ET stating that the fraction of energy not used for evapotranspiration is dissipated as increased sensible heat flux that in turn increases potential ET.

Performance of improved CR formulation in estimating ET
No parameter calibration is required in the improved CR formulation, thus strongly facilitating its application on a global scale. The performance of the improved CR formulation was compared with other commonly adopted CR formulations against the FLUXNET2015 Dataset. Results indicated that the RMSE of the estimated actual ET generally decreased from 1.32 mm d −1 , 1.56 mm d −1 and 1.7 mm d −1 to 0.89 mm d −1 when the improved CR formulation was used instead of the other commonly adopted CR formulations (figure S1, see section 1 in supplementary information for details). The calibrated data-driven methods (i.e. multiple regression and Random Forest methods in figures 1(b) and (c)) performed slightly better than the improved CR formulation in estimating ET in terms of RMSE. However, the application of methods requiring parameter calibration is limited in data-scarce regions, hence only the CR formulation not requiring any parameter adjustment was considered in this study to estimate global terrestrial ET. Using the required meteorological forcing data for CR formulations from ERA-Interim data, the global terrestrial ET (ET CR_ERA-Interim ) over the 1979-2017 period was estimated. Compared to other widely used ET data derived from ERA-Interim, GLDAS

Change in RH from 1979-2017
Changes in RH trends were assessed using the ERA-Interim meteorological data. A piecewise linear regression method was used to assess the changes in the RH trends and detect the presence of mutation points (Muggeo 2003). We found that RH exhibited a slightly increasing trend before 1998 but started to decrease sharply afterwards ( figure 3(a)). A decrease in RH was observed over approximately 68% of the land (except Antarctica), especially in northern parts of South America, central parts of Africa and Asia ( figure 3(b)). In contrast with the expectation that reduced RH would mainly occur in arid regions where soil drying can suppress ET and reduce RH, we also observed reductions in humid regions ( figure S6). Byrne and Gorman (2016) stated that this reduction was related to the relatively greater warming over land than ocean. Thus, the specific humidity of air advected from the ocean to land increased less than the increase in saturation-specific humidity over land.

Change in global terrestrial ET from 1979-2017 and its main driver
The spatiotemporal variations in ET from 1979-2017 based on ET CR_ERA-Interim revealed that global terrestrial annual average ET increased before 1998 and slightly decreased after 1998 ( figure 4(a)). A decline in the trend of global terrestrial ET from 1998-2000 until 2008 was also found by Jung et al (2010) using a data-driven method based on FLUXNET data. The approach proposed here further supports a decline in the global average terrestrial ET since 1998 that extended until 2017. At the regional scale, different trends of annual average ET were obtained (figure 4(b)). From 1998 to 2017, approximately half of the land (excluding Antarctica) was dominated by decreased annual ET, and even more regions with reduced increasing rates of annual ET covered approximately 58% of the land surface ( figure 4(c)). The decreased annual ET was more significant in arid regions ( figure S7).
Before 1998, RH remained almost constant (figure 3(a)); thus, the terrestrial ET was generally energy limited, and the ET increased with increasing temperatures ( figure S8). Afterwards, RH decreased implying an increased VPD (figure S9) inducing the closure of plant stomata (Mcadam and Brodribb 2015). Novick et al (2016) found that VPD limits ET to a greater extent than soil water stress in many biomes. Hence, the overall reduction in terrestrial ET after 1998 was considered to be mostly caused by reduced RH, which induced plant stomatal regulation to prevent excessive water loss. Notedly, this explanation is based on leaf-level (micro-scale) dynamic, that may not entirely transfer at canopy level since a fraction of canopy evaporation is made  ET between 1979ET between -1997ET between and 1998ET between -2017. Trends were scaled by the corresponding multi-year annual average actual ET. The dashed lines in (a) indicate that the linear trend was not significant at the 5% significance level, and the shaded areas represent the 95% confidence interval. by evaporation from the soil and/or other wet surfaces (Perez-Priego et al 2018). Therefore, we estimated both canopy conductance (G s ) and leaf-level stomatal conductance (g s ) and verified whether the variation of normalized G s is consistent with g s . By inverting the Penman-Monteith equation (Allen et al 1998), we computed the canopy-level G s (see section 8 in supplementary information for details). Results showed that the G s indeed sharply decreased after 1998 (figure 5). Since increasing atmospheric CO 2 is also known to stimulate stomatal closure (Ainsworth and Rogers 2007), to disentangle the role of atmospheric CO 2 concentration and RH in reducing G s , we applied the linearized optimal stomatal control model (OptiL) (Katul et al 2010) to approximate the relative variations in leaf-level g s (see section 9 in supplementary information for details). It was expected that the relative variations in normalized G s and g s would remain similar when the regulation of G s was exclusively through stomatal control (Rigden and Salvucci 2017). Figure 5 shows that the relative variations in normalized G s and g s are highly correlated, and the reduced G s is largely controlled by the decreased RH (i.e. increased VPD). This result highlights the critical role of the stomatal response to decreased RH in modulating the transfer of moisture to the atmosphere. Nevertheless, atmospheric demand represents only one of the key constraints to stomata regulation. Further observations are needed to quantify the present and future role of other resistances to atmospheric demand: the soil (Carminati and Javaux 2020) and the xylem (Eller et al 2020).

Discussion
Although we have attributed reduced ET to stomatal regulation, we must also consider that wide regions exist on the Earth without or poorly covered by vegetation. ET is dominated by soil evaporation in some arid regions, such as deserts in Sahara and central Asia (Zhang et al 2016). Here, figure 6(a) shows that the global average surface soil moisture also significantly decreased after 1998 based on the data from GLDAS, ERA-Interim and ESA CCI (Liu et al 2012, Dorigo et al 2017, Gruber et al 2017. The reduction in rootzone (or 0-255 cm soil depth) soil moisture occurred widely across the globe, especially in most parts of Africa, Australia and Europe ( figure 6(b)). Thus, soil drying further suppresses both soil evaporation and plant evapotranspiration, and a continued global soil drying is expected in the future (Gu, et al 2019).
The start time of RH reduction was consistent with the beginning of warming hiatus around 1998, while such RH reduction did not terminate with the end of warming hiatus around 2012 (Medhaug et al 2017). According to ET reduction shown here (figure 4), decreased energy exchange as ET would result in increased energy exchange as sensible heat flux, accelerating warming over the land surface (Vogel et al 2017). This mechanism may therefore have contributed to the end of the warming hiatus itself, however such a connection needs to be further analysed and demonstrated. Besides, the continued greater warming over land than the ocean enhances the faster increase in saturation specific humidity over land, leading to further reduced RH (Byrne and O'Gorman 2016), therefore representing a positive feedback between decreased ET and reduced RH. Another positive feedback mechanism exists between reduced RH and soil drying, and the concurrent soil drying and atmospheric aridity (low RH) are found to be greatly exacerbated by land-atmosphere feedbacks (Zhou et al 2019). Reduced RH (i.e. rising VPD) has been considered as a major contributor in recent drought-induced plant mortality (Grossiord et al 2020), making plant adaptation and species change in response to a decline in RH worth of future attention.

Conclusions
Our results indicate that global terrestrial ET is not increasing since 1998 due to the reduction of available water in the soil-vegetation-atmosphere continuum. We found that the decreased atmospheric RH cancelled the increase in the global terrestrial ET observed until 1998. Physically, the decreased RH induced plant water stress and triggered stomatal closure to prevent excessive water loss. Terrestrial ET is the main source of precipitation on the land surface and a key driver of the intensified water cycle under global warming (Huntington 2006, Donat et al 2016). The continuing decrease in terrestrial ET observed here indicates that the ongoing temperature-driven acceleration of the terrestrial water cycle is getting constrained due to terrestrial vegetation feedbacks. Nevertheless, further research is needed to disentangle the complexity of vegetation-soil-atmosphere interactions in the climate system.

Acknowledgments
This work used FLUXNET data acquired and shared by the global FLUXNET community, including AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, China-Flux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux-TERN, Swiss FluxNet, TCOS-Siberia, and USCCC. The ERA-Interim reanalysis data obtained from ECMWF and processed by LSCE. The FLUXNET eddy covariance data processing and harmonization was carried out by the European Fluxes Database Cluster, AmeriFlux Management Project, and Fluxdata project of FLUXNET, with the support of CDIAC and ICOS Ecosystem Thematic Center, and the OzFlux, ChinaFlux and AsiaFlux offices. In detail, US-Prr is supported by JAMSTEC and IARC/UAF collaboration study and Japanese Arctic Research Program (ArCS); Department of Energy Office of Science, Ameriflux Network Management Project Support for UW ChEAS Cluster (US-Los, US-PFa, US-WCr, US-Syv, 2012-present). CPC Global Unified Precipitation data was provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their Web site at https://www.esrl.noaa.gov/psd/data/gridded/. Monthly GRDC Reference Dataset was provided by the Global Runoff Data Centre, 56068 Koblenz, Germany. The multi-satellite surface soil moisture data was provided by ESA CCI SM project. The authors acknowledge contributions from individual researcher who made sincere efforts to gather these data sets and made them available for synthesis work. The work was supported by National Key R&D Program of China (2016YFC0402706, 2016YFC0402710); the National Natural Science Foundation of China (51909057, 51539003, 41901041); China Postdoctoral Science Foundation (2017M610292). MG acknowledges funding by Swiss National Science Foundation project ICOS-CH Phase 2 20FI20_173691. The authors also acknowledge Professor Changliang Shao, Dr. Mika Aurela for their interest and helpful correspondence. Finally, we thank two anonymous reviewers and the Editorial Board Member for their constructive comments and review of the manuscript.

Data availability statements
The data that support the findings of this study are available from the corresponding author upon reasonable request. The FLUXNET2015 dataset is available from the global FLUXNET community website (https://fluxnet.fluxdata.org/data/fluxnet2015dataset/); the multi-satellite surface soil moisture data are available from the ESA CCI SM project website (https://www.esa-soilmoisturecci.org/); ERA-Interim datasets are available from the NCAR Research Data Archive at https://rda.ucar.edu/datasets/ds627.1/; GLDAS datasets are available from the NASA Goddard Earth Sciences (GES) Data and Information Services Center (DISC) at https://disc.gsfc.nasa. gov/datasets?keywords=GLDAS; GLEAM datasets are available from https://www.gleam.eu/.