Journal cover Journal topic
Hydrology and Earth System Sciences An interactive open-access journal of the European Geosciences Union
Journal topic
Hydrol. Earth Syst. Sci., 23, 925-948, 2019
https://doi.org/10.5194/hess-23-925-2019
Hydrol. Earth Syst. Sci., 23, 925-948, 2019
https://doi.org/10.5194/hess-23-925-2019

Research article 18 Feb 2019

Research article | 18 Feb 2019

# Potential evaporation at eddy-covariance sites across the globe

Potential evaporation at eddy-covariance sites across the globe
Wouter H. Maes1, Pierre Gentine2, Niko E. C. Verhoest1, and Diego G. Miralles1 Wouter H. Maes et al.
• 1Laboratory of Hydrology and Water Management, Ghent University, Coupure Links 653, 9000 Gent, Belgium
• 2Department of Earth and Environmental Engineering, Columbia University, New York City, New York 10027, USA
Abstract

Potential evaporation (Ep) is a crucial variable for hydrological forecasting and drought monitoring. However, multiple interpretations of Ep exist, which reflect a diverse range of methods to calculate it. A comparison of the performance of these methods against field observations in different global ecosystems is urgently needed. In this study, potential evaporation was defined as the rate of terrestrial evaporation (or evapotranspiration) that the actual ecosystem would attain if it were to evaporate at maximal rate for the given atmospheric conditions. We use eddy-covariance measurements from the FLUXNET2015 database, covering 11 different biomes, to parameterise and inter-compare the most widely used Ep methods and to uncover their relative performance. For each of the 107 sites, we isolate days for which ecosystems can be considered unstressed, based on both an energy balance and a soil water content approach. Evaporation measurements during these days are used as reference to calibrate and validate the different methods to estimate Ep. Our results indicate that a simple radiation-driven method, calibrated per biome, consistently performs best against in situ measurements (mean correlation of 0.93; unbiased RMSE of 0.56 mm day−1; and bias of −0.02 mm day−1). A Priestley and Taylor method, calibrated per biome, performed just slightly worse, yet substantially and consistently better than more complex Penman-based, Penman–Monteith-based or temperature-driven approaches. We show that the poor performance of Penman–Monteith-based approaches largely relates to the fact that the unstressed stomatal conductance cannot be assumed to be constant in time at the ecosystem scale. On the contrary, the biome-specific parameters required by simpler radiation-driven methods are relatively constant in time and per biome type. This makes these methods a robust way to estimate Ep and a suitable tool to investigate the impact of water use and demand, drought severity and biome productivity.

1 Introduction

Since its introduction 70 years ago by Thornthwaite (1948), the concept of potential evaporation (Ep), defined as the amount of water which would evaporate from a surface unconstrained by water availability, has been widely used in multiple fields. It has been incorporated in hydrological models dedicated to estimate runoff (e.g. Schellekens et al., 2017) or actual evaporation (Wang and Dickinson, 2012), as well as in drought severity indices (Sheffield et al., 2012; Vicente-Serrano et al., 2013). Long-term changes in Ep have been regarded as a driver of ecosystem distribution and aridity (Scheff and Frierson, 2013) and used to diagnose the influence of climate change on ecosystems based on climate model projections (e.g. Milly and Dunne, 2016). However, many different definitions of Ep exist, and consequently many different methods are available to calculate it. In recent years, there has been an increasing awareness of the impact of the underlying assumptions and caveats in traditional Ep formulations (Weiß and Menzel, 2008; Kingston et al., 2009; Sheffield et al., 2012; Seiller and Anctil, 2016; Bai et al., 2016; Milly and Dunne, 2016; Guo et al., 2017). As such, a global appraisal of the most appropriate method for assessing Ep is urgently needed. Yet, current formulations reflect a disagreement on the mere meaning of this variable, which requires the definition of some form of reference system (Lhomme, 1997). Ep has been typically defined as the evaporation which would occur in given meteorological conditions if water was not limited, either (i) over open water (Shuttleworth, 1993); (ii) over a reference crop, usually a wet (Penman, 1963) or irrigated (Allen et al., 1998) short green grass completely shading the ground; or (iii) over the actual ecosystem transpiring under unstressed conditions (Brutsaert, 1982; Granger, 1989).

A second source of disagreement on the definition of Ep relates to the spatial extent of the reference system and the consideration (or not) of feedbacks from the reference system on the atmospheric conditions. Several authors found it convenient to define Ep taking an extensive area as a reference system, because this reduces the influence of advection and entrainment flows (Penman, 1963; Priestley and Taylor, 1972; Brutsaert, 1982; Shuttleworth, 1993). Such an idealised extensive and well-watered ecosystem evaporating at maximal rate for the given atmospheric conditions can be expected to raise air humidity until the vapour pressure deficit (VPD) tends to zero. In this case, evaporation is only driven by radiative forcing and no longer by aerodynamic forcing. Meanwhile, others have defended the use of reference systems that are infinitesimally small (Morton, 1983; Pettijohn and Salvucci, 2009; Gentine et al., 2011b), in order to avoid the feedback of the reference system on aerodynamic forcing. The effect of the choice of reference system is best exemplified by the complementary relationship framework (Bouchet, 1964), which uses both approaches to link potential and actual evaporation, through $\left(\mathrm{1}+b\right){E}_{{\mathrm{p}}_{\mathrm{0}}}={E}_{\mathrm{pa}}+b{E}_{\mathrm{a}}$, with b an empirical constant (Kahler and Brutsaert, 2006; Aminzadeh et al., 2016), ${E}_{{\mathrm{p}}_{\mathrm{0}}}$ the evaporation from an extensive well-watered surface (i.e. in which the feedback from the ecosystem on the VPD and aerodynamic forcing is considered and where evaporation is only driven by a radiative forcing), Epa the evaporation from a well-watered but infinitesimally small surface (i.e. where evaporation is driven by both radiative and aerodynamic forcing) and Ea the actual evaporation (Morton, 1983).

In light of all this controversy, the net radiation of the reference system remains another point of discussion: some scientists argue that the (well-watered) reference system should have the same net radiation as the actual (water-limited) system (e.g. Granger, 1989; Rind et al., 1990; Crago and Crowley, 2005). Yet, this is inherently inconsistent as the surface temperature reflects the surface energy partitioning; thus a well-watered system transpiring at a potential rate is expected to have a lower surface temperature (Maes and Steppe, 2012) and correspondingly a higher net radiation (e.g. Lhomme, 1997; Lhomme and Guilioni, 2006). Meanwhile, to some extent, the albedo also depends on soil moisture (Eltahir, 1998; Roerink et al., 2000; Teuling and Seneviratne, 2008) and it can be argued that it should be adjusted to reflect well-watered conditions (Shuttleworth, 1993). Finally, extensive reference surfaces can be expected to exert a feedback, not only on the aerodynamic forcing, but also on the incoming radiation (via impacts on air temperature, humidity and cloud formation). Yet, these larger-scale feedbacks are not acknowledged when computing Ep, even when considering extensive reference systems.

As it can be concluded from the above discussion, a unique and universally accepted definition of Ep does not exist, and the most appropriate definition remains tied to the specific interest and application. Nonetheless, as different applications make use of different Ep formulations, a good understanding of the implications of the choice for a specific method is required (Fisher et al., 2011). For terrestrial ecosystems, the use of an open-water reference system is uninformative about the actual available energy and the aerodynamic properties of the actual ecosystem (Shuttleworth, 1993; Lhomme, 1997). The approach of considering an idealised well-watered crop system only takes climate forcing conditions into account and not the actual land cover; as such, it has become the standard to estimate aridity and trends in global drought (Dai, 2011). When the actual ecosystem transpiring at unstressed rates is considered as a reference system, both climate forcing conditions and ecosystem properties are taken into account. This has been the preferred approach when calculating Ep as an intermediate step to estimate actual evaporation, often by applying a multiplicative stress factor (S) varying between 0 and 1, such that Ea=SEp (e.g. Barton, 1979; Mu et al., 2007; Fisher et al., 2008; Miralles et al., 2011; Martens et al., 2017). This S factor can be considered analogous to the β factor used in some land surface models to incorporate the effect of soil moisture in the estimation of gross primary production and surface turbulent fluxes (Powell et al., 2013).

Several studies have attempted to compare and evaluate different Ep methods. Some of these studies have compared the performance of different Ep formulations in hydrological (Xu and Singh, 2002; Oudin et al., 2005a; Kay and Davies, 2008; Seiller and Anctil, 2016) or climate models (Weiß and Menzel, 2008; Lofgren et al., 2011; Milly and Dunne, 2016). Others considered the Penman–Monteith method as the benchmark to test less input-demanding formulations (e.g. Chen et al., 2005; Sentelhas et al., 2010). All these studies have their own merits, yet an evaluation of Ep methods based on empirical data of actual evaporation measurements is to be preferred (Lhomme, 1997). To date, such approaches have been hampered by limited data availability (Weiß and Menzel, 2008). Lysimeters provide arguably the most precise evaporation measurements available (e.g. Abtew, 1996; Pereira and Pruitt, 2004; Yoder et al., 2005; Katerji and Rana, 2011), but are sparsely distributed and not always representative of larger ecosystems. Pan evaporation measurements are more easy to perform and are broadly available (Zhou et al., 2006; Donohue et al., 2010; McVicar et al., 2012) but provide a proxy of open-water evaporation, rather than actual ecosystem potential evaporation; they also exhibit biases related to the location, shape and composition of the instrument (Pettijohn and Salvucci, 2009). Eddy-covariance measurements are an attractive alternative, but, apart from an unpublished study by Palmer et al. (2012), have so far only been used in Ep studies focusing on local to regional scales (Jacobs et al., 2004; Sumner and Jacobs, 2005; Douglas et al., 2009; Li et al., 2016).

The overall purpose of the present work is to identify the most suitable method to estimate Ep at the ecosystem-scale across the globe. Because we are using an empirical dataset of actual evaporation at FLUXNET sites, the reference system considered in this study is the actual ecosystem, so Ep is defined as the evaporation of the actual ecosystem when it is completely unstressed. As mentioned above, this definition is the most suitable for hydrological studies, studies of ecosystem drought and derivations of actual evaporation through constraining Ep calculations. Following this definition, Ep is similar to ${E}_{{\mathrm{p}}_{\mathrm{0}}}$ in the complementary relationship. We used the most recent and complete eddy-covariance database available, i.e. the FLUXNET2015 archive (http://fluxnet.fluxdata.org/, last access: 14 February 2019). The most frequently adopted Ep methods are applied based on standard parameterisations as well as calibrated parameters by biome and are inter-compared in order to gain insights into the most adequate means to estimate Ep from ecosystem to global scales.

2 Material and methods

## 2.1 Selection of Ep methods

Methods to calculate Ep can be categorised based on the amount and type of input data required. In this overview, we will only discuss the ones evaluated in our study, which are arguably the most frequently used.

### 2.1.1 Methods based on radiation, temperature, wind speed and vapour pressure

The well-known Penman–Monteith equation (Monteith, 1965) expresses latent heat flux λEa (W m−2) as

$\begin{array}{}\text{(1)}& \mathit{\lambda }{E}_{\mathrm{a}}=\frac{s\left({R}_{\mathrm{n}}-G\right)+\frac{{\mathit{\rho }}_{\mathrm{a}}{c}_{p}\mathrm{VPD}}{{r}_{\mathrm{aH}}}}{s+\mathit{\gamma }+\mathit{\gamma }\frac{{r}_{\mathrm{c}}}{{r}_{\mathrm{aH}}}},\end{array}$

with λ being the latent heat of vaporisation (J kg−1), Ea the actual evaporation (kg m−2 s−1), s the slope of the Clausius–Clapeyron curve relating air temperature with the saturation vapour pressure (Pa K−1), Rn the net radiation (W m−2), G the ground heat flux (W m−2), ρa the air density (kg m−3), γ the psychrometric constant (Pa K−1), cp the specific heat capacity of the air (J kg−1 K−1), VPD the vapour pressure deficit (Pa), raH the resistance of heat transfer to air (s m−1) and rc the surface resistance of water transfer (s m−1). While λ, cp, s and γ are air-temperature-dependent, raH is a complex function of wind speed, vegetation characteristics and the stability of the lower atmosphere (see Sect. 2.3). In most methods to estimate Ea or Ep, raH is estimated from a simple function of wind speed.

The Penman–Monteith equation is frequently used to calculate Ep by adjusting rc to its minimum value (the value under unstressed conditions). If the reference system is the actual system or a reference crop, rc is usually considered a fixed, constant value larger than zero (even if the soil is well-watered). In this study, both a universal, fixed value of rc for reference crops and a biome-specific constant value are used (see Sect. 2.5). When instead of a well-watered canopy a wet canopy (i.e. a canopy covered by water) is considered, rc=0 and Eq. (1) collapses to

$\begin{array}{}\text{(2)}& \mathit{\lambda }{E}_{\mathrm{p}}=\frac{s\left({R}_{\mathrm{n}}-G\right)+\frac{{\mathit{\rho }}_{\mathrm{a}}{c}_{p}\mathrm{VPD}}{{r}_{\mathrm{aH}}}}{s+\mathit{\gamma }}.\end{array}$

Equation (2) is often referred to as the Penman (1948) formulation and can be conveniently rearranged as $\mathit{\lambda }{E}_{\mathrm{p}}=\frac{s\left({R}_{\mathrm{n}}-G\right)}{s+\mathit{\gamma }}+\frac{{\mathit{\rho }}_{\mathrm{a}}{c}_{p}\mathrm{VPD}}{\left(s+\mathit{\gamma }\right){r}_{\mathrm{aH}}}$ to illustrate that Ep is driven by a radiative (left term) and an aerodynamic (right term) forcing (Brutsaert and Stricker, 1979).

### 2.1.2 Methods based on radiation and temperature

When the reference system is considered an idealised extensive area, or when radiative forcing is very dominant, the aerodynamic component of Eq. (2) may become negligible; thus the whole equation collapses to $\mathit{\lambda }{E}_{\mathrm{p}}=\frac{s\left({R}_{\mathrm{n}}-G\right)}{s+\mathit{\gamma }}$, which is commonly referred to as “equilibrium evaporation” (Slatyer and McIlroy, 1961). Priestley and Taylor (1972) analysed time series of open water and water-saturated crops and grasslands and found that the evaporation over these surfaces closely matched the equilibrium evaporation corrected by a multiplicative factor, commonly denoted as αPT:

$\begin{array}{}\text{(3)}& \mathit{\lambda }{E}_{\mathrm{p}}={\mathit{\alpha }}_{\mathrm{PT}}\frac{s\left({R}_{\mathrm{n}}-G\right)}{s+\mathit{\gamma }}.\end{array}$

This formulation is known as the Priestley and Taylor equation. Because usually a constant value of αPT is adopted, it assumes that the aerodynamic term in the Penman equation (Eq. 2) is a constant fraction of the radiative term. Typically, αPT=1.26 is considered, as estimated by Priestley and Taylor (1972) in their original experiments. In this study, we also include a biome-specific value to extend its applicability to all biomes (see Sect. 2.5). Since this method does not require wind speed or VPD as input, it is widely applied in hydrological models (Norman et al., 1995; Castellvi et al., 2001; Agam et al., 2010), remote sensing evaporation models (Norman et al., 1995; Fisher et al., 2008; Agam et al., 2010; Miralles et al., 2011) and drought monitoring methods (Anderson et al., 1997; Vicente-Serrano et al., 2018).

### 2.1.3 Methods based on radiation

Other studies such as Lofgren et al. (2011), or the more recent Milly and Dunne (2016), further simplified Eq. (3) to make it a linear function of the available energy by defining a constant multiplier here referred to as αMD:

$\begin{array}{}\text{(4)}& \mathit{\lambda }{E}_{\mathrm{p}}={\mathit{\alpha }}_{\mathrm{MD}}\left({R}_{\mathrm{n}}-G\right).\end{array}$

In the case of Milly and Dunne (2016) this equation was applied to climate model outputs based on a constant and universal value of αMD=0.8. On a daily scale, (RnG) expresses the total amount of energy available for evaporation, and the fraction of this energy that is actually used for evaporation is typically referred to “evaporative fraction”, or $\mathrm{EF}=\frac{\mathit{\lambda }{E}_{\mathrm{a}}}{\left(H+\mathit{\lambda }{E}_{\mathrm{a}}\right)}=\frac{\mathit{\lambda }{E}_{\mathrm{a}}}{\left({R}_{\mathrm{n}}-G\right)}$. From Eq. (4), it follows that the parameter αMD can be interpreted as the EF of the unstressed ecosystem. In this study, we test both the general value of αMD=0.8 and a biome-specific constant (see Sect. 2.5).

### 2.1.4 Methods based on temperature

Of the many empirical methods to estimate Ep, temperature-based methods have arguably been the most commonly used because of the availability of reliable air temperature data. For an overview of these methods, we refer to Oudin et al. (2005a). In this study, three methods are included. First, Pereira and Pruitt (2004) formulated a daily version of the well-known Thornthwaite (1948) equation:

$\begin{array}{}\text{(5a)}& & {T}_{\mathrm{eff}}<\mathrm{0}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathit{\lambda }{E}_{\mathrm{p}}=\mathrm{0},\text{(5b)}& & \mathrm{0}<{T}_{\mathrm{eff}}<\mathrm{26}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathit{\lambda }{E}_{\mathrm{p}}={\mathit{\alpha }}_{\mathrm{Th}}{\left(\frac{\mathrm{10}{T}_{\mathrm{eff}}}{I}\right)}^{b}\left(\frac{N}{\mathrm{360}}\right),\text{(5c)}& & \mathrm{26}<{T}_{\mathrm{eff}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathit{\lambda }{E}_{\mathrm{p}}=-c+d{T}_{\mathrm{eff}}-e{T}_{\mathrm{eff}}^{\mathrm{2}},\end{array}$

with Teff the effective temperature, based on maximum and minimal temperatures (see Sect. 2.5); αTh an empirical parameter (see below); I the yearly sum of (Ta_month∕5)1.514, with Ta_month the mean air temperature for each month; N the number of daylight hours; b a parameter depending on I and c; and d and e empirical constants (see Sect. 2.5). The general value of αTh=16 is often adopted; in this study, we will also calculate and apply a biome-specific value.

The second temperature-based formulation is that proposed by Oudin et al. (2005a), selected and developed after comparing 27 physically based and empirical methods with runoff data from 308 catchments:

$\begin{array}{}\text{(6a)}& & {T}_{\mathrm{a}}<\mathrm{5}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathit{\lambda }{E}_{\mathrm{p}}=\mathrm{0},\text{(6b)}& & {T}_{\mathrm{a}}>\mathrm{5}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathit{\lambda }{E}_{\mathrm{p}}=\frac{{R}_{e}}{{\mathit{\rho }}_{\mathrm{a}}}\frac{\left({T}_{\mathrm{a}}-\mathrm{5}\right)}{{\mathit{\alpha }}_{\mathrm{Ou}}},\end{array}$

with Ta the air temperature (C) and Re the top-of-atmosphere radiation (MJ m−2 day−1), depending on latitude and Julian day. Oudin et al. (2005a) suggested the use of αOu=100. This value will be used, in addition to a biome-specific value.

Finally, the third temperature-based method is the Hargreaves and Samani (1985) formulation, which includes minimum (Tmin) and maximum (Tmax) daily temperature, next to Ta and Re:

$\begin{array}{}\text{(7)}& \mathit{\lambda }{E}_{\mathrm{p}}={\mathit{\alpha }}_{\mathrm{HS}}{R}_{e}\left({T}_{\mathrm{a}}+\mathrm{17.8}\right)\sqrt{{T}_{\mathrm{max}}-{T}_{\mathrm{min}}},\end{array}$

with αHS a constant, normally assumed to equal 0.0023. As for the other methods, we additionally apply a biome-specific value. A detailed description of the calibration of all Ep methods is given in Sect. 2.5.

## 2.2 FLUXNET2015 database

The Tier 2 FLUXNET2015 database based on half-hourly or hourly measurements from eddy-covariance sensors is used to evaluate the different Ep formulations (http://fluxnet.fluxdata.org/data/fluxnet2015-dataset/, last access: 14 February 2019). Sites lacking at least one of the basic measurements required for our analysis (i.e. Rn; G; λEa; H; precipitation; wind speed, u; friction velocity, u*; Ta; and relative humidity, RH, or VPD) were not considered further. For latent and sensible heat fluxes, we used the data corrected by the Bowen ratio method. In this approach, the Bowen ratio is assumed to be correct, and the measured λEa and H are multiplied by a correction factor derived from a moving window method; see http://fluxnet.fluxdata.org/data/fluxnet2015-dataset/data-processing/ (last access: 14 February 2019) for a detailed description of this standard procedure. Nonetheless, taking the uncorrected λEa instead did not impact the main findings (not shown). For the main heat fluxes (G, H, λEa), medium and poor gap-filled data were masked out according to the flags provided by FLUXNET. As no quality flag was available for Rn measurements, the flag of the shortwave incoming radiation was used instead. All negative values for H or λEa were masked out, as these relate to periods of interception loss and condensation when accurate measurements are not guaranteed (Mizutani et al., 1997). Similarly, all negative values of Rn were masked out.

Finally, sub-daily measurements were aggregated to daytime composites (i) by applying a minimum threshold of 5 W m−2 of top-of-atmosphere incoming shortwave radiation and (ii) after excluding the first and last (half-)hours from these aggregates. Based on these daytime aggregates, the daytime means of s, γ and ρa were calculated using the parameterisation procedure described by Allen et al. (1998). We used Ta to calculate s. Only days in which more than 70 % of the data were measured directly were retained, and days with rainfall (between midnight and sunset) were removed from the analyses to avoid the effects of rainfall interception. Furthermore, only sites with at least 80 retained days were used for the further analysis. The global distribution of the final selection of sites is shown in Fig. 1 and detailed information about these sites is provided in Table S1 of the Supplement. The IGBP classification was used to assign a biome to each site.

Figure 1Location of the flux sites used in this study per biome. CRO: cropland; GRA: grasslands; DBF: deciduous broadleaf forest; EBF: evergreen broadleaf forest; ENF: evergreen needleleaf forest; MF: mixed forest; CSH: closed shrubland; WSA: woody savannah; SAV: savannah; OSH: open shrubland; WET: wetlands.

## 2.3 Calculation of resistance parameters

Estimates of raH are required for the Penman and Penman–Monteith equations. The resistance of heat transfer to air, raH, was calculated as

$\begin{array}{ll}{r}_{\mathrm{aH}}& =\frac{u}{{u}_{*}^{\mathrm{2}}}+\frac{\mathrm{1}}{k{u}_{*}}\left[\mathrm{ln}\left(\frac{{z}_{\mathrm{0}\mathrm{m}}}{{z}_{\mathrm{0}\mathrm{h}}}\right)+{\mathrm{\Psi }}_{\mathrm{m}}\left(\frac{z-d}{L}\right)-{\mathrm{\Psi }}_{\mathrm{m}}\left(\frac{{z}_{\mathrm{0}\mathrm{h}}}{L}\right)\right\\ \text{(8)}& & -{\mathrm{\Psi }}_{\mathrm{h}}\left(\frac{z-d}{L}\right)+{\mathrm{\Psi }}_{\mathrm{h}}\left(\frac{{z}_{\mathrm{0}\mathrm{h}}}{L}\right)],\end{array}$

in which k=0.41 is the von Kármán constant; z the (wind) sensor height (m); d the zero displacement height (m); z0m and z0h the roughness lengths for momentum and sensible heat transfer (m), respectively; L the Obukhov length (m); and Ψm(X) and Ψh(X) the Businger–Dyer stability functions for momentum and heat for the variable X, respectively. These were calculated based on the equations given by Garratt (1992) and Li et al. (2017) for stable, neutral and unstable conditions. Note that, in neutral and stable conditions, Ψm(X)=Ψh(X) and that ${\mathrm{\Psi }}_{\mathrm{m}}\left(\frac{z-d}{L}\right)-{\mathrm{\Psi }}_{\mathrm{m}}\left(\frac{{z}_{\mathrm{0}\mathrm{h}}}{L}\right)-{\mathrm{\Psi }}_{\mathrm{h}}\left(\frac{z-d}{L}\right)+{\mathrm{\Psi }}_{\mathrm{h}}\left(\frac{{z}_{\mathrm{0}\mathrm{h}}}{L}\right)=\mathrm{0}$. This is not the case for the unstable conditions that usually prevail during the daytime. Daytime averages of all variables were used as input in Eq. (8). The sensor height z was collected individually for each tower through online and literature research or through personal communication with the towers' principal investigator. The Obukhov length L was calculated as

$\begin{array}{}\text{(9)}& L=\frac{-{u}_{*}^{\mathrm{3}}{\mathit{\rho }}_{\mathrm{a}}{T}_{\mathrm{a}}\left(\mathrm{1}+\mathrm{0.61}{q}_{a}\right){c}_{p}}{kgH},\end{array}$

with qa being the specific humidity (kg kg−1) and g=9.81 m s−2 the gravitational acceleration.

The displacement height d and the roughness length for momentum flux z0m were estimated as a function of the canopy vegetation height (VH), as d=0.66 VH and z0m=0.1 VH (Brutsaert, 1982). VH was estimated from the flux tower measurements using the approach of Pennypacker and Baldocchi (2016):

$\begin{array}{}\text{(10)}& \mathrm{VH}=\frac{z}{\mathrm{0.66}+\mathrm{0.1}\mathrm{exp}\left(\frac{ku}{{u}_{*}}\right)}.\end{array}$

This equation was applied to the full (half-)hourly database and only when conditions were near-neutral ($|z/L|<\mathrm{0.01}$) and friction velocities lower than 1 standard deviation below the mean u* at each site. The daily VH was then aggregated by averaging out the half-hourly estimates to daily values, excluding the 20 % outliers of the data, and then calculating a 30-day-window moving average on the dataset. When not enough (half-)hourly vegetation height observations (<160) were available, the site was excluded from the analysis. This gave robust results for all remaining sites; an example of VH temporal development for a specific location is given in Fig. 2a. For homogeneous canopies, the VH calculated this way represents the true vegetation height. For savannah-like ecosystems, it corresponds to the vegetation height that the vegetation would have if it was represented by a single big leaf model.

Figure 2(a) Vegetation height dynamics in time (grey dots: half-hourly measurements; dark grey lines: daily mean vegetation height; red line: 30-day moving average, i.e. the final vegetation height dataset). (b) Relation between Stanton number (kB−1) and Reynolds number (Re). Both plots correspond to the woody savannah site of Santa Rita Mesquite, US-SRM (Arizona, USA).

The Stanton number (defined as $k{B}^{-\mathrm{1}}=\mathrm{ln}\left({z}_{\mathrm{0}\mathrm{m}}{z}_{\mathrm{0}\mathrm{h}}^{-\mathrm{1}}\right)$) was calculated by assuming that the surface aerodynamic temperature T0 (defined by $H={\mathit{\rho }}_{\mathrm{a}}{c}_{p}\frac{\left({T}_{\mathrm{0}}-{T}_{\mathrm{a}}\right)}{{r}_{\mathrm{aH}}}$) is equal to the radiative surface temperature Ts derived from the longwave fluxes. Then, through an iterative approach, an optimal value of z0h was obtained, using the following equations for T0 (Garratt, 1992) and Ts (Maes and Steppe, 2012):

$\begin{array}{ll}{T}_{\mathrm{0}}& ={T}_{\mathrm{a}}+\left(\frac{H}{k{u}_{*}{\mathit{\rho }}_{\mathrm{a}}{c}_{p}}\right)\left[\mathrm{ln}\left(\frac{z-d}{{z}_{\mathrm{0}\mathrm{h}}}\right)-{\mathrm{\Psi }}_{\mathrm{h}}\left(\frac{z-d}{L}\right)\right\\ \text{(11)}& & +{\mathrm{\Psi }}_{\mathrm{h}}\left(\frac{{z}_{\mathrm{0}\mathrm{h}}}{L}\right)],\end{array}$

$\begin{array}{}\text{(12)}& {T}_{\mathrm{s}}=\sqrt[\mathrm{4}]{\frac{{\mathrm{LW}}_{\mathrm{out}}-\left(\mathrm{1}-\mathit{\epsilon }\right){\mathrm{LW}}_{\mathrm{in}}}{\mathit{\sigma }\mathit{\epsilon }}},\end{array}$

with σ the Stefan–Boltzmann constant and ε the emissivity. The (half-)hourly data were used for this calculation. Following the approach of Li et al. (2017), only summertime data were used, and only when H was larger than 20 W m−2 and u* larger than 0.01 m s−1. Summertime was defined as those months in which the maximal daily value of H is at least 85 % of the 98th percentile of H based on the entire tower time series. In addition, (half-)hourly observations with counter-gradient heat fluxes were excluded from the analysis. For each observation, z0h was optimised by minimising the difference between T0 and Ts. Then, kB−1 was calculated at each site based on its relation with the observed Reynolds number (Re) by fitting the following function, based on the work by Li et al. (2017):

$\begin{array}{}\text{(13)}& k{B}^{-\mathrm{1}}={a}_{\mathrm{0}}+{a}_{\mathrm{1}}R{e}^{a\mathrm{2}}.\end{array}$

Note that Eq. (12) requires a value for ε, which is often assumed to be equal to 0.98 for all sites (e.g. Li et al., 2017; Rigden and Salvucci, 2015). Under the assumption that T0=Ts, ε can also be calculated separately per site. If H=0, it follows that T0=Ta and, from Eq. (12),

$\begin{array}{}\text{(14)}& \mathit{\epsilon }=\frac{{\mathrm{LW}}_{\mathrm{out}}-{\mathrm{LW}}_{\mathrm{in}}}{\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}}-{\mathrm{LW}}_{\mathrm{in}}}.\end{array}$

Here, ε was calculated for each site using (half-)hourly data, selecting those measurements where H was close to 0 ($-\mathrm{2} W m−2) and excluding rainy days as well as measurements in which the albedo (calculated as α= SWoutSW${}_{\mathrm{in}}^{-\mathrm{1}}$ with SWin the incoming and SWout the outgoing shortwave radiation) was above 0.4, to avoid influences of snow or ice. Negative estimates of ε were masked out, and the ε of the site was calculated as the mean, after excluding the outlying 20 % of the data. Equation (3) was applied both with a fixed ε of 0.98 and with the observed ε, and the ε value yielding the lowest RMSE in Eq. (12) was retained for each site. An example of such a function between kB−1 and Re is shown in Fig. 2b.

Finally, the surface resistance rc (s m−1) was calculated by inverting the Penman–Monteith equation as

$\begin{array}{}\text{(15)}& {r}_{\mathrm{c}}=\frac{s\left({R}_{\mathrm{n}}-G\right){r}_{\mathrm{aH}}+{\mathit{\rho }}_{\mathrm{a}}{c}_{p}\mathrm{VPD}}{\mathit{\gamma }\mathit{\lambda }{E}_{\mathrm{a}}}-\frac{\left(s+\mathit{\gamma }\right){r}_{\mathrm{aH}}}{\mathit{\gamma }}.\end{array}$

We converted the rc estimates to surface conductance gc (mm s−1) using ${g}_{\mathrm{c}}=\mathrm{1000}{r}_{\mathrm{c}}^{-\mathrm{1}}$; we will continue using gc (rather than rc) for the remainder of this paper. Note that the approach of calculating kB−1 directly requires a separate measurement of LWin and LWout, which was only available in 95 of the 107 selected sites. For the remaining sites, an alternative approach was used to calculate kB−1 (see Supplement).

## 2.4 Selection of unstressed days

We include two different approaches to identify a subset of measurements per eddy-covariance site in which the ecosystem was unstressed and provide the results for both methods. A first approach is based on soil moisture levels. For those sites where soil moisture measurements were available, the maximal soil moisture level for each site was determined as the 98th percentile of all soil moisture measurements. We then split up the dataset of each site in five equal-size classes based on evaporation percentiles. For each class, days having soil moisture levels belonging to the highest 5 % of soil moisture levels within each class were selected as unstressed days, but only if the soil moisture level of these selected days was above 75 % of the maximum soil moisture. This guaranties the sampling of unstressed evaporation during all seasons.

As soil moisture data were not available for a large number of sites, and because using soil moisture data does not exclude days in which evaporation may be constrained by other kinds of biotic or abiotic stress, a second approach for defining unstressed days was applied based on an energy balance criterion. We calculated the EF from the daytime λEa and H values and considered it as a direct proxy for evaporative stress; i.e. we assumed that, under unstressed conditions, a larger fraction of the available energy is used to evaporate water (Gentine et al., 2007, 2011a; Maes et al., 2011). This approach is similar to those of other Ep studies using eddy-covariance or lysimeter data, in which the Bowen ratio (e.g. Douglas et al., 2009) or the ratio of λEa and (SWin+ LWin) (Pereira and Pruitt, 2004) are used to define unstressed days. The unstressed record was comprised of all days with EF exceeding the 95th percentile EF threshold for each particular site or, if fewer than 15 days fulfilled this criterion, the 15 days with the highest EF. Consequently, we assume that at each site during at least 5 % of the days the conditions are such that evaporation is unstressed and Ea reflects Ep. The measured Ea from the identified unstressed days is further referred to as Eunstr (mm day−1) and used as reference data to evaluate the different Ep methods.

To assess whether the atmospheric conditions of the unstressed datasets are representative for the FLUXNET sites as such, a random bootstrap sample having the same number of records as the unstressed dataset was taken from the entire dataset of daily records. The mean, standard deviation, and 2nd and 98th percentile of SWin, Ta, VPD and u were calculated. This procedure was repeated 1000 times per site. A t test comparing the values of the unstressed subsample with those of the 1000 random samples was used to analyse whether the atmospheric conditions of the unstressed subsample were representative for the overall site conditions. This analysis was carried out for both methods to select unstressed days: the soil moisture threshold and the energy balance criterion.

## 2.5 Calculation and calibration of the different Ep methods

An overview of the different methods to calculate Ep is given in Table 1. If possible, three versions of each method were calculated: (1) a reference crop version, (2) a standard version and (3) a biome-specific version. The reference crop version calculates Ep for the reference short turf grass crop, with Rn and other properties of this reference crop as well. The standard version uses the same non-biome-specific parameters of the reference crop but considers Rn and other properties of the actual ecosystem. The biome-specific version of each method applies a calibration of the key parameters per biome (Table 1) and considers Rn and other properties of the actual ecosystem. These calibrated values per biome are based on the mean value of this key parameter for the unstressed dataset for each site, averaged out per biome type.

Table 1Overview of the different Ep methods used in this study and their calculation.

N is the number of daylight hours; Tmax and Tmin are the maximum and minimum daily air temperature; RHmax and RHmin are the maximum and minimum RH; SW* and LW* are the net shortwave and net longwave radiation; SWTOA is the shortwave incoming radiation at the top of the atmosphere; FAO-56 refers to the methodology described by Allen et al. (1998).

To estimate the radiation and crop properties of the reference crop versions, the equations described by Allen et al. (1998) in the FAO-56 method (Food and Agricultural Organization) were used and G was considered to be 0. Rn was calculated as

$\begin{array}{}\text{(16)}& {R}_{\mathrm{n}}={\mathrm{SW}}_{\mathrm{in}}\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{ref}}\right)+{\mathrm{LW}}^{*},\end{array}$

where αref=0.23 (Allen et al., 1998) and LW* is the net longwave radiation, calculated after Allen et al. (1998; Eq. 39, Sect. 3).

In the case of the reference crop version of the Penman–Monteith equation (Eq. 1), the FAO-56 method was used as described by Allen et al. (1998), with gc_ref fixed as 14.49 mm s−1 (corresponding with rc_ref=69 s m−1) and using Eq. (16) to calculate Rn. The standard version of the Penman–Monteith equation used observed (Rn, G, VPD) and calculated (s, γ, ρa, raH) daytime values as described in Sect. 2.2 in Eq. (1), and it also assumed gc_ref=14.49 mm s−1. The biome-specific version was calculated with the same data but used a biome-dependent value of gc. First, for each individual site, the unstressed gc was calculated as the mean of the gc values of the unstressed record (see Sect. 2.4). The mean value per biome was then calculated from these unstressed gc values. Regarding the Penman method (Eq. 2), the reference crop and standard versions were calculated using the same input data as for the Penman–Monteith methods; given Penman's consideration of no surface resistance, no biome-specific version was calculated.

The reference crop version of the Priestley and Taylor method was calculated from Eq. (3) with Rn from Eq. (16), s and γ from the FAO-56 calculations, and with αPT=1.26. The standard version used the same value for αPT but the observed daytime values for Rn and G. The biome-specific version followed a calibration of αPT similar to the gc_ref calculation. For each site, the unstressed αPT was calculated as the average αPT, obtained by solving Eq. (3) for αPT using the unstressed dataset. Finally, the mean per biome was calculated and used in the Ep estimation. Regarding the method by Milly and Dunne (2016) (Eq. 4), the reference crop, standard and biome-specific calculation were calculated accordingly, with Rn from Eq. (16) for the reference crop version, αMD=0.8 for the reference crop and standard version, and a calibrated αMD by biome type for the biome-specific version.

For the Thornthwaite, Oudin, and Hargreaves–Samani formulations (Eqs. 5–7), only the standard and biome-specific versions were calculated. The standard version of Thornthwaite's formulation used αTh=16. In the biome-specific version, this parameter was again calculated per site as the mean value of the unstressed records (e.g. Xu and Singh, 2001; Bautista et al., 2009), and then averaged per biome type. The effective temperature Teff was calculated as ${T}_{\mathrm{eff}}=\mathrm{0.36}\left(\mathrm{3}{T}_{\mathrm{max}}-{T}_{\mathrm{min}}\right)$ (Camargo et al., 1999). The parameter b was calculated as $b=\left(\mathrm{6.75}×{\mathrm{10}}^{-\mathrm{7}}{I}^{\mathrm{3}}\right)\left(\mathrm{7.71}×{\mathrm{10}}^{-\mathrm{7}}{I}^{\mathrm{2}}\right)+\mathrm{0.0179}I+\mathrm{0.492}$ and the parameters ce in Eq. (5c) were 415.85, 32.24 and 0.43, respectively (Pereira and Pruitt, 2004). For Oudin's temperature-based formulation, αOu=100 was taken for the standard version (Eq. 6). In the biome-specific version, this value was recalculated by calculating αOu for the unstressed records through Eq. (6), calculating the mean αOu per site and finally the biome-dependent αOu. Similarly, for the Hargreaves–Samani method, αHS=0.0023 is used for the standard version (Eq. 7), whereas, in the biome-specific version, this value was calculated using the unstressed records. Altogether, this exercise yielded a total of 17 different methods to estimate Ep, whose specificities are documented in Table 1.

The influence of climatic forcing data on Eunstr and on gc_ref, αPT and αMD was investigated. This was done by calculating for each individual site the correlations between the daily estimates of the atmospheric conditions and the daily values of the unstressed datasets. Analyses were then performed on these correlations of all sites.

3 Results

## 3.1 Representativeness of climate forcing data of unstressed datasets

Table S2 provides an overview of the analyses used to verify if the climatic forcing data of the unstressed subsets were representative for the atmospheric conditions of the sites as such. For the subsets of both unstressed criteria, atmospheric conditions were very representative for the site conditions, including for VPD. For the energy balance criterion, for instance, only at one site, the unstressed subset of the 98th percentile was significantly different from the random sampling-based simulations and in only two sites the 2nd percentile was significantly different from the simulations.

## 3.2 Key parameters by biome

We focus here on the parameter estimates of the unstressed record based on the energy balance criterion (Sect. 2.4). Of the full dataset, 107 flux sites meet all the selection criteria (i.e. at least 80 days without rainfall, good quality measurements of radiation and main fluxes, and at least 160 vegetation height observations – see Sect. 2.2 and 2.3). Despite considerable variation, gc_ref does not differ statistically across biomes, in contrast with αPT and αMD. Overall, croplands (CRO) are characterised by a higher measured Eunstr, which translates into the highest gc_ref, αPT and αMD of all biomes. Grasslands (GRA), deciduous broadleaf forest (DBF) and evergreen broadleaf forest (EBF) also have a relatively high gc_ref, αPT and αMD, whereas mixed forest (MF) and savannah ecosystems (closed shrubland, CSH; woody savannah, WSA; open shrublands, OSH; and savannah, SAV, ecosystems) are characterised by lower gc_ref, αPT and αMD. Only five sites (DE-Kli and IT-BCi, CRO; CA-SF3, OSH; AU-Rig, GRA; and AU-Wac, EBF) have mean values of αPT higher than the typically assumed 1.26 (Table 2). In contrast, 27 sites, including 9 croplands, have a mean value of αMD above 0.80, and 42 sites have a mean gc_ref above 14.49 mm s−1. Finally, wetlands (WET) show a large standard deviation of αPT and αMD (Table 2) due to their location in tropical, temperate and in arctic regions. The parameters of the three temperature-based methods differed significantly across biomes, but trends were different for each key parameter and did not always match those for gc_ref, αPT and αMD (Table 2).

Table 2Overview of the difference of the key parameters (gc_ref, αPT, αMD, αTh, αOu and αHS) during unstressed conditions per biome. The energy balance method was used for defining unstressed days (see Sect. 2.4; see Table 1 for definition of key parameters). The p value of the ANOVA test is given, as well as the mean ±1 SD for each biome. Different alphabetic superscripts indicate significantly differing means (Tukey post-hoc test: p<0.05). The number of sites per biome is given between brackets. Different parts are used to group biomes into broader ecosystem types (in descending order: croplands, grasslands, forests, savannah ecosystems and wetlands).

CRO: cropland; DBF: deciduous broadleaf forest; EBF: evergreen broadleaf forest; ENF: evergreen needleleaf forest; MF: mixed forest; CSH: closed shrubland; WSA: woody savannah; SAV: savannah; OSH: open shrubland; GRA: grasslands; WET: wetlands.

Table 3Influence of atmospheric conditions on Eunstr and on selected key parameters (gc_ref, αPT, αMD), based on unstressed days only defined using the energy balance criterion. Left: mean ±1 SD of the correlations of Eunstr, gc, αPT and αMD against the atmospheric conditions. Right: number of sites (out of a total of 107) with significant negative/positive correlations between Eunstr, αPT, gc_ref and αMD and the climate forcing variables.

* significantly different from 0.

Figure 3Histograms of correlations between the climate forcing variables and selected key parameters (a) Eunstr, (b) gc_ref, (c) αPT and (d) αMD measured in all flux tower sites, based on unstressed days only defined using the energy balance criterion.

Figure 4Scatterplot of the measured Eunstr versus Ep calculated with the different methods, based on unstressed days only, defined using the energy balance criterion. The discontinuous line is the 1:1 line.

Next, the effect of the atmospheric conditions on Eunstr and on the key parameters gc_ref, αPT and αMD of the unstressed dataset is investigated. Figure 3 gives the distribution of the correlations between the climatic variables (Rn, Ta, VPD, u and [CO2]) and Eunstr, gc_ref, αPT and αMD of the unstressed records at each site. We did not include αTh, αHS or αOu because the temperature-based methods did not perform well (see Sect. 3.3). Eunstr was strongly positively correlated with Rn, Ta and VPD in most sites, but less with u (Fig. 3a, Table 3). Considering all sites, the correlation between gc_ref and the climatic variables is not significantly different from zero for any climate variable, yet gc_ref is significantly negatively correlated with Tair and with VPD in 40 % and 45 % of the flux tower sites, respectively (Table 3, Fig. 3b). The two other parameters, αPT and αMD, appear less correlated to any of the climatic variables across all sites (Table 3). In the case of αMD, in particular, the distributions of the correlations against all climate forcing variables peak around zero (Fig. 3d): αMD is hardly influenced by Rn and is overall not dependent on u, Ta, [CO2] or VPD in most sites (Table 3).

## 3.3 Evaluation of different Ep methods

We first list the results of the analysis using the energy balance criterion for selecting the unstressed records (Sect. 2.4). The scatterplots of measured Eunstr versus estimated Ep based on the 17 different methods are shown in Fig. 4 for three sites belonging to different biomes. Despite the overall skill shown by the different Ep methods, considerable differences can be appreciated. In general, the methods designed for reference crops (PMr, Per, PTr, MDr) overestimate Eunstr and only two methods, MDB and PTB, match the measured Eunstr closely.

Table 4Mean correlations per biome between the measured Eunstr and the different Ep methods, based on unstressed days only defined using the energy balance criterion. The methods with the highest correlation per biome are highlighted in bold. Different parts are used to group biomes into broader ecosystem types (in descending order: croplands, grasslands, forests, savannah ecosystems and wetlands).

CRO: cropland; DBF: deciduous broadleaf forest; EBF: evergreen broadleaf forest; ENF: evergreen needleleaf forest; MF: mixed forest; CSH: closed shrubland; WSA: woody savannah; SAV: savannah; OSH: open shrubland; GRA: grasslands; WET: wetlands.

Table 4 gives the mean correlation per biome for each method. The results are very consistent and reveal that the highest correlations for nearly all biomes are obtained with the standard and biome-specific radiation-based method (MDs and MDb), closely followed by the standard and biome-dependent Priestley and Taylor method (PTs and PTb). Temperature-based methods have the lowest overall mean correlation as well as lower mean correlations per biome, with the Hargreaves–Samani method performing slightly better than the other two temperature-based methods. Note that the correlations are the same for the standard and biome-specific version in the case of the Priestley and Taylor (PTs and PTb), radiation-based (MDs and MDb), Oudin (Ous and Oub), and Hargreaves–Samani (HSs and HSb) methods (Table 4) – this is to be expected, as the only difference between the standard and biome-specific version of these methods is the value of their key parameters (αPT, αMD, αOu, αHS), which are multiplicative (see Eqs. 3, 4, 6 and 7). Differences are however reflected in the unbiased root mean square error (unRMSE) and mean bias – see Tables 5 and 6. The biome-specific versions of the radiation-based method (MDb) and of the Priestley and Taylor method (PTb) consistently have the lowest unRMSE for all biomes. Though the difference between these two methods is small, MDb performs slightly better. The standard Penman method (Pes) has the highest unRMSE. All reference crop methods (PMr, Per, PTr, MDr) have a mean unRMSE above 1 mm day−1, and the temperature-based methods (Ths, Ous, HSs, Thb, Oub and HSb) also have a relatively high unRMSE. Finally, bias estimates are given in Table 6. Again, MDb is overall the best-performing method (mean bias closest to 0 mm day−1), closely followed by the PTb method. Both methods consistently have the bias closest to zero among all biomes, except for wetlands. Most reference crop methods (PMr, Per, PTr, MDr), as well as Pes, overestimate Ep in all biomes.

Table 5Unbiased root mean square error (unRMSE) (in mm day−1) for the Ep methods per biome, based on unstressed days only defined using the energy balance criterion. The methods with the lowest unRMSE per biome are indicated in bold. Different parts are used to group biomes into broader ecosystem types (in descending order: croplands, grasslands, forests, savannah ecosystems and wetlands).

CRO: cropland; DBF: deciduous broadleaf forest; EBF: evergreen broadleaf forest; ENF: evergreen needleleaf forest; MF: mixed forest; CSH: closed shrubland; WSA: woody savannah; SAV: savannah; OSH: open shrubland; GRA: grasslands; WET: wetlands.

The use of soil water content as the criterion to select unstressed days (see Sect. 2.4) is explored. In total, 62 sites have soil water content data and meet the other selection criteria documented in Sect. 2.2. The results of this analysis are given in Tables S3–S5. To allow for a fair comparison, the same statistics have also been computed for just the same 62 tower sites with the energy balance criterion (Tables S6–S8). Using the soil moisture criterion, the correlations are overall lower and the results of the mean correlation, unRMSE and biases are less consistent. However, the overall performance ranking of the different models remains similar: PTb is the best-performing method with overall the highest mean correlation (R=0.84) and the lowest unRMSE (0.78 mm day−1) and a bias close to zero (−0.04), closely followed by the MDb method, with R=0.81, unRMSE = 0.89 and a mean bias of −0.12. More complex Penman-based models, and especially the empirical temperature-based formulations, show again a lower performance.

So far, all flux sites were used to calibrate the key parameters (Table 2) and those same sites were also used for the evaluation of the different methods. This was done to maximise the sample size. However, to avoid possible overfitting, we also repeated the analysis after separating calibration and validation samples. For each biome, two-thirds of the sites were randomly selected as calibration sites and one-third as validation sites. The key parameters were then calculated from the calibration subset and applied to estimate Ep of the biome-specific methods of the validation subset. This procedure was repeated 100 times and the mean correlation, unbiased RMSE and bias per biome were calculated. Results show no substantial differences in overall correlation, unRMSE and bias of each method, which are provided in Tables S9–S11 for completeness.

Table 6Mean bias (in mm day−1) for the Ep methods per biome, based on unstressed days only defined using the energy balance criterion. The best-performing method per biome is indicated in bold. Different parts are used to group biomes into broader ecosystem types (in descending order: croplands, grasslands, forests, savannah ecosystems and wetlands).

CRO: cropland; DBF: deciduous broadleaf forest; EBF: evergreen broadleaf forest; ENF: evergreen needleleaf forest; MF: mixed forest; CSH: closed shrubland; WSA: woody savannah; SAV: savannah; OSH: open shrubland; GRA: grasslands; WET: wetlands.

4 Discussion

## 4.1 Comparison of criteria to define unstressed days

We prioritised the energy balance over the soil moisture criterion to select unstressed days (see Sect. 2.4), because it can be applied to sites without soil moisture measurements and because it implicitly allows the exclusion of days in which the ecosystem is stressed for reasons other than soil moisture availability (e.g. insect plagues, phenological leaf-out, fires, heat and atmospheric dryness stress, nutrient limitations). In addition, soil moisture at specific depths can be a poor indicator of water stress, as rooting depth can vary and is not accurately measured (Powell et al., 2006; Douglas et al., 2009; Martínez-Vilalta et al., 2014). This is confirmed by our results: sampling unstressed days based on the energy-balance-based criterion resulted in higher correlations between Ep and Eunstr for all methods (Table S6 versus Table S3) and in lower unRMSE, with the exception of the temperature-based methods (Table S7 versus Table S4).

Nonetheless, it could be argued that, because the MD method assumes a constant evaporative fraction, the use of the evaporative fraction as a criterion for selecting unstressed days may favour the MD and even the closely related PT formulation. However, the soil moisture criterion adopted here provides an independent check of the results and confirms the robust and superior performance of the energy-driven PTb and MDb methods, independently of the framework chosen to select unstressed days. In the following discussion, the primary focus is on the results of selecting unstressed days based on the energy balance criterion.

## 4.2 Estimation of key ecosystem parameters

The resulting biome-specific values of the key parameters in Table 2 are within the range of values used in reference crop and standard applications of the models (Table 1), with the exception of αPT, which is typically lower than the frequently adopted value of 1.26. Other studies also found αPT values far below 1.26 but within the range of our study, mainly for forests (e.g. Shuttleworth and Calder, 1979; Viswanadham et al., 1991; Eaton et al., 2001; Komatsu, 2005) but also for tundra (Eaton et al., 2001) and grassland sites (Katerji and Rana, 2011) – see McMahon et al. (2013) for an overview. Our results and these studies demonstrate that the standard level of αPT=1.26 is close to the upper bound and will lead to an overestimation of Ep at most flux sites (Table 5).

## 4.3 Lower performance of the Penman–Monteith method

The poor performance of the PMr, PMs and PMb methods was relatively unexpected. Because the Penman–Monteith method incorporates the effects of Ta, VPD, Rn and u, it is often considered superior (e.g. Sheffield et al., 2012) and is even used as reference to evaluate other Ep formulations (e.g. Xu and Singh, 2002; Oudin et al., 2005b; Sentelhas et al., 2010). However, in studies dedicated to estimate Ea at eddy-covariance sites, in which gc is adjusted so it reflects the actual stress conditions in the ecosystem, the Penman–Monteith method has already been shown to perform worse than other simpler methods, such as the Priestley and Taylor equation (e.g. Ershadi et al., 2014; Michel et al., 2016). It is well known that its performance depends on the reliability of a wide range of input data and on the methods used to derive raH and gc (Singh and Xu, 1997; Dolman et al., 2014; Seiller and Anctil, 2016). In our case, the strict procedure followed to select the samples of 107 FLUXNET sites (see Sect. 2.2) ensured that all relevant variables were available, and the meteorological measurements were quality-checked. Hence, in our analysis, input quality is unlikely to be the cause of low performance.

We believe that the underlying assumption of a constant gc, typically adopted by PM methods (PMr, PMs, PMb) when estimating Ep, is a more likely explanation of the poor performance. PM was the only method in which the biome-specific calibration did not improve the performance. This is partially because of the large variation in gc_ref among the different flux sites of the same biome type (Table 2). In addition, of all the key parameters, gc_ref showed the largest mean relative standard deviation of the unstressed datasets of individual sites (results not shown). Surface conductance of the unstressed dataset was significantly negatively correlated with VPD in 45 % of the sites (Fig. 3b, Table 3). The relationship between gc and VPD for two of these sites is illustrated in Fig. 5. It is clear that gc of unstressed days (red dots) is always high for a given VPD, confirming the validity of the energy balance method. However, for these sites, it is also shown that gc of the unstressed days becomes very high when VPD becomes very low. As a consequence, the mean value of gc of the unstressed records, used to ultimately calculate gc_ref per biome type, is highly influenced by local VPD and is not necessarily a representative ecosystem property.

Figure 5Surface conductance gc as a function of vapour pressure deficit (VPD) of the regular and the unstressed datasets of two flux sites, (a) the evergreen needle forest Niwot Ridge Forest and (b) the open savannah woodland site Santa Rita Creosote.

The dependence of gc on VPD, even when soil moisture is not limiting, has been well studied (e.g. Jones, 1992; Granier et al., 2000; Sumner and Jacobs, 2005; Novick et al., 2016) and incorporated in most conventional stomatal or surface conductance models (e.g. Jarvis, 1976; Ball et al., 1987; Leuning, 1995). Yet, out of practical reasons, gc_ref is usually considered constant in Ep methods using the Penman–Monteith approach, with the PMr as the best illustration. Our data confirm that in unstressed conditions stomata open maximally only when VPD is very low. As such, the VPD dependence of gc smooths the impact of VPD in the Penman–Monteith equation – drops in VPD are compensated for by increases in gc, and vice versa, lowering the impact of VPD on Ea (Eq. 1). As such, assuming a constant gc_ref value overestimates the influence of VPD (and wind speed) on Ep.

Moreover, assuming a constant gc_ref value in the Penman–Monteith method also ignores the influence of CO2 levels on gc. As a result, Milly and Dunne (2016) found that the Penman–Monteith methods with constant gc_ref overpredicted Ep in models estimating future water use. Calibrating the sensitivity of gc_ref to VPD and [CO2] in the Penman–Monteith equation is outside the scope of this study, but could certainly improve Ep calculations – yet, it would further increase the complexity of the model. Finally, we note that the above discussion also applies to Penman's method: taking a wet canopy as reference in the Penman method (gc=∞ or rc=0) may not only overestimate Ep (Table 6) but also the influence of VPD and u on Ep.

## 4.4 Considerations regarding simple energy-based methods

The simpler Priestley and Taylor and radiation-based methods came forward as the best methods for assessing Ep with both criteria to define unstressed days. These observations are in agreement with studies highlighting radiation as the dominant driver of evaporation of saturated or unstressed ecosystems (e.g. Priestley and Taylor, 1972; Abtew, 1996; Wang et al., 2007; Song et al., 2017; Chan et al., 2018). They also agree with Douglas et al. (2009), who found that PT outperformed the PM method for estimating unstressed evaporation in 18 FLUXNET sites.

Both PTb and MDb are attractive from a modelling perspective, as they require minimum input data. However, this simplicity can also hold risks. The Priestley and Taylor method has been criticised for its implicit assumptions, which are also present in the MDb method. For instance, by not incorporating wind speed explicitly, it is assumed that the effect of u on Ep is somehow embedded within αPT (or αMD). Yet, several studies indicate that wind speeds are decreasing (“stilling”) globally (McVicar et al., 2008, 2012; Vautard et al., 2010). McVicar et al. (2012) also reported an associated decreasing trend in observed pan evaporation worldwide as well as in Ep calculated with the PMr method. With PT (or MD) methods, this trend cannot be captured. A similar criticism can be drawn with regards to the effect of [CO2] on stomatal conductance, water use efficiency and thus potential evaporation (Field et al., 1995). A separate question is whether more complex Ep methods that incorporate the effects of u, [CO2] or VPD explicitly do this correctly; the above-mentioned issues about the fixed parameterisation of the Penman–Monteith methods for estimating Ep indicate that this may typically not be the case.

Regarding the non-explicit consideration of u by simpler methods, our records show a limited effect of u on Ea and Ep, even when considering larger temporal scales. Of the 16 flux towers with at least 10 years of evaporation data, we calculated the yearly average Ea as well as the annual mean climatic forcing variables. Yearly averages were calculated from monthly averages, which in turn were calculated if at least three daytime measurements were available. Despite a relatively large mean standard deviation in yearly average u of 7.0 %, yearly average u was not significantly correlated with Ea in any of these sites. In contrast, yearly average Rn was positively correlated with yearly average Ea in 7 of the 16 sites, with comparable mean standard deviation in annual Rn (8.5 %). Moreover, looking at all individual towers and using the daily estimates, neither αMD nor αPT were correlated with u (Fig. 3c and d, Table 3). In fact, since αMD appears hardly affected by any climatic variable, and given the relatively small range in αMD values within each biome (Table 2), it appears that αMD is a robust biome property that can be used in the seamless application of these methods at a global scale. The robustness of αMD as a biome property is furthermore confirmed by the analysis with independent calibration and validation sites, which hardly affected the unRMSE and bias (Tables S10 and S11).

## 4.5 Use of available energy under stress conditions

The best-performing methods rely heavily on measurements of available energy (RnG) (Eqs. 2 and 3). In Sect. 3.2, all Ep calculations used available energy obtained during unstressed conditions. The question is whether Ep can also be calculated correctly using the actual Rn and G when the ecosystems are not unstressed. As mentioned in Sect. 1, there is discussion on whether SWout and LWout should be considered forcing variables or ecosystem responses (e.g. Lhomme, 1997; Lhomme and Guilioni, 2006; Shuttleworth, 1993). Among other considerations, it is clear that Ts will be lower if vegetation is healthy and soils are well-watered (Maes and Steppe, 2012), which results in lower LWout and higher Rn. Therefore, while using the observations of SWin and LWin as forcing variables in the computation of Ep is defendable (despite the potential atmospheric feedbacks that may derive from the consideration of unstressed conditions), we agree that SWout, LWout and G should ideally reflect unstressed rather than actual conditions to estimate Ep. Note that, in that case, Ep deviates from ${E}_{{\mathrm{p}}_{\mathrm{0}}}$ in the complementarity relationship, in which atmospheric feedbacks affecting incoming radiation or VPD are implicitly considered (Kahler and Brutsaert, 2006).

A method to derive the unstressed estimates of SWin (through the albedo, α), LWin (through Ts) and G under stressed conditions is presented in the Sect. S2 and is based on the MDb method and on flux data of the unstressed datasets. We further refer to this method as the “unstressed” Ep. It requires a large amount of input data and is not practically applicable at a global scale. Comparing the mean unstressed Ep with the “actual” Ep (i.e. Ep calculated from the actual Rn and G) for all sites reveals that the actual Ep is 8.2±10.1 % lower than the unstressed Ep. There are no significant differences between biomes, but the distribution of the underestimation is left-skewed and the actual Ep is more than 10 % lower than the unstressed Ep in 22 % of the sites (Fig. 6). The main reason explaining about 65 % of the difference between the actual and the unstressed Ep is the difference in Ts. Assuming that the unstressed Ts can be estimated as the mean of Ta and the actual Ts results in a straightforward alternative to approximate the unstressed Ep with only data of Ta and radiation:

$\begin{array}{ll}{E}_{\mathrm{p}}& ={\mathit{\alpha }}_{\mathrm{MD}}\left(\left(\mathrm{1}-\mathit{\alpha }\right){\mathrm{SW}}_{\mathrm{in}}+\mathit{\epsilon }{\mathrm{LW}}_{\mathrm{in}}-\mathrm{0.5}\mathit{\epsilon }{\mathrm{LW}}_{\mathrm{out}}\right\\ \text{(17)}& & -\mathrm{0.5}\mathit{\epsilon }\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}}-G).\end{array}$

This approach was also tested and resulted in a mean underestimation of Ep of 2.6±5.8 % compared to the unstressed Ep, with a mean median value at −0.1 % (Fig. 6). Given the low error and the straightforward calculation, we recommend this method to calculate Ep at global scales.

Figure 6Comparison of the Ep calculated with a modelled method calculating (RnG) for unstressed conditions (Sect. S2) using flux tower data, and the actually observed (RnG) (red) or a simplified correction of (RnG) using Ta (Eq. 17) (green). Negative y values indicate a lower estimation of Ep compared to the modelled method. For each distribution, the mean and median are indicated with a full and dashed line, respectively.

Figure 7Seasonal evolution of Ea (top, green), Ep (top, red) and the S factor ($S={E}_{\mathrm{a}}{E}_{\mathrm{p}}^{-\mathrm{1}}$, confined between 0 and 1) for two ecosystems: (a) a grassland crop, Sturt Plains in Australia; and (b) a deciduous broadleaf forest, Ohio Oak forest in the USA. Ep was calculated with the MDb method and using the tower-based correction of (RnG) as presented in S2 of the Supplement.

Figure 7 gives an example of the seasonal evolution of Ea and Ep and the S factor ($S={E}_{\mathrm{a}}{E}_{\mathrm{p}}^{-\mathrm{1}}$) in a grassland (Fig. 7a) and a deciduous forest site (Fig. 7b). The short growing season in the grassland site, when Ea is close to Ep and values of S are close to 1, stands in clear contrast with the winter period, when grasses have died off and Ea and consequently also S are very low. In the relatively wet broadleaf forest, Ea and Ep follow a similar seasonal cycle. In winter, when total evaporation is limited to soil evaporation, S is very low; in spring, when leaves are still developing, Ea lags Ep. In summer, S remains high and close to one.

5 Conclusion

Based on a large sample of eddy-covariance sites from the FLUXNET2015 database, we demonstrated a higher potential of radiation-driven methods calibrated by biome type to estimate Ep than of more complex Penman–Monteith approaches or empirical temperature-based formulations. This was consistent across all 11 biomes represented in the database and for two different criteria to identify unstressed days, one based on soil moisture and the other on evaporative fractions. Our analyses also showed that the key parameters required to apply the higher-performance radiation-driven methods are relatively insensitive to climate forcing. This makes these methods robust for incorporation into global offline models, e.g. for hydrological applications. Finally, we conclude that, at the ecosystem scale, Penman–Monteith methods for estimating Ep should only be prioritised if the unstressed stomatal conductance is calculated dynamically and high accuracy observations from the wide palate of required forcing variables are available.

Data availability
Data availability.

The FLUXNET2015 dataset can be downloaded from http://fluxnet.fluxdata.org/data/fluxnet2015-dataset/ (FLUXNET, 2016). The main script for calculating potential evaporation with the different method as well as the daily flux data of one site (AU-How), for which permission of distribution was granted, is available as supplement. For further questions, we ask readers to contact the corresponding author at wh.maes@ugent.be.

Appendix A

Table A1Abbreviation list.

Supplement
Supplement.

The supplement contains a description of the calculation method of kB−1 for sites where radiance fluxes are not separately measured (Sect. S1), a description of the method to estimate unstressed Ep (see Sect. 4.5) from flux data (Sect. S2) and several tables (Sect. S3), including an overview of the FLUXNET sites used in the study (Table S1); a table on the representativeness of the climate forcing conditions of the unstressed datasets for the full dataset (Table S2); and correlation, unbiased RMSE and bias tables for different selection criteria (Tables S3–S11). The supplement related to this article is available online at: https://doi.org/10.5194/hess-23-925-2019-supplement.

Author contributions
Author contributions.

WHM and DGM designed the research. PG and NECV assisted in developing the optimal method for analysing all flux tower data. WHM performed the calculations and analyses. WHM and DGM wrote the manuscript, with contributions from PG and NECV.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

This study was funded by the Belgian Science Policy Office (BELSPO) in the frame of the STEREO III program project STR3S (SR/02/329). Diego G. Miralles acknowledges support from the European Research Council (ERC) under grant agreement no. 715254 (DRY-2-DRY). 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 harmonisation were 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. The Atqasuk and Ivotuk towers in Alaska were supported by the Office of Polar Programs of the National Science Foundation (NSF) awarded to DZ, WCO and DAL (award number 1204263) with additional logistical support funded by the NSF Office of Polar Programs; by the Carbon in Arctic Reservoirs Vulnerability Experiment (CARVE), an Earth Ventures (EV-1) investigation, under contract with the National Aeronautics and Space Administration; and by the ABoVE (NNX15AT74A; NNX16AF94A) program. The OzFlux network is supported by the Australian Terrestrial Ecosystem Research Network (TERN, http://www.tern.org.au, last access: 14 February 2019). We would also like to thank the anonymous reviewers as well as Giovanni Ravazzani for their useful suggestions on an earlier version of the manuscript.

Edited by: Alberto Guadagnini
Reviewed by: Giovanni Ravazzani and one anonymous referee

References

Abtew, W.: Evapotranspiration measurements and modelling for three wetland systems in South Florida, J. Am. Water Resour. Assoc., 32, 465–473, https://doi.org/10.1111/j.1752-1688.1996.tb04044.x, 1996.

Agam, N., Kustas, W. P., Anderson, M. C., Norman, J. M., Colaizzi, P. D., Howell, T. A., Prueger, J. H., Meyers, T. P., and Wilson, T. B.: Application of the Priestley–Taylor Approach in a Two-Source Surface Energy Balance Model, J. Hydrometeorol., 11, 185–198, https://doi.org/10.1175/2009jhm1124.1, 2010.

Allen, R. G., Pereira, L. S., Raes, D., and Smith, M.: Crop Evapotranspiration – guidelines for computing crop water requirements, FAO Irrigation and Drainage Paper No. 56, Food and Agriculture Organization of the United Nations (FAO), Rome, Italy, 1998.

Aminzadeh, M., Roderick, M. L., and Or, D.: A generalized complementary relationship between actual and potential evaporation defined by a reference surface temperature, Water Resour. Res., 52, 385–406, https://doi.org/10.1002/2015WR017969, 2016.

Anderson, M. C., Norman, J. M., Diak, G. R., Kustas, W. P., and Mecikalski, J. R.: A two-source time-integrated model for estimating surface fluxes using thermal infrared remote sensing, Remote Sens.f Environ., 60, 195–216, 1997.

Bai, P., Liu, X. M., Yang, T. T., Li, F. D., Liang, K., Hu, S. S., and Liu, C. M.: Assessment of the influences of different potential evapotranspiration inputs on the performance of monthly hydrological models under different climatic conditions, J. Hydrometeorol., 17, 2259–2274, https://doi.org/10.1175/jhm-d-15-0202.1, 2016.

Ball, J. T., Woodrow, I. E., and Berry, J. A.: A model predicting stomatal conductance and its contribution to the control of photosynthesis under different environmental conditions, in: Progress in photosynthesis research: Volume 4 proceedings of the VIIth international congress on photosynthesis providence, 10–15 August 1986, Rhode Island, USA, edited by: Biggins, J., Springer Netherlands, Dordrecht, 221–224, 1987.

Barton, I. J.: A parameterization of the evaporation from nonsaturated surfaces, J. Appl. Metereol., 18, 43–47, 1979.

Bautista, F., Baustista, D., and Delgado-Carranza, C.: Calibration of the equations of Hargreaves and Thornthwaite to estimate the potential evapotranspiration in semi-arid and subhumid tropical climates for regional applications, Atmósfera, 22, 331–348, 2009.

Bouchet, R. J.: Evapotranspiration réelle, évapotranspiration potentielle et production agricole, L'eau et la production végétale, INRA, Paris, 1964.

Brutsaert, W.: Evaporation into the atmosphere. Theory, history and applications, D. Reidel Publishing Company, Dordrecht, 299 pp., 1982.

Brutsaert, W. and Stricker, H.: An advection-aridity approach to estimate actual regional evapotranspiration, Water Resour. Res., 15, 443–450, https://doi.org/10.1029/WR015i002p00443, 1979.

Camargo, A. P., R., M. F., Sentelhas, P. C., and Picini, A. G.: Adjust of the Thornthwaite's method to estimate the potential evapotranspiration for arid and superhumid climates, based on daily temperature amplitude, Revista Brasileira de Agrometeorologia, 7, 251–257, 1999.

Castellvi, F., Stockle, C., Perez, P., and Ibanez, M.: Comparison of methods for applying the Priestley–Taylor equation at a regional scale, Hydrol. Process., 15, 1609–1620, 2001.

Chan, F. C., Arain, M. A., Khomik, M., Brodeur, J. J., Peichl, M., Restrepo-Coupe, N., Thorne, R., Beamesderfer, E., McKenzie, S., and Xu, B.: Carbon, water and energy exchange dynamics of a young pine plantation forest during the initial fourteen years of growth, Forest Ecol. Manage., 410, 12–26, 2018.

Chen, D., Gao, G., Xu, C.-Y., Guo, J., and Ren, G.: Comparison of the Thornthwaite method and pan data with the standard Penman–Monteith estimates of reference evapotranspiration in China, Clim. Res., 28, 123–132, 2005.

Crago, R. and Crowley, R.: Complementary relationships for near-instantaneous evaporation, J. Hydrol., 300, 199–211, https://doi.org/10.1016/j.jhydrol.2004.06.002, 2005.

Dai, A.: Drought under global warming: a review, Wiley Interdisciplin. Rev.: Clim. Change, 2, 45–65, https://doi.org/10.1002/wcc.81, 2011.

Dolman, A. J., Miralles, D. G., and de Jeu, R. A. M.: Fifty years since Monteith's 1965 seminal paper: the emergence of global ecohydrology, Ecohydrology, 7, 897–902, https://doi.org/10.1002/eco.1505, 2014.

Donohue, R. J., McVicar, T. R., and Roderick, M. L.: Assessing the ability of potential evaporation formulations to capture the dynamics in evaporative demand within a changing climate, J. Hydrol., 386, 186–197, https://doi.org/10.1016/j.jhydrol.2010.03.020, 2010.

Douglas, E. M., Jacobs, J. M., Sumner, D. M., and Ray, R. L.: A comparison of models for estimating potential evapotranspiration for Florida land cover types, J. Hydrol., 373, 366–376, https://doi.org/10.1016/j.jhydrol.2009.04.029, 2009.

Eaton, A. K., Rouse, W. R., Lafleur, P. M., Marsh, P., and Blanken, P. D.: Surface energy balance of the Western and Central Canadian Subarctic: variations in the energy balance among five major terrain types, J. Climate, 14, 3692–3703, https://doi.org/10.1175/1520-0442(2001)014<3692:SEBOTW>2.0.CO;2, 2001.

Eltahir, E. A.: A soil moisture–rainfall feedback mechanism: 1. Theory and observations, Water Resour. Res., 34, 765–776, 1998.

Ershadi, A., McCabe, M. F., Evans, J. P., Chaney, N. W., and Wood, E. F.: Multi-site evaluation of terrestrial evaporation models using FLUXNET data, Agr. Forest Meteorol., 187, 46–61, https://doi.org/10.1016/j.agrformet.2013.11.008, 2014.

Field, C. B., Jackson, R. B., and Mooney, H. A.: Stomatal responses to increased CO2 – Implications from the plant to the global scale, Plant Cell Environ., 18, 1214–1225, 1995.

Fisher, J. B., Tu, K. P., and Baldocchi, D. D.: Global estimates of the land–atmosphere water flux based on monthly AVHRR and ISLSCP-II data, validated at 16 FLUXNET sites, Remote Sens. Environ., 112, 901–919, https://doi.org/10.1016/j.rse.2007.06.025, 2008.

Fisher, J. B., Whittaker, R. J., and Malhi, Y.: ET come home: potential evapotranspiration in geographical ecology, Global Ecol. Biogeogr., 20, 1–18, https://doi.org/10.1111/j.1466-8238.2010.00578.x, 2011.

FLUXNET: FLUXNET2015 Dataset, available at: http://fluxnet.fluxdata.org/data/fluxnet2015-dataset/ (last access: 14 February 2019), 2016.

Garratt, J. R.: The atmospheric boundary layer, Cambridge University Press, Cambridge, New York, 1992.

Gentine, P., Entekhabi, D., Chehbouni, A., Boulet, G., and Duchemin, B.: Analysis of evaporative fraction diurnal behaviour, Agr. Forest Meteorol., 143, 13–29, https://doi.org/10.1016/j.agrformet.2006.11.002, 2007.

Gentine, P., Entekhabi, D., and Polcher, J.: The diurnal behavior of evaporative fraction in the soil–vegetation–atmospheric boundary layer continuum, J. Hydrometeorol., 12, 1530–1546, https://doi.org/10.1175/2011JHM1261.1, 2011a.

Gentine, P., Polcher, J., and Entekhabi, D.: Harmonic propagation of variability in surface energy balance within a coupled soil–vegetation–atmosphere system, Water Resour. Res., 47, W05525, https://doi.org/10.1029/2010WR009268, 2011b.

Granger, R. J.: A complementary relationship approach for evaporation from nonsaturated surfaces, J. Hydrol., 111, 31–38, https://doi.org/10.1016/0022-1694(89)90250-3, 1989.

Granier, A., Biron, P., and Lemoine, D.: Water balance, transpiration and canopy conductance in two beech stands, Agr. Forest Meteorol., 100, 291–308, 2000.

Guo, D., Westra, S., and Maier, H. R.: Sensitivity of potential evapotranspiration to changes in climate variables for different Australian climatic zones, Hydrol. Earth Syst. Sci., 21, 2107–2126, https://doi.org/10.5194/hess-21-2107-2017, 2017.

Hargreaves, G. H. and Samani, Z.: Reference crop evapotranspiration from temperature, Appl. Eng. Agricult., 1, 96–99, https://doi.org/10.13031/2013.26773, 1985.

Jacobs, J. M., Anderson, M. C., Friess, L. C., and Diak, G. R.: Solar radiation, longwave radiation and emergent wetland evapotranspiration estimates from satellite data in Florida, USA, Hydrolog. Sci. J., 49, 476, https://doi.org/10.1623/hysj.49.3.461.54352, 2004.

Jarvis, P. G.: Interpretation of variations in leaf water potential and stomatal conductance found in canopies in field, Philos. T. Roy. Soc. Lond. B, 273, 593–610, 1976.

Jones, H. G.: Plants and microclimate. A quantitative approach to environmental plant physiology, 2nd Edn., Cambridge University Press, Cambridge, UK, 1992.

Kahler, D. M. and Brutsaert, W.: Complementary relationship between daily evaporation in the environment and pan evaporation, Water Resour. Res., 42, W05413, https://doi.org/10.1029/2005WR004541, 2006.

Katerji, N. and Rana, G.: Crop reference evapotranspiration: a discussion of the concept, analysis of the process and validation, Water Resour. Manage., 25, 1581–1600, https://doi.org/10.1007/s11269-010-9762-1, 2011.

Kay, A. L. and Davies, H. N.: Calculating potential evaporation from climate model data: A source of uncertainty for hydrological climate change impacts, J. Hydrol., 358, 221–239, https://doi.org/10.1016/j.jhydrol.2008.06.005, 2008.

Kingston, D. G., Todd, M. C., Taylor, R. G., Thompson, J. R., and Arnell, N. W.: Uncertainty in the estimation of potential evapotranspiration under climate change, Geophys. Res. Lett., 36, L20403, https://doi.org/10.1029/2009GL040267, 2009.

Komatsu, H.: Forest categorization according to dry-canopy evaporation rates in the growing season: comparison of the Priestley–Taylor coefficient values from various observation sites, Hydrol. Process., 19, 3873–3896, https://doi.org/10.1002/hyp.5987, 2005.

Leuning, R.: A critical appraisal of a combined stomatal-photosynthesis model for C3 plants, Plant Cell Environ., 18, 339–355, https://doi.org/10.1111/j.1365-3040.1995.tb00370.x, 1995.

Lhomme, J.-P.: Towards a rational definition of potential evaporation, Hydrol. Earth Syst. Sci., 1, 257–264, https://doi.org/10.5194/hess-1-257-1997, 1997.

Lhomme, J. P. and Guilioni, L.: Comments on some articles about the complementary relationship, J. Hydrol., 323, 1–3, https://doi.org/10.1016/j.jhydrol.2005.08.014, 2006.

Li, D., Rigden, A., Salvucci, G., and Liu, H.: Reconciling the Reynolds number dependence of scalar roughness length and laminar resistance, Geophys. Res. Lett., 44, 3193–3200, https://doi.org/10.1002/2017GL072864, 2017.

Li, S., Kang, S., Zhang, L., Zhang, J., Du, T., Tong, L., and Ding, R.: Evaluation of six potential evapotranspiration models for estimating crop potential and actual evapotranspiration in arid regions, J. Hydrol., 543, 450–461, https://doi.org/10.1016/j.jhydrol.2016.10.022, 2016.

Lofgren, B. M., Hunter, T. S., and Wilbarger, J.: Effects of using air temperature as a proxy for potential evapotranspiration in climate change scenarios of Great Lakes basin hydrology, J. Great Lakes Res., 37, 744–752, https://doi.org/10.1016/j.jglr.2011.09.006, 2011.

Maes, W. H. and Steppe, K.: Estimating evapotranspiration and drought stress with ground-based thermal remote sensing in agriculture: a review, J. Exp. Bot., 63, 4671–4712, 2012.

Maes, W. H., Pashuysen, T., Trabucco, A., Veroustraete, F., and Muys, B.: Does energy dissipation increase with ecosystem succession? Testing the ecosystem exergy theory combining theoretical simulations and thermal remote sensing observations, Ecol. Model., 23–24, 3917–3941, 2011.

Martens, B., Miralles, D. G., Lievens, H., van der Schalie, R., de Jeu, R. A. M., Fernández-Prieto, D., Beck, H. E., Dorigo, W. A., and Verhoest, N. E. C.: GLEAM v3: satellite-based land evaporation and root-zone soil moisture, Geosci. Model Dev., 10, 1903–1925, https://doi.org/10.5194/gmd-10-1903-2017, 2017.

Martínez-Vilalta, J., Poyatos, R., Aguadé, D., Retana, J., and Mencuccini, M.: A new look at water transport regulation in plants, New Phytol., 204, 105–115, https://doi.org/10.1111/nph.12912, 2014.

McMahon, T. A., Peel, M. C., Lowe, L., Srikanthan, R., and McVicar, T. R.: Estimating actual, potential, reference crop and pan evaporation using standard meteorological data: a pragmatic synthesis, Hydrol. Earth Syst. Sci., 17, 1331–1363, https://doi.org/10.5194/hess-17-1331-2013, 2013.

McVicar, T. R., Van Niel, T. G., Li, L. T., Roderick, M. L., Rayner, D. P., Ricciardulli, L., and Donohue, R. J.: Wind speed climatology and trends for Australia, 1975–2006: Capturing the stilling phenomenon and comparison with near-surface reanalysis output, Geophys. Res. Lett., 35, L20403, https://doi.org/10.1029/2008gl035627, 2008.

McVicar, T. R., Roderick, M. L., Donohue, R. J., Li, L. T., Van Niel, T. G., Thomas, A., Grieser, J., Jhajharia, D., Himri, Y., Mahowald, N. M., Mescherskaya, A. V., Kruger, A. C., Rehman, S., and Dinpashoh, Y.: Global review and synthesis of trends in observed terrestrial near-surface wind speeds: Implications for evaporation, J. Hydrol., 416, 182–205, https://doi.org/10.1016/j.jhydrol.2011.10.024, 2012.

Michel, D., Jiménez, C., Miralles, D. G., Jung, M., Hirschi, M., Ershadi, A., Martens, B., McCabe, M. F., Fisher, J. B., Mu, Q., Seneviratne, S. I., Wood, E. F., and Fernández-Prieto, D.: The WACMOS-ET project – Part 1: Tower-scale evaluation of four remote-sensing-based evapotranspiration algorithms, Hydrol. Earth Syst. Sci., 20, 803–822, https://doi.org/10.5194/hess-20-803-2016, 2016.

Milly, P. C. D. and Dunne, K. A.: Potential evapotranspiration and continental drying, Nat. Clim. Change, 6, 946–949, https://doi.org/10.1038/nclimate3046, 2016.

Miralles, D. G., Holmes, T. R. H., De Jeu, R. A. M., Gash, J. H., Meesters, A. G. C. A., and Dolman, A. J.: Global land-surface evaporation estimated from satellite-based observations, Hydrol. Earth Syst. Sci., 15, 453–469, https://doi.org/10.5194/hess-15-453-2011, 2011.

Mizutani, K., Yamanoi, K., Ikeda, T., and Watanabe, T.: Applicability of the eddy correlation method to measure sensible heat transfer to forest under rainfall conditions, Agr. Forest Meteorol., 86, 193–203, https://doi.org/10.1016/S0168-1923(97)00012-9, 1997.

Monteith, J. L.: Evaporation and environment, Symp. Soc. Exp. Biol., 19, 205–234, 1965.

Morton, F. I.: Operational estimates of areal evapotranspiration and their significance to the science and practice of hydrology, J. Hydrol., 66, 1–76, https://doi.org/10.1016/0022-1694(83)90177-4, 1983.

Mu, Q., Heinsch, F. A., Zhao, M., and Running, S. W.: Development of a global evapotranspiration algorithm based on MODIS and global meteorology data, Remote Sens. Environ., 111, 519–536, https://doi.org/10.1016/j.rse.2007.04.015, 2007.

Norman, J. M., Kustas, W. P., and Humes, K. S.: Source approach for estimating soil and vegetation energy fluxes in observations of directional radiometric surface temperature, Agr. Forest Meteorol., 77, 263–293, 1995.

Novick, K. A., Ficklin, D. L., Stoy, P. C., Williams, C. A., Bohrer, G., Oishi, A. C., Papuga, S. A., Blanken, P. D., Noormets, A., Sulman, B. N., Scott, R. L., Wang, L., and Phillips, R. P.: The increasing importance of atmospheric demand for ecosystem water and carbon fluxes, Nat. Clim. Change, 6, 1023–1027, https://doi.org/10.1038/nclimate3114, 2016.

Oudin, L., Hervieu, F., Michel, C., Perrin, C., Andréassian, V., Anctil, F., and Loumagne, C.: Which potential evapotranspiration input for a lumped rainfall–runoff model: Part 2 – Towards a simple and efficient potential evapotranspiration model for rainfall–runoff modelling, J. Hydrol., 303, 290–306, https://doi.org/10.1016/j.jhydrol.2004.08.026, 2005a.

Oudin, L., Michel, C., and Anctil, F.: Which potential evapotranspiration input for a lumped rainfall–runoff model: Part 1 – Can rainfall–runoff models effectively handle detailed potential evapotranspiration inputs?, J. Hydrol., 303, 275–289, https://doi.org/10.1016/j.jhydrol.2004.08.025, 2005b.

Palmer, C., Fisher, J. B., Mallick, K., and Lee, J.: The Potential of Potential Evapotranspiration, Fall Meeting, 3–7 December 2012, American Geophysical Union, San Francisco, 2012.

Penman, H. L.: Natural evaporation from open water, bare soil and grass, P. Roy. Soc. Lond. A, 193, 120–145, https://doi.org/10.1098/rspa.1948.0037, 1948.

Penman, H. L.: Vegetation and Hydrology, in: Soil Science, Vol. 96, 357, 1963.

Pennypacker, S. and Baldocchi, D.: Seeing the Fields and Forests: Application of Surface-Layer Theory and Flux-Tower Data to Calculating Vegetation Canopy Height, Bound.-Lay. Meteorol., 158, 165–182, https://doi.org/10.1007/s10546-015-0090-0, 2016.

Pereira, A. R. and Pruitt, W. O.: Adaptation of the Thornthwaite scheme for estimating daily reference evapotranspiration, Agr. Water Manage., 66, 251–257, https://doi.org/10.1016/j.agwat.2003.11.003, 2004.

Pettijohn, J. C. and Salvucci, G. D.: A new two-dimensional physical basis for the complementary relation between terrestrial and pan evaporation, J. Hydrometeorol., 10, 565–574, https://doi.org/10.1175/2008jhm1026.1, 2009.

Powell, T. L., Bracho, R., Li, J., Dore, S., Hinkle, C. R., and Drake, B. G.: Environmental controls over net ecosystem carbon exchange of scrub oak in central Florida, Agr. Forest Meteorol., 141, 19–34, https://doi.org/10.1016/j.agrformet.2006.09.002, 2006.

Powell, T. L., Galbraith, D. R., Christoffersen, B. O., Harper, A., Imbuzeiro, H. M. A., Rowland, L., Almeida, S., Brando, P. M., da Costa, A. C. L., Costa, M. H., Levine, N. M., Malhi, Y., Saleska, S. R., Sotta, E., Williams, M., Meir, P., and Moorcroft, P. R.: Confronting model predictions of carbon fluxes with measurements of Amazon forests subjected to experimental drought, New Phytol., 200, 350–365, https://doi.org/10.1111/nph.12390, 2013.

Priestley, C. H. B. and Taylor, R. J.: On the assessment of surface heat flux and evaporation using large-scale parameters, Mon. Weather Rev., 100, 81–92, https://doi.org/10.1175/1520-0493(1972)100<0081:otaosh>2.3.co;2, 1972.

Rigden, A. J. and Salvucci, G. D.: Evapotranspiration based on equilibrated relative humidity (ETRHEQ): Evaluation over the continental US, Water Resour. Res., 51, 2951–2973, https://doi.org/10.1002/2014wr016072, 2015.

Rind, D., Goldberg, R., Hansen, J., Rosenzweig, C., and Ruedy, R.: Potential evapotranspiration and the likelihood of future drought, J. Geophys. Res.-Atmos., 95, 9983–10004, https://doi.org/10.1029/JD095iD07p09983, 1990.

Roerink, G. J., Su, Z., and Menenti, M.: S-SEBI: a simple remote sensing algorithm to estimate the surface energy balance, Phys. Chem. Earth Pt. B, 25, 147–157, 2000.

Scheff, J. and Frierson, D. M. W.: Scaling potential evapotranspiration with greenhouse warming, J. Climate, 27, 1539–1558, https://doi.org/10.1175/JCLI-D-13-00233.1, 2013.

Schellekens, J., Dutra, E., Martinez-de la Torre, A., Balsamo, G., van Dijk, A., Weiland, F. S., Minvielle, M., Calvet, J. C., Decharme, B., Eisner, S., Fink, G., Florke, M., Pessenteiner, S., van Beek, R., Polcher, J., Beck, H., Orth, R., Calton, B., Burke, S., Dorigo, W., and Weedon, G. P.: A global water resources ensemble of hydrological models: the eartH2Observe Tier-1 dataset, Earth Syst. Sci. Data, 9, 389–413, https://doi.org/10.5194/essd-9-389-2017, 2017.

Seiller, G. and Anctil, F.: How do potential evapotranspiration formulas influence hydrological projections?, Hydrolog. Sci. J. – J. Sci. Hydrol., 61, 2249–2266, https://doi.org/10.1080/02626667.2015.1100302, 2016.

Sentelhas, P. C., Gillespie, T. J., and Santos, E. A.: Evaluation of FAO Penman–Monteith and alternative methods for estimating reference evapotranspiration with missing data in Southern Ontario, Canada, Agr. Water Manage., 97, 635–644, https://doi.org/10.1016/j.agwat.2009.12.001, 2010.

Sheffield, J., Wood, E. F., and Roderick, M. L.: Little change in global drought over the past 60 years, Nature, 491, 435–438, https://doi.org/10.1038/nature11575, 2012.

Shuttleworth, W. J.: Evaporation, in: Handbook of Hydrology, chap. 4, edited by: Maidment, D. R., McGraw-Hill, New York, USA, 4.1–4.53, 1993.

Shuttleworth, W. J. and Calder, I. R.: Has the Priestley-Taylor equation any relevance to forest evaporation?, J. Appl. Meteorol., 18, 639–646, https://doi.org/10.1175/1520-0450(1979)018<0639:htptea>2.0.co;2, 1979.

Singh, V. P. and Xu, C. Y.: Sensitivity of mass transfer-based evaporation equations to errors in daily and monthly input data, Hydrol. Process., 11, 1465–1473, https://doi.org/10.1002/(SICI)1099-1085(199709)11:11<1465::AID-HYP452>3.0.CO;2-X, 1997.

Slatyer, R. and McIlroy, I.: Practical micrometeorology, CSIRO, Canberra, Australia, 1961.

Song, Q. H., Braeckevelt, E., Zhang, Y. P., Sha, L. Q., Zhou, W. J., Liu, Y. T., Wu, C. S., Lu, Z. Y., and Klemm, O.: Evapotranspiration from a primary subtropical evergreen forest in Southwest China, Ecohydrology, 10, e1826, https://doi.org/10.1002/eco.1826, 2017.

Sumner, D. M. and Jacobs, J. M.: Utility of Penman–Monteith, Priestley–Taylor, reference evapotranspiration, and pan evaporation methods to estimate pasture evapotranspiration, J. Hydrol., 308, 81–104, https://doi.org/10.1016/j.jhydrol.2004.10.023, 2005.

Teuling, A. J. and Seneviratne, S. I.: Contrasting spectral changes limit albedo impact on land-atmosphere coupling during the 2003 European heat wave, Geophys. Res. Lett., 35, L03401, https://doi.org/10.1029/2007GL032778, 2008.

Thornthwaite, C. W.: An approach toward a rational classification of climate, Geogr. Rev., 38, 55–94, https://doi.org/10.2307/210739, 1948.

Vautard, R., Cattiaux, J., Yiou, P., Thepaut, J. N., and Ciais, P.: Northern Hemisphere atmospheric stilling partly attributed to an increase in surface roughness, Nat. Geosci., 3, 756–761, https://doi.org/10.1038/ngeo979, 2010.

Vicente-Serrano, S. M., Gouveia, C., Camarero, J. J., Beguería, S., Trigo, R., López-Moreno, J. I., Azorín-Molina, C., Pasho, E., Lorenzo-Lacruz, J., Revuelto, J., Morán-Tejeda, E., and Sanchez-Lorenzo, A.: Response of vegetation to drought time-scales across global land biomes, P. Natl. Acad. Sci. USA, 110, 52–57, https://doi.org/10.1073/pnas.1207068110, 2013.

Vicente-Serrano, S. M., Miralles, D. G., Domínguez-Castro F., Azorín-Molina, C., El Kenawy, A., McVicar, T. R., Tomás-Burguera, M., Beguería, S., Maneta, M., and Peña-Gallardo, M.: Global assessment of the Standardized Evapotranspiration Deficit Index (SEDI) for drought analysis and monitoring, J. Climate, 31, 5371–5393, https://doi.org/10.1175/jcli-d-17-0775.1, 2018.

Viswanadham, Y., Silva Filho, V. P., and André, R. G. B.: The Priestley–Taylor parameter α for the Amazon forest, Forest Ecol. Manage., 38, 211–225, https://doi.org/10.1016/0378-1127(91)90143-J, 1991.

Wang, K., Wang, P., Li, Z., Cribb, M., and Sparrow, M.: A simple method to estimate actual evapotranspiration from a combination of net radiation, vegetation index, and temperature, J. Geophys. Res.- Atmos., 112, D15107, https://doi.org/10.1029/2006JD008351, 2007.

Wang, K. C. and Dickinson, R. E.: A review of global terrestrial evapotranspiration: observation, modeling, climatology, and climatic variability, Rev. Geophys., 50, RG2005, https://doi.org/10.1029/2011rg000373, 2012.

Weiß, M. and Menzel, L.: A global comparison of four potential evapotranspiration equations and their relevance to stream flow modelling in semi-arid environments, Adv. Geosci., 18, 15–23, https://doi.org/10.5194/adgeo-18-15-2008, 2008.

Xu, C. Y. and Singh, V. P.: Evaluation and generalization of temperature-based methods for calculating evaporation, Hydrol. Process., 15, 305–319, https://doi.org/10.1002/hyp.119, 2001.

Xu, C. Y. and Singh, V. P.: Cross comparison of empirical equations for calculating potential evapotranspiration with data from Switzerland, Water Resour. Manage., 16, 197–219, https://doi.org/10.1023/A:1020282515975, 2002.

Yoder, R. E., Odhiambo, L. O., and Wright, W. C.: Evaluation of methods for estimating daily reference crop evapotranspiration at a site in the humid Southeast United States, Appl. Eng. Agricult., 21, 197–202, https://doi.org/10.13031/2013.18153, 2005.

Zhou, M. C., Ishidaira, H., Hapuarachchi, H. P., Magome, J., Kiem, A. S., and Takeuchi, K.: Estimating potential evapotranspiration using Shuttleworth–Wallace model and NOAA-AVHRR NDVI data to feed a distributed hydrological model over the Mekong River basin, J. Hydrol., 327, 151–173, https://doi.org/10.1016/j.jhydrol.2005.11.013, 2006.