Land-Atmosphere Coupling Constrains Increases to Potential Evaporation in a Warming Climate: Implications at Local and Global Scales

Abstract

The magnitude and extent of runoff reduction, drought intensification, and dryland expansion under climate change are unclear and contentious. A primary reason is disagreement between global circulation models and current potential evaporation (PE) models for the upper limit of evaporation under warming climatic conditions. An emerging body of research suggests that current PE models including Penman-Monteith and Priestley-Taylor may overestimate future evaporation for non-water-stressed conditions. However, they are still widely used for climatic impact analysis although the underlying physical mechanisms for PE projections remain unclear. Here, we show that current PE models diverge from observed non-water-stressed evaporation across site (>1,500 flux tower site years), watershed (>10,000 watershed-years), and global (25 climate models) scales. By not incorporating land-atmosphere coupling processes, current models overestimate non-water-stressed evaporation and its driving factors for warmer and drier conditions. To resolve this, we introduce a land-atmosphere coupled PE model by extending the Surface Flux Equilibrium theory. The proposed PE model accurately reproduces non-water-stressed evaporation across spatiotemporal scales. We find that terrestrial PE will increase at a similar rate to ocean evaporation but much slower than rates suggested by current PE models. This finding suggests that land-atmosphere coupling moderates continental drying trends. Budyko-based runoff projections incorporating our PE model are well aligned with those from coupled climate simulations, implying that land-atmosphere coupling is key to improving predictions of climatic impacts on water resources. Our approach provides a simple and robust way to incorporate coupled land-atmosphere processes into water management tools.

Introduction

Potential evaporation (PE), which sets the upper limit of evaporation, is widely used as a key constraint for estimating actual evaporation and evapotranspiration, as well as runoff, crop water use, aridity, and drought (Ault, 2020; Vicente-Serrano et al., 2020). PE is operationally understood as the maximum evaporation rate when evaporation is not limited by soil moisture although there are several definitions and forms of PE in the literature (Lhomme, 1997). Among various proposed PE models, Penman-Monteith (PM) and Priestley-Taylor (PT) PE are the most widely used since their underlying definitions and formulations are primarily physically based. As such, these PE equations are widely used to predict and analyze changes in drought, aridity, and water availability in relation to a changing climate (Dai, 2013; Greve et al., 2014; Marvel et al., 2019; McEvoy et al., 2020; Milly & Dunne, 2020; Piemontese et al., 2019; Sheffield et al., 2012; Su et al., 2018; Trenberth et al., 2013; Wang et al., 2018).

According to current PE models, PE under climate warming is projected to increase at a greater rate than precipitation over land, leading to increased aridity (Scheff & Frierson, 2014; Sherwood & Fu, 2014). Therefore, many studies have sought to explore the implications of enhanced land surface drying under climate change. However, several recent studies have demonstrated that calculations based on PE lead to overestimation of non-water-stressed evaporation for warmer and drier future climate conditions (Milly & Dunne, 2016, 2017), resulting in overestimates of actual evapotranspiration, soil drying, and runoff reductions compared to direct projections by climate models (Milly & Dunne, 2016, 2017; Roderick et al., 2015; Yang et al., 2019). More recently, the widely accepted dryland expansion trend under climate change (Berdugo et al., 2020; J. Huang et al., 2016; Overpeck & Udall, 2020) has been questioned (Berg & McColl, 2021; Greve et al., 2019; Keenan et al., 2020; Shi et al., 2021), and debate on the magnitudes of past and future trends in drought has intensified (Berg & Sheffield, 2018; Swann et al., 2016; Tomas-Burguera et al., 2020; Yang et al., 2020). These scientific debates originate to a large extent as a result of PE overestimation. That is, if PE models are used for estimating non-water-stressed evaporation, changes in these PE models could result in the overestimation of the changes in non-water-stressed evaporation under changing climatic conditions (Berg & Sheffield, 2018; Vicente-Serrano et al., 2020). Therefore, there is an urgent need to re-evaluate PE to correctly understand and predict climatic impacts on water resources.

Current physically based PE models assume saturated or near saturated surface conditions, and these land surface conditions are assumed to be stationary. Based on this assumption, the Clausius-Clapeyron relationship (i.e., the relationship between temperature and saturation vapor pressure) can be introduced, resulting in a positive, exponential relationship between temperature and PE. As a result, available energy (AE) partitioning by PT PE, which is based on the equilibrium evaporation concept, is largely controlled by temperature. While PM PE also scales with temperature, a decrease in atmospheric relative humidity (RH) also results in increased values for the PM PE due to the increased vertical gradient of RH from the functionally saturated land surface relative to the overlying atmosphere. Consequently, warmer and drier atmospheric future conditions result in substantially increased PE computed using current PE models.

However, the stationary land surface assumptions in PE models are not necessarily robust. For instance, elevated atmospheric CO₂ concentrations can increase the surface resistance, an empirical parameter in the PM model (Swann et al., 2016; Yang et al., 2019, 2020). To resolve the CO₂ fertilization effect, Yang et al. (2019) recently provided a way to correct the surface resistance based on climate simulation outputs. However, this approach was not based on observational evidence nor physical principles and failed to fully resolve the PE overestimation issue noted in more recent studies (Berg & McColl, 2021; Liu et al., 2022; Vicente-Serrano et al., 2020). Furthermore, if we need to estimate surface resistance for calculating PE as suggested by Yang et al. (2019), there may be no practical advantage of the PE concept over the direct estimation of actual evapotranspiration. Direct predictions of actual evapotranspiration by climate simulation and land surface models utilize semi-empirically parameterized surface resistance to account for physiological responses to changes in atmospheric CO₂ concentrations and RH (or vapor pressure deficit). On the other hand, PE models can be simply executed using only meteorological information without considering the land surface conditions, which provides a practical advantage of the PE concept for some analytical contexts. Therefore, the idea of correcting surface resistance for calculating PE may blur the line between PE models and process-based models by giving up the utility of the PE concept.

Maybe a more fundamental problem of the current PE models stems from a paradoxical assumption itself. Assuming a saturated surface for any given meteorological conditions is a contradiction in terms in that the overlying atmospheric conditions are not independent of land surface wetness due to land-atmosphere coupling processes (Bouchet, 1963; Brutsaert & Stricker, 1979; Kim et al., 2021; McColl et al., 2019; Rigden & Salvucci, 2017; Salvucci & Gentine, 2013). Since the warmer and drier future climate over land is already regulated by soil moisture (Berg et al., 2016; Dirmeyer et al., 2021), considering this apparent trend in current PE models as an increased upper limit of evaporation can lead to a “double-counting” of soil drying (Berg & Sheffield, 2018).

Therefore, the rapid increase in PE suggested for warming climate conditions may in fact be a methodological artifact caused by the structure of current PE models that, in effect, ignore land-atmosphere feedback processes (Berg & Sheffield, 2018). Recognizing this problem, Milly and Dunne (2016) suggested an alternative PE model which simply assumes 80% of AE is used for evaporation for non-water-stressed conditions as they speculated the land-atmosphere coupling effect could be embedded in this empirical approach. Although this approach successfully estimated non-water-stressed evaporation (Maes et al., 2019; Milly & Dunne, 2016), there is currently no clear physical explanation for why this empirical equation works.

Here, we directly assess how current PE models diverge from observed non-water-stressed evaporation, and suggest a physically based alternative PE model that constrains the upper limit of evaporation based on land-atmosphere coupling processes. In the following section, we derive the alternative model and discuss the theory behind it. We then compare our novel PE model with the most commonly used PE models including PM (Allen et al., 1998; Monteith, 1965), PT (Priestley & Taylor, 1972), and an empirical model that calculates PE as proportional to AE (hereafter the Milly-Dunne [MD] PE model (Milly & Dunne, 2016)).

Since applications of PE are broad and of great societal importance, we hierarchically evaluate the three current PE models along with our alternative PE model over a range of scales, and present these results at the field, watershed, and global scales. We first use in-situ field-scale observations from 212 eddy covariance tower sites worldwide contained in the FLUXNET2015 data set representing over 1,500 site-years (Pastorello et al., 2020) to test the performance of PE models in reproducing non-water-stressed evaporation and its sensitivity to temperature and RH. This is followed by watershed scale assessment of model performance compared to water balance observations for 338 US watersheds (Duan et al., 2006) for the 1983–2020 period. We then examine global-scale changes in PE models from a historical reference period to the future period using 25 general circulation models (GCMs) that were included in Coupled Model Intercomparison Project Phase 5 (CMIP5) (Taylor et al., 2012). Through these analyses, we evaluate the causes underlying the differing responses of PE models to changing climatic conditions.

PE Models Alternative PE Model Based on Land-Atmosphere Coupling Theory

Terrestrial evaporation is constrained by soil moisture (supply side) as well as climatic conditions (demand-side). Increasing soil moisture (i.e., supply) increases evaporation until evaporation approaches its maximum rate. At the maximum level of evaporation, additional soil moisture cannot increase evaporation further, and thus evaporation and soil moisture become independent of each other. This transition point is known as critical soil moisture (Denissen et al., 2020). Mathematically, this condition may be expressed as the derivative of evaporation (E) with respect to soil moisture (θ) $$\frac{\partial E}{\partial \theta }=0$$.

The soil moisture control on terrestrial evaporation also interacts with the overlying atmosphere since the vertical gradient of RH from the land surface to the atmosphere simultaneously changes with evaporation (Kim et al., 2021; Salvucci & Gentine, 2013). Soil moisture supply to a dry soil (i.e., soil wetting) increases RH at the land surface (RH₀) at a faster rate than the change in atmospheric RH (i.e., $$\frac{\partial {\mathrm{R}\mathrm{H}}_{0}}{\partial \theta } > \frac{\partial \mathrm{R}\mathrm{H}}{\partial \theta }$$), which leads to increased evaporation (Figure 1). As RH₀ approaches saturation, the increasing rate of RH₀ reduces, becoming $$\frac{\partial {\mathrm{R}\mathrm{H}}_{0}}{\partial \theta }=\frac{\partial \mathrm{R}\mathrm{H}}{\partial \theta }$$ at some point, and thus evaporation reaches its maximum. Therefore, the maximum evaporation should satisfy not only $$\frac{\partial E}{\partial \theta }=0$$ but also $$\frac{\partial {\mathrm{R}\mathrm{H}}_{0}}{\partial \theta }=\frac{\partial \mathrm{R}\mathrm{H}}{\partial \theta }$$ if considering the overlying air as a coupled system with the land surface. This conceptual understanding in Figure 1 is supported by empirical data sets using the Soil Moisture Active Passive (SMAP) satellite soil moisture observations and watershed water balance observations (Figure S1 in Supporting Information S1).

Figure 1. Conceptual framework of two mathematical constraints of the alternative potential evaporation model in Equation 1. The vertical profile of relative humidity (RH) evolves with soil moisture (θ). When soil is dry, the sensitivity of surface relative humidity (RH0) to θ change is larger than the sensitivity of atmospheric RH to changes in θ. Once θ approaches critical soil moisture, the two sensitivities become equivalent, and evaporation reaches its maximum.

We provide a detailed derivation of $$\frac{\partial E}{\partial \theta }$$ in Appendix A and its solution in Appendix B. In this derivation, we assume well-coupled land-atmosphere conditions over large area that are comparable to the emerging Surface Flux Equilibrium (SFE) theory where the near-surface atmosphere is much more sensitive to changes in surface conditions than to changes in atmospheric boundary layer processes (McColl et al., 2019). Also, we assume AE is independent of the soil moisture, although AE is a weak function of temperature, albedo, and emissivity that can be affected by soil moisture. Substituting $$\frac{\partial E}{\partial \theta }=0$$ and $$\frac{\partial {\mathrm{R}\mathrm{H}}_{0}}{\partial \theta }=\frac{\partial \mathrm{R}\mathrm{H}}{\partial \theta }$$ into the derivative $$\frac{\partial E}{\partial \theta }$$ and then assuming $$\frac{\partial {\mathrm{R}\mathrm{H}}_{0}}{\partial \theta } > 0$$ yields [Image Omitted. See PDF]where, $$s\left(=\frac{\partial {q}^{\ast }}{\partial T}\right)$$ is the linearized slope of saturation specific humidity versus temperature, γ is psychrometric constant, λ is the latent heat of vapourization, R_n is net radiation, G is soil heat flux, and R_n − G is AE at the land surface.

Equation 1 represents a formulation for the upper limit of evaporation, and as such it can be considered an alternative PE model. Equation 1 hereafter referred to as the PE-LAC model (i.e., PE based on land-atmosphere coupling perspective) as the underlying assumption of the derivation relies on the well-coupled land-atmosphere system. Our proposed alternative has several unique characteristics. First, in deriving Equation 1, the overlying air is considered as a coupled system with the land. Therefore, the physical meaning of the PE-LAC can be stated as “the maximum evaporation that can be supported by the land surface to produce the current atmospheric conditions.” Since RH in the atmosphere is connected with the land surface dryness in our model, a decrease in RH reduces evaporation in Equation 1, which stands in contrast to the PM model. Finally, while the current PE models require empirical parameters, Equation 1 does not, and as such it does not require any parameter calibration. Fourth, the only required variables for using Equation 1 are standard measurements (air temperature, RH, and AE).

Other PE Formulations

We calculated PM PE based on the FAO reference crop method (Allen et al., 1998). [Image Omitted. See PDF]where, $$s\left(=\frac{\partial {q}^{\ast }}{\partial T}\right)$$ is the linearized slope of saturation specific humidity versus temperature, γ is psychrometric constant, λ is the latent heat of vapourization, R_n is net radiation, G is soil heat flux, ρ is air density, c_p is the specific heat capacity of the air, q* is saturation specific humidity, q_a is air specific humidity, T_a is air temperature. r_s is surface resistance and set to 70 s m⁻¹ based on the FAO method (Allen et al., 1998). r_a is aerodynamic resistance and calculated as $$\frac{208}{{u}_{2}}$$, where u₂ is 2 m height wind speed. u₂ was calculated based on FAO method (Allen et al., 1998) from wind speed and measurement height (z) as $${u}_{2}={u}_{z}\frac{4.87}{\mathrm{ln}(67.8z-5.42)}$$.

We then calculated PT PE as follows [Image Omitted. See PDF]where, $${\alpha }_{\mathrm{P}\mathrm{T}}$$ is PT coefficient and set to 1.26 (Priestley & Taylor, 1972).

Finally, we calculated MD PE as follows [Image Omitted. See PDF]where, $${\alpha }_{\mathrm{M}\mathrm{D}}$$ is MD coefficient and set to 0.8 (Milly & Dunne, 2016).

Theoretical Differences Among PE Models

There are various definitions of PE in the literature, with current PE models relying on different concepts or framings. For example, according to the definition of WMO (2012), PE is defined as the “Quantity of water vapor which could be emitted by a surface of pure water under existing conditions.” Here, the existing conditions indicate current meteorological conditions. This definition is aligned with the definition of the PM style PE by Morton (1983) who highlighted the independency of the overlying air to evaporation in defining PE. Evaporation from open water to the current atmosphere could be effectively calculated using the Penman free water evaporation model (Penman, 1948), and thus it can be considered PE in this definition. Indeed, the Penman equation works well for point-scale evaporation from a free water surface such as pan evaporation observations (Roderick et al., 2007). If the free water surface is replaced by a well-watered vegetative surface with the same definition, the PM reference evapotranspiration model (Allen et al., 1998; Monteith, 1965) could be considered PE (Lhomme, 1997). Based on this definition of PE, the Penman or the PM models are widely considered to represent atmospheric evaporative demand (Vicente-Serrano et al., 2020).

However, some authors do not prefer to use this definition of PE as meteorological conditions are not independent of surface water availability (Bouchet, 1963; Brutsaert & Stricker, 1979). For instance, Brutsaert (2015) argued that PM type models yield “apparent” PE instead of PE as they assume evaporation from the potentially saturated surface into non-potential meteorological conditions (i.e., the existing meteorological conditions and typically unsaturated air). In Brutsaert's theory (which originated from the complementary relationship perspective), PE should be defined as evaporation from an “extensive area” or “large area” of a potentially saturated surface into potential meteorological conditions (saturated or near-saturated air). The PT model (Priestley & Taylor, 1972) is suitable for estimating PE in this definition since this model relies on the equilibrium evaporation which is commonly defined as evaporation from a saturated surface into saturated air (Raupach, 2001). While the PM style models yield PE from a saturated (or near-saturated) surface into air for given meteorological conditions, the PT model yields PE from a saturated surface into near-saturated air due to land-atmosphere feedback (Brutsaert, 2015; Kahler & Brutsaert, 2006).

The proposed PE-LAC model (i.e., Equation 1) is more closely related to the definition of PE for the PT PE model than the definition of PE in the PM PE model in that land-atmosphere coupling is considered. However, there are considerable differences between the PT PE and the proposed PE-LAC. Fundamentally, the PT model could be understood as an adjustment of the equilibrium evaporation, and thus the energy partitioning relies only on temperature because of the relationship between temperature and saturation vapor pressure. On the other hand, the PE-LAC model can be understood as an adjustment of SFE evaporation (McColl et al., 2019) which does not assume a saturated surface, and thus RH plays an important role, in particular the vertical gradient of RH. Kim et al. (2021) showed that evaporation could be expressed as a sum of the SFE evaporation and an additional water vapor flux driven by the vertical gradient of the RH. The proposed PE-LAC model builds from Kim et al. (2021) by providing the upper limit of the RH gradient-driven flux. If one rewrites Equation 1 as $$\mathrm{P}\mathrm{E}=\underset{\begin{array}{c}\mathrm{S}\mathrm{F}\mathrm{E}\\ \mathrm{e}\mathrm{v}\mathrm{a}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\end{array}}{\underbrace{\frac{\mathrm{R}\mathrm{H}s}{\mathrm{R}\mathrm{H}s+\gamma }\frac{{R}_{n}-G}{\lambda }}}+\underset{\begin{array}{c}\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}\\ \mathrm{R}\mathrm{H}\,\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\,\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x}\end{array}}{\underbrace{\frac{\gamma }{\mathrm{R}\mathrm{H}s+\gamma }\frac{{R}_{n}-G}{2\lambda }}}$$, this perspective becomes apparent. Therefore, the PE-LAC model is constrained by the maximum gradient of RH at given atmospheric conditions. If the near-surface air at a random location A is drier than a random location B, the maximum evaporation at location A is lower than the maximum evaporation at location B due to the connection between the land-atmosphere system through the maximum RH gradient invoked by the PE-LAC model. This is the reason why low RH corresponds to a low PE-LAC value under a given temperature and AE.

Because of these various underlying definitions and concepts of PE formulations, it is difficult to evaluate PE models using measurements. However, the PE models should provide the upper limit of actual evapotranspiration in the current atmospheric conditions (Lhomme, 1997), and thus actual evaporation should be close to modeled PE when there is no water stress (Maes et al., 2019; Milly & Dunne, 2016). Therefore, maximum and non-water-stress evaporation observations are used as reference points in this analysis evaluating the PE models. Nevertheless, it should be noted that there may be methodological uncertainty due to various definitions of PE.

Materials and Methods

We hierarchically evaluate our alternative PE model and the three current PE models (i.e., Equations 1–4) over a range of scales, and present these results at the field, watershed, and global scales. All resulting figures were generated by using the R statistical language (R Core Team, 2020). In the following subsections, we describe details of the data set we used.

FLUXNET2015 in Situ Observations

The FLUXNET2015 data set, which includes 212 empirical eddy-covariance flux tower sites around globe representing over 1,500 site years (Pastorello et al., 2020), was used for the site-scale analyses. Latent and sensible heat fluxes, net radiation, soil heat flux, air temperature, RH, wind speed, and barometric pressure were obtained at weekly scales from the FLUXNET2015 data set. We selected data only when all required variables for PE calculations are available. Bowen ratio corrected turbulent heat fluxes suggested by Twine et al. (2000) were used for this analysis following methods employed in previous analysis (Maes et al., 2019). Measurement heights for each site were also retrieved to calculate aerodynamic resistance for the PM model.

We only included data for periods for which the quality control flag indicated more than 80% of the half-hourly data were used for generating the daily or weekly data sets (i.e., measured data or good quality gap-filled data). Also, we filtered out data points when AE (i.e., net radiation minus soil heat flux) was negative or local advection was strong (i.e., negative sensible heat flux) (Maes et al., 2019). Also, data in which surface energy imbalance (AE minus sum of sensible and latent heat flux) was greater than 50 W m⁻² were excluded.

Following a recent study (Maes et al., 2019; Tu et al., 2022), we isolated non-water-stressed conditions by selecting data with evaporative fraction (EF) exceeding 95% of each site's EF distribution. This selection strategy was adopted since soil moisture observations are only available for a few FLUXNET2015 sites, and previous research found no significant difference between soil moisture-based and EF-based criteria (Maes et al., 2019). The isolated non-water-stressed evaporation observations were used as a reference to assess the PE models. PE models were calculated based on Equations 1–4 using weekly averaged observed meteorological variables at the flux towers. It should be noted that the PM model was originally developed for shorter time scales, but the FAO PM method was developed for a daily time scale (Allen et al., 1998), and it is widely used for longer time scales (e.g., weekly and monthly) because of the relative insensitivity of PE change to computational time step (Milly & Dunne, 2016). The bigleaf R package (Knauer et al., 2018) was used for analyzing this flux data set.

US Watersheds

For the watershed-scale analysis, watersheds included in Model Parameter Estimation Experiment (MOPEX) were analyzed (Duan et al., 2006). We selected MOPEX watersheds for which more than 30 years of runoff observations are available for the period of 1983–2020. The total number of watersheds meeting the criterion was 328, resulting in over 10,000 watershed-years of runoff observations. We first calculated uncorrected annual evaporation as the difference between observed precipitation based on the PRISM data set (Daly et al., 1994) and USGS runoff observations on a water-year (i.e., 1 October–30 September) basis. Groundwater storage changes for each watershed were then estimated by the average storage changes of the two reanalysis data sets which provide all water balance components (i.e., ERA5-Land (Hersbach et al., 2020) and FLDAS (McNally et al., 2017) data sets). We then corrected annual scale watershed evaporation based on the mean estimated storage change from ERA5-Land and FLDAS.

Maximum annual evaporation for the 1983–2020 period for each watershed was used as a reference to assess the PE models. At the watershed scale, PE models are calculated based on Equations 1–4 using monthly meteorological variables retrieved from ERA5-Land (Hersbach et al., 2020) and FLDAS (McNally et al., 2017) data sets. Here, we assumed that soil heat flux is zero. PE values were calculated by the average of the two reanalysis data sets. In order to analyze relationship between PE and soil moisture, we obtained soil moisture derived by the Soil Moisture Active Passive (SMAP) satellite mission and calculated annual mean percent soil moisture for each watershed using the NASA-USDA Enhanced SMAP data set (Mladenova et al., 2020). Percent soil moisture for the NASA-USDA Enhanced SMAP data set represents plant available water divided by the total water holding capacity of the soil profile (1 m depth). The soil moisture, precipitation, and reanalysis data sets were downloaded from Google Earth Engine (Gorelick et al., 2017) using the rgee R package (Aybar et al., 2020), while the USGS runoff observations were downloaded using the dataRetrieval R package (De Cicco et al., 2018).

It should be noted that this reference point (i.e., maximum annual evaporation for the 1983–2020 period for each watershed) would be lower than the PE models if a watershed is relatively dry as one could expect some period of water stress even during the wettest year. Therefore, in addition to analyzing the full 328 watersheds using the maximum annual evaporation approach, we repeated the analysis for the 125 watersheds remaining after filtering out dry watersheds (watersheds in which the annual mean SMAP percent soil moisture is <75%). This criterion has been chosen as the PE rate can be sustained to 75% of the total water holding capacity in the soil moisture model used in the NASA-USDA Enhanced SMAP data set.

CMIP5 Simulations

For the global-scale analysis, we used 25 publicly available GCMs that participated in CMIP5. Although Coupled Model Intercomparison Project Phase 6 models recently became publicly available, we used CMIP5 models to enable comparison of this study with the relevant previous studies (Berg & McColl, 2021; Milly & Dunne, 2016; Yang et al., 2019). Latent and sensible heat fluxes, air temperature, RH, wind speed, barometric pressure, precipitation, runoff, and evaporation were obtained from the models' outputs. CMIP5 models that provide all required variables for the PE calculations were selected, and the models include: ACCESS1-0, ACCESS1-3, CNRM-CM5, GISS-E2-R-CC, HadGEM2-CC, HadGEM2-ES, bcc-csm1-1, bcc-csm1-1-m, CanESM2, CESM1-CAM5, CSIRO-Mk3-6-0, GFDL-CM3, GFDL-ESM2G, GFDL-ESM2M, GISS-E2-H, GISS-E2-H-CC, GISS-E2-R, inmcm4, IPSL-CM5A-LR, IPSL-CM5A-MR, IPSL-CM5B-LR, MIROC-ESM, MIROC-ESM-CHEM, MIROC5, and MRI-CGCM3. These CMIP5 models' output were obtained from the Columbia University Lamont-Doherty Ocean and Climate Physics Data Library (http://strega.ldeo.columbia.edu:81/CMIP5/).

We used monthly outputs of the historical reference period (1981–2000) and the high emission future period (2081–2100 under RCP 8.5) to calculate the mean values of each period for PE, precipitation, and runoff. In order to calculate the multimodel mean and median, the models were regridded to a 2° × 2° resolution (Berg & McColl, 2021). Runoff output is not available for some models (the first six models in the above model list); in these cases, we estimated runoff as the difference between precipitation and evaporation (Milly & Dunne, 2016).

PE models were calculated based on Equations 1–4 using monthly CMIP5 models' meteorological output, and then aggregated into 20 years average for the historical and future periods, respectively. Since soil heat flux is not available in CMIP5 model outputs, we calculated AE (i.e., net radiation minus soil heat flux) as the sum of latent and sensible heat fluxes following recent studies (Berg & McColl, 2021; Milly & Dunne, 2016). We calculated PE only for the land fraction which does not include Greenland and Antarctica (Berg & McColl, 2021; Milly & Dunne, 2016).

Budyko Model

To evaluate hydrological implications of the varying increasing rates of the different PE models estimated by CMIP5 models, we estimated runoff (Q) based on the Budyko water balance model, which is forced by PE and precipitation (P). The Budyko water balance model can be written following (Milly & Dunne, 2016) [Image Omitted. See PDF]

In this equation, all variables should be understood as 20-year mean values over the historical or future periods. Although there are several functions representing the Budyko model, we select this original equation following Milly and Dunne (2016). We also tested an equation used by Yang et al. (2019), and found similar results, implying choice of the Budyko equation may be a minor issue at global scale although prediction skill can be improved at a regional scale.

Results Site-Level Evaluation 1: Performance

The performances of PE models are evaluated using FLUXNET2015 in situ observations around globe (Pastorello et al., 2020), by isolating non-water-stressed conditions (Maes et al., 2019) (Methods). The observed values of non-water-stressed evaporation (E_unstr) are used as a reference of the upper limit to evaporation to assess the PE models (PM, PT, MD, and PE-LAC) since PE should become equivalent to actual evaporation under non-water-stressed conditions. We found that PE-LAC most accurately reproduces observed E_unstr in terms of root mean square error (RMSE), coefficient of determination (R²), and regression slope relative to current PE models (Figures 2a–2d). The MD PE model yielded similar site-scale performance with the PE-LAC, while the widely used PM PE showed the lowest performance.

Figure 2. Observed non-water-stressed evaporation predicted by potential evaporation (PE) models, and their bias. From (a) to (d), the y-axis is weekly observations of non-water-stressed evaporation (Eunstr), and the x-axis is PE calculated by Penman-Monteith, Priestley-Taylor, Milly Dunne, and PE-LAC, respectively. Shaded points represent all evaporation observations while colored triangles represent Eunstr. The regression lines are based only on the colored triangles. From (e) to (h), biases of each model (y-axis) are depicted as a function of air temperature (T) and relative humidity. Color in panels (e–h) represents data density.

Importantly, the PE-LAC does not show significant bias regarding temperature and RH, unlike biases present in all current PE models (Figures 2e–2h). For example, the PM and PT models overestimate observed E_unstr when the temperature is high and/or RH is low. Our findings for the PM and PT models are consistent with the PE overestimation bias first reported by Milly and Dunne (Milly & Dunne, 2017) for models used to predict PE under future climate conditions. In contrast, our PE model exhibits the smallest bias with respect to temperature and RH, making it more appropriate for PE calculations when evaluating future climate scenarios in frameworks that use PE as a parameter.

Site-Level Evaluation 2: Climate Sensitivity of PE

We then use multiple-linear regression to determine the sensitivity of PE models to temperature ($$\frac{\partial \mathrm{P}\mathrm{E}}{\partial T}$$) and to RH ($$\frac{\partial \mathrm{P}\mathrm{E}}{\partial \mathrm{R}\mathrm{H}}$$). Regression slopes for each PE model and E_unstr are considered as the sensitivity in Figures 3a and 3b. It should be noted that since we only use temperature and RH as independent variables, these sensitivities represent not only direct effects (e.g., saturation vapor pressure) but also indirect effects (e.g., net radiation).

Figure 3. Climate sensitivities of potential evaporation (PE) models and measured evaporation for non-water-stressed periods (Eunstr) (a, b) and evaporative fraction (EF) of PE models as a function of climate variables (c–e). From (a) to (b), the y-axis of each panel is the PE sensitivities to temperature (T) and to relative humidity, respectively. The error bar indicates 2 standard error, and the dotted lines represent the Eunstr error bar range. From (c) to (e), the y-axis of each panel is EF of Penman-Monteith (PM), Priestley-Taylor (PT), and PE-LAC, respectively. The dashed line in panels (c–e) represents Milly Dunne EF (fixed at 0.8). The short black lines on the right-side margins of panels (c–e) shows EF distributions for PM EF, PT EF, and Equation 1 EF, respectively.

The PT model yielded the largest overestimates of observed E_unstr sensitivity to increasing temperature. Of the four PE models, the PM model was found to be the most largely biased for RH $$\left(\frac{\partial \mathrm{P}\mathrm{E}}{\partial \mathrm{R}\mathrm{H}}\right)$$ compared to the observed sensitivity of E_unstr to RH. These results show how the PM and PT models overestimate evaporative demand in warmer and drier future climates. On the other hand, the PE-LAC exhibits good performance in reproducing observed sensitivity to temperature and RH. It should be noted that the negative sensitivity of the PE-LAC model to RH is because only temperature and RH were used as independent variables in this analysis. Higher RH commonly corresponds to lower AE, and this may cause negative sensitivity. Likewise, the negative sensitivities to RH for the PT and MD models could be explained this way.

To further understand the influence of temperature and RH on each PE model, we depict each PE model's EF in Figures 3c–3e, where EF represents the ratio of evaporation to AE. Here, EF for MD 's model is fixed at a constant 0.8 by definition (dashed line). Increases in temperature result in increasing EF for the remaining three PE models (PE-LAC, PM, and PT) due to the saturation vapor pressure effect, but EF computed with our proposed model increases only modestly compared to the PM and PT models. Also, the EF determined with the PE-LAC decreases as RH declines, since RH reflects land surface dryness. As a result, EF determined with the PE-LAC covers a much narrower range of values than the PM and PT models, with values centered on MD 's 0.8 fixed value of EF. In contrast, the PM model's EF rapidly increases for declining RH due to the constant surface resistance assumption, and it is probably the key reason for the relatively poor performance of the PM model.

Watershed-Level Evaluation 1: Performance

We then evaluated the performance of the four PE models (PM, PT, MD, and PE-LAC) using runoff observations from 338 US watersheds for the period of 1983–2020 (Duan et al., 2006). Annual evaporation was estimated as the difference between observed precipitation (PRISM data set (Daly et al., 1994)) and USGS runoff observations considering groundwater storage change (Materials and Methods). Maximum annual evaporation for the 1983–2020 period for each watershed was selected and used as a validation criterion for the PE models which were parameterized using two reanalysis data sets: ERA5-Land (Hersbach et al., 2020) and FLDAS (McNally et al., 2017). As depicted in Figure 4, the PE-LAC and the PT model most accurately reproduce observed maximum watershed evaporation, as an indicator of upper limit of evaporation, in terms of RMSE, R², and regression slope, while the PM model leads to the least accurate results. These watershed-scale results are consistent with the results we obtained for the site-level analysis.

Figure 4. Observed maximum annual evaporation at US watersheds predicted by potential evaporation (PE) models. From (a) to (d), the y-axis is annual water balance evaporation for each of 338 watersheds, and the x-axis is PE calculated by Penman-Monteith, Priestley-Taylor, Milly Dunne, and our proposed model, respectively. Shaded points represent annual water balance evaporation while colored triangles represent maximum annual evaporation for each watershed. The regression lines are based only on the colored triangles. The color bar indicates the annual mean percent soil moisture of each watershed derived from Soil Moisture Active Passive satellite observation.

Nevertheless, it is possible that the maximum annual evaporation for the 1983–2020 period for each watershed would be smaller than PE if a watershed is relatively dry. Indeed, all PE models overestimate the maximum annual evaporation for drier watersheds in Figure 4 (red colored points). Interestingly, the PE-LAC overestimates to a lesser extent, implying the PE-LAC can be a practical upper limit of evaporation even for drier watersheds. In order to validate the performance of the PE models further, we repeated the analysis using only wet watersheds by selecting watersheds in which the annual mean SMAP percent soil moisture is higher than 75%. This additional analysis using only wet watersheds still shows that the PE-LAC model (RMSE = 75.9 mm, R² = 0.55, slope = 0.98) and the PT model (RMSE = 74.3 mm, R² = 0.58, slope = 0.79) more accurately reproduce observed maximum watershed evaporation compared to the PM model (RMSE = 89.4 mm, R² = 0.45, slope = 0.65) and the MD model (RMSE = 107.6 mm, R² = 0.52, slope = 1.1).

Watershed-Level Evaluation 2: Evaporation Sensitivity to PE

Using the long-term watershed observations, we evaluate the sensitivity of annual scale evaporation to each PE model ($$\frac{\partial E}{\partial \mathrm{P}\mathrm{E}}$$). We use multiple-linear regression to determine the sensitivity of evaporation to PE and to precipitation (P) assuming evaporation is constrained by precipitation (supply side) as well as PE (demand-side). Regression slopes are considered as the sensitivity. Theoretically, dry watershed ($$\frac{P}{\mathrm{P}\mathrm{E}}< 1$$) evaporation is largely determined by precipitation while wet watershed ($$\frac{P}{\mathrm{P}\mathrm{E}} > 1$$) evaporation is primarily determined by PE. Therefore, $$\frac{\partial E}{\partial \mathrm{P}\mathrm{E}}$$ should be around zero to one for dry watersheds while $$\frac{\partial E}{\partial \mathrm{P}\mathrm{E}}$$ should be close to one for wet watersheds.

For wet watersheds, where the role of PE is important in controlling evaporation, $$\frac{\partial E}{\partial \mathrm{P}\mathrm{E}}$$ determined by MD model and the PE-LAC are aligned with this theoretical expectation while the PM model is not (Figure 5a). PM model's $$\frac{\partial E}{\partial \mathrm{P}\mathrm{E}}$$ is significantly less than unity for most wet watersheds, meaning changes in PE are always larger than changes in evaporation. This finding implies that future evaporation can be overestimated if evaporation is constrained by the PM PE model even for wet watersheds.

Figure 5. Sensitivities of watershed annual evaporation to potential evaporation (PE) models (a) and relationship between PE models and soil moisture (b–e). (a) Each watershed sensitivity is group by ratio between precipitation and the Penman-Monteith (PM) PE. From (b) to (e), the y-axis is PM, Priestley-Taylor, Milly Dunne, and our proposed model, respectively, and the x-axis is annual mean percent soil moisture retrieved by the NASA-USDA enhanced Soil Moisture Active Passive product. Each regression line represents one of the 328 watershed. The inset shows the slope of each regression line.

Is there any theoretical reason underlying this result? Fundamentally, terrestrial evaporation is constrained by soil moisture limitations (supply) and PE (demand), and thus one can write evaporation as $$E={f}_{\mathrm{s}\mathrm{m}}\mathrm{P}\mathrm{E}$$, where $${f}_{\mathrm{s}\mathrm{m}}$$ represents the soil moisture constraint ranging from zero to one. Therefore, the sensitivity of evaporation to PE can be written as follows. [Image Omitted. See PDF]

If $${f}_{\mathrm{s}\mathrm{m}}$$ and PE are independent, the first term becomes negligible, and thus one can write $$\frac{\partial E}{\partial \mathrm{P}\mathrm{E}}={f}_{\mathrm{s}\mathrm{m}}$$. For wet watersheds, $${f}_{\mathrm{s}\mathrm{m}}$$ is close to one and thus $$\frac{\partial E}{\partial \mathrm{P}\mathrm{E}}$$ becomes one in principle. On the other hand, if $${f}_{\mathrm{s}\mathrm{m}}$$ and PE have a negative correlation, $$\frac{\partial E}{\partial \mathrm{P}\mathrm{E}}$$ cannot approach one even for wet watersheds due to the first term in Equation 6. As depicted in Figures 5b–5e, the PM model shows the most apparent negative relationship with satellite-derived soil moisture (Mladenova et al., 2020), which explains why $$\frac{\partial E}{\partial \mathrm{P}\mathrm{E}}$$ does not generally approach one using the PM model, even for wet watersheds. In contrast, the MD model and the PE-LAC do not show an apparent dependency on soil moisture.

Projections of the PE Models to Anthropogenic Climate Change

Following Milly and Dunne (2016), we compared century-scale changes in PE models from a historical reference period (1981–2000) to a future scenario (2081–2100) using 25 CMIP5 models under Representative Concentration Pathway 8.5 (Methods and Table S1 in Supporting Information S1). The four tested PE models suggest increasing PE over most of the terrestrial regions for the future relative to the reference period, but the magnitude of the changes vary substantially. Consequently, projected median PE changes over the global land surface vary, from smallest to largest: MD (60), PE-LAC (83), PT (136), and PM (227 mm yr⁻¹) (Figure 6). These differences originate from the EF responses to future climatic conditions represented by the individual PE models. Both rising temperatures and declining RH result in increased PM EF, and hence the largest increase in PE is projected by the PM PE model. For the PT model, EF is not directly affected by RH, while rising temperatures increase EF. Thus, the projected mean increase in PT PE is lower than for the PM PE model. In contrast, declining RH projected for future climate conditions results in reduced EF with PE-LAC, and thus PE changes projected in the present study are lower than those obtained using PT and PM models.

Figure 6. Global scale changes in potential evaporation (PE) in the future period 2081–2100 (RCP 8.5) relative to the historic period 1981–2000. The bars represent the ensemble median of 25 Coupled Model Intercomparison Project Phase 5 models while the points and error bars indicate the ensemble mean ±1 standard deviation. Changes in different PE models and their components are presented and change in ocean evaporation is presented for reference. Here, EF is evaporative fraction and AE is available energy. The yellow areas represent changes in PE due to changes in AE while the orange areas represent changes in PE due to changes in EF.

Interestingly, we found that changes in ocean evaporation projected by GCMs are most closely matched by our PE-LAC model (Figure 6). Over the ocean, temperatures are increasing at a slower rate than for the land. Also, RH is roughly steady over the ocean, unlike the declining trend in RH over land, a difference known as the “land-ocean contrast” (Byrne & O’Gorman, 2016, 2018). “Land-ocean contrast” effects on evaporative demand can be reconciled using the PE-LAC because the influence of temperature on EF is conditioned by changes in RH. This may be a reason why our proposed PE model is well-matched with the ocean evaporation change. On the other hand, in the PM PE model, the “land-ocean contrast” increases PE compared to ocean evaporation, resulting in a projected change in PE for the land surface that is nearly 3 times that projected for evaporation from the ocean.

Varying Runoff Projections Resulting From the Different PE Models

In order to evaluate hydrological implications of the varying increasing rates of the different PE models, we compared runoff changes estimated using the Budyko water balance approach forced by PE and precipitation (Figure 7), following Milly and Dunne (2016) (Methods). The Budyko-estimated runoff change forced with our proposed PE model (17 mm yr⁻¹ at global scale) most closely matches the direct CMIP5 projections (20 mm yr⁻¹). The MD PE based Budyko runoff change (25 mm yr⁻¹) slightly overestimates the direct CMIP5 projections. In contrast, when the PM model is used, the Budyko estimated runoff change (−1 mm yr⁻¹) largely underestimates the direct CMIP5 output. The PT model (7 mm yr⁻¹) is a better comparator with CMIP5 than PM, but it still underestimates the direct CMIP5 projections.

Figure 7. Multi Coupled Model Intercomparison Project Phase 5 (CMIP5) models median of the relative change of the annual mean runoff (a–e), and global scale changes runoff (f) in the future period 2081–2100 (RCP 8.5) relative to the historic period 1981–2000. From (a) to (d), runoff change estimated by Budyko model forced by each potential evaporation (PE) model while (e) represents direct CMIP5 models' output. Here, relative change indicates future minus historical divided by average of historical and future. (f) Changes in runoff estimated by the Budyko model using different PE models and change in runoff directly projected by CMIP5 is presented for reference. The bars represent the ensemble median of 25 CMIP5 models while the points and error bars indicate the ensemble mean ±1 standard deviation.

In terms of spatial patterns, the Budyko-estimated runoff forced with the PM PE model shows apparent bias while other models show reasonable agreement with the direct CMIP5 projections. Especially, the PM PE model shows a large negative bias in wet regions such as the tropics. As we demonstrated in the watershed scale analysis (Figure 5 and Section 4.4), the PM PE model overestimates evaporation increase particularly for wet watersheds, where the role of PE is important in controlling evaporation. Therefore, evaporation increases could be largely overestimated in wet regions if one applies the PM PE model, which results in a large negative bias in runoff projections for the wet regions. This runoff projections bias can be reduced using the PT PE model and further reduced using the MD PE model and the PE-LAC.

Discussion Shortcomings of the Constant Surface Parameter in the PM PE Model

The overall results from our hierarchical analyses consistently indicated that the widely used PM PE model showed the worst performance while the proposed PE-LAC model showed the best performance in reproducing the non-water-stress evaporation and its changes in relation to climatic drivers. The lower performance of the PM PE model may be due to several reasons. First, one possible problem is the usage of the FAO PM equation that is designed for a short height reference crop although empirical parameters of this model (e.g., surface resistance) may vary by ecosystem. However, this may be not a major reason since the PM model showed poorer performance than the PT and MD models in a prior study despite using ecosystem-specific parameter tuning for the PM model (Maes et al., 2019). This means that the poor performance of the PM model may not be limited to the reference crop-specific parameters, but could be due to the sensitivity of surface resistance to atmospheric conditions (e.g., vapor pressure deficit and CO₂) even if one selects non-water-stressed evaporation observations. Surface resistance cannot be assumed to be constant in time at the plant stand scale, even in non-water-stressed conditions (e.g., McNaughton & Black, 1973). Therefore, if one estimates evaporation even for non-water-stressed conditions by using process-based models like the PM equation, surface resistance should be correctly modeled as commonly conducted in land surface models using semi-empirical stomatal models. As such, the constant surface resistance assumption in the PM PE model is likely to be problematic in reproducing non-water-stressed evaporation.

Because vapor pressure deficit and soil moisture are negatively correlated with each other, the PM PE model with constant surface resistance should have a negative correlation with soil moisture (Figure 5). This unique feature of the PM PE model can be useful if the PM PE model is used for revealing drier conditions (e.g., wildfire risk). However, the negative correlation between the PM PE and soil moisture even in wet watersheds is problematic as changes in PM PE is not equivalent to changes in actual evaporation anymore even for wet watersheds (Section 4.4 and Figure 5). This could be a main reason for the low performance of the PM PE model in reproducing changes in non-water-stress evaporation and runoff for future climatic conditions (Sections 4.5 and 4.6). Therefore, the constant surface resistance in the PM PE model should be carefully evaluated in an analysis using the PM PE model, particularly under warming and drying climatic conditions with elevating atmospheric CO₂ concentrations.

Theoretical Justification of the MD PE Model

Another important finding is that the proposed PE-LAC behaves very similarly to the MD PE model. The EF of the PE-LAC is close to 0.8 for a wide range of temperature and RH (Figure 3e). Therefore, the proposed PE-LAC model could be understood as a theoretical justification for the MD PE model. As a result, the PE-LAC and the MD PE model act similar to each other in other analyses. These results imply that the proposed PE-LAC is dominantly controlled by AE as effects of temperature and RH commonly cancel each other out.

Nevertheless, the PE-LAC model shows slightly better performance than the MD PE model, particularly in the watershed scale analysis (Figure 4). Maes et al. (2019) showed that the empirical parameter of the MD PE model (i.e., $${\alpha }_{\mathrm{M}\mathrm{D}}$$ in Equation 4) slightly varies with land cover types. They found that $${\alpha }_{\mathrm{M}\mathrm{D}}$$ is lower than 0.8 for semi-dry regions (e.g., savanna and shrubland) or cold regions (e.g., evergreen needleleaf forest). The PE-LAC model results are consistent with the empirical adjustment of the $${\alpha }_{\mathrm{M}\mathrm{D}}$$ parameter as the EF of the PE-LAC decreases when RH or temperature decreases (Figure 3e).

Potential Caveats

We should point out potential caveats in interpreting this analysis. One should be carefully when interpreting our results as the definition of PE for the PM model is quite different from the definition of the PT PE model and our PE-LAC model. For example, the PM type models with constant surface resistance (including the Penman equation) could be better than other PE models for determining open water evaporation (e.g., pan evaporation), wetland evaporation, or well-watered reference crop evaporation, although its performance was low for reproducing maximum evaporation at an ecosystem and watershed scales.

Furthermore, there are several methodological limitations. The FLUXNET2015 data set includes systematic uncertainty as eddy covariance observation is subject to the well-known energy balance closure problem (Mauder et al., 2020; Wilson et al., 2002). We provide additional analyses and discussion regarding the energy balance problem in Text S1 and Figures S2–S4 in Supporting Information S1. Also, the method of selecting non-water-stressed evaporation using EF introduced by Maes et al. (2019) includes methodological uncertainty as they explicitly discussed in their article. For instance, filtering out local advection conditions (i.e., negative sensible heat flux) could be unfavorable to the PM PE model as the aerodynamic term of the PM PE model could be important in the local advection conditions. As for the watershed scale analysis, the water balance evaporation is subject to the uncertainty of the storage change although we tried to resolve this problem using reanalysis data sets. Furthermore, we used changes in ocean evaporation as a reference point under future climatic conditions in the analysis using CMIP5 models, and this reference point implicitly assumes that non-water-stressed terrestrial evaporation acts similarly to ocean evaporation under the warming climate. In spite of these limitations of our analyses, the consistent results across spatiotemporal scales support the key findings of this study.

Conclusions

Although warmer and drier air is widely recognized to correspond to high atmospheric evaporative demand, the terrestrial supply side mechanism constraining atmospheric aridity tends to be overlooked. That is, soil moisture represents a supply constraint to warmer and drier air through the land-atmosphere feedback. This feature is ignored in current PE models which assumes a saturated surface for any given meteorological condition. As such, current PE models untenably overestimate the upper limit of evaporation for warmer and drier conditions and thus overestimate evaporation change, even for non-water limited conditions. The systematic biases of current PE models have serious implications that could lead to inappropriate planning in relation to needed climate change mitigation and adaptation. Arguably, to fundamentally resolve this problem, one should consider the overlying atmosphere as a coupled system with the land surface instead of solely as a source of evaporative demand that is independent from terrestrial processes.

Perhaps, current PE models would still be essential tools for some purposes (Vicente-Serrano et al., 2020). For instance, the PM PE model could be a useful indicator of wildfire risk in that high PM PE values represent drier air (Y. Huang et al., 2020; McEvoy et al., 2020). It should be noted that our analyses are not intended to deny or replace the various applications of these PE models. However, if the current PE model is used as the upper limit of evaporation that controls actual evaporation and evapotranspiration or other water balance components, the systematic biases toward drying are unavoidable under anthropogenic climate change. Therefore, care should be taken while applying and interpreting the PE models.

Internally consistent climate simulations that incorporate coupled land-ocean-atmosphere processes such as GCMs can be a solution to this issue. However, these sophisticated simulations and low-resolution outputs cannot fully replace widely used operational approaches based on PE such as watershed hydrological models, crop growth and crop water use models, drought and aridity analysis, and global satellite-based evaporation products (e.g., MOD16, GLEAM or PT-JPL). By presenting a land-atmosphere coupled PE model that can be easily implemented in established hydrologic approaches using readily measurable parameters, we believe that the land-atmosphere coupling perspective can be effectively implemented into a wide range of hydrological planning tools, particularly those focused on evaluating responses to changing climatic conditions.

A AppendixThe Derivative of Evaporation With Respect to Soil Moisture

Vertical water vapor flux from the Earth's surface (i.e., evaporation) is constrained by the difference in specific humidity (q) between the land surface and the atmosphere. If we express specific humidity as the product of saturation specific humidity and relative humidity (RH) following Kim et al. (2021), a flux gradient expression of evaporation can be written as follows. [Image Omitted. See PDF]where, E is evaporation, λ is the latent heat of vapourization, ρ is air density, q* is saturation specific humidity, T is temperature, RH is relative humidity, r_a is aerodynamic resistance to water vapor transfer (s m⁻¹). The subscript a indicates the atmospheric state near the land surface, and the subscript 0 indicates the land surface.

In this derivation, we assume well-coupled land-atmosphere conditions that are similar to the SFE state where the near-surface atmosphere is much more sensitive to changes in surface fluxes compared to changes in atmospheric boundary layer processes (McColl et al., 2019). Therefore, the derivative of E with respect to soil moisture (θ) can be expressed as follows by assuming the atmospheric state, as well as surface state, change with soil moisture because of this land-atmosphere coupling. [Image Omitted. See PDF]where, $$s\left(=\frac{\partial {q}^{\ast }}{\partial T}\right)$$ is the linearized slope of saturation specific humidity versus temperature (kg water vapor (kg moist air)⁻¹ K⁻¹). We assume identical s for the land surface and the atmosphere as is typically assumed in evaporation models. Also, we assume that r_a is independent to soil moisture.

The land surface state and the atmospheric state can be related as follows. [Image Omitted. See PDF] [Image Omitted. See PDF]

Substituting Equations A3 and A4 into Equation A2 yields [Image Omitted. See PDF]

If we approximate $${q}^{\ast }\left({T}_{0}\right)-{q}^{\ast }\left({T}_{a}\right)\approx s\left({T}_{0}-{T}_{a}\right)$$ and assume identical r_a for water vapor and sensible heat transfer, the second and the third terms of the right-hand side of Equation A5 can be expressed using sensible heat flux (i.e., $$H=\rho {c}_{p}\frac{{T}_{0}-{T}_{a}}{{r}_{a}}$$). [Image Omitted. See PDF]where, $$\gamma \left(=\frac{{c}_{p}}{\lambda }\right)$$ is the psychrometric constant and c_p is the specific heat capacity of the air. We then substitute the energy balance equation (i.e., $$H=\left({R}_{n}-G\right)-\lambda E$$) into the second and the third terms of the right-hand side of Equation A6, and then arrange it as follows. [Image Omitted. See PDF]

As Equation A7 is too complicated to analyze, it would be useful to introduce additional assumptions to simplify the equation. Incoming solar energy to the land surface can be considered as independent to soil moisture. Although other components of AE (i.e., $${R}_{n}-G$$) such as outgoing longwave radiation are affected by surface temperature, albedo, and emissivity that can be affected by soil moisture, their effects on AE are relatively smaller than incoming solar radiation. Therefore, we assume $$\frac{\partial \left({R}_{n}-G\right)}{\partial \theta }$$ is relatively negligible. This assumption is empirically supported in Figure 5d as the correlation between soil moisture and AE is weak. Substituting $$\frac{\partial \left({R}_{n}-G\right)}{\partial \theta }=0$$ into Equation A7 yields: [Image Omitted. See PDF]

Next, we replace T₀ with moist static enthalpy to further simplify Equation A8. Here, moist static enthalpy is known as $$d{h}_{0}={c}_{p}d{T}_{0}+\lambda d{q}_{0}$$ and the vertical flux of moist static enthalpy multiplying air density is equivalent to AE. Since we assume there are no changes in AE, moist static enthalpy at the land surface also does not change with soil moisture (i.e., $$\frac{\partial {h}_{0}}{\partial \theta }=0$$). This assumption is supported by climate simulation results in Berg et al. (2016). If we express specific humidity as a multiplication of q* and RH, moist static enthalpy can be written as $$d{h}_{0}=\left({c}_{p}+\lambda {\mathrm{R}\mathrm{H}}_{0}s\right)d{T}_{0}+\lambda {q}^{\ast }\left({T}_{0}\right)d{\mathrm{R}\mathrm{H}}_{0}$$. Thus, temperature change can be written as $$d{T}_{0}=\frac{1}{\lambda \left({\mathrm{R}\mathrm{H}}_{0}s+\gamma \right)}d{h}_{0}-\frac{{q}^{\ast }\left({T}_{0}\right)}{{\mathrm{R}\mathrm{H}}_{0}s+\gamma }d{\mathrm{R}\mathrm{H}}_{0}$$. Substituting this equation into the first term of the right-hand side of Equation A8 and assuming $$\frac{\partial {h}_{0}}{\partial \theta }=0$$ yield [Image Omitted. See PDF]

Kim et al. (2021) found that $$\frac{{\mathrm{R}\mathrm{H}}_{0}s+\gamma }{{\mathrm{R}\mathrm{H}}_{a}s+\gamma }\approx 1.1\frac{{q}^{\ast }\left({T}_{a}\right)}{{q}^{\ast }\left({T}_{0}\right)}$$ for most inland regions. This finding implies that $$\frac{{\mathrm{R}\mathrm{H}}_{0}s+\gamma }{{\mathrm{R}\mathrm{H}}_{a}s+\gamma }$$ should be between 1 and 1.1 for wet surfaces (i.e., RH₀ > RH_a) and for advection free (i.e., T₀ > T_a) conditions. Based on this reasoning, we approximate $$\frac{{\mathrm{R}\mathrm{H}}_{0}s+\gamma }{{\mathrm{R}\mathrm{H}}_{a}s+\gamma }\approx 1$$ in the first term of Equation A9, and thus Equation A9 becomes as follows: [Image Omitted. See PDF]

Next, we substitute the PM_rh actual evaporation model (Kim et al., 2021) to the first term of the right-hand side of Equation A10. Here, the PM_rh evaporation model provides an equation for RH flux (i.e., $$\frac{\rho {c}_{p}{q}^{\ast }\left({T}_{0}\right)}{{\mathrm{R}\mathrm{H}}_{a}s+\gamma }\frac{{\mathrm{R}\mathrm{H}}_{0}-{\mathrm{R}\mathrm{H}}_{a}}{{r}_{a}}=\lambda E-\frac{{\mathrm{R}\mathrm{H}}_{a}s\left({R}_{n}-G\right)}{{\mathrm{R}\mathrm{H}}_{a}s+\gamma }$$). Therefore, Equation A10 becomes [Image Omitted. See PDF]

Equation A11 can be solved for E by substituting $$\frac{\partial E}{\partial \theta }=0$$ and $$\frac{\partial {\mathrm{R}\mathrm{H}}_{0}}{\partial \theta }=\frac{\partial \mathrm{R}\mathrm{H}}{\partial \theta }$$ and then assuming $$\frac{\partial {\mathrm{R}\mathrm{H}}_{0}}{\partial \theta } > 0$$, which yields [Image Omitted. See PDF]

Equation A12 is exactly equivalent to Equation 1 in the main text. It should be noted that we simplify RH_a to RH in the main text.

B AppendixThe Solution of Equation A11

In Appendix A, we derived Equation A12 by introducing two constraints ($$\frac{\partial E}{\partial \theta }=0$$ and $$\frac{\partial {\mathrm{R}\mathrm{H}}_{0}}{\partial \theta }=\frac{\partial \mathrm{R}\mathrm{H}}{\partial \theta }$$) to Equation A11. However, in order to demonstrate that Equation A12 is the maximum instead of minimum, one may need to formally solve Equation A11. In this appendix, we suggest a mathematical solution to Equation A11.

Although it is difficult to derive the exact solution to Equation A11, it can be approximated with two additional simplifications. First, we assume $$\frac{2{\mathrm{R}\mathrm{H}}_{a}s+\gamma }{2{\mathrm{R}\mathrm{H}}_{a}s+2\gamma }$$ is constant, while recognizing that it actually varies slightly with soil moisture. This assumption could be justified using Figure 3e, as $$\frac{2{\mathrm{R}\mathrm{H}}_{a}s+\gamma }{2{\mathrm{R}\mathrm{H}}_{a}s+2\gamma }$$ is close to 0.8 for broad ranges of temperature and RH.

Second, we parameterize the last term of the right-hand side of Equation A11 as follows for the simplification of the equation. [Image Omitted. See PDF]where $$0< \alpha < 1$$. The justification of this parameterization is as follows. According to the PM_RH equation and the conceptual Figure 1, changes in the vertical RH gradient can partly explain changes in $$\lambda E$$. That is, the partial derivative of the PM_RH equation with respect to soil moisture has several terms including $$\frac{\rho {c}_{p}{q}^{\ast }\left({T}_{a}\right)}{\left({\mathrm{R}\mathrm{H}}_{a}s+\gamma \right){r}_{a}}\frac{\partial \left({\mathrm{R}\mathrm{H}}_{0}-{\mathrm{R}\mathrm{H}}_{a}\right)}{\partial \theta }$$. Therefore, this term partly explains changes in $$\lambda E$$. Furthermore, the evaporative fraction, which is largely affected by water availability instead of available energy, is dominantly controlled by the second term of the PM_RH equation (Kim et al., 2021), implying α is close to 1. This theoretical reasoning is empirically supported by Figure S1 in Supporting Information S1.

Substituting Equation B1 into Equation A11 yields [Image Omitted. See PDF] [Image Omitted. See PDF]

The integrating factor (F) of the above equation can be calculated as follows. [Image Omitted. See PDF]

Multiplying F to both sides of the above equation and rearranging it using $$\frac{\partial }{\partial \theta }\left[{\left({\mathrm{R}\mathrm{H}}_{0}s+\gamma \right)}^{\tfrac{2}{1-\alpha }}\right]=\frac{2s}{\left({\mathrm{R}\mathrm{H}}_{0}s+\gamma \right)(1-\alpha )}{\left({\mathrm{R}\mathrm{H}}_{0}s+\gamma \right)}^{\tfrac{2}{1-\alpha }}\frac{\partial {\mathrm{R}\mathrm{H}}_{0}}{\partial \theta }$$ yields [Image Omitted. See PDF] [Image Omitted. See PDF]

One can now solve Equation B4 by integrating with respect to soil moisture from is 0 to $$\theta $$, with assumption of $$\lambda E{\vert }_{0}=0$$: [Image Omitted. See PDF] [Image Omitted. See PDF]where $${\mathrm{R}\mathrm{H}}_{0}{\vert }_{\theta }$$ is the land surface RH when soil moisture is $$\theta $$ and $${\mathrm{R}\mathrm{H}}_{0}{\vert }_{0}$$ is the land surface RH when soil moisture is 0, which is close to zero. Thus, we can anticipate $$\frac{{\mathrm{R}\mathrm{H}}_{0}{\vert }_{0}s+\gamma }{{\mathrm{R}\mathrm{H}}_{0}{\vert }_{\theta }s+\gamma }< 1$$ for well-coupled land-atmosphere conditions. Since α is close to 1, $${\left(\frac{R{H}_{0\vert 0}s+\gamma }{R{H}_{0\vert \theta }s+\gamma }\right)}^{\tfrac{2}{1\mbox{-}\alpha }}$$ can quickly approach zero for wet conditions. Therefore, this solution clearly shows that Equation A12 is the upper limit of evaporation, that is, PE.

Acknowledgments

We express our thanks to the data providers, site investigators, and technicians without who this effort would not have been possible. We want to thank Brenda D’Acunha, Stephen Chignell, and Emily Mistick for helping with watershed data processing. We acknowledge the support of the Canadian Space Agency (CSA) Grant 21SUESIELH.

Conflict of Interest

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

Data Availability Statement

All data described in the main text are available. The FLUXNET2015 data set is available from FLUXNET (https://fluxnet.org/data/fluxnet2015-dataset/). The US watershed runoff data set is publicly available by USGS and can be downloaded using dataRetrieval R package (https://doi.org/10.5066/P9X4L3GE). The PRISM precipitation, SMAP soil moisture and the two reanalysis data sets (ERA5-Land and FLDAS) are publicly available and can be accessed using tools such as Google Earth Engine (https://doi.org/10.1002/2015WR017031). All climate model simulations are from CMIP5 and are publicly available, and are hosted on various servers including the Columbia University Lamont-Doherty Ocean and Climate Physics Data Library (http://strega.ldeo.columbia.edu:81/CMIP5/).