Estimating Turbulent Fluxes in the Tropical Andes

The correct estimation of Sensible Heat Flux (H) and Latent Heat Flux (LE) (i.e., turbulent fluxes) is vital in the understanding of exchange of energy and mass among hydrosphere, atmosphere, and biosphere in an ecosystem. One of the most popular methods to measure these fluxes is the Eddy Covariance (EC) technique; however, there are a number of setbacks to its application, especially in remote and topographically complex terrain such as the higher altitudes of the Andes. Efforts have been made by the scientific community to parameterise these fluxes based on other more commonly measured variables. One of the most widespread methods is the so-called bulk method, which relates average temperature, humidity, and wind vertical profiles to the turbulent fluxes. Another approach to estimate LE is the Penman-Monteith (PM) equation which uses meteorological measurements at a single level. The objective of this study was to validate these methods for the first time in the Tropical Andes in Southern Ecuador (in the páramo ecosystem at 3780 m a.s.l.) using EC and meteorological measurements. It was determined that the bulk method was the best to estimate H, although some adjustments had to be made to the typical assumptions used to estimate surface meteorological values. On the other hand, the PM equation yielded the best LE estimations. For both fluxes, the error in the estimations was within the uncertainty range of the EC measurements. It can be concluded that it is possible to accurately estimate H and LE using the methods described in this paper in this ecosystem when no direct measurements are available.


Introduction
The turbulent fluxes (e.g., the Latent Heat Flux (LE) and Sensible Heat Flux (H)) are crucial elements of an ecosystem's energy balance. The estimation of these atmospheric processes at the land surface has long been recognized as vital in the understanding of exchange of energy and mass among hydrosphere, atmosphere, and biosphere [1]. One of the most popular methods to measure turbulent fluxes is the eddy covariance technique, which consists of measuring the 3D field of the wind, the air temperature, pressure, humidity, and its water vapour content at high frequencies (typically 10-50 Hz). Once these variables are known, the latent heat flux and the sensible heat flux can be derived [2,3].
One of the major advantages is that this method provides continuous measurements that are spatially averaged [4]. However, the eddy covariance method has some restrictions that limit its applicability: (i) it requires constant supervision and maintenance, (ii) measurements in complex terrain are challenging because of the theoretical assumptions to transform the high-frequency measurements

Study Area
The study area is located in the Zhurucay Ecohydrological Observatory at 3870 m a.s.l. in Southern Ecuador, on the divide of waters draining into the Pacific Ocean and through the Amazon River to the Atlantic Ocean on the Andean Cordillera (03 • 03'45.12" S, 79 • 14'04.96" W). Annual rainfall is 1345 mm and displays low seasonality, with a relatively dry period from June to August [23]. Maximum seasonal rainfall values are related to the Inter Tropical Convergence Zone [24]. The mean annual temperature is 6 • C (max. = 14.2 • C, min. = 0.4 • C), mean wind speed is 2 m s −1 [25], the relative air humidity is 93.6%, and the prevalent wind direction is from the northeast.
The vegetation within the catchment is typical of páramo grasslands: composed of tussock grass (mainly Calamagrostis sp.), which covers more than 70% of the surface with heights of between 30 and 80 cm [24], and cushion plants (such as Plantago rigida, Xenophyllum humile, and Azorella spp.), which cover around 25% of the surface [26]. The sensors are located in an area where the tussock grass is predominant. More information about the study site can be found in the cited studies in this subsection and several others [27][28][29]. Figure 1 shows the study location and a landscape view of the observatory.

Instrumentation
The turbulent fluxes were measured by means of an enclosed path LI-COR Eddy Covariance system (LI-7200) and a Gill WindMaster 3D sonic anemometer. These sensors are installed at 3.6 m above the surface. At a short distance, a separate automatic weather station was installed to measure the basic meteorological variables. Air temperature and relative humidity were measured by means of a Campbell Scientific CS-215 combined probe with radiation shield. A Met-One 034B Windset anemometer was used to measure wind speed and direction and the atmospheric pressure was measured with a Vaisala PTB110 sensor. The measurement height of all sensors was 1.8 m above the surface.
The four radiation components were measured with a Kipp & Zonen CNR4 radiometer installed at 0.5 m above the surface, which is associated with a 3 m radius footprint. This height was chosen in order to cover vegetation only within the sensor footprint. Therefore, both the fence surrounding the Atmosphere 2020, 11, 213 3 of 15 site and the small creek down-slope did not disturb the measurements. Another CNR4 installed also at 0.5 m a few meters away was used to double check the measurements, and they were comparable (data not shown).

Data
The turbulence data were collected at a sampling rate of 20 Hz and averages, variances, and fluxes were calculated every 30 minutes using the EddyPro software (version 6.2.0, LI-COR, Lincoln, NE, USA). This also provides a quality control of the data using quality flags given by the sensors, plausibility limits, and spike removal. Specific corrections were also applied, such as planar fit adjustment, co-spectra filtering, time-lag compensation, and humidity dependent spectrum correction. More information about these corrections can be found in [24].
Measurements from the weather station were sampled every 15 seconds and averages were recorded every 5 minutes, and these data were used to calculate the 30-min averages after quality control was performed. Both data sets were one-year long (March 2017-Feb 2018), with a total of 17,502 half-hourly observations.
Atmosphere 2020, 11, x FOR PEER REVIEW 3 of 15 at 0.5 m a few meters away was used to double check the measurements, and they were comparable (data not shown).

Data
The turbulence data were collected at a sampling rate of 20 Hz and averages, variances, and fluxes were calculated every 30 minutes using the EddyPro software (version 6.2.0, LI-COR, Lincoln, NE, USA). This also provides a quality control of the data using quality flags given by the sensors, plausibility limits, and spike removal. Specific corrections were also applied, such as planar fit adjustment, co-spectra filtering, time-lag compensation, and humidity dependent spectrum correction. More information about these corrections can be found in [24].
Measurements from the weather station were sampled every 15 seconds and averages were recorded every 5 minutes, and these data were used to calculate the 30-min averages after quality control was performed. Both data sets were one-year long (March 2017-Feb 2018), with a total of 17,502 half-hourly observations. Before performing any calculations, we filtered the data to obtain a high-quality data set. From the eddy covariance data-set, missing, low-quality (based on the flag given by the EddyPro software for H and LE, values with a flag value higher than 1 were not considered), or implausible values (LE > 600 W m −2 , T < −50 °C, T > 50 °C) were filtered out. The same was done for the weather station data set. After filtering, we selected only measurements for which data was available for both datasets, leaving 11,492 observations (66%), which is an acceptable value for eddy-covariance measurements when compared to other flux sites [30].

Methods
The sensible and latent heat flux (H and LE, respectively) are defined as: Before performing any calculations, we filtered the data to obtain a high-quality data set. From the eddy covariance data-set, missing, low-quality (based on the flag given by the EddyPro software for H and LE, values with a flag value higher than 1 were not considered), or implausible values (LE > 600 W m −2 , T < −50 • C, T > 50 • C) were filtered out. The same was done for the weather station data set. After filtering, we selected only measurements for which data was available for both datasets, leaving 11,492 observations (66%), which is an acceptable value for eddy-covariance measurements when compared to other flux sites [30].

Methods
The sensible and latent heat flux (H and LE, respectively) are defined as: where ρ is the near-surface air density, c p is the specific heat capacity of air, (ω'θ') s is the vertical kinematic turbulent heat flux, (ω'q') s is the vertical kinematic turbulent moisture flux, and L v is the latent heat of vaporization of water. The subscript s stands for surface values. The fluxes are defined positive towards the surface.

Bulk Method
The bulk method is based on the assumption that the turbulent fluxes (H and LE) can be related to the vertical gradients of temperature, specific humidity, and wind speed between the surface and one level in the atmosphere. The advantage of using the surface as one of the levels is that the gradients become relatively large. Large gradients reduce the sensitivity and measurement errors when compared to calculating the gradient from measurements made at two levels in the atmosphere. The turbulent fluxes for this method are expressed as follows: where u* is the friction velocity, θ* friction temperature, and q* friction specific humidity, and are defined as: where V is wind speed, z is the height of the measurements (i.e., 1.8 m), θ is potential temperature, κ the dimensionless Von Karman constant (0.41), z om is the momentum or aerodynamic roughness length, which is defined as the height above the surface where wind speed theoretically becomes zero, and z oh is the roughness length for heat. We assumed that the roughness lengths for heat and moisture are the same, which is commonly done [31]. We excluded the displacement height from the calculations (Equations (6) and (7)) as we assumed it to be zero. To calculate density, we used the calculations of saturated vapour pressure and water vapour pressure, which are discussed below.
Ψ m and Ψ h , in Equation (5) to (7), are the vertically integrated stability functions for momentum and heat. The turbulence flux values can then be found via iteration assuming for the first iteration neutral conditions so that z/L = 0. For the definition of the stability functions and the method of iteration we followed [32].
The surface temperature was estimated from the outgoing long-wave radiation data using Equation (8): where OLR is outgoing long-wave radiation, σ is the Stefan-Boltzmann constant, ε is the emissivity of the surface, and T s is the surface temperature in Kelvin.

of 15
We calculated specific humidity (q) based on temperature, relative humidity, and pressure as follows: where, p sat is saturation vapour pressure, T is air temperature, RH is relative humidity, p v is actual vapour pressure, p d is the pressure of dry air, p is pressure, R d is the gas constant for dry air, R v is the gas constant for water vapour, p s0 is the saturation vapour pressure at T 0 = 273.16 K, and w is the mass mixing ratio of water vapour to total air. This method was used to calculate specific humidity at 1.8 m and at the surface. For surface humidity estimations we considered two options. First, we assumed that the surface is always saturated, hence relative humidity RH s = 100% at all times. The second assumption was to limit the relative humidity gradient to a maximum value during daytime. For this second approach, the assumption that RH s = 100% was retained during the night, because the air is more humid as there is no incoming shortwave radiation, and during most nights the RH measured at 1.8 m is also close to saturation (mean RH at night during the study period was about 97%).
To calculate the surface specific humidity a second method was applied as well, following the study from [33]. In this relation, the surface specific humidity can be calculated from the surface temperature: where r c is the stomatal resistance of the vegetation. We calculated the roughness length for momentum (z o ) using only data with near neutral conditions so that the wind speed profiles are theoretically logarithmic. Neutrally stable conditions were assumed to occur when the temperature difference between T and T s was smaller than or equal to +/− 0.2 • C. The momentum roughness length was then calculated assuming a logarithmic wind profile: which yielded an average value of z om = 0.077 m.
To calculate z oh we used two different approaches. First, we assumed it to be constant throughout the entire study period. In this case, it was defined as z oh = z om /7.4 (following the equations from [34]). When z 0h is calculated in this manner, it is assumed that the roughness length for heat only depends on the land surface characteristics, as is the case for the momentum roughness length.
In the second approach, the dependence of the value of z oh on u * (and thus on the conditions of the flow) is considered, following the recent work by [35] and [36]. The calculation was performed as follows: with kB −1 defined as: where the roughness Reynolds number is Re * = z om u * v −1 , in which v is the kinematic viscosity.

The Penman-Monteith Equation
The latent heat flux can also be calculated by means of the Penman-Monteith equation. This approach does not use the vertical profiles as input but combines the energy balance closure and aerodynamic characteristics. As formulated in Equation (17), the method includes all the parameters that govern energy exchange and the corresponding latent heat flux from a uniform extension of the so-called reference crop. The reference crop resembles 0.12 m, well-watered, actively growing grass. We calculated reference evapotranspiration (ET o ) following the procedure described by [37]: where R n is net radiation, G is soil heat flux (estimated to be 10% of R n during daytime and 50% during night-time [37]), and V is wind speed. Details about the calculation of γ, ∆, e 0 (T) and e a can be found in [37]. This method therefore only necessitates routine meteorological measurements as an input (temperature, solar radiation, relative humidity, wind speed, and atmospheric) and avoids the need for assumptions about the values of the variables at the surface. After this, we calculated actual evapotranspiration (ET a ) by multiplying ET 0 by the value of k c , which serves to transform ET 0 to ET a . The value of k c is first taken as a constant (0.9) according to [24]. However, in the same paper it is reported to vary between 0.87 in the drier period and 0.92 in the wetter period for this study site. Finally, we multiplied this value by L v to obtain the latent heat flux.

The Bulk Method
In this section, the first results shown correspond to H and LE calculated by means of the bulk method using the parameters ε = 1, q s = q sat (T s ), and z om = 0.077 m from the logarithmic wind profile assumption) and constant z oh equal to z om /7.4. Note that the iterations did not converge for wind speed values smaller than 1.5 m s −1 . Therefore, we only made the calculations for the time steps with measured wind speed higher than 1.5 m s −1 . These results are shown in the scatter plots in Figure 2, where we compare the bulk flux results to those from the turbulence measurements (correlations were found to be significant with p < 0.01).
Atmosphere 2020, 11, x FOR PEER REVIEW 6 of 15 so-called reference crop. The reference crop resembles 0.12 m, well-watered, actively growing grass. We calculated reference evapotranspiration (ETo) following the procedure described by [37]: where Rn is net radiation, G is soil heat flux (estimated to be 10% of Rn during daytime and 50% during night-time [37]), and V is wind speed. Details about the calculation of γ, Δ, e 0 (T) and ea can be found in [37]. This method therefore only necessitates routine meteorological measurements as an input (temperature, solar radiation, relative humidity, wind speed, and atmospheric) and avoids the need for assumptions about the values of the variables at the surface. After this, we calculated actual evapotranspiration (ETa) by multiplying ET0 by the value of kc, which serves to transform ET0 to ETa. The value of kc is first taken as a constant (0.9) according to [24]. However, in the same paper it is reported to vary between 0.87 in the drier period and 0.92 in the wetter period for this study site. Finally, we multiplied this value by Lv to obtain the latent heat flux.

The Bulk Method
In this section, the first results shown correspond to H and LE calculated by means of the bulk method using the parameters ε = 1, qs = qsat(Ts), and zom = 0.077 m from the logarithmic wind profile assumption) and constant zoh equal to zom/7.4. Note that the iterations did not converge for wind speed values smaller than 1.5 m s −1 . Therefore, we only made the calculations for the time steps with measured wind speed higher than 1.5 m s −1 . These results are shown in the scatter plots in Figure 2, where we compare the bulk flux results to those from the turbulence measurements (correlations were found to be significant with p < 0.01).  Absolute Bias (MAB) is close to 100 W m −2 . These errors might be related to the assumptions made for the calculation of surface temperature leading to an underestimation of the temperature gradient during daytime and an overestimation of the gradient during night-time. This combination would result in the observed systematic underestimation in sensible heat flux calculations.
In contrast, LE estimations ( Figure 2b) are less strongly (R 2 = 0.60 compared to R 2 = 0.81 for H) correlated to the observations. There is heavy overestimation of the higher values of LE as reflected in the slope of the regression line being much higher than 1. The poor quality of the calculations is also reflected in the value of MAB close to 200 W m −2 . These large overestimation for higher values of LE could be related to the assumptions made for the calculation of surface humidity, leading to either an overestimation of the temperature gradient or of the specific humidity gradient. The implications of the assumptions for the estimation of H and LE are addressed in the next subsection.

Alternative Assumptions in the Application of the Bulk Method
As found in the previous subsection, the bulk method results deviate considerably from those of the Eddy Covariance (EC) method. Most likely this is related to the assumptions that we initially made to estimate the values of T s and q s at the surface. Therefore, we study a variety of different approaches for the estimation of different parameters in this subsection.

Emissivity (ε)
With ε = 1, as was assumed in the previous subsection, the surface is a black body. Values in the range 0.95-0.975 have been found to be reasonable in the past for natural surfaces [38][39][40]. In this study we found that the optimum value for the páramo is 0.95, as the surface temperature shows the most realistic values (higher values led to freezing temperatures during most nights, which is not the case in the area). Taking a lower ε value will result in higher values of estimated T s . Therefore, ∆T will be smaller during the night, when negative H is observed, and this will result in less negative H estimations at night. On the other hand, ∆T will be higher during the day and thus H will also be higher. This solves the systematic underestimation problem.
The second row of Table 1 shows the result of applying the bulk method for ε = 0.95. It is observed that H estimations improved significantly, MAB decreased from 99.74 W m −2 to 32.05 W m −2 (this is below the uncertainty range of the H measurements in the study site, as reported by [24]), and R 2 increased from 0.81 to 0.90. However, the quality of LE estimations decreased (MAB increased and R 2 slightly increased). This was expected as higher temperatures lead to higher specific humidity values in an exponential relationship. Hereafter, we used ε = 0.95 for our calculations, i.e., every next row in Table 1 also show the results using ε = 0.95.

Surface Humidity
We used two alternatives to calculate the surface specific humidity estimations, and thus ∆q. First, we used the approach presented in Equation (13) [33]. This method was applied only during daytime (06H00-18H00), as the assumption of saturation at the surface seems realistic during night-time. The major drawback of this approach is that it requires a value for the stomatal resistance of the vegetation cover (r c ) as input, which is not known for the páramo tussock grasslands. However, we used two values found in the literature to test the potential impact that this refinement in the calculations could have in the estimation of LE.
First, we used 10 s m −1 , as representative of green well-watered grasslands [37]. As observed in Table 1, the values for H did not change because this only impacts the humidity gradient, whereas the MAB for LE decreased by a factor of more than 2; however, R 2 decreased. Secondly, we used the value estimated by [41] for tussock grasslands in New Zealand (200 s m −1 ). Even though the vegetation is similar, those estimations were based on measurements of only two days during daytime under dry conditions, which are not representative of our humid study site. The results of this estimation were very poor, LE was largely overestimated, and the correlation was almost zero. Even though these estimations did not yield accurate results, it is important to note that LE appears to be quite sensible to r c , and hence, it shows potential for when the stomatal resistance for this landscape is known.
The second approach we used is based on the idea that the assumption of a permanently saturated surface is unrealistic, especially during daytime (06H00-18H00). This is supported by the observations, which show that RH values as low as 20% (data not shown) are possible and that the soil can get so dry that even soil moisture measurements below the wilting point have been recorded, as is shown later in this section. Thus, we hypothesize that during the daytime the assumption of a saturated surface might not be realistic, causing unrealistically high values of ∆q and hence overestimation of LE. In order to mimic this, we limited the daytime relative humidity gradient (∆RH) to a maximum of 5% and 10% (Table 1) and kept the assumption of saturation at the surface during night-time. In both cases, the value of MAB was considerably reduced and R 2 also improved. The best results were obtained when ∆RH max was set to 5%.

Thermal Roughness Length
Alternatively, we followed [35], who found a relation between the values of z 0h and u * . The only unknown value to apply this method is the constant C 1 (Equation (16)), which depends on the emissivity of the surface and on the vegetation cover. We calculated the turbulent fluxes using values for C 1 in a range between grass and shrub lands, which is representative of our study site. The results in Table 1 show that both H and LE calculations improved via this method. The MAB for H decreased substantially and did not vary too much for the chosen values. It also vastly improved the LE estimations, and the respective MAB values showed a wider range for this flux; R 2 was not affected by changes in C 1 .

The Best Possible Estimation
Based on the results from this section, optimising for best possible results for H and LE was achieved by applying the bulk method using ε = 0.95, ∆RH max = 5%, and C 1 = 0.7. The results of these calculations are summarized in Table 1 and shown in Figure 3. The results for H are quite good for both positive and negative values, with the only exception being a few positive observations that were estimated to be close to zero (Figure 3). In the case of LE, although it was still systematically overestimated, the MAB was for the first time below 100 W m −2 . In addition, the correlation was the best among all the methods used (all correlations are found to be significant with p < 0.01).

The Penman-Monteith (PM) Equation
As shown in Figure 3, the estimation of LE (Figure 3b) was not as good as the one for H (Figure 3a) using the bulk method even when the alternatives to the assumptions were implemented. The main issue when calculating LE via the bulk method seems to be the estimation of the q gradient. In contrast to this, the PM equation requires only measurements from a weather station and thus has the strength that no assumptions are needed to estimate the equation inputs. This approach outperformed the previously shown estimations of LE: MAB is significantly lower (30.52 W m −2 ) and R 2 is higher than in the previous cases ( Figure 4).

Monthly Analysis
In order to gain more insight into our results and to determine which variables control and lead to the best estimation of H and LE, we analysed the calculations on a monthly time scale. The summary is shown in Table 2, in which a seasonal pattern can be recognised. When compared to precipitation and soil moisture during the study period ( Figure 5), estimations seem to be better (i.e., the best combinations of high R 2 and low MAB) in general for the wet months and worse during the dry months. In particular, the MAB is lower for wet months compared to dry months. When we focus on the turbulent fluxes as estimated by means of the bulk method, March is the month when both estimations have a low MAB, while the correlation is high. In contrast, the estimations are especially erroneous for LE in November (highest MAB of all estimations), when soil moisture reaches its minimum value ( Figure 5). Interestingly, for LE estimations based on the Penman-Monteith equation, the estimations were also the poorest (Table 2) during the months of November and December when the soil is drier ( Figure 5). It should be noted that soil moisture from June to October is nearly constant, whereas precipitation decreases compared to the March-April and January-February periods. This could be explained by the fact that precipitation in the form of drizzle (a common precipitation type in the páramo [23]) and fog are not measured by the tipping bucket rain gauge. This "occult" precipitation could contribute to maintain the soil moisture at the levels observed in

Monthly Analysis
In order to gain more insight into our results and to determine which variables control and lead to the best estimation of H and LE, we analysed the calculations on a monthly time scale. The summary is shown in Table 2, in which a seasonal pattern can be recognised. When compared to precipitation and soil moisture during the study period ( Figure 5), estimations seem to be better (i.e., the best combinations of high R 2 and low MAB) in general for the wet months and worse during the dry months. In particular, the MAB is lower for wet months compared to dry months. When we focus on the turbulent fluxes as estimated by means of the bulk method, March is the month when both estimations have a low MAB, while the correlation is high. In contrast, the estimations are especially erroneous for LE in November (highest MAB of all estimations), when soil moisture reaches its minimum value ( Figure 5). Interestingly, for LE estimations based on the Penman-Monteith equation, the estimations were also the poorest (Table 2) during the months of November and December when the soil is drier ( Figure 5). It should be noted that soil moisture from June to October is nearly constant, whereas precipitation decreases compared to the March-April and January-February periods. This could be explained by the fact that precipitation in the form of drizzle (a common

Monthly Analysis
In order to gain more insight into our results and to determine which variables control and lead to the best estimation of H and LE, we analysed the calculations on a monthly time scale. The summary is shown in Table 2, in which a seasonal pattern can be recognised. When compared to precipitation and soil moisture during the study period ( Figure 5), estimations seem to be better (i.e., the best combinations of high R 2 and low MAB) in general for the wet months and worse during the dry months. In particular, the MAB is lower for wet months compared to dry months. When we focus on the turbulent fluxes as estimated by means of the bulk method, March is the month when both estimations have a low MAB, while the correlation is high. In contrast, the estimations are especially erroneous for LE in November (highest MAB of all estimations), when soil moisture reaches its minimum value ( Figure 5). Interestingly, for LE estimations based on the Penman-Monteith equation, the estimations were also the poorest (Table 2) during the months of November and December when the soil is drier ( Figure 5). It should be noted that soil moisture from June to October is nearly constant, whereas precipitation decreases compared to the March-April and January-February periods. This could be explained by the fact that precipitation in the form of drizzle (a common precipitation type in the páramo [23]) and fog are not measured by the tipping bucket rain gauge. This "occult" precipitation could contribute to maintain the soil moisture at the levels observed in Figure 5b, as hypothesised by [42]. Furthermore, the period when fog occurs more frequently is from June through August (while almost no fog is observed during November and December). Additionally, the highest values of solar radiation were observed in November and December (Figure 5c). The combination of these factors could explain the patterns observed in Figure 5a,b. Table 2. Monthly values for the mean absolute bias (MAB is in W m −2 ) and R 2 when applying the best possible estimation and the Penman-Monteith method.

Best Estimation
Penman-Monteith We found that the differences in LE might be explained best by precipitation ( Figure 5a) and soil moisture (Figure 5b). March was the wettest month, as can be seen from the precipitation and soil moisture data. The opposite was the case for November, which was the driest month in our study period. The systematic overestimation of the bulk method LE throughout the year stands out in November. Based on precipitation and soil moisture, we can hypothesize that these variables are related to this larger error during this month. As observed in Figure 5b, November is the month where the soil is the driest and below the wilting point (which was found to be 0.5 [44] for the páramo grasslands). The wilting point is the humidity at which the roots of the plants cannot uptake water from the soil pores anymore. Therefore, the larger errors in this month could be explained by the fact that there is no water available to evaporate. Therefore, the bulk method overestimated the latent heat flux, as it does not consider water availability. The systematic overestimation of ∆q was also more pronounced during this month because it was the warmest one (data not shown), and the error in the specific humidity estimations is larger for higher temperatures.
However, during the drier months a higher bias was also observed using the Penman-Monteith method. To discriminate between the wet and dry period, Figure 7 shows a scatter plot with the same data shown in Figure 4 but now for separate datasets for the dry and wet period. The brown points and fitted line correspond to data from the dry conditions (November, December, and January). It is observed that they are clustered above the 1:1 line and that the overestimation of LE is higher for this period (MAB = 60.66 W m −2 ). The data for the wet period is shown using green points and a green fitted line; the line overlaps the 1:1 line. This indicates that the systematic overestimation is restricted to the dry period, and hence, LE is more accurately estimated using this method under wet conditions (MAB decreased to 21.52 W m −2 ).   Figure 6 shows the time series for the first 12 days of March and November. In the case of H (Figure 6a,b), it was observed that the fluxes are larger in November, primarily as a result of the increased solar radiation in the study area during this month [43]. The main errors in the H calculations are made when large negative values occur (i.e., when H < 0), as observed in Figure 6a,b. This causes the MAB to be higher in November than in March (Table 2). LE is shown in Figure 6c,d for March and November respectively. In Figure 6c,d, we observe the characteristic overestimation of LE calculations; nonetheless, this overestimation is more pronounced during November (Figure 6d), and this is why MAB is so much higher in this month compared to the March value.  Figure 6 shows the time series for the first 12 days of March and November. In the case of H (Figure 6a,b), it was observed that the fluxes are larger in November, primarily as a result of the increased solar radiation in the study area during this month [43]. The main errors in the H calculations are made when large negative values occur (i.e., when H < 0), as observed in Figure 6a,b. This causes the MAB to be higher in November than in March (Table 2). LE is shown in Figure 6c,d for March and November respectively. In Figure 6c,d, we observe the characteristic overestimation of LE calculations; nonetheless, this overestimation is more pronounced during November (Figure 6d), and this is why MAB is so much higher in this month compared to the March value. We found that the differences in LE might be explained best by precipitation ( Figure 5a) and soil moisture (Figure 5b). March was the wettest month, as can be seen from the precipitation and soil moisture data. The opposite was the case for November, which was the driest month in our study period. The systematic overestimation of the bulk method LE throughout the year stands out in November. Based on precipitation and soil moisture, we can hypothesize that these variables are  The reason for this overestimation under dry conditions may be related to the fact that we used a constant value for kc. By keeping the value constant, we assume that the proportion of reference evapotranspiration that is actually realized is also constant throughout the year. This might not be the case in drier periods, when water may not be available for different reasons (i.e., dry soil, limited plant transpiration, plant stress). To test this hypothesis, we used the two different values of kc proposed by [24] for the wet and the dry seasons in the same study site. However, this did not improve our estimations significantly (data not shown). This might be due to the fact that the distinction between the wet and dry seasons in the aforementioned publication was based on monthly precipitation, whereas the differences for LE seem to be more related to soil moisture conditions.

Conclusions
The main goal of this project was to determine if Sensible Heat Flux (H) and Latent Heat Flux (LE) can be accurately estimated in the páramo ecosystem located in the high Tropical Andes using data from an automatic weather station. The results were evaluated by comparing them to the measurements from an eddy covariance tower located a few meters away.
The best estimations for H were obtained via the bulk method and a combination of our modifications, i.e., emissivity equal to 0.95 instead of 1, considering zoh as dependent on the conditions of the flow, and limiting the relative humidity gradient to 5%, to reduce the error on dry days. The results of these calculations were significantly and strongly correlated (R 2 = 0.9) with the eddy covariance observations and the mean absolute bias (MAB) below the uncertainty range of the H measurements in the study site [24].
The best estimations of LE were obtained via the Penman-Monteith equation. This method is only dependent on variables which are directly available from the weather station measurements. Hence, assumptions for surface temperature and relative humidity are avoided. The results of these calculations were significantly and highly correlated (R 2 = 0.69) with the eddy covariance observations and the Mean Absolute Bias (MAB) was close to the uncertainty range of the LE measurements in the study site [24]. However, LE was still overestimated, especially under dry conditions. The best results were obtained in March (the wettest month).
All in all, the evidence from this study suggests that the modified bulk method and the Penman -Monteith equation are suitable to estimate the turbulent fluxes in this ecosystem. Future studies are encouraged to focus on dry periods, when the largest errors were found. Furthermore, one method that could improve the results is to consider the stomatal resistance (rc) dependency on the weather The reason for this overestimation under dry conditions may be related to the fact that we used a constant value for kc. By keeping the value constant, we assume that the proportion of reference evapotranspiration that is actually realized is also constant throughout the year. This might not be the case in drier periods, when water may not be available for different reasons (i.e., dry soil, limited plant transpiration, plant stress). To test this hypothesis, we used the two different values of kc proposed by [24] for the wet and the dry seasons in the same study site. However, this did not improve our estimations significantly (data not shown). This might be due to the fact that the distinction between the wet and dry seasons in the aforementioned publication was based on monthly precipitation, whereas the differences for LE seem to be more related to soil moisture conditions.

Conclusions
The main goal of this project was to determine if Sensible Heat Flux (H) and Latent Heat Flux (LE) can be accurately estimated in the páramo ecosystem located in the high Tropical Andes using data from an automatic weather station. The results were evaluated by comparing them to the measurements from an eddy covariance tower located a few meters away.
The best estimations for H were obtained via the bulk method and a combination of our modifications, i.e., emissivity equal to 0.95 instead of 1, considering z oh as dependent on the conditions of the flow, and limiting the relative humidity gradient to 5%, to reduce the error on dry days. The results of these calculations were significantly and strongly correlated (R 2 = 0.9) with the eddy covariance observations and the mean absolute bias (MAB) below the uncertainty range of the H measurements in the study site [24].
The best estimations of LE were obtained via the Penman-Monteith equation. This method is only dependent on variables which are directly available from the weather station measurements. Hence, assumptions for surface temperature and relative humidity are avoided. The results of these calculations were significantly and highly correlated (R 2 = 0.69) with the eddy covariance observations and the Mean Absolute Bias (MAB) was close to the uncertainty range of the LE measurements in the study site [24]. However, LE was still overestimated, especially under dry conditions. The best results were obtained in March (the wettest month).
All in all, the evidence from this study suggests that the modified bulk method and the Penman -Monteith equation are suitable to estimate the turbulent fluxes in this ecosystem. Future studies are encouraged to focus on dry periods, when the largest errors were found. Furthermore, one method that could improve the results is to consider the stomatal resistance (r c ) dependency on the weather conditions. L was found to be very sensitive to different values of r c , which was assumed to be a constant within this study. Research (e.g., its value, its seasonality, and its dependence on the weather and soil conditions) on this parameter is encouraged as this could lead to better quality estimations of q s , and hence LE.