Testing a maximum evaporation theory over saturated land: implications for potential evaporation estimation

. State-of-the-art evaporation models usually as-sume net radiation ( R n ) and surface temperature ( T s ; or near-surface air temperature) to be independent forcings on evaporation. However, R n depends directly on T s via outgoing longwave radiation, and this creates a physical coupling between R n and T s that extends to evaporation. In this study, we test a maximum evaporation theory originally developed for the global ocean over saturated land surfaces, which explicitly acknowledges the interactions between radiation, T s , and evaporation. Similar to the ocean surface, we ﬁnd that a maximum evaporation ( LE max ) emerges over saturated land that represents a generic trade-off between a lower R n and a higher evaporation fraction as T s increases. Compared with ﬂux site observations at the daily scale, we show that LE max corresponds well to observed evaporation under non-water-limited conditions and that the T s value at which LE max occurs also corresponds with the observed T s . Our results suggest that saturated land surfaces behave essentially the same as ocean surfaces at timescales longer than a day and further imply that the maximum evaporation concept is a natural attribute of saturated land surfaces, which can be the basis of a new

Abstract. State-of-the-art evaporation models usually assume net radiation (R n ) and surface temperature (T s ; or near-surface air temperature) to be independent forcings on evaporation. However, R n depends directly on T s via outgoing longwave radiation, and this creates a physical coupling between R n and T s that extends to evaporation. In this study, we test a maximum evaporation theory originally developed for the global ocean over saturated land surfaces, which explicitly acknowledges the interactions between radiation, T s , and evaporation. Similar to the ocean surface, we find that a maximum evaporation (LE max ) emerges over saturated land that represents a generic trade-off between a lower R n and a higher evaporation fraction as T s increases. Compared with flux site observations at the daily scale, we show that LE max corresponds well to observed evaporation under non-water-limited conditions and that the T s value at which LE max occurs also corresponds with the observed T s . Our results suggest that saturated land surfaces behave essentially the same as ocean surfaces at timescales longer than a day and further imply that the maximum evaporation concept is a natural attribute of saturated land surfaces, which can be the basis of a new approach to estimating evaporation.

Introduction
Potential evaporation (E P ), defined as the rate of evaporation (E) that would occur under non-water-stressed conditions, determines the upper boundary of E over a specific land surface for a given meteorological forcing. Although E P is more of a hypothetical variable and is generally very difficult to observe, it is often the starting point for partitioning rainfall between E, runoff, and soil moisture changes in hydrological, agricultural, ecological, and other related studies (Maes et al., 2019;Milly and Dunne, 2016;Scheff and Frierson, 2014;Schellekens et al., 2017;Sheffield et al., 2012;Vicente-Serrano et al., 2013;Wang and Dickinson, 2012). Over the years, numerous mathematical models have been proposed with varying structures and complexities to quantify E P (e.g., Allen et al., 1998;Priestley and Taylor, 1972;Penman, 1948;Shuttleworth, 1993;Thornthwaite, 1948). Among them, the Penman-Monteith type models (e.g., either the Open Water Penman model; Shuttleworth, 1993; or the Food and Agriculture Organization Penman-Monteith model; Allen et al., 1998) are most widely used, given their explicit consideration of the radiative and aerodynamic components of E, and are hence generally considered as a physical-based and accurate approximation of the real E processes.
Nevertheless, recent empirical evidence shows that the Penman-Monteith type models perform unsatisfactorily in estimating E P compared with eddy covariance observations (i.e., the observed E under non-water-stressed conditions; Maes et al., 2019). Instead, the energy balance-based approaches work better in reproducing E P in both observations (Maes et al., 2019) and climate model simulations (Milly and Dunne, 2016). From an energy balance point of view, the magnitude of E (or in its energy form -latent heat flux or LE) is determined by the energy balance equation: with R n the net radiation (W m −2 ) and G the ground heat flux (W m −2 ), which is often negligibly small over land for timescales longer than a day. In Eq. (1), β is the Bowen ratio and represents the ratio of sensible heat flux (H ) over LE (Bowen, 1926). As a result, LE is determined by the available energy at an evaporating surface (i.e., R n − G) and the ability of that evaporating surface to convert the available energy into LE, which is represented by the 1/(1 + β) term and is often known as the evaporative fraction. With no restriction on water supply, β is known to be a decreasing function of temperature at an evaporating surface (T s ; Aminzadeh et al., 2016;Andreas et al., 2013;Guo et al., 2015;Philip, 1987;Priestley and Taylor, 1972;Slatyer and McIlroy, 1961;. This implies that when water is not limiting, both T s and the available energy determine the rate of E. Hence, with fixed available energy, a higher T s corresponds to a lower β (or a higher evaporative fraction) and therefore a larger LE. This line of reasoning has directly led to the development of energy balancebased evaporation models, including the classic equilibrium evaporation approach (Slatyer and McIlroy, 1961) and the Priestley-Taylor evaporation model (Priestley and Taylor, 1972). Compared with Penman-Monteith type models, the energy balance-based approach simplifies the representation of the aerodynamic component of E and usually takes the aerodynamic component of E as a fixed fraction of its radiative counterpart (e.g., 0.26 in the Priestley-Taylor model).
However, a key issue in the above energy balance-based approach is that it takes R n to be an independent forcing of E. A similar idea was also adopted in Penman-Monteith type models (Penman, 1948;Monteith, 1965). Nevertheless, it is clear that R n cannot be physically independent of either E or T s . On one hand, a higher T s corresponds to higher outgoing longwave radiation and therefore a lower R n . On the other hand, a higher E is associated with larger evaporative cooling, which lowers T s and ultimately feed backs to R n . This latter process confirms that T s is not independent of E. Consequently, the intrinsic interdependence between R n , E, and T s has long been ignored in the state-of-the-art evaporation models that require R n as model input .
To deal with the above issue, a recent study by  explicitly considered the interdependence between radiation, T s , and evaporation and tested the new approach over global ocean surfaces. They found that with the increase of T s , R n decreases while evaporative fraction increases (since β decreases as T s increases), in agreement with a number of previous studies (Aminzadeh et al., 2016;Andreas et al., 2013;Guo et al., 2015;Philip, 1987;Priestley and Taylor, 1972;Slatyer and McIlroy, 1961). This generic and explicit trade-off between a lower R n and a higher evaporative fraction with the increase of T s directly yields a maximum evaporation along the T s gradient according to Eq. (1) ; also see Sect. 2.2). This maximum evaporation emerges naturally from the R n -T s -E inter-actions and does not require a priori knowledge of T s thereby alleviating the need for the assumption that R n and T s are independent of E in traditional evaporation models. As a result, the maximum evaporation theory does not consider R n to be an independent forcing of E. Instead, it only requires the incoming and reflected solar radiation and the assumption that β decreases with the increase of T s (see Sect. 2.2). Compared with observations of ocean surface evaporation and temperature,  demonstrated the validity of the maximum evaporation theory over global ocean surfaces. Here, we test this new maximum evaporation theory over land by asking and answering two questions: (i) does the theory recover the observed E and (ii) does it recover the observed T s ? By recovering, we mean that the maximum E as per theory corresponds to the observed E and that the T s at which the maximum E occurs corresponds to the observed T s under non-water-stressed conditions. Testing the maximum evaporation theory over land is important, as vegetation transpiration generally dominates the total evaporative flux over land (Jasechko et al., 2013;Lian et al., 2018), which is essentially different from ocean surfaces where the evaporative flux only consists of evaporation from open water surfaces. In addition, land surfaces usually have a larger surface roughness than ocean surfaces, which may result in different energy partitioning (into sensible heat and latent heat) between the ocean and the land. Therefore, it is crucial to test the maximum evaporation theory over land to determine whether saturated land behaves like the ocean surface and whether the maximum evaporation theory can be the basis of a new approach to estimating E P over land.

Flux site observations
Observations of actual daily evaporation (or latent heat flux), sensible heat flux, and ground heat flux along with relevant meteorological variables, radiative fluxes and soil moisture were originally obtained from 212 flux sites collected in the FLUXNET2015 database (http://fluxnet.fluxdata.org/ data/fluxnet2015-dataset/, last access: 2 November 2020). Only days with the data quality metric for LE and H higher than 0.9 (on a scale of 0-1) were used. The daily scale variables were obtained based on 15 min/30 min observations using the standard approach (Pastorello et al., 2020). The residual approach (i.e., assuming the observed H is correct and considering LE as the residual of the energy balance equation) was used to recalculate the fluxes based on a forced energy balance closure at each flux site (Ershadi et al., 2014). We also used the Bowen ratio approach (Twine et al., 2000) to force the flux-site energy balance closure and this resulted in similar model performance (Fig. S1 in the Supplement). The surface temperature for each site-day combination was calculated based on the observed longwave radiation following: where R lo and R li are, respectively, the outgoing and incoming longwave radiation, σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W m −2 K −4 ), and ε is the surface emissivity, which is acquired from the MODIS (Moderate Resolution Imaging Spectroradiometer) emissivity product (i.e., MOD11A1 Version 6; https://lpdaac.usgs.gov/products/ mod11a1v006, last access: 2 January 2021). The MOD11A1 surface emissivity has a daily temporal resolution and a 1 km spatial resolution. To obtain the emissivity for each eddy covariance (EC) flux site, we center on the pixel where the site is located and take the mean value of the 81 neighboring pixels (9 × 9 pixels) as the emissivity value of the site. For conditions when the MOD11A1 emissivity are not available, we deleted these site-days.
To select a subset of observations at each flux site in which the actual evaporation is not limited by water availability, the energy balance criterion and the soil moisture criterion used by Maes et al. (2019) were adopted. Specifically, at each flux site, the evaporative fraction (EF; i.e., EF = LE/(LE + H )) was first calculated and the unstressed measurements consisted of all days with EF exceeding the 95th percentile EF threshold at each site. Following that, we removed days with soil moisture (averaged over all measured depths) lower than 50 % of the maximum soil moisture (taken to be soil moisture at the 98th percentile) at each site. In addition, any remaining site-days with daily EF values lower than 0.6 were also removed. Finally, we removed days having a negative H value (accounting for ∼ 5 % of the total daily data) to avoid dealing with strongly advective conditions when accurate measurements are not guaranteed (Paw U et al., 2000;Wilson et al., 2002). As a result, a total of 1128 non-waterstressed site-days from 86 sites passed the above criteria and were used in this study ( Fig. 1 and Table S1 in the Supplement).

Overview of the maximum evaporation model
The maximum evaporation model calculates evaporation from a wet surface based essentially on the surface energy balance (Eq. 1), with R n and β both explicitly represented as functions of T s : In the above equation, the first term on the right-hand side (i.e., 1/[1 + β(T s )]) is the evaporative fraction, which is the ratio of the latent heat flux over the total available energy. Over wet surfaces, since the Bowen ratio decreases with T s (Aminzadeh et al., 2016;Andreas et al., 2013;Guo et al., 2015;Philip, 1987;Priestley and Taylor, 1972;Slatyer and McIlroy, 1961;, evaporative fraction increases with T s . On the other hand, the second term on the right-hand side of Eq. (3) is the total available energy, which decreases with the increase of T s as a higher T s directly leads to a higher outgoing longwave radiation and hence a lower R n . As a result, the trade-off between a higher evaporative fraction and a lower R n with the increase of T s would naturally lead to a maximum LE along the T s gradient according to Eq. (3). A previous study by  demonstrated that this naturally emergent maximum LE corresponds well to the actual LE over global ocean surfaces and the T s at which the maximum LE occurs also corresponds to the observed sea surface temperature. Here we will test whether this maximum evaporation approach is also valid over land under non-water-stressed conditions.

Parameterization of R n and β as a function of T s
To explicitly acknowledge the dependence of R n on T s , R n (T s ) is expressed as where R sn is the net shortwave radiation (W m −2 ) and is taken to be unchanged with T s . T is the temperature dif-ference between T s and the effective radiating temperature of the atmosphere (T rad ; assuming blackbody radiation, T rad = 4 √ R li /σ ), and is parameterized as a function of atmospheric transmissivity and geographic latitude : where τ is the atmospheric transmissivity for shortwave radiation (dimensionless) and is calculated as the ratio of incoming shortwave radiation at the Earth's surface to that at the top of the atmosphere. The parameter "lat" is the geographic latitude (in decimal degrees), which is considered here to account for a longer pathway of shortwave radiation going through the atmosphere in higher latitudes compared to lower latitudes. n 1 , n 2 , and n 3 are the fitting coefficients. Using extensive data over the global ocean (n = 202 794), Yang and Roderick (2019) determined the values of these coefficients to be n 1 = 2.52, n 2 = 2.38, and n 3 = 0.035, respectively. Here, we directly adopt these same coefficient values over land for two reasons: (i) the key processes governing the interactions between incoming and outgoing longwave radiation are essentially the same for ocean and land (mainly greenhouse gases that affect the vertical temperature structure of the atmosphere, and water vapor and aerosols that affect the formation of clouds), and (ii) there were many more samples available for parameterizing Eq. (5) over the ocean than that over land. Validation against observations from all 1128 non-water-limited site-days demonstrates an overall good performance of Eq. (5) in estimating T over land under saturated conditions (Fig. S2). The Bowen ratio (β) is expressed as a function of T s : where m is a fitting coefficient. γ is the psychrometric constant (kPa K −1 ), and is the slope of the saturation vapor pressure curve (kPa K −1 ), both of which are functions of T s : where C P is the specific heat of air at constant pressure (1.01 kJ kg −1 K −1 ), P a is the air pressure (kPa), and e s is the saturated vapor pressure (kPa). L is the latent heat of vaporization (kJ kg −1 ) and is calculated as a weak function of temperature: L (T s ) = 2.51 × 10 3 − 2.32 × (T s − 273.15) .
To apply the maximum evaporation model, an array of T s (e.g., from 250 to 330 K at an interval of 0.1 K) is generated along with the observed R sn and G; these are applied to Eqs. (4) and (6) and then Eq.
(3) to estimate LE at each corresponding T s . The maximum evaporation is then located in array as well as the surface temperature at which this maximum occurs (see Fig. 3 for an example).

Results
The maximum evaporation theory is tested at 86 flux sites globally, covering a wide range of bio-climates ( Fig. 1 and Table S1). By pooling daily observations of H , LE, and T s across all 1128 site-days, we first obtain a generic β-T s relationship as β = 0.27γ / . Similarly, we also obtained a β-T s relationship for each separate biome type as shown in Fig. 2. By comparison, Yang and Roderick (2019) reported a β-T s relationship over ocean of β = 0.24γ / . This means that for the same T s , β over land is generally larger than that over ocean. Interestingly, the ocean surface β-T s relationship is identical to that in wetlands obtained here. These β-T s relationships will be used in the following calculations of LE using the maximum evaporation approach.
To get an overview of how each of the energy fluxes varies with T s , we first examine the maximum evaporation theory using the pooled data over all 1128 site-days (Fig. 3). Under this condition, the mean observed net shortwave radiation (R sn ) over all site-days is 176.6 W m −2 and G is 1.0 W m −2 . Since R sn is not directly dependent on T s and G is negligibly small, the term R sn minus G is held constant across the entire T s range. With the increase of T s , it is readily apparent that both outgoing and incoming longwave radiation (R lo and R li ) steadily increase (see Sect. 2.2 for details about the coupling between R lo and R li ), with R lo increasing slightly faster than R li , leading to a decreased net longwave radiation and thus a decreased R n as T s increases (Fig. 3). With this and the observed generic dependence of β on T s (β = 0.27γ / , Fig. 2), a maximum LE emerges along the T s gradient that represents the interaction between decreasing R n and increasing evaporative fraction as T s increases. For the pooled dataset used here, the maximum LE (LE max ) is found to be 105.6 W m −2 and the corresponding T s is 294.7 K, both of which are very close to the averages computed from all daily flux site observations (i.e., LE obs = 102.4 W m −2 and T s_obs = 292.3 K) (Fig. 3).
Having demonstrated the overall concept, we next perform the detailed calculations using data for all individual site-days (Figs. 4-6). Using the same generic β dependence on T s (β = 0.27γ / ), the LE max estimated from the maximum evaporation model agrees very well with flux site observations, yielding an R 2 of 0.92, a root-mean-squared error (RMSE) of 14.6 W m −2 , and a mean bias of 1.6 W m −2 (Fig. 4a). The performance of the maximum evaporation model improves slightly when the biome-specific model parameters are used (RMSE decreases to 14.1 W m −2 and mean bias decreases to 1.4 W m −2 ; Fig. 4b). This result demonstrates that LE max corresponds to the observed evaporation under well-watered conditions across a broad range of bio- Figure 2. Relationship between Bowen ratio (β) and surface temperature (T s ) over saturated land surfaces. The thick black curve represents the fitted β-T s relationship across all data points (i.e., n = 1128, β = 0.27γ / , R 2 = 0.11, p < 0.001), and the colored lines represent different biome types with the number of data points (n site-days) and fitted β-T s relationship for each biome type shown in the legend. climates. In fact, when the previously identified ocean surface β-T s relationship is adopted, the maximum evaporation approach performs only slightly worse than those based on the calibrated β-T s relationship over saturated lands, yielding an R 2 of 0.91, an RMSE of 14.8 W m −2 , and a mean bias of 2.8 W m −2 (Fig. 4c).
We next test whether the maximum evaporation approach could recover T s over the same saturated land surfaces. Similar to the test of LE, the three β-T s relationships are respectively used. Results show that when the generic β-T s relationship over land is used (i.e., β = 0.27γ / ), the T s at which LE max occurs corresponds reasonably well to the observed T s , with an R 2 of 0.62, an RMSE of 4.3 K, and a mean bias of 0.3 K (Fig. 5a), indicating that the maximum evaporation approach is also able to recover T s under saturated conditions. Again, the model's performance in recovering T s increases slightly when the biome-specific β-T s relationships are used (Fig. 5b). When the ocean surface β-T s relationship is used, the model performs similarly in estimating the variability of T s to that of the generic land β-T s relationship (Fig. 5c). However, the ocean surface β-T s relationship ( Fig. 5c) results in a higher T s mean bias compared to the β-T s relationships obtained over land (Fig. 5a).
Different from most state-of-the-art evaporation models, the maximum evaporation approach does not rely on observed R n (or independent R n estimates) as model input but estimates R n as a result of the R n -T s -E interaction. Here, we also test the estimated R n calculated using the maximum evaporation approach as the discrepancy between LE max and LE obs is mainly caused by the slight difference between T s_max and T s_obs that leads to different R lo and R li (and thus a different R n ) being used in the calculation of LE max . It should be noted that since the observed shortwave radiation is used in the maximum evaporation model, validation of R n is essentially the same as the validation of net longwave radiation. We find that the maximum evaporation model could satisfactorily reproduce the observed R n when the generic land β-T s relationship is used, as indicated by an R 2 of 0.93, an RMSE of 14.4 W m −2 , and a mean bias of 2.3 W m −2 (Fig. 6a). Using biome-specific β-T s relationships or the ocean surface β-T s relationship does not considerably increase or decrease the model's performance in estimating R n (Fig. 6b and c).

Discussion
Taking R n and/or T s (or near-surface air temperature) to be independent forcings has long been identified as a scientific concern in the use of evaporation models (Milly, 1991;Monteith and Unsworth, 2013;Philip, 1987). Here, we test a maximum evaporation theory developed over the global ocean surface that addresses this concern by explicitly acknowledging the interdependence between radiation, surface temperature, and evaporation . Our new results show that there exists a maximum evaporation along the T s gradient that corresponds to the observed evaporation under saturated conditions over land (Figs. 3 and 4). In addition, the T s at which LE max occurs also corresponds reasonably well to the observed T s (Figs. 3 and 5). These results mirror those found previously over the global ocean (Yang    Fig. 2). (c) Ocean surface β-T s relationship (β = 0.24γ / , . The colors indicate different biome types (as provided in Fig. 1). The dashed black line indicates the 1 : 1 line. and . This is not a surprise since the basic principles are the same for a wet land surface and the ocean surface. These results suggest that saturated land surfaces behave essentially the same as ocean surfaces and imply that LE max is a natural attribute of the land surface when water availability does not limit evaporation.
A key assumption involved in the maximum evaporation model is that β decreases with the increase of T s under saturated conditions. Nevertheless, this key assumption that β decreases with the increase of T s under saturated conditions has been extensively validated in previous studies based on theoretical relationships (Philip, 1987;Priestley and Taylor, 1972;Slatyer and McIlroy, 1961;Lhomme, 1997) and in situ observations (Andreas et al., 2013;Guo et al., 2015;; also see Fig. S3). Moreover, our results also found this held over saturated lands (Fig. 2). The original maximum evaporation study reported that β = 0.24γ / over global ocean surfaces . Here, we find the generic land surface coefficient increases to 0.27 (i.e., β = 0.27γ / , Fig. 2), which indicates a slightly higher β over wet vegetated land than that over the ocean surface for the same T s . This is biophysically reasonable, as the stoma of plant leaves represents an additional resistance to vapor transfer between the land and the atmosphere (Swann et al., 2016;, which lowers the ability of a generic vegetated surface to convert available energy into LE for a given T s , compared to open water surfaces. In addition, different surface roughness can also lead to different β-T s relationships between the land and the ocean. Compared with the ocean surface that shows a tight β-T s relationship , the β-T s relationships over saturated vegetated land are relatively weak with considerable scatter (Fig. 2). This data scatter could be caused by several reasons. First, the observations by EC towers can be a source of uncertainty. This is threefold, including (i) the quality of the observations, (ii) the footprint within each EC tower (possibly) being heterogeneous (Lee et al., 2004;Paw et al., 2000), and (iii) whether or not the selected days are truly non-water-limited (however, see Fig. S4). Second, as is seen in Fig. 2, different biome types exhibit different β-T s relationships. This can be caused by different surface resistance and roughness between biome types and even between sites. Nevertheless, these data-based limitations only have limited impacts on the model performance, as similar performance is obtained using both the generic β-T s relationship (i.e., β = 0.27γ / ) and biome-specific β-T s relationships (Fig. 4). Third, wind speed could be another factor that leads to the scatter. For the same surface roughness, a different wind speed will lead to a different aerodynamic resistance and therefore a different β. However, this effect is usually very small, as demonstrated by the long-standing similarity theory (the transfer of mass and heat share the same aerodynamic process in the lower atmospheric boundary layer; Monin and Obukhov, 1954). In fact, our findings that one can make a reasonable estimate of LE using a generic land or ocean β-T s relationship instead of a site-specific relationship (Fig. 4) imply that R n is the primary determinant of LE over saturated surfaces. As evaporation tends to operate at its maximum strength, sensible heat (and β) are usually very small over warm saturated land surfaces. As a result, once R n can be accurately determined, any reasonable β-T s relationship (Fig. 2) would result in a satisfactory LE estimate (Figs. 4 and S5). Our result highlighted in Fig. 3 shows that R n (and hence LE) is only a weak function of T s and this explains why one can obtain an accurate estimate of LE using a generic β-T s relation. However, the same logic also leads to the conclusion that an accurate β-T s relationship will be necessary to estimate T s , since T s is very sensitive to changes in LE (Fig. 3). In this regard, using the land β-T s relationships (preferably site-specific relations) is preferable to a generic ocean surface relation (Fig. 5). To further demonstrate the above points, we conduct an uncertainty test by varying the coefficient m in the β-T s relationship in Fig. S6. We find that when m ranges from 0.18 to 0.36 (all other forcings as per Fig. 3), the change in estimated LE max is only 9 W m −2 (101.7-110.7 W m −2 ), whereas the change in estimated T s_max is as high as 11.6 K (287.9-299.5 K).
Besides the data scattering that leads to an uncertainty in the β-T s relationships, there are also uncertainties associated with (i) parameterization of the longwave coupling and (ii) selection of non-water-stressed observations in the current study. In the maximum evaporation approach, the coupling between outgoing and incoming longwave radiation is calculated using the temperature difference between the surface and an effective radiating height in the atmosphere ( T ) and is parameterized as a function of shortwave atmospheric transmissivity and geographic latitude. However, the shortwave atmospheric transmissivity is primarily affected by aerosols, while the longwave transmissivity is mainly affected by the concentration of greenhouse gases. Nevertheless, here we only deal with wet conditions, under which the vapor concentration of the atmosphere is also relatively high and more aerosols would favor the development of more clouds that simultaneously affect both shortwave and longwave radiation. We suspect that this underlies the excellent performance of Eq. (5) in estimating T at the flux sites (Fig. S2). To further evaluate that conclusion, we additionally evaluate the estimated longwave radiation against four global products (i.e., ERA5, Hersbach et al., 2020;CERES, Kato et al., 2018; the Princeton global forcing data, Sheffield et al., 2006; the GLDAS global forcing data, Rodell et al., 2004) and compare our longwave estimates with other two semi-empirical models (i.e., Brutsaert, 1975;Shakespeare and Roderick, 2021). The results show our T -based approach to be the best performer across a wide range of conditions when the surface is wet (Fig. S7). In addition, we further note that our maximum evaporation model is only tested at the daily timescale  and longer (Fig. 3). In particular, for timescales shorter than that (e.g., hourly), the diurnal cycle of E can be very different for ocean and land surfaces (Kleidon and Renner, 2017). In addition, the parameterization of the coupling between incoming and outgoing longwave radiation via Eq. (5) requires a timescale that is long enough to allow the surface heat fluxes to be fully redistributed through the atmospheric column . At sub-daily scales, Eq. (5) is likely invalid because R lo usually exhibits a larger diurnal range than R li during a typical cloudless day (Monteith and Unsworth, 2013). For even longer periods, especially for assessing the impacts of climate change, the relationship between shortwave and longwave radiation used herein may also be invalid, as we expect this relationship to evolve with anthropogenic climate change.
As for the selection of non-water-stressed evaporation observations from global EC towers, we rely largely on the same selection criteria used in a previous study (Maes et al., 2019). However, these selection criteria are somewhat subjective and represent a compromise between better data quality and more data samples. As a result, the selected site-days are not necessarily non-water-limited. Nevertheless, varying the selection criteria (changing the thresholds) of non-waterstressed evaporation only resulted in minor changes in the overall model performance (Fig. S4), which suggests that the uncertainties in the selection of non-water-stressed evaporation observations would not materially change our conclusion.
The ability of the maximum evaporation model to recover LE and T s over vegetated lands under saturated conditions has an important implication for the estimation of potential evaporation, which is a central concept in hydrology and agriculture (and especially in irrigation). The underlying idea of E P is straightforward -it is the evaporation that would occur with an unlimited supply of water. However, the formal physical definition of E P has been widely debated in the literature (Brutsaert, 2005;Donohue et al., 2010;Granger, 1989;Nash, 1989) and the calculation of E P using conventional evaporation models is problematic (Aminzadeh et al., 2016;Roderick et al., 2015). The key scientific issue is that the meteorological forcing variables observed over actual surfaces are generally not equivalent to the meteorological variables that would be measured over a hypothetical surface with an unlimited water supply. Compared with existing evaporation models, the maximum evaporation model presented here requires fewer meteorological variables than existing approaches (but performs similarly with existing approaches under wet conditions; see Fig. S8 for details). This new approach only requires the incoming and reflected solar radiation, a relationship that describes a decreasing dependence of β on T s , and a relation for the coupling of the incoming and outgoing longwave radiation. With these modest requirements, LE max naturally follows from the physical interdependence between radiation, surface temperature, and evaporation. These features suggest that the maximum evaporation model can be used to make a strictly independent estimate of E P . In fact, the maximum evaporation formulation directly maps to one particular definition of E P that was proposed by Brutsaert (2015) as "the maximum evaporation that would occur over real surfaces with the actual solar forcing and a prescribed Bowen ratio".

Conclusions
In this study, we test a maximum evaporation theory that explicitly acknowledges the interdependence between radiation, surface temperature, and evaporation over saturated land surfaces. Validated against flux site observations, we show that the maximum evaporation approach could recover the observed evaporation across a broad range of bioclimates. In comparison, although the model is also able to reasonably recover the observed T s , the model's performance in recovering T s is not as good as that for LE. Nevertheless, this does not materially lead to larger errors in LE estimates, as we additionally demonstrate that LE is not sensitive to T s changes. The overall good performance of the maximum evaporation approach over saturated surfaces implies the method's great potential for use in estimating potential evaporation. To calculate E P in practice using the maximum evaporation approach, a detailed site-specific (or biome-specific) β-T s relationship (e.g., Fig. 2) would be favorable; otherwise, a generic default β-T s relationship (β = 0.27γ / or even β = 0.24γ / ) can also lead to a reasonable E P estimate that remains consistent with the E P definition by Brutsaert (2015). Table S2 gives a working example of applying the maximum evaporation model for E P estimation.
Data availability. All data for this paper are properly cited and referred to in the reference list. The FLUXNET2015 dataset is acquired from http://fluxnet.fluxdata.org/data/fluxnet2015-dataset/ (FLUXNET, 2016).
Author contributions. YY and MLR conceived the idea and designed the research. ZT performed the calculation. YY drafted the manuscript. All authors contributed to the results, discussion, and manuscript writing.
Competing interests. The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Acknowledgements. Michael L. Roderick acknowledges the support of the Australian Research Council (DP190100791). The FLUXNET community is greatly appreciated for making the eddy covariance data publicly available.
Financial support. This study is financially supported by the National Natural Science Foundation of China (grant nos. 42041004, 42071029, 41890821) and the Ministry of Science and Technology of China (grant no. 2019YFC1510604).
Review statement. This paper was edited by Adriaan J. (Ryan) Teuling and reviewed by Miriam Coenders-Gerrits and two anonymous referees.