Practical and Theoretical Benefits of an Alternative to the Penman‐Monteith Evapotranspiration Equation

The Penman‐Monteith equation is used widely to estimate evapotranspiration (E) and to understand its governing physics. I present an alternative to the Penman‐Monteith equation that has both practical and theoretical advantages, at no appreciable cost. In particular, the new equation requires no additional assumptions, empiricism, or computational cost compared with the Penman‐Monteith equation. Practically, the new equation is consistently more accurate over a wide range of conditions when compared with eddy covariance observations: The new equation has lower errors compared with Penman‐Monteith estimates of ET at all of the 79 eddy covariance sites available for the analysis. Using the new equation reduces errors, on average, by 67%, from 8.55 to 2.81 [W m−2]. At night, the improvement is even greater (92% reduction in error; from 1.26 to 0.097 [W m−2]). This improvement is achieved without calibration. Theoretically, the new equation corrects a conceptual error in the Penman‐Monteith equation, in which the Penman‐Monteith equation incorrectly implies that E from a saturated surface into a saturated, turbulent atmosphere (“equilibrium” E) is exactly equivalent to E from an unsaturated surface into an unsaturated, laminar atmosphere. The conceptual error is traced back to the failure of the Penman‐Monteith equation in important limiting cases; these errors are eliminated by the new equation. I use the new equation to revise an existing theory of land‐atmosphere coupling affected by the conceptual error in the Penman‐Monteith equation and to reassess several common but incorrect definitions of equilibrium E.


Introduction
Evapotranspiration (E [kg m −2 s −1 ]) is the second largest flux of water, after precipitation, in the terrestrial water cycle. It is a key component of the surface energy budget, in which the difference between net radiation and ground heat flux (R n −G [W m −2 ]) balances turbulent fluxes of sensible and latent heat (H+ E [W m −2 ], where is the latent heat of vaporization [J kg −1 ]). During the daytime, when typically R n − G > 0, E cools and moistens the lower atmosphere. E can also remain significant at night (Dawson et al., 2007;Groh et al., 2019;Novick et al., 2009), when typically R n − G < 0. E is governed by the surface energy budget, where T s is surface temperature [K], T a is air temperature at screen-level [K], q * (T s ) is saturated specific humidity at the surface [-] specified by the Clausius-Clapeyron relation, q a is specific humidity of air at screen level [-]), is air density [kg m −3 ], and c p is the specific heat capacity of air at constant pressure [J kg −1 K −1 ]. E is constrained by water availability and plant physiology, modeled by the bulk parameter g s [m s −1 ], often referred to as the "ecosystem conductance." It is also constrained by turbulent transport, modeled by the bulk parameter g a [m s −1 ], or "aerodynamic conductance." Equation (1) gives the "radiatively uncoupled" surface energy budget (Raupach, 2001), in which R n and G are treated as known and do not vary explicitly with surface temperature T s ; this applies, for instance, if direct observations of R n and G are available. If R n and G are modeled rather than observed, their dependence on T s can be retained in the analysis by introducing a "radiative conductance," g r , and "storage conductance," g g (the "radiatively coupled" case, considered in Appendix B).  (8): (a) approximations to the Clausius-Clapeyron relation, where q * (T s ) is the saturated specific humidity [-] and T s is surface temperature [K]. In this example, the vertical gray line is the air temperature T a [K] used in each approximation. The linear approximation, used in the PM equation, is a first-order Taylor expansion around T s = T a . The quadratic and quartic approximations are second-and fourth-order expansions, respectively. Comparison of equation (8) (b, c) and the PM equation (d, e) with half-hourly observations from 79 FLUXNET sites. The red dashed line is the 1:1 line. "PM" means "Penman-Monteith." "RMSE" means "root-mean-squared error." Given observations or model estimates of R n , G, q a , T a , g s and g a , equation (1) is an implicit equation for T s . Since T s is present in the definition of E, it is also an implicit equation for E. Equation (1) cannot be solved for T s analytically, due to the presence of the nonlinear q * (T s ). An approximate solution can be obtained by linearizing q * (T s ) around T = T a (Figure 1), that is, Using (i) the definition of H to replace T s − T a with H∕( c p g a ), (ii) equation (1) to replace H with R n − G − E, (iii) substituting the resulting equation into the definition of E, and (iv) solving for E yields the explicit equation where = Δ c p [-]. This is the radiatively uncoupled Penman-Monteith (PM) equation (Monteith, 1965;Penman, 1948) (see Appendix B for the radiatively coupled version). The major benefit of the PM equation is its removal of explicit dependence on T s . While the surface energy balance could be solved numerically for T s and E, the PM equation provides an approximate but explicit solution that has proven theoretically useful in understanding the coupled land-atmosphere system (e.g., Jarvis & McNaughton, 1986;Konings et al., 2011;Raupach, 2001;Scheff & Frierson, 2013;Stap et al., 2014;van Heerwaarden et al., 2010). It has also been widely applied in estimating E from field data, used within macroscale hydrology models (Hamman et al., 2018;Liang et al., 1994) and land surface models (Kumar et al., 2017) and has been inverted to estimate ecosystem-scale g s from observations of E and other variables (e.g., Baldocchi et al., 1991;Lin et al., 2018;Novick et al., 2016).

MCCOLL 2 of 15
Unfortunately, the linearization of the Clausius-Clapeyron relation in the PM equation introduces both empirical and conceptual errors into estimates of E. Empirically, the linearization of the Clausius-Clapeyron relation can be quite inaccurate (Milly, 1991;Paw U & Gao, 1988), particularly at night and in cold environments. Conceptually, the PM equation is incorrect in important limiting cases, which has led to common conceptual mistakes in the literature, related to E under well-watered conditions and other limiting cases (Paw U & Gao, 1988).
To remedy these problems, I propose an alternative to the PM equation-based on an approximation of the Clausius-Clapeyron relation proposed by Vallis et al. (2019)-that substantially reduces its empirical errors and eliminates its conceptual errors. Like the PM equation, the new equation is an explicit solution for E based on equation (1). Unlike the PM equation, it reproduces the correct limiting behavior, eliminating conceptual errors. Empirically, the new equation is more accurate than the PM equation across a wide range of real-world conditions, at no additional cost in terms of required assumptions, inputs, parameters, empiricism, or computation. Many challenges in modeling ET are not addressed by this work, including the modeling and estimation of g a and g s ; rather, the focus is on improving the PM equation without introducing additional costs.
This manuscript is structured as follows. In section 2, an alternative to the PM equation is presented. In section 2.1, the alternative equation is shown to outperform the PM equation in a range of real-world cases using eddy covariance observations. Readers who are primarily interested in seeing empirical evidence of the accuracy of the new equation should feel free to skip to this section. In section 2.2, the theoretical benefits of equation (8) are introduced. One benefit is its accuracy in limiting cases, which contrasts with the PM equation; the behavior of both equations in various limiting cases is considered in section 2.3. In section 2.4, the new equation is used to revise an existing theory of land-atmosphere coupling. In light of results in the previous sections, five different definitions of equilibrium E are critically reassessed in section 2.5. Conclusions are presented in section 3. For readers who are interested in immediately using the more accurate radiatively uncoupled equation for E, it is given in equation (8). Code for applying the equation is available from the author's website.

An Alternative to the Penman-Monteith Equation
In this section, I introduce the alternative to the PM equation. The new equation is analogous to the PM equation: It uses an approximation of the Clausius-Clapeyron relation to provide an equation for E that does not have any explicit dependence on T s . However, compared with the PM equation, I use a much less severe approximation of the Clausius-Clapeyron relation, resulting in a more accurate expression for E at no additional cost: no extra assumptions, parameters, inputs, or computational expense are required. Several previous studies have used higher-order Taylor expansions of the Clausius-Clapeyron relation ( Figure 1a) to derive quadratic and quartic versions of the PM equation that are also more accurate (Baldocchi et al., 2005;Milly, 1991;Paw U & Gao, 1988). While useful in daytime conditions (R n − G > 0), the quadratic and quartic polynomial equations on which these solutions are constructed can be undefined in some cases when R n − G < 0. This work focuses exclusively on general solutions of equation (1), which are defined for all values of R n − G, like the PM equation itself. For notational simplicity, throughout this manuscript, I use the term "daytime" to mean "R n − G > 0" and "nighttime" to mean "R n − G < 0," while recognizing that this is a simplification.
The Clausius-Clapeyron relation can be written as where R v is the gas constant for water vapor [J kg −1 K −1 ]. The relation can be written in terms of q * rather than saturated vapor pressure e * [Pa] because these terms are proportional to a very good approximation in Earth's lower atmosphere. While is a function of T, its dependence on T is relatively weak under standard atmospheric conditions, and it is conventionally treated as constant. Applying this approximation and integrating equation (4) between the screen-level air temperature T a and the surface temperature T s gives

10.1029/2020WR027106
Assuming T s − T a is small, as assumed in the PM equation, Vallis et al. (2019) proposed that this could be approximated as ) .
(6) Figure 1a compares this approximation to the approximations made in the linear PM equation (equation (2)) and two higher-order approximations. While all approximations are reasonable when |T s − T a | is small, equation (6) is much more robust to increases in |T s −T a |. If the exponential term in equation (6) is linearized around T a (using the relation exp(x) ≈ 1 + x near x = 0), and the Clausius-Clapeyron relation (equation (4)) is substituted in, the resulting expression is identical to equation (2). Therefore, the approximation in equation (6) is already assumed implicitly in the derivation of the PM equation and is not a new assumption.
Substituting equation (6) into the radiatively uncoupled surface energy balance (equation (1)) gives This equation can be solved exactly for T s ; see Appendix A for a derivation of the solution. Substituting this solution into equation (6) and then into the definition of E, an expression for E is obtained that has no explicit dependence on T s , analogous to the PM equation: Equation (8) uses the principal branch of the Lambert W function W(x), an analytic function defined by W(x) is multivalued for −1 e ≤ x < 0, with a principal branch W 0 (x) and a negative branch W −1 (x), and single valued for x ≥ 0. Since x ≥ 0 in equation (8), it is safe to restrict our attention to W 0 (x) ( Figure  S1). While less common than other analytic functions (e.g., the natural logarithm, which is defined similarly as log(exp(x)) = x), the Lambert W function has been used for centuries. (Corless et al., 1996). Applications of the Lambert W function in the environmental sciences include an exact expression for the lifting condensation level (Romps, 2017;Yin et al., 2015), a closed-form solution for convective available potential energy (Romps, 2016), a solution for the temperature profile in an idealized model of moist convection (Vallis et al., 2019), a solution of Richards' equation for unsaturated soil water transport (Barry et al., 1993), and various applications in ecology (Lehtonen, 2016).
Equation (8) applies to the radiatively uncoupled surface energy budget. An equivalent equation is derived for the radiatively coupled case in Appendix B.

Practical Benefits of Equation (8)
In this section, eddy covariance observations are used to assess the accuracy of equation (8) across a wide range of conditions. The performance of equation (8), in terms of root-mean-squared error (RMSE), is shown to be better than that of the PM equation at all sites available for the analysis.

Data
Observations of E, R n , G, T a , wind speed (u [m s −1 ]), relative humidity (RH = q a ∕q * (T a ) [-]), air pressure (P [Pa]), and friction velocity (u * [m s −1 ]) were obtained from the FLUXNET database (fluxnet.ornl.gov). The FLUXNET data set includes half-hourly observations from eddy covariance sites around the world. Sites that did not include all of these observations were excluded. Observations with quality control flags corresponding to poor-quality gapfill were removed. Estimates of vegetation height at each site were obtained from Lin et al. (2018). Both daytime and nighttime observations were used in this analysis. Nighttime observations were filtered using site-specific friction velocity thresholds provided with the FLUXNET database (Barr et al., 2013;Papale et al., 2006). Seventy-nine eddy covariance sites met these requirements and were used in the analysis; information on these data is summarized in Table S1.
The aerodynamic conductance g a was estimated at each site using the standard relation (Garratt, 1994;Thom & Oliver, 1977) where k = 0.41 is the von Karman constant [-]; z is the measurement height [m]; u is the mean wind speed [m s −1 ] at height z; d = 2 3 h is the assumed zero-plane displacement height [m]; h is the vegetation height [m]; z 0h and z 0m are the thermal and momentum roughness heights [m], respectively (chosen to be 0.01h and 0.1h, respectively, consistent with previous studies (Lin et al., 2018)); and M and H are the stability correction functions for momentum and heat transfer, respectively [-]. Standard relations are used for M and H for unstable (Paulson, 1970) and stable conditions (Holtslag & De Bruin, 1988).
For eddy covariance observations, the ecosystem-scale surface conductance g s is often estimated using a rearranged form of the PM equation (e.g., Lin et al., 2018;Medlyn et al., 2017;Novick et al., 2016;Wullschleger et al., 2002). The aim of this analysis is to quantify errors in this equation, so a different approach is required: instead, the radiatively uncoupled surface energy budget (equation (1)) is solved numerically for g s at each point in time and space. Specifically, observations of H (estimated as R n − G − E) and T a , combined with estimates of g a , are used to estimate T s , based on the definition of H in equation (1). The estimated T s is then combined with observations of E and q a to estimate g s , using the definition of E in equation (1).

Validation Against Observations
Errors in equation (8) and the PM equation are estimated using observed values of E, R n − G, T a , q a , and estimated values of g a and g s , as described in the previous section. Overall, equation (8) consistently outperforms the PM equation. At all 79 sites, compared with the PM equation, equation (8) has lower RMSE with respect to observed E. Across sites, and using both daytime and nighttime observations, equation (8) results in a 67% reduction in RMSE, reducing it from 8.55 to 2.81 W m −2 . At night (R n − G < 0), absolute values of RMSE are lower, since absolute values of E are lower. However, across sites, the average reduction in RMSE using equation (8) at night is 92%, reducing from 1.26 to 0.0974 W m −2 (Figures 1b and 1d). Restricting observations to daytime only (R n − G > 0) results in a reduction in RMSE using equation (8) (2)) systematically underestimates q * (T s ) (Figure 1a), causing the PM equation to systematically underestimate E. In contrast, the approximation of q * (T s ) used in the derivation of equation (8) (equation (6)) slightly overestimates q * (T s ) (Figure 1a), causing equation (8) to overestimate E.
At the half-hourly timescale, the PM equation occasionally performs slightly better than equation (8), but these cases are rare and errors are small. More specifically, for 87% of the half-hourly observations across all sites, equation (8) performs better than the PM equation in terms of RMSE (the RMSE is just the absolute value of one residual at the half-hourly timescale, since there is one observation per half hour). In these cases, the median absolute difference between the half-hourly RMSE for equation (8) (8) performs better than the PM equation at the half-hourly timescale in a broad majority of cases. When it does not perform better, the difference in performance is very small. When aggregated to longer timescales, the performance of equation (8) is consistently better than that of the PM equation.
These results do not appear to be an artifact of errors in the observations or in the methods used to estimate g a and g s from the data. To check this, first, all analyses were repeated using E estimated as the residual of the observed energy balance at each site, rather than the directly observed value. The results were qualitatively similar (not shown). Second, the equations were also tested on synthetic "observations" (including synthetic observations of g a and g s ), which are not subject to observation error, or other estimation errors (supporting information Text S1). The results are qualitatively similar to those obtained using real data ( Figures S4-S5).

10.1029/2020WR027106
In summary, equation (8) substantially and consistently reduces errors in estimates of E based on real-world observations, compared with the PM equation. It does so without calibration and without requiring any additional assumptions, inputs, or computational cost.

Theoretical Benefits of Equation (8)
Important theoretical insights can be gained from analyzing explicit solutions of equation (1), such as the PM equation or equation (8). One major theoretical benefit of an explicit solution over numerically solving equation (1) is that the explicit solution can be differentiated. For example, van Heerwaarden et al. (2010) differentiated the PM equation to analyze forcings and feedbacks in the coupled land-atmosphere system. Another benefit of an explicit solution is that limiting cases-that is, cases in which a governing parameter, such as g a or g s , approaches zero or infinity-can be studied analytically. The study of limiting cases of evapotranspiration has a long history (e.g., Paw U & Gao, 1988;Raupach, 2001;Yang & Roderick, 2019). For example, a common example of a limiting case in the study of E is the concept of "potential evapotranspiration": one of many definitions of potential ET is E in the limit of g s → ∞. Another definition is the case in which, in addition to g s → ∞, the near-surface air is also saturated. Under this definition, the radiatively uncoupled PM equation converges to the "equilibrium" value, Raupach (2001) provides a comprehensive review of equilibrium ET. Conveniently, this relation does not depend on wind speed or surface conditions, including surface temperature (this applies for the radiatively uncoupled case (Raupach, 2001)). The relation has been empirically adapted to the more common case in which the air is unsaturated, by multiplying the equilibrium value by a constant, often taken to be 1.26 (Priestley & Taylor, 1972). This variant of equilibrium ET is at the core of modern ET estimation algorithms, such as the Priestley-Taylor Jet Propulsion Lab (PT-JPL) algorithm (Fisher et al., 2008)  Given the considerable importance of equilibrium ET and other limiting cases in understanding and modeling ET, any approximate solution of E based on equation (1) should reproduce the correct limiting behavior. In the following section, I show that, in contrast to equation (8), the PM equation does not reproduce the correct limiting behavior in all cases.

Limiting Behavior of the PM Equation and Equation (8)
In this section, I compare predictions of E in several limiting cases based on the PM equation and equation (8), with the true limiting behavior implied by the radiatively uncoupled surface energy budget (equation (1)). The radiatively uncoupled case is of primary interest in this section because, as we will see, two common definitions of equilibrium ET are based on limits of the radiatively uncoupled surface energy budget. I will show that one of these definitions is incorrect.
For readers not familiar with limits, the expression " E → X as g a → 0" can be interpreted as "X is a reasonable approximation of E when g a is sufficiently small." The expression " E → X as g a → ∞" can be interpreted as "X is a reasonable approximation of E when g a is sufficiently large." When considered this way, a limit has clear practical relevance as a useful approximation in real-world cases.
While various limiting cases of the surface energy budget could be considered, this study focuses exclusively on the limiting cases in which g a and g s approach both zero and infinity, consistent with previous studies (Raupach, 2001).

The Wet Limit: g s → ∞
The wet limit can be a reasonable approximation of E from a saturated surface, such as a lake or saturated soil. We consider two variants of this case: one in which the air is also saturated (D = 0 and g s → ∞) and one in which it is not (D > 0 and g s → ∞).
For the case where D = 0, the surface energy budget can be rewritten exactly as which does not contain explicit dependence on g a . However, there is implicit dependence, since T s is an implicit function of g a , and so the limiting value varies with g a . The PM equation approximates sa ≈ , removing all dependence on T s and thus g a . While this turns out to be a reasonable approximation in this limit, technically, there is still some dependence on g a . Equation (8) is a function of g a in this limit and is numerically more accurate, compared to the PM equation, as shown numerically in Figure S7.
For the case where D > 0, the PM equation systematically underestimates the true value in this limit. This result is illustrated numerically using synthetic observations, since there is no closed-form solution to equation (1) in this case. Synthetic "truth" observations are generated by solving equation (1)

The Dry Limit: g s → 0
The dry limit can be a reasonable approximation of E in an arid environment, where water availability limits E. Physically, E should be zero, since the surface conductance is limiting. Both the PM equation and equation (8) give the correct limiting behavior: E → 0 as g s → 0.

The Rough Limit: g a → ∞
The rough limit can be a reasonable approximation of E over a rough surface, such as a forest. Physically, the temperature gradient approaches zero (T s → T a ) to maintain finite sensible heat flux. This implies that E = g s g s ∕g a +1 (q * (T s ) − q a ) → g s (q * (T a ) − q a ) as g a → ∞. Both the PM equation and equation (8) correctly reproduce this limit (see Appendix C for a derivation).

The Calm Limit: g a → 0
This limit corresponds to a case in which turbulent diffusion becomes small. The radiatively uncoupled PM equation (equation (3)) implies that the calm limit (g a → 0) is exactly equivalent to the "wet" limit (g s → ∞ and D = 0), with E converging to the "equilibrium" value (equation (10)) in both cases. It is still common to find references to equilibrium E as the radiatively uncoupled limit in which g a → 0 (e.g., Jones, 2014;Raupach, 2001). From this, a puzzle arises: For a given observed value of R n − G, why would E from a saturated surface into a saturated, turbulent atmosphere be exactly equivalent to E from an unsaturated surface into an unsaturated, laminar atmosphere, in general? If this were true, it would require a deep and fundamental connection between turbulent and laminar processes.
In this section, I show that this apparent equivalence is incorrect. Paw U and Gao (1988) showed the limiting behavior produced by the PM equation for the calm limit is incorrect during daytime conditions, and I extend that analysis here to show that it is also incorrect at night.
The PM equation converges to equilibrium E (equation (10)) in two cases. The first case is a wet limit, in which g s → ∞ and D = 0. This case corresponds to the original definition of equilibrium E, is physically justified, and predates the PM equation. Figure S3 shows two example cases, corresponding to daytime (R n − G > 0, Figure S3b) and nighttime (R n − G < 0, Figure S3d), comparing the solution of the PM equation (equation (3)) to the exact solution obtained from numerically solving the surface energy budget (equation (1)), for different values of g s . All solutions correctly converge to the equilibrium value as g s → ∞ when D = 0.
The second case is a calm limit, in which g a → 0. Physically, this limit is approached when there is little atmospheric turbulence (i.e., the atmosphere approaches a laminar state). Diffusion of water vapor from the land surface in the absence of turbulence is slow. This is one of several textbook definitions of equilibrium E (e.g., Jones, 2014), first proposed by Thom (1975). It is derived directly from the PM equation. However, there is no obvious physical reason why this case should be equivalent to the first case. In fact, it is not and is an artifact of the PM equation. To illustrate this, Figures S3a and S3c compare the PM solution to the true solution for different values of g a , for daytime (R n − G > 0) and nighttime (R n − G < 0) conditions, respectively. As g a → 0, diffusion of heat and water vapor from the surface becomes very slow and inefficient.

Radiatively uncoupled
Radiatively coupled when RH = 1 E PM underestimates when E PM underestimates when RH < 1 and R n − G < 0 RH < 1 and R * n − G * < 0 Note. True values are given in black text. Equation (8) is correct in all limiting cases (see Appendix C for analytical derivations). Where E PM incorrectly diverges from the true limit, it is noted with red text.
The land surface is unable to turbulently transport away the incoming energy, and as a result, |T s | increases substantially. For R n − G > 0, the nonlinear Clausius-Clapeyron relation q * (T s ) in the definition of E increases much more rapidly than T s , the equivalent term in the definition of H. In the limit of T s → ∞, E dominates H and consumes all available energy at the surface (Bateni & Entekhabi, 2012;Yang & Roderick, 2019), resulting in the limit E → R n − G as g a → 0. For R n − G < 0, as g a → 0, T s decreases substantially, rather than increasing. In this case, T s decreases much more rapidly than q * (T s ), which asymptotes to zero; therefore, H dominates E and consumes all available energy at the surface, resulting in the limit E → 0 as g a → 0. In comparison, as g a → 0, the PM equation incorrectly approaches the equilibrium value, for both positive and negative R n − G.
Why do these artifacts occur? The major assumption behind the PM equation is that a linear approximation of the Clausius-Clapeyron relation (equation (2)) is accurate. This approximation is reasonable when |T s − T a | is small but can be extremely inaccurate when this assumption is violated. Figure 1a gives an illustration of this. For small |T s − T a |, the linear approximation used in the PM equation is quite accurate. However, as the difference between T s and T a grows, the accuracy of the linear approximation (and even higher-order quadratic and quartic approximations) degrades. In particular, for the case where T s ≪ T a , none of the linear, quadratic, or quartic approximations approach the correct limiting value of q * (T s ) → 0; in fact, errors in the higher-order quadratic and quartic approximations grow more rapidly than those for the linear approximation. Since |T s − T a | grows very large in the calm limit, this suggests that the PM equation should fail to reproduce the correct limiting behavior.

Summary
In summary, the apparent equivalence between E in the radiatively uncoupled calm (g a → 0) and wet (g s → ∞ and D = 0) limits is purely an artifact of the PM equation. A major theoretical advantage of equation (8) is that it avoids this problem: In particular, Figure S3 shows that, unlike the PM equation, equation (8) correctly predicts that E → R n −G as g a → 0 when R n −G > 0 and that E → 0 when R n −G < 0. It also reproduces the correct limiting behavior when g a → ∞, g s → 0, and g s → ∞. The correct limiting behavior for each case is summarized in Table 1, and inaccurate limits in the PM equation are highlighted in red. It is shown analytically in Appendix C that equation (8) converges to the correct limits in each case.

A Revised Decoupling Parameter
By correctly representing limiting behavior, equation (8) can be used to correct a conceptual mistake in a theory of land-atmosphere coupling. Jarvis and McNaughton (1986) argued that E could be scaled between two limiting cases: an "equilibrium" case, in which the atmospheric state was set by the evaporating surface, and an "imposed" case, in which the atmospheric state was independent of the evaporating surface. To quantify this scaling, they proposed a "decoupling parameter" Ω, defined as 10.1029/2020WR027106 or equivalently as where E eq is equilibrium latent heat flux and E imp is the latent heat flux in the limit of g a → ∞. Based on this definition, they proposed, using the PM equation, that Ω could be estimated aŝ According to this definition, Ω → 1 (and, therefore, E → E eq ) when g s → ∞ or when g a → 0. As shown previously, the latter case is an artifact of the PM equation and is incorrect (Paw U & Gao, 1988).
Using equation (8), this problem can be resolved. I define E eq as the latent heat flux in the wet limit (g s → ∞ and D = 0) in equation (8). E imp is defined as the latent heat flux in the rough limit (g a → ∞), the definition given by Jarvis and McNaughton (1986), in equation (8). Equation (8) is correct in both limiting cases, unlike the PM equation. Substituting these expressions along with equation (8) into equation (12) gives the following estimate: ) . (14) While complicated, this expression does not require any additional information or assumptions beyond those made already by Jarvis and McNaughton (1986). As g a → 0, both E and E eq approach R n − G; therefore,Ω → 1, by equation (12), as for the formulation of Jarvis and McNaughton (1986). The key difference is, in the new formulation, E eq goes to the correct limit (R n − G) as g a → 0. This remedies a significant conceptual weakness in the theory.
This discussion has been focused on the radiatively uncoupled case, since this case is most relevant to previous definitions of equilibrium E (equilibrium E is not a limiting value in the radiatively coupled case) and has also been used most widely (e.g., Fisher et al., 2009;De Kauwe et al., 2017). The behavior of the decoupling parameter is different in the radiatively coupled case (Paw U & Gao, 1988), and the decoupling parameter framework has since been generalized to the radiatively coupled case (Martin, 1989;McNaughton & Jarvis, 1991). I provide a derivation of a radiatively coupled equivalent expression to equation (14) in Appendix B.1.

Revisiting Equilibrium E
In this section, I critically reexamine previous definitions of equilibrium E in light of the presented results. This work implies that several conventional definitions of equilibrium E are incorrect (i.e., are not equivalent to, or even approximations of, equation (10)). Raupach (2001) provides a comprehensive review of equilibrium E, listing five different definitions prevalent in the literature. Here, I critically reassess these definitions, finding that two definitions are clearly incorrect, one is probably incorrect, one is probably correct, and one is correct when considered to be a useful approximation. The definitions are as follows: • E in the limit of g a → 0: This is incorrect, as discussed in previous sections.
• E in the limit of complete decoupling (Ω = 1): This definition is also based on the PM equation and is also incorrect (Paw U & Gao, 1988). Even ifΩ is appropriately redefined, as in equation (14),Ω = 1 does not imply equilibrium E. For example, in the radiatively uncoupled case given above,Ω → 1 as g a → 0, but E does not go to the equilibrium value. • E that is independent of g a : This relation, proposed by Monteith (1965) and Thom (1975), was obtained by differentiating the PM equation with respect to g a , setting the resulting expression to zero, and rearranging to solve for E. While equation (8) is differentiable, the same procedure does not yield a closed-form expression when applied to equation (8). Furthermore, it is not clear when a minimum exists in E for finite g a , in general: For example, one exists in Figure S3c but not Figure S3a. It therefore seems likely that this definition is also an incorrect artifact of the PM equation. • E of a closed system forced with incoming energy and allowed to evolve over a sufficiently long period of time: Previous studies that established this relation used either the PM equation explicitly (Raupach, 2001) or other relations based on the linearization of the Clausius-Clapeyron relation (McNaughton, 1976a(McNaughton, , 1976bSlatyer & McIlroy, 1961). This provides some cause for concern. However, since the results of these studies hold for any value of g a and g s , and problems with the PM equation mainly arise in limiting cases, this suggests that the definition and results of these studies are likely to be broadly correct. However, the problem requires further consideration, which is left to future work. • E in the limit of g s → ∞ and D = 0: This definition is exactly correct when sa is used in the definition of equilibrium E (Raupach, 2001). When is used instead, it is typically an accurate approximation over wet surfaces (Milly, 1991).

Summary and Conclusions
This study has introduced an alternative to the PM evapotranspiration equation (equation (8) Theoretically, the PM equation is shown to be incorrect in several important limiting cases, which has led to incorrect definitions of equilibrium E: In particular, the definition of equilibrium E as the limiting value of E in the radiatively uncoupled surface energy budget as g a → 0, which is an incorrect artifact of the PM equation. The new equation does not suffer from this problem. I use the new relation to show that several other common definitions of equilibrium E are incorrect and to remedy a related conceptual error in the decoupling parameter theory of Jarvis and McNaughton (1986).
Beyond the practical and theoretical applications considered here, this work opens up opportunities for more accurately studying the surface energy budget at night, and in cold environments, where the PM equation is particularly inaccurate.

Appendix A: Derivation of Equation (8)
In this section, I solve equation (7) for T s and obtain equation (8). Equation (7) can be rearranged to the form (A4) Based on the definition of the Lambert W function (equation (2)), this gives

10.1029/2020WR027106
Substituting in equations (A2, A3, A4) yields By the surface energy balance (equation (1)), E = R n − G − c p g a (T s − T a ). Substituting equation (A5) into this yields equation (8). Figure S3 is reproduced for the radiatively coupled case in Figure S6. For R * n − G * > 0, the limiting behavior for g s is similar, although the limit as g s → ∞ is E Figure S6b, Raupach, 2001). The limiting behavior as g a → 0 is quite different for the radiatively coupled case. This is because, in the calm limit of the radiatively coupled case, T s remains finite, with outgoing longwave radiation balancing incoming radiation, and consequently, H and E both approach zero ( Figure S6a). Equation (B5) is accurate in all cases.

B.1. Decoupling Parameter
For the radiatively coupled case, the decoupling parameter is given bŷ (g a + g s ) Based on a similar derivation to that in Appendix C for the limiting behavior of E, it can be shown that, as g a → 0,Ω → 1 and E → 0.
Appendix C: Limiting Behavior of Equation (8) in the Limits g s → ∞, g s → 0, g a → ∞, and g a → 0 Raupach (2001) comprehensively characterizes the limiting behavior of both the radiatively coupled and radiatively uncoupled PM equation. In this section, I characterize the limiting behavior of equation (8) in both radiatively coupled and radiatively uncoupled cases.

C.1. g s → ∞ and D = 0
In this case, for both the radiatively uncoupled and radiatively coupled cases, where p = 1 gives the expression for the radiatively uncoupled case, and 0 ≤ p < 1 describes the radiatively coupled case.
Simulations conducted in this limit ( Figure S7) demonstrate that, while the assumed form in PM gives a reasonable first-order estimate, it is entirely insensitive to variability due to varying g a , as expected. In contrast, equation (8) is more accurate and captures reasonably the sensitivity of E to variation in g a , even in the limit of g s → ∞ and D = 0.

C.3.1. Radiatively Uncoupled Case
As x → 0, W 0 (x) ∼ x. In addition, g s g a g s +g a = g s g s ∕g a +1 → g s as g a → ∞. These results will be used in this section. As g a → ∞, equation (8) goes to − g s q a = g s (q * (T a ) − q a ).

C.3.2. Radiatively Coupled Case
As g a → ∞, p → 1, and so the radiatively coupled case is identical to the radiatively uncoupled case.
As g a → 0, equation (8) goes to Since both the numerator and the denominator diverge to infinity as g a → 0, L'Hôpital's rule is used to evaluate the limit. Differentiating the numerator and the denominator and simplifying give lim g a →0 Term III = c p g a ( (R n − G) + ( g a g s +g a ) 2 (g s 2 q * (T a ) + c p (g a + g s )R v T 2 a ) R v T 2 a ( (R n − G) + c p g a R v T 2 a log(

C.4.2. Radiatively Coupled Case
As g a → 0, p → 0 and g a ∕p → g g + g r . Therefore, equation ( ∕(R v T 2 a ) − g s g a g s + g a q a → c p (g g + g r ) W 0 (0) − 0 = 0 since W 0 (0) = 0. The key difference for the radiatively coupled case is that p varies with g a , whereas in the radiatively uncoupled case, p = 1 and is invariant to changes in g a . Thanks to Pierre Gentine and Daniel Short Gianotti for providing feedback on a draft manuscript, Angela Rigden for assistance with processing the FLUXNET data, Changjie Lin for providing vegetation height information at each FLUXNET site, and Trevor Keenan for providing FLUXNET citations and site information via his FLUXNET_citations package on GitHub. This work was supported by funding from a Winokur Seed Grant in Environmental Sciences from the Harvard University Center for the Environment. The eddy covariance data used in this study can be accessed at https://daac.ornl.gov/cgi-bin/ dataset_lister.pl?p=9. This work used eddy covariance data acquired and shared by the FLUXNET community, including these networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux-TERN, TCOS-Siberia, and USCCC. The ERA-Interim reanalysis data are provided by ECMWF and processed by LSCE. The FLUXNET eddy covariance data processing and harmonization was carried out by the European Fluxes Database Cluster, AmeriFlux Management Project, and Fluxdata project of FLUXNET, with the support of CDIAC and ICOS Ecosystem Thematic Center, and the OzFlux, ChinaFlux, and AsiaFlux offices.