remote sensing Medium-Resolution Mapping of Evapotranspiration at the Catchment Scale Based on Thermal Infrared MODIS Data and ERA-Interim Reanalysis over North Africa

: Accurate quantiﬁcation of evapotranspiration (ET) at the watershed scale remains an important research challenge for managing water resources in arid and semiarid areas. In this study, daily latent heat ﬂux (LE) maps at the kilometer scale were derived from the two-source energy budget (TSEB) model fed by the MODIS leaf area index (LAI), land surface temperature (LST) products, and meteorological data from ERA-Interim reanalysis from 2001 to 2015 on the Tensift catchment (center of Morocco). As a preliminary step, both ERA-Interim and predicted LE at the time of the satellite overpass are evaluated in comparison to a large database of in situ meteorological measurements and eddy covariance (EC) observations, respectively. ERA-Interim compared reasonably well to in situ measurements, but a positive bias on air temperature was highlighted because meteorological stations used for the evaluation were mainly installed on irrigated ﬁelds while the grid point of ERA-Interim is representative of larger areas including bare (and hot) soil. Likewise, the predicted LE was in good agreement with the EC measurements gathered on the main crops of the region during 15 agricultural seasons with a correlation coefﬁcient r = 0.70 and a reasonable bias of 30 W/m 2 . After extrapolating the instantaneous LE estimates to ET daily values, monthly ET was then assessed in comparison to monthly irrigation water amounts provided by the local agricultural ofﬁce added to CRU precipitation dataset with a reasonable agreement; the relative error was more than 89% but the correlation coefﬁcient r reached 0.80. Seasonal and interannual evapotranspiration was analyzed in relation to local climate and land use. Lastly, the potential use for improving the early prediction of grain yield, as well as detecting newly irrigated areas for arboriculture, is also discussed. The proposed method provides a relatively simple way for obtaining spatially distributed daily estimates of ET at the watershed scale, especially for not ungauged catchments.


Introduction
Sustainable water resources management is a key issue in regions such as North Africa suffering historically from water scarcity and, in recent years, from an increase in drought frequency [1] and drastic changes in agricultural practices, leading to an increase in irrigation water demand [2][3][4][5][6]. Evapotranspiration (ET), which includes soil evaporation and plant transpiration, is a key component for the understanding of global water cycle [7], for catchment hydrology and streamflow prediction [8], as well as for agronomical studies, as it represents the crop water use [9,10]. However, vegetation ET and its partition is largely unknown at a global scale [11], particularly for highly anthropized catchments composed of heterogeneous agricultural landscapes mixing a large range of crops with complex structure [12] with various development stage and hydric conditions [13]. There is, thus, an imperious need for the development of methods to map ET at the catchment scale that are easy to implement in order to help water managers in their day-to-day decision-making processes.
In situ measurements are usually carried out using eddy covariance stations [14] or scintillometers [15]. While such instruments provide accurate measurements of ET at an infra-daily timestep, they require very specific skills, they are costly in human force and funding, and their measurements have a limited spatial representativity. By contrast, significant progress has been made in the last three decades on the mapping of ET based on remote sensing observations and modeling [16]. While ET as a water flux (or its energy-equivalent, the latent heat flux) cannot be directly observed from space, several approaches combining surface parameters and variables derived from remote sensing data with land surface models have been developed following basically two directions [17][18][19]: (1) the first is hereafter called the "water budget" approach, aimed at predicting the vertical profile of water availability for the plant within the root zone by solving a soil water budget. Several methods combine remote sensing products with land surface models of various complexity from concepts such as the FAO-56 double coefficient [3,20] to complex physically based soil-vegetation-atmosphere transfer (SVAT) models [21,22]. Nevertheless, solving a water budget model on a complex agricultural and irrigated landscape relies on a good knowledge of the irrigation water inputs, information that is largely unknown for large areas at the field scale [23][24][25][26]; (2) the second approach hereafter called "energy budget" is based on surface energy budget models fed by remotely sensed land surface temperature (LST) observations taken as a proxy of the crop hydric conditions when water availability is limited [27]. Indeed, when the plant undergoes stress, the stomata start closing and the energy dissipation through the transpiration processes is limited, leading to an increase in the crop temperature. In this context, several approaches based on LST have been developed to retrieve ET [17,[28][29][30][31]: (1) contextual approaches mapping evapotranspiration by projecting each pixel of an image in a two-dimensional space, consisting of surface temperature and the NDVI or the albedo [27,32]; (2) residual approaches which estimate the evapotranspiration as the residual of the surface energy balance [33][34][35]. These latter approaches showed good performance for reliable estimation of evapotranspiration over semiarid areas [36,37] and crop water stress [18,38,39]. Among the numerous residual models (SEBAL [40]; TSEB [35]; ALEXI [41]; DTD [42]; SPARSE [38]; METRIC [43]; ETEML [44]), the two-source energy balance (TSEB) model [35] has proven to be fairly robust for a wide range of landscapes [36,45,46] and weather conditions [9,27,28] and has already been calibrated, evaluated, and validated at our study site [18,19,47,48].
The recent advance in the free availability of satellite land surface products at different spatial and temporal resolutions, together with the future launch of sensors combining high resolution and revisit time in the thermal infrared domain (the Copernicus "Land Surface Temperature Monitoring" mission-LSTM onboard Sentinel-8 to be launched in 2028 and the Indo-French "Thermal Infrared Imaging Satellite for High-Resolution Natural Resource Assessment"-TRISHNA mission to be launched in 2024), has further the its potential for water management application. Among those products, mediumresolution Moderate-Resolution Imaging Spectroradiometer (MODIS) products provide globally, since the beginning of the 2000s, most of the physical parameters for surface energy balance models including leaf area index (LAI), albedo, and LST [29][30][31]. They are also freely available from two sensors (Terra and Aqua), thus enhancing temporal resolution. Likewise, meteorological forcings are currently available at the global scale with a good accuracy thanks to the ongoing progress of reanalysis systems. Among them, ERA-Interim [49] reanalysis, a single 32 year global land surface reanalysis dataset covering the period 1979-2010, has proven to be a good candidate for climate monitoring in various part of the world. Jones et al. [50] showed excellent agreement for global-and continental-scale trends in surface air temperatures over land with conventional stationbased datasets. Simmons et al. [51] demonstrated the quality of ERA-Interim near-surface fields by comparing with observation-only climatic data records. The incoming solar radiation provided by the ERA-Interim reanalysis with ground-based measurements was evaluated over France by Szczypta [52]. They showed a slight positive bias, with a modest impact on land surface simulations. These studies demonstrate that reanalyses have the specific advantage of being spatially and temporally complete through the process of physical/dynamic representation of the climate system, which provides internally consistent fields across most surface atmospheric variables, as well as in the atmospheric column up to the stratosphere [53].
The main objective of this work is to present a regional mapping evapotranspiration system based on forcing data available at the global scale and freely distributed (MODIS products and the ERA-Interim reanalysis data). This system is evaluated over the 2001-2015 period thanks to the longtime meteorological and micrometeorological network of the Tensift Observatory [54,55]. As a preliminary step, the quality of the ERA-Interim meteorological forcing is assessed compared to the Tensift Observatory meteorological stations. The predicted components of the energy budget are then compared to eddy covariance station measurements at the time of the satellite overpass. Two complementary indirect assessment are conducted: (1) by comparing predicted monthly ET to precipitation + irrigation water amounts provided at the monthly scale by the office in charge of the agricultural water management in the region (Office Régional de Mise en Valeur du Haouz-ORMVAH) at the scale of the main irrigated perimeter of the catchment; (2) by comparing integrated ET data during the growing season to crop yield statistics. Lastly, the seasonal and interannual variability of the ET in relation to the climatic variability of the region and the agricultural practices are analyzed.

Study Area
The Tensift Al Haouz region, which is considered as a typical watershed of the Southern Mediterranean, is a semiarid watershed (20,450 km 2 ) located in the center of Morocco (region of Marrakech-Safi) ( Figure 1). Located between 30 • 50 and 32 • 1 north and 7 • 12 and 9 • 25 west, the mountainous part of the catchment culminating at 4000 m with the "Jbel Toubkal" is located to the south. The northern border is composed of lower mountains less than 1000 m in elevation called "Jbilet". In the east, the little marked water parting separates the ouadi of Tensift from that of Tessaout and Oum ErR'bia and borders the Atlantic Ocean in the west, known as its release place. Annual crops (wheat, barley, and alfalfa), olive groves, and citrus are the dominant crops (approximately 80% of arable land) in the Tensift area [56][57][58]. The cereal production, which occupies approximately 15,000 km 2 , presents an efficiency of 15 quintals per hectare (qt/ha) on average [59][60][61], a relatively low and variable production from 1 year to the next [62][63][64] with strong discrepancies between the rainfed and irrigated areas, where the production can reach 75 qt/ha. The climate of the watershed is mainly continental semiarid, with specificities on the coastal fringe (influenced by the cold flow form Canary Islands) and on the mountain range of the High Atlas in the south, where most rainfall is concentrated. The rainfall is globally weak in the plain with a strong spatiotemporal variability. The average annual rainfall in the plain is 250 mm, concentrated between November and April, and it can exceed 600 mm in the High Atlas [54]. Please note that R3 stations are not installed on the same fields and crops (see Table 1); however, for ease of visualization, only one point is reported on the map. Table 1. List, location, and acquisition period of the Tensift observatory eddy covariance tower and meteorological stations used in this study for evaluating the predicted energy budget components and for the evaluation of ERA-Interim analysis dataset. * Please note that each equipment and each variable are found at each site. However, for ease, they are reported only once in the

Fields Measurements
The measured data fluxes and meteorological datasets used in this study were provided by the Tensift Observatory composed of a network of meteorological, micrometeorological, and hydrological stations covering the Tensift watershed since 2002. Concerning micrometeorological measurements, most of the eddy covariance stations are located in irrigated areas (named Agafay, Agdal, R3, and Saada; Figure 1), and the instruments (Table 1) are compliant with the FLUXNET standards [65]. They are equipped with commercially available instrumentation: (1) a Krypton hygrometer (model KH2O-Campbell Scientific Inc., Logan, UT, USA) measuring the water vapor sampled at a rate of 20 Hz, and a 3D sonic anemometer (CSAT3, Campbell Scientific Ltd.) measuring the wind speed at high frequency; both data are used to derive the sensible H and latent LE heat fluxes (W/m 2 ) with a half hourly time; (2) air temperature (K) and relative humidity (%) are measured using an HMP45C, Vaisala, Finland; (3) incoming and outgoing solar radiation in the shortand longwave domain is measured using a CNR1 (Kipp & Zonen, Delft, The Netherlands). Details related to the processing, including the quality compared to lysimeter and sap-flow measurement, were given by Diarra et al. [18], Duchemin et al. [66], Elfarkh et al. [47], Ezzahar et al. [67], and Simonneaux et al. [68]. To maintain the integrity of the observations, no gap-filling was performed for these data. Concerning meteorological observations, the network (Table 1) is equipped with standard Campbell Scientific measurement sensors for the measurement of hourly rainfall, global radiation, temperature, relative humidity, and wind speed. The daily latent heat flux observations are computed daily with no missing data by aggregating the 48 half-hourly observations.

MODIS Data
Collection 6 products of Moderate-Resolution Imaging Spectroradiometer (MODIS) onboard the Terra satellite were used to feed the TSEB model. NDVI is provided by MOD13Q1 product, with 250 m/16 days spatiotemporal resolution. Albedo is provided by MCD43A3 product, with 500 m/8 days spatiotemporal resolution. LST from the Terra sensors (MOD11A1 product) is provided at daily and 1 km spatiotemporal resolution at the time over the satellite overpass. Albedos and vegetation indices are interpolated at the daily timestep. The leaf area index (LAI), height of the vegetation, and the fractional vegetation cover (FVC) were obtained from the NDVI thanks to empirical relations. FVC was computed on the basis of the methodology proposed by Gutman and Ignatov [69] and further improved by Zeng et al. [70]. All the variables were finally oversampled at 250 m using the nearest neighbor resampling method ( Figure 2). MODIS products are freely available from Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC): https://ladsweb.modaps.eosdis.nasa.gov/ (accessed on 31 July 2022).

ERA-Interim Dataset
The meteorological variables for the resolution of the surface energy budget are the temperature and the humidity of the air at 2 m, the wind speed, and the incoming radiation in the short-and longwave [33]. In this study, the meteorological variables (air temperature, incoming radiation in the shortwave domain, air humidity, and wind speed) came from ERA-Interim reanalysis [49] for the period of availability of MODIS products (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). The spatial resolution of the dataset is approximately 80 km with a 3 h temporal resolution, and the last version of the ERA-Interim dataset published by the European agency website was used (e.g., Beck et al. [71]).
The air temperature at 2 m and wind speed are key input variables for TSEB (see, among others, the sensitivity study of Bigeard et al. [19]. In particular, the main factors governing the variability of the air temperatures are the elevation and the water and energy exchanges between the surface and the atmosphere, which both depend on soil and vegetation conditions, as well as on advection. In this study, air temperature was disaggregated using the digital elevation model from the Shuttle Radar Topography Mission (SRTM, 30 m) to correct for the difference between the elevation of the last level of the European center model and the real elevation of the surface (Figure 2), following Gao et al. [72].
Likewise, the wind speed was corrected for the difference in elevation assuming a logarithmic profile. The incoming radiation in the longwave was computed using the formula of Brutsaert [73]. The relative humidity was then computed from the temperature at the dew point and from the disaggregated air temperature. Lastly, all the variables were temporally interpolated at the time of the satellite overpass from the ERA-Interim three hourly basis using a simple linear interpolation, and gridded on the MODIS grid at 250 m resolution using the nearest neighbor resampling method.
Long-term estimates of precipitation named climatic research unit (CRU) timeseries are also used in version 3.24 released in 2017 [74]. The CRU precipitation dataset is based on monthly observational data. The data are available at a monthly temporal resolution and at a spatial scale of 0.5 • for the period 1901-present. This dataset was used to indirectly assess the predicted ET at the monthly timescale by comparison to rainfall plus irrigation water amounts. It is preferred to ERA-Interim reanalysis rainfall estimates for two reasons: (1) it is an independent dataset from the ERA-Interim that is used as input of TSEB; (2) significant bias and errors on ERA-Interim rainfall have been highlighted by comparison to observational datasets [75].

Irrigation Water Amounts and Yields Statistics
Ancillary data including monthly irrigation water amounts and cereal grain yields statistics were also gathered to assess indirectly the quality of the ET estimates, particularly their monthly and interannual variabilities, respectively. Irrigation water amounts are provided at the monthly timescale by our partner in charge of the management of the agricultural water for the year 2011. They gather information at the monthly timescale for each irrigation perimeter equipped with concrete channel. The considered irrigated perimeters are reported in Figure 1. Data on cereal productions over the study period 2001-2011 at the administrative provincial scale of the Tensift were kindly made available by the Economic Services of the Ministry of Agriculture in Morocco (https://www.agriculture. gov.ma/, accessed on 31 July 2022), which gathers production of the main cereals cropped in the region (Soft wheat, barley, and durum wheat). More details about the cereal yield data are provided in Balaghi et al. [59], Bouras et al. [60,76], and Jarlan et al. [61].

The Two-Source Energy Balance Model
The TSEB model was described in detail in Norman et al. [35] and Kustas and Norman [77]. Latent and sensible heat fluxes arise from two sources: soil and canopy. The model estimates the turbulent fluxes by simultaneously solving the energy balance composed of two energy balance equations for soil and canopy and two equations of continuity for latent and sensible heat fluxes to the aerodynamic level: where R n is net radiation (W/m 2 ), G is the soil heat flux (W/m 2 ), H is the sensible heat flux (W/m 2 ), and LE is the latent heat flux (W/m 2 ); superscripts s and c denote "soil" and "canopy", respectively. The TSEB model uses a single-time observation of radiometric surface temperature T rad to solve the double-energy balance at the time of the satellite overpass. Soil and vegetation temperatures contribute to T rad in proportion to the fraction of the radiometer view occupied by each component as follows: where T c and T s are the canopy and soil temperatures (K), respectively, f c (θ) is the fractional vegetation cover observed at the radiometer view angle θ. For a canopy with a spherical leaf angle distribution, f c (θ) can be expressed as where Ω, the Kappa factor, indicates the degree to which vegetation is clumped. The Ω value ranges from 0.4 to 0.9 [78]. A value of 0.45 corresponding to no clumping conditions was used in this study. This value was shown to provide the best estimates of LE in our study region by Diarra et al. [18], after a sensitivity analysis. The net radiation R n is computed using where R g and R atm are the short-and longwave incident radiation (W/m 2 ), respectively, the emissivity ε is equal to 0.98, and the albedo α evolves seasonally with the vegetation cover fraction through a weighted average of a soil and a vegetation albedo (taken equal to 0.15 and 0.2, respectively in this study); σ is the Stefan-Boltzman constant (5.67 The net radiation at the soil surface is computed assuming an extinction coefficient of 0.45 [35]. R s n = R n e (−0.45 LAI) .
The leaf area index (LAI) is derived from an empirical relationship using the NDVI (see Section 2.3).
There are two versions of the TSEB model, depending on the connection scheme of the resistance network. The parallel version assumes that both components are located side by side and that the ET from both sources is computed separately. The series version represents the crops as a layer of vegetation above the soil with interactions in terms of radiative and convective fluxes [31]. In this study, the parallel version was selected. With this resistance network, the sensible heat fluxes from the soil H s and canopy H c were computed as follows: where ρ Cp is the volumetric heat capacity (J/m 3 /K) of air, T a is the air temperature (K) measured at 2 m, r ah is the aerodynamic resistance (s/m) to heat transport (e.g., Brutsaert and Reidel [79]), and r sh is the resistance (s/m) to heat flow in the boundary layer immediately above the soil surface. The second formulation for computing H s and H c , by permitting the soil and vegetated canopy fluxes to interact, is the series resistance arrangement (see Appendix A in Norman et al. [35]). The conduction flux in the soil is expressed as a function of the net radiation at the soil surface R s n : were Γ is an empirical coefficient equal to 0.35 [80]. Lastly, the resolution of this set of equations relies on the (strong) assumption that, most of the time, the vegetation transpires at a potential rate. The Priestley-Taylor equation [81], with the parameter α PT nominally equal to 1.26 for unstressed conditions, gives a first estimation of LE c : where R c n is the net radiation to the canopy (W/m 2 ), α PT is the Priestley-Taylor parameter (α PT ≈ 1.26), ∆ is the slope of the saturation vapor pressure-temperature relation (kPa/ • C), Υ is the psychrometric constant (kPa/ • C), and the sign convention is positive toward the canopy.
H c is then computed as a residual of the vegetation energy balance (Equation (2)). T c is derived from H c ; T s is derived from Equation (5); H s is computed from T s and LE s as a residual of the soil energy balance (Equation (1)). If LE s is found negative, indicating condensation of water on the soil surface, which is very unlikely during the day, the assumption of vegetation transpiration at the potential rate is bypassed. LE s is set to zero, and a new value of LE c is computed. Likewise, if LE c is found negative, LE is set to zero. Further details on the resolution are given in French [82] (see Figure 2.6 p.39).

Extrapolation to Daily Value and Monthly Values
In order to extrapolate LE estimates at the time of the satellite overpass to the daily value of ET in mm/day, several approaches have been developed. Delogu et al. [83] among others intercompared different methods of daily extrapolation. Most methods rely on the daily preservation of the evaporative fraction scaling up evapotranspiration with the daily cycle of available energy or are based on the extrapolation of a stress factor with the help of instantaneous potential evapotranspiration value. In this study, a method simply based on the daily course of incoming solar radiation was chosen. The approach was proposed and described in Ryu et al. [84]. This method, particularly simple to implement because incoming solar radiation has a low spatial variability and can be measured easily, is particularly interesting in the view of a practical application by managers. It was evaluated on our region of study by [18], who highlighted a limited RMSE of 0.52 mm in comparison to eddy covariance measurements. In practice, the 48 half-hourly values are cumulated and divided by the water latent heat of vaporization. Monthly values are obtained thanks to a simple average.

Results
As a preliminary step, ERA-Interim meteorological variables were evaluated in comparison to the data acquired by the network of meteorological stations located in the Tensift watershed. The components of the energy budget predicted by TSEB were then compared to the eddy covariance measurements. Table 2 reports the statistical metrics related to the comparison between ERA-Interim reanalysis data after disaggregation and meteorological stations data from the Tensift Observatory at the time of the Terra satellite overpass. The wind speed was unsurprisingly the worst reproduced variable. Indeed, this is one of the meteorological variables, with rainfall, that is most difficult to be retrieved by the weather forecast and climate models, even if the models are driven by observations as for reanalysis systems, because of the strong influence of the local environment (i.e., presence of high vegetation, buildings, a marked orography, and possible advective flux). This difference between ERA-Interim wind speed and weather station data has already been reported by other authors as Liu et al. [85] and Srivastava et al. [86]. In our study, the coefficient of correlation was low (although significant at 99%), but the bias (0.7 m/s) was reasonable and demonstrated the added value of the correction of elevation applied. Indeed, the simplified representation of the orography by ERA-Interim at the 80 km resolution has been shown to be one of the main causes of discrepancy with local observations. Szczypta et al. [87], in particular, who evaluated the ERA-Interim meteorological variables in comparison to the SAFRAN reanalysis over France, highlighted a bias slightly higher (0.9 m/s) except for the mountainous regions.

Evaluation of ERA-Interim
Despite a relatively high RMSE with 116 W/m 2 (Figure 3), the incoming radiation was also reasonably biased (+6 W/m 2 ). This result is in agreement with the results of Szczypta et al. [87], who highlighted a bias of 10 W/m 2 over France. By contrast, the relative humidity is in excellent agreement with the observations with low RMSE and bias. Lastly, the air temperature dynamic, closely linked to that of the air humidity, was also correctly reproduced with r = 0.97 for all the stations ( Table 2). Simmons et al. [88] also showed high levels of correlation (r > 0.94) between ERA-Interim and the CRU dataset for air temperature at the global scale [89]. Szczypta et al. [87] showed that ERA-Interim highlighted a lower diurnal amplitude over France with an overestimation at night (at 12:00 and 6:00 a.m. UTC) and an underestimation around midday. The positive bias of +3 • highlighted here is closely linked to the underlying the network of meteorological station. Indeed, as the Tensift Observatory stations are used to compute the FAO-56 reference evapotranspiration [10], they are installed on irrigated fields with fully developed vegetation cover. This leads to observed temperatures cooler than the average values of the ERA-Interim grid point, which is mostly occupied by hotter areas of bare soils. Several points strengthen this hypothesis. First, this hot bias is as strong as the plant cover is dense. Indeed, on the site of Agafay, which is occupied by fully developed irrigated citrus while the surroundings are mainly bare soil, the bias was 5.2 • C, while it was only 2.4 • C for R3_2012, a station with a much more homogeneous surrounding mainly composed of wheat. Likewise, the Chichaoua station surrounded by bare soil also exhibited a hot bias (4 • C). For the relative humidity, biases were much more limited with a maximum value of +7% on Agdal, and it was reasonably estimated on mostly all stations. However, the errors raised in this study are in line with the literature [72,87,88,90] and remain acceptable.
As a conclusion, the quality of the ERA-Interim dataset compared to meteorological measurement is reasonable, although some discrepancies can be observed in specific conditions.

Energy Budget Component from TSEB
To estimate the quality of the predictions of TSEB forced by ERA-Interim and MODIS products, we evaluated the instantaneous energy components predicted by TSEB at the time of the satellite overpass in comparison to the EC measurements. Firstly, our focus was on the seasonal variation of the four components of the energy budget over two contrasted sites (Figure 4). These were an irrigated olive orchard in 2003 and a rainfed wheat in 2015. On both sites, the seasonal cycle of Rn was very well reproduced with a coefficient r = 0.81 on olive and r = 0.74 on rainfed wheat, although a slight positive bias was observed on both sites. The conduction flux G, in contrast, was rather poorly simulated, especially on the olive orchard with almost no correlation. On the rainfed wheat site, despite a strong overestimation of TSEB with a bias of 111 W/m 2 , the correlation was quite good with r = 0.72. This can be explained by several reasons. In the TSEB model, G is directly related to the energy available at the ground level. The partition of this energy between the ground and the vegetation depends on the fraction of the vegetation cover and a divergence coefficient which is itself linked to the distribution of leaves and the clumping effect of the vegetation. Thus, differences can be significant, especially when the cover fraction is not representative of the surrounding conditions. Indeed, the MODIS NDVI has a resolution of 250 m. In addition, it should be stressed that the divergence coefficient assumption is not always verified, as discussed by Anderson et al. [91], especially over tree crop. Moreover, the soil heat flux is the most difficult variable of the energy budget to estimate. Several studies have pointed out the difficulty of comparing predicted G to point-scale measurements, especially in the study area, with row tree crops leading to contrasted illumination soil conditions [18,48,92].
The performance of the model in reproducing the partition of energy between sensible and latent heat fluxes was very consistent. This partition is strongly dependent on the surface water status, which was very variable from one site to another, between the irrigated perimeters where most experiments were carried out, and the surrounding dry soils. Between the two sites, the model performed better on the rainfed wheat site. The switch between latent and sensible heat fluxes at harvest time (between Q2 and Q3) was, in particular, rather well reproduced. On the olive orchards, we observed an overestimation of LE with a bias of 43 W/m 2 . The overestimation can be explained by the fact that LE is a residual term affected by errors in estimated Rn, G, and H. In addition, the raw EC data, which are known to underestimate LE fluxes by about 20% [14], were used without correction for the energy balance closure. The application of the Bowen ratio method would have certainly reduced this bias by increasing the observed fluxes by about 50 W/m 2 . Zhang et al. [30] reported errors of the same order of magnitude on the LE without using the energy balance closure. In contrast, the sensible heat flux was remarkably well simulated on rainfed wheat. This can be explained by the relative homogeneity of the site around the instrumented plot considering the difference in scale between the footprint of the eddy covariance (around 100 m), on one hand, and the size of the MODIS pixel (250 m) and a spatial representativity larger for ERA-Interim reanalysis even if disaggregation was applied for some variables, on the other hand. Indeed, the rainfed field is located within a rainy area covering about 1 km 2 , while the irrigated olive perimeter is surrounded by dry soils. Likewise, regarding the bias, choosing a lower value of the Priestley-Taylor coefficient could most certainly improve these results because it has been shown that the value of the Priestley-Taylor coefficient should be reduced under dry conditions from its nominal value of 1.26 to a value of 0.8 [44]. An encouraging result, however, is the energy shift observed in Q4 at the Agdal site in 2003. This site was affected by a strong water stress at this time of the year; a technical problem preventing irrigation. This switch between latent and sensible heat flux, although a little too strong compared to the observations, was well observed by the TSEB model. These different points allow us to confirm that the TSEB model can follow the phenology of crops as already mentioned by Chirouze et al. [93], Diarra et al. [18], and Zhang et al. [30].
Considering the good performance of TSEB over the rainfed site, we benchmarked the four components of the energy balance over the period 2001-2015 on all wheat sites ( Figure 5). Results evidence large discrepancies between modeled and observed components, but reasonable statistics were found on average, with a particularly good reproduction of the available energy with a correlation coefficient r = 0.72 and a very slight positive bias of 3 W/m 2 . The G-flux was always the least well reproduced. We also always observed an overestimation of the latent heat flux with a positive bias of 37 W/m 2 and an underestimation of the sensible heat flux with a negative bias of -34 W/m 2 . This is a well-known drawback. As mentioned above, a large part of this difference can be explained by the significant scale mismatch between the EC footprint (on the order of 100 m) and the TSEB predictions. Nevertheless, as reported by Castelli et al. [94], a sensitivity analysis performed on TSEB demonstrated that additional sources of error may be related to inaccuracies in remote sensing and meteorological input data.

Comparison to the Water Inputs + Precipitation
Cumulated CRU precipitation + irrigation water amounts were compared to monthly evapotranspiration estimates for some irrigated perimeters of the Tensift watershed. The monthly evapotranspiration was obtained from the extrapolation of instantaneous estimates at the time of the satellite overpass to daily values converted to mm/day. Table 3 lists the statistics per perimeter and the total statistics for year 2011.  Figure 6 represents the comparison of the monthly ET cumulated at the perimeter scale and the monthly water supply composed of irrigation + CRU rainfall throughout all perimeters during the year 2011. It emerges from the analysis of these results that if the dynamics were correctly reproduced, and a systematic overestimation was highlighted, indicating that ET was higher than the water supply. This systematic overestimation of 60% on average for all irrigated perimeters is further discussed in Section 4.  Table 3.
Another striking result is the good performance of our approach when mapping ET to match water supply for the irrigated perimeters dominated by cereal crops such as the Z1 or R3 perimeter. This result is of course not new because the TSEB model has shown good performance on wheat (see the study case in this manuscript, together with [18,95] among others). Of course, clustering effects affect TSEB performance in areas dominated by tree crop production. The coefficients of correlation are on the same order of magnitude as those obtained by Saadi et al. [3], who led a similar study on a Tunisian irrigated perimeter but using the FAO-56 double crop coefficient approach. The scatter plot of Figure 6 shows a remarkable agreement, but the gaps seem to be more important for large perimeters such as the RG or H2 perimeter. In general, the RMSE appears to be high, but this is largely attributed to the positive bias, as already mentioned. The coefficient of correlation is greatly satisfactory. These initial and very encouraging results were confirmed by investigations in Nassah et al. [96]. However, they deserve to be expanded over the years and analyzed much more finely in collaboration with our management partners regarding the dominant types of crops and known seeding densities.

Comparison to Grain Yield
The evapotranspiration is closely related to crop growth and production, as shown by Semmens et al. [97] for vineyards or Abtew and Melessa [98] and Melesse and Nangia [99] for wheat, among others. This tight linkage is also the basic principle of the Aquacrop model to predict above ground biomass and grain yields [100]. In order to analyze the interannual variability of the ET maps, Figure 7 shows the annual cereal production over the Marrakech agricultural province and the timeseries of calculated ET anomalies over the cereal growth period cumulated from January to April in mm/day over the 2001-2011 period. With a coefficient of correlation r = 0.92, this work shows that the interannual variability of ET integrated over the wheat season is in close agreement with the yield statistics. In particular, the stronger anomalies of cereal production match the stronger anomalies of cumulated ET such as the below-normal season 2007-2008 and the above-normal conditions in 2008-2009 further analyzed in Section 3.3.2. A narrow link exists between evapotranspiration and the water consumed by plants to produce and maintain their raw material and the production in grain. This close relationship has been studied and confirmed by several authors in our study site [60,63,101]. Significant correlations were also highlighted with yields and productions of three cereals taken separately (not shown).

Climatology
On the basis of this performance assessment, we mapped the ET over the entire Tensift watershed, over the period 2001-2015. The instantaneous values at the time of the satellite overpass were aggregated to the daily scale. The daily values were then averaged monthly. Figure 8 shows the monthly climatology of simulated evapotranspiration on the Tensift (2001-2015) watershed during April and August; April represents the peak of plant development for the annuals including the wheat, while August is the driest month during which only the irrigated tree crops are green. Observed spatial variability perfectly fits what we know of the hydrological functioning of the region. The most watered oceanic coasts show significantly higher values of evapotranspiration than the inland region. Likewise, foothills of the Atlas in the southeast part of the Tensift receive higher rainfall than the plain areas, and irrigation can be applied throughout the year by deriving the water of nearby ouadi, which also exhibited higher evapotranspiration throughout the year, as reported by Elfarkh et al. [47] and Ouassanouan et al. [102]. The region on both sides of Marrakech which groups many of the irrigated perimeters and presents high values of evapotranspiration thanks to irrigation. Lastly, for the central region of Tensift, around the city of Chichaoua, which is much less watered than Marrakech, evapotranspiration is much lower. The temporal variability also agrees well with the climate of the region, where most rainfall is concentrated in the winter months. Annual ET is correlated to this cycle, with significant evapotranspiration values in winter spring corresponding to the annual development peak of crops such as wheat and many other crops in the region, as reported by Kharrou [103] and Diarra et al. [18]. The pattern located at the center of the maps, which presents highly persistent evapotranspiration in August, corresponds to irrigated tree crops which are very present over this irrigated perimeter.

Interannual Variability of Evapotranspiration
These monthly ET maps allow us to calculate monthly anomalies of ET over the Tensift watershed. Figure   These favorable conditions were well reproduced on the anomaly maps of February 2009, with values reaching 1.5 mm/day. The most important anomalies are observed in rainfed zones (the west of Marrakech) where the farmers practice opportunist agriculture according to the distribution of rains in the beginning of season. In November, heavy rainfall extensively encouraged the farmers to sow, and the development was regular afterward thanks to the good hydric conditions. This above-normal ET in the rainfed areas persisted up to the peak of development (at the end of March). It was a little less spectacular in the irrigated zones where the contributions of water are guaranteed well enough, even during the loss-making years. By contrast, year 2012 was marked by a deficit of precipitation in spring. Following rather neutral conditions in January, the hydric deficit becomes more marked in reaching its maximum at the time of the peak of development, which often coincides with the phenological filling stage of the grains, during which hydric and thermal stress must be avoided. These strong anomalies at this key stage of the development of the wheat very certainly penalized the production, because it was very low for the 2011-2012 campaign (Economic Services of the Ministry of Agriculture in Morocco: https://www.agriculture.gov.ma/, accessed on 31 July 2022).

Detection of New Irrigated Area
The irrigated zones developed considerably during the last decades mainly because of easy access to underground water, aggravated by periods of recurring droughts, which led to an extension of the amount of drilling in the southern part of the Mediterranean [104]. Especially in the northern zone of Marrakech, we observed a change in the landscape with the appearance of private farms exclusively based on extensive tree crops irrigated through groundwater pumping. The location and the maps of these new irrigated zones are essential information for the study of the watershed hydrological functioning, as well as for the local managers of the water resource. Nevertheless, the maps of irrigated zones resulting from international organizations such as MIRCA2000 [105] or the world map of irrigated zones [106] do not report the fast dynamics of changes in irrigated zones, and they provide coarse-scale information, which limits their utility. Several studies have already mentioned the potential of using ET maps to analyze the evolution of agricultural areas [107] and oases [108]. Our daily regional mapping of ET could allow strengthening these systems of irrigated zone mapping by bringing further information to the existing tools mostly based on spectral indices such as the NDVI. Figure 10  , for Moroccan agriculture that was fostered through incentive funding with conversion to drip irrigation as the main driver, the southeast region of the area was replaced by citrus crop characterized by much higher water needs than more traditional wheat crop [102]. The maps of evapotranspiration allow not only detecting the changes in land use but also quantifying the changes in water consumption, as the average daily consumption in August increased from 2.5-3 to 4-4.5 mm/day.

Discussion
As demonstrated by several sensitivity analysis studies [18,19], meteorological variables including air temperature and wind speed are key inputs of the TSEB model. The new ERA5 [109] and ERA-land reanalysis [110] recently released, as well as the general ongoing improvement of reanalysis datasets, open perspectives for more accurate ET mapping. Indeed, they present a better spatio-temporal resolution (31 km vs. 80 km) leading to significant improvements in the representation of the land surface variables (LSVs) linked to the terrestrial water cycle such us turbulent fluxes [22].
Lastly, this work could be of interest for some important finalized issues. Considering the importance of the production of wheat for the Moroccan economy, it is important for decision makers to be able to have precise information on the production of wheat early in the season to anticipate importation needs and to take up their position on markets [59]. The methodological developments for agricultural yield seasonal forecast follow two main approaches: (1) empirical models relating the yields to environmental indicators (temperature, precipitation, and vegetation index), and (2) crop models fed by the seasonal climatic forecast [111,112]. The latter tools nevertheless suffer from (1) the high demand of input data in the models [112], (2) the uncertainties about the seasonal climatic forecasts [113], and (3) the disparity of the spatial and temporal scales between the climatic forecasts and the expectations of model growth [114]. The limitations evoked above have led most countries to adopt the empirical approach. Most of these studies are based on a series of local predictors such as temperature, rainfall, or information about the crop cover activity through vegetation indices measured from space. Balaghi et al. [59] demonstrated over Morocco that good-quality forecasted yields can be obtained approximately 1 month before harvest. Nevertheless, these approaches, in the south of the Mediterranean area, suffer from a low density of the meteorological network. This work illustrates the potential of evapotranspiration estimates for the early forecast of cereal yields. These daily maps could, for instance, be integrated into key crop phenological stages to develop key predictors closer to the physiological functioning of plants than traditional vegetation indices in order to strengthen early warning systems in the region such as the MEDI system (Mediterranean Drought Index) developed by CESBIO (http://osr-cesbio.ups-tlse.fr/medi/, accessed on 31 July 2022 [115]).
Considering the growing number of farms relying on groundwater pumping in the whole mediterranean region, the daily estimates of ET on a continuous basis provided within this study could be an important information for water managers to quantify groundwater pumping that is largely unknown by lack of accurate drill inventory as most of them are dig outside any legal framework.
In addition, considering the 1 km scale of ET estimates provided in this study, focus was put on the water management at the catchment scale as basin hydraulic agency doesn't need fine information at the field scale. By contrast, another important application of evapotranspiration maps could be irrigation scheduling. Such application is limited considering the sensors acquiring LST into orbit nowadays that provides either data at coarse scale such as MODIS or with a low revisit time such as LANDSAT. Nevertheless, several thermal infrared missions are currently under preparation including TRISHNA and LSTM that will provide both high spatial resolution and revisit time (about 50m and 3 days for TRISHNA). This will open perspective for the application of "energy budget" approaches such as the one presented in this study to irrigation scheduling at the field scale on a day-to-day basis as an alternative of the traditional "water budget" approach based on the FAO-56 method usually used to this objective. Further steps would be to mix both approaches through data assimilation of snapshot estimate of ET into the FAO-56 model.

Conclusions
In this study, we presented a method that allows for daily, spatially distributed, and contiguous estimates of ET. On the basis of the results of Diarra et al. [18], which showed the good performances of the TSEB model in reproducing the energy budget components on the cereal crops in the south of the Mediterranean area, we developed a prototype for the regional cartography of evapotranspiration only based on MODIS and ERA-Interim reanalysis products. These data are freely distributed and are, thus, accessible to all: farmers, watershed managers, and policymakers. This system was directly compared to the Tensift observatory flux tower data (acquired since 2002 on the main crops of the region) and indirectly compared to consumer water on the irrigated surface, provided by the Haouz plain agricultural office (ORMVAH), as well as yields statistics. The main conclusions are as follows: (a) Despite some differences, particularly a positive bias on the air temperature, the ERA-Interim meteorological data constitutes a true solution of forcing of the models compared to the stations of measurements for the areas where meteorological station networks are scarce. (b) TSEB, forced by a low-resolution remote sensing dataset, can reproduce four components of the energy budget reasonably well. Rn is correctly reproduced (r = 0.63 and RMSE = 73 W/m 2 ), while the partition between sensible and latent heat seems very difficult to evaluate given the scale mismatch between the "satellite" predictions and the eddy correlation footprint. Nevertheless, the main trends of the latent heat fluxes seem to be properly reproduced, including the switch of energy occurring after harvest or the occurrence of a stress period. (c) For monthly water consumption of irrigated crops, the approach presents a strong overestimation (almost a factor of 2) compared to the water supply (equal to the irrigation water amounts provided by ORMVAH, the office in charge of the agricultural water in the region plus rainfall from an independent dataset). The first explanation is linked to the trend over the last few years of uncontrolled drilling that has affected the region. Farmers are supplementing ORMVAH water supplies with their own wells in order to secure the harvest, regardless of the crop's actual needs. A similar result was shown in Tunisia by Saadi et al. [3]. (d) Daily mapping allows for anomaly analysis of intra-and interannual climate variability at the watershed scale. This allows understanding the distribution of cultivation practices. (e) This simple tool provides information, supported by yield data, to better understand the variability of grain production. This is information of the highest importance in the current climatic context and also provides a simple approach to detect new irrigated areas.
This modeling framework that does not rely on any in situ data also has potential for application to similar agricultural regions over North Africa, which would facilitate regional water budget assessments, drought monitoring, and other hydrologic research and applications such as identifying the areas where irrigation water management is better conducted than in their poorly efficient counterparts. Lastly, as a wider objective, such simple approaches quite easy to implement could contribute to analyzing the efficiency of irrigation water use and providing objective information to managers to help them rationalize agricultural water. Funding: This work was carried out within the framework of the Joint International Laboratory TREMA (http://lmi-trema.ma, accessed on 31 July 2022). Financial support for the experiment was provided by the Institut de Recherche pour le Développement (IRD), the ERANET-MED 03-62 CHAAMS project, and the H2020 RISE ACCWA project (grant agreement no. 823965). H2020 PRIMA IDEWA is also acknowledged for addition funding. M. Diarra received a PhD grant from the Programme Doctoral International "Modélisation des Systeme complexe" or (PDI MSC) associated with the IRD and the Pierre et Marie Curie University or UPMC (Paris, France).

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.