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

Abstract

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 find 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² = 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.

1. 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 vegetational characteristics that influence it and make its estimations more complex [7,8]. Major challenges in ET_a estimation are those that make the process even more dynamic over time and variable in space. For example, when water is a limiting factor, plants decrease their transpiration through physiological adaptations such as stomata closure. Furthermore, regional advection can bring additional energy to the system, resulting in higher ET_a rates that can even exceed potential ET [8]. As arid lands are vastly different from irrigated farms, typically having native vegetation with high resistance to transpiration and low ground cover, traditional ET_a estimation methods are not suitable for these environments. The most common approach to determine ET_a is the crop coefficient method, in which ET_a is estimated from the reference evapotranspiration (ET_o), computed from the Penman–Monteith formula [7] and using a crop coefficient (K_c), i.e., ET_a = K_c × ET_o, where K_c represents four effects that distinguish the crop from reference grass: aerodynamic resistance, albedo, surface resistance, and soil evaporation [7,8]. Nonetheless, it has been demonstrated that the crop coefficient method in general is not suitable for determining ET_a of native vegetation adapted to arid conditions, because transpiration is overestimated when plants encounter suboptimal conditions of soil water as a result of not considering stomatal regulation and plant adaptations to drought [9]. Only in few exceptions has the crop coefficient method been suitable to estimate ET_a in crops cultivated in arid or semi-arid regions [10,11].

Remote sensing methods have been developed to estimate ET_a and have been positioned as the only feasible approach for wide areas of mixed landscapes, allowing for an improvement in water balance estimations over basin and regional scales [12,13,14,15]. One of the most used operational global ET_a satellite products is MOD16, which is based on Moderate Resolution Imaging Spectroradiometer (MODIS) information. Moreover, the Satellite Application Facility on Land Surface Analysis (LSA-SAF) is widely used, but it only covers Europe, Africa, and most of South America [16,17]. Hu et al. [16] demonstrated that LSA-SAF has a better performance than MOD16, but neither products capture ET_a in water-limited regions. ET_a data can also be obtained from the Agricultural Research Service of the United States Department of Agriculture (USDA-ARS) ET dataset, and the data provided from the European Centre for Medium-Range Weather Forecasts (ECMWF) or the Global Land Data Assimilation System (GLDAS) models [18,19]. However, their spatial resolutions are even coarser than that of the MOD16 and LSA-SAF datasets [17]. In 2019, the European Space Agency (ESA) released the Sen-ET open-source software application for ET_a modeling at high (tens of meters) and medium (1 km) spatial resolutions, and at a temporal resolution of ~5 to 10 days. The Sen-ET software uses observations of Sentinel-2 and Sentinel-3 for field-scale applications. The first validation procedure in latent heat flux was performed in the Skjern river basin (Denmark) and resulted in a correlation of 0.76 when comparing data obtained from three eddy covariance flux towers, with the best performance estimated in croplands [17]. The most common remote sensing ET_a approaches are based on the surface energy balance (SEB) equation, where sensible heat (H) is estimated using land surface temperature (LST) derived from thermal infrared (TIR) sensor on satellites [12]. Although these methods have been cataloged as operational, there are difficulties on their implementation: small errors in the estimation of the LST translate into large errors in H estimates, and only few sensors offer open source TIR data [12,20].

Vegetation index (VI)-based methods to estimate ET_a were developed to take advantage of remote sensing, avoiding the disadvantages associated with the methods based on SEB. VIs were developed for vegetation monitoring due to spectral reflectance signature revealing information about the state, biogeochemical composition, and structure of a leaf and canopy, but VIs can also provide information about water and carbon cycles [21]. ET_a estimation methods based on VIs depend on an estimate of the density of green vegetation over the landscape, as measured by VIs or related products that combine the VIS and IR bands [12]. For example, the Normalized Difference Vegetation Index (NDVI) captures the contrast in light reflection from green leaves between the red and near infrared (NIR) bands, because red light is strongly absorbed by chlorophyll and nearly all the NIR is 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 physics-constrained 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.

2. Materials and Methods

2.1. 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. Site name, location (latitude and longitude), elevation, International Geosphere-Biosphere Programme (IGBP) land cover classification, and mean annual precipitation and temperature of the study sites. The country of each site is coded as follows: Chile (CH), Australia (AU), and the United States (US).

Site	Location		Elevation (m ASL)	IGBP Land Cover Classification
	Latitude (°)	Longitude (°)			Mean Annual Precipitation (mm)	Mean Annual Temperature (°C)
CH-AT1	22.02 S	68.05 W	4182	Barren or sparsely vegetated	78	5.8
CH-AT2	22.01 S	68.05 W	4330	Barren or sparsely vegetated	78	5.8
CH-AT3	22.52 S	68.02 W	4255	Barren or sparsely vegetated	106	1.7
AU-Cpr	34.00 S	140.59 E	166	Open shrubland	250	18
AU-Ync	34.99 S	146.29 E	125	Grassland	419	17.3
US-Cop	38.09 N	109.39 W	1520	Grassland	216	12
US-SRG	31.79 N	110.83 W	1290	Grassland	377	19
US-SRM	31.82 N	110.87 W	1116	Open shrubland	377	19
US-Whs	31.74 N	110.05 W	1370	Grassland	285	17.6
US-Wkg	31.74 N	109.94 W	1530	Grassland	294	17.3

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

2.2. 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. Period of measurements, sensor height, dominant vegetation type, and approximate vegetation height of each site [39,40,41,42,43,44,45,51]. The country of each site is coded as follows: Chile (CH), Australia (AU), and the United States (US).

Site	Time Period		Sensor Height (m)	Dominant Vegetation Type	Vegetation Height (m)
	Start Date	End Date
CH-AT1	18-01-2018	29-05-2019	2.11	Oxychloe andina and grass Deyeuxia sp.	0.7
CH-AT2	22-02-2018	25-04-2019	2.11	Festuca genera	0.2
CH-AT3	19-04-2018	28-05-2019	2.49	Oxychloe andina, Festuca and Deyeuxia genera grass	0.2
AU-Cpr	01-01-2010	31-12-2014	20	Several species of Eucalyptus	4.0
AU-Ync	01-01-2012	31-12-2014	8	perennial tussock grasses	0.3
US-Cop	01-01-2001	31-12-2007	1.85	Hilaria jamesii, Stipa hymenoides bunchgrasses and Coleogyne ramosissima shrub	0.3
US-SRG	01-01-2008	31-12-2014	14	South African warm-season bunchgrass, Eragrostis lehmanniana, and Prosopis velutina	0.3
US-SRM	01-01-2004	31-12-2014	8	Prosopis velutina and native and nonnative perennial grasses, subshrubs, and scattered succulents	1.5
US-Whs	01-01-2007	31-12-2014	5	Parthenium incanum, Acacia constricta, Larrea tridentata, and Flourensia cernua	4.3
US-Wkg	01-01-2004	31-12-2014	5	Eragrostis lehmanniana, Bouteloua eripoda, and Aristida spp.	0.3

presents the period where the data were collected in each site as well as the height at which the ET_a measurements were performed.

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:

$$x_{max}=\frac{u}{u_{*}}\frac{\left(z-d\right)}{2\kappa }$$

(1)

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

2.3. 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].

2.3.1. 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].

$$NDVI=\frac{{\rho }_{\mathrm{NIR}}-{\rho }_{\mathrm{RED}}}{{\rho }_{\mathrm{NIR}}+{\rho }_{\mathrm{RED}}}$$

(2)

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

2.3.2. 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].

$$SAVI=\frac{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}-{\rho }_{\mathrm{RED}}}{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}+{\rho }_{\mathrm{RED}}+L}\left(1+L\right)$$

(3)

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

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

$$EVI=G_f\frac{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}-{\rho }_{\mathrm{RED}}}{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}+C1{\rho }_{\mathrm{RED}}+C2{\rho }_{\mathrm{BLUE}}+L_c}$$

(4)

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

2.3.4. 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].

$$NDWI=\frac{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}-{\rho }_{\mathrm{SWIR}}}{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}+{\rho }_{\mathrm{SWIR}}}$$

(5)

where ρ_SWIR is the reflectance of the short-wave infrared (SWIR) band (1.55–1.75 µm). 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.

2.3.5. 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]:

$$NDGI=\frac{\epsilon {\rho }_{\mathrm{GREEN}}+\left(1-\epsilon \right){\rho }_{\mathrm{NIR}}-{\rho }_{\mathrm{RED}}}{\epsilon {\rho }_{\mathrm{GREEN}}+\left(1-\epsilon \right){\rho }_{\mathrm{NIR}}+{\rho }_{\mathrm{RED}}}$$

(6)

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

2.4. 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. Variables analyzed in this study and their respective symbols and units.

Variable	Symbol	Units	Variable	Symbol	Units
Available energy	Rn − G	MJ m−2 d−1	Volumetric water content	WVC	cm3 cm−3
Precipitation	PPT	mm	Water vapor deficit	VPD	kPa
Mean temperature	T	°C	Wind speed	WS	m s−1
Minimum temperature	Tmin	°C	Reference evapotranspiration	ETo	mm
Maximum temperature	Tmax	°C	Normalized difference vegetation index	NDVI	-
Soil temperature	Ts	°C	Soil-adjusted vegetation index	SAVI	-
Minimum soil temperature	Tsmin	°C	Enhanced vegetation index	EVI	-
Maximum soil temperature	Tsmax	°C	Normalized difference water index	NDWI
Relative humidity	RH	-	Normalized difference greenness index	NDGI	-

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

A linear equation to estimate ET_a was constructed using the four main variables found with the EFS algorithm:

$${\mathrm{ET}}_{\mathrm{a}}=a_0+{\displaystyle \sum }_{i=1}^4{a}_{i}Va{r}_i$$

(7)

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²) 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.

3. Results

3.1. Remote Sensing Information

For all the Chilean sites, a high correlation was found (coefficient of determination, R² > 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. Comparison between vegetation index (VI) values obtained with the footprints calculated with the Kljun et al. [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²) and root mean square error (RMSE) between the results of the Kljun et al. [57] and the Schuepp et al. [55] approaches are shown.

Site	Footprint Model, R2 and RMSE	NDVI	SAVI	EVI	NDWI	NDGI
CH-AT1	Kljun mean	0.09	0.14	0.20	−0.03	−0.04
	Schuepp mean	0.08	0.12	0.19	−0.03	−0.04
	R2	0.91	0.92	0.94	0.88	0.84
	RMSE	0.01	0.02	0.03	0.03	0.01
CH-AT2	Kljun mean	0.04	0.07	0.08	0.00	−0.06
	Schuepp mean	0.05	0.07	0.09	0.01	−0.06
	R2	0.96	0.96	0.95	0.99	0.93
	RMSE	0.00	0.01	0.01	0.01	0.00
CH-AT3	Kljun mean	0.15	0.22	0.27	0.16	0.00
	Schuepp mean	0.22	0.33	0.41	0.25	0.05
	R2	0.99	0.99	0.99	0.87	0.99
	RMSE	0.09	0.14	0.17	0.10	0.06

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.

The temporal evolution of ET_a, precipitation, and VIs on the validation sites is shown in Figure 5, Figure 6 and Figure 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.

3.2. ET_a Estimation Formulae

The global and the site-specific ET_a estimation formulae are presented in

Table 5. Global and site-specific ET_a estimate formulae obtained with the training data set, and their respective coefficient of determination. ET_a is expressed in mm/day in all cases.

Site	ETa Estimate Formula	R2	Notes
Global	$E{T}_a=-2.86\times E{T}_o+10.25\times \left(R_n-G\right)+5.17\times T_{min}-4.18\times T{s}_{max}-3.31$	0.60	Daily
	$E{T}_a=4.26\times \left(R_n-G\right)-3.46\times E{T}_o+3.21\times T_{min}-2.26\times T{s}_{max}-0.07$	0.70	Monthly
	$E{T}_a=1.17\times \left(R_n-G\right)+2.36\times PPT+2.41\times NDGI+1.95\times NDWI-2.14$	0.67	Monthly with VI
CH-AT1	$E{T}_a=-4.19\times \left(R_n-G\right)-1.18\times T+2.05\times Ts+5.72\times E{T}_o+1.39$	0.90	Daily
	$E{T}_a=0.81\times E{T}_o-0.36\times T_{max}+0.89\times Ts+0.15\times T{s}_{max}+1.14$	0.99	Monthly
	$E{T}_a=0.41\times PPT+1.43\times Ts-0.62\times RH+0.02\times NDWI+1.28$	0.97	Monthly with VI
CH-AT2	$E{T}_a=-0.43\times T{s}_{max}+0.45\times RH+0.21\times VWC+0.49\times E{T}_o+0.08$	0.25	Daily
	$E{T}_a=0.13\times PPT+0.6\times RH-0.03\times VPD+0.09\times WS+0.07$	0.83	Monthly
	$E{T}_a=0.12\times PPT+0.62\times RH+0.09\times WS-0.04\times NDWI+0.06$	0.83	Monthly with VI
CH-AT3	$E{T}_a=-1.19\times T+3.81\times Ts-0.88\times VWC+1.73\times E{T}_o+0.73$	0.82	Daily
	$E{T}_a=-2.02\times PPT-1.24\times T_{min}+1.97\times Ts+2.83\times T{s}_{min}+0.77$	0.98	Monthly
	$E{T}_a=-1.75\times PPT+4.49\times Ts-1.4\times T{s}_{max}+0.85\times NDWI+0.59$	0.99	Monthly with VI
AU-Cpr	$E{T}_a=3.09\times \left(R_n-G\right)+3.49\times T-2.83\times T{s}_{max}-1.86\times E{T}_o-0.34$	0.58	Daily
	$E{T}_a=1.31\times \left(R_n-G\right)+3.15\times T-3.08\times Ts-0.67\times WS+0.31$	0.45	Monthly
	$E{T}_a=3.54\times \left(R_n-G\right)-2.76\times E{T}_o+0.25\times VPD-0.01\times NDGI-0.1$	0.34	Monthly with VI
AU-Ync	$E{T}_a=2.0\times \left(R_n-G\right)-2.45\times Ts+1.37\times VPD+0.56\times WS+0.32$	0.30	Daily
	$E{T}_a=1.62\times \left(R_n-G\right)-3.14\times Ts+1.4\times VPD+0.26\times WS+0.7$	0.55	Monthly
	$E{T}_a=1.55\times \left(R_n-G\right)-2.46\times Ts+1.26\times VPD+0.32\times NDWI+0.43$	0.67	Monthly with VI
US-Cop	$E{T}_a=1.81\times \left(R_n-G\right)+2.71\times T-1.15\times T_{max}-2.22\times VPD-0.43$	0.33	Daily
	$E{T}_a=1.23\times \left(R_n-G\right)+0.64\times PPT+0.22\times T_{min}-0.59\times VPD-0.03$	0.55	Monthly
	$E{T}_a=0.37\times PPT+0.97\times T-0.6\times VPD-0.06\times NDWI+0.33$	0.21	Monthly with VI
US-SRG	$E{T}_a=6.5\times \left(R_n-G\right)+5.56\times T_{min}-4.6\times T{s}_{max}+1.69\times VWC-3.39$	0.78	Daily
	$E{T}_a=2.97\times \left(R_n-G\right)+3.59\times T-4.5\times Ts+1.89\times RH-0.6$	0.89	Monthly
	$E{T}_a=4.28\times \left(R_n-G\right)-4.49\times E{T}_o+1.57\times T_{min}+0.87\times NDWI+0.25$	0.85	Monthly with VI
US-SRM	$E{T}_a=1.98\times \left(R_n-G\right)+5.01\times T-3.16\times VPD+2.44\times VWC-2.76$	0.76	Daily
	$E{T}_a=3.24\times \left(R_n-G\right)-4.86\times E{T}_o+2.76\times T+0.89\times VWC-0.12$	0.89	Monthly
	$E{T}_a=4.19\times T-3.05\times VPD+1.54\times VWC+0.43\times NDWI-0.63$	0.90	Monthly with VI
US-Whs	$E{T}_a=2.08\times \left(R_n-G\right)+1.27\times T{s}_{min}+3.12\times VWC-1.61\times E{T}_o-0.65$	0.77	Daily
	$E{T}_a=2.62\times \left(R_n-G\right)-3.37\times E{T}_o+1.85\times T+0.91\times VWC-0.05$	0.91	Monthly
	$E{T}_a=2.61\times T-2.5\times VPD+1.23\times VWC+0.44\times NDWI-0.27$	0.90	Monthly with VI
US-Wkg	$E{T}_a=5.86\times \left(R_n-G\right)+2.43\times T_{min}+1.92\times VWC-3.17\times E{T}_o-1.91$	0.71	Daily
	$E{T}_a=4.09\times \left(R_n-G\right)-3.78\times E{T}_o+5.36\times T-4.35\times Ts-0.2$	0.90	Monthly
	$E{T}_a=3.34\times \left(R_n-G\right)-4.36\times E{T}_o+1.96\times T+1.37\times NDGI-0.78$	0.89	Monthly with VI

. The daily global equation developed with the seven sites of the training data reached an R² of 0.60 and an RMSE of 0.64 mm/day (Figure 8a). On monthly timescales and only considering meteorological information, the global equation developed with the training sites reaches an R² of 0.70 and an RMSE of 0.47 mm/day (Figure 8b). The monthly global equation that does include VIs reached an R² of 0.67 and an RMSE of 0.49 mm/day (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² 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² 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² 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).

3.3. 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. Number of times that every variable is selected in a site-specific equation for daily, monthly, and monthly with VI estimates.

Variable	Daily	Monthly	Monthly with VI	Variable	Daily	Monthly	Monthly with VI
Rn − G	8	7	4	Tsmax	3	1	1
VPD	3	3	5	PPT	0	3	4
VWC	6	2	2	ETo	6	4	3
RH	1	2	2	WS	1	3	1
T	5	5	4	NDVI	-	-	0
Tmin	2	2	1	NDWI	-	-	8
Tmax	1	1	0	SAVI	-	-	0
Ts	3	6	3	EVI	-	-	0
Tsmin	1	1	0	NDGI	-	-	2

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.

4. 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² 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² = 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² 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² = 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² = 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² = 0.03, 0.00, 0.16 and R² = 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 regional-scale 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.

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