At which timescale does the complementary principle perform best in evaporation estimation?

The complementary principle has been widely used to estimate evaporation under different conditions. However, it remains unclear at which timescale the complementary principle performs best. In this study, evaporation estimations were conducted at 88 eddy covariance (EC) monitoring sites at multiple timescales (daily, weekly, monthly, and yearly) by using sigmoid and polynomial generalized complementary functions. The results indicate that the generalized complementary functions exhibit the highest skill in estimating evaporation at the monthly scale. The uncertainty analysis shows that this conclusion is not affected by ecosystem type or energy balance closure method. Through comparisons at multiple timescales, we found that the slight difference between the two generalized complementary functions only exists when the independent variable (x) in the functions approaches 1. The results differ for the two models at daily and weekly scales. However, such differences vanish at monthly and annual timescales, with few high x values occurring. This study demonstrates the applicability of generalized complementary functions across multiple timescales and provides a reference for choosing a suitable time step for evaporation estimations in relevant studies.

1 Introduction

Terrestrial evaporation (E), including soil evaporation, wet canopy evaporation, and plant transpiration, is one of the most important components in the global water cycle and energy balance (Wang and Dickinson, 2012). The evaporation process affects the atmosphere through a series of feedbacks involving humidity, temperature, and momentum (Brubaker and Entekhabi, 1996; Neelin et al., 1987; Shukla and Mintz, 1982). Quantifying evaporation is crucial for a deep understanding of water and energy interactions between the land surface and the atmosphere. Generally, meteorological studies focus on evaporation changes at hourly and daily scales; hydrological applications require evaporation data at weekly, monthly, or longer timescales (Morton, 1983); and climate change studies focus more on interannual variations. The observation of E can occur at different timescales. For example, the eddy covariance, lysimeter, and scintillometer can measure evaporation at the half-hour scale, and water balance methods can observe evaporation at monthly to yearly scales (Wang and Dickinson, 2012). However, in most situations, an observation is unavailable, and the estimation of E is necessary. There are several types of methods for evaporation estimations, for example, the Budyko-type methods (Budyko, 1974; Fu, 1981), the Penman-type methods (Penman, 1948; Monteith, 1965), and the complementary-type methods (Bouchet, 1963; Brutsaert and Stricker, 1979). The Budyko-type methods perform well at annual or longer timescales, and the Penman-type methods can be applied at hourly and daily scales, while the complementary-type methods are used at multiple timescales (Crago and Crowley, 2005; Han and Tian, 2018; Ma et al., 2019) without explicit consideration of the timescale issue.

Recently, the complementary principle, as one of the major types of E estimation methods, has drawn increasing attention because it can be implemented with standard meteorological data (radiation, wind speed, air temperature, and humidity) without complicated underlying surface properties. Based on the coupling between the land surface and the atmosphere, the complementary principle assumes that the limitation of the wetness state in the underlying surface on evaporation can be synthetically reflected by atmospheric wetness (Han and Tian, 2020). Bouchet (1963) first proposed the “complementary relationship” (CR), which suggested that apparent potential evaporation (E_pa) and actual E depart from potential evaporation (E_po) in equal absolute values but opposite directions ($E_{\mathrm{pa}}-E_{\mathrm{po}}=E_{\mathrm{po}}-E$). According to the advection–aridity approach (AA; Brutsaert and Stricker, 1979), E_pa is formulated by Penman's (1948) equation (E_pen), and E_po is formulated by Priestley and Taylor's (1972) equation (E_PT). Subsequently, the CR was extended to a linear function with an asymmetric parameter (Brutsaert and Parlange, 1998). Further studies have found that the linear function underestimates E in arid environments and overestimates E in wet environments (Han et al., 2008; Hobbins et al., 2001; Qualls and Gultekin, 1997). To address this issue, Han et al. (2011, 2012) and Han and Tian (2018) proposed a sigmoid generalized complementary function (SGC; see Eq. 1 for details). As a modification to the AA approach, the SGC function illustrates the relationship between two dimensionless terms, E/E_pen and E_rad∕E_pen, where E_pen is Penman evaporation (Penman, 1948) andE_rad is the radiation term of E_pen. The SGC function shows higher accuracy in estimating E (Han and Tian, 2018; Ma et al., 2015b; Zhou et al., 2020) and outperforms the linear functions, especially in dry desert regions and wet farmlands (Han et al., 2012). Obtaining the impetus from Han et al. (2012), Brutsaert (2015) proposed a quartic polynomial generalized complementary function (PGC; see Eq. 5 for details). The PGC function describes the relationship between E∕E_pa and E_po∕E_pa, where E_pa and E_po are formulated in the manner of the AA approach. The PGC function has also been frequently used in recent years (Brutsaert et al., 2017; Hu et al., 2018; Liu et al., 2016; Zhang et al., 2017).

The prerequisite of the complementary principle is adequate feedback between the land surface and the atmosphere, which results in an equilibrium state. In this situation, the wetness condition of the land surface can be largely represented by the atmospheric conditions. Therefore, the timescales used in the complementary principle need to satisfy the adequate feedback assumption. However, this issue involves the complex processes of atmospheric horizontal and vertical motion, and these processes are difficult to explain theoretically. Morton (1983) noted this problem earlier and suggested that the complementary principle is not suitable for short timescales (e.g., less than 3 d), mainly because of the potential lag times associated with the response of energy and water vapor storage to disturbances in the atmospheric boundary layer. However, there is no solid evidence or theoretical identification to support this inference. The original complementary relationship and the AA function are not limited by applicable timescales. In the derivation of the advanced generalized complementary functions (SGC of Han and Tian, 2018, and PGC of Brutsaert, 2015), no specific timescale is defined. In practice, the complementary principle has been widely adopted to estimate E at multiple timescales, including hourly (Crago and Crowley, 2005; Parlange and Katul, 1992), daily (Han and Tian, 2018; Ma et al., 2015b), monthly (Ma et al., 2019; Brutsaert, 2020), and annual scales (Hobbins and Ramirez, 2004). The accuracy of the results has varied in different studies. Crago and Crowley (2005) found that the linear complementary function performs well in estimating E at small timescales of less than half an hour using the data from several well-known experimental projects (e.g., International Satellite Land Surface Climatology Project). The correlation coefficient between simulated E and observed E ranges from 0.87 to 0.92 in different experiments. The results of Ma et al. (2015b) indicated that the SGC function (root mean square error, RMSE = 0.39 mm d⁻¹) performs well in estimating E in an alpine steppe region of the Tibetan Plateau at the daily scale. Han and Tian (2018) applied the SGC function to the daily data of 20 eddy covariance (EC) sites from FLUXNET and found that it performed well in estimating E, with a mean Nash–Sutcliffe efficiency (NSE) value of 0.66. Crago and Qualls (2018) evaluated the PGC function and their rescaled complementary functions using the weekly data of seven FLUXNET sites in Australia, and the results showed that all the functions performed adequately, with a correlation coefficient between simulated E and observed E of higher than 0.9. Ma et al. (2019) also validated an emendatory polynomial complementary function at the monthly scale, and the NSE values of 13 EC sites in China were higher than 0.72. At the annual scale, Zhou et al. (2020) found that the mean NSE of the SGC function was 0.28 for 15 catchments in the Loess Plateau. Since these results were derived with different functions under varied conditions, it is difficult to determine at which timescale the performance is the best, and it is more difficult to explain theoretically how long the land–atmosphere feedback needs to achieve equilibrium.

In previous studies, the model validations were mostly completed at the daily scale (Brutsaert, 2017; Han and Tian, 2018; Wang et al., 2020), and the datasets of evaporation estimation were often established at the monthly scale (Ma et al., 2019; Brutsaert et al., 2020). However, each study only focused on a single timescale. In this study, we assessed the performance of the complementary functions in evaporation estimation at multiple timescales (daily, weekly, monthly, and yearly). The assessment was carried out at 88 EC monitoring sites with > 5-year-long observation records. In view of the fact that the complementary principle has developed to the nonlinear generalized forms, we selected two nonlinear complementary functions in the literature, i.e., the SGC function (Han et al., 2012, 2018) and the PGC function (Brutsaert, 2015). The key parameters of the complementary functions need to be determined by calibration. We chose the uniform database and the uniform parameter calibration methods for the optimization of the two complementary functions. We aimed to determine the most suitable timescale for the complementary functions through a comparison of the performances at different timescales. It is important for not only a deep understanding of the application of the complementary principle but also time step selection in evaporation database establishment and evaporation trend analysis.

This paper is organized as follows: Sect. 1 briefly describes the development of the complementary theory and our motivations to investigate the timescale issue. Section 2 describes the two functions, the parameter calibration method, and the data sources and processing. Section 3 shows and discusses the performance of the complementary functions at multiple timescales, the dependence of the key parameters on timescales, and the uncertainties in the analysis. The conclusions are given in Sect. 4.

2 Methodology

2.1 Sigmoid generalized complementary function

Han et al. (2012, 2018) proposed a generalized form of the complementary function that expresses E∕E_pen as a sigmoid function (SGC) of E_rad∕E_pen:

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

where x_max corresponds to the certain maximum value of x under extremely wet environments, and x_min corresponds to the certain minimum value of x under extremely arid environments. In this study, x_max and x_min were set as 1 and 0, respectively, for convenience. The E_pen term is defined by Penman's equation (Penman, 1948, 1950), which can be expressed as

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{E}_{\mathrm{pen}}& ={\displaystyle \frac{\mathrm{\Delta }(R_{\mathrm{n}}-G)}{\mathrm{\Delta }+\mathit{\gamma }}}+{\displaystyle \frac{\mathit{\rho }{c}_p}{\mathrm{\Delta }+\mathit{\gamma }}}{\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}u}{\ln \left(\frac{z-d_{\mathrm{0}}}{z_{\mathrm{0}\mathrm{m}}}\right)\ln \left(\frac{z-d_{\mathrm{0}}}{z_{\mathrm{0}\mathrm{v}}}\right)}}\\ & (e_{\mathrm{a}}^{*}-e_{\mathrm{a}}),\end{array}\end{array}$$

where Δ (kPa C⁻¹) is the slope of the saturation vapor curve at air temperature; R_n is the net radiation; G is the ground heat flux; γ (kPa C⁻¹) is a psychrometric constant; ρ is the air density; c_p is the specific heat; κ= 0.4 is the von Karman constant; u is the wind speed at measurement height; $e_{\mathrm{a}}^{*}$ and e_a are the saturated and actual vapor pressures of air, respectively; z is the measurement height (Table S1); d₀ is the displacement height; z_0m and z_0v are the roughness lengths for momentum and water vapor, respectively, which are estimated from the canopy height (h_c; Table S1), d₀=0.67h_c, z_0m=0.123h_c, and z_0v=0.1z_0m (Monin and Obukhov, 1954; Allen et al., 1998). E_rad is the radiation term of Penman evaporation:

$$\begin{array}{}\text{(3)}& E_{\mathrm{rad}}={\displaystyle \frac{\mathrm{\Delta }(R_{\mathrm{n}}-G)}{\mathrm{\Delta }+\mathit{\gamma }}}.\end{array}$$

The two parameters m and n of Eq. (1) can be determined by the Priestley–Taylor coefficient α and the asymmetric parameter b (Han and Tian, 2018).

$$\begin{array}{}\text{(4)}& \left\{\begin{array}{l}n=\mathrm{4}\mathit{\alpha }\left(\mathrm{1}+b^{-\mathrm{1}}\right)x_{\mathrm{0.5}}\left(\mathrm{1}-x_{\mathrm{0.5}}\right)\\ m={\left(\frac{x_{\mathrm{0.5}}}{\mathrm{1}-x_{\mathrm{0.5}}}\right)}^n,\end{array}\right.\end{array}$$

where x_0.5 is a variable that corresponds to y= 0.5 and equals $\frac{\mathrm{0.5}+b^{-\mathrm{1}}}{\mathit{\alpha }(\mathrm{1}+b^{-\mathrm{1}})}$.

2.2 Polynomial generalized complementary function

Brutsaert (2015) proposed the polynomial generalized complementary (PGC) function, which describes the relationship between E/E_pa and E_po∕E_pa. We uniformed the independent variable as E_rad∕E_pen to compare the two functions conveniently, and the polynomial function can be expressed as

$$\begin{array}{}\text{(5)}& y=\left(\mathrm{2}-c\right){\mathit{\alpha }}^{\mathrm{2}}{x}^{\mathrm{2}}-\left(\mathrm{1}-\mathrm{2}c\right){\mathit{\alpha }}^{\mathrm{3}}{x}^{\mathrm{3}}-c{\mathit{\alpha }}^{\mathrm{4}}{x}^{\mathrm{4}},\end{array}$$

where c is an adjustable parameter. When c= 0, Eq. (5) reduces to

$$\begin{array}{}\text{(6)}& y=\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}{x}^{\mathrm{2}}-{\mathit{\alpha }}^{\mathrm{3}}{x}^{\mathrm{3}}.\end{array}$$

2.3 Parameter optimization method

Typically, α has a default value of 1.26 (Priestley and Taylor, 1972). Since some studies have shown that a constant α may cause illogical results and biases in estimating E, it is suggested to specify α for diverse scenarios (Hobbins et al., 2001; Ma et al., 2015a; Sugita et al., 2001; Szilagyi, 2007). According to the complementary principle, under wet conditions, E is close to E_pen and the Priestley–Taylor's evaporation (E_PT=αE_rad). Specifically, when E∕E_pen is larger than a threshold (0.9 is commonly adopted), E_PT can be considered to be approximately equal to the observed E; thus, α can be calculated by E∕E_rad (Kahler and Brutsaert, 2006; Ma et al., 2015a). In this study, α was calculated using this method based on the mean value of E∕E_rad under wet conditions (E∕E_pen > 0.9). When all the E∕E_pen values are less than 0.9, α was set as the default value of 1.26. The key parameter b in SGC was calibrated by an optimization algorithm, with the objective function as the minimization of the mean absolute error (MAE) between the estimated E (by Eq. 1) and the observed E. Similarly, the key parameter c in PGC was calibrated by an optimization algorithm, with the objective function as the minimization of the MAE between the estimated E (by Eq. 5) and the observed E. Since we used the optimization algorithm to determine the parameter b in the SGC function, it is a fair choice to use the optimal c value instead of a constant value (c=0) in the PGC function.

To make the model parsimonious, we gave one value for the parameters (α, b and c) at each site for every different timescale. If the parameter was alterable, for example, it was month-dependent, and we would have to calibrate 12 parameters instead of one value for the whole study period. The purpose of this study is to determine the most suitable timescale for the complementary functions, and the variances of the key parameter within a timescale will introduce extra uncertainties. The accuracy will increase when an alterable parameter (that means a higher number of parameters) is used; however, the probability of overfitting risk will increase at the same time. In addition, in comparison to a group of parameters, a general representation of the parameter is more helpful in detecting its overall trend as the change in the timescale.

2.4 Data sources and data processing

The eddy flux data analyzed in this study were obtained from the FLUXNET database (http://fluxnet.fluxdata.org; Baldocchi et al., 2001). Observations from a total of 88 sites around the world were analyzed. Detailed information on these sites is listed in Table S1. These sites were selected from the FLUXNET database because they have observations from a period longer than 5 years. The 88 sites include 11 IGBP (International Geosphere-Biosphere Programme) land cover classes: ENF, evergreen needleleaf forests (27 sites); EBF, evergreen broadleaf forests (8); DBF, deciduous broadleaf forests (13); MF, mixed forests (5); OSH, open shrublands (4); CSH, closed shrublands (1); WSA, woody savannas (3); SAV, savannas (4); GRA, grasslands (15); CRO, croplands (6); and WET, permanent wetlands (2). The climates of the 88 sites range from arid to humid. Among the 88 sites, 11 sites have mean annual precipitation levels lower than 200 mm, 47 sites have precipitation levels between 200–500 mm, and 30 sites have precipitation levels above 500 mm. A total of 11 sites are located in the Southern Hemisphere (i.e., Australia, Brazil, and South Africa), and the others are located in the Northern Hemisphere.

Variables including net radiation, sensible heat flux, latent heat flux, ground heat flux, wind speed, air temperature, air pressure, precipitation, relative humidity, and vapor pressure deficit were acquired from the daily, weekly, and monthly datasets on the FLUXNET website. We analyzed the observations in the growing seasons from April to September for the Northern Hemisphere and from October to March for the Southern Hemisphere. These study periods were selected to avoid the high biases caused by the low level of solar radiation or extremely low evaporation (≈ 0) during the non-growing season. The seasonal and annual data were acquired by averaging the monthly data of the growing seasons. Following Ershadi et al. (2014), the energy-residual-corrected latent heat fluxes were used, which means the residual term in the energy balance is attributed to the latent heat to force the energy balance closure. To investigate the influence of different residual correction methods, the Bowen ratio energy balance method was also adopted in the uncertainty analysis. In the Bowen ratio method, the residual term is attributed to sensible heat and latent heat by preserving the Bowen ratio (Twine et al., 2000). The latent heat, sensible heat, and available energy (R_n – G) are restricted to positive values (Han and Tian, 2018). The energy balance residual (W m⁻²) and energy balance closure ratio for each site are shown in Table S1.

The Nash–Sutcliffe efficiency (NSE; Legates and McCabe, 1999) is used to evaluate the efficiency of estimating E by the two generalized complementary functions:

$$\begin{array}{}\text{(7)}& \mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{\sum (E-E_{\mathrm{est}}{)}^{\mathrm{2}}}{\sum (E-\stackrel{\mathrm{‾}}{E}{)}^{\mathrm{2}}}},\end{array}$$

where E_est (W m⁻²) is the estimated evaporation according to Eqs. (1) or (5), and $\overline{E}$ is the mean value of E (W m⁻²).

3 Results and discussion

3.1 Performance of the SGC function at multiple timescales

The relationship between the estimated E_est (site mean values) based on the SGC function (Eq. 1) and the observed E at the 88 sites at multiple timescales is shown in Fig. 1. The regression equations and determination coefficients (R²) were calculated using the site mean results. Each dot in Fig. 1 represents the site mean result averaged by daily (Fig. 1a), weekly (Fig. 1b), monthly (Fig. 1c), and yearly (Fig. 1d) results, and the total observation number is 88 (sites) at each timescale. Most of the results are near the 1:1 line, and all the regression slopes are close to 1 with high R² (0.95–0.99), which means the sigmoid function exhibits a good performance in estimating E at multiple timescales. The evaluation merits show that the performance varies at each timescale. The NSE values of the SGC functions for each site at different timescales are listed in Table S2. For the 88 sites, nearly half of the sites (40) have the highest NSE at the monthly scale, 12 sites have the highest NSE at the daily scale, 13 sites have the highest NSE at the weekly scale, and 23 sites have the highest NSE at the annual scale. The mean results of NSE_H, $R_{\mathrm{H}}^{\mathrm{2}}$, and RMSE_H (the subscript H corresponds to the sigmoid function proposed in Han and Tian, 2018) of these sites are shown in Table 1. $R_{\mathrm{H}}^{\mathrm{2}}$ represents the mean value averaged by the determination coefficients within each site. When the timescale changes from day to month, the mean NSE_H increases from 0.33 to 0.55, and $R_{\mathrm{H}}^{\mathrm{2}}$ also increases from 0.61 to 0.75 (Table 1). However, they both decrease at the annual scale (NSE_H= 0.18 and $R_{\mathrm{H}}^{\mathrm{2}}=$ 0.61). These results indicate that the SGC function exhibits the highest skill at the monthly scale. We inferred that there is a trade-off between the random error and the number of observations. RMSE_H values decrease from 24.56 W m⁻² at the daily scale to 7.33 W m⁻² at the annual scale, which means that the random error decreases as the timescale increases. At the same time, the fewer observations at the annual scale result in decreased variabilities of x and y, which affect the performance of the SGC function. On the other hand, Morton (1983) did not suggest using the complementary principle for short time intervals (e.g., less than 3 d), mainly considering the lag times associated with heat and water vapor change in the atmosphere, which may provide a possible inference for the weak performance at the daily scale.

Figure 1 The estimated evaporation based on the SGC function (Eq. 1) vs. the observed site mean evaporation at the daily scale (a), weekly scale (b), monthly scale (c), and yearly scale (d). Each dot represents the site mean result (N=88 in each panel). The regression equations and determination coefficients (R²) were calculated using the site mean results of the 88 EC sites.

Table 1 The evaluation merits (NSE, R², and RMSE in W m⁻²) of the two generalized complementary functions using the energy residual (ER) closure correction method. The subscripts H and B correspond to the SGC function proposed in Han and Tian (2018) and the PGC function proposed in Brutsaert (2015), respectively.

In previous studies, the SGC function was mainly applied at the daily scale. For example, the results of Ma et al. (2015b) in the alpine steppe region showed that the NSE of the sigmoid function is 0.73 at the daily scale, which is equal to our mean value in the grassland (0.73±0.08). The RMSE (11.06 W m⁻²) is smaller than ours (16.36±1.48 W m⁻²). The mean NSE of the 20 EC sites from FLUXNET is 0.66 at the daily scale in Han and Tian (2018), approximately 2 times the result in this study, and the RMSE (18.6 ± 0.94 W m⁻²) is lower than our mean result of 88 sites (24.56 ± 0.95 W m⁻²).

The SGC function for the five selected sites of different ecosystem types is shown in Fig. 2 to show the performance at multiple timescales (red lines in Fig. 2). These five EC monitoring sites were selected because they have long-term observations (> 10 years). The five sites include an evergreen needle forest (CA-TP1, Fig. 2a to d), a deciduous broad forest (US-UMB; Fig. 2e to h), a woody savanna (US-SRM; Figure 2i to l), a cropland (US-Ne2; Fig. 2m to p), and a grassland (US-Wkg; Fig. 2q to t). As observations decrease from the daily to the annual scale, the results converge on the middle part of the sigmoid curves and lie closer to the fitted lines. For some sites, the annual results concentrate on a narrow range with lower annual variabilities (e.g., Fig. 2h, l, and t). Generally, the key parameter (b) of the SGC function at these sites increases from the daily scale to the annual scale, which indicates that the sigmoid curves in the two-dimensional space of $E_{\mathrm{rad}}/E_{\mathrm{pen}}-E/E_{\mathrm{pen}}$ move upwards. A detailed discussion about the variation in the parameters is provided in Sect. 3.4.

Figure 2 Plots of E∕E_pen with respect to E_rad∕E_pen for five selected sites at multiple timescales. The black dots represent the observations; the red lines represent the SGC function; the green lines represent the PGC function; the blue lines are the P–T and Penman boundary lines. ENF, evergreen needleleaf forests; DBF, deciduous broadleaf forests; WSA, woody savannas; CRO, croplands; GRA, grasslands.

3.2 Performance of the PGC function at multiple timescales

The relationship between the estimated E_est (site mean values) based on the PGC function (Eq. 5) and the observed E at the 88 sites at multiple timescales is shown in Fig. 3. The slopes of the regression increase from 0.9 to 1 as the timescale changes from day to month and further increase to 1.01 at the annual scale. The intercept terms decrease from 13.06 W m⁻² at the daily scale to 0.01 W m⁻² at the monthly scale and further decrease to −0.25 W m⁻² at the annual scale. The R² values increase from 0.83 to 0.99 as the timescale increases. These coefficients of the regression show that the PGC function exhibits the highest skill at the monthly scale. The NSE values of the PGC functions for each site at different timescales are listed in Table S2. For the 88 sites, 42 sites have the highest NSE at the monthly scale, 7 sites have the highest NSE at the daily scale, 14 sites have the highest NSE at the weekly scale, and 25 sites have the highest NSE at the annual scale. The mean values of NSE_B, $R_{\mathrm{B}}^{\mathrm{2}}$, and RMSE_B (the subscript B corresponds to the polynomial function proposed in Brutsaert, 2015) of these sites are shown in Table 1. When the timescale changes from day to month, NSE_B increases from 0.19 to 0.50, and $R_{\mathrm{B}}^{\mathrm{2}}$ increases from 0.61 to 0.75. They decrease at the annual scale (NSE = 0.25 and $R_{\mathrm{H}}^{\mathrm{2}}=$ 0.63). Again, these evaluation merits indicate that the PGC function also exhibits the highest skill at the monthly scale, which is the same as for the SGC function.

Figure 3 As in Fig. 1 except for the PGC function (Eq. 5).

The PGC function has been applied at multiple timescales in previous studies. Zhang et al. (2017) evaluated the performance of the PGC function in estimating evaporation at four EC flux sites located across Australia, and their results showed that the mean RMSE (24.67 W m⁻²) and R² (0.65) are close to our results (RMSE = 26.83 ± 1.16 W m⁻² and R²= 0.61) at the daily scale. In Crago and Qualls (2018), the mean RMSE of seven EC sites at the weekly scale was 20.6 W m⁻², and the mean R² was 0.81, which are close to our mean results (RMSE = 19.17 ± 0.95 W m⁻² and R²= 0.7).

The PGC functions for the five selected sites are also shown in Fig. 2 (green lines). The fitted lines are almost the same as those of the SGC function in most situations when x is not too high. However, they diverge from each other when x becomes larger. Finally, y exceeds 1 when x is larger than 1∕α. Generally, the key parameter (c) of the PGC function at these sites decreases from the daily scale to the annual scale, which also indicates that the fitted curves move upwards.

3.3 Performance comparison of the SGC and PGC functions

The results from the 88 sites (Figs. 1, 3 and Table 1) show that the performances of the two functions are similar at monthly and annual timescales, while the SGC function performs slightly better than the PGC function at daily and weekly timescales. According to the results in Fig. 2, the two functions with calibrated parameters are approximately identical under non-humid environments, but their difference increases as x(E_rad∕E_pen) increases. We found that the values of α for all sites are greater than 1.0 in our study, which means that the PGC model cannot work properly under the condition of $\mathrm{1}/\mathit{\alpha }\mathit{<}{E}_{\mathrm{rad}}/E_{\mathrm{pen}}\mathit{<}\mathrm{1.0}$. At daily and weekly timescales, a substantial number of ecosystems can produce very high E_rad∕E_pen values. Specifically, 63 of the 88 sites have high E_rad∕E_pen values (x > 1∕α) at the daily scale, and 24 sites have high values at the weekly scale. However, there are only three sites with an x > 1∕α at the monthly scale, and no site has that value at the yearly scale. For the SGC function, in super-humid conditions, the upper part of the sigmoid curve is nearly flat and closer to the observations (e.g., Fig. 2a, m, and n). However, for the PGC function, theoretically, it cannot be applied when x is over 1∕α because the estimated E_est will be higher than E_pen, which is illogical. Thus, the sigmoid function performs slightly better at daily and weekly timescales than the polynomial function. However, the difference vanishes at the monthly scale as few high E_rad∕E_pen values occur.

According to the results, the performance of the PGC function is more sensitive to the time step than that of the SGC function. On the one hand, the regression relationship between E_est and the observed E of the 88 sites shows that the performance of the SGC function remains more stable (Fig. 1), while the regression results of the PGC function have higher variation when the timescale changes (Fig. 3). On the other hand, the estimation merits (Table 1) further confirm the sensitivity of the PGC function. From the daily scale to the monthly scale, the increase in NSE_H is 0.22, while the increase in NSE_B is 0.31; RMSE_H decreases by 11.36 W m⁻² (46 %), and RMSE_B decreases by 13.13 W m⁻² (49 %). At the daily scale, quite a few ecosystems (63 of 88 sites) can experience frequent high E_rad∕E_pen (> 1∕α) values, and the PGC function does not have the ability to simulate E accurately in this situation (E_est > E_pen), resulting in lower efficiency. We have carried out an additional analysis that adopts E=E_pen for $\mathrm{1}/\mathit{\alpha }\mathit{<}{E}_{\mathrm{rad}}/E_{\mathrm{pen}}\mathit{<}\mathrm{1.0}$ in the PGC function, and the resultant NSE_B (0.19 vs. 0.19) and RMSE_B (26.83 W m⁻² vs. 26.68 W m⁻²) present very similar results. As the timescale increases, the results converge on the middle part of the fitted line, and the number of high x greatly decreases (Fig. 2). Thus, the efficiency of the PGC function obviously increases. This is the reason that the polynomial function acts more sensitively to the time step.

In addition, we found that the two complementary functions perform reasonably well at shorter timescales (i.e., day and week), with relatively high R² values. Additionally, the estimations of site mean evaporation at shorter timescales are accurate (Figs. 1 and 3), especially for the SGC function. These results suggest that the generalized complementary functions have the ability to estimate evaporation accurately, even at shorter timescales.

3.4 Dependence of the key parameters of the SGC and PGC functions on timescales

The key parameters of the two complementary functions (b of the SGC function and c of the PGC function) vary at multiple timescales (Fig. 2). To explore their changes, the values of 1∕b and c at the 88 sites were averaged at each timescale. To take into account the situation in which b is equal to infinity, we used 1∕b instead of b in this analysis. Figure 4 shows the change in the two complementary functions with varied parameters at multiple timescales. The averaged 1∕b decreases from 0.45 ± 0.05 at the daily scale to 0.24 ± 0.03 at the annual scale (Fig. 4a), and the averaged c decreases from 0.98 ± 0.19 at the daily scale (Fig. 4b) to −0.37 ± 0.22 at the annual scale. The sign of c changes from positive to negative at the monthly scale.

Figure 4 Plots of the SGC Eq. (1) with α= 1.26 and varying 1∕b values at multiple timescales (a). Plots of the PGC Eq. (5) with α= 1.26 and varying c values at multiple timescales (b). The blue lines are the P–T and Penman boundary lines.

We show the histograms of 1∕b and c at multiple timescales in Figs. 5 and 6, respectively. At the daily scale, half of the 1∕b values are lower than 0.3, and the mean value is 0.45 ± 0.05. At the weekly scale, the peak of the distribution moves left, and almost half of the 1∕b values are lower than 0.2, with a mean value of 0.36 ± 0.04. At the monthly scale, the mean value is 0.29 ± 0.04, and the 1∕b values continue to decrease. At the annual scale, the mean value decreases to 0.24 ± 0.03, and 61 % of the 1∕b values are lower than 0.2. According to Fig. 6, at the daily scale, c follows a normal distribution (p value = 0.17, Kolmogorov–Smirnov test), with a mean value of 0.98 ± 0.21. Nearly 1∕3 of the c values are lower than 0. At the weekly scale, the center of the distribution moves left, with a mean value of 0.43 ± 0.24. Half of the c values are lower than 0. At the monthly scale, the mean value is −0.04 ± 0.23, and 58 % of the c values are lower than 0. At the annual scale, the mean value decreases to −0.37 ± 0.25, and 63 % of the c values are lower than 0. These results support our conclusion that 1∕b and c decrease as the timescale increases. Generally, the distribution of 1∕b and c also moves left within each ecosystem type according to Figs. 5 and 6.

Figure 5 Distribution of the key parameter 1∕b at the daily scale (a), weekly scale (b), monthly scale (c), and yearly scale (d). EBF, evergreen broadleaf forests (8); ENF, evergreen needleleaf forests (27); DBF, deciduous broadleaf forests (13); MF, mixed forests (5); shrub (12), closed shrubland, open shrublands, woody savannas, and savannas; GRA, grassland (15); CRO, croplands (6); WET, permanent wetlands (2).

Figure 6 Distribution of the key parameter c at the daily scale (a), weekly scale (b), monthly scale (c), and yearly scale (d). EBF, evergreen broadleaf forests (8); ENF, evergreen needleleaf forests (27); DBF, deciduous broadleaf forests (13); MF, mixed forests (5); shrub (12), closed shrubland, open shrublands, woody savannas, and savannas; GRA, grassland (15); CRO, croplands (6); WET, permanent wetlands (2).

The reduction in 1∕b and c indicates that the curves of the complementary functions move upwards as the timescale increases. Under non-humid conditions, the sigmoid function is a concave function, which means

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}[f\left(x_{\mathrm{1}}\right)+f\left(x_{\mathrm{2}}\right)]>f\left({\displaystyle \frac{x_{\mathrm{1}}+x_{\mathrm{2}}}{\mathrm{2}}}\right),\end{array}$$

where f is the concave function, and x₁ and x₂ represent any two values on the x axis. Since most of the results follow the fitted line, the averaged results of the longer time step will move upwards in the two-dimensional space of $E_{\mathrm{rad}}/E_{\mathrm{pen}}-E/E_{\mathrm{pen}}$, as well the new fitted curve. Although under super humid conditions, the SGC function is a convex function, there are fewer data under this condition as the timescale increases, and the shape of this part is almost unchanged (Fig. 4a). For the PGC function, when x is in the range of 0 to 1∕α, most of it is a concave function. For example, in the situation where c is equal to 0, the second derivative is higher than 0 as long as x is lower than 2∕3.

Furthermore, we found that the two key parameters b and c present a significant correlation, which provides additional evidence that the two functions can substitute each other in a sense. In other words, the two functions with calibrated parameters substantially provide similar descriptions of the distribution of the results in the state space ($x=E_{\mathrm{rad}}/E_{\mathrm{pen}}$, $y=E/E_{\mathrm{pen}}$). They can covert to each other in most situations since the two functions are generally equivalent to the linear asymmetric function when x is neither excessively large nor excessively small. The relationship can be described as follows: $\mathrm{1}/b=$ 0.01$c^{\mathrm{2}}+\mathrm{0.11}c+\mathrm{0.24}$, with R² being higher than 0.96 at the monthly scale (Fig. 7). The relationship remains at other timescales with a slight difference in the regression coefficients. At the daily scale, when c is equal to 0, the corresponding b is equal to 4.5, which is the same as that of the theoretical derivation in Brutsaert (2015).

Figure 7 Relationships between 1∕b and c at the monthly scale.

In this study, the physical meaning of the Priestley–Taylor coefficient α, which represents the ratio of E_PT (the Priestley–Taylor evaporation) and E_rad with the default value of 1.26 (Priestley and Taylor, 1972; Brutsaert and Stricker, 1979), was retained. This fundamental definition of α may result in a smaller range of E_rad∕E_pen in the PGC function. Liu et al. (2016) suggested that α_e (the calibrated α with c=0) in the PGC function is only a weak analog of the Priestley–Taylor coefficient, and Brutsaert (2019) directly considered α_e to be an adjustable parameter, which can be equal to or smaller than 1. We added the analysis that c is fixed to 0, and α is calibrated as α_e. This analysis showed that the two methods provide similar results (mean RMSE = 14.99 W m⁻² for α_e vs. 16.67 W m⁻² for α), and the conclusion of the timescale issue is consistent by adopting either α or α_e in the analysis. The optimal α_e has a significantly negative linear relationship with the optimal c, and the Pearson correlation coefficient is −0.8. This scenario suggests that calibrating either of the two parameters (α_e and c) is equivalent (Han et al., 2012).

3.5 Uncertainty analysis

3.5.1 Influence of ecosystem types

The evaluation merits of the generalized complementary functions may differ among ecosystem types. However, our results show that such variation generally does not affect our conclusion that the complementary functions perform best at the monthly scale. We show the performance of the two functions at multiple timescales for each ecosystem type in Table S3. Generally, the SGC function and the PGC function perform best at the monthly scale in most ecosystem types (9 of 11), with the highest NSE and R², which is consistent with the overall results. The exceptions include a closed shrubland site (CSH, N= 1) and evergreen broadleaf forests (EBF, N = 8), in which the complementary functions do not perform as well as in other ecosystem types. The CSH site (IT-Noe) has the highest NSE_H (0.11) and NSE_B (0.12) at the annual scale. In the EBF group, the highest NSE_H (0.15) and NSE_B (0.03) occur at the weekly scale, but the R² values at the weekly scale ($R_{\mathrm{H}}^{\mathrm{2}}=$ 0.64; $R_{\mathrm{B}}^{\mathrm{2}}=$ 0.62) and those at the monthly scale ($R_{\mathrm{H}}^{\mathrm{2}}=$ 0.62; $R_{\mathrm{B}}^{\mathrm{2}}=$ 0.61) are similar. The RMSEs at the weekly scale are 14.95 W m⁻² and 16.08 W m⁻² for the sigmoid function and polynomial function, respectively, and those values at the monthly scale are 12.36 W m⁻² (RMSE_H) and 12.93 W m⁻² (RMSE_B). We inferred that the abnormal results of these two exceptions are related to the lower NSE values in these ecosystem types. The mean NSE values at multiple timescales of CSH (−0.75) and EBF (−0.66) are negative, while the values of the other ecosystem types are all positive.

3.5.2 Performance at the seasonal scale

In consideration of the substantial discrepancy between the monthly results and the annual results, we added an analysis at the seasonal scale, which is between the two time steps. The relationship between the estimated E_est (site mean values) and the observed E of the 88 sites at the seasonal scale is shown in Fig. S1. For the SGC function, the regression result at the seasonal scale is similar to that at the monthly scale (Figs. S1a and 1c). The values of NSE_H (0.33), $R_{\mathrm{H}}^{\mathrm{2}}$ (0.61), and RMSE_H (10.16 W m⁻²) at the seasonal scale are between the monthly results and the yearly results (Table 1). For the PGC functions, the regression result at the seasonal scale is extremely close to that at the yearly scale (Fig. S1b and 3d). The evaluation merits (NSE_B = 0.31; $R_{\mathrm{B}}^{\mathrm{2}}=$ 0.63; RMSE_B= 9.94 W m⁻²) also range between the monthly results and the yearly results (Table 1). These results indicate that the decline in model efficiency has already occurred at the seasonal scale and support our conclusion that the complementary functions perform best at the monthly scale.

In addition, we also tested the influence of the different energy balance closure methods. The results based on both the energy residual (ER) closure correction (e.g., Ershadi et al., 2014; Han and Tian, 2018) and the Bowen ratio (BR) closure correction support our conclusion that the generalized complementary functions perform best at the monthly scale (Table S4).

4 Conclusions

In this study, evaporation estimations were assessed at 88 EC monitoring sites at multiple timescales (daily, weekly, monthly, and yearly) by using two generalized complementary functions (the SGC function and the PGC function). The performances of the complementary functions at multiple timescales were compared, and the variation in the key parameters at different timescales was explored. The main findings are summarized as follows:

1. The sigmoid and polynomial generalized complementary functions exhibit higher skill in estimating evaporation at the monthly scale than at the other evaluated scales. The highest evaluation merits were obtained at this timescale. The accuracy of the complementary functions highly depends on the calculation time step. The NSE increases from the daily scale (0.26, averaged by NSE_H and NSE_B) to the weekly scale (0.37) and monthly scale (0.53), while it decreases at the seasonal scale (0.32) and the annual scale (0.22). The regression parameters between estimated E_est and observed site mean E also support this conclusion for the PGC function. The variations among the different ecosystem types or between different energy balance closure methods generally have no effect on this conclusion. Further evaporation estimation studies with complementary functions can choose the monthly time step to achieve the most accurate results.

2. The SGC function and the PGC function are approximately identical under non-humid environments, while the SGC function performs better under super humid conditions implied by high values of x (> 1∕α) when the PGC function is theoretically useless (E_est > E_pen). At daily and weekly timescales, a substantial number of ecosystems can experience frequent high x values, and thus, the SGC function performs slightly better than the PGC function at these timescales. However, both functions perform very similarly at monthly and annual timescales, with few high x values. In addition, the performance of the PGC function is more sensitive to the time step than that of the SGC function.

3. The key parameter b of the SGC function increases and the key parameter c of the PGC function decreases as the timescale increases. The value of 1∕b is a quadratic function of c, with a higher R² (> 0.96). The relationship at the monthly scale can be described as $\mathrm{1}/b=\mathrm{0.01}{c}^{\mathrm{2}}+\mathrm{0.11}c+\mathrm{0.24}$. This relationship indicates that the two functions serve as substitutes to some extent.

   In this study to determine the most suitable timescale for applying the complementary principle, the key parameters (b and c) were calibrated to achieve the best model performance at each timescale. Further studies on the prognostic application of the complementary principle could focus on the reasonable prediction of the key parameters, and with the predictable flexible parameters at different timescales, the complementary principle could be integrated into hydrological models to reduce the uncertainty associated with evaporation estimations.