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.

Abstract. 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.

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 complementarytype 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 meteo-rological 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 . 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 pa − E po = E 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. (2011Han et al. ( , 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 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 func-tion are not limited by applicable timescales. In the derivation of the advanced generalized complementary functions (SGC of Tian, 2018, andPGC 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 −1 ) 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 landatmosphere feedback needs to achieve equilibrium.
In previous studies, the model validations were mostly completed at the daily scale 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(Han et al., , 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. Han et al. (2012Han et al. ( , 2018 proposed a generalized form of the complementary function that expresses E/E pen as a sigmoid function (SGC) of E rad /E pen :

Sigmoid generalized complementary function
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(Penman, , 1950, which can be expressed as where (kPa C −1 ) 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 −1 ) 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 * a and e a are the saturated and actual vapor pressures of air, respectively; z is the measurement height (Table S1); d 0 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 = 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: 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).

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 where c is an adjustable parameter. When c = 0, Eq. (5) reduces to

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.

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 nongrowing 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 −2 ) and energy balance closure ratio for each site are shown in Table S1.
The Nash-Sutcliffe efficiency (NSE; Legates and Mc-Cabe, 1999) is used to evaluate the efficiency of estimating E by the two generalized complementary functions: where E est (W m −2 ) is the estimated evaporation according to Eqs.
3 Results and discussion

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 2 ) 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)  When the timescale changes from day to month, the mean NSE H increases from 0.33 to 0.55, and R 2 H also increases from 0.61 to 0.75 (Table 1). However, they both decrease at the annual scale (NSE H = 0.18 and R 2 H = 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 −2 at the daily scale to 7.33 W m −2 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.
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 −2 ) is smaller than ours (16.36 ± 1.48 W m −2 ). 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 −2 ) is lower than our mean result of 88 sites (24.56 ± 0.95 W m −2 ).
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 twodimensional space of E rad /E pen − E/E pen move upwards. A detailed discussion about the variation in the parameters is provided in Sect. 3.4.

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.  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 2 B increases from 0.61 to 0.75. They decrease at the annual scale (NSE = 0.25 and R 2 H = 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. 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.
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 −2 ) and R 2 (0.65) are close to our results (RMSE = 26.83 ± 1.16 W m −2 and R 2 = 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 −2 , and the mean R 2 was 0.81, which are close to our mean results (RMSE = 19.17 ± 0.95 W m −2 and R 2 = 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.

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 1/α < E rad /E pen < 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 −2 (46 %), and RMSE B decreases by 13.13 W m −2 (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 1/α < E rad /E pen < 1.0 in the PGC function, and the resultant NSE B (0.19 vs. 0.19) and RMSE B (26.83 W m −2 vs. 26.68 W m −2 ) 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 2 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.

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. 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 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.
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 where f is the concave function, and x 1 and x 2 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 rad /E pen − E/E 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 exam-  ple, 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 rad /E pen , y = E/E 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: 1/b = 0.01c 2 +0.11c+0.24, with R 2 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).
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 −2 for α e vs. 16.67 W m −2 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).

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 2 , 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 2 values at the weekly scale (R 2 H = 0.64; R 2 B = 0.62) and those at the monthly scale (R 2 H = 0.62; R 2 B = 0.61) are similar. The RMSEs at the weekly scale are 14.95 W m −2 and 16.08 W m −2 for the sigmoid function and polynomial function, respectively, and those values at the monthly scale are 12.36 W m −2 (RMSE H ) and 12.93 W m −2 (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.

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 2 H (0.61), and RMSE H (10.16 W m −2 ) 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 2 B = 0.63; RMSE B = 9.94 W m −2 ) 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).

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 2 (> 0.96). The relationship at the monthly scale can be described as 1/b = 0.01c 2 + 0.11c + 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.
Code and data availability. All the data used in this study are from FLUXNET (http://fluxnet.fluxdata.org; U.S. Department of Energy, 2020; Baldocchi et al., 2001). The information for the EC data used in this study is listed in the Supplement. The intermediate data are available on request from the corresponding author (tianfq@mail.tsinghua.edu.cn).