Annual evapotranspiration retrieved from satellite vegetation indices for the eastern Mediterranean at 250 m spatial resolution

We present a model to retrieve actual evapotranspiration (ET) from satellites’ vegetation indices (Parameterization of Vegetation Indices for ET estimation model, or PaVI-E) for the eastern Mediterranean (EM) at a spatial resolution of 250 m. The model is based on the empirical relationship between satellites’ vegetation indices (normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) from MODIS) and total annual ET (ETAnnual) estimated at 16 FLUXNET sites, representing a wide range of plant functional types and ETAnnual. Empirical relationships were first examined separately for (a) annual vegetation systems (i.e. croplands and grasslands) and (b) systems with combined annual and perennial vegetation (i.e. woodlands, forests, savannah and shrublands). Vegetation indices explained most of the variance in ETAnnual in those systems (71 % for annuals, and 88 % for combined annual and perennial systems), while adding land surface temperature data in a multiple-variable regression and a modified version of the Temperature and Greenness model did not result in better correlations (p > 0.1). After establishing empirical relationships, PaVI-E was used to retrieve ETAnnual for the EM from 2000 to 2014. Models’ estimates were highly correlated (R = 0.92, p < 0.01) with ETAnnual calculated from water catchment balances along rainfall gradient of the EM. They were also comparable to the coarser-resolution ET products of the Land Surface Analysis Satellite Applications Facility (LSA-SAF MSG ETa, 3.1 km) and MODIS (MOD16, 1 km) at 148 EM basins with R of 0.75 and 0.77 and relative biases of 5.2 and −5.2 %, respectively (p < 0.001 for both). In the absence of high-resolution (< 1 km) ET models for the EM the proposed model is expected to contribute to the hydrological study of this region, assisting in water resource management, which is one of the most valuable resources of this region.


Introduction
Actual evapotranspiration (ET) is a primary component of the global water cycle.Its assessment at global and regional scales is essential for forecasting future atmospheric feedback (Jung et al., 2010;Oki and Kanae, 2006;Zemp et al., 2014).Estimating ET at such scales, however, is not straightforward and requires the use of models (Chen et al., 2014;Hu et al., 2015;Jung et al., 2009;Trambauer et al., 2014).Data-driven models using satellite information benefit from a continuous spatio-temporal direct observation of the land surface (Ma et al., 2014;Shi and Liang, 2014).
Although physically based models are much more common their performance is comparable to that of the empirically based models (Glenn et al., 2010).The accuracy of both approaches is within that of the eddy covariance mea-Published by Copernicus Publications on behalf of the European Geosciences Union.

D. Helman et al.: High-resolution ET from satellite vegetation index
surements (70-90 %) used for their calibration or validation (Kalma et al., 2008).Yet the empirical approach is simpler than the physically based model and requires less additional information.
The basis for the empirical model is the resource optimization theory.This theory suggests that plants adjust their foliage density to the environmental capacity to support photosynthetic activity and transpiration (Glenn et al., 2010).Accordingly, changes in vegetation foliage cover (and VIs) would affect ET, resulting in high ET-VI correlations.Then, the empirical equation could be used to retrieve ET in space and time.
This approach is mostly used in vegetation systems with an annual cycle of growth and drying where VIs define well the phenological stages (Glenn et al., 2011;Senay et al., 2011).However, in complex systems comprised of annual (i.e.herbaceous) and perennial (i.e.woody) vegetation the model must be adjusted with additional meteorological data (Maselli et al., 2014).
The main drawback of the empirically based approach is that it is limited to a specific site and vegetation type (Glenn et al., 2010;Maselli et al., 2014;Nagler et al., 2012).No common relationship was found between ET and VIs for different sites and climatic conditions.
Here we used MODIS VIs and land surface temperature (LST) products and eddy covariance ET from 16 FLUXNET sites with different plant functional types to establish empirical relationships between VIs (and/or LST) and ET in Mediterranean vegetation systems.We first examined those relationships in annual vegetation systems and complex systems comprising both annual and perennial vegetation.Three empirical models were examined with 16-day and mean annual data: (1) simple and (2) multiple-variable regressions, and (3) a modified version of the Temperature and Greenness (TG) model .We used a performance-simplicity criterion to choose the best model to retrieve ET for the EM.Models' estimates were compared with two ET operational products: MODIS (MOD16) and the Land Surface Analysis Satellite Applications Facility product based on MSG satellites' data (LSA-SAF MSG ETa, hereafter MSG).Finally, we evaluated our model against ET calculated from water catchment balances in the EM.

Evapotranspiration from eddy covariance towers
In situ ET was derived from eddy covariance towers that constitute the international flux towers net (FLUXNET).Two open FLUXNET sources were used to acquire the data sets: the Oak Ridge National Laboratory Distributed Active Archive Centre (available online at http://fluxnet.ornl.gov from ORNL DAAC, Oak Ridge, Tennessee, USA) and the European fluxes database (http://gaia.agraria.unitus.it/home).Half-hourly level 4 ET data were checked for acceptable quality (Reichstein et al., 2005) and gap-filled using methods described in Reichstein et al. (2005) and Moffat et al. (2007).Then, data were aggregated to 16-day means (mm d −1 ) and total annual ET (mm yr −1 ).Only ET data from the time period for which MODIS VI products are available were used (i.e.since 2000).

Satellite products
We used 16-day normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) at a spatial resolution of 250 m (MOD13Q1) and 8-day LST at 1 km spatial resolution (MOD11A2) from MODIS on board the Terra satellite.Although Terra provides LST twice a day (around 10:30 and 22:30 local time) here we used only daytime LST, which is the relevant for ET processes.The 8-day LST was averaged to match the 16-day temporal resolution of the VI product.
The MODIS 16-day VI product is a composite of a singleday value selected from 16-day period based on a maximumvalue criterion (Huete et al., 2002).It represents the vegetation status of the entire 16-day period because of the gradual development of the vegetation.This enables regressing the MODIS VI product against 16-day averages of ET.NDVI is defined as (Rouse et al., 1974 and EVI as (Huete et al., 2002) where R 0.8 , R 0.6 and R 0.5 are the reflectance at near infrared (0.8 µm), red (0.6 µm) and blue (0.5 µm) bands, respectively.NDVI suffers from asymptotic problems (saturation) over high density of vegetation biomass, while EVI is more sensitive in such cases (Huete et al., 2002).For model development, time series of NDVI, EVI and LST at each FLUXNET site were obtained from MODIS Land Product Subsets (http://daac.ornl.gov/MODIS/modis.html) for the years when ET data were available since 2000 (see "Period" column in Table 1).NDVI and EVI time series were smoothed using the local weighted scatter plot smoothing technique (LOWESS) as in Helman et al. (2014aHelman et al. ( , b, 2015)).For model implementation, tiles h20v05, h21v05, h20v06 and h21v06 of the MOD13Q1 product were downloaded for 2000-2014 using the USGS EarthExplorer tool (http://earthexplorer.usgs.gov).These tiles fully cover the eastern Mediterranean (EM) region.
Model results were compared with two satellite operational ET products in 2011 at 148 main basins in the eastern Mediterranean.These operational MODIS and MSG ET products are based on different physical models and have different spatial and temporal resolutions (1 km/8 day for MOD16, and 3.1 km/daily for MSG) (Hu et al., 2015).The annual MODIS (MOD16A3) and daily MSG (LSA-SAF MSG ETa daily) ET products were downloaded for 2011 for the EM region.The basins map was taken from Hy-droSHEDS, a mapping product based on a high-resolution elevation layer developed by the Conservation Science Program of the World Wildlife Fund (http://hydrosheds.cr.usgs.gov).Only main basins with an area greater than 10 km 2 were selected (Fig. S1 in the Supplement).

Evapotranspiration from water catchment balances for validation
We evaluated our model with mean annual ET calculated from water balances at six catchments along a north-south rainfall gradient (130-800 mm yr −1 ) in the eastern Mediterranean (Fig. S2).The calculation follows the classical water balance equation: where P and Q are the total annual precipitation and discharge measured in the catchment, and dS / dt is the change in water storage.Precipitation data (P ) were collected for 2000-2013 from a total of 30 stations of the Israel Meteorological Service: 5 in Kziv, 2 in HaShofet, 21 in the Mountain Aquifer (north, centre and south) and 2 stations in the Mamashit catchment.Data were interpolated for the entire catchment area using ArcGIS and the inverse-distance weighting (IDW) methodology (Lu and Wong, 2008).Discharges (Q) were measured for the same period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) for Kziv, Hashofet and Mamashit catchments using runoff gauges of the Israel Hydrological Service (IHS) at Gesher Haziv hydrometric station for Kziv, HaShofet-Hazorea for HaShofet and Mamashit station for the Mamashit catchment.Annual runoffs at the upper parts of the Mountain Aquifer (drainage areas without hydrometric stations at the Hedera, Alexander, Yarkon, Ayalon, Soreq and Lachish basins) were calculated using the HEC-HMS (Hydrologic Engineering Centre -Hydrologic Modelling System) model (Feldman, 2000) run by the IHS (http://www.water.gov.il).For timescales of several years dS / dt is assumed to be negligible (dS / dt ≈ 0), so the mean annual ET could be simply calculated from P minus Q (Conradt et al., 2013).Following this assumption, we averaged water components over the 14 years of data (i.e. 2000-2013) to calculate mean annual ET.Water balances components (P and Q) and calculated mean annual ET for the six catchments are presented in Table 2.
Calculating ET from water balances has some drawbacks like the difficulty to properly estimate the spatial distribution of precipitation over the entire catchment and uncertainties of catchment boundaries (Conradt et al., 2013).However, this is the best existing approach to compare in situ ET with satellite-derived ET at a basin scale.

Site selection
Perennial and annual vegetation in Mediterranean regions have distinct phenology, contributing differently to the VI signal (Helman et al., 2015;Karnieli, 2003;Lu et al., 2003).
Here we examined VI-ET relationships in vegetation systems comprising both annual and perennial vegetation (i.e.forests, woodlands, savannah and shrublands, hereafter PA) separately from those comprising only annual vegetation (i.e.croplands and grasslands, hereafter AN).
We found that annual vegetation in the understory of PA systems might contribute significantly to VIs while having very small contribution to the total ecosystem ET.In some cases, this results in an apparent phase shift between ET and VIs (Fig. 1), leading to a negative correlation or a lack of correlation.Moreover, we found that AN sites exhibit one single ET-VI relationship under a wide range of rainfall conditions whereas similar types of PA systems have significantly different ET-VI relationships (i.e.different slopes) under different climatic regimes (unpublished results).
Therefore, AN sites (FLUXNET sites in AN systems) were selected from wide range of climatic regimes, while PA sites (FLUXNET sites in PA systems) were selected only from Mediterranean climate regions.Selection of the FLUXNET sites had to fulfil the following criteria: (1) at least 3 years of satellite and eddy covariance data in the FLUXNET site; (2) missing data less than 30 days yr −1 for ET and 15 % for VIs; and (3) homogeneous vegetation cover near the FLUXNET tower within at least the 250 m spatial resolution of the MODIS VI product.The last criterion was manually assured using Google Earth ™ .These led us to select 16 FLUXNET sites that represent a wide range of plant functional types and ET rates (Table 1, Figs.S3 and S4).

Empirical ET models using VIs and LST
Three regression models were examined using the satellitederived NDVI, EVI, LST and eddy covariance ET data: 1. Simple regressions of ET against VIs or LST with 16day and annual data.
2. Multiple-variable regressions using NDVI (or EVI) and LST as independent variables and ET as the dependent variable.Regressions were conducted with both 16-day averages and annual data.
3. Modified version of the TG model proposed by Sims et al. (2008) using LST as a proxy for radiation and potential ET (Maeda et al., 2011) with 16-day data alone.
We used all models with 16-day ET averages and 16-day VIs and/or LST data but only the first two models with total annual ET and mean annual VIs and/or LST because the TG model was designed to work only with 16-day data (Sims et al., 2008).In AN, we subtracted the annual minimum VIs before integrating them over the growing season instead of using the original 16-day VI data (see Helman et al., 2014aHelman et al., , b, 2015)).The integral over the VIs during the growth season was used in the two first models against total annual ET.Multiple-variable regressions were applied only on NDVI and LST data or EVI and LST data, but not on NDVI with EVI data because NDVI and EVI were highly correlated (R > 0.95, p < 0.001).
The original TG model is based on the observed correlations between MODIS-EVI and gross primary production (GPP) measured at several FLUXNET sites, which were further refined by incorporating LST data (Sims et al., 2008): where EVI scaled is the scaled EVI set to 0 at EVI = 0.1 (i.e.EVI scaled = EVI − 0.1) due to absence of photosynthetic activity at this value (Sims et al., 2006); a is the slope of the relationship that enables parameterization of the model; and LST scaled is daytime LST scaled to 1 at an optimum temperature for leaf photosynthetic response around 30 • C, decreasing towards 0 at lower and higher temperatures as follows (Sims et al., 2008): Note that LST scaled in Eq. ( 5) is negative at LST higher than 50 • C. In this case, LST scaled is set to 0 in Eq. ( 4) assuming no photosynthetic activity at those high temperatures due to stomata closure (Sims et al., 2008).
Here, we modified the TG model by using ET instead of GPP in Eq. ( 4): The rationale is that GPP and ET are correlated through the trade-off of carbon gain and water loss by the stomata opening during photosynthetic activity.We used the modified TG model with EVI and NDVI alternatively in Eq. ( 6).

Model evaluation
Pearson's correlation coefficient (R) and mean absolute error (MAE) were chosen as accuracy metrics to evaluate the VI-based ET models.The best model was considered as the one with the highest |R| and lowest MAE or at least lower than the eddy covariance accuracy (< 30 %).If two (or more) models fulfil these requirements, the one with the best performance with respect to its complexity i.e. with respect to the number of variables and operations needed, was preferred.A two-tailed Student's t test was applied to examine statistical differences between the models' p values.

Land cover map for model implementation
ET was assessed for the eastern Mediterranean using the best models for AN and PA systems separately.To produce the required land cover map, we classified pixels as AN and PA based on their NDVI during the year.Low NDVI during the dry season (< 0.25) implies absent or dry vegetation typical for AN systems (Lu et al., 2003).Yet some PA systems (e.g.open shrublands) also have low NDVI during this period but differ from AN systems by smaller NDVI change (< 0.4) during the growth season (Lu et al., 2003;Roderick et al., 1999).Hence, we classified pixels with minimum NDVI < 0.25 as AN only if their NDVI increased by more than 0.4 during the growth season.To account for the high NDVI in agricultural fields of the Nile Delta, pixels with minimum NDVI smaller or equal to 0.35 were also classified as AN only if their NDVI increased by more than 0.35.All remaining pixels were classified as PA (Fig. S5).Although this classification procedure might be coarse, we preferred it to the MODIS land cover product for two reasons.First, a significant discrepancy was found between the MODIS-based land cover product and actual land cover type distribution in the eastern Mediterranean (Sprintsin et al., 2009a).Second, this procedure produces the required mask layer at the spatial resolution of the model (250 m), while the MODIS-derived land cover product is available at a coarser resolution (500 m).
The produced AN/PA land cover map showed the general pattern known for this region (Fig. S5).Moreover, the total AN area estimated for Israel not considering the Golan Heights grasslands (i.e.considering mainly Israel's croplands) was 255 × 10 3 ha.This agreed well with the total cropland area reported by the Israel Central Bureau of Statistics for the same years (220 × 10 3 ha, http://www1.cbs.gov.il).

ET-VI simple relationships in complex systems comprising annual and perennial vegetation
On average, the absolute correlation coefficient (|R|) for the ET-VI linear regressions using annual data was higher by 60 % (for NDVI) and 40 % (for EVI) than the |R| for the 16day regressions at PA sites.Total annual ET was highly correlated with mean annual NDVI at PA sites: 0.85 < R < 0.93 (Table 3; Fig. 2).In contrast, 16-day ET averages were only poorly correlated with 16-day NDVI (0.17 < R < 0.63).The same was for total annual ET and mean annual EVI, with 0.66 < R < 0.94 compared to 0.28 < R < 0.70 when using 16day EVI and ET data.The year-to-year changes in mean annual NDVI and EVI were significant enough to detect even small interannual changes in ET of 20-40 mm yr −1 (see e.g.ES-Amo site in Fig. 2).
Correlation coefficients from the cross-site comparisons were as high as those from site-specific regressions when using annual data in PA sites (Fig. 3).Correlations were equally high for both linear and exponential functions (R = 0.94, p < 0.05 for both VIs and estimating functions).The linear functions were ET = 1277 NDVI -189 and ET = 2844 EVI -300 (mm yr −1 ).Exponential functions were ET = 85 e 3.12 NDVI and ET = 65 e 6.31 EVI (mm yr −1 ).
Although a linear regression function is usually preferred to explain simple relationships between two parameters, the exponential relationship is more realistic in the case of ET-VIs.This is because VIs exhibit exponential relationships with the leaf area index (LAI; Baret et al., 1989;Duchemin et al., 2006), which is directly related to water consumption and ET.Also, ET is usually greater than 0 in places with low vegetation cover (VIs ≤ 0.1) due to soil evaporation.The mean annual NDVI and EVI explained 71 and 88 % of the variability (R 2 ) in the total annual ET using these functions, respectively.This is within the accuracy of the eddy covariance technique for estimating ET (Glenn et al., 2010;Kalma et al., 2008).Cross-site correlations of annual ET and LST in PA were also high with a negative relationship (R = −0.89,p < 0.05, Fig. 3).
The contribution of annual and perennial vegetation to VIs at the sub-pixel level is most difficult to distinguish in PA systems.In some cases, one of those components might have a dominant contribution to VIs, albeit insignificant for the ecosystem flux exchange (Fig. 1).This is probably one of the reasons that VIs could not be used to assess ET at a seasonal timescale (i.e. using 16-day data) in such systems.However, at interannual timescales (i.e. using the annual mean) relationships between ET and VIs were strong and might be used to retrieve total annual ET in PA systems.

Comparison between empirical VI-based ET models
In AN, correlation coefficients from the cross-site regressions of ET against VIs (i.e. the integrals over the growing season period) using the annual data were comparable to those achieved when using the 16-day data (Table 4).The R from the linear regression using 16-day data was as high as 0.86 for both indices (p < 0.001).When using the annual data, R was even higher for ET-NDVI (R = 0.88, p < 0.001), but lower for ET-EVI (R = 0.79, p < 0.001).The mean relative error (i.e.MAE/mean) was substantially lower for regressions using annual data (12-16 %) than for those using the 16-day data (32-33 %, Table 5).The relatively high R for the 16-day ET-VI regressions in AN supports the biomass-ET-VI relationship in those systems, which is described elsewhere (Glenn et al., 2010).
Correlations did not significantly improve (p > 0.1) when LST was added in a multiple-variable regression at the AN sites (Tables 4 and 5).The R from the multiple-variable regressions of LST, VIs and ET was 0.87 when using 16-day data (for LST with each one of the VIs).The R from the multiple-variable regressions on the annual data was 0.89 and 0.79 (ET against LST with NDVI or EVI, respectively, with p < 0.001 for both VIs).
In PA, correlation coefficients from the multiple-variable regressions were substantially higher (p < 0.05 using both VIs) than those obtained from simple ET-VI regressions.R from multiple-variable regressions were 0.71 or 0.73 for 16day ET against LST with NDVI or EVI, respectively, compared to 0.51 and 0.61 for ET against NDVI and EVI.R from the single-and multiple-variable regressions were not statistically different (p > 0.1) in PA when using annual data.The R was 0.94 and 0.96 for multiple-variable models with NDVI and EVI, respectively, and 0.94 for simple regression of ET against VIs (both VIs).
The modified TG model resulted in significantly higher R (p < 0.05 for both indices) only for PA with 16-day data (R = 0.80 and 0.78 using NDVI or EVI in Eq. 6).However, it was still significantly lower (p < 0.05 for both VIs) than the R obtained from simple ET-VI regressions when using the annual data (Table 4 and Fig. S6b).In AN, the R from the modified TG with 16-day data were not significantly different than those obtained from simple ET-VI regressions (p > 0.1, Table 4 and Fig. S6a).

PaVI-E model
NDVI and EVI explained most of the interannual changes in ET in both AN and PA systems (Table 4).This means that a single ET-VI regression function could be used to estimate total annual ET in those systems.Multiple-variable regressions and the modified TG model had higher R and lower MAE in some cases (Tables 4 and 5), but differences were not significant (p > 0.05).Hence, following the performancesimplicity criterion, we chose to use the simple regression functions.The functions obtained from ET-NDVI and ET-EVI regressions were averaged for PA systems as and AN systems as where ET Annual is the total annual ET in millimetres per year (mm yr −1 ).NDVI and EVI in Eq. ( 7) are the mean annual NDVI and EVI.NDVI GSI and EVI GSI in Eq. ( 8  respectively.We used exponential functions because VIs exhibit exponential relationships with LAI, which is directly related to ET, and because ET is greater than 0 in areas with low vegetation cover due to soil evaporation.Finally, we named this model the Parameterization of Vegetation Indices for ET estimation model (PaVI-E).The mean relative error of PaVI-E was 13 and 12 % for AN and PA, respectively.This is within the accuracy of the eddy covariance measurements that were used for calibration and much lower than those reported for more complex models (Glenn et al., 2010;Kalma et al., 2008).PaVI-E was used to assess total annual ET at a spatial resolution of 250 m for the EM after using the land cover map created for AN and PA as a mask layer (Sect.3.3 and Fig. S5).
Figure 4 shows the mean annual ET at the EM for the period of 2000-2014.The annual products of PaVI-E is available by request at 1 km spatial resolution for the entire EM and at 250 m for Israel (http://davidhelman.weebly.com).

Comparison with MOD16 and LSA-SAF MSG ETa
ET estimates from PaVI-E were compared with two operational remote sensing ET products in 148 large basins (> 10 km 2 ).The spatial patterns of annual ET for 2011 from PaVI-E, MOD16 and MSG were generally similar over the EM (Fig. 5).The three models show a general west-to-east and a south-to-north ET gradient along the eastern coastline, matching the rainfall gradients of this region (Ziv et al., 2014).Also, all three models show higher ET estimates over agricultural fields in the Nile Delta compared to the surrounding desert.However, some discrepancies also exist.MOD16 estimates were lower along the EM coast compared to PaVI-E and MSG.ET estimates from MSG were higher along the eastern coast especially to the east of the Galilee Sea (mean ET of ∼ 800 mm yr −1 ).Differences between models were particularly noted over the Nile Delta.Annual ET for 2011 over the Nile Delta was on average 160 mm yr −1 from MSG, 530 mm yr −1 from MOD16, and 680 mm yr −1 from PaVI-E.While MSG estimates seem extremely low for such a highly productive area, PaVI-E and MOD16 estimates agreed well with the high ET reported from in situ measurements (Elhag et al., 2011).Besides the advantage of an improved spatial resolution (250 m compared to 1 km and 3.1 km of MOD16 and MSG) PaVI-E also has the ability to produce spatially continuous ET compared to MSG and MOD16 products (Fig. 5).
Comparing the three models at a basin scale resulted in good agreement between them (R = 0.77 and 0.75 for PaVI-E vs. MOD16 and MSG, respectively, p < 0.001 for both; Fig. 6).MOD16 and MSG had small biases with respect to PaVI-E, with relative biases (i.e.bias/mean) of −5.2 % and 5.2 % and slopes of 0.76 and 1.17 for MOD16 and MSG ET products, respectively.
The relatively higher (lower) MOD16 estimates in xeric (mesic) Mediterranean areas (Fig. 6) was already pointed out by Trambauer et al. ( 2014), who compared this product with several independent ET models.Furthermore, comparison of MOD16 and MSG ET products in Europe showed that correlations with in situ ET (from 15-eddy covariance sites) were better for MSG (Hu et al., 2015), and that MOD16 underestimates ET in Mediterranean dry regions similarly to the observation in this study (Fig. 5).

Evaluation against ET calculated from water catchment balances along rainfall gradient
ET estimates for PaVI-E were evaluated against ET calculated from six water catchments along rainfall gradient in the EM.PaVI-E estimates were highly correlated with the ET calculated from water balances (R = 0.92, p < 0.01) at six catchments along the north-south rainfall gradient in the EM (Fig. 7).ET from MOD16 and MSG was also significantly correlated with the water-balance-derived ET (p < 0.05, Fig. S7).All three models had very similar ET estimates in the Mountain Aquifer catchments (MA-N, MA-CS, and MA-S), lower than those calculated from water balances (Fig. 7) but still within the accuracy of the models (∼ 12 %) and gauging/rainfall distribution uncertainties (∼ 10-15 %; Conradt et al., 2013).
As shown in Fig. 5, ET estimates derived from PaVI-E are significantly higher than those from MOD16 and MSG in the dry areas of the EM.This is due to the exponential functions used in PaVI-E (Eqs.7 and 8).It derived a comparable ET to the calculated from the water balance equation at the dry catchment of Mamashit with a slight overestimation of 15 mm (< 15 %, Fig. 7).MSG largely underestimated the calculated ET in Mamashit (by more than 85 %), while MOD16 had no data for this area.

Conclusions
Three empirical VI-based ET models using only eddy covariance ET and MODIS vegetation indices and land surface temperature data for Mediterranean vegetation systems were tested.VIs in vegetation systems comprising mostly annual vegetation (i.e.grasslands and croplands) had strong relationships with intra-annual (16-day ET averages) and interannual (total annual ET) ET estimates.The mean relative error was larger for intra-annual relationships compared to interannual relationships (32 compared to 12 %).In complex systems with annual and perennial vegetation (i.e.forests, woodlands, savannah and shrublands) ET-VI relationships were strong only at interannual timescales (i.e. using annual data).Intra-annual relationships were poor probably due to the mixed VI signal contributed by annual and perennial vegetation that constitute different vertical layers in those sys-tems.While annual vegetation (mostly herbaceous vegetation in the understory) is the main contributor to the intraannual VI change, it constitutes only a minor contributor to the total ecosystem ET in complex Mediterranean systems.Multiple-variable regression and a modified version of the TG model with VIs and LST were not significantly better than the simple ET-VI model for both PA and AN vegetation systems (p > 0.1).
The empirical ET-VI model, PaVI-E, had comparable estimates to MOD16 and LSA-SAF MSG ET models in the eastern Mediterranean.PaVI-E also agreed well with ET calculated using the water balance equation at six catchments along the south-north EM rainfall gradient.PaVI-E is the first ET model with such high resolution (250 m) for this region.Its advantage is in its simplicity and spatial resolution compared to the coarser resolutions of MOD16 and LSA-SAF MSG ETa products (1 and 3.1 km, respectively).We are confident that using PaVI-E will enhance the hydrological study in this region, where ET plays a major role in the hydrological cycle.
The Supplement related to this article is available online at doi:10.5194/acp-15-12567-2015-supplement.

Figure 1 .
Figure 1.Sixteen-day eddy covariance ET averages (black dots) and MODIS-derived NDVI (green line) at two vegetation systems: (top) PA, i.e. comprising perennial and annual vegetation (evergreen coniferous forest), and (bottom) AN, i.e.only annual vegetation (corn and soybean cropland).Note: in the cropland site (bottom) is the NDVI during the growing season after the annual minimum was subtracted.

Figure 2 .
Figure 2. Relationships between annual ET (mm yr −1 ) from eddy covariance towers and mean annual MODIS-derived NDVI, EVI and LST ( • C) at PA sites (perennials and annuals vegetation systems, i.e. forests, woodlands, savannah and shrublands).

Figure 3 .
Figure 3. Same as Fig. 2 but for all PA sites together.The linear (dashed line) and exponential (solid line) functions are presented for the ET-VI relationships, and the R is for the exponential function.

Figure 5 .
Figure 5.Total annual ET for the eastern Mediterranean from PaVI-E, MODIS (MOD16) and LSA-SAF MSG ETa (MSG) for 2011.Grey colouring in MOD16 and MSG indicates missing data.

Figure 7 .
Figure 7. (Left) Scatter plot of the mean annual ET (2000-2013) retrieved from PaVI-E and calculated using the water balance equation at six catchments along the EM north-south rainfall gradient (Fig. S2).(Right) Comparison between mean annual ET estimates from PaVI-E, MOD16, MSG and the water balances in those six water catchments.MA-N, MA-CS and MA-S stand for the northern, central-southern, and southern parts of the Mountain Aquifer of Israel, respectively.

Table 1 .
Description of the 16 selected FLUXNET sites.Horizontal line divides between the six FLUXNET sites in PA systems (top) and the nine FLUXNET sites in AN systems (bottom).Plant functional types (PFTs) include CSH: closed shrublands; WDL: woodland; SAV: savannah; ENF: evergreen needle-leaved forest; WSA: woody savannah; CRO: croplands; and GRA: grasslands.Mean annual precipitation (P ) is in millimetres per year (mm yr −1 ) for the years in which ET data were used ("Period").

Table 3 .
Correlation coefficients (R) from the linear regression between eddy covariance ET and MODIS NDVI, EVI and LST using 16-day and annual data at six FLUXNET sites in PA systems (perennials and annuals vegetation systems, i.e. forests, woodlands, savannah and shrublands).Statistically significant correlations at p < 0.05 are indicated by * while * * indicates p = 0.06 and * * * indicates p = 0.07.

Table 4 .
Correlations coefficients (R) of three empirical VI-based ET models using MODIS-derived NDVI, EVI and LST.Results are for models using 16-day/annual data in AN (annual vegetation systems i.e. croplands and grasslands), and PA (perennials and annuals vegetation systems i.e. forests, savannah and shrublands) systems.All R were significant at p < 0.05 except for the 16-day ET-LST simple regression in PA.Mean annual NDVI and EVI were regressed against annual ET using linear and exponential functions.

Table 5 .
The mean absolute error (MAE) for Table4.The 16-day MAE is in millimetres per day (mm d −1 ), while annual MAE is in millimetres per year (mm yr −1 ).In parenthesis is the mean relative error (MAE/mean) in percent (%).