Evaluation of the Performance of Di ﬀ erent Methods for Estimating Evaporation over a Highland Open Freshwater Lake in Mountainous Area

: Lake evaporation is an important link connecting the water cycle and the surface energy cycle and remains one of the most uncertain terms in the local catchment’s water balance. Quantifying lake evaporation and its variability is crucial to improve water resource management and understand the response of the lake system towards climate change. In this study, we evaluated the performances of nine evaporation methods at di ﬀ erent timescales and calibrated them by using the continuous eddy covariance (EC) observation data during 2015–2018 over Erhai Lake, a highland open freshwater lake situated in the Dali valley, China. The nine evaporation methods could be classiﬁed into combination methods (Bowen-ratio energy budget, Penman, Priestley–Taylor, DeBruin–Keijman and Brutsaert–Stricker), solar radiation-based methods (Jensen–Haise and Makkink) and Dalton-based method (mass transfer and Ryan–Harleman) based on their parameterization schemes. The Dalton-based Ryan–Harleman method is most suitable for estimating evaporation at daily to weekly scales, while the combination methods and solar radiation-based method had good estimates at monthly timescale. After calibration, the biases of the Jensen–Haise and Ryan–Harleman method were slightly reduced, while the biases of the Makkink and mass transfer methods were reduced substantially. The calibrated Jensen–Haise method with small annual bias ( − 2.2~2.8%) and simple input variables was applied to estimate the long-term trend of evaporation during 1981–2018. The annual total evaporation showed an insigniﬁcant increasing trend of 0.30 mm year − 1, mainly caused by the signiﬁcant rising air temperature. This study showed the performance of evaporation methods over water bodies had large discrepancies on di ﬀ erent time scales, which indicated the importance of the choice of evaporation methods and provided instruction for water resource management of this region under climate change.


Introduction
Evaporation is a major water loss term in the water budget of inland water bodies, such as lakes and reservoirs [1]. As an important source of atmospheric moisture, nearly 60% of annual precipitation is returned to the atmosphere via evaporation on a global scale over land [2], and the proportion has increased to about 74% under climatic variability [3]. An increase in air and lake water surface temperature [4,5] and a decrease in the extent and duration of lake ice cover [6] have been observed under the background of global warming, therefore leading to an acceleration in lake evaporation [7].

Observation Site and Instruments
Erhai Lake is located in the city of Dali, Yunnan Province in China. It is surrounded by the Cangshan Mountains to the west and the Yu'an Mountains to the east (Figure 1). The Lake covers an area of 252.2 km 2 , with a south-north length of about 42.6 km and an east-west length of about 3-9 km. The average and maximum depths are 10 m and 20.7 m, respectively. This region lies in a subtropical monsoon climate zone, which is characterized by a warm-wet season from May to October and a cold-dry season from November to April [35]. According to the local climatology record [36], the monthly mean air temperature ranges from 8.4 • C in December to 20.3 • C in June, and the annual mean precipitation is 1055 mm, with 85% of precipitation concentrated in the wet season. There are mainly 18 rivers and streams from the foothills of Cangshan mountain flowing into the lake, and exists only one natural outlet-Xi'er River in the south. The water level is controlled artificially due to the agricultural irrigation for cropland during dry season, and ranges from 1971.1 m to 1974.1 m [37]. The water body is ice-free all year round and well mixed without obvious thermal stratification [38]. The platform (25.46 • N, 100.10 • E) is established approximately 70 m away from the west bank of Erhai Lake ( Figure 1). The height between the platform and the water surface ranged from 1.5 m to 3 m with the changing water level. The EC instruments were mounted at 2 m above the platform. The fluctuations of wind velocity components (u, v and w) and sonic temperature are captured by a three-dimensional sonic anemometer (CSAT3, Campbell Scientific, Logan, UT, USA). The water vapor and carbon dioxide density fluctuation are measured by an open-path infrared gas analyzer (IRGA, LI-7500, LI-COR, Lincoln, NE, USA). The data sampling frequency is 10 Hz. Sensors measured incoming and outgoing shortwave/longwave radiation (CNR1, Kipp and Zonen, Delft, The Netherlands) and active photosynthetic radiation (PAR, LI-190SB, Cambell Scientific, Logan, UT, USA) were mounted at 1.5 m above the platform. Air temperature and relative humidity are measured by HMP45C (Vaisala, Vantaa, Finland), and wind speed/direction was measured by 034B (Met One Instruments, Grants Pass, OR, USA) at 2.5 m above the platform. The water temperature profile at 0.05, 0.2, 0.5, 1, 2, 4, 6 and 8 m under the water surface are measured by water temperature sensors (CS616, Campbell Scientific, Logan, UT, USA). Measurement of basic meteorological variables and water temperature profile are taken every minute and recorded as 30-min averaged values. More information on the observation site and instruments is already documented in our previous study [39].

EC Data PostProcessing
The EC measurement dataset during the period from 2015 to 2018 was analyzed in this study. The raw 10 Hz, time-series data were processed into 30-min average values with the EddyPro software, Version 6.1 (LI-COR, Inc. 2015, Lincoln, NE, USA). The raw data point was considered as a spike when its magnitude exceeded the mean values 3.5 times the standard deviations within a certain moving time window and was further replaced through linear interpolation [40]. The double rotation method was adopted to align the x-axis wind component to the local streamline and turn the vertical and crosswind component into zero [41]. The block average based on the Reynolds decomposition rule was applied to remove a trend in the raw time series and obtain the fluctuations. Therefore, Latent and sensible heat fluxes were calculated from 30-min to mean covariances between the fluctuations of the scalars and the vertical wind speed component. Correction for density fluctuations was applied to latent heat and CO 2 fluxes using the method developed by [42]. The data quality control was applied to EC 30-min fluxes using the 0-1-2 flag system [43], which the fluxes with flag = 0 have the best quality, with flag = 1 is suitable for general analysis, and with flag = 2 do not pass the steady-state test and the developed turbulent conditions test and needed to be abandoned. Moreover, the data points were also removed when active gain control (AGC) values larger than 40 indicated the low-quality data caused by precipitation, dust, or other contamination on the sensor optics. Because the EC observation site is located near the west side of the lake, the flux footprint area is calculated by the footprint model [44]. As shown in Figure 1c, direction 225-315 • contains fluxes from the land surface is removed to eliminate the influence of the surrounding land area. After the data filtering as described above, the remaining 30-min data for sensible and latent heat fluxes are 69.8% and 74.3% in 2015, 67.1% and 70.8% in 2016, 69.7% and 72.0% in 2017, and 68.9% and 75.2% in 2018, respectively. To estimate the annual budgets of evaporation, an artificial neural network (ANN) method [38] is applied to fill the data gap in the daily timescale. The network was trained with four observed meteorological variables, including wind speed, air temperature, vapor pressure difference, and net radiation. The annual total evaporation losses are 1264.0 mm, 1214.6 mm, 1266.1 mm, and 1299.6 mm during 2015, 2016, 2017, and 2018, respectively.
To obtain the interannual variation of lake evaporation, a long term dataset (1981-2018) from measurements in the neighboring Dali National Climatic Observatory (DNCO, [36]) site and gridded China Meteorological Forcing Dataset (CMFD, [45,46]) are selected in this study. The DNCO site is about 15 km away from the Erhai Lake site (Figure 1a). The CMFD dataset was made through a fusion of remote sensing products, a reanalysis dataset and in situ observation data at the weather station, with a spatial resolution of 0.1 • and temporal resolution of three hours. Variables of wind speed, air temperature, relative humidity and precipitation are obtained from the DNCO measurements dataset. Incoming solar radiation is obtained from the long-term CMFD dataset. Using the land-based forcing data will likely introduce errors in calculating lake evaporation. Corrections of input variables were made to reduce systematic errors by comparing the land-based data with in situ measurements of Erhai Lake.

Introduction of Evaporation Methods
Nine evaporation methods based on different physical constraints and data requirements were selected for the estimation of open-water evaporation and classified into three categories, including the combination methods, solar radiation-based methods and Dalton-based methods.

Combination Methods
The Bowen ratio energy budget (BREB) method was first introduced by [18] and was generally regarded as a standard method for estimating evaporation over water surfaces in the absence of in situ observation [1,25]. The formula of the BREB method is given as follows: where E is the estimated evaporation rate (mm d −1 ), R n is the net radiation (W m −2 ), ∆Q is the heat storage change in the water body (W m −2 ), ρ w is the water density (1000 kg m −3 ), λ is the latent heat of vaporization (2.45 × 10 6 J kg −1 ), β is the Bowen ratio (dimensionless), T s is the water surface temperature ( • C), T a is the air temperature ( • C), P is the atmospheric pressure (kPa), ε is the ratio of molecular weight of the water to that of the dry air (0.622, dimensionless), e * s is the saturated vapor pressure at the water surface temperature (kPa), e * a is the saturated vapor pressure at the air temperature (kPa), e a is the vapor pressure at the air temperature (kPa), RH is the relative humidity (%), c pw is the specific heat of water at constant pressure (4192 J kg −1 • C), c p is the specific heat of air at constant pressure (1004 J kg −1 • C), D is the multipliers that convert the unit of evaporation into mm d −1 (86.4 × 10 6 ), ∆T w ∆t is the mean water temperature change during the time interval ∆t ( • C s −1 ), z d is the largest depth at the profile measurement (m).
The Penman (PM) method [20] combined the surface energy budget and aerodynamic components, which was developed to calculate the potential evaporation from an open water surface. The PM equation is expressed as: where ∆ is the slope of saturated vapor pressure-air temperature curve (kPa • C −1 ), γ is the psychometric constant (kPa • C −1 ), f (u) is the wind function (W m −1 hPa −1 ), U is the wind speed at the measurement height (m s −1 ). The first term of the PM method is usually defined as equilibrium evaporation; the lower limit of evaporation when the air above the water surface is saturated. The second term of the PM method is generally called the aerodynamic component correlated to the drying power of the air. The Priestley-Taylor (PT) method [47] was originally developed for estimating potential evaporation over open water surfaces and saturated land surfaces under advection-minimum conditions. The PT method is a simple expression of the PM equation, which introduced a constant to account for the aerodynamic component of the PM method. The PT formula is expressed as: where α is the Priestley-Taylor empirical constant and taken as a constant of 1.26. The DeBruin-Keijman (dBK) method [48] is derived from the PT equation and determined as: The following empirical relationship between β and γ/∆ [49] is applied in the derivation of the dBK method: The Brutsaert-Stricker (BS) method [50] also called the advection-aridity model, is based on the concept of a symmetric complementary relationship that depends on the feedback between actual evaporation (E a ) and potential evaporation (E p ), where E w is the wet environment evaporation and is often calculated by the PT method, E p is generally calculated by the PM method. E a is the evaporation defined by the BS method, and therefore the formula can be written as:

Solar Radiation-Based Methods
Two different generalized forms of solar radiation-based methods, including the Jensen-Haise (JH) method [51] and Makkink (Mak) method [52], were evaluated in this study. The JH method was originally developed for the arid environment in the western United States and determined by the following expression: where R s is the incoming solar radiation (W m −2 ). The Mak method was first developed for estimating evapotranspiration for grasslands in The Netherlands and is currently applied to open water. The equation is described as follows:

Dalton-Based Methods
The Dalton-based methods, consisting of mass transfer (MT) method [21,53] and Ryan-Harleman (RyH) method [54], estimated the evaporation through wind speed (U) and vapor pressure difference (VPD) between the lake surface and the atmosphere. A simple formula of the MT method is given as follows: where N is the mass transfer coefficient and parameterized as a function of lake surface area A s (acre) [21]. The RyH method considered the two main evaporation-driven processes: free convection caused by buoyancy and forced convection caused by wind. The formula of the RyH method is given as follows:

Calibration of Evaporation Methods
Due to the difference in physical and climatic conditions, some parameters or empirical coefficients can be replaced or calibrated to achieve a relatively better agreement with observations. The details of the adjustments are summarized in Table 1. Since long term continuous observation of temperature profiles sometimes is absent because of the instrument failure, the determination of heat storage in the water body needs alternate expression related to meteorological factors. For the combination methods, the heat storage is determined using the measured water temperature profiles in the default condition (V0) and a hysteresis function of net radiation [55] in the calibrated condition (V1), respectively. The initial climate conditions assumed for the two solar radiation-based methods are totally different from the Erhai Lake, and therefore its empirical coefficients are calibrated site-specific through multiple regression against measured evaporation. The transfer coefficient N for the MT method, which reflects the efficiency of vertical transport of water vapor by turbulent eddies induced from wind shear, is calibrated by using the aerodynamic formula [22]. The empirical constants of the RyH method are also modified by forcing multiple regression against measured data under the condition V1. Table 1. Comparisons between the default condition V0 and calibrated condition V1 for the evaporation equations.

V0 V1
Combination methods BREB ∆Q is calculated from measured water temperature profiles.

Evaluation Criteria of Evaporation Methods
To evaluate the performance of evaporation methods from daily to monthly timescale, statistical measures including linear correlation coefficient (R), root mean square error (RMSE) and mean absolute error (MAE) between the observed and estimated evaporation were selected and calculated as follows: where n denotes the number of data points, OBS i and MOD i are the observed and modeled evaporation, respectively, OBS and MOD are mean values of the observed and modeled evaporation, respectively.

Assessment of Relative Contribution
In order to quantify the relative contribution of meteorological factors to the long-term trend of annual evaporation, the detrend method [56] and the sensitivity method [57] were used. In the detrend method, the relative contribution of a factor xi is determined by the following formula: where RC detrend (xi) is the relative contribution of factor xi to evaporation changes based on the detrend method, dE original dt is the trend of lake evaporation calculated by using the original inputs, is the trend of lake evaporation calculated by using the linear detrend xi and original data for other variables, nt is the total number of affecting factors. In the sensitivity method, the relative contribution of a factor xi is determined by the following formula: where RC sensitivity (xi) is the relative contribution of factor xi to evaporation changes based on the sensitivity method, ∂E original ∂xi is the partial derivatives of xi, dxi dt is the trend of xi during the long-term period.

Meteorological Conditions and Surface Energy Budget
The general characteristic of meteorological conditions during 2015-2018 is shown in Figure 2.  17.69 • C, which were slightly larger than the annual mean T a . The daily average temperature difference (∆T) between the lake surface and the atmosphere was almost negative from February to June and was positive during other months in a year. Annual mean ∆T is always positive during 4-year period. The daily average vapor pressure difference (VPD) gradually increased with seasons but slightly dropped in summer due to a large amount of rainfall. The daily average incoming solar radiation (R s ) generally peaked in mid-summer and reached its minimum in winter. Small daily mean values of R s were also observed in summer because of the influence of rainfall events. The precipitation was mainly concentrated on the monsoon season from May to October, accounted for 76% to 89% of the annual precipitation during 2015-2018. The annual total precipitation in 2016 (1153.6 mm) was larger than the long-term mean annual precipitation (1055 mm), while the annual precipitation in 2015 (1016.2 mm), 2017 (909.7 mm), and 2018 (1000.8 mm) were slightly lower than the long-term average. The maximum LE was about one month to three months lag behind the maximum R n over Erhai Lake. For large and deep lakes with the high heat capacity of lake water, a two to five months delay between maximum LE and R n was observed [17,[58][59][60]. The seasonal variation of the heat storage change (∆Q) can be divided into storaging period (∆Q > 0) and releasing period (∆Q < 0). The heat was absorbed into the water body during spring to summer and used to warm the lake water. The peak value of heat absorption was observed in May. A slight fluctuation of heat absorption occurred during summer due to a large amount of rainfall. The heat released into the atmosphere since September. The maximum heat release usually occurred from October to November. The seasonal variation of ∆Q has a significant influence on the estimation of seasonal lake evaporation by weakening evaporation during spring and summer and enhancing it largely during the autumn and winter [61]. The sensible heat flux (H) was relatively small when compared with LE, with monthly average values ranged from −11.9 W m −2 to 19.5 W m −2 . Monthly mean H was usually remained negative from February to June when T a is greater than T s and became positive afterward when T a is smaller than T s . Based on the surface energy partitioning, most of the R n absorbed by the water body was consumed for lake evaporation with annual mean LE/R n ranged from 0.85 to 0.88 during 2015-2018. Monthly LE/R n exceeded 1 during September (only for 2018)/October to January, indicated the energy stored in the water body served as an additional energy source of evaporation during these periods. A similar result was also found in Ross Barnett Reservoir under humid and subtropical climate, with an annual average LE/R n of 0.81 and monthly average LE/R n > 1 during September to January [62]. For a hypersaline terminal lake under an arid environment, LE accounted for a lower portion of R n with an average value of 65% over four seasons [34]. The energy balance closure (EBC) is quantified by the ratio of the turbulent flux (LE + H) to the available energy (R n − ∆Q). Annual mean energy EBC ranged from 0.86 to 0.97 during 2015-2018 and was comparable to values obtained from other lakes [15,31,62] and terrestrial sites [63].

Evaluation of Combination Methods
The performances of five combination methods were improved as the timescale is extended from daily to monthly. Under the default condition V0, the mean correlation coefficient (R) of the five methods increased from 0.53 for daily to 0.85 for monthly timescale, the RMSE decreased from 1.85 to 0.83, and the MAE reduced from 1.30 to 0.67, respectively (Table 2, V0). Under the calibrated Condition V1, the heat storage change ∆Q is considered as a hysteresis function of R n . The estimate of evaporation using hysteresis ∆Q was also comparable to the default condition using measured ∆Q as input. The coefficient of determination R 2 between measured and modeled ∆Q is 0.84 at monthly timescales during 2015-2018, which is identical to the average results of 22 lakes [55]. Using hysteresis relation is able to provide reasonable ∆Q estimates when measurements of water temperature profiles are not available, and produce good estimates of monthly evaporation for the combination methods. Reasonable evaporation estimates can also be achieved by taking heat storage as a simple linear function of R n in Poyang lake [64] with a mean depth of 8 m comparable to the mean depth of Erhai Lake (10 m). The difference in monthly average evaporation between five combination methods and EC observations showed small seasonal bias under the default and calibrated conditions. The combination methods overestimated evaporation during spring and summer and underestimated it during late autumn to winter (Figure 4). Variation of seasonal bias was almost consistent among the five combination methods. Under the default condition V0, the departure was within the range of −1.5 mm d −1 to 1.5 mm d −1 , except for the BS method with bias beyond this range. Calibrated condition V1 showed similar seasonal bias related to default condition V0. More positive bias during April (2017 and 2018) and May (2015 and 2016) and less negative bias during December were found in V1 (Figure 4, V1). Such seasonal bias of combination methods has also been found in other lakes [25,30,34,65]. The combination methods are based on the assumption of closed surface energy balance. However, the turbulent fluxes measured by the EC system are generally underestimated, and therefore, the surface energy balance is not always closed [66], eventually leading to seasonal deviation of combination methods. the Andes mountain range [67]. More favorable results from the PT than the dBK method coincided with the analysis in a small highland lake, which pointed out the PT method is more suitable for estimating evaporation over high altitude lakes than the dBK method developed from the sea-level environment [32]. The P-T coefficient taken as a constant of 1.26 without consideration its seasonal variability (larger than 1.26 in winter and smaller than 1.26 in summer, not shown) might introduce uncertainties and cause larger RMSE and MAE of the PT method in comparison with the PM method. However, in Poyang lake, a large shallow subtropical lake in China, no obvious seasonal trend in α was observed [68], and the PT method performed slightly better than the PM method [64]. The BS method displayed the largest bias for lake evaporation estimation among the combination methods from daily to monthly timescales, where its MAE value was almost twice as much as the value from the PM method.

Evaluation of Solar Radiation-Based Methods
The performances of the solar radiation-based methods were similar to the combination methods. The radiation-based methods yielded poor estimates at short timescales (daily to weekly) and relatively good estimates at monthly timescale. Under the condition V0, the default version of two radiation-based methods compared less favorably with most of the combination methods, except for the BS method. The difference in monthly average evaporation between the default JH method   Among all the combination methods, the PM method has the best performance with large R and small RMSE and MAE under the range of daily to monthly timescales, followed by the PT, dBK, BREB, and BS method. Good performance of the PM method is also found in other lakes, e.g., a large deep lake in Egypt [29], the Dead sea in an arid environment [34], and the Laja Lake in the middle of the Andes mountain range [67]. More favorable results from the PT than the dBK method coincided with the analysis in a small highland lake, which pointed out the PT method is more suitable for estimating evaporation over high altitude lakes than the dBK method developed from the sea-level environment [32]. The P-T coefficient taken as a constant of 1.26 without consideration its seasonal variability (larger than 1.26 in winter and smaller than 1.26 in summer, not shown) might introduce uncertainties and cause larger RMSE and MAE of the PT method in comparison with the PM method. However, in Poyang lake, a large shallow subtropical lake in China, no obvious seasonal trend in α was observed [68], and the PT method performed slightly better than the PM method [64]. The BS method displayed the largest bias for lake evaporation estimation among the combination methods from daily to monthly timescales, where its MAE value was almost twice as much as the value from the PM method.

Evaluation of Solar Radiation-Based Methods
The performances of the solar radiation-based methods were similar to the combination methods. The radiation-based methods yielded poor estimates at short timescales (daily to weekly) and relatively good estimates at monthly timescale. Under the condition V0, the default version of two radiation-based methods compared less favorably with most of the combination methods, except for the BS method. The difference in monthly average evaporation between the default JH method and EC observations also had a seasonal cycle, in which lake evaporation was generally overestimated during spring and summer and was underestimated during late autumn and winter, resulting in an annual bias of −2.6% to 4.3% during 2015-2018 ( Figure 5, V0). The default Mak method underestimated the monthly average evaporation almost all year, leading to a large negative annual bias of −13.4% to −17.8% during 2015-2018 ( Figure 5, V0). Under the calibrated condition V1, the simulation bias of the default radiation-based methods was reduced after calibration with smaller RMSE and MAE, while the correlation between observation and modified JH method declined slightly compared to the default JH method ( Table 2, V1). The large positive bias of the default JH method in summer was reduced or became negative after calibration ( Figure 5, V1). The difference in annual bias caused by the JH method before and after calibration was not significant. The modified Mak method reduced the large negative deviation of default condition and overestimated evaporation during spring, consequently resulting in less negative annual bias ranged from −4.4% to 0.65% during 2015-2018 ( Figure 5, V1). For the sake of simplicity and accuracy, the modified solar radiation-based methods have an advantage over the combination methods and were applied in estimating the long-term lake evaporation [32,69].

Evaluation of Dalton-based Methods
In contrast to the two former methods based on the surface energy budget and solar radiation, the Dalton-based methods had a larger correlation with observations at short timescales (daily to weekly) and relatively poor correlation at a monthly timescale. Under the default condition V0, the RyH method had much better performance than the MT method in terms of smaller RMSE and MAE values among the entire range of timescales (Table 2, V0). With the transfer coefficient N parameterized as a function of the lake surface area [21], the default MT method strongly overestimated the monthly average evaporation during all the year except for August 2018, thus leading to overestimation of the annual total evaporation by 25.7% to 40.6% during period 2015-2018 ( Figure 5, V0). A similar result was also found in Lake Taihu [31] that the MT method caused an overestimation of annual evaporation up to 87%. For the default RyH method, the monthly departure of evaporation was within the scope of −1 to 1 mm d −1 . The underestimation generally occurred during summer to early autumn. The annual bias was relatively small and varied from −2.9% to 5.2% during 4-yr period ( Figure 5, V0). Under the calibrated condition V1, two modified Dalton-based methods both reduced the bias in terms of smaller RMSE and MAE values from daily to monthly scales ( Table 2, V1) compared to the default condition. Continuous strong overestimation of monthly average evaporation of default MT method was reduced substantially after using the aerodynamic formula of mass transfer coefficient , and generally shifted into negative bias during the second half the year, resulting in a small positive annual bias in 2015 and 2017, and negative annual bias in 2016 and 2018 ( Figure 5, V1). The difference in an annual bias of the RyH method was not significant before and after calibration, indicated the RyH method without site-specific calibration performed well over Erhai Lake. The good behavior of the default RyH method is also found in the study of a variety of lakes and ponds in Minnesota [70] and small mountain lakes in southern Wyoming, USA [71] because of the large variability in wind speed and relative humidity.

Inter-Comparisons between Evaporation Methods at Different Timescales
Nine evaporation methods are related to different main drivers of lake evaporation. The correlation between meteorological variables and lake evaporation obtained from observations and default evaporation methods across daily to yearly timescales during 2015-2018 is shown in Figure  6. On the basis of the measurement results, lake evaporation had the largest correlation with from daily to weekly timescale, and the relationship gradually decreased as the time interval is extended, which correlation dropped from 0.76 on a daily scale to 0.01 on a yearly timescale. Lake evaporation had a greater correlation with and on a monthly timescale. On an annual timescale, evaporation was highly correlated to VPD and radiation ( and ) with a correlation higher than

Evaluation of Dalton-based Methods
In contrast to the two former methods based on the surface energy budget and solar radiation, the Dalton-based methods had a larger correlation with observations at short timescales (daily to weekly) and relatively poor correlation at a monthly timescale. Under the default condition V0, the RyH method had much better performance than the MT method in terms of smaller RMSE and MAE values among the entire range of timescales ( Table 2, V0). With the transfer coefficient N parameterized as a function of the lake surface area [21], the default MT method strongly overestimated the monthly average evaporation during all the year except for August 2018, thus leading to overestimation of the annual total evaporation by 25.7% to 40.6% during period 2015-2018 ( Figure 5, V0). A similar result was also found in Lake Taihu [31] that the MT method caused an overestimation of annual evaporation up to 87%. For the default RyH method, the monthly departure of evaporation was within the scope of −1 to 1 mm d −1 . The underestimation generally occurred during summer to early autumn. The annual bias was relatively small and varied from −2.9% to 5.2% during 4-year period ( Figure 5, V0). Under the calibrated condition V1, two modified Dalton-based methods both reduced the bias in terms of smaller RMSE and MAE values from daily to monthly scales ( Table 2, V1) compared to the default condition. Continuous strong overestimation of monthly average evaporation of default MT method was reduced substantially after using the aerodynamic formula of mass transfer coefficient N, and generally shifted into negative bias during the second half the year, resulting in a small positive annual bias in 2015 and 2017, and negative annual bias in 2016 and 2018 ( Figure 5, V1). The difference in an annual bias of the RyH method was not significant before and after calibration, indicated the RyH method without site-specific calibration performed well over Erhai Lake. The good behavior of the default RyH method is also found in the study of a variety of lakes and ponds in Minnesota [70] and small mountain lakes in southern Wyoming, USA [71] because of the large variability in wind speed and relative humidity.

Inter-Comparisons between Evaporation Methods at Different Timescales
Nine evaporation methods are related to different main drivers of lake evaporation. The correlation between meteorological variables and lake evaporation obtained from observations and default evaporation methods across daily to yearly timescales during 2015-2018 is shown in Figure 6. On the basis of the measurement results, lake evaporation had the largest correlation with U from daily to weekly timescale, and the relationship gradually decreased as the time interval is extended, which correlation dropped from 0.76 on a daily scale to 0.01 on a yearly timescale. Lake evaporation had a greater correlation with T a and R n on a monthly timescale. On an annual timescale, evaporation was highly correlated to VPD and radiation (R n and R s ) with a correlation higher than 0.8, while it had a negligible relationship with U and T a . Apparently, VPD, R n , and R s both had an increasing correlation with lake evaporation as the timescale is extended from daily to annual. The combination and solar radiation-based methods failed at short timescales (daily to weekly), due to the mismatched controlling factors as lake evaporation was highly correlated with U while these methods were controlled by T a and radiation (R n and R s ), and showed relative good performances at monthly timescale due to the consistency of large relationship with T a and radiation (R n and R s ). The Dalton-based methods showed large correlations with U and VPD, make it possible to yield good estimates of evaporation from daily and weekly timescale, while produced less satisfactory results at monthly timescale because lake evaporation is largely correlated to T a and radiation. Due to the difference in weather and climatic conditions, lake evaporation dynamics are driven by different atmospheric forcing and processes, which will influence the suitability of evaporation methods among lakes. For instance, the good performance of calibrated Dalton-based methods during a shorter timescale over Erhai Lake is different from the results from lake Okeechobee in subtropical South Florida [72], because lake Okeechobee experienced a large amount of rainfall and small VPD condition, lake evaporation is mainly governed by solar radiation. Therefore, the Dalton-based methods driven by U and VPD is not suitable, and the solar radiation-based methods can provide more realistic estimates of evaporation. The solar radiation-based methods also have an advantage over the Dalton-based methods in estimating evaporation from lakes and reservoirs in a semi-arid environment [30].
The combination methods agreed well with the solar radiation-based methods across the entire range of timescale as they have similar controlling factors-T a and radiation (R n and R s ) (Figure 7). These two methods showed obvious disagreement with the Dalton-based methods at shorter timescales because the latter is mainly driven by U and VPD, indicated the importance of choosing the appropriate method for lake evaporation estimation at a smaller timescale. The correlations between the methods increased with the extended timescale, indicated that the deviation of methods based on different physical parameterization would be reduced on a longer timescale. Despite the more consistency at the annual timescale, large differences could still occur between these methods.

Estimation of Long-Term Lake Evaporation
To investigate the interannual variation of lake evaporation during the past decades (1981-2018), the calibrated JH method was forced with historical measurements dataset from the DNCO station and the CMFD dataset. The calibrated JH method had a relatively small annual bias varied from −2.6% to 2.8% during 2015-2018 (Figure 4, V1), and only required T a and R s as input, which both easily obtained, therefore made it become a practical and reliable option for estimating the long-term trend of evaporation.
The estimated annual evaporation displayed an insignificant increasing trend of 0.30 mm year −1 during 1981-2018 and ranged from 1146.2 mm to 1362.0 mm with a mean value of 1249.0 mm (Figure 8a). The low annual evaporation in the early 1990s resulted from the small incoming solar radiation. The increasing trend of lake evaporation has also been reported in lakes with different climate background, including a small highland lake in the Tibetan Plateau [32], Lake Taihu in a subtropical region of China [73], Lake IJssel in the Netherlands [69] and a small reservoir in the Brazilian savannah [74]. The long-term trend of evaporation of Erhai Lake calculated by the modified JH method is opposite to the results obtained from evaporation pan E601 situated in the surrounding land site [75], in which pan evaporation showed an obvious decreasing trend since the 1980s with rising temperature. This inverse trend is known as the "evaporation paradox", which is in agreement with the results from reservoirs in the contiguous United States [76]. This discrepancy is likely attributed to different environmental and meteorological conditions and the limited measured area of the evaporation pan. During period 1981-2018, T a increased significantly with a trend of 0.046 • C year −1 , while U declined significantly at a rate of −0.016 m s −1 year −1 , R s and annual total precipitation both decreased slightly at rates of −0.037 W m −2 year −1 and −1.4 mm year −1 , and specific humidity q a changed insignificantly (Figure 8). The increasing lake evaporation accompanied by decreasing precipitation is likely to reduce the water resources of Erhai Lake, which is supported by the decreasing trend of water resources observed in this basin [77] and in the mainland of China [78]. The climate over Erhai Lake basin was mainly cold and wet in the 1960s and 1970s, but getting warm after into 1980s and warmer and dry into the 21st century [77]. This changing climatic condition may cause an increase in lake evaporation since the 1980s despite the incoming solar radiation decreased during this period. The relative contributions of T a and R s were quantified by the detrend method and sensitivity method to find out what caused the evaporation trend during the past decades, as shown in Figure 9. A large difference occurred between the original evaporation and the recalculated evaporation under the detrended T a condition. The increasing trend (0.30 mm year −1 ) of evaporation based on the original inputs was largely reduced to −0.24 mm year −1 when using the detrended T a and original R s as input (Figure 9c). Under the detrended R s condition, the trend of original annual evaporation was slightly increased to 0.54 mm year −1 (Figure 9c). Based on the detrend method, the increasing T a contributed to 69.2% of evaporation variation, while the decreasing R s resulted in −30.8% of evaporation change (Figure 9f). Relative contributions of these two factors based on the sensitivity method were consistent with the results of the detrend method. The change in T a contributed 69.5% to evaporation variation (Figure 9f). Therefore, the rising T a contributed to the increasing trend of evaporation during 1981-2018, indicating the impact of changing climate on the hydrological cycle over Erhai Lake. Because the variation of R s was relatively small and insignificant during 1981-2018, its contribution to the long-term trend of evaporation was not the largest, although the magnitude of annual lake evaporation changed largely to the disturbance of R s . It has been reported that annual evaporation for lakes in low-altitude areas is projected to increase about 200 mm by 2091-2100 in comparison with 2006-2015, in spite of little changes in incoming solar radiation [8]. According to a study about the characteristics of climate change around the Erhai Lake area, the annual average air temperature will keep a warming trend in the future [79], which may continuously accelerate the evaporation of the lake and, in turn, alter the local water balance and energy budget allocation. The "E_original" line means lake evaporation is calculated by using original data. The "E_detrend R s " line means lake evaporation is calculated by using detrend R s and original T a . The "E_detrend T a " line means lake evaporation is calculated by using detrend T a and original R s . Variations of partial derivatives of (d) T a and (e) R s . (f) Relative contributions of R s and T a to the annual total evaporation trend based on the results of the detrending method and sensitivity method.

Conclusions
This study is aimed to evaluate the nine evaporation methods at different timescales on the basis of the comprehensive measurements conducted by an EC system over the Erhai Lake and revealed the long-term trend of evaporation. The accuracy of methods varied with timescales resulted from the (dis)agreement between its parameterization scheme and the drivers of lake evaporation. Under the default condition, the Dalton-based RyH method is the most accurate method at short timescales (daily to weekly), while the MT method resulted in a large overestimation of lake evaporation. On a monthly scale, the variation of lake evaporation is better described by the combination and solar radiation-based methods since evaporation is mainly correlated to air temperature and radiation (R n and R s ). Under the calibrated condition, the combination methods using the hysteresis heat storage produced comparable results to the default condition using measured heat storage as input, indicated the hysteresis model could be considered as a good alternative for estimating heat storage without the measurement of water temperature profiles. The JH and RyH method improved slightly after optimization. The large negative bias of the default Mak method and positive bias of the default MT method were substantially reduced after optimization. The calibrated JH method had a relatively small annual bias, and a simple data requirement was applied in the estimates of long-term evaporation. The annual total evaporation showed an insignificant increasing trend of 0.30 mm year −1 during 1981-2018, with a minimum occurred in the early 1990s. Relative contributions of meteorological factors were quantified by the detrend method and sensitivity method. Rising air temperature contributed to the increasing trend of evaporation by 69.2% and 69.5% based on the detrend method and sensitivity method, respectively. This result indicated the amplified impact of a warming climate on the water cycle over the Erhai Lake basin, which will help to better understand the response of the lake system to climate change. In this study, we only concentrated on the estimate of historical evaporation. The extent of lake evaporation over Erhai Lake in the future and its variability are underexplored and needed to discuss further.