Actual Evapotranspiration Estimates in Arid Cold Regions Using Machine Learning Algorithms with In Situ and Remote Sensing Data

: Actual evapotranspiration (ET a ) estimations in arid regions are challenging because this process is highly dynamic over time and space. Nevertheless, several studies have shown good results when implementing empirical regression formulae that, despite their simplicity, are comparable in accuracy to more complex models. Although many types of regression formulae to estimate ET a exist, there is no consensus on what variables must be included in the analysis. In this research, we used machine learning algorithms—through implementation of empirical linear regression formulae—to ﬁnd the main variables that control daily and monthly ET a in arid cold regions, where there is a lack of available ET a data. Meteorological data alone and then combined with remote sensing vegetation indices (VIs) were used as input in ET a estimations. In situ ET a and meteorological data were obtained from ten sites in Chile, Australia, and the United States. Our results indicate that the available energy is the main meteorological variable that controls ET a in the assessed sites, despite the fact that these regions are typically described as water-limited environments. The VI that better represents the in situ ET a is the Normalized Difference Water Index, which represents water availability in plants and soils. The best performance of the regression equations in the validation sites was obtained for monthly estimates with the incorporation of VIs (R 2 = 0.82), whereas the worst performance of these equations was obtained for monthly ET a estimates when only meteorological data were considered. Incorporation of remote-sensing information results in better ET a estimates compared to when only meteorological data are considered.


Introduction
Arid and semi-arid regions cover approximately 41% of the world's land and are inhabited by more than 2500 million people [1]. These regions are expected to expand because of unsustainable land and water use, as well as a result of climate change, which is exacerbating desertification [1]. In this context, an accurate quantification of evapotranspiration (ET), a relevant hydrological process in arid regions, is important for managing water resources to ensure their availability for human and environmental needs [2][3][4]. Although there have been many efforts to quantify actual evapotranspiration (ET a ) in arid regions [5], few ET a direct observations exist in cold desert climates (also known as arid cold regions). For instance, from a total of 267 sites in the FLUXNET 2015 dataset, only seven sites located in arid cold regions have more than one year of records [6]. This lack of data hinders the understanding of the main processes that drive ET in these environments. Hence, the motivation of this work is to further explore if ET a in arid cold regions is driven by similar or different variables than other climates. Broadly speaking, ET a is mainly driven by energy exchange and water availability, but there are plenty of meteorological and transmitted [22], whereas the Normalized Difference Water Index (NDWI) is sensitive to other properties, e.g., leaf water content [23], and it is also capable of representing both canopy and soil water content [24]. Thus, it is able to represent plant water stress [25]. VIs have several advantages for use in ET a algorithms: they are available from multiple sensors, and they usually change on time scales of weeks to months, so it is feasible to interpolate VI values with observations obtained several days apart, especially if they represent the activity of the vegetation. In addition, VI methods are usually simple and resilient in the presence of data gaps [4]. Moreover, VIs have been applied to natural ecosystems and achieved good results [26,27].
Several studies have demonstrated the potential of combining site-specific ET a data with remotely sensed and meteorological parameters to develop empirical models based on statistical correlations for regional-scale ET a estimates [3,12,14,15,28]. Despite their simplicity, empirical regression formulae can produce ET a values that are comparable in accuracy to more complex models, without as much computational requirements for specific expertise [4]. However, the estimation of ET a with a higher degree of accuracy and over extended time scales has forced researchers to look for techniques such as machine learning [28][29][30][31][32][33][34][35][36]. Elbeltagi et al. [37] applied five machine learning techniques to develop the Combined Terrestrial Evapotranspiration Index. Zhao et al. [30] constructed a physicsconstrained machine learning model to estimate ET a , which conserves the surface energy balance and successfully reproduces extreme values. Torres et al. [29] forecasted potential ET using the multivariate relevance vector machine and limited climatic data. Granata [28] provides some examples of machine learning applications in hydrology and mentions some investigations related to ET a . However, he states that these investigations are limited and that the knowledge on the topic is still partial and fragmented. Moreover, studies that use empirical regression formulae and basic machine learning concepts usually focus on the form of the formulae that predict ET a instead of the factors that drive ET a [4,20].
The aim of this research is to investigate the main variables that control ET a in arid cold regions through implementation of empirical regression formulae using machine learning algorithms. The machine learning algorithms, based on the exhaustive feature selection (EFS) approach, which has not been used before in ET a applications, were formulated first with meteorological data and then with remote sensing data. Consequently, a secondary objective of this work is to analyze if the inclusion of remote sensing data improves ET a estimations. The scope of this work is restricted to spatial extents within the field and landscape scales (between a few hundred and a few thousand square meters) and in arid cold regions, as our review revealed that these locations are underrepresented in the scientific literature.

Study Sites
In this study, we used 10 sites located around the world that, according to the Köeppen climate classification system, correspond to arid cold climate (BSk and BWk) [38]. Arid cold climates are characterized by little precipitation with warm/dry summers and cold/dry winters, as opposed to hot desert climates. Arid cold climates have an annual average temperature below 18 • C, and a threshold that depends on both precipitation and temperature is used to define a region as an arid cold climate [38]. Three of the study sites are located in the Chilean Altiplano, two in Australia, and five in the United States. Figure 1 presents the location of the study sites, and Table 1 shows their main characteristics. Chilean sites are classified as desert cold climates, while the other sites are cold steppe. Additionally, the Chilean sites are located above 4000 m ASL, whereas the sites in Australia and United States are located between 125 and 1530 m ASL. The study sites represent different ecosystems of arid cold environments, which include grasslands, savannah, and shrubland (see Appendix A for a general description of the sites).

ET a Fluxes and Meteorological Data
Actual evapotranspiration data in the Chilean sites were obtained from three eddy covariance systems (IRGASON, Campbell Sci., Utah, United States), each one having a meteorological station that allowed measuring net radiation (R n ) (CNR4, Kipp & Zonen, The Netherlands), soil heat flux (G) (HFP01SC, Hukseflux, The Netherlands), precipitation (PPT) (TE525, Campbell Sci., Logan, UT, USA), atmospheric pressure (P) (PTB110, Vaisala, Helsinki, Finland), air temperature (T) and relative humidity (RH) (CS215, Campbell Sci., Logan, UT, USA), soil temperature (Ts) (TCAV, Campbell Sci., Logan, UT, USA), and soil's volumetric water content (VWC) (CS655, Campbell Sci., Logan, UT, USA). Vapor pressure deficit (VPD) was estimated using the previous meteorological data, and the wind speed (WS) in these sites was calculated using the measurements of the eddy covariance sonic anemometer. On the other hand, the data from Australia and the United States were obtained from the FLUXNET 2015 dataset [5,[39][40][41][42][43][44][45]. The sampling frequency of the data was 30 min in all sites, which was then integrated into hourly, daily, and monthly timescales. The processing/correction methods of the eddy covariance data of all the sites are described in detail in [5,30,[46][47][48][49][50]. Briefly, the following data filtering methods were applied: (i) only high-quality infilled data were chosen; (ii) ET a data samples collected in rainy periods were removed; and (iii) daytime data were utilized to avoid stable boundary layer conditions. Moreover, ET o was estimated in each site with the Penman-Monteith equation [7]. Table 2 presents the period where the data were collected in each site as well as the height at which the ET a measurements were performed. timescales. The processing/correction methods of the eddy covariance data of all the sites are described in detail in [5,30,[46][47][48][49][50]. Briefly, the following data filtering methods were applied: (i) only high-quality infilled data were chosen; (ii) ETa data samples collected in rainy periods were removed; and (iii) daytime data were utilized to avoid stable boundary layer conditions. Moreover, ETo was estimated in each site with the Penman-Monteith equation [7]. Table 2 presents the period where the data were collected in each site as well as the height at which the ETa measurements were performed.  [51], (c) CH-AT2 [51], (d) CH-AT3 [51], (e) AU-Cpr [48], (f) AU-Ync [49], (g) US-Cop [40], (h) US-SRG [42,52], (i) US-SRM [43,53], (j) US-Whs [44,54] , and (k) US-Wkg [45]. The latitude, longitude, and elevation of each study site are presented in Table 1.  [51], (c) CH-AT2 [51], (d) CH-AT3 [51], (e) AU-Cpr [48], (f) AU-Ync [49], (g) US-Cop [40], (h) US-SRG [42,52], (i) US-SRM [43,53], (j) US-Whs [44,54] , and (k) US-Wkg [45]. The latitude, longitude, and elevation of each study site are presented in Table 1. To incorporate the remote sensing data, as described below, it is important to estimate the footprint of the ET a measurements. Here, we approximated the footprint to a circle whose radius corresponds to the position of the footprint peak, x max , following the Schuepp et al. [55] approach: where u is the average wind speed (m/s), u * is the average friction velocity (m/s), z is the measuring height (m), d is the displacement height (m), and k is the von Kármán constant [56]. The footprint in each site was calculated using the sampling frequency of the data, i.e., 30 min, and then aggregating the footprint at a monthly level [57]. This approximation was chosen instead of a more complex approach, such as the Kljun et al. [57] model, for one principal reason: the required information for more complex footprint models is not available in the FLUXNET dataset, e.g., the crosswind distance standard deviation (σ y ) [57]. The performance of the Schuepp et al. [55] approximation was assessed by comparing the area and VI values obtained with this model and with the 80% flux footprint calculated with the Kljun et al. [57] approach in the Chilean sites, where all the required information was available ( Figure 2). Note that this approach has recently been used to investigate landscape transformation processes in urban areas located in arid regions [58].
Water 2021, 13, x FOR PEER REVIEW 6 of 23 one principal reason: the required information for more complex footprint models is not available in the FLUXNET dataset, e.g., the crosswind distance standard deviation (σy) [57]. The performance of the Schuepp et al. [55] approximation was assessed by comparing the area and VI values obtained with this model and with the 80% flux footprint calculated with the Kljun et al. [57] approach in the Chilean sites, where all the required information was available ( Figure 2). Note that this approach has recently been used to investigate landscape transformation processes in urban areas located in arid regions [58].

Remote Sensing and Vegetation Indices
Reflectance images were obtained from the Level-2 science products of the Landsat 7 satellite mission and then analyzed via Google Earth Engine [59] to estimate different VIs to be incorporated into the ETa estimates. Every selected image corresponded to the less cloudy image of each month.
For every selected image, the following VIs were calculated: the NDVI, the Soil Adjusted Vegetation Index (SAVI), the Enhanced Vegetation Index (EVI), the NDWI, and the Normalized Difference Greenness Index (NDGI). Then, the average of each VI was obtained in the footprint area approximated by the Schuepp et al. [55] approach described above. These VIs were selected as they are easy to implement, and they correctly represent the vegetation state over time periods of weeks, months, and years [60]. The main disadvantages of these VIs are that (i) they have shown errors when used in bare soils; (ii) they do not have a physical meaning, which complicates a direct comparison between VI among different sites; and (iii) they are not good at representing stress effects on vegetation in the short term (hours to days) [27,60].

Remote Sensing and Vegetation Indices
Reflectance images were obtained from the Level-2 science products of the Landsat 7 satellite mission and then analyzed via Google Earth Engine [59] to estimate different VIs to be incorporated into the ET a estimates. Every selected image corresponded to the less cloudy image of each month.
For every selected image, the following VIs were calculated: the NDVI, the Soil Adjusted Vegetation Index (SAVI), the Enhanced Vegetation Index (EVI), the NDWI, and the Normalized Difference Greenness Index (NDGI). Then, the average of each VI was obtained in the footprint area approximated by the Schuepp et al. [55] approach described above. These VIs were selected as they are easy to implement, and they correctly represent the vegetation state over time periods of weeks, months, and years [60]. The main disadvantages of these VIs are that (i) they have shown errors when used in bare soils; (ii) they do not have a physical meaning, which complicates a direct comparison between VI among different sites; and (iii) they are not good at representing stress effects on vegetation in the short term (hours to days) [27,60].

Normalized Difference Vegetation Index (NDVI)
The NDVI is the most utilized VI because it is strongly correlated with several biophysical characteristics and physiological processes of plants, including ET a [22]. The NDVI ranges between −1 and 1, where negative values correspond to water pixels, positives values but near 0 correspond to bare soil, and values near 1 are related to dense canopy. The NDVI is calculated as [12,22].
where ρ NIR corresponds to the reflectance of the NIR band (0.77-0.90 µm) and ρ RED (0.63-0.69 µm) is the reflectance of the visible red band. As the flux footprint areas are located at latitudes lower than 40 • and have steepness lower than 5 • , topographic illumination bias in the NDVI is expected to be negligible [61].

Soil-Adjusted Vegetation Index (SAVI)
The SAVI is a VI derived from the NDVI that includes a correction factor L, which minimizes the variations produced by the soil presence in heterogeneous surfaces. This index is calculated as [22].
The optimal value of L decreases as vegetation cover increases, i.e., L = 1 when the density is low, L = 0.5 for intermediate vegetation cover, and L = 0.25 for high density. For this investigation, L = 0.5 was used, because this value has shown good performance in reducing the noise produced by the presence of bare soil in a great range of vegetation cover densities [62].

Enhance Vegetation Index (EVI)
The EVI was developed to improve the sensitivity of the signal in high-biomass regions and to reduce the atmosphere influence. The EVI responds better than the NDVI to structural changes in plants and extends the range over which the NDVI responds to increases in foliage density [12,20]. The EVI is calculated as where C1 and C2 are area correction coefficients used to account for aerosol resistance, using the blue band to correct the influence of the aerosol in the red band. ρ BLUE (0.45-0.52 µm) is the reflectance of the blue band, G f is the gain factor (set as 2.5 [22]), and L c is the canopy background adjustment (set to 1 [22]). C1 and C2 were set as 6 and −7.5, respectively [22].

Normalized Difference Water Index (NDWI)
The NDWI, unlike the others VIs, focuses on identifying trends in the humidity of the studied surface, combining the water content of bare soil and vegetation. The NDWI is defined as [63].
Ji et al. [23] suggested naming this index the Normalized Difference Infrared Index (NDII), because the NDWI was first used in MODIS, whose SWIR band is between 1.23 and 1.25 µm. However, in this study, we prefer to call it the NDWI.

Normalized Difference Greenness Index (NDGI)
The NDGI is a VI developed to minimize variations between background reflectance of different surfaces and to maximize the contrast between vegetation and other background components in order to prevent the effects of snow in phenology estimation [64]. The NDGI is calculated as follows [64]: where ρ GREEN is the reflectance of the green band (0.52-0.60 µm) and ε is a coefficient that depends on the satellite (ε = 0.63 for Landsat 7).

Determination of Main Variables and ET a Estimates Using Machine Learning
The general procedure to generate ET a estimates is shown in Figure 3. Remote sensing, meteorological, and flux data were used as an input in the exhaustive feature selection (EFS) algorithm [65] to determine the main variables that control ET a . The EFS algorithm selects the subset of the original features or variables that better achieves an objective, usually finding the high value of a performance metric given by an arbitrary linear regressor or classifier [65]. The EFS algorithm is the most computationally expensive feature selection method because it needs to evaluate all possible combinations of the original factors considering that only a certain amount of these should be selected [66]. However, the EFS is the optimal feature selection method as the size of the dataset and the number of required features allow this method to be computationally feasible [65]. In this research, a subset of four features was selected, and a maximum of 18 features that are typically used in ET studies were evaluated (Table 3). It was decided to use subsets of four factors because, in preliminary tests, no significant improvement in formulae performance was achieved when more factors were added. The NDGI is a VI developed to minimize variations between background reflectance of different surfaces and to maximize the contrast between vegetation and other background components in order to prevent the effects of snow in phenology estimation [64]. The NDGI is calculated as follows [64]: where ρGREEN is the reflectance of the green band (0.52-0.60 μm) and ε is a coefficient that depends on the satellite (ε = 0.63 for Landsat 7).

Determination of Main Variables and ETa Estimates Using Machine Learning
The general procedure to generate ETa estimates is shown in Figure 3. Remote sensing, meteorological, and flux data were used as an input in the exhaustive feature selection (EFS) algorithm [65] to determine the main variables that control ETa. The EFS algorithm selects the subset of the original features or variables that better achieves an objective, usually finding the high value of a performance metric given by an arbitrary linear regressor or classifier [65]. The EFS algorithm is the most computationally expensive feature selection method because it needs to evaluate all possible combinations of the original factors considering that only a certain amount of these should be selected [66]. However, the EFS is the optimal feature selection method as the size of the dataset and the number of required features allow this method to be computationally feasible [65]. In this research, a subset of four features was selected, and a maximum of 18 features that are typically used in ET studies were evaluated (Table 3). It was decided to use subsets of four factors because, in preliminary tests, no significant improvement in formulae performance was achieved when more factors were added.
where Var i is the main variables that explain ET a , which are found using the EFS algorithm; a i represents the regression coefficients of the linear equation. The regression coefficients were found with ordinary least squares (OLS) [67]. To identify the main variables and the regression coefficients, the input data were normalized; i.e., each one of the Var i ranged between 0 and 1. This normalization ensured that the EFS chose the main variables for their contribution to the ET a variability and not because of their magnitude. The coefficient of determination (R 2 ) of the linear regression was chosen as the performance metric. The combination of EFS and OLS was chosen to create the ET a estimation equations instead of non-linear machine learning models like tree-base models, because linear equations allow a deeper interpretation of the results. Data were separated into two groups: the training data (CH-AT2, CH-AT3, AU-Cpr, US-Cop, US-SRG, US-SRM, and US-Wkg) and the validation data (CH-AT1, AU-Ync, and US-Wkg). The training data were used to generate a global estimation that could fit all of the sites. The performance of this equation was evaluated with the validation data. Furthermore, site-specific formulae were constructed with the data of each site with the aim of determining the most important factors minimizing the influence of the amount of data that each of the sites has on the results. Global and site-specific equations were found for daily and monthly time scales, both with ET a expressed as mm/day, with only meteorological data. Then, VIs were incorporated in monthly estimations to evaluate the relevance of incorporating remote sensing data into estimations that consider places with different cover types but the same climate.

Remote Sensing Information
For all the Chilean sites, a high correlation was found (coefficient of determination, R 2 > 0.84 and root mean square error, RMSE < 0.17) between VI values calculated with both footprint approaches, despite the difference in the footprint areas (Table 4 and Figure 4). Moreover, when increasing x max in more than one order of magnitude, the high correlation between VI values calculated with the Schuepp et al. [55] and the Kljun et al. [57] approaches was maintained and only decreased when the footprint reached a surface that had different characteristics. Hence, although the Schuepp et al. [55] approach may not precisely represent the footprint, it allows the estimation of VIs that agree with those calculated with a more sophisticated footprint method, such as the Kljun et al. [57] approach.  [57] and the Schuepp et al. [55] approaches. The mean values of each VI calculated with both approaches in the study period are presented. In addition, the coefficient of determination (R 2 ) and root mean square error (RMSE) between the results of the Kljun et al. [57] and the Schuepp et al. [55] approaches are shown. The temporal evolution of ETa, precipitation, and VIs on the validation sites is shown in Figures 5-7. At the US-Wkg validation site, ETa responded to water availability determined by the amount of precipitation ( Figure 5). Additionally, the VIs responded in the same way as ETa, except for the SAVI and NDVI, whose values decreased drastically in the presence of precipitation events. This behavior is not common for all the sites. For example, at the AU-Ync validation site (Figure 6), the relationship between precipitation and ETa is weak. However, it seems that ETa responded to water availability, represented as the NDWI. Furthermore, the EVI explains some of the ETa peaks. In the CH-AT1 validation site, no relationship between the Vis and ETa was found (Figure 7), most likely because groundwater is shallow in this site. Hence, there is water available for evapotranspiration throughout the year.  The temporal evolution of ET a , precipitation, and VIs on the validation sites is shown in Figures 5-7. At the US-Wkg validation site, ET a responded to water availability determined by the amount of precipitation ( Figure 5). Additionally, the VIs responded in the same way as ET a , except for the SAVI and NDVI, whose values decreased drastically in the presence of precipitation events. This behavior is not common for all the sites. For example, at the AU-Ync validation site (Figure 6), the relationship between precipitation and ET a is weak. However, it seems that ET a responded to water availability, represented as the NDWI. Furthermore, the EVI explains some of the ET a peaks. In the CH-AT1 validation site, no relationship between the Vis and ET a was found (Figure 7), most likely because groundwater is shallow in this site. Hence, there is water available for evapotranspiration throughout the year.

ETa Estimation Formulae
The global and the site-specific ETa estimation formulae are presented in Table 5.  (Figure 8c). In all cases, monthly estimations were more accurate than daily estimations, especially because monthly averages can mask outliers. In general, the monthly global equation that only considers meteorological information performed better than the equation that includes VIs. However, in general, the site-specific equations that include VIs resulted in better outcomes (Table 5).

ET a Estimation Formulae
The global and the site-specific ET a estimation formulae are presented in Table 5 (Figure 8c). In all cases, monthly estimations were more accurate than daily estimations, especially because monthly averages can mask outliers. In general, the monthly global equation that only considers meteorological information performed better than the equation that includes VIs. However, in general, the site-specific equations that include VIs resulted in better outcomes (Table 5).

ETa Estimation Formulae
The global and the site-specific ETa estimation formulae are presented in Table 5.  (Figure 8c). In all cases, monthly estimations were more accurate than daily estimations, especially because monthly averages can mask outliers. In general, the monthly global equation that only considers meteorological information performed better than the equation that includes VIs. However, in general, the site-specific equations that include VIs resulted in better outcomes (Table 5).    As shown in Figure 9, in the validation sites, the best results were obtained at US-Wkg, whereas the results at CH-AT1 and AU-Ync were not satisfactory. As most of the training data came from a site located near US-Wkg, these results imply that accurate ET a estimations are obtained when the training and validation sites have similar characteristics. Figure 9 shows the R 2 and RMSE values for the three validation sites. Daily estimations were usually less accurate than monthly estimates, excluding that for the AU-Ync site. Moreover, estimations that included a VI performed better than those that only considered meteorological information. The case with the best performance corresponds to the monthly estimate that includes a VI in US-Wkg (R 2 of 0.82 and RMSE of 0.42 mm/day).
In the validation sites, only acceptable results were obtained for the US-Wkg cases, most likely because a large amount of the training data came from a site located near US-Wkg. Figure 9 shows the R 2 and RMSE values for the three validation sites. Daily estimations were usually less accurate than monthly estimates, excluding that for the AU-Ync site. Additionally, estimations that include a VI performed better than those that only considered meteorological information. The case with the best performance corresponds to the monthly estimate that includes a VI in US-Wkg (R 2 of 0.82 and RMSE of 0.42 mm/day). Water 2021, 13, x FOR PEER REVIEW 14 of 23 In the validation sites, only acceptable results were obtained for the US-Wkg cases, most likely because a large amount of the training data came from a site located near US-Wkg. Figure 9 shows the R² and RMSE values for the three validation sites. Daily estimations were usually less accurate than monthly estimates, excluding that for the AU-Ync site. Additionally, estimations that include a VI performed better than those that only considered meteorological information. The case with the best performance corresponds to the monthly estimate that includes a VI in US-Wkg (R² of 0.82 and RMSE of 0.42 mm/day).

Variables Controlling ETa
For both the daily and monthly global equations, the main variables that influence ETa, obtained as a result of applying the EFS algorithm, are Rn − G, ETo, Tmin, and Tsmax (Table 5). In both cases Rn − G is the variable that has the highest regression coefficients and, hence, is the variable that can represent better temporal evolution of ETa. Note that although ETo is among the variables that are typically used to determine ETa [10], and that it includes Rn − G, RH, VPD, and WS, it is not always among the best ranked main variables. This is most likely due to limited water availability in many of the selected study Figure 9. Performance of the global formulae in the validation sites. CH-AT1, AU-Ync, and US-Wkg are shown from left to right. Daily, monthly, and VI monthly formulae are shown from the top to the bottom. Each panel includes the main variable selected for the construction of each formula, the RMSE, and R 2 . The blue band corresponds to the 95% confidence interval.

Variables Controlling ET a
For both the daily and monthly global equations, the main variables that influence ET a , obtained as a result of applying the EFS algorithm, are R n − G, ET o , T min , and Ts max ( Table 5). In both cases R n − G is the variable that has the highest regression coefficients and, hence, is the variable that can represent better temporal evolution of ET a . Note that although ET o is among the variables that are typically used to determine ET a [10], and that it includes R n − G, RH, VPD, and WS, it is not always among the best ranked main variables. This is most likely due to limited water availability in many of the selected study sites; hence, other variables, such as volumetric water content, could be controlling ET a . In the case where remote sensing information is used, the main variables are R n − G, PPT, the NDGI, and the NDWI. The variable with the greatest contribution to this equation is the NDGI and that with the lowest was the R n − G. Table 6 shows the occurrence of the main variables found for all the sites for the site-specific equations. For the daily estimates, the most important variables are R n − G, VWC, and ET o , which show that daily ET a depends on both energy and water availability. In the case of the monthly estimations, the main variables are R n − G, T, and Ts. The monthly estimates that include VIs have VPD and the NDWI as the principal variables.
When only meteorological information is considered, the available energy is the principal variable that explains ET a in both daily and monthly scales. In the case when remote sensing information is incorporated the main variable is NDWI, however, the available energy is still equally important as other variables that represent the availability of water, such as PPT and VPD.

Discussion
A comparison of ET a estimation formulae between studies is difficult due to the many differences between them: (1) calibration and validation procedures; (2) data selection and processing; (3) temporal scale of estimates; (4) number and characteristic of the variables used; and (5) number and location of field sites considered [20]. However, in this section, the results obtained in other studies that have used regression formulae or machine learning algorithms to estimate ET a are discussed and compared to our results. Carter and Liang [4] evaluated seven regression algorithms for daily ET a estimations with meteorological and/or remote sensing data of different cover types, reaching R 2 between 0.43 to 0.52 for all sites. The performance obtained by Carter and Liang [4] in the sites with climate different than arid cold desert is slightly lower than that obtained in this study for daily estimations considering the training data (R 2 = 0.6). The algorithms evaluated by Carter and Liang [4] correspond to simple linear equations, such as the Yebra et al. [20] formula, and to more complex equations, such as that developed in [68] Granata [28] fitted three daily ET a estimations models that include different meteorological data with four different machine learning algorithms in a subtropical humid site located in Florida. All of them reached R 2 values of over 0.92, because only one site was considered and the availability of water is not limited, as opposed to what occur in arid cold regions. However, better results were obtained in the model with a greater number of variables.
Studies in natural arid zones landscapes are scarce compared to studies performed in agricultural lands located in mesic environments. Investigations performed in the western of the U.S. have provided the basis for improved estimation of ET a in arid and semi-arid environments [2,3,27,60,69]. Bunting et al. [3] evaluated three regression equations that estimates ET a in a period of 16 days in riparian and upland sites in California. One of the equations is a multiple linear regression that includes the MODIS EVI and precipitation data (R 2 = 0.74). Nagler et al. [2,27] developed two different regression equations that require meteorological and MODIS EVI information to estimate ET a in riparian environments of the Colorado, Rio Grande, and San Pedro rivers in Colorado, U.S. Both equations are based on the relationship between the leaf area index (LAI) and light absorption by the canopy, and the linear relationship between the EVI and LAI. Both equations have good predictive capability (R 2 = 0.73 and 0.74, respectively). Performances of the best results obtained in this research are comparable to the studies reviewed.
In our work, different arid cold climate sites were used to generate linear regression formulae to estimate ET a . Although other works show better results in sites that only have the same vegetation cover [20], we obtained very different performances in daily, monthly, and monthly with VIs estimations in the validation sites with similar vegetation cover, such as in the case of AU-Ync and US-Wkg (R 2 = 0.03, 0.00, 0.16 and R 2 = 0.69, 0.78, 0.82, respectively). In regard to this, note that the linear formulation of the regression formulae should not be an important source of error, as Carter and Liang [4] demonstrated that different regression formulae, with different theoretical bases and the same input data, have similar performance.
According to Allen et al. [7], the main meteorological variables affecting ET a are radiation, air temperature, air humidity, and wind speed. Several investigations have studied the relative importance of these variables in ET a processes in arid regions. However, these studies generally focus on the behavior of ET o instead of ET a , so they do not consider the effects of water stress. For example, Adnan et al. [70] and Eslamian et al. [71] studied the influence of meteorological variables on ET o estimations in semi-arid, arid, and hyper-arid climates (Pakistan and Iran) using the Penman-Monteith formulation. Both studies concluded that air temperature and humidity are the most important meteorological variables affecting ET o . One of the few studies that analyze the sensitivity of ET a estimations to variations in meteorological and remote sensing data in a semi-arid region is that conducted by Mokhtari et al. [72]. They analyzed the sensitivity of METRIC (mapping evapotranspiration at high resolution with internalized calibration) and concluded that it is highly sensitive to surface temperature, net radiation, and air temperature, and it is less sensitive to the LAI, SAVI, and WS (except for WS at low level of vegetation cover).
Our results indicate that available energy is the main variable that controls ET a , which agrees with previous studies that investigated ET a components in as many climate and vegetation cover types as possible [4,[73][74][75]. For example, Wang et al. [73] correlated ET a measurements with radiation, air, and land surface temperature, the EVI and NDVI, and soil moisture. They concluded that correlation coefficients between R n and ET a are the highest, followed by Ts and VIs.
In this research, of the four most important meteorological variables, only wind speed was not decisive in estimating ET a in all of the cases studied. This fact agrees with the findings of Granata [28], who proved that it is possible to generate accurate and precise estimates of daily ET a through machine learning algorithms only with mean temperature, net solar radiation, and relative humidity data, pointing out that the incorporation of wind speed does not improve the ET a estimations compared to the case when it is not accounted for. However, he analyzed ET a in a subtropical humid climate, where the number of sunshine hours is considered to be the more dominant variable as opposed to arid climates, where wind speed is an important variable [46]. Irmak et al. [76] compared 11 ET a models in a crop field in Nebraska, USA, to study their complexity on hourly, daily, and seasonal scales. They concluded that wind speed, and other meteorological variables such as temperature, gained importance in daily and hourly calculations, while on seasonal scales, radiation was the dominant variable. As shown in these studies, it was expected that wind speed would be an important variable in daily ET a estimations. However, the method and the number of variables chosen in this research could mask its effects: EFS selects the most important meteorological variable or variables that explain ET a , in this case R n − G and the NDWI, accompanied by variables whose unique objective is to make the equation work numerically; furthermore, the WS influence could be well represented in ET o , so R n − G, VWC, and Ts are variables that provide more information about ET a variability than WS itself. In arid regions, WS plays an important role when advection of dry air enhances evaporation and affects the energy balance by horizontal transport of latent heat [46,77].
Our findings support that regardless of the climate type, atmospheric demand and available energy determine ET a when water supply is sufficient, whereas soil moisture becomes an important factor controlling ET a when soil water supply is deficient [75]. Bunting et al. [3] proved that ET a estimations in semi-arid upland sites using multiple linear regression improve with the incorporation of a moisture input. However, variables such as precipitation and soil moisture are not usually used as the moisture input for several reasons: (1) surface precipitation and soil moisture measurements are point measurements, limiting the possibilities for upscaling; (2) a lag effect must be considered with precipitation; and (3) soil moisture remote sensing products are difficult to process, and their resolution is of several kilometers [4]. In this context, we expected that few of the ET a estimation formulae presented in Table 5 would incorporate PPT and VWC in the four most important variables found by the EFS algorithm. However, a novel and non-expected result that we found when remote sensing information was added, is the appearance of the NDWI as the main variable that explains the moisture input (see ET a estimation formulae shown in Table 5). Moreover, in studies performed in climates different than arid cold deserts, NDVI and EVI are the remote sensing VIs that are typically used to determine ETa because they represent the vegetation activity [4,26,75,78]. Unlike other VIs, NDWI is capable of indicating trends in soil and vegetation wetness [25,79], so it is a valid water availability input that does not have the same disadvantages of PPT and VWC, as mentioned above. Hence, NDWI is a VI that improves the estimation of ET a in arid cold regions, and that has not received too much attention in studies conducted in other climates.
The use of remote sensing information is fundamental in estimating ET a for regionalscale and heterogeneous landscapes [12]. This research proved that the incorporation of VIs helps to extrapolate global equations to each one of the sites. However, it has been proven that VIs are not sufficient to accurately estimate ET a [20]. Carter and Liang [4] noted that, at minimum, ET a estimates with Vis require the inclusion of radiation data. However, it is recommended to increase the number of input variables. Our results demonstrate that acceptable results were achieved with four variables.
Although the contribution of the VIs to the improvement in ET a estimations at the regional level is indisputable, there are several sources of error that must be addressed. One of the most important is the influence of bare soil on the reflectance response, especially in high-resolution satellites, such as Landsat. Jarchow et al. [69,80] compared Landsat 5 and Landsat 8 EVI values to the MODIS EVI in a riparian zone of the Colorado River, Mexico, finding low correlations over bare soil and sparsely vegetated areas. Additionally, they suggest being cautious when high-resolution Landsat EVI data are analyzed over heterogeneous areas with low vegetation densities, such as those commonly encountered in cold arid and semi-arid environments, because soil presence contributes to increased variability in the response of the NIR and red bands.
The poor correlations obtained in this study between VIs and ET a could be explained by several factors. Firstly, as mentioned before, the presence of bare soil can perturb the calculation of VIs [80]. Secondly, in this research, only ET a outliers were extracted, whereas other studies selected data that accomplished some characteristics. For example, Yebra et al. [20] selected data of days where only transpiration was expected to be dominant, and Scott et al. [81] excluded data from precipitation events and outliers of meteorological variables. In the presence of important rainfall events, most of the VIs considered in this study, except for the NDWI and EVI, have negative values. The values of the VIs indicate that there should be a lower ET a rate when it rains, since they actually increase. VI values obtained in this research are different to those reported in previous studies [82][83][84]. However, they are different from each other, highlighting the importance of satellite selection.

Conclusions
In this study, we generated linear regression formulae to estimate daily and monthly ET a in arid cold sites. Different performances were obtained for every site, and the following trends were identified: (1) better results were obtained for monthly than for daily estimates; (2) incorporation of remote sensing information allows one to extrapolate formulae to other sites in order to obtain better results than estimations with only meteorological data; (3) the available energy is the most important meteorological variable in ET a estimations for the sites evaluated in this research; and (4) in arid regions, it is important to incorporate estimations of water availability. As precipitation and soil moisture are point measurements that do not allow one to extrapolate estimations in wide areas, the NDWI could be incorporated as a proxy for water availability in the heterogeneous landscapes located in arid cold regions. Furthermore, more studies that analyze variables controlling ET a in arid natural landscapes are needed, because ET a in drylands is exposed to different factors than in more humid environments, such as water stress, advection, and vegetation with adaptations to drought. Global ET o investigations cannot study the complexity of ET a in arid cold regions in depth.