Estimating Evapotranspiration of Mediterranean Oak Savanna at Multiple Temporal and Spatial Resolutions. Implications for Water Resources Management

: Mediterranean oak savanna is composed of a mixture of scattered oak trees, crops, pasture, and shrubs. It is the most widespread agroforestry landscape in Europe, and its conservation faces multiple threats including water scarcity, which has been exacerbated by global warming and greater climate variability. Evapotranspiration (ET) can be used as a proxy of the vegetation water status and response to water shortage conditions, providing relevant information about the ecosystem stability and its hydrological dynamics. This study evaluates a framework to estimate ET at multiple spatial and temporal scales and applies it to the monitoring of the oak savanna vegetation water consumption for the years 2013–2015. We used a remote sensing-based energy balance model (ALEXI/DisALEXI approach), and the STARFM data fusion technique to provide daily ET estimates at 30 m resolution. The results showed that modeled energy balance components compared well to ground measurements collected by an eddy covariance system, with root mean square error (RMSE) values ranging between 0.60 and 2.18 MJ m − 2 d − 1 , depending on the sensor dataset (MODIS or Landsat) and the ﬂux. The daily 30 m ET series generated by STARFM presented an RMSE value of 0.67 mm d − 1 , which yielded a slight improvement compared to using MODIS resolution or more simple interpolation approaches with Landsat. However, the major advantage of the high spatio-temporal resolution was found in the analysis of ET dynamics over different vegetation patches that shape the landscape structure and create different microclimates. Fine-scale ET maps (30 m, daily) provide key information difﬁcult to detect at a coarser spatial resolution over heterogeneous landscapes and may assist management decisions at the ﬁeld and farm scale.


Introduction
In water-controlled ecosystems with limited water resources, soil moisture dynamics play a central role in the existence and spatial distribution of the different vegetation functional types [1,2]. The discontinuities in the functioning of these ecosystems are related to an alternation in the dry and wet periods [3]. A feedback relationship is observed in these environments, where the vegetation water consumption strongly conditions the hydrological balance of the system, while plants are impacted by water stress situations.
Higher variability in climate patterns coupled with an increase in water demand due to anthropogenic factors is expected to cause a considerable reduction in the quality Within this context, the monitoring of ET with high resolution in both time and space can help to assess landscape structure and to perform a disaggregated evaluation of areas mostly covered by grasslands, and areas with higher tree or bush coverage. The analysis of the spatial and temporal variations in the hydrological conditions controlling the production of pastures and acorns, the two primary sources of livestock feed in the dehesa, can provide information critical for adjusting management practices at the farm scale, with an important impact on the ecosystem conservation and profitability.
Space agencies and companies are making major efforts to distribute high-resolution products, but the current satellite programs present limitations in offering thermal information with high spatial and temporal resolution (e.g., 30 m/daily) simultaneously. An additional problem is the presence of clouds that often result in longer periods without useful images. Multi-scale remotely sensed data fusion techniques are a viable alternative to improve the spatio-temporal resolution of ET estimates derived from thermal-based modeling. The Spatial and Temporal Adaptive Reflectance Fusion Model, STARFM, was developed by Gao et al. [37] to initially integrate surface reflectance data from multiple sensors [38,39], but more recently it has been applied to fuse ET retrievals with satisfactory results in homogeneous agricultural and forested/mosaic areas [16,17,[40][41][42][43]. These studies demonstrated that data fusion modeling improved the accuracy in the estimation of the water consumed by the vegetation, identifying small areas with deficits or excessive irrigation.
This research monitors the vegetation water use in the oak-grass savanna landscape of a small Mediterranean watershed of Southern Spain at different spatial and temporal scales. The specific objectives of the study were: (i) To evaluate the utility of a surface energy balance model (ALEXI/DisALEXI) and the STARFM data fusion technique, using multiple remote sensing platforms (Landsat 7/8 and MODIS), to estimate high-resolution ET in time and space over the complex canopy structure of Mediterranean oak savannas. (ii) To analyze the opportunities offered by this high-resolution product to provide information that is useful to improve the water and vegetation management of this agroforestry system at a field scale. To do that, we evaluated the water use patterns of the herbaceous stratum and other small heterogeneous vegetation patches typical of the dehesa (e.g., scrubs, humid areas, creek shore), which shape the landscape structure and reflect the existence of different micro-ecosystems and climates. Finally, the cumulative monthly ET generated by the different approaches with different spatial resolutions (1 km and 30 m) was quantified over the same vegetation patches.

Description of the Study Area and Experimental Site
The study was conducted over the Martin Gonzalo watershed (48.4 km 2 ), part of the larger Guadalquivir River basin located in Southern Spain (Figure 1a,c). The elevation range varies from 760 m in the north to 280 m at the outlet of the watershed, where there is a dam. The continental Mediterranean climate of this area is highly seasonal, with moderately cold winters alternating with hot and long dry summers. The rainfall presents intra-and inter-year variability, with an annual average of 895 mm (from 1990 to 2015) concentrated during spring and fall.
The landscape is mainly occupied by homogeneous dehesa, along with some conifer forests and olive groves. Dehesa is a Mediterranean oak savanna, whose canopy structure is composed of sparse clumped trees and a grassland or crops understory. In this ecosystem, extensive livestock production fed with acorns and grass is the main economic activity, but they also provide other uses, such as cork, cereal production, hunting, mushroom harvesting, and beekeeping. They offer as well multiple environmental services, such as biodiversity hotspots, water provisioning, CO 2 fixation and high diversity of habitats [28,44]. The landscape is mainly occupied by homogeneous dehesa, along with some conife forests and olive groves. Dehesa is a Mediterranean oak savanna, whose canopy structur is composed of sparse clumped trees and a grassland or crops understory. In this ecosys tem, extensive livestock production fed with acorns and grass is the main economic activ ity, but they also provide other uses, such as cork, cereal production, hunting, mushroom harvesting, and beekeeping. They offer as well multiple environmental services, such a biodiversity hotspots, water provisioning, CO2 fixation and high diversity of habitat [28,44].
Ground validation measurements were taken at an experimental site (Santa Clotild 38°12′ N, 4°17′ W, 736 m a.s.l, Figure 1c), located in a dehesa farm within the studied wa tershed. The setup included an eddy covariance tower (ECT) over the combined tree grassland system and two grazing exclusion enclosures (over open grassland and unde an oak tree respectively) to take into account the heterogeneity of the area (Figure 1b). A energy balance components: the turbulent fluxes of sensible heat (H) and latent heat (LE net radiation (Rn) and the heat flux transport across the surface soil (G), were measure continuously.
The ECT system was installed on a tower at 18 m above ground level in April 2012 registering the ecosystem response as a whole (Figure 2a). The system included a 3D soni anemometer (model CSAT3, Campbell Scientific Inc. Mention of trade names or commer cial products in this publication is solely for the purpose of providing specific informatio and does not imply recommendation or endorsement by the U.S. Department of Agricu ture) that measures horizontal and vertical fluctuations of temperature and wind speed and a hygrometer (model KH20, Campbell Scientific Inc.) to estimate water vapor fluctu ations, the latter being replaced in 2015 by an open Path CO2/H2O Gas Analyzer LICOR 7500, taking simultaneous measurements of carbon dioxide and water vapor in turbulen air structures (Figure 2b). These eddy covariance measurements were recorded with a fre quency of 10 Hz and corrected for density effects due to heat and water vapor transfe [45]. Ground validation measurements were taken at an experimental site (Santa Clotilde, 38 • 12 N, 4 • 17 W, 736 m a.s.l, Figure 1c), located in a dehesa farm within the studied watershed. The setup included an eddy covariance tower (ECT) over the combined tree + grassland system and two grazing exclusion enclosures (over open grassland and under an oak tree respectively) to take into account the heterogeneity of the area (Figure 1b). All energy balance components: the turbulent fluxes of sensible heat (H) and latent heat (LE), net radiation (Rn) and the heat flux transport across the surface soil (G), were measured continuously.
The ECT system was installed on a tower at 18 m above ground level in April 2012, registering the ecosystem response as a whole (Figure 2a). The system included a 3D sonic anemometer (model CSAT3, Campbell Scientific Inc. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the U.S. Department of Agriculture) that measures horizontal and vertical fluctuations of temperature and wind speed, and a hygrometer (model KH20, Campbell Scientific Inc.) to estimate water vapor fluctuations, the latter being replaced in 2015 by an open Path CO 2 /H 2 O Gas Analyzer LICOR-7500, taking simultaneous measurements of carbon dioxide and water vapor in turbulent air structures (Figure 2b). These eddy covariance measurements were recorded with a frequency of 10 Hz and corrected for density effects due to heat and water vapor transfer [45].
A relative humidity and air temperature probe (model HMP155, Vaisala) and a net radiometer (model NR-Lite, Kipp&Zonen, Delft, Netherlands) to obtain the net radiation were also installed at the tower. In 2015, the net radiometer was replaced by a four-component-net-radiation sensor (model NR01, Hukseflux Thermal Sensors, Delft, Netherlands) (Figure 2b). More information on the experimental site and the equipment can be found in Andreu et al. [20] and Carpintero et al. [46].
The quality of the energy balance fluxes measured with the ECT was tested by Andreu et al. [20] for the 2012 summer season, resulting in an average closure balance (Rn − G = LE + H) of 86% of available energy. For the period 2014-2015, the closure balance was 91%. For the comparison with modeled LE, the latent heat fluxes obtained by forcing the closure of the energy balance with the residual method was preferred to direct LE measurement. This method assumes that the sensible heat, H, is correctly measured and LE is obtained by solving the surface energy balance equation. Daytime-integrated energy fluxes were calculated by averaging half-hourly values focusing on daytime fluxes and not considering overnight fluxes. A relative humidity and air temperature probe (model HMP155, Vaisala) and a ne radiometer (model NR-Lite, Kipp&Zonen, Delft, Netherlands) to obtain the net radiatio were also installed at the tower. In 2015, the net radiometer was replaced by a four-com ponent-net-radiation sensor (model NR01, Hukseflux Thermal Sensors, Delft, Nethe lands) (Figure 2b). More information on the experimental site and the equipment can b found in Andreu et al. [20] and Carpintero et al. [46].
The quality of the energy balance fluxes measured with the ECT was tested by An dreu et al. [20] for the 2012 summer season, resulting in an average closure balance (Rn G = LE + H) of 86% of available energy. For the period 2014-2015, the closure balance wa 91%. For the comparison with modeled LE, the latent heat fluxes obtained by forcing th closure of the energy balance with the residual method was preferred to direct LE mea urement. This method assumes that the sensible heat, H, is correctly measured and LE obtained by solving the surface energy balance equation. Daytime-integrated energ fluxes were calculated by averaging half-hourly values focusing on daytime fluxes an not considering overnight fluxes. Figure 3 shows the flowchart of the modeling framework applied for the period 2013 2015. The process starts with a global ET product (5 km, daily) developed using the ALEX model based on MODIS day-night temperature differences [47]. In the second step, th ALEXI ET fluxes are disaggregated to higher resolutions over the study area, using th flux disaggregation scheme DisALEXI applied to both MODIS (1 km, daily) and Landsa 7/8 (60-100 m, 16 days) images. Finally, both types of ET maps, the infrequent Landsat a 30 m resolution and the daily images at MODIS-scale generated by DisALEXI, are com bined using the data fusion technique STARFM to provide ET estimates with both fin spatial (30 m) and temporal (daily) resolution.  Figure 3 shows the flowchart of the modeling framework applied for the period 2013-2015. The process starts with a global ET product (5 km, daily) developed using the ALEXI model based on MODIS day-night temperature differences [47]. In the second step, the ALEXI ET fluxes are disaggregated to higher resolutions over the study area, using the flux disaggregation scheme DisALEXI applied to both MODIS (1 km, daily) and Landsat 7/8 (60-100 m, 16 days) images. Finally, both types of ET maps, the infrequent Landsat at 30 m resolution and the daily images at MODIS-scale generated by DisALEXI, are combined using the data fusion technique STARFM to provide ET estimates with both fine spatial (30 m) and temporal (daily) resolution.
The TSEB model partitions the surface blackbody radiance into the soil/substrate (T S ) and vegetated canopy (T C ) blackbody temperatures, being weighted by the cover fraction of each component at the sensor view angle, f (ϕ), and resulting in: The surface energy balance is solved for the whole soil-canopy-atmosphere system, and for each individual component, with Equations (2) and (3): where the subscripts "s" and "c" represent fluxes from the soil and canopy, respectively. Norman et al. [48] present a description of the original model, and further improvements can be found in Kustas et al. [18].
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 24 Figure 3. Flowchart of modeling framework applied to generate high resolution ET estimations.
The TSEB model partitions the surface blackbody radiance into the soil/substrate (TS) and vegetated canopy (TC) blackbody temperatures, being weighted by the cover fraction of each component at the sensor view angle, f (φ), and resulting in: The surface energy balance is solved for the whole soil-canopy-atmosphere system, and for each individual component, with Equations (2) and (3): where the subscripts "s" and "c" represent fluxes from the soil and canopy, respectively. Norman et al. [48] present a description of the original model, and further improvements can be found in Kustas et al. [18]. The ALEXI model was designed to reduce the use of ancillary meteorological data while maintaining a physically realistic representation of land-atmosphere exchange over a wide range of vegetation coverages. ALEXI implements the TSEB in a time-differential mode, applying this scheme two times during the morning, using surface radiometric temperature data generally provided by geostationary satellites. In this mode, the sensibility of the model to absolute temperature biases is reduced [25]. In ALEXI, TSEB is coupled with an atmospheric boundary layer model to internally simulate the effect of land-atmosphere feedback on the near-surface air temperature [49]. One limitation of this procedure is the dependence on geostationary datasets, with different calibrations and temporal extensions. For this reason, in this research, a global ET product developed using the non-geostationary MODIS sensor day-night temperature differences has been used. As a result, ALEXI generates surface energy fluxes at the continental scale, but at the coarse spatial resolution, of several km.
The DisALEXI scheme runs the TSEB using higher resolution thermal data from polarorbiting satellites, such as MODIS or Landsat, for mapping finer spatial resolution fluxes.
It iteratively adjusts the air temperature boundary conditions such that the disaggregated daily ET flux field reaggregates up to the coarse-resolution ALEXI regional baseline. This spatial disaggregation technique facilitates consistent flux evaluations at local to continental scales [49]. This ALEXI/DisALEXI framework has been widely evaluated across the United States and Europe [49][50][51].

Remote Sensing Data Fusion Method
The STARFM model [37] is a remotely sensed data fusion algorithm that allows an improvement of the temporal representation of ET variations between clear Landsat dates. It was originally designed to blend surface reflectance data, but the fusion of ET maps has been successfully conducted over a number of land covers by Cammalleri et al. [40,41] over rain-fed and irrigated agricultural areas in the Midwestern United States and over irrigated crops in Texas, by Semmens et al. [42] and Knipper et al. [17] over California vineyards, by Anderson et al. [16] over the California Delta and by Yang et al. [43][44][45][46][47][48][49][50][51][52] over forested landscapes. It combines information at a high-temporal frequency from MODIS and high-spatial-resolution from Landsat. STARFM compares one or two Landsat/MODIS image pairs acquired on the same day to predict maps at Landsat spatial scale on other MODIS dates.
Following the weighting function (Equation (4)), STARFM model predicts ET for the central pixel of a selected moving window at generic date (t 0 ): where w is the searching window size; (x w/2 , y w/2 ) is the central pixel of this moving window; n is the number of Landsat and MODIS pairs used (in this case only one pair) and W is the weighting factor. The weighting function (W) integrates the spatial differences between the Landsat (L) and MODIS (M) images on the acquisition date t k and also the temporal differences between MODIS images from observed and predicted dates, t k and t 0 respectively. The prediction for the central pixel only uses spectrally (or ET) similar pixels within the searching window. Further details about the STARFM model are provided by Gao et al. [37].

ET Data Gap Filling
Prior to the fusion algorithm, MODIS and Landsat ET maps produced by DisALEXI were preprocessed to fill the gaps created by clouds or instrument issues (e.g., the failure of scan-line corrector in Landsat-7 ETM+). The gaps in MODIS images were filled to obtain a full daily coverage at 1 km resolution. To accomplish this, the ratio between MODIS and ALEXI ET was computed on MODIS days (f ALEXI = ET MODIS/ET ALEXI), then smoothed and filled using the Savitsky-Golay method described by Sun et al. [53]. The smoothed and filled f ALEXI time series maps were then multiplied by daily ALEXI ET to obtain daily MODIS ET ( Figure 3). The technique used to ensure optimal spatial coverage in Landsat maps employing the STARFM model was described by Yang et al. [43] (Figure 3). It combines the gapped ET image with a prediction from STARFM on the target date using a nearby MODIS-Landsat date pair.

Simple ET Interpolation Methods
In parallel to the application of the STARFM framework to obtain a daily 30 m ET series, daily ET values obtained with a simpler data interpolation method were also generated for comparison purposes. It used the potential ET and MODIS ET as scaling fluxes. The objective was to evaluate the advantages of using the STARFM technique with respect to these simpler and less demanding methods. In this interpolation approach, the ratio between the Landsat ET and a daily scaling flux (potential and MODIS ET) was computed on clear dates. These ratios (FPET = Landsat ET/Potential ET and FMOD = Landsat ET/MODIS ET) were then linearly interpolated in time and fused in the daily ET estimation.

Landsat Data
This work employed a set of 82 Pre-collection and Collection 1 scenes (path 200/row 33) from the Landsat 7 and 8 satellites, acquired from January 2013 to December 2015. Scenes with high cloud coverage were excluded. The Landsat surface reflectance climate data record (SR CDR) (http://espa.cr.usgs.gov/, accessed on 2 September 2017) (distributed atmospherically corrected) was used to compute the leaf area index (LAI) and albedo parameters. A simple approach proposed by Liang et al. [54] was applied to calculate albedo at a 30 m resolution using six surface reflectance bands. To estimate LAI at 30 m on Landsat overpass dates, the regression tree approach developed by Gao et al. [55] was used, which was trained with samples from the MODIS LAI products and Landsat reflectance. This method has been successfully applied to annual and woody crops [42,55]. Landsat thermal band data were corrected for atmospheric and surface emissivity effects using the atmospheric radiative transfer model MODTRAN [56], and subsequently sharpened to the shortwave bands' spatial resolution (30 m) using the Data Mining Sharpener (DMS) method [57].

MODIS Data
MODIS images from Collection 5 (at 1 km spatial resolution) and Collection 6 (500 m) with variable temporal resolution were also used as inputs in the DisALEXI application, including data of land surface temperature (LST; MOD11_L2), geolocation (MOD03), and LAI (MCD15A3). The radiometric temperature with a daily frequency was derived from the 5-min swath LST product, MOD11_L2. The 4-day LAI product, integrated into the TIMESAT algorithm [58], allowed for creating a smoothed and filled daily LAI time series.
As with other DisALEXI application inputs, the global albedo product from Boston University (ftp://rsftp.eeos.umb.edu/data02/Gapfilled/, accessed on 20 September 2017), in addition to the yearly and global MODIS product of land cover (MCD12Q1, 500 m) were included. The albedo data is a global gap-filled snow-free albedo product, which is temporally smoothed and composed over a global grid.

Meteorological Input Data and Vegetation Properties
The energy balance-based modeling framework requires several regional meteorological datasets: a series of surface atmospheric pressure, wind speed at 30 m, and air temperature and specific humidity at 2 m. These were provided every 3 h at a 25 km spatial resolution by the Climate Forecast System Reanalysis (CFSR) [59]. Hourly insolation data at a spatial resolution of 25 km were also obtained from the CFSR. As part of the regional information, a global digital elevation model at a 30 m resolution from ASTER (a product of METI and NASA) and the Andalusian land cover map 2010 were also employed.
A set of biophysical parameters associated with the vegetation in Mediterranean savanna ecosystems were used. The vegetation ground cover fraction, clumping factor and canopy height were calculated following Andreu et al. [20].

Input Data Filtering
In order to isolate errors caused by poor input data quality, the insolation and LST data were filtered for fluxes comparison using ground information. It is known that the quality of insolation data is the primary driver of the total energy available for the ET process [16], and TSEB is mostly driven by thermal information. Both variables were compared with the data measured at the experimental site, showing a mean relative error (mean absolute error/observed mean value; RE) of 9.4% in insolation and 14.5% in surface temperature. The days in which the RE value was higher than 25% in insolation data or higher than 50% in surface temperature were eliminated from the flux comparison analysis. Remotely sensed Remote Sens. 2021, 13, 3701 9 of 23 LST was compared with the half-hourly infrared temperature value corresponding to the time of the satellite overpass. In this case, given the variations between both values and the mismatch representativeness (IR thermometer is mostly measuring canopy- Figure 2b), the threshold considered for filtering the data is less restrictive than for radiation.

Global Remotely Sensed ET Product
The 8-day total ET provided by a product at a 500 m pixel resolution (MOD16A2, Collection 6) was used for comparison purposes. This data collection uses Mu et al.'s [60] improved algorithm. It is based on the logic of the Penman-Monteith equation which includes inputs of daily meteorological reanalysis data along with MODIS data, such as albedo, LAI, fraction of photosynthetically active radiation (FPAR) and land cover.

Evaluation of Surface Energy Fluxes at the Flux Tower Site
The frequency distribution of the wind direction at the ECT for 2015 showed the southwest as the predominant fetch. A footprint analysis was performed by estimating the contribution areas to daily fluxes for 132 days of 2015 as described by Hsieh et al. [61]. An extensive description of the methodology can be found in Kustas et al. [62]. Over the period analyzed, it is observed that 78% of the system ECT fluxes were collected from an area within 130 m upwind (Figure 1b), with the maximum contribution to the measured energy fluxes at approximately 33 m upwind. According to these results, average values of a 3 × 3 grid cell (with 30 × 30 m size-cell and located to the southwest of the ECT) of the modeled fluxes were evaluated.
The relationship between the daily energy fluxes (daytime integrated fluxes) observed at the ECT and modeled by the DisALEXI application is plotted in Figure 4, using MODIS ( Figure 4a) and Landsat images (Figure 4b). The comparison was made on days in which both in-situ measurements and input information for the model were available, with a larger set of usable field data corresponding to 2014 and 2015. In addition, 22 days were filtered from the comparison with DisALEXI/MODIS and 5 days for DisALEXI/Landsat, due to the lack of input data quality (for further information about the input data filtering, see Section 2.3.4). A total of 313 days were used in the comparison for MODIS data and 44 images for Landsat. In all three years, there was a general shortage of clear images and/or field measurements during the last part of the winter and early spring. This is the rainy season for this region and the absence of available data hindered incorporating enough observations representative of this period in the statistical analysis. Statistical metrics, such as mean absolute error (MAE), root mean square error (RMSE), mean bias error (MBE) and coefficient of determination (R 2 ) are presented for all the daily energy fluxes generated by DisALEXI application for the period 2013-2015 (Table 1).     ET measurements acquired by the dehesa ECT system and the rainfall observations. It can be observed in Figure 5 that in this ecosystem (with ~25% fraction tree cover) the annual ET curve is bimodal, with two distinct peaks of different sizes. The largest peak occurs during the spring, reaching maximum values of around 4-5 mm d −1 . A second and smaller peak appears in autumn, with maximum values of around 2-2.5 mm d −1 in 2013-2014 and 1.5 mm d −1 in 2015. This pattern is closely linked to the distribution of the annual rainfall throughout the year, defining the general pattern of soil water availability, and the vegetation growth dynamics, as well as by the energy available for evaporation. The modeled ET reasonably reproduced the temporal dynamics of consumptive water use observed by the flux tower instruments. In general terms, the maximum peaks reached in the spring seasons, as well as the gradual decrease of ET during the end of the spring and summer seasons (periods with higher data available) were accurately reproduced by the STARFM and DisALEXI approaches ( Figure 5).  The modeled ET reasonably reproduced the temporal dynamics of consumptive water use observed by the flux tower instruments. In general terms, the maximum peaks reached in the spring seasons, as well as the gradual decrease of ET during the end of the spring and summer seasons (periods with higher data available) were accurately reproduced by the STARFM and DisALEXI approaches ( Figure 5).

Analysis of ET Time Series from DisALEXI and STARFM
The annual ET modeled by STARFM was equal to 600 mm in the hydrological year 2013/2014 (1 October to 30 September), and 578 mm in the year 2014/2015, with rainfall values of 704 mm and 511 mm for each hydrological year, respectively. These results reflect the evolving vegetation water needs and how the ecosystem adapts to differing water availability conditions, with 193 mm less rainfall during the second year. It also highlights the major role that ET plays in the water balance of the system.
Looking at the details of seasonal water use more closely, we see that the water consumption in the spring of 2015 was similar to that of 2014, as well as the vegetation growth (NDVI values were similar for both years), even though the rainfall was considerably lower (spring rainfall of 120 mm in 2014 and 71 mm in 2015). In this case, the vegetation water consumption in the spring-2015 may have been tapping soil water storage from an unusually high rainfall in the previous autumn (311 mm in autumn-2014 compared to 103 mm of autumn-2013), pointing out the importance of the antecedent soil moisture conditions and the recharge of the subsurface layers.
During the summers, observed ET fell to low values (around 1.1 and 1.4 mm d −1 on average for the months of July and August), with minima around 0.7 mm d −1 at the end of the dry season (in early September before the first rainfalls). In addition to this summer-time soil-moisture limited minimum, a similar shut-down of the vegetation was observed during the winter (average ET around 0.7 in January and February), due to the lack of insolation.
Despite the generally good fit of ET estimations with ECT observations during the study period, some seasonal mismatches between modeled and observed fluxes can be identified. For example, modeled ET increased at the end of the summer in 2013 (on September 3) after a light rainfall event (<2 mm), while ET observations continued to decrease. While the second peak in the observed flux series is missing due to a gap in the field data from September 28, it appears that the rise in ET outlined by Landsat and STARFM preceded the observed rise following rainfall on September 27 (36 mm event). At the end of the dry season of 2014, similar behavior was observed with a gradual increase of modeled estimates starting on September 1, although it was not until September 17 when a significant rainfall event happened (approx. 46 mm). During the wet season in mid-May of 2015, the model showed a lagged response to a drydown that caused the vegetation to dry quickly. This may be related to smoothing/gap-filling in the baseline ALEXI time series, which can tend to smooth out abrupt changes in fPET.
A comparative analysis between the different daily ET series and the observed measurements at the ECT is shown in Table 2 including the quantitative performance of the following daily ET retrievals: (a) from DisALEXI 1 km gap-filled MODIS, (b) STARFM from Landsat-MODIS, and (c) from a simple Landsat-only interpolation technique (using the potential ET as a daily scaling input, through the ratio FPET). The analysis was carried out for a period of 584 days. Results obtained with the Landsat interpolated ET series using MODIS as a daily scaling input (not shown) are similar to those obtained using FPET. The data fusion algorithm for estimating daily 30 m ET outperformed the simpler interpolation technique based on FPET and using MODIS 1-km ET, yielding an RMSE value of 0.67 mm d −1 . Figure 6 shows data of the 8-day total ET layer provided by the global product MOD16A2 (500 m spatial resolution) at the tower pixel, together with the data observed at the ECT and the water consumption modeled by the STARFM approach. The MODIS ET product agreed poorly with the ECT observations (MAE = 0.79 mm d −1 and RMSE = 1.02 mm d −1 ), with a significant underestimation (MBE = −0.56 mm d −1 ). The error was unevenly distributed throughout the year, being more pronounced in the spring and summer seasons, and lessened in autumn and winter. This trend is observed in the three years and it highlights that the use of this product will underestimate ET rates over this semiarid ecosystem. . The error was unevenly distributed throughout the year, being more pronounced in the spring and summer seasons, and lessened in autumn and winter. This trend is observed in the three years and it highlights that the use of this product will underestimate ET rates over this semiarid ecosystem.

Water Resources Management at Field Scale Using High-Resolution ET Maps
At the field scale, the dehesa ecosystem presents variations in the distribution of vegetation strata (oak trees, grassland, crops, or scrubs) that affect the provisioning of ecosystem services and must be considered in the management of the livestock. Within this con-

Water Resources Management at Field Scale Using High-Resolution ET Maps
At the field scale, the dehesa ecosystem presents variations in the distribution of vegetation strata (oak trees, grassland, crops, or scrubs) that affect the provisioning of ecosystem services and must be considered in the management of the livestock. Within this context, the evolution of daily ET at four characteristic micro-ecosystems/patches of the dehesa within the study area (Figure 1c), with different vegetation canopy structures and livestock feed production strategies has been analyzed at two spatial resolutions (30 m and 1 km) to evaluate the potential of Landsat-scale ET retrievals (Figure 7).
Zone A (1700 m 2 ) is an open grassland without trees. It is a grazing area, which exhibits high temporal and spatial variability in vegetation growth. Zone B (2400 m 2 ) is a pasture area where there is a land depression accumulating water through most of the year and maintaining high humidity conditions. Zone C (4300 m 2 ) is a mixture of tree/grass, with an average tree fractional cover of around~25%, common in dehesa production systems and similar to the ECT surrounding area. Acorn production to feed Iberian pigs is usually the main economic use of these landscapes. The final area, zone D (3300 m 2 ), is located near a stream with a high ground fraction covered by trees and an evergreen underlayer of scrubs.
Significant differences were observed between the daily ET provided by DisALEXI/ MODIS at a 1 km resolution and modeled ET with STARFM approach at a 30 m resolution, especially pronounced over zones A and D (Figure 7). Despite the relatively small area of these sites, the 30 m spatial resolution better captured the scale of heterogeneity in this landscape, whereas at 1 km, multiple vegetation covers are contributing, resulting in similar ET patterns for all the sites at this scale.
In the grassland area (Figure 7a), the STARFM 30 m ET decreases more rapidly in the spring than the MODIS-DisALEXI 1 km ET. This decrease corresponds to the drying of the herbaceous stratum, and thus the reduction in the available pasture for forage. Because there are no trees in the STARFM extraction area, transpiration flux is not appreciable during the summer months. The MODIS estimates include contributions from surrounding trees and therefore, maintain significant ET rates during the dry seasons. In each year, the pasture drying process started in May and ended in early to mid-July.
In contrast to zone A, in zone B the late spring decrease in ET is more gradual due to the high soil moisture conditions. The behavior during this season was accurately reproduced by the 30 m resolution model data. In this case, the 30 m and 1 km ET timeseries are more similar to the wet grassland (absent moisture limitations), mimics the transpiration curve of the surrounding trees. The main differences are during springtime peaks, where STARFM reaches values around 4.5-5 mm d −1 , and about 4 mm d −1 with MODIS.
In zone C, the differences between DisALEXI/MODIS and STARFM were not pronounced due to the homogeneity of the surrounding landscape. This type of landscape is representative of the most common structure of the dehesa ecosystem, following a behavior similar to that described in Section 3.2.
Finally, in zone D we see that the 1 km MODIS information is not able to resolve the higher ET focused in this riparian forested area. ET estimates at 30 m suggest the use of groundwater (even shallow water from the creek) by the vegetation during summer, where the ET was maintained over 1.5 mm d −1 , while MODIS ET incorporating surrounding grasslands drops to near 0. A more stable ET trend is seen over this river/scrub zone, where the ET reached similar maximums during the 3 years. These scrub species may have a high nutritional interest both for the domestic livestock, as well as for game species [63].
Cumulative annual water use curves, relevant for water resources management purposes, at 30 m and 1 km resolution are contrasted between relatively wet and dry hydrological years (2013/2014 and 2014/2015, respectively) over the zones A and D in Figure 8. These two areas were chosen because they showed higher scale-dependent differences (    At both resolutions, cumulative ET is higher in the wet year (2013/14) as expected. While cumulative ET at 1 km resolution was similar between zones (around 450-500 mm/year), the 30 m resolution data showed substantial variations between sites, with values ranging between 350 and 900 mm/year.
In the open grassland area, zone A (Figure 8a), the DisALEXI/MODIS approach produced an overestimation of the annual ET in the order of 120-180 mm when compared with the fusion model. Differences were more pronounced during the dry season when the STARFM curve remained stable from May on due to grass senescence. In zone D (Figure 8b), with transpiring vegetation throughout the year, the cumulative ET continuously increased at 30 m resolution. Neither intra-nor inter-annual significant variations in water consumption were observed with STARFM in this area close to a stream, probably due to higher use of groundwater by trees and scrubs and a lesser coupling of water consumption with rainfall events. MODIS cumulative ET showed a slight plateau during the months of July and August, capturing the behavior of the surrounding grassland area in the analysis. An underestimation of approximately 450 mm/year was quantified in this area using 1 km resolution data.
In addition, spatially distributed maps of ET over the Martin Gonzalo watershed and Santa Clotilde dehesa farm (Figure 1c) on three days, representative of different seasons of the period 2013-2015 are presented in Figure 9, highlighting the area with a dehesa ecosystem. They were produced by the ALEXI model (5 km), the DisALEXI/MODIS application (1 km) and the STARFM technique (30 m). In the open grassland area, zone A (Figure 8a), the DisALEXI/MODIS approach produced an overestimation of the annual ET in the order of 120-180 mm when compared with the fusion model. Differences were more pronounced during the dry season when the STARFM curve remained stable from May on due to grass senescence. In zone D (Figure 8b), with transpiring vegetation throughout the year, the cumulative ET continuously increased at 30 m resolution. Neither intra-nor inter-annual significant variations in water consumption were observed with STARFM in this area close to a stream, probably due to higher use of groundwater by trees and scrubs and a lesser coupling of water consumption with rainfall events. MODIS cumulative ET showed a slight plateau during the months of July and August, capturing the behavior of the surrounding grassland area in the analysis. An underestimation of approximately 450 mm/year was quantified in this area using 1 km resolution data.
In addition, spatially distributed maps of ET over the Martin Gonzalo watershed and Santa Clotilde dehesa farm (Figure 1c) on three days, representative of different seasons of the period 2013-2015 are presented in Figure 9, highlighting the area with a dehesa ecosystem. They were produced by the ALEXI model (5 km), the DisALEXI/MODIS application (1 km) and the STARFM technique (30 m).
At the watershed scale, the more detailed representation of ET at a fine resolution provides a better understanding of the variability of the dehesa vegetation water consumption. At 30 m the variations in topography and soils, affecting the vegetation species composition and abundance, can be observed and are largely lost at the 1 km resolution and not existent at the 5 km scale. At the farm scale, the STARFM model allowed for the observation of ET variability in the different vegetation patches previously described, highlighting the existence of micro-ecosystems and climates at the local scale. At the watershed scale, the more detailed representation of ET at a fine resolution provides a better understanding of the variability of the dehesa vegetation water consumption. At 30 m the variations in topography and soils, affecting the vegetation species composition and abundance, can be observed and are largely lost at the 1 km resolution and not existent at the 5 km scale. At the farm scale, the STARFM model allowed for the observation of ET variability in the different vegetation patches previously described, highlighting the existence of micro-ecosystems and climates at the local scale.

DisALEXI Model Validation
The estimation of Rn flux presented an error (RMSE = 1.15 MJ m −2 d −1 ) slightly higher compared to other applications over agricultural areas [40,41,49], but it was aligned with the errors over sparse vegetation cover crops, such as grapevines [17,42] and forested/mosaic areas [16,43]. Ideally, future filtering of insolation inputs needs to be avoided in order to build an automated application, which should include a proxy to account for the quality of the input data. An alternative could be to migrate to an MSG Land-SAF satellite insolation product. Anderson et al. [64] tested different satellite insolation products generated by geostationary satellites, which significantly improved insolation performance over CFSR, although it did not translate into comparable improvement in the ET retrieval accuracy. Likewise

DisALEXI Model Validation
The estimation of Rn flux presented an error (RMSE = 1.15 MJ m −2 d −1 ) slightly higher compared to other applications over agricultural areas [40,41,49], but it was aligned with the errors over sparse vegetation cover crops, such as grapevines [17,42] and forested/mosaic areas [16,43]. Ideally, future filtering of insolation inputs needs to be avoided in order to build an automated application, which should include a proxy to account for the quality of the input data. An alternative could be to migrate to an MSG Land-SAF satellite insolation product. Anderson et al. [64] tested different satellite insolation products generated by geostationary satellites, which significantly improved insolation performance over CFSR, although it did not translate into comparable improvement in the ET retrieval accuracy.
Likewise, the errors in H flux (RMSE value of 1.99 MJ m −2 d −1 when using MODIS data and 2.18 MJ m −2 d −1 with Landsat) were higher than other applications over agricultural areas (in the range of 1.2-1.7 MJ m −2 d −1 observed by Cammalleri et al. [40,41] and Semmens et al. [42]), but in the same range as those obtained in more complex vegetation [17,43].
In relation to LE flux modeled by DisALEXI, a similar behavior over another dehesa site (fluxnet code ES-LMA) using LST data measured in the field and MODIS images, was observed by Andreu et al. [20,21]. These authors used a lower Priestley-Taylor coefficient and stated that even integrating the green fraction and reducing the coefficient, during the dry season the flux was so low as to approach zero. Furthermore, the same trend had been previously observed by Carpintero et al. [65,66] in the study area, where the scheme was applied using a different meteorological dataset source and study period. In this regard, Burchard-Levine et al. [22] implemented the TSEB model over a dehesa area with modifications, considering phenological dynamics of this tree-grass ecosystem and assuming a dominant vegetation structure and cover for different seasons, and different seasonal vegetation parameters. This new adaptation improved the model performance, decreasing the overestimation of LE flux and RMSD errors compared to the application of TSEB assuming a single vegetation source.

Temporal Patterns in ET Curves
The ET time series estimated for the hydrological years 2013/2014 and 2014/2015 reflected the evolution of the vegetation water consumption, depending on water availability. Both hydrological years have similar ET rates, despite the lower precipitation rates in the second one. In 2014/15 water use exceeded the rainfall, suggesting that the vegetation tapped water from other sources to sustain transpiration rates within a threshold that plants could survive. It is known that oak tree roots are able to explore a large volume of soil, with a high dependence on deep water reserves [67,68]. Moreover, oaks form both arbuscular and ectotrophic mycorrhizae, and create mycorrhizal symbioses with partners, adapting to different conditions and accessing many different resources [69].
The similarity of water consumption in spring 2015 with respect to 2014 (with a significantly lower rainfall but an unusually high rainfall in the previous autumn) highlighted the importance of significant antecedent rainfall events to recharge the subsurface layers. In this sense, Fernández [70] found a significant correlation between grassland species abundance in the dehesa and autumn rainfall.
The ET values observed during the summers ( Figure 5; minima values of 0.7 mm d −1 at the end of the dry season) were comparable with David et al.'s results [71], who found that oak trees maintained transpiration rates above 0.7 mm d −1 during the dry season with more than 70% of the transpired water being taken from groundwater sources. In seasons with low LE values (summers and winters), the uncertainty in the data observed by the eddy covariance system [72], added to the fact that LE flux was obtained by closing the balance with the residual method must be considered in the error assessment.
The accuracy of the different modeled daily ET series, shown in Table 2, is below the target error of 0.8 mm d −1 suggested by Seguin et al. [73]. The errors (RMSE = 0.67 mm d −1 with STARFM approach) were close to those found by other authors for woody sparse semiarid crops [21,74,75], and slightly higher than the RMSE value (~0.5 mm d −1 ) found by Campos et al. [76] and Carpintero et al. [46] in a dehesa ecosystem using a locally calibrated remote sensing-based soil water balance.

Performance of MOD16A2 ET
Sriwongsitanon et al. [77] showed a strong agreement of MOD16A2 values with the bulk ET computed using a water balance framework in the forested humid tropics of Thailand, on a monthly and annual scale. The estimations from that study were consistent with the soil moisture conditions and land use classes. An acceptable accuracy of this product was also found by Niyogi et al. [78] and Aguilar et al. [79] for water consumption monitoring over agricultural areas in Indiana (Midwestern United States) and Northwestern Mexico. Comparisons of this product with in-situ measurements at 15 flux tower sites with different climates and biome types ranging from croplands, grasslands, shrublands, savannas, to forests over Europe were evaluated by Hu et al. [80]. The global ET was consistent over most of Europe, with the best results over crops and meadows located in a temperate and humid climate. However, the seasonal performance of MOD16A2 was worse over vegetation under a semiarid climate, with RMSE equal to 1.17 mm d −1 , comparable to the error found in this work ( Figure 6; RMSE = 1.02 mm d −1 ). Global product ET was always underestimated for a dehesa ecosystem in Spain (similar study area and results to those of this work), demonstrating the weakness of this algorithm over limited water availability conditions [80]. The same behavior of ET underestimation was observed over a semiarid region in Iran [81], highlighting the high dependence of MOD16A2 performance with the climate type.
As shown in Figure 6, through the STARFM framework application we were capable of estimating with higher accuracy the temporal trend of ET (8 mm d −1 ) for this semiarid ecosystem, when compared with the MOD16A2 product.

Variability of Dehesa Vegetation Water Use at Field Scale
Dehesa landscape structure features, with variations in the distribution of vegetation strata, are associated with the creation of different microclimates, directly influencing ecological processes, such as plant regeneration and growth, soil respiration, nutrient cycling, and wildlife habitats [35].
The trends observed in the grassland area (zone A; Figure 7a), where the pasture drying process started in May and ended in early to mid-July, were consistent with the pasture production cycle observed in a similar dehesa farm in Spain (Las Majadas). The maximum production was obtained in the spring (around 60-70%) and autumn (15-20%), while the pasture growth was at a minimum in winter (5-15%) and zero in the summer [28]. Information about the pasture drying date at Landsat-scale, for which a high temporal frequency is crucial, can benefit the assessment of the nutritional quality for livestock feed [28]. The ET behavior in zone B (Figure 7b), where high humidity conditions existed, demonstrates the greater transpiration and production capacity of these moist grasslands, important for managing grazing rotations.
Due to the great expansion of the zone C vegetation structure (Figure 7c; a mixture of grass and trees with a low coverage fraction), previous applications of models for estimating water consumption at a low spatial resolution (from 1 km to 5 km) generally have worked well over this ecosystem. These models have provided very useful information at the basin or regional scale [20,21,23]. In zone D, the high spatial resolution ET reflects the presence of dense evergreen vegetation transpiring throughout the year. These areas provide essential ecosystem services, such as shelter for livestock and wildlife (rabbits, Iberian Lynx, etc.), or the creation of specific conditions under the canopy, with different radiation interception, and soils with increased fertility and water retention [82].
The spatial variability in the ET maps is represented in Figure 9 over the Martin Gonzalo watershed (left) and Santa Clotilde farm (right). From a hydrological modeling viewpoint, the 30 m resolution ET product may be very useful to identify and delineate hydrological zones with different water storage capacities and runoff processes within the basin. Nevertheless, the coarse resolution may be enough for regional managing purposes, such as drought monitoring. At the field scale, regular availability of this type of map can lead to developing better agricultural management decisions related to grazing rotations, delineating areas containing fragile ecosystems and important ecological riparian areas with higher water holding capacity and vegetation cover that can sustain a wide diversity of plant and animal species.

Conclusions
The ALEXI/DisALEXI model and the STARFM data fusion technique adequately estimated the ET dynamics of a Mediterranean oak savanna at fine spatial and temporal resolution (30 m, daily) for the period 2013-2015. The energy fluxes provided by DisALEXI using MODIS images (1 km, daily) and Landsat images (60-100 m, 16 days) were compared with in-situ measurements from an eddy covariance flux tower system. The modeled fluxes compared reasonably well to field observations, with RMSE values ranging between 0.60 and 2.18 MJ m −2 d −1 depending on the sensor resolution and the flux component, a similar range to those obtained by other authors over complex landscapes. In addition, the daily 30 m ET series from STARFM model obtained a RMSE value of 0.67 mm d −1 , which is considered an acceptable error for dehesa management purposes.
The ET time series generated by both approaches, DisALEXI and STARFM with different resolutions, showed an annual bimodal behavior of the vegetation water consumption over the dehesa structure (grass and trees, with a tree cover fraction of 25%), with two marked peaks of different magnitudes (around 4-5 mm d −1 during the spring, and between 1.5 and 2.5 mm d −1 in the autumn). Modeled ET accurately reproduced the temporal dynamics of the water consumed by the vegetation, clearly linked to the distribution of the annual rainfall and the energy availability. Nevertheless, some difficulties of the model to reproduce certain abrupt changes in the ET behavior were identified, showing some limitations of the modeling scheme without more frequent observations and the importance of high-quality input data.
The results showed that when the data fusion model was applied, a slight improvement of RMSE, compared to using MODIS or simpler linear interpolations with Landsat data, was observed. This can be due to the combination of daily temporal frequency capturing the changes in ET and the high spatial resolution providing a better representation of the footprint of the flux tower. Modeled annual ET, close to 600 mm for both hydrological years despite their different rainfall regimes, provides an indicative value of current vegetation water needs and how the system adapts to meet these needs under the different water availability conditions. Holm oak's ability to explore a large volume of soil and access groundwater resources under higher than typical water-limited conditions can explain a lack of a difference in ET between years and the minimum modeled values of transpiration of 0.7 mm d −1 at the end of the dry season. However, for periods of long-term and sustained drought, the ecosystem is likely to have a different response.
The major advantage of the data fusion technique and the high spatio-temporal resolution was found in the analysis of ET dynamics over different vegetation microbiomes, characteristics of the dehesa landscape. The results showed that high-resolution ET maps (daily, 30 m) can provide key information to better understand the hydrological functioning of different vegetation distributions difficult to detect at a 1 km pixel resolution and assist in decision making at the farm or local level. However, 1-5 km resolution ET information may suffice for regional monitoring and water planning purposes, such as regional drought detection and impacts on water resources of different biomes.