Complementary principle of evaporation: From original linear relationship to generalized nonlinear functions

The complementary principle is an important methodology for estimating evaporation. Throughout the 56-year development, related studies have shifted from adopting a symmetric linear complementary relationship (CR) to employing generalized nonlinear functions. Studies based on the linear CR have been maintained for a long time by rationally formulating 10 the potential (Epo) and "apparent" potential evaporation (Ep) and/or employing an asymmetric parameter. These works have also advanced two types of generalized nonlinear complementary functions by invoking the boundary conditions. The first type inherits the concepts of three types of evaporation yet still requires the prognostic modelling of Epo. Polynomial functions are derived and tested for this type of function. Meanwhile, the second type does not involve Epo, yet requires a diagnostic modelling of actual evaporation by using the radiation and aerodynamic components of the Penman (1948) equation as inputs. 15 A sigmoid function is derived by satisfying the boundary conditions based on physical considerations. The generalized nonlinear functional approach has improved the understandings on the complementary principle, and shows potential in advancing the evaporation research. Further studies may cover several topics including boundary conditions, analytical forms, parameterization, and application.


Introduction 20
The complementary principle provides a framework for estimating terrestrial land surface evaporation by adopting routinely observed meteorological variables, and offers strong potential applications (Brutsaert and Stricker, 1979;Morton, 1983;McMahon et al., 2016). In this paper, the terms "evaporation" and "evapotranspiration" are considered equivalent. As its underlying physical basis, this principle describes the feedback of areal evaporation on evaporation demand (Bouchet, 1963;Brutsaert, 2015) as illustrated by the fact that reducing areal evaporation can make the overpassing air hotter and drier (Morton, 25 1983). Based on the complementary principle, Bouchet (1963) first proposed a complementary relationship (CR) among three types of evaporation (Brutsaert, 2015), namely, the actual evaporation (E) from an extensive landscape under natural conditions by relating the apparent potential evaporation ( pa E ) of a small saturated surface inside the landscape that does not affect the overpassing air and the natural evaporation process, and the potential evaporation ( po E ) that occurs from the same large-size ( 1 ) pa E and po E should be specified in Eq. (1). Bouchet (1963) assumed po E to be half the input solar radiation. Morton (1976) calculated pa E by using the modified Penman's (1948) equation proposed by Kohler and Parmele (1967)( KP Pen E ), in which a constant vapor transfer coefficient was used to replace the wind function, and calculated po E by using the Priestley-Taylor's 55 (1972) Morton (1971) is physically unrealistic and added that the complementarity of the negative relationship is not supported by any proof.
The physical basis of the CR has been further advanced by using climate models. McNaughton and Spriggs (1989) Kim and Entekhabi (1997) added the surface energy balance and atmospheric thermal radiation fluxes into the model to extend the study of McNaughton and Spriggs 105 (1989). By using the Penman-Monteith equation to govern the areal latent heat flux at the surface, Lhomme (1997a) proposed a closed-box model with an impermeable lid at a fixed height while Lhomme (1997b) used a realistic open-box model of the ABL with entrainment to assess the CR. Sugita et al. (2001) tested the CR by using a modified version of Lhomme (1997b) model, which was calibrated by using a dataset obtained from the Hexi Corridor desert area in Northwest China. But a strict symmetric CR wss hardly confirmed by these studies. 110

Rational formulation of and/or
The imperfect asymmetric CR has inspired researchers to apply a rational formulation of pa E and/or po E for retaining the symmetry of CR. One direct method is to improve the formulations of and/or based on the AA approach through calibration. For Pen E , the empirical wind function was calibrated to improve the CR (Hobbins, 2001). However, Penman's 115 wind function cannot work under wet and dry conditions simultaneously (Pettijohn and Salvucci, 2006). The wind function derived from Monin-Obukhov's similarity theory was then employed (Crago and Crowley, 2005;Parlange and Katul, 1992b;Pettijohn and Salvucci, 2006;Ma et al., 2015a). The surface roughness and surface albedo were also calibrated to improve the CR (Lemeur and Zhang, 1990). Meanwhile, for PT E , the Priestley-Taylor coefficient (  ) is regarded varying, thereby leaving a range for calibration (Han et al., 2006;Yang et al., 2012;Xu and Singh, 2005). In addition to Pen E and PT E , the mass-120 transfer type potential evaporation (van Bavel, 1966) ( MT E ) was considered another formulation of potential evaporation (Granger, 1989 Crago and Crowley (2005). The local re-parameterizations and/or calibrations have significantly improved the evaporation estimation (Hobbins et al., 2001b;Xu and Singh, 2005;Ma et al., 2015a). However, the calibration approach and trial-and-error process are deemed ineffective because of their high computation demand, which is a key 125 stumbling block when applying the CR in large-scale (e.g., continental or global) E modelling (Ma et al., 2019).  (Morton, 1983). To address these limitations, Morton (1983) 130 combined the energy balance and water vapor transfer equations by using an equilibrium temperature (Tp) and derived a Morton-type potential evaporation to denote pa E . By attributing the asymmetry to the assumption that pa E conceptually includes a transpiration component, Pettijohn and Salvucci (2006) (Han et al., 2014d;Han et al., 2017). 135 PT E was proposed by Priestley and Taylor (1972) to represent evaporation from extensive saturated surfaces and has been widely used (Brutsaert, 1982;Priestley and Taylor, 1972). This way it could be used to represent po E (Brutsaert and Stricker, 1979). The AA approach calculates PT E by using the atmospheric variables that correspond to the current natural landscape. However, the atmosphere in contact with the land surface will change if the land surface was brought into saturated (Morton, 1983;Brutsaert, 2015). Therefore, PT E cannot represent the "true" po E since the slope of the saturation vapor 140 pressure at the current air temperature ( () a T  ) is imperfect because the temperature corresponding to po E is different from that corresponding to a non-saturated environment (Morton, 1983;Szilagyi and Jozsa, 2008). Moreover, PT E does not fully consider the effects of advection, which are inevitable in reality (Morton, 1983(Morton, , 1975Parlange and Katul, 1992a). To this end, Morton (1983) derived po E by using a modified Priestley-Taylor equation with net radiation and the slope of the saturation vapor pressure that is calculated at equilibrium temperature Tp ( p T PT E ). The effects of advection were considered by an empirical 145 correction factor in p T PT E (Morton, 1975(Morton, , 1983. Parlange and Katul (1992a) attributed the asymmetry to the horizontal advection of dry air, which would make Pen E larger than the available energy ( n RG  ) (i.e., 0 n Pen to improve the CR on an hourly basis. Szilagyi and Jozsa (2008) argued that Δ in PT E should be calculated at the air temperature corresponding to the wet environment (Twa) instead of actual air temperature. While it is not straightforward to derive Twa, Szilagyi and Jozsa (2008) proposed an iterative approach based on the Bowen 150 ratio method for a small wet patch to estimate the surface temperature under wet environments ( ws T ). By assuming a negligible temperature gradient over such a small wet area, Tibetan Plateau (Ma et al., 2015a). The evaporation estimations across the US were also improved by applying the modified 155 AA approach (Szilagyi and Jozsa, 2008;Szilagyi et al., 2009;Szilagyi, 2015).

Asymmetric linear complementary relationship
With pa E and po E denoted by MT E and Pen E , respectively, Granger (1989) proposed an alternative CR as follows: where  is the psychrometric constant. Despite being identical to the surface energy balance, Eq.
(2) has inspired researchers 160 to examine whether the CR should be symmetric (Szilagyi, 2007;Pettijohn and Salvucci, 2006). By using pan evaporation to denote pa E , Brutsaert and Parlange (1998) extended the symmetric CR as follows: where b is the coefficient that denotes asymmetry. Kahler and Brutsaert (2006) clarified and tested the validity of Eq.
(3) at a daily timescale and attributed the asymmetry to the nature of the heat transfer between the pan and its surroundings, 165 which made the changes in pan E larger than those in E . Szilagyi (2007) showed that the asymmetry is not limited only to pan E but is inherently linked to the definition of pa E . The asymmetric CR is widely used, and Brutsaert (2015) stated that asymmetry is an inherent characteristic of the CR. The parameter b is often considered a calibrated parameter in evaporation estimation (Ma et al., 2015a;Szilagyi, 2015), and its controlling factors were also detected (Szilagyi, 2015;Lintner et al., 2015;Szilagyi, 2007). 170 The asymmetric CR can be illustrated in a dimensionless form ( Figure 2) (Kahler and Brutsaert, 2006 The scaled pa E and E are both functions of the dimensionless variable po EE , while po EE serves as the evaporative surface moisture index. Compared with the original form (Eq. (1) and Figure 1), the CR here is illustrated without the appearance of the water availability explicitly.

Normalized complementary functions
Unlike the normalization by po E (Kahler and Brutsaert, 2006), Han (2008)   , the AA approach can be expressed as where Pen EE is regarded as a linear function of rad Pen EE. The bias of the AA function under arid and wet environments can be easily understood in its dimensionless form, but the AA approach with a tuned b still underestimated evaporation in arid environments . The work of Crago and Brutsaert (1992) in Kansas during FIFE 1987 revealed that the parameter b are obviously different for days with differing degrees of soil moisture. These studies imply that the CR may deviate from its linear characteristics. 190 The CR model proposed by Granger (1989) based on Eq.
(2) has demonstrated promising application across different land covers and regional climate conditions (Carey et al., 2005;Granger, 1999;Granger and Gray, 1989b;Pomeroy et al., 1997;Xu and Singh, 2005). In fact, the relationship between relative evaporation and relative drying power plays a key role in reflecting the dryness of the surface (Granger and Gray, 1989a). Normalized by Pen E , Granger's model is similar to the AA function in that Pen EE is expressed as a function of the relative magnitude of drying power to net radiation (Han et al., 2011). 195 By synthesising the dimensionless forms of the AA function and the Granger's model, Han et al. (2011) proposed the following logistic function as an alternative: (1 ) where 1 c and d are the parameters. Eq. (6) approximates the linear AA function under normal conditions neither too wet nor too dry but amends its bias (Han et al., 2011). Actual evaporation can be estimated using routinely measured meteorological data by using the climatological resistance to parameterize the bulk surface resistance in the Penman-Monteith equation (Liu et al., 2012;Rana et al., 1997;Katerji and Perrier, 1983;Ma et al., 2015b). A linear relationship between the ratio of surface resistance to aerodynamic resistance and the ratio of climatological resistance to aerodynamic resistance was proposed by Katerji and Perrier (1983). Han et al. (2014c) integrated this linear relationship into the Penman-Monteith equation and derived a dimensionless form via normalization by 205 where k and l are the empirical calibration parameters. With similar variables yet different mathematical formulations, Eq. (7) can also be considered a complementary function (Han et al., 2014c).

Sigmoid function relating E/EPen to Erad/EPen 210
By synthesizing the aforementioned three functions, Han et al. (2012)  indirectly denoted by the drying power of air with a constant radiation energy input (Brutsaert, 1982). Accordingly, Eq. (8) is considered a "general form" of the CR (Han et al., 2014b) (hereinafter referred to as H12). The existing analytical forms of the function can be classified into linear, concave, or sigmoid (Table 1 in Han and Tian (2018a)). Studies on the complementary principle can be advanced by formulating a proper analytical form for H12. 220 The exact analytical form of H12 is inadequately understood at present. However, some of its characteristics can be detected from its boundary conditions under extremely arid and completely wet environments. Han et al. (2012) where m and n are parameters. The linear AA and nonlinear H2012 have been compared in a 2D space ( rad Pen EE , Pen EE ) (Han et al., 2012). The results obtained from an extremely dry desert and a wet farmland reveal that the sigmoid H2012 corrects the bias of the linear AA and Equation (6) (Han et al., 2012); the application of this sigmoid function has also been 230 recommended for an alpine meadow region of the Tibetan Plateau (Ma et al., 2015b).  Szilagyi et al., 2017;Ková cs, 1987 Based on the boundary conditions, Han and Tian (2018a) found that the growth of where rad Pen EE adopts the feasible domain ( min x , max x ), which is a subdomain of (0, 1). Both the linear AA function and H2012 are special cases of H2017. Han where c is a parameter. Brutsaert (2015)   13 Ai et al., 2017). In this case, we refer to Eq. (14) in the manner of the AA approach as B2015 to avoid confusion. Although x  of B15 may not hold in the manner of the AA approach (Ková cs, 1987;Szilagyi et al., 2017;Crago et al., 2016;Han and Tian, 2018a). To address these challenges, Szilagyi et al. (2017); Crago et al. (2016) (Crago and Qualls, 2018). After rescaling, Crago et al. (2016) proposed a new linear version of the generalized complementary function (C2016) (i.e., BC yx  ; Table 2), while Szilagyi et al. (2017) used the third order polynomial function (S2017) by replacing B15 with c=0. With the same independent variable yet different functions (Table 2), C2016 and S2017 demonstrate improvements in their evaporation estimation performance 295 (Szilagyi et al., 2017;Crago et al., 2016;Crago and Qualls, 2018).

Comparisons between different analytical forms
The linear AA, the polynomial B2015, and the sigmoid H2017 are three analytical forms of H12 (Table 1 in Han and Tian (2018a)). Han and Tian (2018a) (Brutsaert, 2015) while considering that the wet regional 300 evaporation must always be smaller than the wet patch evaporation (  OMN (Figure 3). The limits of Pen E and PT E on E can be approximately satisfied by H2017 with the parameters transformed from the linear AA function (Han and Tian, 2018a). By contrast, B2015 with 1 c  produces a physically unreasonable E that is larger than PT E (Han and Tian, 2018a and polynomial B2015 have one and two stages respectively. As it is difficult for one site to cover all the three stages with a wide range of wetness, the linear AA can effectively represent the complementary curve under normal conditions falling in the middle stage. The polynomial B2015 is effective if the first two stages exist. Given that the third stage under a wet environment is uncommon, the polynomial B2015 performs well with calibrated parameters Liu et al., 2016;Zhang et al., 2017;Han and Tian, 2018a). However, the sigmoid H2017 shows the best performance in estimating 315 evaporation as validated by using data from FLUXNET (Han and Tian, 2018a).

Improved understanding on the correlation between actual and potential evaporation
Interpreting the changes in E based on the trends in Pen E (or pan evaporation) greatly relies on the understanding of whether the correlation between E and Pen E is positive or negative. The corresponding confusion has resulted in a discrepancy between the Penman hypothesis and the complementary principle (Yang et al., 2006) and encouraged debates on whether the 320 increasing or decreasing trend in E corresponds to reductions in the observed pan evaporation in the past (Brutsaert and Parlange, 1998;Roderick and Farquhar, 2002;Roderick et al., 2009;Wang et al., 2017). According to the symmetric CR, Pen E would be negatively correlated with E when the energy input is constant (Morton, 1983). Based on the asymmetric linear CR, Brutsaert and Parlange (1998) stated that the decreasing pan E can be used to indicate an increasing E in water-limited regions. However, the interpretation is not general (Roderick et al., 2009) because of the inherent weakness of the linear CR (Han et 325 al., 2014b). conditions where the variations in rad E and aero E are comparable or when aero E obviously varies (which tends to occur on a daily or annual basis), the influence of aero E comes to force. As a result, the correlation between E and Pen E changes from 335 negative to positive along with increasing water availability. The theoretical results were validated in a grassland site in Northeast China, and can rationally interpreted the trends in E over China (Han et al., 2014b).  Han and Tian (2018a) found that the first-order wet boundary conditions of H12 and B2015 are two possible solutions of the assumption that actual evaporation proceeds at Pen E , and they discarded each other. However, there was not perfect instance to demonstrate H12's first-order wet boundary condition. Thus, the controversies on the first-order wet boundary condition Han and Tian, 2019) 350 require further studies from both the theoretical and practical aspects. Thus, the flux data of sites over lakes or wetlands need to be well examined.

Parameterizations of generalized complementary functions
Determining the parameters of the generalized complementary functions is the urgent work for the application of B2015 and H2017 for evaporation estimation, as well as the development of the generalized complementary principle. Given the 355 variations in  , the AA, B2015 and H2012 all have two parameters. The linear AA with a default value of b=1 has achieved a great success in evaporation estimation. For the B2015, c was thought to be only applied to accommodate unusual situations (Brutsaert, 2015). In practice, c = 0 is adopted and the Priestley-Taylor coefficient is calibrated Brutsaert et al., 2017;Liu et al., 2016;Brutsaert, 2015). But the calibrated  is smaller than the widely accepted constant 1.26 or even smaller than the unit at several sites, which is physically unrealistic. Han and Tian (2018a)  to b=4.16, not b=1. The consistency suggests that c needs to be calibrated. By calibrating both  and c, the B2015 performed well in estimating evaporation for 20 FLUXNET sites, and the value of  were more rational (Han and Tian, 2018a).
By contrast, two more parameters ( min x and max x ) are added to H2017. Because the sigmoid complementary curve are insensitive to min x and max x , Han and Tian (2018a) suggested that they could be treated them as constant parameters for 365 application convenience. min x and max x may change along with rad E , and were thought to vary with the time scales (Han and Tian, 2018a), min 0 x  and max 1 x  are appropriate at a daily scale for convenience, as have be evidenced by the well performances when compared to the flux measurements (Han and Tian, 2018a;Han et al., 2012). min x and max x are expected to be calculated by applying certain approaches to reduce the number of parameters of H2017 to two .
Although  would vary in theory (Assouline et al., 2016), it is widely used with a constant value of 1.26 in practice 370 (Priestley and Taylor, 1972). After calibrating, the variations of  is much less significant than those of the other parameters.
Moreover, the calibrated  approaches 1.26, especially for the H2017. Thus, the constant 1.26

 
was suggested with acceptable weakening of the accuracy of E estimation (Han and Tian, 2018a;Han et al., 2012). In practice,  was also determined from the observed E values in wet condition when E is close to Pen E and/or PT E (Kahler and Brutsaert, 2006;Ma et al., 2015a;Wang et al., 2019). A novel method by using observed air temperature and humidity data under wet environment 375 was proposed by Szilagyi et al. (2017) when measured E is lacking, and was successfully used for large-scale CR model applications (Ma and Szilagyi, 2019;Ma et al., 2019).
After determining  in advance, only a single parameter in the generalized complementary functions needs to be calibrated. As the parameters of the B2015 and H2017 can be transferred from the asymmetric parameter b of the original CR (Han and Tian, 2018a), the former studies on the characteristics of b could help its parameterization. The b in the desert was 380 much smaller than those in the oases or irrigated farmlands (Han et al., , 2012. b was thought to be related to the characteristics of the atmosphere, i.e., the atmospheric humidity (Szilagyi, 2015), the Clausius-Clapeyron relationship between saturation-specific humidity and temperature (Lintner et al., 2015), or the characteristics of the land surface, i.e., the surface temperature (Szilagyi, 2007), the water availability of the evaporating surface (Han and Tian, 2018b;Lhomme and Guilioni, 2010), or the ecosystem types (Wang et al., 2019). Szilagyi (2015) applied a sigmoid function of relative humidity to 385 parameterize b -1 . Wang et al. (2019) used the ecosystem mean b values of 217 sites around the world in the B2017 with litter weakening of the evaporation estimation accuracy. However, the characteristics and determination methods of b need further studies toward a calibration-free evaporation estimation model.

Applications of generalized complementary principle
The generalized complementary functions have been validated or applied in evaporation estimation for many sites (Ai 390 et al., 2017;Brutsaert et al., 2017;Zhang et al., 2017;Han and Tian, 2018a;Crago and Qualls, 2018), and several basins in China (Liu et al., 2016;Gao et al., 2018). It should be noted that most, if not all, above mentioned CR applications need "prior" knowledge in E (either ground-measured or water-balance-derived) to calibrate the parameters. Recently, the calibration-free CR model of S2017 was applied for monthly evaporation estimation across the conterminous China (Ma et al., 2019) and United States (Ma and Szilagyi, 2019). A wide range of model evaluations against the plot-scale flux measurements and basin-395 scale water balance results suggested that the generalized complementary functions could serve as a benchmark tool for validating the large-scale E results simulated by those Land Surface models and Remote Sensing models (Ma and Szilagyi, 2019). However, further applications over the world are still needed to develop more long-term, high-resolution E datasets for use in hydrological and atmospheric communities. Morton (1983) thought that the ability of the complementary principle to estimate actual evaporation by using 400 meteorological variables only can significantly influence the science and practice of hydrology. However, the attempts in using the complementary principle for hydrological modelling (Oudin et al., 2005;Barr et al., 1997;Nandagiri, 2007) have been suspended, while those attempts in applying such principle in drought assessment (Kim and Rhee, 2016;Hobbins et al., 2016) are still in their infancy. Moreover, the potential applications in agriculture water management are limited in the sense that the irrigation-induced changes in potential evaporation at an annual timescale (Ozdogan et al., 2006;Han et al., 2014a;Han et al., 405 2017). Apparently, the complementary principle did not develop to its full capacity via the linear CR, which leaves a broad space for further applications of the generalized complementary functions.

Conclusion
The complementary principle conceptualizes the feedbacks of surface evaporation on potential evaporation and offers advantages in evaporation research. In the CR, both E and pa E are thought to converge to po E . An elaborated prognostic 410 simulation of po E is crucial in studying the CR (Parlange and Katul, 1992a;Szilagyi and Jozsa, 2008). Several efforts have attempted to retain the linear CR by rationally formulating pa E and/or po E and by employing an asymmetric parameter.
Inheriting the concepts of three types of evaporation, the linear CR has evolved into the generalized nonlinear function B15, ()