Daily Actual Evapotranspiration Estimation in a Mediterranean Ecosystem from Landsat Observations Using SEBAL Approach

: Quantifying actual evapotranspiration (ET a ) over natural vegetation is crucial in evaluating the water status of ecosystems and the water-use patterns in local or regional hydrological basins. Remote sensing-based surface energy balance models have been used extensively for estimating ET a in agro-environments; however, the application of these models to natural ecosystems is still limited. The surface energy balance algorithm for land (SEBAL) physical-based surface energy balance model was applied to estimate the actual evapotranspiration over a heterogeneous coverage of Mediterranean maquis in a natural reserve in Sardinia, Italy. The model was applied on 19 Landsat 5 and 8 images from 2009 to 2014, and the results were compared to the data of a micrometeorological station with eddy covariance ﬂux measurements. Comparing the SEBAL-based evaporative fraction ( Λ S ) to the corresponding tower-derived evaporative fractions ( Λ T ) showed good ﬂux estimations in the Landsat overpass time (Coefﬁcient of determination R 2 = 0.77, root mean square error RMSE = 0.05 and mean absolute error MAE = 0.076). Three methods were evaluated for upscaling instantaneous latent heat ﬂux ( λ E) to daily actual evapotranspiration (ET a,D ). The upscaling methods use the evaporative fraction ( Λ ), the reference evapotranspiration fraction (EF r ) and the ratio of daily to instantaneous incoming shortwave radiation (Rs 24 /Rs i ) as upscaling factors under the hypothesis of diurnal self-preservation. A preliminary analysis performed using only in-situ measured data demonstrated that the three factors were relatively self-preserved during the daytime, and can yield good ET a,D estimations, particularly when obtained at near the Landsat scene acquisition time ( ≈ 10:00 UTC). The upscaling factors obtained from SEBAL retrieved instantaneous ﬂuxes, and some ancillary measured meteorological data were used to upscale SEBAL-estimated instantaneous actual λ to daily ET. The Λ EFr and Rs 24 /Rs i methods on average overestimated the measured ET a,D by nearly 20, 61 and 18%, respectively. The performance of the Λ and Rs 24 /Rs i methods was considered satisfactory, bearing in mind the high variable ground cover and the inherent variability of the biome composition, which cannot be properly represented in the Landsat moderate spatial resolution. In this study, we tested the potential of the SEBAL model application in a complex natural ecosystem. This modeling approach will be used to represent the spatial dynamics of ET, which will be integrated into further environmental and hydrological applications. Regarding single acquisition days, some estimation errors can be noticed. In such cases it is difﬁcult to give a deﬁnitive explanation for the source of errors because the results are strongly dependent on the instantaneous energy balance closure, the diurnal variability of this imbalance and the performance of the correction technique used. Nevertheless, we think that such uncertainties cannot be avoided, and the effects of SEB closure and the correction process are beyond the scope of the current study.


Introduction
Understanding eco-hydrological processes is critical for evaluating ecosystem resilience and vulnerability [1]. Evapotranspiration (ET) is one of the most critical components of the land-surface water and energy budget, particularly in arid and semi-arid Mediterranean environments. Chen et al. [2] demonstrated the importance of biotic hydrological processes, particularly the interactions between terrestrial vegetation and hydrology (i.e., canopy interception and ET). The ET variability is a key indicator of ecosystem complexity. Quantifying the spatial and temporal dynamics of ET can help in the assessment of ecosystem sustainability and climate variability. Vincente-Serrano et al. [3] proposed the use of ET and precipitation data as a multi-scalar drought index in the context of global warming. Spano et al. [4] highlighted that, in forest ecosystems, managers need a practical method for continuous, historical and near real-time ET estimation to guide forest management.
Several ground-based ET measurement techniques are available and are widely used, including micrometeorological techniques such as eddy covariance [5][6][7], scintillometry [8,9], the Bowen ratio [10,11] and surface renewal [12,13]. Other ground-based methods measure ET directly by soil-weighing lysimeters [14][15][16], or estimate it by soil water balancing techniques [17,18]. Some methods measure plant transpiration via sap flow [19,20]. Even though field-based methods provide accurate ET measurements, these methods are labor-and cost-intensive [21,22], and require skilled operational staff. Moreover, these measurements are representative of a local spatial scale, from single plant to footprint area. Interpolating or extrapolating such measurements over larger spatial scales can misrepresent the ET variability, particularly for natural environments, where eco-hydrologic studies revealed considerable spatial and temporal variability in plant-water interactions due to microclimatic, land and vegetation heterogeneity [23][24][25]. In this context, there is a need for different methodologies that allow the estimation of ET of natural vegetation over a broader scale.
Knowing that the ET process is governed by the energy exchanges in the soil-plantatmosphere (SPA) system and is limited by the available energy at the land surface, ET can be retrieved by applying the energy conservation principle to the land surface energy balance (SEB). Several models are used to provide an instantaneous estimation of actual ET as a residual term of the land surface energy balance [26][27][28][29][30][31][32][33][34][35]. Such models apply the SEB approach to available free or low-cost earth observation data [21,36]. Satellite-based SEB models have been validated on agricultural land in various parts of the world [26,28,[37][38][39]. The SEB approach for modeling ET spatial distribution was extensively used in irrigation management [26,[40][41][42], water accounting [43], assessing irrigation system performance [36,44,45], agricultural water productivity [46,47], groundwater management [48][49][50], and hydrological modeling [51][52][53][54]. Although the remotely sensed surface energy balance ET approach is widely accepted and used, and despite its ability to provide spatially distributed ET estimation, most of the applications were aimed at hydrological studies and water management in agriculture. According to the author's knowledge, a gap still exists in the application of this approach to natural ecosystems.
Mediterranean forest ecosystems provide multiple environmental and socio-economic goods and services. The Mediterranean maquis is a formation of a complex multistrata structure that represents a dynamic inter-equilibrium, including small trees, tall bushes, dwarf bushes, grasses and microorganisms [55], adapted and resistant to semi-arid condition. In the last few decades Mediterranean forests have been subjected to numerous threats, such as forest fires, over-exploitation, deforestation, and degradation, which are accentuated in the context of climate and land-use changes [56]. Pirastru et al. [57] demonstrated the effects of Mediterranean maquis clearing on soil properties and near-surface hydrology. Folton et al. [58] demonstrated the changes in the hydrological response of a Mediterranean sub-catchment as a result of changes in land physical condition due to wildfire, and Niedda and Pirastru [59] showed the severe effects of land-use and climate changes on the water resources availability in a small Mediterranean catchment in northern Sardinia, Italy.
In this paper, we applied the surface energy balance algorithm for land (SEBAL) model [27,60,61] using Landsat moderate-resolution imagery to estimate ET over a natural Mediterranean ecosystem. SEBAL is a single-source model, and it does not differentiate between components within a single pixel (e.g., soil and vegetation) [62,63]. The model requires surface radiometric temperature derived from thermal infrared radiation on a cloud-free image scene [64,65]. In addition to remote sensing (RS) data, SEBAL relies on a few ground-based ancillary meteorological data [26,60]. The model has been previously and extensively validated under a range of environmental conditions in agricultural land [37,61]. However, there is a need to further verify its potential in natural environments.
SEBAL initially retrieves the net radiation, Rn, soil heat flux at the surface, G 0 , and sensible heat flux, H, at the satellite overpass time, then estimates latent heat flux, λE, as the residual term of the energy budget. The λE represents the energy attributed to the instantaneous actual evapotranspiration. The SEBAL model uses the evaporative fraction (i.e., the ratio of the latent heat flux to the available energy, Rn-G 0 ) to upscale instantaneous estimations of the actual evapotranspiration at the satellite overpass time to daily values [28,36], and operates under the hypothesis of diurnal self-preservation of the evaporative fraction [66][67][68][69]. Alternative upscaling methods that also rely on daytime preservation use the reference crop evapotranspiration [70,71] and the incoming shortwave radiation [72,73] as upscaling parameters.
The general objective of this paper is to evaluate the potential of the SEBAL model combined with Landsat 5 and 8 images in estimating actual ET over a natural ecosystem typical to the Mediterranean coasts. The study is performed over a Mediterranean maquis ecosystem in northwest Sardinia, Italy.
Specifically, this study is carried out with the following objectives: i. testing the diurnal self-preservation hypothesis and different ET upscaling methods based on in-situ flux measurements; ii. evaluating the SEBAL model's performance in estimating instantaneous surface energy fluxes over Mediterranean maquis using the evaporative fraction; iii.
evaluating the upscaling methods for retrieving daily actual ET values from instantaneous SEBAL evapotranspiration estimates of Mediterranean maquis.

Study Area
The study site is located inside a natural reserve called "Le Prigionette", Northwest Sardinia, Italy. The climate is Mediterranean, semi-arid with a warm summer, mild winter, and a high water deficit from May through September. The mean annual temperature is 15.9 • C, the minimum and maximum temperature are 7.0 • C and 28.6 • C, respectively. Precipitation is mainly concentrated from autumn to spring and the annual mean is 588 mm. The soil is 0.3-0.4 m deep Lithic Xerorthent. The ecosystem is a typical coastal Mediterranean maquis (schlerophyll species) mainly consisting of juniper (Juniperus phoenicea L.), lentisk (Pistacia lentiscus L.), tree phyllirea (Phyllirea angustifolia L.), and dwarf fan palm (Chamaerops humilis L.). These species form sparsely vegetated shrub land, in which juniper and lentisk cover 53% and 22% of the vegetated surface, respectively, and they are aggregated into variable-sized patches. Tree phillyrea and palm are found only as isolated elements inside the main patches. Other species present in the experimental site are rosemary (Rosmarinus officinalis L.), Genista corsica (Loisel) DC., Daphne gnidium L., Smilax aspera L., Euphorbia characias L., Helichrysum microphyllum DC., Asphodelus microcarpus Salzm., and Ferula communis L. The vegetation is a secondary succession following a fire event in 1963 and agricultural abandonment in 1970. The mean vegetation height ranges between 0.93 and 1.43 m. Ground cover varies between 42% and 91%.

Micrometeorological and Eddy Covariance Measurements
Eddy Covariance (EC) data were collected by the micro-meteorological station (IT-Noe), installed in the study area as part of the CARBOEUROPE flux monitoring network [74]. The monitoring site is located at the Le Prigionette natural reserve (lat: 40.61, lon: 8.15, 25 m of elevation above sea level) ( Figure 1). The IT-Noe station consists of an EC system with a sonic anemometer (CSAT3, Campbell Scientific, Logan, UT, USA) and an open path gas analyzer (LiCOR 7500, Lincoln, NE, USA), placed within the maquis at 3.5 m above the ground for measuring λE and H. The station also measures G 0 by Data from the IT-Noe station were obtained for the days that correspond to the selected Landsat satellite acquisitions (Table 1). Several studies have reported imbalances in the EC-measured SEB components [75][76][77]. Considering that RS surface energy balance modeling approaches estimate the instantaneous ET under the assumption of the surface energy budget closure, the measured data were corrected by forcing the energy closure using the Bowen ratio. Data quality was evaluated by the energy balance closure (CR) representing the ratio of available energy to the turbulent flux components [75]. As suggested by Prueger et al. [78], CR was assessed for Rn greater than 100 W m −2 : Wilson et al. [75] reported a typical energy budget closure error of 20%. Allen et al. [64] considered that measurements with ±15% energy budget closure error are reliable. Here the energy budget closure was deemed to be satisfactory when CR > 0.85 [75,78]. However, knowing that the SEBAL model is based on energy balance closure and it is more reasonable to compare its results with the closure corrected ground-based data, force closure was applied to all data. The energy budget closure was forced by attributing Data from the IT-Noe station were obtained for the days that correspond to the selected Landsat satellite acquisitions (Table 1). Several studies have reported imbalances in the ECmeasured SEB components [75][76][77]. Considering that RS surface energy balance modeling approaches estimate the instantaneous ET under the assumption of the surface energy budget closure, the measured data were corrected by forcing the energy closure using the Bowen ratio. Data quality was evaluated by the energy balance closure (C R ) representing the ratio of available energy to the turbulent flux components [75]. As suggested by Prueger et al. [78], C R was assessed for Rn greater than 100 W m −2 : Wilson et al. [75] reported a typical energy budget closure error of 20%. Allen et al. [64] considered that measurements with ±15% energy budget closure error are reliable. Here the energy budget closure was deemed to be satisfactory when C R > 0.85 [75,78]. However, knowing that the SEBAL model is based on energy balance closure and it is more reasonable to compare its results with the closure corrected ground-based data, force closure was applied to all data. The energy budget closure was forced by attributing available energy Rn-G 0 to H and λET with the preservation of the Bowen ratio (β), defined as H/λET [78], as follows: where λET c and H c are the corrected values of the latent heat flux and of the sensible heat flux, respectively.   (Table 1). The image selection was driven by both the presence of scene clear-sky acquisitions and good-quality EC tower data. Landsat images and metadata were obtained from the United States Geological Survey (USGS) Earth Explorer (https://earthexplorer.usgs.gov/), last accessed on the 10th of June 2020, and the USGS Global Visualization Viewer (GloVis) (https://glovis.usgs.gov/), last accessed on the 10th of June, 2020. We processed the surface reflectance from visible and infrared bands (near-infrared, mid-infrared and short-wave-infrared) (TM bands 1-5 and 7 and OLI bands 2-7) obtained at 30 m spatial resolution and Landsat thermal bands (TM 6th band and TIRS 10th band). The thermal infrared (TIR) bands were acquired at 120 m and 100 m spatial resolution for TM and TIRS, respectively. TIR data were resampled by the data providers using the cubic convolution method to a 30 m pixel size to match the multispectral bands. Landsat TM and OLI image subsets (cloud-free over study area) were created for the area of interest. Subset images were calibrated for atmospheric scattering and absorption by the "Second Simulation of a Satellite Signal in the Solar Spectrum code" (6S code) [79].

The Surface Energy Balance Algorithm for Land (SEBAL) Model
The SEBAL algorithm is implemented in the ERDAS software Model Maker. SEBAL simulates the radiative and turbulent fluxes within an image pixel without differentiating between soil and vegetation (single-source approach). The SEBAL process flowchart is shown in Figure 2. The Landsat images, accompanied by a 10 m digital elevation model (DEM) (http://www.sardegnageoportale.it/) last accessed on the 15 June 2020 and half-hourly meteorological inputs, were used to produce the SEBAL data layers (i.e., land surface temperature, surface albedo, normalized difference vegetation index, and The SEBAL algorithm is implemented in the ERDAS software Model Maker. SEBAL simulates the radiative and turbulent fluxes within an image pixel without differentiating between soil and vegetation (single-source approach). The SEBAL process flowchart is shown in Figure 2. The Landsat images, accompanied by a 10 m digital elevation model (DEM) (http://www.sardegnageoportale.it/) last accessed on the 15 June 2020 and halfhourly meteorological inputs, were used to produce the SEBAL data layers (i.e., land surface temperature, surface albedo, normalized difference vegetation index, and surface emissivity). SEBAL retrieves Rn, G0 and H and estimates λE as a residue of the energy budget. Figure 2. Flowchart of SEBAL model, main inputs, and output. (U is the wind speed, RH is the relative humidity and Ta is air temperature).
The instantaneous net surface radiation flux Rn (W m −2 ) at the land surface is estimated from Landsat and meteorological data as a balance of incoming and outgoing radiation fluxes [27,61], as follows: where α is surface albedo (-), RS ↓ is incoming short-wave radiation (W m −2 ), RL ↓ is incoming long-wave radiation (W m −2 ), RL ↑ is outgoing long-wave radiation (W m −2 ), and ɛo is the surface emissivity (-). The term (1−ɛo) RL ↓ represents the long-wave radiation reflected from the land surface. The instantaneous value of soil heat flux, G0 (W m −2 ), is estimated from the normalized difference vegetation index, NDVI, the surface radiometric temperature, Trad (°K), α and Rn, as proposed by Bastiaanssen [60]: Usually the difference between the aerodynamic temperature and the air temperature above the canopy is needed for sensible heat flux calculations. In the SEBAL model, a near-surface temperature gradient (ΔT) is used to obtain H [27,61]: The instantaneous net surface radiation flux R n (W m −2 ) at the land surface is estimated from Landsat and meteorological data as a balance of incoming and outgoing radiation fluxes [27,61], as follows: where α is surface albedo (-), R S↓ is incoming short-wave radiation (W m −2 ), R L↓ is incoming long-wave radiation (W m −2 ), R L↑ is outgoing long-wave radiation (W m −2 ), and ε o is the surface emissivity (-). The term (1−ε o ) R L↓ represents the long-wave radiation reflected from the land surface. The instantaneous value of soil heat flux, G 0 (W m −2 ), is estimated from the normalized difference vegetation index, NDVI, the surface radiometric temperature, T rad ( • K), α and Rn, as proposed by Bastiaanssen [60]: Usually the difference between the aerodynamic temperature and the air temperature above the canopy is needed for sensible heat flux calculations. In the SEBAL model, a near-surface temperature gradient (∆T) is used to obtain H [27,61]: where ρ is air density (kg m −3 ), c p is the specific heat of air (J kg −1 K −1 ), r ah is the aerodynamic resistance to heat transport (s m −1 ) and ∆T ( • K) is the temperature gradient between two reference heights, near-surface and air, that governs the transfer of heat. The temperature gradient ∆T is scaled as a function of T rad . The model assumes that an image pixel reproduces a homogenous transfer layer and that ∆T varies linearly with T rad in an image scene [27,61]. This approach captures the relative variability in land surface temperature and reduces the need for absolute accuracy in thermal infrared imagery [80]. Initially, SEBAL solves the surface energy balance equation and obtains the value of ∆T at two hydrological extremes, the anchor pixels, where reliable values of the sensible heat flux can be estimated. The first extreme is referred to as the cold pixel and represents non-stressed full ground cover vegetation, where ET is at the potential level. The second, the hot pixel, is usually a dry bare agricultural soil, where ET is near to zero. The temperature gradient at the cold pixel (∆T cold , Equation (7)) is determined by solving the energy balance assuming the ET is near potential and G 0 is negligible. The same procedure is done for the temperature gradient at the hot pixel (∆T hot , Equation (8)) where the evaporative flux λE hot is considered negligible.
The obtained ∆T values at the anchor pixels are used to compute the correlation coefficient a and b of Equation (9) and estimate the T values from the radiometric temperature (T rad ) across the whole image scene: SEBAL then retrieves the H at each pixel based on ∆T and estimates the instantaneous latent heat flux λE (W m −2 ) as a residual of the surface energy balance:

Upscaling Instantaneous to Daily Evapotranspiration
Three different methods are used to upscale the instantaneous λE at the satellite overpass time to the daily actual evapotranspiration values.
The first method upscales λE based on the self-preservation of the evaporative fraction [41,61,72]. The instantaneous evaporative fraction (Λ S ) is obtained by SEBAL as the ratio of latent heat flux to the available energy at the land surface: The evaporative fraction's self-preservation means that the ratio of the daily energy fluxes is consistent and invariable during daytime. The daily upscaling is a straightforward approach, as the Λ instantaneous value is supposed to be representative of the average daytime value under fair weather conditions (the so-called daytime self-preservation) [28,68,69,81]. The Λ S is used to upscale instantaneous λE to the actual evapotranspiration on the image acquisition day, as follows: where ET D,Λ (mm·d −1 ) is the actual daily evapotranspiration upscaled by Λs, R n24 (MJ·m −2 ·d −1 ) is the net daily radiation estimated by the procedure outlined in the FAO Irrigation and Drainage Paper No. 56 [82], G 24 is the daily soil heat flux considered to be zero on a daily basis by SEBAL, λ (MJ·m −3 ) is the latent heat of vaporization and ρ w (kg·m −3 ) is the water density. The second method uses the standardized reference evapotranspiration as an upscaling variable [70,71]. Both instantaneous reference evapotranspiration, ET r,i (mm·h −1 ), and daily reference evapotranspiration, ET r,24 (mm·d −1 ), were computed by meteorological data following the FAO-56 Penman-Monteith method [82]. This method uses an approach similar to the crop coefficient (K c ) concept [82] and assumes proportionality between the ratio of daily to instantaneous actual and reference ET. This proportionality yields a relatively invariable reference evaporative fraction (EF r ), defined as the ratio of the hourly instantaneous actual evapotranspiration, ET a,i (mm·h −1 ), to the corresponding ET r,i . The upscaling method uses EF r to obtain the daily actual evapotranspiration, as follows: where ET D, EFr (mm d −1 ) is the actual daily evapotranspiration upscaled by EF r . The third upscaling method uses shortwave radiation as the upscaling variable. The incoming shortwave radiation, R S↓ , which is the main source of energy for evapotranspiration at the daily scale [72,73], is used as a reference integration variable, and under the hypothesis of proportionality between changes in ET and R S↓ at the daily scale. This approach also assumes the preservation of the ratio of daily to instantaneous actual ET and the Rs 24 /Rs i ratio during daytime hours: where Rs i is instantaneous incoming shortwave radiation, and Rs 24 is the daily incoming shortwave radiation measured by the micrometeorological station.

SEBAL Data Extraction
The procedure for comparing SEBAL surface fluxes and ground observations considers the scale mismatch between the EC flux footprint and Landsat pixels. The SEBAL values near the acquisition time outputs were extracted for the study area with elliptical polygons (GIS shapefiles) that extend 100 m from the eddy covariance station and fall into line with the measured footprint orientation (upwind direction). The extracted value is the average of a few pixels (from 4 to 10) falling within the polygons. The polygon extension was chosen to partially match the coarser spatial resolution of Landsat TIR and to exclude the interference of the pine forest and the sea near the tower, which exist respectively at 125 and 200 m from the EC tower ( Figure 1).

Measured Energy Fluxes
The wind direction, measured clockwise from the north, ranged from 17.28 • to 355 • , with an average of 218.4 • , during the Landsat image acquisition days (Table 1), indicating that measured turbulent fluxes mainly originated from the North-East. The footprint length representing 90% of the total fluxes (FP 0.90 ) was, on average, 217 m, with a standard deviation of ±58.29 m. The average length of the footprint peak, representing the source area where the majority of turbulent fluxes originates, was approximately 11.43 m with a standard deviation of ±3.07 m. As such, most of the flux contribution originated from an area smaller than a single Landsat image pixel and covered by Mediterranean maquis, which adds confidence to the footprint-pixel data comparison.
The 30 minutes' energy balance components Rn, H, G 0 and λE were analyzed using the flux tower data in all the Landsat image acquisition days (Table 1) (Figure 3). These dates were selected based on the absence of significant data gaps. As shown in Figure 3, Rn exhibited a smooth sinusoidal shape, indicating clear sky for the selected days, with the peak between 12:00 and 14:00. The other fluxes relatively followed the variations of Rn. The sensible heat flux H was the dominant energy component after Rn during the daytime, while λE showed relatively lower values, indicating a limited water availability for the ET process. Negative soil heat fluxes in the afternoon were mainly due to vegetation shading on some soil heat plates. The energy balance closure on the three considered days was acceptable on a daily basis, with an average closure error of less than 15% [64]. The energy balance closure at the near-acquisition time (at 10:00 UTC) was also satisfactory, with C R of 0.835, 1.04 and 1.05 for the DOY 67-2009,  Figure 3b,d,f). Table 2 shows the near-acquisition time C R , the slope of linear regression equation and the coefficient of determination, describing the diurnal variation of turbulent fluxes with respect to available energy on all the selected image acquisition days. The results were similar to those in the three selected days presented in Figure 3. As shown in Figure 3 and Table 2, the 30 min balance closure was sometimes less satisfactory. The partial closure failure in the open field energy balance measurements could be due to different causes. Some terms of the balance, such as the metabolic and the storage terms (i.e., the energy used for biomass production, metabolic activity, and heat storage in vegetation), are usually not accounted for, and therefore contribute to the global error [83]. Another contributing factor is the land heterogeneity [84]. Although the instrumentation functioning and calibration were checked before measurements, errors due to the partial malfunctioning of sensors and/or human causes can contribute to the partial energy balance unclosure. Closure failure also depends on the scale mismatch of the source areas of the energy components and the flux averaging period [85]. Masseroni et al. [86] showed that an averaging period of 30 min may miss the temporal scale variability, particularly in the case of turbulent fluxes.

Diurnal Self-Preservation and Performance of Tower-Derived Upscaling Factors
In this section, we evaluate the validity of the daytime preservation of the upscaling factors, Λ EFr and Rs24/Rsi. In-situ flux observations were used to obtain the diurnal variation of the three considered upscaling factors. The 30 min diurnal variation of these factors is shown in Figure 4 for three Landsat acquisition days, DOY 67-2009, 118-2010 and

Diurnal Self-Preservation and Performance of Tower-Derived Upscaling Factors
In this section, we evaluate the validity of the daytime preservation of the upscaling factors, Λ EF r and Rs 24 /Rs i . In-situ flux observations were used to obtain the diurnal variation of the three considered upscaling factors. The 30 min diurnal variation of these factors is shown in Figure 4 for three Landsat acquisition days, DOY 67-2009, 118-2010 and 257-2014, and similar results were found in the rest of the selected image acquisition days. Regarding the tower-derived evaporative fraction (Λ T ), the diurnal variations showed a typical concave-up shape with two peaks, close to sunrise and sunset. Although λET and the available energy, Rn-G 0 , varied significantly during the day (Figure 3), Λ T showed relative conservation during the daytime, which supports the self-preservation hypothesis approximately from 9:00 to 18:00. Here the Λ T self-preservation can be attributed to several factors, including the relative magnitude of turbulent heat fluxes, and moreover the daytime shape of Λ T was mainly linked to the amplitude and the phase differences between H and λE ( Figure 3). This was confirmed by Gentine et al. [81], who found Λ T mostly independently of the major forcing factors, namely incoming solar radiation and wind speed. The EF r upscaling factor trend was generally in line with the Λ T behavior, as EFr showed relative stability during the daytime, and varied near sunset and sunrise. This is mainly due to the fact that both EF r and Λ T take into account the partitioning of latent energy to the available energy. However, the EF r factor includes additional instantaneous (half-hourly) meteorological variables that are incorporated in the Penman-Monteith equation. The trend of Rs 24 /Rs i also showed low variability and relative consistency during the daytime. Similar to Λ T , Rs 24 /Rs i showed a concave-up shape with peaks near sunrise and sunset. However, compared to Λ T and EF r , the Rs 24 /Rs i showed more considerable diurnal variability, with higher Rs 24 /Rs i values in the early morning (from 7:30 to 8:30 UTC) and late afternoon (from 17:00 to 20:00 UTC) due to relatively low Rs i values at these time intervals. equation. The trend of Rs24/Rsi also showed low variability and relative consistency during the daytime. Similar to ΛT, Rs24/Rsi showed a concave-up shape with peaks near sunrise and sunset. However, compared to ΛT and EFr, the Rs24/Rsi showed more considerable diurnal variability, with higher Rs24/Rsi values in the early morning (from 7:30 to 8:30 UTC) and late afternoon (from 17:00 to 20:00 UTC) due to relatively low Rsi values at these time intervals. With the aim of analyzing the potential of the upscaling factors in retrieving daily ET values, the ΛT, EFr and Rs24/Rsi obtained from the ground-based measurements were used to upscale each of the 30 min measured latent energy values, λE or ETa,i, to the daily actual evapotranspiration (ETa,D). In this analysis, the use of the measured data aims to isolate the upscaling errors from the SEBAL model-specific instantaneous flux estimation errors. Figure 5 shows the diurnal performance of the upscaling methods on the three considered Landsat acquisition days. Each of the estimated ETa,D values was compared to the corresponding measured daily actual evapotranspiration, ETID, EC (black line in Figure 5). The upscaling methods yielded good ETa,D results near the satellite overpass time (at 10:00 UTC) with errors lower than a threshold of ±15% corresponding to EC closure error tolerability [64] (red dashed lines in Figure 5). For example, the near acquisition time ETa,D estimation errors were 9.6, 5.3 and 3.5% for the ΛT method, 0.07, 1.2 and 7.9% for the EFr method and 11.3, 11.3 and 9% for the Rs24/Rsi method on DOY 67-2009, 257-2014 and 118-2010, respectively ( Figure 5). The performance of the three upscaling methods on all the selected image acquisition days (Table 1) was also demonstrated by box and whisker plots ( Figure 6). The results show that beyond the satellite overpass time, the upscaled ETa,D values derived by the EFr and T methods behave similarly, remaining within or near the tollerance limits for most of the daytime, and significantly falling out of those limits near dawn and sunset. Moreover, some outliers were possible during the daytime. An example can be seen at 8:30 and 9:00 of DOY 67-2009 ( Figure 5) due to the significant instantaneous energy closure error and the possible data correction failure where the Bowen ratio could not reproduce the actual energy balance partitioning. Another outlier is shown at 16:00 of With the aim of analyzing the potential of the upscaling factors in retrieving daily ET values, the Λ T , EF r and Rs 24 /Rs i obtained from the ground-based measurements were used to upscale each of the 30 min measured latent energy values, λE or ET a,i , to the daily actual evapotranspiration (ET a,D ). In this analysis, the use of the measured data aims to isolate the upscaling errors from the SEBAL model-specific instantaneous flux estimation errors. Figure 5 shows the diurnal performance of the upscaling methods on the three considered Landsat acquisition days. Each of the estimated ET a,D values was compared to the corresponding measured daily actual evapotranspiration, ET ID, EC (black line in Figure 5). The upscaling methods yielded good ET a,D results near the satellite overpass time (at 10:00 UTC) with errors lower than a threshold of ±15% corresponding to EC closure error tolerability [64] (red dashed lines in Figure 5). For example, the near acquisition time ET a,D estimation errors were 9.6, 5.3 and 3.5% for the Λ T method, 0.07, 1.2 and 7.9% for the EF r method and 11.3, 11.3 and 9% for the Rs 24 /Rs i method on DOY 67-2009, 257-2014 and 118-2010, respectively ( Figure 5). The performance of the three upscaling methods on all the selected image acquisition days (Table 1) was also demonstrated by box and whisker plots ( Figure 6). The results show that beyond the satellite overpass time, the upscaled ET a,D values derived by the EF r and Λ T methods behave similarly, remaining within or near the tollerance limits for most of the daytime, and significantly falling out of those limits near dawn and sunset. Moreover, some outliers were possible during the daytime. An example can be seen at 8:30 and 9:00 of DOY 67-2009 ( Figure 5) due to the significant instantaneous energy closure error and the possible data correction failure where the Bowen ratio could not reproduce the actual energy balance partitioning. Another outlier is shown at 16:00 of DOY 118-2010, and can be attributed to the unexpected shift of phase differences between H and λE observable at that time ( Figure 3). Regarding the Λ T method, the results showed a wide time window where ET upscaling could be accepted. The applicability is mainly limited to daylight hours (i.e., about between 8:30 and 17:00 for the shown dates), as also indicated by other researchers who showed the method's limitations at dawn or sunset, and at night when the energy fluxes are relatively low, and the advective fluxes are strong [87,88]. The EF r method performed similarly; however, EF r showed a wider temporal window of applicability especially in the late afternoon (i.e., between 17:30 and 18:30 DOY 67-209 and 257-2014).This was in line with other studies showing that the EF r method yields good ET a,D estimations even in the afternoon under advective conditions [26,71]. Using the Rs 24 /Rs i values determined during the midday hours slightly underestimated ET a,D . This could be due to the more significant increase in Rs i in this period compared to λE (Equation (14)). On average, the upscaling factors Λ T , EF r and to a lower extent Rs 24 /Rs i yielded good ET a,D estimations. The upscaling methods based on the preservation hypothesis were valid only during daytime, with some exceptions near dawn and sunrise. However, such a limitation generally does not hinder the applicability because at night, except for rare cases, the latent heat flux is negligible [89,90]. Figure 7 compares the daily actual evapotranspiration estimated by upscaling the near overpass time measured latent energy (at 10:00 UTC) to the eddy covariance measured actual daily evapotranspiration (ET ID, EC ). The results showed that the Λ T , EF r and Rs 24 /Rs i scaling approaches yielded good ET a,D estimates. The three methods performed very similarly with coefficients of determination (R 2 ) greater than 0.75. On average, all the three methods yielded ET a,D estimations that fall within the tolerability range of ±15%. The Λ T and Rs 24 /Rs i methods outperformed the EF r method, providing results closer to the equity line, while EF r tends to slightly overestimate ET a,D . The results gave more confidence to the applicability of these upscaling methods near the Landsat image acquisition time (between 9:52 and 10:06 UTC). The performed analysis further supported the hypothesis of the upscaling factors' daytime self-preservation. Moreover, the results indicated that the upscaling variables retrieved near the Landsat overpass time could estimate the daily actual evapotranspiration in this natural Mediterranean ecosystem. Our findings agree with other studies that successfully validated the considered upscaling methods in agricultural lands [66,69,70,72,73]. Regarding single acquisition days, some estimation errors can be noticed. In such cases it is difficult to give a definitive explanation for the source of errors because the results are strongly dependent on the instantaneous energy balance closure, the diurnal variability of this imbalance and the performance of the correction technique used. Nevertheless, we think that such uncertainties cannot be avoided, and the effects of SEB closure and the correction process are beyond the scope of the current study.  [26,71]. Using the Rs24/Rsi values determined during the midday hours slightly underestimated ETa,D This could be due to the more significant increase in Rsi in this period compared to λE (Equation (14)). On average, the upscaling factors ΛT, EFr and to a lower extent Rs24/Rs yielded good ETa,D estimations. The upscaling methods based on the preservation hypoth esis were valid only during daytime, with some exceptions near dawn and sunrise. How ever, such a limitation generally does not hinder the applicability because at night, excep for rare cases, the latent heat flux is negligible [89,90].    can be noticed. In such cases it is difficult to give a definitive explanation for the source o errors because the results are strongly dependent on the instantaneous energy balance closure, the diurnal variability of this imbalance and the performance of the correction technique used. Nevertheless, we think that such uncertainties cannot be avoided, and the effects of SEB closure and the correction process are beyond the scope of the current study

Validation of the SEBAL model
The SEBAL model estimates the instantaneous surface energy balance components at the Landsat overpass time. Generally, the measured fluxes and particularly the H and λE are integrated over time. The IT-Noe station-measured fluxes were integrated and pro vided on a half-hourly basis, which makes the comparison with instantaneous RS-re trieved surface energy fluxes improper. To overcome this inconsistency, similarly to the approach proposed by Bastiaanssen et al. [37], the surface flux ratios that were proven to have diurnal stability were selected as a basis to validate the RS retrieved instantaneous surface energy fluxes. The evaporative fraction estimated by SEBAL (ΛS) was compared to the EC tower-derived evaporative fractions (ΛT) for all the selected satellite acquisitions The results showed that the model estimated Λ with reasonable accuracy, with an R 2 o 0.77, an RMSE of 0.08 and an MAE of 0.05 ( Figure 8).

Validation of the SEBAL Model
The SEBAL model estimates the instantaneous surface energy balance components at the Landsat overpass time. Generally, the measured fluxes and particularly the H and λE are integrated over time. The IT-Noe station-measured fluxes were integrated and provided on a half-hourly basis, which makes the comparison with instantaneous RSretrieved surface energy fluxes improper. To overcome this inconsistency, similarly to the approach proposed by Bastiaanssen et al. [37], the surface flux ratios that were proven to have diurnal stability were selected as a basis to validate the RS retrieved instantaneous surface energy fluxes. The evaporative fraction estimated by SEBAL (Λ S ) was compared to the EC tower-derived evaporative fractions (Λ T ) for all the selected satellite acquisitions. The results showed that the model estimated Λ with reasonable accuracy, with an R 2 of 0.77, an RMSE of 0.08 and an MAE of 0.05 (Figure 8). Considering single acquisition days, we can observe few cases where the ΛS is significantly overestimated. An example is the outlier in Figure 8 that refers to DOY 161 of 2014 (ΛS = 0.6; ΛT = 0.34). On this particular day, in addition to the SEBAL estimation error, the uncertainties can also be attributed to the energy flux measurements at the near-acquisition time and the energy balance closure errors where the Bowen ratio failed in representing real flux partitioning. In fact, even if the daytime energy balance closure was satisfac- Considering single acquisition days, we can observe few cases where the Λ S is significantly overestimated. An example is the outlier in Figure 8 that refers to DOY 161 of 2014 (Λ S = 0.6; Λ T = 0.34). On this particular day, in addition to the SEBAL estimation error, the uncertainties can also be attributed to the energy flux measurements at the near-acquisition time and the energy balance closure errors where the Bowen ratio failed in representing real flux partitioning. In fact, even if the daytime energy balance closure was satisfactory (Cr = 1.1) on DOY 161 of 2014, it was not the case at the near-acquisition time (Cr = 1.36). Excluding this outlier from the data comparison gives R 2 , MAE and RMSE values of 0.87, 0.04 and 0.05, respectively.
The Λ S , EF r and Rs 24 /Rs i factors were used to upscale the SEBAL-estimated instantaneous actual λE to the daily ET. In this case the EF r fraction was obtained from the SEBAL-retrieved instantaneous available energy, near overpass time and daily ancillary measured meteorological data. The Rs 24 was measured by the micrometeorological station and Rs i was estimated from the SEBAL at overpass time. The performance of the upscaling methods is also evaluated here. The upscaled ET a,D was compared to the EC measured daily actual evapotranspiration (ET ID, EC ) ( Figure 9). The Λ and Rs 24 /Rs i methods showed similar performances, as they slightly overestimated the daily actual evapotranspiration. The EFr scaling factor significantly overestimated the ET data, particularly for values greater than 1.5 mm d −1 .
The red dashed line represents the 1:1 line.
Although the SEBAL model showed good instantaneous surface energy flux estima tions at the overpass time (as suggested by the Λ comparison done in Figure 8), we ob serve a decline in performance when comparing the daily upscaled ETa,D to ETID, EC. Here it is important to consider that the ETa,D embeds both instantaneous flux estimations error and upscaling methods errors. In fact, comparing the performance of upscaling factor obtained from SEBAL ( Figure 9) to those obtained from ground-based data (Figure 7 showed that the model-derived factors introduced ETa,D overestimation and increased er ror indices. Besides the errors arising from the model estimations and the upscaling ap proaches, uncertainties in the comparison between simulated data and measured data are dependent on the energy closure errors and the failure of the SEB budget correction. Thi can be manifested in the analysis done by using only in-situ measured data (Figure 7) wherein errors were attributed to the self-preservation of upscaling factors and to both instantaneous and accumulated daily surface energy balance closure errors. In fact, alt hough comparing the overall performance of tower-derived upscaling factors in retriev ing actual daily evapotranspiration (Figure 7) showed a satisfactory performance, single day basis errors cannot be avoided.
We consider positively the results obtained by the Λ and Rs24/Rsi upscaling factors taking into consideration the particular complications of applying an SEB model in such Although the SEBAL model showed good instantaneous surface energy flux estimations at the overpass time (as suggested by the Λ comparison done in Figure 8), we observe a decline in performance when comparing the daily upscaled ET a,D to ET ID, EC . Here, it is important to consider that the ET a,D embeds both instantaneous flux estimations errors and upscaling methods errors. In fact, comparing the performance of upscaling factors obtained from SEBAL ( Figure 9) to those obtained from ground-based data (Figure 7) showed that the model-derived factors introduced ET a,D overestimation and increased error indices. Besides the errors arising from the model estimations and the upscaling approaches, uncertainties in the comparison between simulated data and measured data are dependent on the energy closure errors and the failure of the SEB budget correction. This can be manifested in the analysis done by using only in-situ measured data (Figure 7), wherein errors were attributed to the self-preservation of upscaling factors and to both instantaneous and accumulated daily surface energy balance closure errors. In fact, although comparing the overall performance of tower-derived upscaling factors in retrieving actual daily evapotranspiration (Figure 7) showed a satisfactory performance, single-day basis errors cannot be avoided.
We consider positively the results obtained by the Λ and Rs 24 /Rs i upscaling factors, taking into consideration the particular complications of applying an SEB model in such a complex natural ecosystem in dry conditions. These include difficulties in simulating the actual ET values of surfaces characterized by high heterogeneity (i.e., high variable ground cover), which is typical to the maquis, particularly the maquis of the studied site that is characterized by variable sized patches of mixed species. This surface heterogeneity adds complexity and can be misrepresented in the Landsat 5 and 8 single pixels, especially in the course TIR spatial resolution (120 and 100 m). Some SEB models were developed to partition the energy fluxes and the surface temperature between canopy and soil; an example is the two-source energy balance model (TSEB) model [30,[91][92][93]. The two-source modeling approach may overcome the single-source SEB model, which deals with a pixel as one transfer layer, especially in areas characterized by non-full vegetation cover. However, we think that the complexity and the heterogeneity of the maquis ecosystem is far beyond the partitioning fluxes between soil and vegetation. In fact, the land surface heterogeneity complication even manifests in ground-based measurements. Stoy et al. [77] analyzed energy balance closure across 173 ecosystems in the FLUXNET database and showed that the worst average energy balance closures (0.70-0.78) were in deciduous broadleaf forests, mixed forests and wetlands, due to landscape heterogeneity.

Conclusions
Surface energy balance models combined with free or low-cost remote sensing data were widely used in estimating evapotranspiration in agricultural landscapes. However, an information gap in the application of these models in natural ecosystems still exists. The general objective of this study was to assess the potential of the one-source SEBAL model coupled with Landsat images in simulating the daily actual ET over a coastal Mediterranean maquis.
The SEBAL model firstly estimates the instantaneous latent heat fluxes at the satellite overpass-time, then upscales these instantaneous values to daily values. In this study we used the Λ, ET r and Rs 24 /Rs i upscaling factors under the assumption of diurnal preservation. An analysis performed with in-situ measured data confirmed that the self-preservation of Λ, ET r and Rs 24 /Rs i was assured during daytime. Moreover, the comparison between tower-derived and upscaled daily actual ET with eddy covariance measured daily actual evapotranspiration indicated a significant time window wherein the upscaling approaches were effective, and these were mainly limited to daylight hours. This supported the application of the upscaling approaches in the selected Landsat scene acquisition time (near 10:00 UTC). The comparison of the simulated instantaneous evaporative fractions with the ones obtained from half-hour integrated tower data, provided at the near-overpass time, showed that the SEBAL model could efficiently reproduce the instantaneous energy fluxes. However, when the upscaled daily ET values obtained from the SEBAL instantaneous latent heat fluxes were compared to the measured ET values, the efficiency of the model tended to decrease. This was due to the superimposition of the model instantaneous flux estimation errors with the upscaling methods errors. Moreover, uncertainties in the comparison can be attributed to both instantaneous (at the near overpass time) and daily surface energy balance closure errors. Generally, the Λ and Rs 24 /Rs i were more efficient in estimating daily ET compared to the ET r upscaling approach, which significantly overestimated the daily evapotranspiration. Moreover, compared to the first two methods the EF r approach requires additional meteorological inputs both at the instantaneous (half-hourly or hourly) and the daily scale. Based on our results, the Λ and Rs 24 /Rs i methods are recommended for ET estimation over heterogeneous natural landscapes such as the maquis environment investigated in this study. Furthermore, future applications could also benefit from the remote sensing of incoming shortwave radiation (Rs ↓ ) required in the Λ and Rs 24 /Rs i upscaling methods. Currently, Rs ↓ data can be obtained with relatively high accuracy by geostationary satellites [94,95]. This reduces the dependence on ground data and provides further opportunities for SEBAL large-scale applications in natural areas characterized by data scarcity.
The study demonstrated that the SEBAL model with the Landsat moderate spatial resolution data can estimate the actual evapotranspiration of maquis in a semiarid natural Mediterranean environment with relatively acceptable reliability, taking into consideration the complexity of the studied ecosystem. Such application represents a valuable tool, which better represents the spatial and temporal variability of ET. This spatialized data by itself, or integrated into hydrological modeling, can be a great resource for environmental studies and water resources management.
Author Contributions: H.A. developed the hydrological analysis. M.P., S.M. and C.S. acquired and provided the field data that were processed by H.A., S.D.P. and F.G.; M.P. and D.S. coordinated the team and supervised the research. All the authors analyzed the results and contributed to write and to revise the manuscript. All authors have read and agreed to the published version of the manuscript.