Mountain Lake Evaporation: A Comparative Study between Hourly Estimations Models and In Situ Measurements

Abstract

The quantification of evaporation is very important for water resource management, determination of water availability and environmental modeling. Evaporation (E) is a key component of the hydrological cycle, which becomes more relevant under a changing climate scenario, where increases in temperature and E are projected on a regional scale. This study evaluated the hourly applicability of 21 methods to estimate E in a mountain lake. Models were grouped as combination, solar radiation-temperature and mass transfer methods. Estimated E by all models were compared to three measurement campaigns with an eddy covariance system in Lake Laja during the month of January (22 to 28), March (17 to 23) and November (09 to 15) of 2016. The average evaporated water during the measurement sessions were 3.40, 3.38 and 1.89 mm d⁻¹, respectively. Best model performance was obtained with models whose main E principles are heat flow and water vapor flow. The best performance in this group was the Penman model using a calibrated wind function with a determination coefficient (R²) of 0.91, Nash–Sutcliffe coefficient (NS) of 0.9, and index of agreement (W) of 0.98. Comparing daytime and nighttime hours, it was found that the daytime E in the lake is best explained by the product of the vapor pressure deficit and wind speed, while the nighttime E was explained by the flow of heat in the water. The results highlight the importance of the analysis of diurnal dynamics of water flux and energy stored in water to better understand the E in water bodies.

1. Introduction

The effects of climate change and rising temperatures throughout the world, and particularly in arid and semi-arid regions, are threatening to reduce available surface water through enhanced surface evaporation, especially in reservoirs [1]. Mountains and agricultural regions located in snowmelt-dominated Mediterranean climate watersheds have been identified as highly vulnerable to climate change. Temperature levels in mountain regions are expected to increase by 3–4 °C [2,3], thus increasing E rates in mountain lakes.

Evaporation is an important source of water loss in lakes and reservoirs in arid and semi-arid regions. The Mediterranean climate or dry summer climate exists in the Mediterranean basin, some parts of South Africa, California, Australia, and Chile [4]. Loss of water by evaporation can exceed industrial [5] and municipal water use [6]. Martínez et al. [5] determined that water loss due to evaporation from farm dams in the southeast of Spain was greater than the industrial water demand and equivalent to 27% of the domestic water use for a population of two million inhabitants. Similarly, Gökbulak y Özhan [7] estimated that annual evaporation from lakes and dams in Turkey is greater than the human and industrial consumption and in some regions of Australia, 40% of the water stored in reservoirs is lost through evaporation [8,9]. Bouin et al. [10] found that extreme wind events in a French lake can triple the evaporation, reaching values higher than 350 W m⁻². Miranda et al. [11] reported evaporation measurements in a lake with Mediterranean climate obtaining average daily values of 3.7 mm d⁻¹, 4.0 mm d⁻¹, 4.5 mm d⁻¹ y 2.5 mm d⁻¹ during June, July, August and September, respectively.

A sequence of uninterrupted dry years has affected central Chile (30–38° S) since 2010 with annual rainfall deficits ranging from 25% to 45% [12]. This phenomenon, which has been referred to as the mega-drought due to its unprecedented duration and large spatial extent, has coincided with the warmest decade on record in this area of the country [13], increasing the severity of the drought due to higher E loss from the reservoirs [12,14].

The Laja Lake is a natural reservoir located at 1360 m above sea level (37°10′ S, 71°22′ W), in the area affected by the mega-drought. It has been operated as an interannual regulation system since 1960. At that time, the El Toro Hydroelectric Plant was built, beginning their activities in 1973, transforming the Laja Lake into a multipurpose reservoir system. More recently, a new demand was added to the system by the construction of an irrigation canal (Laja-Diguillín canal), located downstream of the lake in the Laja River. In general, water demand has increased over time not only due to the incorporation of the Laja-Diguillín Canal, which started transferring water from the Laja River in 2001, but also due to climate variability and population growth. In this sense, precipitation in the study area has decreased due to climate change [15] and, therefore, supply has decreased, while increasing temperatures are associated with increased evaporation and, in turn, greater water needs for agriculture. Similarly, population growth has increased energy demand, resulting in an increased need for lake water [16]. The lake has a temperate snowmelt and rainfall-dominated headwater system (the Laja Lake basin). In addition, it has had a tunnel that supplies water from the Polcura River since 1978 [17].

Since El Toro began operating in 1973, there were no important conflicts between users in the first three decades, mainly because the lake levels were adequate to meet water demands. However, conflicts have occurred between users over the past ~15 years and have increased in number over the past decade, probably enhanced due to the 2010–2019 drought in Chile [12,16]. Due to the conflicts between the users of the lake water and the impact of climate change, it has been necessary to reformulate the Laja’s operating rules. Establishing a maximum volume for both irrigation and power generation according to the available water volume [18].

One of the main problems in water resource management, such as water reservoirs, lakes and lagoons, and particularly in the case of the Laja Lake, is the estimation of all the components of the water budget. Such knowledge is necessary for the regulation of water supply, the estimation of potential percolation and for the water quality analysis of the reservoir. In water reservoirs, generally, it may be relatively easy to measure inflow, outflow and precipitation. However, to identify percolation through the reservoir floor requires an accurate estimation of E loss. Therefore, E from water bodies has been an important topic in many studies during the years [19].

In Chile, agricultural water demand has been commonly studied by the determination of potential evapotranspiration. In this sense, some studies have been carried out for specific crops [20,21,22,23], but E from water bodies has been poorly studied. Under climate change scenarios, the ability to accurately measure and estimate surface evaporation becomes even more crucial for adequate water management. However, there might be considerable uncertainty in hydrological projections associated with the model structures and parameterization [24,25].

The eddy covariance (EC) technique is one of the few methods to obtain direct E measurements and it is considered one of the most accurate and reliable methods to determine E rates and the main components of the surface energy balance from water bodies [19,26,27]. However, its use is limited due to the expensive and complex nature of its operation. Therefore, efficient modeling approaches are required.

A wide variety of methods to estimate E have been reported in the literature. They can be classified into several types, including mass balance, energy budget models, the aerodynamic approach (based on water vapor flux), solar radiation and temperature models (based on solar radiation and air temperature), and combination equations (based on heat flux and water vapor flux). Most of the equations were developed for specific studies and areas. Therefore, the methods are applicable to areas with similar conditions to those in which they were developed [28].

Several studies have evaluated the performance of E models at daily (or larger) scales [19,28,29,30,31], but to our knowledge, hourly performance of E models has been less studied. In this sense, it is important to note that normally all the environmental variables required for these models are easily measurable on sub-daily time intervals, which would allow the models to be applied at smaller time scales. Granger and Hedstrom [32] developed relationships between hourly E rates and parameters that affect E and concluded that wind speed was one of the most important factors controlling E. McGloin et al. [33] evaluated three methods for modeling sub-daily E in a small reservoir in southeast Queensland, Australia, determining that the theoretical mass transfer model showed the best correlation.

Reliable estimates of surface energy balance fluxes are often needed for representing the effect that different land surfaces have on regional climate [34]. To accurately model diurnal dynamic variability in the atmospheric boundary layer (ABL) it is crucial that environmental processes are simulated frequently (i.e., time-steps of one hour) ([33,34,35,36]). Estimating this variable in nearly real time is important for both agricultural and hydrological management, while monitoring long-term trends allows the identification of the climatological impacts on the global hydrosphere [37]. The analysis of daytime E allows the association of this phenomenon with specific events that occur during the day, e.g., wind gusts, heat waves, cloudiness, among others. Hourly E rates can also be compared to surface energy balance models that use remotely sensed images (satellite or unmanned aerial vehicles) to estimate latent heat for land and water surfaces.

The main objectives of this study were to evaluate the performance of 21 models to estimate hourly E in water bodies and to identify the main environmental variables that control its diurnal dynamics. These estimates are compared with measurements from an EC system during three field campaigns in a mountain lake. Models are evaluated according to their original approach, and then adjusted to the conditions of the Chilean Laja Lake by calibrating wind functions.

2. Materials and Methods

2.1. Study Area

The Laja Lake is located in central-southern Chile, in the middle of the Andes mountain range (37°10′–37°27′ S and 71°22′–71°14′ W) at a height of 1360 m above sea level. It is 33 km long and 3 km wide, with an average depth of 75 m, and an approximate surface of 70 km². The climatic patterns in the area are generated mainly from frontal rainfall and are enhanced due to the orographic effects of the Andes mountains [38]. In this area, snowfall is concentrated mainly during June to August and accounts for 47% of the total annual precipitation. The Laja Lake is a snow-dominated hydrological system and has a snow-dominated hydrological regime where runoff is strongly regulated by the lake.

The Laja Lake sub-basin (975 km²) is conformed by Andean volcanic soils with low permeability [39]. The soils are mainly volcanic and sandy with the presence of steep slope topography. The lake has two outlets: a natural seepage through fractured rocks at the west side of the lake, and a Hydroelectric Power Plant tunnel (El Toro), which connects the lake with the Polcura River. The flows are controlled and managed according to energy and irrigation demands [16].

The Laja Lake is a natural reservoir operated as an interannual regulation system. It supplies the Abanico, El Toro, Antuco, Rucúe and Quilleco hydroelectric plants, as well as the irrigation system Laja sur, Zañartu, Laja Diguillin and a series of other canals [16]. The Laja Lake is a natural water body but regulated in a controlled way by the El Toro Hydroelectric Plant. It has a volume of approximately 5600 Hm³, which is used for agricultural purposes and as the cornerstone of hydroelectric supply from the interconnected power system in Chile.

The measurement site was located on the southern shore of the lake (Figure 1), 4.5 km from the Antuco ski center (37°21′43.86″ S, 71°20′24.60″ O a 1360 m.a.s.l). To the south of the measurement site is the Antuco volcano, and to the north is a hilly area with little vegetation. The soils in the area adjacent to the lake are sandy and stony, not suitable for agriculture, livestock or forestry, and with limited use for wildlife, recreation or the protection of river basins. The site location was selected considering previous information on the direction of the prevailing winds in the lake. This previous analysis was considered to ensure that the fluxes of the surface energy balance measured by the EC system comes mainly from the body of water. Wind speed in the measurement time period reached 8 m s⁻¹, predominantly blowing in a direction from west to east.

2.2. Equipment and Measurements

An EC system was installed in the study area to measure the main components of the energy budget and the main environmental variables that affect E at 30 min intervals in January (from 22 to 28), March (from 17 to 23) and November (from 9 to 15) of 2016. Specific days during field campaigns were selected considering the weather forecast to find high probabilities of clear and sunny days during measurements.

A pyranometer (LI-200R, Licor, Lincoln, NE, USA) installed at a height of 3.0 m, was used to measure solar radiation (Q_s). Net radiation (Q*) was measured with a net radiometer (NR-Lite2, Kipp & Zonen, Delft, The Netherlands) installed at a height of two meters, using a 2.0 m extension arm. Sensible heat flux (H_s) was measured with a 3D sonic anemometer (CSAT3, Campbell Scientific, Logan, UT, USA) installed at a height of 2.0 m and placed horizontally to the wind direction and a fine thermocouple (FW3, Campbell Scientific, Logan, UT, USA). Latent heat flux (λE) was measured by a sonic anemometer and an infrared gas analyzer (EC150, Campbell Scientific, Logan, UT, USA).

Air temperature and relative humidity were measured with a probe (HMP45C, Campbell Scientific, Logan, UT, USA) installed at a height of 2.5 m. Water surface temperature was measured with 8 temperature recorders (iButton DS1921Z, Maxim Integrated, San Jose, CA, USA) at 0.1, 0.2, 0.3, 0.4, 0.5, 1.5, 3 and 4.5 m depth, all attached to a marine rope with a weight and affixed with a buoy. Wind speed was measured with the same 3D sonic anemometer. Signals from these instruments were sampled every 5 s and stored as the mean value after 30 min in a data logger (CR3000, Campbell Scientic, Logan, UT, USA). Two 12 V batteries connected in parallel supplied the power.

To determine water heat flux (Q_x), the reduced form of the energy budget was used following [27,31,33,40,41]:

$$Q_x=Q^{*}-\lambda E-H_s$$

(1)

where, $Q_x$, $Q^{*}$, $\lambda E$ and $H_s$ is in W m⁻².

In order to ensure that latent heat and sensible heat fluxes came mainly from the water body, a footprint analysis was performed using the method developed by Hsieh et al. [42].

2.3. Evaporation Estimation Methods

Twenty-one models were used to estimate evaporation on an hourly basis, as indicated in the following groups:

Combination equations: The combined method generally includes energy transfer and aerodynamic components. For practical applications, the formulation requires measurement of radiation fluxes on the water body, water surface temperature, and meteorological data such as wind speed, air temperature, and dew point temperature at a given height.

Solar radiation and temperature: This group consists of four models, and has the advantage of requiring fewer physical parameters, such as solar radiation, air temperature, atmospheric pressure, and relative humidity. However, as in all the other methods, measurements should be representative of the conditions of the water body.

Mass transfer: This group consists of three models. Mass transfer is based on physical principles. Evaporation calculated from this method is considered proportional to a function of the mean wind speed, and the difference between air vapor pressure at the evaporation surface and at some level above the surface of the water body.

Table 1. Evaporation estimation methods (mm h⁻¹).

Method	Equation
Combination group
Priestley-Taylor	$E=\alpha \left(\frac{\Delta }{\Delta +\gamma }\right)\left(\frac{Q^{*}-Q_x}{\lambda {\rho }_w}\right) 3.6$
deBrauin-Keijman	$E=\left(\frac{\Delta }{0.85\Delta +0.63\gamma }\right)\left(\frac{Q^{*}-Q_x}{\lambda {\rho }_w}\right) 3.6$
Penman-Monteith	$E=\frac{\Delta \left(Q^{*}-Q_x\right)+{\rho }_a{C}_p\left(e_s-e_a\right)/r_a}{\Delta +\gamma \left(1+\frac{r_s}{r_a}\right)}\frac{3.6}{\lambda {\rho }_w}$
Penman	$E=\left(\frac{\Delta }{\Delta +\gamma }\left(Q^{*}-Q_x\right)+\frac{\gamma }{\Delta +\gamma }{K}_u·f\left(U\right)\left(e_s-e_a\right)\right)\frac{3.6}{\lambda {\rho }_w}$
Penman-Kimberly	$E=\left(\frac{\Delta }{\Delta +\gamma }\left(Q^{*}-Q_x\right)+\frac{\gamma }{\Delta +\gamma }{K}_u\left(e_s-e_a\right)W_f\right)\frac{3.6}{\lambda {\rho }_w}$
Brutsaert-Stricker	$E=\left(\left(2\alpha -1\right)\frac{\Delta }{\Delta +\gamma }\left(Q^{*}-Q_x\right)-\frac{\gamma }{\Delta +\gamma }{K}_u·f\left(U\right)\left(e_s-e_a\right)\right)\frac{3.6}{\lambda {\rho }_w}$
deBruin	$E=\left(\frac{\alpha }{\alpha -1}\right)\left(\frac{\gamma }{\Delta +\gamma }\right)\frac{f\left(U\right)\left(e_s-e_a\right)}{\lambda {\rho }_w}\ast 3.6$
Solar radiation—temperature group
Jensen-Haise	$E=C_r\left(T_a-T_x\right)Q_s\frac{3.6}{\lambda {\rho }_w}$; $C_r=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\left(C_1+C_2·C_H\right)$}\right.$; $C_H=\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$\left(e_s-e_a\right)$}\right.$; $C_1=38-\left(\frac{2Elev}{305}\right)$; $T_x=-2.5-1.4\left(e_s-e_a\right)-Elev/550$; $C_2=7.3{ }^{\circ }\mathrm{C}$
Makkink	$E=\left(0.61\left(\frac{\Delta }{\Delta +\gamma }\right)Q_s-0.12\right)\frac{3.6}{\lambda {\rho }_w}$
Stephens-Stewart	$E=\left(0.0082T_a-0.19\right)\left(Q_{s}x3.495x{10}^{-2}\right)/24$
Turc	$E=\left(0.013\frac{T_a}{T_a+15}\left(Q_s+50\right)\right)/24$ for $RH>50\%$ and for $RH<50\%$ $E=\left(0.013\frac{T_a}{T_a+15}\left(Q_s+50\right)\left[1+\frac{\left(50-RH\right)}{70}\right]\right)/24$
Mass transfer group
Ryan-Harleman	$E=\left(\frac{\left(2.7{\left|{\theta }_v\right|}^{0.333}+3.1U_2\right)\left(e_{s*}-e_a\right)}{\lambda {\rho }_w}\right)3.6$
Trivett	$E=0.024\left(e_{s*}-e_{a*}\right)U_2$/24
Quinn	$E=\left(0.052+0.0066U_3\right)\left(e_{s*}-e_{a*}\right)U_3/24$
(A)	$E=a·\left(e_{s*}-e_{a*}\right)$
(B)	$E=a·U_2\left(e_{s*}-e_{a*}\right)$
(C)	$E=a·\left(1-e^{-U_2}\right)\ast \left(e_{s*}-e_{a*}\right)$
(D)	$E=a·\left(1+b{U}_2\right)·\left(e_{s*}-e_{a*}\right)$
(E)	$E=a·U_2·\left(e_{s*}-e_{a*}\right)·\left(1-b·\left(T_a-T_0\right)\right)$
(F)	$E=a·{\left(T_a+25\right)}^2·\left(100-HR\right)$
(G)	$E=a·\left(1+b·U_2\right)·\left(e_{s*}-e_{a*}\right)·\left(1-c·\left(T_a-T_0\right)\right)$
3.6 and 1/24 multipliers appear in the equations to convert output to mm h−1
∆	Slope saturated vapor pressure-temperature curve (Pa °C−1).
α	1.26, Priestley-Taylor empirically derived constant, dimensionless.
a,b y c	Adjustment parameters for mass models.
γ	Psychrometric constant (depends on temperature and atmospheric pressure (Pa °C−1)).
θv	Difference in temperature between surface water and air (°C).
Λ	Latent heat of vaporization (MJ kg−1).
${\rho }_w$	Water density (998 kg m−3).
ea	Vapor pressure at the air temperature (kPa) (deBruin in mb).
ea*	Vapor pressure at the dewpoint temperature (mb).
es	Saturation vapor pressure at the air temperature (kPa) (deBruin in mb).
es*	Saturation vapor pressure at the water surface temperature (mb).
$Elev$	Elevation (m).
$f\left(U\right)$	Wind function, $f\left(U\right)=a+b U_2$; with a = 0.5 y b = 0.54 for Penman, a = 1 and b = 0.54 for Brutsaert-Stricker, a = 2.9 and b = 2.1 for deBruin model.
$J$	Day of the year.
$K_u$	Unit conversion factor = 74.44 W m−2 kPa−1 = 0.268 MJ m−2 h−1
RH	Relative humidity (%)
Q*	Net radiation (W m−2).
Qs	Solar radiation (W m−2) (Turc in cal cm−2 d−1).
Qx	Water heat flux (W m−2).
$r_s$	Bulk surface resistance (s m−1); $r_s=0$
$r_a$	Aerodynamic resistance (s m−1)
T0	Water surface temperature (°C)
Ta	Air temperature in °C (Stephens-Steward in °F).
$U_2$	Wind speed at 2 m above water surface (m s−1) (km hr−1 for Trivett).
U3	Wind speed at 3 m above water surface (m s−1); $U_2=U_1{\left(\raisebox{1ex}{$Z_2$}\!\left/ \!\raisebox{-1ex}{$Z_1$}\right.\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$7$}\right.}$, $U_2$ is wind speed at height level 2, $U_1$ is wind speed at height level 1, $Z_2$ is measurement height level 2 and $Z_1$ is measurement height level 1
$W_f$	Wind function, $W_f=a_w+b_w{U}_2$; $a_w=0.4+1.4exp\left\{-{\left[\left(J-173\right)/58\right]}^2\right\}$; $b_w=0.605+0.345exp\left\{-{\left[\left(J-243\right)/80\right]}^2\right\}$; For latitudes south of the equator, one should use $J^{\prime }$ in place of $J$, where $J^{\prime }=\left(J-182\right)$ for $J\ge 182$ and $J^{\prime }=\left(J+182\right)$ for $J<182$

shows the equations involved for each of the 21 models. The last 7 models of the mass transfer group were proposed by Singh and Xu [29], and require the estimation of the parameters a, b and c. The parameters were estimated using a method based on the least squares criterion between measured and estimated evaporation. For empirical models, the values of the parameters need to be calibrated from a sample of input and output values of the model, and from an objective function whose value must be minimized. One of the simplest objectives is a linear regression model that uses an input variable x, and an output variable y. The model to be calibrated is an equation of type y’ = A + BX in which parameters A and B can minimize the objective function ${\sum }_{i=1}^n{(y^{\prime }-y)}^2$, where $n$ is the number of data pairs.

The Penman, Brutsaert–Stricker and deBruin models consider the following wind function:

$$f\left(U\right)=a+b{U}_2$$

(2)

where a and b are adjustment constants, which will be calibrated to maximize the function of the coefficient of determination (R²).

2.4. Evaluations

All models were evaluated using the root mean square error (RMSE), the coefficient of determination (R²) (0–1), the Nash–Sutcliffe (NS) index (−∞–1) and the index of agreement (W), (0–1).

Table 2. Equations and indicators to validate the models.

Indicator	Equation
Coefficient of determination	$R^2=\frac{cov{\left(o,m\right)}^2}{sd\left(o\right)sd\left(m\right)}$
Nash-Sutcliffe	$NS=1-\frac{{\sum }_{i=1}^n{\left(o_i-m_i\right)}^2}{{\sum }_{i=1}^n{\left(o_i-\overline{o}\right)}^2}$
Index of Agreement	$W=1-\frac{{\sum }_{i=1}^n{\left(o_i-m_i\right)}^2}{{\sum }_{i=1}^n{\left(\left|m_i-\overline{o}\right|+\left|o_i-\overline{o}\right|\right)}^2}$
Root mean square error	$RMSE=\sqrt{\frac{{\sum }_{i=1}^n{\left(o_i-m_i\right)}^2}{n}}$

cov = covariance, sd = typical deviation, o = observed values, m = values of the model, n = number of data.

shows the equations corresponding to these indicators.

3. Results

The components of the water surface energy budget and the environmental variables that affect evaporation were measured for seven days in January (from 22 to 28), March (from 17 to 23) and November (from 9 to 15) of 2016. The data were filtered, eliminating the values that were attributable to noise in the data, sensor failure or rainy periods, resulting in a total of 747 half hour data (259, 255 and 233 in January, March and November, respectively).

3.1. Environmental Components

Figure 2 shows the environmental variables measured.

Table 3. Daily maximum (Max), minimum (Min) and mean ($\overline{X}$) values of the environmental variables measured.

DOY	T (°C)			T0 (°C)			U2 (m s−1)			HR (%)
	Max.	Min.	$\overline{\mathit{X}}$	Max.	Min.	$\overline{\mathit{X}}$	Max.	Min.	$\overline{\mathit{X}}$	Max.	Min.	$\overline{\mathit{X}}$
22	24.3	14.5	18.7	18.1	16.9	17.5	7.3	1.4	4.5	82.5	22.5	51.4
23	21.4	13.5	17.7	18.3	16.5	17.3	9.1	1.0	5.0	88.5	31.5	50.0
24	12.6	9.7	10.7	17.6	15.0	16.3	7.6	4.2	5.5	98.8	74.0	87.2
25	14.1	10.1	12.3	16.8	16.3	16.5	6.7	1.3	4.4	88.9	62.2	70.3
26	16.1	11.0	13.2	17.1	16.3	16.9	7.3	2.5	5.2	81.7	50.8	66.0
27	13.5	9.3	12.1	17.4	16.5	16.7	7.1	4.0	5.2	97.9	62.2	72.5
78	21.2	3.4	120	17.5	15.5	16.6	5.2	0.8	3.1	93.9	9.2	43.5
79	21.6	7.1	14.3	17.6	15.8	16.3	4.0	1.4	2.7	32.5	16.5	24.8
80	22.5	7.4	14.7	17.4	15.9	16.3	6.2	0.9	3.1	46.3	13.7	30.4
81	20.7	5.3	13.3	16.5	15.6	15.9	6.6	0.5	2.7	74.5	14.3	48.0
82	23.9	7.2	14.2	17.6	15.5	15.9	6.2	0.4	2.3	65.7	15.0	42.7
315	14.5	2.8	9.7	-	-	-	7.4	0.8	3.7	71.5	23.0	51.1
316	17.9	5.1	11.2	14.3	10.4	11.6	7.9	0.8	3.5	70.3	17.4	47.4
317	15.0	4.3	11.4	12.6	10.1	11.2	8.6	0.7	4.8	71.1	30.9	50.6
318	12.5	7.5	10.1	12.0	9.7	10.6	7.5	3.3	5.6	85.9	42.4	59.5

shows the maximum, minimum and daily mean values of the variables. Maximum wind speed ($U_2$) values were recorded during the afternoon hours (between 4 and 7 p.m.), reaching 9.1 m s⁻¹, 6.6 m s⁻¹ and 8.6 m s⁻¹ in January, March and November, respectively (Figure 2a). Vapor pressure deficit ($VPD$) (Figure 2c) and air temperature ($T_a$) (Figure 2b) reached their lowest values during the early morning hours of the day and their highest values in the afternoon hours (between 3 and 5 p.m.). Relative humidity ($RH$) (Figure 2e) ranged between 10% and 90%, reaching its maximum and minimum values in morning and afternoon hours, respectively, showing a behavior contrary to the vapor pressure deficit (VPD) curve. Maximum values for relative humidity were recorded in January. The maximum temperature difference between water surface ($T_0$) and air temperature (Figure 2f) occurred between six and eight in the morning, reaching up to 14 °C. The daily variation of water surface temperature (Figure 2d) did not exceed ±2 °C, reaching its maximum during the afternoon hours. In November, there was a failure in the sensor that recorded water surface temperature, so there are no data for the 315 day of the year (DOY) and half of the day 316.

Figure 3 shows the measured solar radiation ($Q_s$). The maximum values were recorded between 2 and 3 p.m., reaching 1200 W m⁻², 900 W m⁻² and 1050 W m⁻² in January, March and November, respectively.

Footprint

The footprint analysis shows that 90% of the fluxes, under unstable conditions, came from an area 500 m away from the station, while the most important contributions came from less than 100 m of distance. The fluxes in the study area predominantly come from the northwest, this assures that all the surface energy balance measurements mainly came from the surface of interest. Figure 4 shows the flux direction under unstable conditions (north west direction), distance (Xpeak) to the point of maximum influence from the measurement point (generally less than 100 m) and the flux footprint distribution function.

3.2. Surface Energy Budget Components

Figure 5 shows the daily curves for the surface energy budget components.

Table 4. Daily maximum (Max) and mean ($\overline{X}$) values of the surface energy budget components.

DOY	$\mathit{\lambda }\mathit{E} (\mathbf{W} {\mathbf{m}}^{-2})$		${\mathit{Q}}^{*} (\mathbf{W} {\mathbf{m}}^{-2})$		${\mathit{H}}_{\mathit{s}} (\mathbf{W} {\mathbf{m}}^{-2})$		${\mathit{Q}}_{\mathit{x}} (\mathbf{W} {\mathbf{m}}^{-2})$
	Máx	$\overline{\mathit{X}}$	Máx	$\overline{\mathit{X}}$	Máx	$\overline{\mathit{X}}$	Máx	$\overline{\mathit{X}}$
22	184	91	722	207	37	2	591	113
23	217	117	862	181	56	2	674	62
24	183	104	221	6	92	57	99	−155
25	159	102	534	78	55	40	345	−64
26	178	120	669	139	60	42	472	−23
27	189	113	508	59	87	46	315	−100
78	209	121	608	100	128	38	324	−58
79	182	114	593	109	97	14	479	−20
80	205	99	602	114	41	3	556	12
81	202	85	599	113	130	24	473	4
82	191	71	580	108	145	10	549	28
315	104	45	665	190	157	12	559	133
316	129	56	686	190	54	−2	587	136
317	117	61	704	202	42	6	576	134
318	103	63	606	151	62	22	458	67

shows the daily maximum and mean values of all components. The net radiation reached a maximum of 862 W m⁻² at 2:30 p.m. in January, recording negative values at nighttime. Furthermore, 73% of the observations of the sensible heat flux were positive, indicating that there is heat transfer from the lake to the atmosphere. Negative values of sensible heat flux were observed during the afternoon and evening hours, indicating that there is heat transfer from the atmosphere to the lake. Water heat fluxes were calculated according to Equation (1), and had positive values during daytime measurements, while nighttime measurements recorded negative values. The latent heat flux reached a maximum value of 217 W m⁻² at 3:00 p.m. on January 23.

The following variables were correlated with the measured latent heat flux to identify those that most significantly affect E: the relationships between latent heat flux (measured every 30 min) and wind speed above the lake surface (Figure 6a), vapor pressure deficit (Figure 6b), the product between wind speed and vapor pressure deficit (Figure 6c), net radiation over the lake surface (Figure 6d) and water heat flux (Figure 6e), recording R² values of 0.32, 0.29, 0.43, 0.20 and 0.06, respectively. Figure 6 shows that the lowest relationship is given by $\lambda E$ and $Q_x$, followed by $\lambda E$ and $Q^{*}$. It can be observed that λE is more influenced by the product of $VPD$ and $U_2$ (R² = 0.43).

During these campaigns, the average evaporated water depths reached 3.40, 3.38 and 1.89 mm per day in January, March and November, respectively.

Table 5. Evaporation (mm d⁻¹) and estimated evaporated flow (m³ s⁻¹).

Date	Evaporated Water Height (mm d−1)	Evaporated Flow (m3 s−1)
22-01-2016	3.20	2.6
23-01-2016	4.02	3.3
24-01-2016	2.44	2.0
25-01-2016	3.56	2.9
26-01-2016	3.95	3.2
27-01-2016	3.23	2.6
18-03-2016	4.14	3.4
19-03-2016	4.01	3.3
20-03-2016	3.46	2.8
21-03-2016	2.97	2.4
22-03-2016	2.33	1.9
10-11-2016	1.53	1.2
11-11-2016	1.93	1.6
12-11-2016	1.88	1.5
13-11-2016	2.21	1.8

shows the evaporated water depths for the days measured, and its equivalence as E flow (m³ s⁻¹), which was calculated based on the area flooded by the lake (approximately 70 km²). The lowest E rate was observed in November, while the highest value was recorded in January.

Figure 7 shows the volume of daily water evaporated according to the elevation of the water level of the lake for different depths of evaporated water. This indicates that the higher the elevation of the lake, the greater the volume of water lost by evaporation. An important change occurs between elevation 1360 and 1370 m because one of its branches is active, which implies a significant jump in the curves of evaporated water volume.

On the other hand, if the lake were at its maximum level (1390 m), at an evaporation rate of 4 mm d⁻¹, the volume of water evaporated could reach 480,000 m³ d⁻¹.

3.3. Model Comparisons with Measurements

Figure 8, Figure 9 and Figure 10 show the relationship between the E obtained with the hourly models and the E measured with the EC system.

As observed in Figure 8, the Penman and Penman–Monteith models show similar relationships (R²), indicating the similarity between the models. The same occurs with the Priestley–Taylor and deBrauin–Keijman models. The best relationship between measured and estimated evaporation in this group was obtained with the Priestley–Taylor model.

The wind functions of the Penman, deBruin, Brutsaert–Stricker models were adjusted according to equation (2). The adjustment was performed to maximize the function of the coefficient of determination (R²) as:

$$f{\left(U\right)}_{Penman}=0.395+0.171U_2$$

(3)

$$f{\left(U\right)}_{deBruin}=0.020+0.491U_2$$

(4)

$$f{\left(U\right)}_{Brutsaert–Striker}=0.0+0.0U_2$$

(5)

With the adjustment, R² changed from 0.84 to 0.91 in Penman, 0.43 to 0.44 in deBruin and 0.14 to 0.84 in Brutsaert–Stricker.

As for the solar radiation and temperature group, the data were filtered only considering the daytime ranges where solar radiation exceeded 100 W m⁻².

All models in this group overestimate E (Figure 10). The efficiency indicator (NS) for all models in this group is negative, indicating that the mean of the observed data has a higher predictive power than the model. A similarity of measured and estimated E relationships is observed between the Stephens–Stewart and Jensen–Haise models. This is explained by the similarity in their equations. Based on this, the best model of this group is the Stephens–Stewart model because it presents the lowest overestimation.

Similarly, the mass transfer models present an acceptable linear relationship, same as most of their indicators, except for the NS indicator. As mentioned before, water surface temperature data of two days were not recorded and, therefore, these models were evaluated with less data (n = 659).

Figure 11 shows the relationship between estimated values with three models of the mass transfer group and measured values. The best linear relationship is given by the Ryan–Harleman model. However, this does not mean that it is the best model, since it presents an overestimation of E values (72%). Of these three models, the best is Trivett’s as it has the best indicators and the least overestimation.

The parameters a, b and c of the generalized form of the mass transfer model developed by Sing and Xu [29], were adjusted through a least-squares fitting

Table 6. Parameters a, b and c of the mass transfer model and goodness-of-fit indicators.

Model	Parameter			Indicator
	a	b	c	R2	NS	W	RMSE
(A)	0.2471	-	-	0.23	0.51	0.77	0.05
(B)	0.0012	-	-	0.64	0.70	0.92	0.04
(C)	0.0057	-	-	0.24	0.45	0.77	0.05
(D)	0.0019	0.4570	-	0.69	0.79	0.93	0.03
(E)	0.0012	0.0002	-	0.63	0.70	0.92	0.04
(F)	7.648·10−7	-	-	0.27	0.48	0.87	0.05
(G)	0.0018	0.4549	0.0009	0.69	0.80	0.93	0.03

. The best adjustments are given by the models B, D, E and G (Figure 12).

According to the analyzed indicators, the G model is the best model of the generalized form of the mass transfer models. It is important to note that the D model requires the estimation of one parameter less, which is an advantage when processing the data. In addition, the difference obtained in the NS index was not significant.

As many of the equations for estimating E have constants or coefficients that were determined by those who developed or modified them to meet their particular needs, adjustment factors ($K_{LL}$) for each of the hourly models were obtained with the least squares model. The models to estimate E are expressed as follows:

$$E=f\left(E\right)\ast K_{LL}$$

(6)

where $f\left(E\right)$ is the E function of each of the models evaluated.

Table 7. Statistical indicators prior to and after the application of adjustment parameters for the Laja Lake ($K_{LL}$) (* indicates that the displayed factor is the one obtained after applying $K_{LL}$).

Group	Method	N	KLL	R2	NS	NS *	W	W *	RMSE (mm h−1)	RMSE * (mm h−1)
Combination	Priestley-Taylor	747	0.98	0.84	0.82	0.83	0.96	0.96	0.03	0.03
	deBrauin-Keijman	747	0.95	0.83	0.80	0.82	0.95	0.95	0.03	0.03
	Penman-Monteith	747	1.23	0.84	0.68	0.83	0.91	0.96	0.04	0.03
	Penman	747	0.69	0.84	−0.03	0.83	0.84	0.96	0.07	0.02
	Penman WF	747	0.96	0.91	0.90	0.91	0.98	0.98	0.02	0.02
	Penman-Kimberly	747	0.53	0.79	−2.49	0.77	0.67	0.94	0.13	0.03
	Brutsaert-Stricker	747	0.91	0.14	−1.61	−1.59	0.56	0.56	0.11	0.11
	Brutsaert-Stricker WF	747	0.81	0.84	0.60	0.83	0.92	0.96	0.04	0.03
	deBruin	747	0.47	0.43	−4.18	0.15	0.56	0.79	0.16	0.10
	deBruin WF	747	2.50	0.44	−1.19	0.02	0.54	0.78	0.10	0.07
Solar radiation-temperature	Jensen-Haise	348	0.36	0.22	−6.33	0.06	0.45	0.59	0.28	0.10
	Makking	348	0.39	0.10	−4.96	0.29	0.45	0.60	0.25	0.09
	Stephens-Stewart	348	0.50	0.16	−1.83	0.28	0.57	0.64	0.17	0.09
	Turc	348	0.35	0.12	−7.28	0.23	0.41	0.42	0.29	0.09
Mass transfer	Ryan-Harleman	659	0.61	0.67	−0.11	0.77	0.85	0.93	0.08	0.04
	Trivett	659	0.73	0.69	0.71	0.77	0.95	0.94	0.04	0.04
	Quinn	659	0.64	0.64	0.54	0.65	0.94	0.91	0.05	0.04

shows the statistical indicators evaluated for conditions prior to and after the application of K_LL. For both conditions, it is observed that the best model is Penman including the wind function $f\left(U\right)$ (Penman WF). For the models of the group of solar radiation and temperature, the best model was the Stephens–Stewart model even after applying $K_{LL}$. Finally, the best model of the mass transfer models was Trivett, both prior to and after the application of $K_{LL}$.

4. Discussion

4.1. Measurements

In the southern hemisphere, January and March correspond to summer months. January is the month with the highest evaporative demand in the year, while November corresponds to a spring month. Days 22 and 23 were clear days, and days 24 to 27 were cloudy, which caused a drop in the air temperature recorded on those days compared to previous days (Figure 2). The variation in water surface temperature between a cloudy and a clear day for the period between days 22 and 27 reached an average of 0.8 °C. Between days 24 and 27, it was observed that the water surface temperature was higher than the air temperature, indicating that there was heat transfer from the lake surface to the near atmosphere. The decrease in air temperature caused a decrease in $e_a$, which reduced the capacity of the atmosphere to store water vapor, as reflected in the decrease in E around those days. The water surface temperature was lower in November, considering snow-rain driven regime. The decrease in water temperature in November can be associated with the beginning of thaws.

Air temperature dropped sharply during nighttime due to the absence of solar radiation. Water temperature had an average daytime variation of 1.3, 1.6 and 2.9 °C in January, March and November, respectively. Therefore, water temperature was higher than that of air during the night, causing heat transfer from the lake water to the atmosphere. For the night period, water heat flux was negative. There was a positive but not significant correlation between ∆T and Hs (R² = 0.22).

The latent heat flux for a water surface is expected to be high and even greater than that of a crop. However, the effect of water depth is considerable and influences the heat storage capacity of the water body. The maximum latent heat flux (217 W m⁻²) was recorded on a clear day in January at 3:00 p.m., a similar value was found in [10]. This peak coincides with the highest wind speed recorded during the measurement time period (9.1 m s⁻¹) but occurred later than the maximum solar radiation (2:30 p.m.). Daytime measurements of latent heat flux had a similar pattern, reaching the highest values in the afternoon hours. The $\lambda E$ reaches average values of 108, 98 and 56 W m⁻² during the months of January, March and November, respectively. Miranda et al. [11] measured E using an EC system in a Mediterranean climate reservoir (Alqueba, Portugal) during summer months, they found for the hottest month a $\lambda E$ value of 113 W m⁻² similar to those found for the Laja Lake.

Nighttime E rates (between 8:00 p.m. and 7:00 a.m.) corresponded to an average of 36.9, 40.7 and 26.4% of the daytime evaporation rates recorded in January, March and November, respectively. This can be explained by the energy release from the water body. Similar rates were found by Liu et al. [41] at the Ross Barnett dam, Mississippi. This lake has an area of 133 km² and water depths between 4–8 m. Measurements with an EC system showed that nighttime evaporation represents 45% of total E.

The heat flux in the water presented the lowest correlation when compared to latent heat. However, when analyzing the daytime and nighttime ranges separately (Figure 13), it can be observed that nighttime measurements showed a significant correlation (R² = 0.64) with $\lambda E$, while variables were inversely proportional (between 8:00 a.m. and 7:00 p.m.), i.e., as water heat flux decreases, latent heat flux increases. As for the daytime range (between 7:30 a.m. and 7:30 p.m.), the relationship between latent heat flux and heat flux in water was not significant (R² = 0.01). On the other hand, when analyzing the other components for the daytime range, it is observed that E is better correlated during these hours of the day, due to the product between vapor pressure deficit and wind speed (R² = 0.51).

Large and high-altitude lakes have a high capacity to store energy. In summer, most of the solar radiation is absorbed by the water body and the greatest heat storage capacity is recorded in late summer and early fall. Afterwards, the stored energy begins to release. In late fall and early winter, when solar radiation is small, this energy is released to the atmosphere as latent and sensible heat fluxes; the lowest stored heat is recorded in spring, when the heating process of the water body starts again [43]. This explains the decrease in the nighttime E rate recorded in November (26.4%).

In the Laja Lake, there is a low relationship between latent heat flux and net radiation. This is mainly due to the depth of the lake; short-wave radiation penetrates the lake surface where the energy is absorbed. Granger and Hedstrom [32] studied this relationship for a 25 m deep lake with hourly measurements and found an R² coefficient of 0.01, an even lower value than that in the present study.

The relationship between wind speed and latent heat flux was greater than the relationship between net radiation and evaporation fluxes. Same as in [32], it is observed that E in the Laja Lake is controlled by wind speed above the lake, similar observations were found in [10]. In addition, there is a directly proportional relationship between the latent heat flux and the product between vapor pressure deficit and wind speed. This occurs because evaporation rate depends directly on vapor pressure deficit (VPD) and is enhanced by higher wind speeds.

During measurement in January and March, evaporated water reached an average of 3.4 mm per day, equivalent to 266 million L. If we consider a daily water consumption of 200 L per person, the evaporated volume is equivalent to the consumption of 3644 people in a year. By expressing E in terms of flux, E flux averaged 2.76 m³ s⁻¹ for the measurement days in January and March, and 1.52 m³ s⁻¹ for those in November. Considering that the Laja Lake receives contributions from the Polcura River with flow rates that averaged 3.47 m³ s⁻¹, 2.34 m³ s⁻¹ and 8,9 m³ s⁻¹ in January, March and November, respectively, it can be indicated that at least half of the water entering through the Polcura tunnel is lost by E.

4.2. Models

The evaluation of the combined group models showed that the best estimate was given by the Priestley–Taylor model. This model shows similarity in its equation with the deBrauin–Keijman model. Both models incorporate the variable γ, which is dependent on elevation. However, the deBrauin–Keijman model was developed at sea level, so the Priestley–Taylor model was more suitable for a mountain lake. Similar results were obtained by Wang et al. [31]. The Penman–Kimberly model overestimates E in the Laja Lake. This model is based on Penman and presents a modification of the wind function. Therefore, the wind function proposed by the Kimberly model was not adequate to estimate the hourly evaporation in the Laja Lake.

One of the lowest yields was given by the deBruin model. The author combined the Penman and Priestley–Taylor equations, and obtained a formula that predicts E. The model does not include net radiation, but the drying power of air (a function of the wind associated with vapor pressure gradient), and an energy component represented by the latent heat. deBruin’s experimental data showed that the model performed well in a shallow body of water and in periods greater than 10 days [43]. The poor performance of this model in the Laja Lake can be associated with the exclusion of $Q^{*}$ and $Q_x$ in its formula. As abovementioned, $Q_x$ plays an important role in nighttime E.

The Brutsaert–Stricker model also presented a low fit with the measured data. Same as the models in its group, it bases its estimates on a radiative and an aerodynamic component. Brutsaert–Stricker [44] determined that the model was relatively insensitive to the choice of wind function, which would result in a poor performance of the model when applied to the Laja Lake. The wind function was adjusted for the Penman, deBruin and Brutsaert–Stricker models in order to maximize R². The adjusted R² increased from 0.84 to 0.91; from 0.43 to 0.44; from 0.14 to 0.84 in Penman, deBruin and Brutsaert–Stricker, respectively. The adjustment of this last model resulted in the elimination of the aerodynamic component (a and b = 0 in the wind function).

The models of the solar radiation and temperature group overestimated E in the lake. These models base predictions on solar radiation and temperature, excluding net radiation, water heat flux, and the aerodynamic component. The overestimation can be explained by the effect of not including lake depth since part of the short-wave radiation is absorbed at the water surface, but the fraction not absorbed by the surface can penetrate deeper layers [41]. Due to the high heat capacity of the water, an important amount of energy is stored in the water during spring and summer months. Normally the water temperature increases and the heat storage peaks occurs in late summer or early fall [41,45].

The Ryan–Harleman, Trivett and Quinn models have adjustment constants calibrated for daily estimates. Due to this, it was necessary to divide the output of these models into 24 parts to obtain mm h⁻¹. The results obtained show that the Trivett model had the best fit. Additionally, the least squares method was applied to determine the constants a, b and c in the generalized form of the mass transfer models proposed by Sing and Xu [29]. Once the constants were determined and the model applied, it was obtained that the best fit was given by the G model. However, it is important to note that the D model presents similar fit indicators and requires fewer parameters. Although the models in this group tend to overestimate evaporation, they produced a good fit with the measured data. This occurred due to the influence of the wind on the E of the lake.

Figure 14 shows the distribution of the model groups with respect to the indicators used. For RMSE versus R², the closer to zero on the abscissa axis and to one on the ordinate axis, the better fit of the model. As for NS and W versus R², a better fit is close to one on both axes. The models of the combined group had the best fit, followed by the mass transfer models. Conversely, the models based on solar radiation and temperature exhibited the worst performance.

Figure 15 shows maximum, minimum, and average absolute errors between modeled and observed values for hourly estimates. It is observed that the mean value of most of the combination models and the mass transfer group was found within an error range of ± 0.04 mm h⁻¹, while the models of the solar radiation and temperature group present a significant overestimation. The mean values of the absolute errors of the generalized form of mass transfer models exhibited the smallest distance to the observed value.

The use of the adjustment factors (K_LL) resulted in a general improvement in all the models for the RMSE, W and NS indicators. However, the best fit was obtained using the combined group, followed by the mass transfer and the temperature and solar radiation group. Rosenberry et al. [30] evaluated 15 models to estimate evaporation in a small mountain lake on a daily scale, and determined that the best comparisons with the measured E were achieved with the methods that include available energy and aerodynamic components. Similar results were found in [31].

Figure 16 shows maximum, minimum and average absolute errors after applying the adjustment factor K_LL. The mean value of the absolute errors improved, achieving an error range of ±0.02 mm h⁻¹ for all the models. As in the previous case, the maximum amplitude of the error is found in the group of the solar radiation and temperature models, while the minimum amplitudes are achieved with the combined group. This shows the effectiveness of the use of an adjustment constant for the particular conditions of the lake.

5. Conclusions

The main goals of this study were to evaluate the performance of 21 models to estimate hourly E rates in water bodies and to identify the main environmental variables that control its diurnal dynamics. Hourly E and environmental variables that affect it were measured in a mountain lake for 18 days using an EC system. Average E rates measured during three field campaigns reached 3.40, 3.38 and 1.89 mm d⁻¹ in January, March and November 2016, respectively.

Best model performance was obtained with the models that include the available energy and aerodynamic components (combined group). The best fit was obtained with the Priestley–Taylor equation; however, the Penman model showed the best performance when the wind function was adjusted. The models based on solar radiation and air temperature exhibited the worst performance because the aerodynamic and available energy components were not considered, the Stephens–Stewart model produced the best estimates of the group. The mass transfer models presented a good agreement with the E measured. The greatest advantage of this group is that it does not require water heat flux measurements, the G model had the best performance in this group.

In order to adjust the results to the particular conditions of the Laja Lake, an adjustment factor (K_LL) was applied to all models, with the exception of those proposed by Sing and Xu [29] as they present their own adjustment parameters. After applying the adjustment factor, the Penman model (with the wind function) presented the best performance. In general, all the models improved their indicators after applying the K_LL.

Nighttime E represented 36.9 and 40.7% of the total daily rates in January and March, respectively, which is explained by the energy release from the water body to the atmosphere. The water heat flux during nighttime accounted for 64% of E. In November, the nighttime component reduced E rates to 26.4%, mainly due to the energy loss from the lake in winter. Daytime E is aerodynamically driven in the Laja Lake and hourly fluxes are strongly correlated with the product of wind speed above the lake and the vapor pressure deficit.

The analysis of the nighttime component of latent heat flux and the better correlation with wind speed found during daytime highlights the importance of the analysis of diurnal dynamics of water flux and energy stored in water. Therefore, hourly measurements of latent heat flux, net radiation, sensible heat flux, and water heat flux are recommended to better understand the E variability and energy balance components of water bodies.