A review of the complementary principle of evaporation: from the original linear relationship to generalized nonlinear functions

The complementary principle is an important methodology for estimating actual evaporation by using routinely observed meteorological variables. This review summaries its 56-year development, focusing on how related studies have shifted from adopting a symmetric linear complementary relationship (CR) to employing generalized nonlinear functions. The original CR denotes that the actual evaporation (E) and “apparent” potential evaporation ($E_{p_{\mathrm{a}}}$) depart from the potential evaporation ($E_{p_{\mathrm{0}}}$) complementarily when the land surface dries from a completely wet environment with constant available energy. The CR was then extended to an asymmetric linear relationship, and the linear nature was retained through properly formulating $E_{p_{\mathrm{a}}}$ and/or $E_{p_{\mathrm{0}}}$. Recently, the linear CR was generalized to a sigmoid function and a polynomial function. The sigmoid function does not involve the formulations of $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ but uses the Penman (1948) potential evaporation and its radiation component as inputs, whereas the polynomial function inherits $E_{p_{\mathrm{0}}}$ and $E_{p_{\mathrm{a}}}$ as inputs and requires proper formulations for application. The generalized complementary principle has a more rigorous physical base and offers a great potential in advancing evaporation estimation. Future studies may cover several topics, including the boundary conditions in wet environments, the parameterization and application over different regions of the world, and integration with other approaches for further development.

1 Introduction

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; McMahon et al., 2016; Morton, 1983). In this review paper, the terms “evaporation” and “evapotranspiration” are considered equivalent. As its underlying physical basis, this principle originated from the negative 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, 1983). By contrast, the Penman approach neglects the abovementioned feedback (Morton, 1983) and relies heavily on land surface variables such as soil moisture content or stomatal resistance (Monteith, 1965; Penman, 1950). Compared to the Budyko framework (Budyko, 1974; Turc, 1954), which is often used for partitioning catchment evaporation from precipitation at mean annual or annual timescales, actual evaporation from the landscape can be estimated within the complementary framework from nearly hourly to annual timescales. However, as a less popular hydro-climatic framework to estimate evaporation, the complementary principle needs to be advanced not only for its own development but also for the integration with other approaches for further development of evaporation research.

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, the apparent potential evaporation ($E_{p_{\mathrm{a}}}$) of a small saturated surface inside the landscape and the potential evaporation ($E_{p_{\mathrm{0}}}$) that occurs from the same large-size surface of E when it is saturated. In practice, $E_{p_{\mathrm{a}}}$ corresponds to current atmosphere in contact with the unsaturated evaporating surface as the overpassing air is not affected by the small saturated surface, whereas the atmosphere corresponding to $E_{p_{\mathrm{0}}}$ is in contact with the “potential” saturated surface. Thus, the surface water availability can be detected from the relative magnitude of $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ because of the land surface–atmosphere interaction, and E can be estimated without the explicit knowledge of the surface. The original symmetric linear “complementary” relationship (Bouchet, 1963; Brutsaert and Stricker, 1979; Morton, 1983) evolved into an asymmetric linear relationship (Brutsaert and Parlange, 1998; Pettijohn and Salvucci, 2006; Szilagyi, 2007). However, its development and applications are hindered by the use of complex formulations of $E_{p_{\mathrm{0}}}$ and $E_{p_{\mathrm{a}}}$ to retain the linear CR, which will be reviewed in more detail in the following sections.

Recent studies have adopted the “generalized” complementary principle, which employs nonlinear functions instead of the linear CR (Brutsaert, 2015; Han et al., 2012; Han and Tian, 2018a). The generalized complementary function comes in two ways. The first attempt abandons the concept of $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ yet uses a sigmoid function to describe the relationship among E, Penman's potential evaporation (E_Pen) and its radiation component (E_rad) (Han et al., 2012; Han and Tian, 2018a). By contrast, the other attempt adopts a polynomial function to describe the relationship between E, $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$. However, $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ still need to be formulated before applying the polynomial function to practical problems (Brutsaert, 2015).

The generalized complementary principle with earlier linear CRs as special cases has a more rigorous physical base (Brutsaert, 2015; Han and Tian, 2018b), and its methodology based on nonlinear functions is robust and effective. The generalized complementary principle has received much attention for its promising applications in estimating evaporation upon its proposal (Ai et al., 2017; Brutsaert et al., 2017, 2020; Han and Tian, 2018a; Liu et al., 2016; Szilagyi et al., 2016; Zhang et al., 2017). However, the boundary conditions and proper mathematical forms of the generalized complementary functions are still under study (Crago et al., 2016; Han and Tian, 2018a; Ma and Zhang, 2017; Szilagyi et al., 2017). In this review, we summarize the 56-year development of the complementary principle with a specific focus on its evolution from a symmetric linear CR to generalized nonlinear functions. We also compare the two types of generalized complementary functions and discuss their future development.

2 Linear complementary relationship

2.1 Concept of the symmetric complementary relationship

The concept of CR is illustrated in Fig. 1. When the water availability of the landscape is not limited, E is assumed to proceed at $E_{p_{\mathrm{a}}}$ and $E=E_{p_{\mathrm{a}}}=E_{p_{\mathrm{0}}}$. Given that the surface dries with constant available energy, E and $E_{p_{\mathrm{a}}}$ depart from $E_{p_{\mathrm{0}}}$ with equal yet opposite changes in fluxes and exhibit a complementary relation as follows:

$$\begin{array}{}\text{(1)}& E_{p_{\mathrm{a}}}-E_{p_{\mathrm{0}}}=E_{p_{\mathrm{0}}}-E.\end{array}$$

The formulations of $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ should be specified in Eq. (1). Bouchet (1963) assumed $E_{p_{\mathrm{0}}}$ to be half the input solar radiation. Morton (1976) calculated $E_{p_{\mathrm{a}}}$ by using the modified Penman (1948) equation proposed by Kohler and Parmele (1967) ($E_{\mathrm{Pen}}^{\mathrm{KP}})$, in which a constant vapor transfer coefficient was used to replace the wind function, and calculated $E_{p_{\mathrm{0}}}$ by using the Priestley and Taylor (1972) equation (E_PT) for an extensive saturated surfaced with minimal advection. This method has been used to calculate monthly evaporation in large areas.

Table 1 Formulations of $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ in different linear CR formulations.

^* Please refer to the main text and the Appendix for details on the symbols.

Figure 1 Schematic illustration of symmetric CR following Bouchet (1963).

Brutsaert and Stricker (1979) proposed the advection–aridity (AA) approach at a daily timescale, where $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ are directly formulated by E_Pen and E_PT, respectively. Although various combinations of $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ exist (Table 1), $E_{p_{\mathrm{0}}}$ is widely accepted to reflect the energy input while $E_{p_{\mathrm{a}}}$ includes the drying power of air simultaneously (Bouchet, 1963; Lhomme and Guilioni, 2006; Morton, 1983). Therefore, the AA approach seems logical and convincing (Lhomme and Guilioni, 2006). This approach has been validated based on hourly (Crago and Crowley, 2005; Parlange and Katul, 1992a), daily (Ali and Mawdsley, 1987; Brutsaert and Stricker, 1979; Qualls and Gultekin, 1997), monthly (Hobbins et al., 2001; Lemeur and Zhang, 1990; Xu and Singh, 2005) and annual (Ramirez et al., 2005; Yu et al., 2009) data from either plot-scale lysimeters and eddy covariance measurements or basin-wide water balance-derived results. By calculating E_Pen and E_PT using the standard meteorological data, the AA approach has been applied to estimate evaporation in regions with various land cover and climatic features (Hobbins et al., 2001; Liu et al., 2006; Ozdogan and Salvucci, 2004; Wang et al., 2011). For instance, this approach has been applied and validated in China from the Gobi Desert with mean annual precipitation of less than 150 mm (Han et al., 2008; Lemeur and Zhang, 1990; Liu and Kotoda, 1998) to humid eastern China with annual precipitation of approximately 1800 mm (Xu and Singh, 2005). Note, however, that the AA approach tends to overestimate E in wet environments but underestimate E in arid environments. Measurement error, imperfect formulations of $E_{p_{\mathrm{a}}}$ and/or $E_{p_{\mathrm{0}}}$, external energy sources and even the nonlinear nature of the complementary principle were considered as potential causes of this bias (Han et al., 2008, 2012; Hobbins et al., 2001; Qualls and Gultekin, 1997).

2.2 Proof of the complementary relationship

Bouchet (1963) and Morton (1965, 1970) approximately validated the CR by using annual and monthly data, respectively. At an annual scale, E and $E_{p_{\mathrm{a}}}$ (which are represented by E_Pen or pan evaporation E_Pan) were plotted against annual precipitation, and their negative relationship was used as a piece of evidence to support the reliable probability of the complementarity (Morton, 1983). Ramirez et al. (2005) tested the CR by using a composite of 192 data pairs from 25 basins across the USA and claimed a direct observational evidence. Yu et al. (2009) examined the CR at 102 observatories across China and found the CR at low elevations. Su et al. (2015) also showed a negative correlation between E from atmospheric reanalysis data and E_Pan in the non-humid regions of China. The large-scale irrigation development in an arid environment provides a large “natural” experimental area for validating the CR by the opposite changes in E and $E_{p_{\mathrm{a}}}$ (Roderick et al., 2009). A study from Turkey revealed that the warm-season $E_{p_{\mathrm{a}}}$ decreased progressively along with an increasing irrigated area (Ozdogan and Salvucci, 2004). Similar results were obtained from arid irrigation districts in northwest China, where an increasing irrigation water consumption reduces $E_{p_{\mathrm{a}}}$ (Han et al., 2014d) whereas a decreasing irrigation water consumption increases $E_{p_{\mathrm{a}}}$ (Han et al., 2017). However, although these studies showed that E and $E_{p_{\mathrm{a}}}$ move in opposite directions in most cases, there was not solid evidence to support the symmetric nature of CR.

The plausibility of CR has also been validated on theoretical bases and has been mathematically rationalized by Bouchet (1963), Morton (1969, 1971) and Seguin (1975). The rationalization proposed by Morton (1969, 1971) considers the governing changes in the humidity and temperature of the equilibrium sublayer of the atmospheric boundary layer (ABL) by assuming that (1) the net radiation will not change with the surface and (2) the heat and vapor eddy transfer characteristics are identical for E and $E_{p_{\mathrm{a}}}$. Relaxing the second assumption of Morton (1983), Szilagyi (2001) derived the CR by using the mass conservation equation for water vapor. However, LeDrew (1979) argued that Morton's two assumption do not necessarily hold and pointed out that the symmetric CR is physically unrealistic by using a diagnostic model of the energy fluxes within a closed system.

The physical basis of the CR has been further explored by using climate models. McNaughton and Spriggs (1989) tested the CR by using a simple model of the atmospheric mixed layer with entrainment in which the latent heat of the surface is simulated by using the bulk mass transfer equation with bulk resistance. During the validation, $E_{p_{\mathrm{a}}}$ is calculated via Penman's equation, which uses the temperature and humidity obtained from the results of the mixed-layer model corresponding to certain resistance ($E_{\mathrm{Pen}}^{r_{\mathrm{s}}}$), while $E_{p_{\mathrm{0}}}$ is calculated with the surface resistance set to zero ($E_{\mathrm{Pen}}^{r_{\mathrm{s}}=\mathrm{0}}$). 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 (1989). By using the Penman–Monteith equation to govern the areal latent heat flux at the surface, Lhomme (1997b) proposed a closed-box model with an impermeable lid at a fixed height, while Lhomme (1997a) used a more 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's (1997a) model, which was calibrated by using a data set obtained from the Hexi Corridor desert area in northwest China. But a strict symmetric CR was hardly confirmed by these studies.

2.3 Asymmetric linear complementary relationship

With $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ denoted by the mass-transfer-type potential evaporation E_MT and E_Pen, respectively, Granger (1989) proposed an alternative CR as follows:

$$\begin{array}{}\text{(2)}& \left(E_{\mathrm{MT}}-E_{\mathrm{Pen}}\right)={\displaystyle \frac{\mathrm{\Delta }\left(T_{\mathrm{a}}\right)}{\mathit{\gamma }}}\left(E_{\mathrm{Pen}}-E\right),\end{array}$$

where γ is the psychrometric constant, Δ(T_a) is the slope of the saturation vapor pressure at air temperature T_a. Despite being identical to the surface energy balance, Eq. (2) has inspired researchers to examine whether the CR should be symmetric or not (Pettijohn and Salvucci, 2006; Szilagyi, 2007). By using pan evaporation to denote $E_{p_{\mathrm{a}}}$, Brutsaert and Parlange (1998) extended the symmetric CR as follows:

$$\begin{array}{}\text{(3)}& \left(E_{p_{\mathrm{a}}}-E_{p_{\mathrm{0}}}\right)=b\left(E_{p_{\mathrm{0}}}-E\right),\end{array}$$

where b is the coefficient that denotes asymmetry and the original symmetric CR is characterized with b=1. Kahler and Brutsaert (2006) clarified and tested Eq. (3) at a daily timescale and attributed the asymmetry to the nature of the heat transfer between the pan and its surroundings, which made the changes in pan evaporation larger than those in E. Szilagyi (2007) showed that the asymmetry is not limited only to the evaporation pan but is also inherently linked to the definition of $E_{p_{\mathrm{a}}}$. Brutsaert (2015) stated that asymmetry is an inherent characteristic of the CR.

The asymmetric CR can be illustrated in a dimensionless form (Fig. 2) (Kahler and Brutsaert, 2006). Normalized by $E_{p_{\mathrm{0}}}$, $E_{p_{\mathrm{a}}}$ and E can be scaled as

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{E}{E_{p_{\mathrm{0}}}}}={\displaystyle \frac{(\mathrm{1}+b)E/E_{p\mathrm{a}}}{\mathrm{1}+bE/E_{p_{\mathrm{a}}}}}\mathrm{and}{\displaystyle \frac{E_{p_{\mathrm{a}}}}{E_{p_{\mathrm{0}}}}}={\displaystyle \frac{\mathrm{1}+b}{\mathrm{1}+bE/E_{p_{\mathrm{a}}}}}.\end{array}$$

The scaled $E_{p_{\mathrm{a}}}$ and E are both functions of the dimensionless variable $E/E_{p_{\mathrm{a}}}$, while $E/E_{p_{\mathrm{a}}}$ serves as the evaporative surface moisture index. Compared with the original form (Eq. 1 and Fig. 1), the CR here is illustrated without the presence of the water availability explicitly. The asymmetric CR has been validated via the opposite changes of $E/E_{p_{\mathrm{0}}}$ and $E_{p_{\mathrm{a}}}/E_{p_{\mathrm{0}}}$ against $E/E_{p_{\mathrm{a}}}$ at several locations across the world. However, the wet conditions were seldom explored, which may hide the true correlation as the two curves of $E/E_{p_{\mathrm{0}}}$ and $E_{p_{\mathrm{a}}}/E_{p_{\mathrm{0}}}$ approach one another.

Figure 2 Scaled $E_{p_{\mathrm{a}}}$ and E, which serve as functions of the evaporative moisture index $E/E_{p_{\mathrm{a}}}$ and are calculated on the basis of the asymmetric CR, according to method of Kahler and Brutsaert (2006).

The asymmetric CR is a significant improvement of the symmetric CR, and the opposite changes of $E/E_{p_{\mathrm{0}}}$ and $E_{p_{\mathrm{a}}}/E_{p_{\mathrm{0}}}$ against $E/E_{p_{\mathrm{a}}}$ were treated as an enhanced illustration of the CR (Brutsaert et al., 2020; Hu et al., 2018; Ma et al., 2015a; Szilagyi, 2007; Zhang et al., 2017). The performances on evaporation estimation are improved by calibrating the asymmetry parameter b (Han et al., 2008; Huntington et al., 2011; Kahler and Brutsaert, 2006; Ma et al., 2015a). Efforts have also made to calculate b by using the meteorological variables, which enhance the predictability of the CR (Aminzadeh et al., 2016; Szilagyi, 2007, 2015). However, the changes in b imply a potential nonlinear characteristic of the CR (Han, 2008; Lintner et al., 2015). The observed values of $E/E_{p_{\mathrm{0}}}$ and $E_{p_{\mathrm{a}}}/E_{p_{\mathrm{0}}}$ can even exhibit a positive correlation under wet conditions at several flux sites, which challenges the linear CR (Han and Tian, 2018a).

2.4 Efforts to retain the linear nature of the complementary relationship through properly formulating $E_{p_{\mathrm{a}}}$ and/or $E_{p_{\mathrm{0}}}$

The imperfect linear CR has inspired researchers to apply rational formulations of $E_{p_{\mathrm{a}}}$ and/or $E_{p_{\mathrm{0}}}$ to retain it. One direct method is to revise the formulations of E_Pen and/or E_PT based on the AA approach through calibration. For E_Pen, the empirical wind function was calibrated to improve the CR (Hobbins et al., 2001). However, Penman's wind function cannot work under the wet and dry conditions simultaneously (Pettijohn and Salvucci, 2006). The wind function derived from the Monin–Obukhov similarity theory was then employed (Crago and Crowley, 2005; Ma et al., 2015a; Parlange and Katul, 1992b; Pettijohn and Salvucci, 2006). The surface roughness and surface albedo were also calibrated to improve the CR (Lemeur and Zhang, 1990). Meanwhile, for E_PT, the Priestley–Taylor coefficient (α) is regarded varying, thereby leaving a range for calibration (Han et al., 2006; Xu and Singh, 2005; Yang et al., 2012). In addition to E_Pen and E_PT, the mass-transfer-type potential evaporation (van Bavel, 1966) (E_MT) was considered as another formulation of $E_{p_{\mathrm{a}}}$ (Granger, 1989). Different combinations of $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$, (i.e., E_Pen, E_PT and E_MT) were tested through the trial-and-error method to retain the linear nature of CR (Anayah and Kaluarachchi, 2014; Crago and Crowley, 2005).

Given the conceptual inadequacy in using E_Pen and E_PT to denote $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ (Morton, 1983; Szilagyi and Jozsa, 2008), a better CR must be obtained by modifying the formulations of E_Pen and/or E_PT on a physical basis. For $E_{p_{\mathrm{a}}}$, the net long-wave radiation depends on the land surface temperature; meanwhile, adjusting surface temperature with air temperature to calculate solar radiation in E_Pen may be problematic (Morton, 1983). To address these limitations, Morton (1983) combined the energy balance and water vapor transfer equations by using an equilibrium temperature (T_p) and derived a Morton-type potential evaporation E_Mor to denote $E_{p_{\mathrm{a}}}$. By attributing the asymmetry to the assumption that $E_{p_{\mathrm{a}}}$ conceptually includes a transpiration component, Pettijohn and Salvucci (2006) improved the asymmetry by replacing E_Pen with the Penman–Monteith equation with a minimum surface resistance ($E_{\mathrm{PM}}^{r_{s_{\min }}}$). Similarly, the reference evapotranspiration (ET₀) was also used to replace E_Pen (Han et al., 2017, 2014d).

In theory, $E_{p_{\mathrm{0}}}$ is the potential evaporation when the land surface is saturated and should be calculated with a proper formula by using meteorological variables corresponding to the “potential” saturated surface. The Priestley–Taylor equation has been widely accepted to represent evaporation from extensive saturated surfaces by using meteorological variables corresponding to these saturated surfaces (Brutsaert, 1982; Priestley and Taylor, 1972). This way it was suggested to represent $E_{p_{\mathrm{0}}}$ (Brutsaert and Stricker, 1979). However, in the AA approach, E_PT is calculated by the Priestley–Taylor equation using the atmospheric variables that correspond to the current unsaturated surface. But the atmosphere in contact with the land surface will change if the land surface is saturated (Brutsaert, 2015; Morton, 1983). Thus, E_PT is in reality a variable dependent on the meteorological variables at the time of calculation and does not represent the “true” $E_{p_{\mathrm{0}}}$.

Obviously, calculating the slope of the saturation vapor pressure at the current air temperature (Δ(T_a)) for E_PT is imperfect because the temperature corresponding to $E_{p_{\mathrm{0}}}$ is different from current T_a corresponding to an unsaturated environment (Morton, 1983; Szilagyi and Jozsa, 2008). Thus, predicting the air temperature corresponding to the extensive saturated surface is critical for properly formulating $E_{p_{\mathrm{0}}}$. Morton (1983) derived $E_{p_{\mathrm{0}}}$ by using a modified Priestley–Taylor equation with net radiation and the slope of the saturation vapor pressure that is calculated at equilibrium temperature T_p ($E_{\mathrm{PT}}^{T_{\mathrm{p}}}$). Szilagyi and Jozsa (2008) argued that Δ in E_PT should be calculated at the air temperature corresponding to the wet environment instead of actual air temperature, which is not straightforward to derive. Thus, Szilagyi and Jozsa (2008) proposed an iterative approach based on the Bowen ratio method to estimate the surface temperature in wet environments (T_ws) and replaced Δ(T_a) with the slope of the saturation vapor pressure curve at T_ws (Δ(T_ws)) in the Priestley–Taylor equation ($E_{\mathrm{PT}}^{T_{\mathrm{ws}}}$) by assuming a negligible temperature gradient over such a small wet area. $E_{\mathrm{PT}}^{T_{\mathrm{ws}}}$ was used to improve the symmetry of the CR in arid shrubland environments (Huntington et al., 2011) and in an alpine steppe of the Tibetan Plateau (Ma et al., 2015a). The evaporation estimations across the USA were also improved by applying the modified AA approach (Szilagyi, 2015; Szilagyi et al., 2009; Szilagyi and Jozsa, 2008). Aminzadeh et al. (2016) derived a steady-state surface temperature via the surface energy balance at which the sensible heat flux is zero, and calculated $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ using a mass-transfer-type reference evaporation corresponding to current and saturated surface water content.

Advection is another factor influencing $E_{p_{\mathrm{0}}}$, which could not be adequately considered by E_PT with an assumption of a minimal advection effect (Morton, 1975, 1983; Parlange and Katul, 1992a). The effects of advection were considered by an empirical correction factor in $E_{\mathrm{PT}}^{T_{\mathrm{p}}}$ (Morton, 1975, 1983). Parlange and Katul (1992a) attributed the asymmetry to the horizontal advection of dry air, which would make E_Pen larger than the available energy (R_n−G) (i.e., $R_{\mathrm{n}}-G-E_{\mathrm{Pen}}<\mathrm{0}$) and proposed to replace E_PT with $E_{\mathrm{PT}}+|R_{\mathrm{n}}-G-E_{\mathrm{Pen}}|$ to improve the CR on an hourly basis.

The efforts of reformulating $E_{p_{\mathrm{a}}}$ and/or $E_{p_{\mathrm{0}}}$ through the calibration, the trial-and-error process and the physical improvement have significantly improved the evaporation estimation (Hobbins et al., 2001; Ma et al., 2015a; Szilagyi, 2015; Xu and Singh, 2005). However, it is always impossible to find formulations of $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ completely rational at present, and these approaches are deemed ineffective because of their high computation demand, which is a key stumbling block when applying the CR on a large scale (e.g., continental or global) (Ma et al., 2019).

3 Generalized complementary principle via nonlinear functions

3.1 Normalized complementary functions

Unlike the normalization by $E_{p_{\mathrm{0}}}$ (Kahler and Brutsaert, 2006), Han (2008) normalized Eq. (3) by using $E_{p_{\mathrm{a}}}$ and found that $E/E_{p_{\mathrm{a}}}$ is expressed as a linear function of $E_{p_{\mathrm{0}}}/E_{p_{\mathrm{a}}}$. Normalized by E_Pen (Han et al., 2008), the AA approach can be expressed as

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{E}{E_{\mathrm{Pen}}}}=\mathit{\alpha }\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{b}}\right){\displaystyle \frac{E_{\mathrm{rad}}}{E_{\mathrm{Pen}}}}-{\displaystyle \frac{\mathrm{1}}{b}},\end{array}$$

where E∕E_Pen is regarded as a linear function of E_rad∕E_Pen. The bias of the AA function in arid and wet environments can be easily understood in its dimensionless form. Also, the AA approach with a tuned b still underestimated evaporation in arid environments (Han et al., 2008), which implies that the CR may deviate from its linear characteristics.

Based on the examination of the CR using a model of the convective boundary layer with entrainment (Lhomme, 1997a), Lhomme and Guilioni (2006, 2010) recommended a form of the CR through the effective surface resistance of the region. Integrating this relationship into the Penman–Monteith equation and the normalization by E_Pen lead to

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{E}{E_{\mathrm{Pen}}}}=(\mathrm{1}+\mathit{\omega }){\displaystyle \frac{E_{\mathrm{rad}}}{E_{\mathrm{Pen}}}},\end{array}$$

where ω is a coefficient accounting for the entrainment of dry air within the atmospheric boundary layer. Equation (6) is a linear function without intercept but was not verified and applied using observed data.

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, 1989; 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, 1989). This relationship was integrated to an asymmetric CR to improve the performance on evaporation estimation (Anayah and Kaluarachchi, 2014). Normalized by E_Pen, Granger's model is similar to the AA function in that E∕E_Pen is expressed as a function of the relative magnitude of drying power to net radiation (Han et al., 2011). By synthesizing the dimensionless forms of the AA function and the Granger model, Han et al. (2011) proposed the following function as an alternative:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{E}{E_{\mathrm{Pen}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+c_{\mathrm{1}}{e}^{d\left(\mathrm{1}-\frac{E_{\mathrm{rad}}}{E_{\mathrm{Pen}}}\right)}}},\end{array}$$

where c₁ and d are the parameters. Equation (7) approximates the linear AA function under normal conditions neither too wet nor too dry but amends its bias (Han et al., 2011); thus it can be regarded as an enhanced nonlinear version of the linear CR.

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 (Katerji and Perrier, 1983; Liu et al., 2012; Ma et al., 2015b; Rana et al., 1997). 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 E_Pen:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{E}{E_{\mathrm{Pen}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+k\left(\frac{E_{\mathrm{Pen}}}{E_{\mathrm{rad}}}-\mathrm{1}\right)+l}},\end{array}$$

where k and l are the empirical calibration parameters. With similar variables yet different mathematical formulations, Eq. (8) can also be considered a complementary function (Han et al., 2014c).

3.2 Sigmoid function relating E∕E_Pen to E_rad∕E_Pen

By synthesizing the aforementioned functions (Table 2), Han et al. (2012) generalized the CR as a function that relates E∕E_Pen to E_rad∕E_Pen:

$$\begin{array}{}\text{(9)}& E/E_{\mathrm{Pen}}=f\left(E_{\mathrm{rad}}/E_{\mathrm{Pen}}\right).\end{array}$$

Equation (8) shares the same form of Penman's approach for estimating evaporation. The function of surface wetness that denotes the reduction of E to E_Pen is replaced by the function of E_rad∕E_Pen, which is termed “atmospheric wetness” (Han and Tian, 2018b). Despite not explicitly exhibiting a CR, Eq. (9) holds the complementary principle that the land surface wetness is indirectly denoted by the drying power of air with a constant radiation energy input (Brutsaert, 1982). Accordingly, Eq. (9) is considered a “general form” of the CR (Han et al., 2014b) (hereinafter referred to as H12, whereas the other type of generalized complementary function first proposed by Brutsaert (2015) is referred to as B15 for comparison). The existing analytical forms of the function can be classified as linear, concave/convex or sigmoid (Table 2). Studies on the complementary principle can be advanced by formulating a proper analytical form for H12.

Table 2 Different analytical formulas for generalized complementary function H12.

^* $x=E_{\mathrm{rad}}/E_{\mathrm{Pen}}$ and $y=E/E_{\mathrm{Pen}}$; the other symbols are parameters.

The exact analytical form of H12 is inadequately understood at present. However, some of its characteristics can be detected from its boundary conditions in extremely arid and completely wet environments. Han et al. (2012) derived the zeroth-order and first-order boundary conditions for H12 as

$$\begin{array}{}\text{(10)}& \left\{\begin{array}{lcc}{y}_{\mathrm{H}}=\mathrm{0},x_H& \to & \mathrm{0}\\ y_{\mathrm{H}}=\mathrm{1},x_H& \to & \mathrm{1}\\ {\displaystyle \frac{\mathrm{d}{y}_{\mathrm{H}}}{\mathrm{d}{x}_{\mathrm{H}}}}=\mathrm{0},x_{\mathrm{H}}& \to & \mathrm{0}\\ {\displaystyle \frac{\mathrm{d}{y}_{\mathrm{H}}}{\mathrm{d}{x}_{\mathrm{H}}}}=\mathrm{0},x_{\mathrm{H}}& \to & \mathrm{1}\end{array}\right.,\end{array}$$

where $x_{\mathrm{H}}=E_{\mathrm{rad}}/E_{\mathrm{Pen}}$ and $y_{\mathrm{H}}=E/E_{\mathrm{Pen}}$. Han et al. (2012) proposed the following sigmoid function (hereinafter this specific analytical form of H12 is referred to as sigmoid H2012):

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{E}{E_{\mathrm{Pen}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+m{\left(\frac{E_{\mathrm{Pen}}}{E_{\mathrm{rad}}}-\mathrm{1}\right)}^n}},\end{array}$$

where m and n are parameters. 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 Eq. (7) (Han et al., 2012); the application of this sigmoid function has also been recommended for an alpine meadow region of the Tibetan Plateau (Ma et al., 2015b).

The zeroth-order arid boundary condition of H12 adopted in H2012 may be imperfect in the sense that the aerodynamic term (E_aero) of E_Pen may not reach infinity under an arbitrary E_rad (Crago et al., 2016; Kovács, 1987; Szilagyi et al., 2017). Moreover, E_rad∕E_Pen cannot easily approach unity because of advection (Kovács, 1987; Priestley and Taylor, 1972). Therefore, Han and Tian (2018a) brought in the minimum and maximum limits of E_rad∕E_Pen (x_min and x_max) under an assumed constant E_rad and re-derived the boundary conditions of H12 by adopting two widely accepted assumptions following Penman's combination theory, namely, $\partial E/\partial E_{\mathrm{Pen}}=\mathrm{0}$ in extremely arid environments and E=E_Pen in completely wet environments. The boundary conditions are set as follows:

$$\begin{array}{}\text{(12)}& \left\{\begin{array}{lcc}{y}_{\mathrm{H}}=\mathrm{0},x_{\mathrm{H}}& \to & x_{\min }\\ y_{\mathrm{H}}=\mathrm{1},x_{\mathrm{H}}& \to & x_{\max }\\ {\displaystyle \frac{\mathrm{d}{y}_{\mathrm{H}}}{\mathrm{d}{x}_{\mathrm{H}}}}=\mathrm{0},x_{\mathrm{H}}& \to & x_{\min }\\ {\displaystyle \frac{\mathrm{d}{y}_{\mathrm{H}}}{\mathrm{d}{x}_{\mathrm{H}}}}=\mathrm{0},x_{\mathrm{H}}& \to & x_{\max }\end{array}\right..\end{array}$$

Based on the boundary conditions, Han and Tian (2018a) speculated that the growth of E∕E_Pen upon E_rad∕E_Pen exhibits a sigmoid feature, which is a three-stage pattern in which E∕E_Pen gradually increases along with E_rad∕E_Pen, rapidly increases along with E_rad∕E_Pen in the following stage and then demonstrates a decelerated growth in the final stage. The sigmoid feature can be detected from the study by Han et al. (2012) at the arid Gobi site and the humid Wudaogou site in China. Han and Tian (2018a) further validated the sigmoid feature according to the much larger regression slopes of E∕E_Pen upon E_rad∕E_Pen in the middle stage than those in the other two stages with smaller or larger values of E_rad∕E_Pen by using 22 eddy covariance towers from the FLUXNET (Baldocchi et al., 2001) data set, which includes representative biomes of grasslands, croplands, shrublands, evergreen needleleaf forests, deciduous broadleaf forests and wetlands.

In 2017, Han and Tian (2018a) proposed the following new sigmoid function in accordance with the boundary conditions (hereinafter referred to as sigmoid H2017):

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{E}{E_{\mathrm{Pen}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+m{\left(\frac{x_{\max }-E_{\mathrm{rad}}/E_{\mathrm{Pen}}}{E_{\mathrm{rad}}/E_{\mathrm{Pen}}-x_{\min }}\right)}^n}},\end{array}$$

where E_rad∕E_Pen adopts the feasible domain (x_min, x_max), which is a subdomain of (0, 1). Both the linear AA function and sigmoid H2012 are special cases of sigmoid H2017. Han and Tian (2018a) performed a first-order Taylor expansion of Eq. (13) at the point where $E/E_{\mathrm{Pen}}=\mathrm{0.5}$ and set the linear equation equivalent to the linear AA function. Afterward, the parameters m and n of sigmoid H2017 can be transformed from the Priestley–Taylor coefficient α and parameter b of the AA function.

Han et al. (2008) was probably the first to plot the AA function as a linear in the state space (E_rad∕E_Pen, E∕E_Pen), in which the biases of the AA function in arid and wet environments can be understood easily. The analytical forms of the generalized complementary function H12 listed in Table 2 can be plotted as curves in a 2-D space (E_rad∕E_Pen, E∕E_Pen) (Fig. 3), and the upper limits of E_Pen and E_PT are illustrated as the curve of OMN. The sigmoid H2012 was compared to the linear AA in the state space (E_rad∕E_Pen, E∕E_Pen) to demonstrate its improvement (Han et al., 2012). Observed E∕E_Pen can be plotted against E_rad∕E_Pen and fitted by the analytical functions of H12 in the state space (E_rad∕E_Pen, E∕E_Pen), which is an obvious improvement compared to the schematic illustrations of the CR in Figs. 1 and 2.

Figure 3 Generalized complementary functions in the state space (E_rad∕E_Pen, E∕E_Pen): linear AA, polynomial B2015 and sigmoid H2017, with α=1.26 and $b=\mathrm{1}\cdot x_{\min }$, and x_max are set to 0 and 0.9, respectively, as revised from Han and Tian (2018a). OM is the edge at which E=E_PT, MN is the edge where E=E_Pen, M corresponds to the condition of the minimal advection evaporation where E_PT=E_Pen and N corresponds to the condition of the equilibrium evaporation where E_Pen=E_rad.

3.3 Polynomial function relating $E/E_{p_{\mathrm{a}}}$ to $E_{p_{\mathrm{0}}}/E_{p_{\mathrm{a}}}$

Inspired by Hanet al. (2012), Brutsaert (2015) reformulated another general dimensionless form of the CR, $E/E_{p_{\mathrm{a}}}=f(E_{p_{\mathrm{0}}}/E_{p_{\mathrm{a}}})$, and proposed its boundary conditions as follows:

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{lcc}{y}_{\mathrm{B}}=\mathrm{0},x_{\mathrm{B}}& \to & \mathrm{0}\\ y_{\mathrm{B}}=\mathrm{1},x_{\mathrm{B}}& \to & \mathrm{1}\\ {\displaystyle \frac{\mathrm{d}{y}_{\mathrm{B}}}{\mathrm{d}{x}_{\mathrm{B}}}}=\mathrm{0},x_{\mathrm{B}}& \to & \mathrm{0}\\ {\displaystyle \frac{\mathrm{d}{y}_{\mathrm{B}}}{\mathrm{d}{x}_{\mathrm{B}}}}=\mathrm{1},x_{\mathrm{B}}& \to & \mathrm{1}\end{array}\right.,\end{array}$$

where $x_{\mathrm{B}}=E_{p_{\mathrm{0}}}/E_{p_{\mathrm{a}}}$ and $y_{\mathrm{B}}=E/E_{p_{\mathrm{a}}}$. The following fourth-order polynomial function was also derived to satisfy the boundary conditions:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{E}{E_{p_{\mathrm{a}}}}}=(\mathrm{2}-c){\left({\displaystyle \frac{E_{p_{\mathrm{0}}}}{E_{p_{\mathrm{a}}}}}\right)}^{\mathrm{2}}-(\mathrm{1}-\mathrm{2}c){\left({\displaystyle \frac{E_{p_{\mathrm{0}}}}{E_{p_{\mathrm{a}}}}}\right)}^{\mathrm{3}}-c{\left({\displaystyle \frac{E_{p_{\mathrm{0}}}}{E_{p_{\mathrm{a}}}}}\right)}^{\mathrm{4}},\end{array}$$

where c is a parameter. Brutsaert (2015) regarded Eq. (15) (hereinafter referred to as B15) as a generalization of the linear CR and referred to the corresponding methodology as the “generalized complementary principle”.

Table 3 Different forms of the generalized complementary function, y=f(x).

^a H12 denotes the t generalized complementary function $E/E_{\mathrm{Pen}}=f(E_{\mathrm{rad}}/E_{\mathrm{Pen}}$), while H2012 and H2017 are the sigmoid analytical formulas for H12. ^b B15 denotes the generalized complementary function $E/E_{p_{\mathrm{a}}}=f(E_{p_{\mathrm{0}}}/E_{p_{\mathrm{a}}})$, while B2015, C2016 and S2017 are the analytical formulas for B15 in application.

The application of Eq. (15) depends on specific formulations of $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$. In the manner of the AA approach, Eq. (15) has been applied to estimate evaporation (Ai et al., 2017; Brutsaert et al., 2017; Liu et al., 2016; Szilagyi et al., 2016; Zhang et al., 2017). In this case, we refer to Eq. (15) in the manner of the AA approach as B2015 to avoid confusion. Although estimating $E_{p_{\mathrm{a}}}$ by using E_Pen is widely accepted by the research community, prognostically predicting $E_{p_{\mathrm{0}}}$ based on E_PT remains a huge challenge considering the theoretical problems of the Priestley–Taylor coefficient. In addition, the lower limit of x_B→0 of B15 may not hold in the manner of the AA approach (Crago et al., 2016; Han and Tian, 2018a; Kovács, 1987; Szilagyi et al., 2017). To address these challenges, Crago et al. (2016) and Szilagyi et al. (2017) used the maximum value of $E_{p_{\mathrm{a}}}$ to rescale x_B and replaced E_PT with $E_{\mathrm{PT}}^{T_{\mathrm{ws}}}$; the latter is based on the air temperature in a wet environment. Crago et al. (2016) applied a mass transfer approach to calculate the maximum value of $E_{p_{\mathrm{a}}}$ ($E_{\mathrm{MT}}^{\max }$) and rescaled x_B as

$$\begin{array}{}\text{(16)}& x_{\mathrm{C}}={\displaystyle \frac{E_{\mathrm{PT}}^{T_{\mathrm{ws}}}/E_{\mathrm{Pen}}-E_{\mathrm{PT}}^{T_{\mathrm{ws}}}/E_{\mathrm{MT}}^{\max }}{\mathrm{1}-E_{\mathrm{PT}}^{T_{\mathrm{ws}}}/E_{\mathrm{MT}}^{\max }}}.\end{array}$$

Szilagyi et al. (2017) employed the Penman equation to calculate the maximum value of $E_{p_{\mathrm{a}}}$ ($E_{\mathrm{Pen}}^{\max }$) and proposed the following rescaled version:

$$\begin{array}{}\text{(17)}& x_{\mathrm{S}}={\displaystyle \frac{E_{\mathrm{Pen}}^{\max }-E_{\mathrm{Pen}}}{E_{\mathrm{Pen}}^{\max }-E_{\mathrm{PT}}^{T_{\mathrm{ws}}}}}{\displaystyle \frac{E_{\mathrm{PT}}^{T_{\mathrm{ws}}}}{E_{\mathrm{Pen}}}}.\end{array}$$

x_C and x_S are essentially same (Szilagyi et al., 2017) except for the different formulations for the maximum value of $E_{p_{\mathrm{a}}}$. However, $E_{\mathrm{MT}}^{\max }$ in Eq. (16) may became invalid under conditions with relatively strong available energy yet weak winds (Ma and Zhang, 2017) and was replaced with $E_{\mathrm{Pen}}^{\max }$ (Crago and Qualls, 2018) in the latest version. After rescaling, Crago et al. (2016) proposed a new linear version of the generalized complementary function (hereinafter referred to as C2016) (i.e., y_B=x_C; Table 2), while Szilagyi et al. (2017) used the third-order polynomial function (hereinafter referred to as 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 (Crago and Qualls, 2018; Crago et al., 2016; Szilagyi et al., 2017).

3.4 Comparisons between the two generalized complementary approaches

The two generalized complementary approaches, H12 and B15, are essentially different, with completely different normalized variables (Table 3). The differences in the analytical forms, sigmoid and fourth-order polynomial, mainly result from their wet boundary conditions. B15 inherits the concept of the three types of evaporation dated from the original CR, and its boundary conditions and analytical form are derived for $x_{\mathrm{B}}=E_{p_{\mathrm{0}}}/E_{p_{\mathrm{a}}}$ and $y_{\mathrm{B}}=E/E_{p_{\mathrm{a}}}$. The original CR adopts the limits of $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ on E in a serial manner ($E\le E_{p_{\mathrm{0}}}\le E_{p_{\mathrm{a}}}$) (Brutsaert, 2015) while considering that the wet regional evaporation must always be smaller than the wet patch evaporation ($E_{p_{\mathrm{0}}}\le E_{p_{\mathrm{a}}}$). Under wet conditions, B15 adopts $\mathrm{d}{y}_{\mathrm{B}}/\mathrm{d}{x}_{\mathrm{B}}=\mathrm{1}$ as x_B→1 by considering that any change in E is the same as the change in $E_{p_{\mathrm{0}}}$, which results in a concave polynomial type function. The limits and boundary conditions of B15 would be appropriate in theory. However, $E_{p_{\mathrm{0}}}$ and $E_{p_{\mathrm{a}}}$ should be formulated before B15 is applied to practical problems. Thus, B15 still faces one of the difficulties that the original CR has, that is, appropriately formulating $E_{p_{\mathrm{0}}}$ and $E_{p_{\mathrm{a}}}$, which determines the validity and application of B15. So future studies can be conducted towards more proper formulations of $E_{p_{\mathrm{a}}}$ and $E_{p_{\mathrm{0}}}$ to satisfy the boundary conditions of B15.

By contrast, H12 goes further from the original CR. The boundary conditions and the analytical form of H12 are derived for $x_{\mathrm{H}}=E_{\mathrm{rad}}/E_{\mathrm{Pen}}$ and $y_{\mathrm{H}}=E/E_{\mathrm{Pen}}$. The knowledge on $E_{p_{\mathrm{0}}}$ is unnecessary, and only the mostly accepted E_Pen and its radiation term appear in H12. Thus, the corresponding theoretical and practical difficulties of formulating $E_{p_{\mathrm{0}}}$ and $E_{p_{\mathrm{a}}}$ are eliminated. H12 adopts E_Pen as the upper limit E≤E_Pen during the derivation and introduce the limit of E≤E_PT by considering that E_PT is widely used as an upper limit of E in practice. Han and Tian (2018a) showed that the upper limits of E_Pen and E_PT on evaporation must be in parallel, that is, $\left\{\begin{array}{l}E\le E_{\mathrm{Pen}}\\ E\le E_{\mathrm{PT}}\end{array}\right.$, and the complementary curves should be constrained by the limits of OMN as illustrated in Fig. 3. The limits of E_Pen and E_PT on E can be approximately satisfied by the sigmoid function H2017 with the parameters transformed from the linear AA function (Han and Tian, 2018a). Furthermore, H12 adopts $\mathrm{d}{y}_{\mathrm{H}}/\mathrm{d}{x}_{\mathrm{H}}=\mathrm{0}$ as y_H→1 by considering that E approaches E_Pen under wet conditions, which results in a sigmoid type function.

In the manner of the AA approach of formulating $E_{p_{\mathrm{0}}}$ and $E_{p_{\mathrm{a}}}$, B15 evolves to one of its analytical forms, the polynomial B2015. Taking the Priestley–Taylor coefficient as a parameter, B2015 can also be regarded as a polynomial analytical form of H12 (Table 2) and can be compared with the sigmoid H2017 in the state space (E_rad∕E_Pen, E∕E_Pen) (Fig. 3). In the polynomial B2015, the limits on the E are specified to $E\le \mathit{\alpha }{E}_{\mathrm{rad}}\le E_{\mathrm{Pen}}$. In practice, a constant α is widely used, and the polynomial curves of B2015 are required to be constrained by the triangle domain OMP (Fig. 3), which discards the domain out of OMP. However, the Priestley–Taylor coefficient varies with several factors, such as the relative transport efficiency of turbulence or the surface/air temperature (Assouline et al., 2016; Szilagyi, 2014). Thus, E_rad∕E_Pen may be larger than 1∕α, revealing that the trapezoidal domain adopted by the sigmoid H2017 may be more accurate. In the state space (E_rad∕E_Pen, E∕E_Pen), the curve of the sigmoid H2017 exhibits a three-stage pattern, whereas the linear AA 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 is uncommon, the polynomial B2015 performs well with calibrated parameters (Brutsaert et al., 2017; Han and Tian, 2018a; Liu et al., 2016; Zhang et al., 2017). However, observed points are located in the domain out of OMP at several flux sites, and the sigmoid H2017 shows the best performance in estimating evaporation as validated by using data from FLUXNET (Han and Tian, 2018a; Wang et al., 2019).

4 Current applications and future developments of the generalized complementary principle

4.1 Current applications of the generalized complementary functions for evaporation estimation

Morton (1983) thought that the ability of the complementary principle to estimate actual evaporation by using meteorological variables only can significantly influence the science and practice of hydrology. However, attempts at using the complementary principle for evaporation estimation in hydrological modeling (Barr et al., 1997; Nandagiri, 2007; Oudin et al., 2005) have been suspended, while attempts at applying the principle in drought assessment (Hobbins et al., 2016; Kim and Rhee, 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 have mainly been evaluated at an annual timescale (Han et al., 2014a, 2017; Ozdogan et al., 2006). Apparently, the complementary principle did not develop to its full capacity via the linear CR, which leaves a broad space for applying the generalized complementary functions for evaporation research.

For example, the generalized complementary functions have been validated or applied in evaporation estimation for many sites (Ai et al., 2017; Brutsaert et al., 2017; Crago and Qualls, 2018; Han and Tian, 2018a; Zhang et al., 2017) and several basins in China (Gao et al., 2018; Liu et al., 2016). B2015 was applied to estimate global terrestrial evaporation with calibrated α as a function of the aridity index (Brutsaert et al., 2020). The modified Granger model was also applied for estimating global evaporation with 0.5^∘ spatial resolution and monthly time steps (Anayah and Kaluarachchi, 2019). It should be noted that most, if not all, abovementioned CR applications need “prior” knowledge in E (either ground measured or water balance derived) to calibrate the parameters. Recently, the Szilagyi et al. (2017) model was applied for monthly evaporation estimation without calibration across China (Ma et al., 2019) and the conterminous USA (Ma and Szilagyi, 2019). A wide range of model evaluations against the plot-scale flux measurements and basin-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 across the world are still needed to develop more long-term, high-resolution E data sets for use in hydrological and atmospheric communities.

4.2 Parameterizing generalized complementary functions for future applications

Determining the parameters of the generalized complementary functions is urgent for the application of B2015 and H2017 to evaporation estimation, as well as the development of the generalized complementary principle. Given the variations in α, the linear AA, polynomial B2015 and sigmoid H2012 all have two parameters. The linear AA with a default value of b=1 was applied at first 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, 2015; Brutsaert et al., 2017; Liu et al., 2016; Zhang et al., 2017). But the calibrated α is smaller than the widely accepted constant 1.26 and even smaller than the unit at several sites, which is physically unrealistic. Han and Tian (2018a) found that c corresponds to b in the AA by setting the B2015 approximately equal to the AA in the middle stage. However, the default value of c=0 corresponds to b with a value around 4.5, not the early default value of b=1. By calibrating both α and c, the B2015 performed well in estimating evaporation for 20 FLUXNET sites, and the value of α was more rational (Han and Tian, 2018a).

By contrast, two more parameters (x_min and x_max) are added to the sigmoid H2017. Because the sigmoid complementary curve is insensitive to x_min and x_max, Han and Tian (2018a) suggested that they could be treated as constant parameters for application convenience. x_min and x_max may change along with E_rad, and they were thought to vary with the timescales (Han and Tian, 2018a). x_min=0 and x_max=1 are appropriate at a daily scale for convenience, as has been evidenced by the good performances when compared to the flux measurements (Han et al., 2012; Han and Tian, 2018a). x_min and x_max are expected to be calculated by applying certain approaches to reduce the number of parameters of H2017 to two (Han and Tian, 2019).

Although α would vary in theory (Assouline et al., 2016), it is widely used with a constant value of 1.26 in practice (Priestley and Taylor, 1972). After calibrating, the variations of α are much less significant than those of the other parameters. Moreover, the calibrated α approaches 1.26, especially for the sigmoid H2017. Thus, the constant α=1.26 was suggested with acceptable weakening of the accuracy of E estimation (Han et al., 2012; Han and Tian, 2018a). In practice, α was also determined from the observed E values in wet condition when E was close to E_Pen and/or E_PT (Kahler and Brutsaert, 2006; Ma et al., 2015a; Wang et al., 2019). A novel method of using observed air temperature and humidity data in wet environments when measured E is lacking was proposed by Szilagyi et al. (2017) 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 much smaller than those in the oases or irrigated farmlands (Han et al., 2008, 2012). b was thought to be related to the characteristics of the atmosphere, i.e., the atmospheric humidity (Szilagyi, 2015); the Clausius–Clapeyron relation 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 land 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 parameterize b⁻¹. Wang et al. (2019) used the ecosystem mean b values of 217 sites around the world in the B2017 with little weakening of the evaporation estimation accuracy. However, the characteristics and determination methods of b need further studies toward a calibration-free evaporation estimation model.

4.3 Integrating with other approaches for further development

Actual evaporation is widely estimated as a reduction of the evaporation demand. The reduction factor was first taken as a function of soil moisture (Penman, 1950; Shuttleworth, 1993) or canopy resistance (Monteith, 1965). This Penman approach or Penman–Monteith approach has played a great role in parameterizing the evaporation process in hydrological models and land surface models. The canopy or surface temperature has also been widely used as a water stress indicator (Jackson et al., 1981, 1988), and the approach based on land surface temperature from remote sensing data has generated increasing attention. At annual or long-term timescales, the reduction factor is taken as a function of the humidity index represented by the ratio of precipitation to potential evaporation, and this method is known as the Budyko approach (Budyko, 1974; Yang et al., 2006; Zhang et al., 2001). In the above approaches, the evaporation demand is assumed to be independent of the land surface (Lhomme, 1997c; Morton, 1983). But in a large area where the land surface significantly interacts with the atmosphere, the evaporation demand will be altered by the changes of the land surface and the independent assumption does not hold. Although problems may not arise in diagnostic modeling as the current evaporation demand can be observed, they should be considered if these approaches are applied to a large area and used for future prediction or management in prognostic modeling (Han and Tian, 2018b).

As opposed to the above approaches that relied on the land surface properties, the reduction factor is determined from the atmospheric wetness in the generalized complementary functions (Table 3). The changes in evaporation demand due to the land surface properties are conceptually considered in the complementary principle, which is a theoretical improvement and would be helpful in predicting evaporation with land use changes. In addition, under conditions of the land surface properties being difficult to get, it is an obvious advantage of the complementary principle using the routinely observed meteorological variables in evaporation estimation. However, the complementary principle assumes that the changes in land surface properties can be accurately and timely detected from the changes of the atmospheric conditions. This assumption requires that the effects of regional or large-scale advection are negligible (Morton, 1983). Outside these situations, the generalized complementary functions may not work well because land surface properties are inadequately involved. Furthermore, the components of evaporation from different patches of the spatially heterogeneous surfaces, especially the evaporation from bare soil and the transpiration from vegetation, cannot be separated in the complementary principle, which is its disadvantage compared to the other approaches.

Considering the above disadvantages, Han and Tian (2018b) proposed a framework to integrate the complementary principle with other approaches for the advancement of evaporation research, which expresses E∕E_Pen as a function of both the land surface properties and the atmospheric wetness. Actually, both the land surface characteristics (e.g., soil moisture and vegetation) and atmospheric variables (e.g., radiation, humidity and temperature) have been used in the Jarvis–Stewart model (Jarvis, 1976; Stewart, 1988) to parameterize the canopy resistance. In fact, several attempts were conducted by integrating the complementary principle with other approaches to derive some of the land surface variables by using the meteorological variables (Han et al., 2015; Mallick et al., 2013; Szilagyi and Jozsa, 2009). A unified formulation of the Penman approach and the linear AA function was proposed by Crago and Brutsaert (1992). The integrated approach is a more rational conceptualization of the evaporation process from the unsaturated surface into the unsaturated atmosphere and is expected to increase the accuracy of evaporation estimation while reducing the burdens of parameterization. The findings of Liu et al. (2018) and Wang et al. (2019) that the parameters of the generalized functions significantly depends on the wetness of the land surface have demonstrated that the integrated approach has bright prospects. However, proper manners of integrating them need further study.

5 Conclusions

The complementary principle conceptualizes the feedbacks of land surface evaporation on atmospheric evaporation demand and offers advantages in evaporation estimation. In this study, the historical development of the complementary principle during the past half century was reviewed, and the two types of generalized complementary functions were focused on. In addition, future development for the generalized complementary principle was summarized based on the review. The concluding remarks are as follows:

1. The studies on the complementary principle first adopted a symmetric CR and then extended to an asymmetric CR. In the meantime, the original CR has evolved to the generalized complementary principle, which employs nonlinear functions as generalizations of the original linear relationship. The generalized complementary principle has a more rigorous physical base and offers potential in advancing actual evaporation estimation by using simple and standardized procedures.

2. Two types of generalized complementary functions were derived based on different understandings of the boundary conditions in completely wet environments: the sigmoid H12 and polynomial B15. The B15 inherits the concepts of “potential evaporation $E_{p_{\mathrm{0}}}$” and “apparent potential evaporation $E_{p_{\mathrm{a}}}$” from the original CR, and uses a polynomial function relating $E/E_{p_{\mathrm{a}}}$ to $E_{p_{\mathrm{0}}}/E_{p_{\mathrm{a}}}$. By contrast, H12 goes further from the original CR without involving the difficulties in formulating $E_{p_{\mathrm{0}}}$ and $E_{p_{\mathrm{a}}}$. Instead, a sigmoid function relating the ratio of actual evaporation to the Penman potential evaporation E_Pen and the proportion of the radiation component in E_Pen was derived. Nevertheless, further validation and application of the two types of generalized complementary functions are required with multiple data sets from different parts of the world.

3. Further studies on both the theoretical and practical aspects are still required before the generalized complementary principle achieves its potential. The generalized complementary principle requires a bold attempt for the practice of hydrology through enhancing its ability of evaporation estimation while reducing the burdens of parameterization. Thus, it should be carefully examined for its physical base of the boundary conditions in a completely wet environment and be integrated with other approaches to include the information of the land surface properly.

Appendix A: List of symbols

Table A1 List of symbols.