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 ETa 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 ETa 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 ETa 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.
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
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
Rn 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
Rn −
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
Rn −
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.