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

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²) when a predictor variable was removed from the multiple regression was calculated. The change in R² is also called partial R², and it represents the amount of unique variance that each predictor variable explains. Thereby, the partial R² 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² and root mean square error (RMSE).

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 well-watered 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):

$$\begin{equation}E{T_p} - E{T_w} = b\left( {E{T_w} - E{T_a}} \right)\end{equation} \tag{ 1 }$$

Therefore, ET_a can be estimated based on equation (1) as:

$$\begin{equation}E{T_a} = \frac{{\left( {1 + b} \right)E{T_w} - E{T_p}}}{b}\end{equation} \tag{ 2 }$$

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

$$\begin{equation}b = \frac{1}{\gamma }\frac{{e_s^* - e_w^*}}{{{T_s} - {T_w}}}\end{equation} \tag{ 3 }$$

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⁻³P (Allen et al 1998) is the psychrometric constant (kPa °C⁻¹) 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:

$$\begin{equation}{e^*}\left( T \right) = 0.6108\exp \left( {\frac{{17.27T}}{{T + 237.3}}} \right)\end{equation} \tag{ 4 }$$

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

$$\begin{equation}\frac{{\partial H}}{{\partial LE}} = - 1 = \gamma \frac{{{T_w} - {T_s}}}{{e_w^* - {e_s}}}\end{equation} \tag{ 5 }$$

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⁻¹) is estimated using the approach proposed by Penman (1948):

$$\begin{equation}E{T_p} = \frac{{0.408\Delta \left( {{T_s}} \right)\left( {{R_n} - G} \right)}}{{\Delta \left( {{T_s}} \right) + \gamma }} + \frac{{\gamma f\left( {{U_2}} \right)e_s^*\left( {1 - RH} \right)}}{{\Delta \left( {{T_s}} \right) + \gamma }}\end{equation} \tag{ 6 }$$

where Δ(T_s) is the slope of the saturation vapor pressure curve at air temperature T_s (kPa °C⁻¹) (Allen et al 1998); R_n and G are net radiation and soil heat flux (MJ m⁻² d⁻¹); RH is relative humidity; f(U₂) (mm d⁻¹ kPa⁻¹) is an empirical wind function to characterize the aerodynamic conditions of the atmosphere (Penman 1948), and it is defined as:

$$\begin{equation}f\left( {{U_2}} \right) = 2.6\left( {1 + 0.54{U_2}} \right)\end{equation} \tag{ 7 }$$

where U₂ is the wind speed (m s⁻¹) at a height of 2 m.

Based on the concept of equilibrium evapotranspiration defined by Bouchet (1963), ET_w (mm d⁻¹) 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:

$$\begin{equation}E{T_w} = \frac{{0.408\left( {{R_n} - G} \right)}}{{1 + {\beta _{wet}}}}\end{equation} \tag{ 8 }$$

Similar to equation (5), the β_wet for the hypothetical wet environment from temperature T_w to T_s is:

$$\begin{equation}{\beta _{wet}} = \frac{{\partial H}}{{\partial LE}} = \gamma \frac{{{T_s} - {T_w}}}{{e_s^* - e_w^*}}\end{equation} \tag{ 9 }$$

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

The FLUXNET2015 Dataset (see section 2 in Supplementary Information for details) was used to evaluate the influence of RH on ET, and the performance of the improved CR formulation at site-level and weekly time scale. ERA-Interim meteorological data (Berrisford et al 2011) have been extensively used and verified for accuracy across the globe (Albergel et al 2013, Dee et al 2014, Vuichard and Papale 2015). To estimate global terrestrial ET (ET_{CR_ERA-Interim}), input data to the improved CR formulation were obtained from gridded monthly ERA-Interim data over the 1979–2017 period (see section 3 in Supplementary Information for details). The spatial resolution of ERA-Interim data is approximately 0.7 degrees in longitude and latitude. Since such coarse spatial scale did not completely match the site-level FLUXNET2015 data, a dataset of the multi-year annual average basin-level ET was obtained and used as an additional independent observation to further verify the accuracy of ET_{CR_ERA-Interim} at a different spatial scale. In this study, 236 basins with long time series data (≥9 years) were selected (figure S3, see section 4 in supplementary information for details) (available online at stacks.iop.org/ERL/15/094066/mmedia). To reduce the uncertainty in selecting the soil moisture dataset, three widely used datasets were analysed, including the root-zone ERA-Interim soil moisture dataset (Berrisford et al 2011), the root-zone soil moisture dataset from Global Land Data Assimilation System (GLDAS) (Beaudoing and Rodell 2016) and the satellite-based surface soil moisture dataset released by European Space Agency, Climate Change Initiative (ESA-CCI) (Liu et al 2012, Dorigo et al 2017) (see section 5 in supplementary information for details).

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.

Zoom In Zoom Out Reset image size

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.56 mm d⁻¹ and 1.7 mm d⁻¹ to 0.89 mm d⁻¹ 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 (Rodell et al 2004, Beaudoing and Rodell 2016) and Global Land Evaporation Amsterdam Model (GLEAM) (Miralles et al 2011, Martens et al 2017), results show that the basin-level ET estimated by ET_{CR_ERA-Interim} outperforms the other methods by approximately 25.9%, 14.4% and 13.3% (as reductions in RMSE), respectively (figure 2), confirming the overall reliability of the ET_{CR_ERA-Interim} proposed in this study.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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 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₂ is also known to stimulate stomatal closure (Ainsworth and Rogers 2007), to disentangle the role of atmospheric CO₂ 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).

Zoom In Zoom Out Reset image size

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 root-zone (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).

Zoom In Zoom Out Reset image size

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.