## 1 Introduction  
The global water cycle is changing due to the combined effects of climate change and human interventions during the 21 st century [1]. One of the greatest challenges is keeping water consumption at sustainable levels, which is more complex due to the increasing population in a context of climate uncertainty [2; 3] and 3.5-4.4 billion people estimated under water scarcity conditions in 2050 [4]. Many regions of the world can expect a combination of increasing temperatures (largely increasing evaporative demand) and decreasing precipitation patterns, which leads to increased stress on tackling water demand [5]. A prime example of this is the Mediterranean region, which is consistently projected as a \"hotspot\" of drying trends and prolonged water scarcity conditions [6; 7]. The Iberian Peninsula is predicted to be among the most affected areas by severe droughts by the end of the 21 st century [8]. In this region, where irrigated agriculture represents over 80% of the total extracted water [9], land use shifts towards higher market-valued crops represent a major driver of change, which will markedly increase water withdrawals [10]. In the Guadalquivir basin in Spain, irrigation water requirements are expected to increase between 15% and 20% by 2050 [11]. This may cause a redistribution of water between the surface and groundwater [12]. Monitoring the variations in the hydrological cycle components more closely, among them evapotranspiration (_ET_), is of major importance in countries facing intensified drought spills [13; 14]. After precipitation, _ET_ is the major component in the terrestrial water budget [15]. Evapotranspiration at specific locations can be estimated by using methods such as lysimeters, the Bowen-ratio, or Eddy Covariance (EC) [16]. However, _ET_ estimations at large-scales are complex to achieve due to spatial heterogeneity in the land surface [17; 18; 19]. As the only connecting term between the water cycle and the land surface energy budget [20; 21; 22], _ET_ can be estimated by hydrological models (HMs) relying on the water balance or by land surface models (LSMs) computed through the mass-transfer equations between the land surface and the atmosphere [23]. Conceptual lumped HMs consider an undivided entity with individual values of input variables and parameters [24]. Spatially distributed HMs or prognostic LSMs enable _ET_ to be estimated at large scales in a spatially explicit manner [25]. This type of models requires a priori information on land use and surface properties, which are not always up-to-date and may be complex for retrieving in situations of fast land use changes such as new irrigation developments. Diagnostic or historical LSMs using remote sensing (RS) datasets to prescribe vegetation and other surface variables can provide spatially distributed fluxes and fast updating on the land surface [26; 27]. Diagnostic RS-LSMs need limited a priori data on soil and vegetation parameters to estimate _ET_ and can rely on remotely sensed land surface temperature (_LST_) as a proxy to the surface moisture status [28]. The complexity of RS-based models to retrieve _ET_ depends on the process considered [29; 30]. _ET_ can be directly estimated as a residual of the surface energy balance equation or can be based on potential evapotranspiration adjusted to actual rates using different constraints factors from remote sensing. The Penman-Monteith (PM) equation [31; 32] is the base to compute _ET_ in some global RS-_ET_ products such as the first widely available _ET_ product MOD16 [33; 34] whose main challenge relies on estimating surface and aerodynamic resistances to the water vapor [35]. [PERSON] et al. [36] found that MOD16 best performs in temperate and fully humid climates and produces a consistent underestimation of _ET_ in semiarid regions such as Spain, Portugal, the neighborhood of the Black Sea and the Caspian Sea, and semiarid climatic zones in the conterminous United States [37].