Emergent Simplicity of Continental Evapotranspiration

Evapotranspiration (ET) is challenging to model because it depends on heterogeneous land surface features—such as soil moisture, land cover type, and plant physiology—resulting in rising model complexity and substantial disagreement between models. We show that the evaporative fraction (ET as a proportion of available energy at the surface) can be estimated accurately across a broad range of conditions using a simple equation with no free parameters and no land surface information; only near‐surface air temperature and specific humidity observations are required. The equation performs well when compared to eddy covariance measurements at 76 inland continental sites, with prediction errors comparable to errors in the eddy covariance measurements themselves, despite substantial variability in surface conditions across sites. This reveals an emergent simplicity to continental ET that has not been previously recognized, in which land‐atmosphere coupling efficiently embeds land surface information in the near‐surface atmospheric state on daily to monthly time scales.


Introduction
Evapotranspiration (ET, the sum of evaporation and transpiration) is a key flux in the terrestrial water, energy, and carbon cycles. As a component of the water cycle, it is of fundamental importance to water resource management, agriculture, and ecosystem health. As a substantial flux of latent heat in the land surface energy budget, it effectively cools the near-surface atmosphere. Plants make trade-offs between transpiration and photosynthesis through stomatal regulation, which directly impacts the terrestrial carbon budget (Friedlingstein et al., 2013;Green et al., 2019). Since ET coregulates the water, energy, and carbon cycles, changes in any one cycle that impact ET can potentially propagate through the other cycles. However, ET has only been measured on a routine basis since the 1990s, and measurements are relatively sparse due to the high cost of installing and maintaining eddy covariance instruments. As a result, ET models remain essential for characterizing ET at regional and global scales.
The land surface strongly constrains ET (Monteith, 1965) but is challenging to model due to its heterogeneity. In an attempt to represent this heterogeneity, land surface models have become increasingly complex (Sellers et al., 1997), at the cost of introducing additional parameters, which are difficult to constrain at global scales. For example, one key parameter in modeling ET-the stomatal conductance g s -is often itself parameterized (Ball et al., 1987) as a function of relative humidity, atmospheric CO 2 concentration, and photosynthetic assimilation rate (A). This scheme includes two parameters and also requires a separate parameterization of A; one common scheme (Farquhar et al., 1980) for A requires at least another four parameters. In addition, the impact of water limitation on ET is often imposed by multiplying g s by an empirical function of soil moisture, varying between 0 and 1, requiring at least another two parameters (Trugman et al., 2018). Most of these parameters further vary between, and within, plant functional types. Constraining these and other parameters relevant to ET is a major challenge in land surface modeling. As a result, global models struggle to represent ET with high accuracy (Mueller & Seneviratne, 2014), particularly over inland continental regions (Ma et al., 2018).
A growing body of evidence suggests that modeling continental ET may be simpler than it first appears. Over many land surfaces, strong coupling between the land and atmosphere mean the land surface state becomes embedded in the atmospheric state at daily or longer time scales (Bouchet, 1963;Morton, 1969;McColl et al., 2019;Novick et al., 2016;Salvucci & Gentine, 2013;Zhou et al., 2019), implying that atmospheric observations alone should contain sufficient information to estimate ET. This is potentially useful because it implies that ET can be measured or modeled without explicitly measuring or modeling the land surface, Geophysical Research Letters 10.1029/2020GL087101 the source of many modeling errors. Various methods have exploited this relation and have been empirically successful in estimating ET while using relatively little land surface information (Brondani et al., 2019;Gentine et al., 2013Gentine et al., , 2016Rigden & Salvucci, 2015;Salvucci & Gentine, 2013. While theoretically possible, extracting the ET signal solely from atmospheric observations has proven more challenging. Methods for estimating ET over water-limited land surfaces almost always require at least some land surface information: At minimum, semiempirical surface roughness and stratified turbulence parameterizations are required, which typically require observations of vegetation height and wind speed (Gentine et al., 2016;Rigden & Salvucci, 2015;Salvucci & Gentine, 2013). The semiempirical parameterizations include parameters that are substantially uncertain and vary with plant functional type (Rigden et al., 2018), so their elimination is also desirable. Schemes for estimating ET from satellite observations also require at least some land surface inputs, such as visible, near-infrared, and/or thermal infrared radiances (Anderson et al., 2011;Bastiaanssen et al., 1998;Fisher et al., 2008;Mu et al., 2011) or observations of soil moisture and vegetation water content (Miralles et al., 2011).
A recent study (McColl et al., 2019) hypothesized that, in many inland continental regions, the near-surface atmosphere is in a state of "surface flux equilibrium" (SFE), in which the surface moistening and surface heating terms in the near-surface relative humidity budget approximately balance. On daily to monthly time scales, the theory predicts that the Bowen ratio, B = H E , can be estimated as where E and H are the latent and sensible heat fluxes (W m −2 ), respectively, q a is screen-level specific humidity (-), T a is screen-level air temperature (K), = 2.5008×10 6 (J kg −1 ) is the latent heat of vaporization of water, c p = 1, 005 (J kg −1 K −1 ) is the specific heat capacity of air at constant pressure, and R v = 461.5 (J kg −1 K −1 ) is the gas constant for water vapor. Here, all quantities that vary in time are daily (i.e., a 24-hr period) or monthly averages. The Bowen ratio can be converted to latent heat flux using the relation E = (1 + B) −1 (R n − G), where (1 + B) −1 is the "evaporative fraction" (EF), R n is net radiation at the surface (W m −2 ), and G is ground heat flux (W m −2 ). G is typically small compared to R n at daily or longer time scales and is often ignored.
Unlike previous approaches, this model of water-limited B has zero free parameters, does not require any land surface information as inputs (including vegetation height), and does not require wind speed. Two lines of evidence were given for the validity of this theory. First, a simple steady-state box model of the atmosphere, coupled to a land surface, reproduced SFE across a wide range of plausible conditions. Second, it was shown analytically that, under steady conditions, equation (1) is exactly equivalent to a recent method-the Evapotranspiration from Relative Humidity at Equilibrium (ETRHEQ) method-that has been extensively evaluated against observations and shown to be empirically successful (Rigden & Salvucci, 2015;Salvucci & Gentine, 2013).
The results of earlier work (McColl et al., 2019), which were primarily about gaining theoretical insight, were focused exclusively on steady conditions. In the real world, the diurnal cycle ensures that steady conditions are never achieved. Is equation (1) accurate in real-world, unsteady conditions? In response to this question, in this study, we provide empirical and analytical evidence that equation (1) also applies reasonably well in transient, real-world conditions. Specifically, we evaluate equation (1) against global eddy covariance measurements of ET, showing that it performs well over a wide range of conditions. This is significant because it shows that water-limited B over inland continental regions can be estimated solely using atmospheric observations at daily to monthly time scales. Since equation (1) is simple and contains no free parameters, it demonstrates an emergent simplicity to continental ET that has not been previously recognized.

Surface Flux Equilibrium (SFE)
We briefly review the concept of "surface flux equilibrium" proposed in McColl et al. (2019). The near-surface atmospheric state (specifically, near-surface air temperature and specific humidity) is sensitive, to some degree, to turbulent surface fluxes (sensible heat flux and evapotranspiration, respectively). If this sensitivity is sufficiently strong, then higher atmospheric specific humidity is a signal of higher evapotranspiration in the recent past; similarly, higher temperature is a signal of higher sensible heat flux in the recent past. One consequence of strong sensitivity is that, to a good approximation, a balance exists between the surface moistening and surface heating terms in the near-surface relative humidity budget; this balance defines SFE ( Figure 1 and supporting information Text S1). Furthermore, it leads directly to equation (1) and is exactly equivalent to a previous method (Rigden & Salvucci, 2015;Salvucci & Gentine, 2013) for estimating ET under steady conditions (Text S2). In this study, we test the hypothesis that SFE is a reasonable representation of real-world (i.e., transient) inland continental regions. Regions with substantial moisture or heat convergence (such as land regions near coasts) are not expected to be in SFE.

Data
Collocated measurements of half-hourly latent heat flux E (W m −2 ), sensible heat flux H (W m −2 ), net radiation R n (W m −2 ), ground heat flux G (W m −2 ), near-surface air temperature T a (K), relative humidity r (-), wind speed u (m s −1 ) and atmospheric pressure P (Pa) were obtained from the FLUXNET (fluxnet.ornl.gov) and AmeriFlux (ameriflux.lbl.gov) networks. For FLUXNET measurements, data with poor quality gap filling, as labeled in the quality control flags, were excluded. In addition, half-hourly measurements were excluded if the half-hourly surface energy imbalance was greater than 300 W m −2 , consistent with Rigden and Salvucci (2015).
At each site, measurements were aggregated to daily and monthly averages. Consistent with previous studies (Michel et al., 2016;Rigden & Salvucci, 2015;Salvucci & Gentine, 2013), we retained nighttime data in estimating daily and monthly averages (throughout this study, a "daily" average refers to an average over a full 24-hr period, rather than just daylight hours). In estimating daily averages, days with any missing half-hourly data were excluded. Furthermore, days in which the daily surface energy imbalance was greater than 50 W m −2 were also excluded, consistent with Rigden and Salvucci (2015). In estimating monthly averages, months with fewer than 15 complete days of data were excluded.

Evaluation at All Sites
In addition to the direct measurements of E provided by eddy covariance, the daily and monthly averaged observations at the 76 continental sites were used as inputs to calculate E in four different ways: (i) the residual of the surface energy budget, (ii) the Priestley-Taylor equation (Priestley & Taylor, 1972), (iii) an Advection-Aridity equation based on the complementary relationship (Brutsaert & Stricker, 1979), and (iv) equation (1).
First, latent heat flux was calculated by solving for E as the residual of the surface energy budget: The difference between this estimate of latent heat flux and the direct measurement from eddy covariance data is a measure of the surface energy budget imbalance and the inherent uncertainty in the eddy covariance measurements.

Second, the Priestley-Taylor equation
, was used similarly to calculate E. The Priestley-Taylor equation provides an estimate of daytime ET over saturated surfaces.
Third, an Advection-Aridity equation from Brutsaert and Stricker (1979) (their equation (20)) was used to calculate latent heat flux. This equation requires wind speed as an additional input. The formulation used is based on the original wind function proposed by Penman (1948), as given in equation (18) of Brutsaert and Stricker (1979). We use the wind speed reported at each FLUXNET and AmeriFlux site, even though the (unreported) measurement height may deviate from the 2 m reference height for which the equation was originally designed.
Fourth, equation (1) was used to estimate B and converted to latent heat flux using the relation E = R n −G

1+B .
We compare SFE to the Priestley-Taylor and Advection-Aridity equations because they are similar in not requiring calibration and in not requiring land surface information as inputs. The Priestley-Taylor equation does not require specific humidity as an input, meaning that it is simpler than SFE. The Advection-Aridity equation requires the same inputs as SFE, and also near-surface wind speed, meaning that it is more complicated than SFE. Various other forms of the complementary relationship exist, but many of these require land surface inputs (e.g., soil moisture or vegetation height). We further discuss differences between SFE and the complementary relationship in section 4.1.
The four different latent heat flux estimates were compared to the latent heat flux directly measured from eddy covariance at each site. The root-mean-square error (RMSE) and mean bias were calculated for all four estimates with respect to the directly measured latent heat flux. We choose the RMSE as a primary comparison metric of interest in this study because it is more appropriate for model selection: Since neither equation (1) nor the Priestley-Taylor equation nor the Advection-Aridity equation require any parameters to be fit to data, the Akaike information criterion (AIC) is a monotonic increasing function of RMSE. Therefore, comparing the models with respect to RMSE is equivalent to comparing them with respect to the AIC, a common approach to choosing between competing models.

Climatological Daily Mean Time Series at Focus Sites
Multiyear daily averaged data were further averaged over multiple years to create climatological daily mean time series of eddy covariance E (both directly measured and calculated as the residual of the surface energy budget), and E estimated using equation (1), the Priestley-Taylor equation, and the Advection-Aridity equation. The resulting time series were then smoothed using a 5-day moving average. We discarded sites in which there were missing values in any of the smoothed climatological daily mean E time series. Twenty-eight sites met this criterion. Six of these sites, denoted "focus sites," were chosen as broadly representative and are presented in  The calculated errors in the eddy covariance measurements and in the SFE estimates at each of these 28 sites were compared with the annual mean evaporative fraction. The evaporative fraction is a ratio of noisy measurements, meaning it is particularly sensitive to measurement errors. To mitigate this problem, the annual mean evaporative fraction was calculated as the annual average of the ratio of monthly mean latent heat flux to monthly mean available energy. (1) is relatively accurate across a wide range of conditions, at both daily and monthly time scales. For comparison, we calculate E as the residual of the surface energy budget components obtained from eddy covariance measurements and compare this estimate to the latent heat flux directly measured by eddy covariance (Figures 2a and 2e). If there were no errors in the eddy covariance measurements, there would be no difference between these two values; however, since there are well-known energy balance closure errors in eddy covariance measurements (Aubinet et al., 2012), the values differ. The errors in this comparison (Figures 2a and 2e) provide an approximate upper bound on the performance of any ET model when compared with eddy covariance measurements: Since the eddy covariance data are subject to errors of 10-30% themselves (Aubinet et al., 2012), even a model with perfect performance would exhibit errors of similar magnitude to those shown in Figures 2a and 2e when compared to eddy covariance measurements. In this light, the error statistics for equation (1) presented in Figures 2d and 2h are particularly impressive. The RMSE values also compare well to equivalent values from a recent intercomparison of substantially more complicated ET estimation algorithms, which reported values in the range 21-56 W m −2 based on daily eddy covariance measurements (Michel et al., 2016).

Figures 2d and 2h show that equation
On the other hand, if errors in eddy covariance measurements are large enough, then even an extremely inaccurate ET model might yield comparable error statistics. To mitigate this concern, we also present equivalent plots, except replacing equation (1) with two alternative estimates of ET of comparable complexity: the Priestley-Taylor equation (Figures 2b and 2f) and the Advection-Aridity equation (Figures 2c and 2g). The Priestley-Taylor equation, which is used to estimate daytime ET over saturated surfaces, has substantially poorer error statistics compared to equation (1) and the eddy covariance measurements, as expected, since most land surfaces are not saturated. The Advection-Aridity equation, which is based on the complementary relationship (Bouchet, 1963;Brutsaert & Stricker, 1979;Morton, 1969), also has poorer error statistics than  Tables S1 and S2. "Obs" is the latent heat flux directly measured from the eddy covariance data. "Obs EB" is the latent heat flux indirectly estimated (as the residual of the surface energy budget) from the eddy covariance data. "Eq 1" is the latent heat flux estimated using equation (1). "PT" is the latent heat flux estimated using the Priestley-Taylor equation (Priestley & Taylor, 1972). "AA" is the latent heat flux estimated using the Advection-Aridity equation (Brutsaert & Stricker, 1979). both equation (1) and the eddy covariance measurements. This is surprising, since the Advection-Aridity equation requires wind speed as an additional input, beyond those required by SFE. The poorer performance may be due, in part, to the fact that the wind speed observations used in this study are not necessarily measured at a reference height of 2 m. If so, this demonstrates the sensitivity of the Advection-Aridity equation to a quantity that is often not reported and to the details of uncertain parameterizations of turbulence. These results demonstrate that the comparable error statistics between SFE and eddy covariance measurements are more than just a coincidence. Overall, our results demonstrate that equation (1) performs considerably better than comparable methods and even exhibits errors similar to those in eddy covariance measurements. Figure 3 shows that equation (1) is also typically quite accurate at individual sites in inland continental regions. We highlight this performance at six focus sites. The focus sites are chosen to be hydrologically and climatologically diverse. They are also chosen to be broadly representative of the accuracy of equation (1), although there is considerable variability in accuracy between sites; equivalent plots for the full set of sites are presented in Figures S2-S5. Equation (1) broadly reproduces the observed climatology at each site, although it slightly overestimates ET under particularly dry conditions. Figure S8 shows that there is some weak structure to the errors in the SFE estimates: SFE errors are greatest at the driest sites, where annual mean evaporative fraction is particularly low.
Equation (1) outperforms the Priestley-Taylor and Advection-Aridity equations (in terms of errors estimated by comparing each method with E directly measured using eddy covariance) at the vast majority of sites ( Figure S6). Overall, compared to both the Priestley-Taylor and Advection-Aridity equations, it has both lower RMSE and mean bias at 93% of sites. Remarkably, it also has lower RMSE and mean bias compared with the estimated error in the eddy covariance measurements at a majority of sites (54% and 57% for RMSE and mean bias, respectively), even after substantially filtering the eddy covariance measurements to remove sites where energy balance closure errors are particularly high.

Relation to Previous Work
It has been shown previously that, under steady conditions, equation (1) predicted by SFE is exactly equivalent to ETRHEQ (McColl et al., 2019). In Text S2, we extend this comparison to provide analytical arguments for why SFE and ETRHEQ remain similar under unsteady diurnally varying conditions. The derivation shows that under unsteady conditions, the right-hand side of equation (1) is multiplied by a random variable, which is constrained by air temperature and specific humidity. We characterize this random variable (denoted ) at all eddy covariance sites with sufficient observations, showing that its mean is very close to one-both at all sites individually and across all sites-with relatively little variability around this value ( Figure S7). Therefore, SFE (equation (1)) is a reasonable approximation of ETRHEQ even under real-world, diurnally varying conditions (note that is not a calibration parameter, since this analysis shows that setting = 1 everywhere is reasonable). Since ETRHEQ has been extensively evaluated at a variety of temporal and spatial scales and shown to perform well (Rigden & Salvucci, 2015;Salvucci & Gentine, 2013), this lends further credibility to SFE, beyond our direct comparison with eddy covariance measurements in this study. The major advantages of SFE over ETRHEQ are that it does not require any land surface information (ETRHEQ requires vegetation height) or wind speed observations as inputs and avoids using uncertain semiempirical parameterizations of surface roughness and stratified turbulence. Furthermore, SFE consists of a single explicit equation, while ETRHEQ couples the energy balance and diffusion equations and needs to be solved using iterative methods.
Equation (1) is similar to, but distinct from, previous work on "equilibrium" ET, which predicts B ≈ R v c p T 2 a 2 q * (T a ) over idealized saturated surfaces (Raupach, 2001;Slatyer & McIlroy, 1961). In contrast to equilibrium estimates, SFE is accurate across a wide range of surface conditions, including substantially water-limited conditions (see McColl et al., 2019, for a detailed comparison between SFE and equilibrium ET).
SFE is also conceptually similar to the "complementary relationship" (Bouchet, 1963;Brutsaert & Stricker, 1979;Morton, 1969), which has been used to estimate ET from mostly meteorological data. In contrast to SFE, some implementations of the complementary relationship require a free parameter, which varies considerably between studies (Kahler & Brutsaert, 2006;Pettijohn & Salvucci, 2009). Mechanistic models of this parameter still require at least some land surface information, such as soil moisture (Aminzadeh et al., 2016). Other implementations of the complementary relationship, such as the Advection-Aridity equation (Brutsaert & Stricker, 1979) used in this study, do not require a free parameter. However, these implementations, when applied to daily time scales, still require observations of wind speed and, in many cases, vegetation height to parameterize turbulent transport. None of this is required by SFE.
Equation (1) predicts that, all else being equal, higher ET will be associated with higher relative humidity on daily to monthly time scales. Since drier atmospheres are known to promote more ET, this prediction might appear puzzling. In fact, both predictions are correct but apply to different time scales. Instantaneously, ET increases with increasing atmospheric dryness. However, over daily or longer time periods, ET moistens and cools the lower atmosphere, resulting in higher ET being associated with higher relative humidity (Betts, 2000;Novick et al., 2016;Zhou et al., 2019). Similarly, recent work has suggested that near-surface atmospheric humidity is largely set by the ocean and atmosphere, with the land surface playing a relatively passive role (Byrne & O'Gorman, 2016, a result that might appear to contradict SFE. However, that work applies to annual and longer time scales and global spatial scales, over which B will be substantially determined by atmospheric conditions. For example, while a wet surface with large ET will lead to a moist atmosphere at daily to monthly time scales, if anomalously high atmospheric drying persists long enough, the surface will dry out and eventually ET and near-surface humidity will reflect this. In this case, SFE still applies, but B is substantially determined by longer time scale processes, resulting in correlations with slower processes in the atmosphere and ocean. Furthermore, a recent study (Yin et al., 2019) providing a means for estimating the evaporative fraction (and, therefore, B) solely from atmospheric observations applies to multiyear time scales, rather than the daily to monthly time scales applicable to SFE.

Limitations
Equation (1) is substantially simpler than most methods for estimating ET, including the ETRHEQ method, but has comparable performance. This is particularly noteworthy, since it has zero parameters and exhibits good performance without tuning or calibration. Like ETRHEQ, SFE systematically overestimates ET when ET is very low and underestimates it when it is very high (Salvucci & Gentine, 2013). The positive bias at low values may be due to the longer time required to achieve equilibrium when ET is small, reducing the likelihood that the system achieves equilibrium in the real world (McColl et al., 2019;Raupach, 2001). The negative bias at high values may be due to enhanced sensitivity to boundary layer drying by entrainment, advection, or cloud base mass flux, particularly under strongly convective conditions or at irrigated sites in otherwise dry regions. Strong sensitivity to these fluxes would violate the assumptions of SFE. Entrainment drying has been shown mechanistically to enhance ET over saturated surfaces on time scales of several hours as the boundary layer grows during the day (De Bruin, 1983;Raupach, 2000). However, its contribution to the full diel (i.e., 24 hr) cycle average, relevant to SFE, is likely smaller, given that entrainment is much smaller at night. In any case, most continental regions fall somewhere between the dry and saturated limiting cases where these biases occur.
Furthermore, like ETRHEQ, SFE is not expected to hold in coastal areas, where substantial atmospheric convergence dominates surface fluxes in setting the atmospheric state. It is also not clear that SFE will hold over oceans, since relative contributions to the near-surface diel-averaged state from radiation, surface fluxes, entrainment, and cloud base mass flux are fundamentally different compared with land surfaces; nor is it clear that SFE will hold in future climates as the world warms. Despite these caveats, our results demonstrate that SFE is sufficient to explain most of the observed signal in global eddy covariance measurements over inland continental regions, spanning a wide range of conditions, without calibration or tuning of parameters.

Future Applications
By drastically simplifying the representation of continental ET, equation (1) opens up considerable new opportunities for constraining terrestrial water, energy, and carbon budgets. First, since observations of air temperature and specific humidity have much greater spatial and temporal coverage compared to existing observations of B, SFE substantially increases available observations of continental evapotranspiration. Second, equation (1) provides a constraint on models that might prove useful in evaluating existing models and developing new models that faithfully represent land-atmosphere coupling. Third, it opens up a new opportunity for remote sensing of ET, using remotely sensed estimates of near-surface atmospheric temperature and specific humidity, which is substantially less model-based compared to existing methods. to Dan Chavas for his publicly available dist_from_coast Matlab code; and to the contributors to the FLUXNET (fluxnet.ornl.gov) and AmeriFlux (ameriflux.lbl.gov) databases. K. A. M. acknowledges funding from a Winokur Seed Grant in Environmental Sciences from Harvard University's Center for the Environment. A. J. R. acknowledges funding from The Rockefeller Foundation Planetary Health Fellows program at Harvard University. This work used eddy covariance data acquired and shared by the FLUXNET community, including these networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux-TERN, TCOS-Siberia, and USCCC. The ERA-Interim reanalysis data are provided by 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. Funding for AmeriFlux data resources was provided by the U.S. Department of Energy's Office of Science.