Integrated Methodology for Estimating Water Use in Mediterranean Agricultural Areas

Agricultural use is by far the largest consumer of fresh water worldwide, especially in the Mediterranean, where it has reached unsustainable levels, thus posing a serious threat to water resources. Having a good estimate of the water used in an agricultural area would help water managers create incentives for water savings at the farmer and basin level, and meet the demands of the European Water Framework Directive. This work presents an integrated methodology for estimating water use in Mediterranean agricultural areas. It is based on well established methods of estimating the actual evapotranspiration through surface energy fluxes, customized for better performance under the Mediterranean conditions: small parcel sizes, detailed crop pattern, and lack of necessary data. The methodology has been tested and validated on the agricultural plain of the river Strimonas (Greece) using a time series of Terra MODIS and Landsat 5 TM satellite images, and used to produce a seasonal water use map at a high spatial resolution. Finally, a tool has been designed to implement the methodology with a user-friendly interface, in order to facilitate its operational use.


Introduction
The world is facing a freshwater crisis and irrigated agriculture creates a high pressure as it is the single largest freshwater consumer on a global scale [1,2].Although driven by the need for increased production, the freshwater exploitation is frequently unsustainable, leading to depletion and deterioration of water resources [3].To alleviate these threats, the European Commission [4] has initiated the Water Framework Directive (WFD), with the general objectives of improving water quality, guaranteeing provision of sufficient quantities of good quality water for all purposes, and avoiding any long-term deterioration.Among others, the WFD demands identification of pressures from agriculture and estimation of water abstractions, with the overall aim of recovering the costs of water services, including environmental and resource costs.The WFD dictates that monitoring should be performed at the basin level.Although scale specifications are not defined, detailed mapping is required in the fragmented Mediterranean agricultural environment, thereby increasing the cost of monitoring [5].
Despite exceeding 80% of the total use in the Mediterranean, agricultural water use is difficult to estimate in large irrigation systems, especially due to the presence of private pumps with no water meters.Numerous methods have been developed to estimate this important, yet often unknown, part of the water cycle.Ground based measurement methods employ specialized instruments such as lysimeters, scintillometers, and eddy covariance meters to estimate, directly or indirectly, the water loss through evapotranspiration [6,7].They provide an accurate estimation of evapotranspiration at a specific location, but because of their high cost they are limited for use in local validation of other methods.The water balance methods use hydrological models, such as Soil Water Atmosphere Plant -SWAP [8], Semi-distributed Land Use-based Runoff Processes-SLURP [9], and Simulation of the Water Balance of a Cropped Soil -SWBACROS [9][10][11], to estimate the components of the water cycle.However, the substantial data requirements for such methods can be prohibitive in ungauged basins.An important milestone in the category of meteorological based methods is the publication of FAO 56 method [12], which has since been considered the standard method for estimating crop evapotranspiration and subsequently crop water use.However, the spatial component of the FAO 56 method was originally missing and thus the crop-coefficient remote sensing methods have combined FAO 56 with satellite earth observation [13,14].Despite their simplicity, the crop-coefficient remote sensing methods may overestimate evapotranspiration at the basin level, as they frequently assume optimal and homogeneous growing of each crop, because it is not feasible to account for the specific conditions (growth stage, nutrient deficiencies, stresses, etc.) at the field level.
Extensive overviews of remote sensing methods for estimating evapotranspiration have been published [15][16][17], within which the energy balance methods stand out because of their advantages: the Surface Energy Balance Index-SEBI [18], Surface Energy Balance System-SEBS [17] and Surface Energy Balance Algorithms for Land-SEBAL [19] take into account the crop growing conditions and do not require a satellite-based crop classification, while the Two-Source Model-TSM [20], Disaggregated Atmosphere Land Exchange Inverse model-DisALEXI [21] and Dual-Temperature Difference-DTD [22] do not require local data or the identification of extreme conditions.These algorithms have been designed to calculate the energy balance at the regional scale, however, local parameterization can improve the accuracy of the model [23,24].Customizations of the above mentioned methods offering solutions to specific conditions of Mediterranean irrigated basins have been suggested in western Turkey [25], France [26][27][28] and Spain [29].But an integrated top-to-bottom remote sensing approach is missing.The aim of this work is to present an integrated methodology for the improved estimation of agricultural water use in Greek basins using satellite earth observation.The specific objectives are: (i) to improve existing methods for operation in Greek river basins and integrate them into an operational methodology, (ii) evaluate the methodology during its application in an irrigated river basin, and (iii) design and implement a user-friendly tool for further use of the methodology by the water managers.

Description of Study Area
The study area is the Greek part of the Strimonas river basin (Figure 1).The total basin of the transboundary river covers an area of 17,250 km 2 , with its upstream areas situated in Bulgaria (65% of total area).A large reservoir (artificial Lake Kerkini) was constructed in the 1930s immediately upstream of the Strimonas plain in Greece, to store water for irrigation and to prevent flood events.With progressive drainage networks, agricultural land has been reclaimed from the wetlands and temporary lakes that had dominated the flood plain of Strimonas.By the 1960's, an extended irrigation system (the largest in Greece) had been developed in the plain, which taps water from the reservoir and river diversion dams.Water is distributed mainly by open canals, which are predominantly lined but often badly maintained, thus offering low conveyance efficiencies.Local groundwater pumping stations provide supplementary water to account for the inefficient operation of the irrigation system.The large quantities of drainage water are reused for irrigating downstream regions.The total irrigated area is estimated in excess of 800 km 2 , including the areas irrigated by privately operated pump wells (upper reaches) [30].The main agricultural crops are irrigated maize, cotton, rice, beetroot, alfalfa, and rainfed winter cereals.Agricultural and livestock production are the main economic activities in the site.Water is managed by 14 Local Water Organizations, which are divided according to canal command areas.They undertake the tasks of water distribution, fees collection from farmers based on their irrigated area, and canal maintenance.Several environmental problems are reported to be connected to the unwise management of irrigation water.Due to Lake Kerkini's operation as an irrigation reservoir, its water level fluctuates by 5 m during the irrigation season (May to September) and its surface decreases from 75 km 2 in spring to 50 km 2 in early autumn.As a result, the riparian forest and other key habitats of this wetland of international importance have decreased drastically.Also the river's estuary suffers from low water quality due to agricultural drainage and low discharge during the irrigation season, which is not sufficient to sustain the coastal ecosystem.Finally, salinization and sodification have been reported in areas where upstream drainage water is the main source of irrigation or where the water table is high.

Data Description
The changing spatial pattern of irrigated crops during the irrigation season in the study area (from beginning of May to end of September) was provided by frequent, free of charge acquisitions from the Terra/MODIS sensor.Three Terra/MODIS products were required for SEBAL implementation: Surface Reflectance (MOD09Q1 at 250 m and MOD09A1 at 500 m resolution) and Land Surface Temperature (LST) and Emissivity products (MOD11A2 at 1 km resolution).All products have a level 3 of processing and provide 8-day composite images.MODIS images are acquired daily and the data quality may be affected by the cloud cover and variable viewing angle, thus decreasing the quality of the outputs.For this reason, low quality pixels (e.g., cloud affected) are substituted with other unaffected pixels of the same location but from another date within the 8-day period.The day of acquisition of each pixel is recorded in the product files.Altogether, 23 layers were downloaded for each of the 14 dates necessary to cover the irrigation season using the Earth Observing System (EOS) data gateway of the Land Processes Distributed Active Archive Center [31].
Daily meteorological data were required for water use estimation.These were extreme air temperature (Tmin and Tmax) and wind speed and relative humidity at the satellite overpass time acquired from three automatic meteorological stations located in the irrigated area.Instantaneous air temperatures (Ta) were also collected from the same stations in order to derive the LST to Ta relation and thus produce a Ta map of the area.Daily ETo values were estimated using the FAO 56 method [12] to describe the fluctuation of evapotranspiration from day to day.Due to the close proximity of the stations it was not possible to perform spatial interpolation, thus the readings were averaged and meteorological conditions were considered uniform across the study area.Data from these meteorological stations were available only for the period starting from June 2007, due to late installation of the systems.
One Landsat 5 TM image (acquired on 19/07/2007) was required to cover the spatial extent of the study area.The Landsat TM sensor provides data at high multispectral content (blue, green, red, near infrared, middle infrared, thermal infrared and short wave infrared bands), at a spatial resolution of 30 m (120 m for the thermal band).Being based on archived data collected every 16 days, it was possible to select an image strategically acquired at the peak of the irrigation season, therefore representative of the spatial distribution of the highest water consumption.

Integrated Methodology Customized for Greek Conditions
An integrated methodology for the production of a high-resolution water use map of Greek river basin is presented, incorporating existing methods, after customization to the local conditions of Greece.The methodology consists of the following components.

Estimating surface energy fluxes
The Surface Energy Balance Algorithm for Land (SEBAL) model aims at the empirical determination of the actual evapotranspiration (ETa) by solving the Energy Balance equation: The latent heat flux (λE) once divided by the latent heat of vaporization (λ) provides an evapotranspiration approximation.
In SEBAL, the net radiation R n is computed from satellite-measured broad-band reflectance and surface temperature, while the soil heat flux G 0 is estimated from R n , the surface temperature, and vegetation indices.The sensible heat flux H, is considered proportional to the ratio between the surface-air temperature difference (dT) and the bulk aerodynamic resistance (r ah ) [19].
SEBAL uses the partitioning of H and λE as described in [19], where the evaporative fraction (Λ) is calculated as: The underlying assumption is that Λ is constant over the day, or said differently: the instantaneous partition of λE and H is equal to the average diurnal partitioning ratio.The difference between Λ at the moment of satellite overpass and Λ derived from the 24-hour integrated energy balance has been considered as non significant [32].
The latent heat flux was calculated based on Λ, multiplied by the diurnal net radiation at earth surface (R nday ) that is itself retrieved spatially from the surface broadband Αlbedo (α) and the singleway transmissivity of the atmosphere (τ sw ).
The implementation of SEBAL requires multi-spectral remote sensing data (blue, green, red, near infrared, shortwave infrared, and thermal infrared bands), meteorological data (instantaneous wind speed, instantaneous relative humidity, and extreme air temperatures), meteorological station characteristics (height of wind measurement and vegetation height around met station), and topological data (elevation).Although SEBAL has been developed for instantaneous data at the time of image overpass, 8-day composite images can be used as long as the date and time of each pixel's acquisition is incorporated into the algorithm, and SEBAL is applied at the pixel level.
SEBAL provides a good estimation of λE in flat agricultural areas.It has been tested extensively under several climatic conditions at both field and regional scales.High accuracy of results has been reported [33].Regarding Mediterranean areas, SEBAL has been validated in the South-East of France [28] and in the Gediz river basin in western Turkey [25].Several customizations are proposed hereafter to provide optimum performance of SEBAL in the local conditions.

Customization of Albedo (α) estimation
Albedo, which is input to SEBAL and represents the shortwave surface reflectance in the range of 0.3-3 μm, was calculated from band-wise surface reflectance (SR).Then, a specific histogram stretching method was designed for each type of input imagery conditions and characteristics.The histogram stretching method was different in calibration for Landsat 5 TM and Terra MODIS, since spatial resolution, atmospheric and geographical conditions are different from one to another.For Terra MODIS, although surface reflectance products were already available by the MODIS Science Team [31], it was found useful to standardize their response inside image time series acquired under different atmospheric conditions/corrections performed by the automated algorithm [34].Available Landsat 5 TM data [35] were top-of-atmosphere readings without any atmospheric correction.We used a specific histogram stretch to produce surface albedo without an atmospheric correction model.
Considering Figure 2, the algorithm first identifies the three highest points (Top: 1, 2, 3) and the six bottom points (bottom: 1a, 1b, 2a, 2b, 3a, 3b).It then stretches the Landsat 5 TM histogram in order to achieve an output Top1 albedo value to 0.05 and Bottom2b albedo value to 0.36, thus forcing the median of water histogram to become 0.05 and the sand median to be 0.36, both values being helpful to SEBAL iteration stability response.For MODIS, the algorithm sets Bottom1a albedo value to 0.05 and Bottom2b albedo value to 0.36, thus stretching with a linear function the whole histogram, thus forcing the water histogram to values higher than or equal to those that correspond to the expected range of water albedo (>0.05).This configuration is necessary for correcting surface albedo of the two input datasets, taking into account their pre-processing steps.

Customization of air temperature (Ta) estimation
Air temperature is required at various processing stages for SEBAL implementation (e.g., R n , H), either as Ta or dT input, dT being the difference between the surface temperature (Ts) and Ta.R n is estimated from the incoming longwave radiation L in that is a function of Ta: Therefore, the accuracy of the Ta estimation is essential for an optimum accuracy in solving the energy balance equation.However, for the H term: An initial value of dT is first inserted (use of the radiometric surface temperature as a surrogate for the aerodynamic surface temperature To), which is further adjusted with a set of iterations of estimating H.
A linear relation between dT and Ts has been proposed [36], while others [37] used a step function to retrieve Ta from Ts values.In this study, a regression analysis was performed using the Ta values available from the three meteorological stations and the radiometric LST provided by the MODIS product.The resulting linear equation for the whole irrigation season is the following: According to Figure 3, the Pawan [36] equation overestimates the LST to Ta temperature difference, while the equation proposed by Bastiaanssen et al. [19] provides closer estimations of the air temperatures.Since ground data are available, the use of the empirical linear relation is advantageous to allow modelling of the local weather conditions in a more realistic way.
As a result, the pixel based assessment of Ta with the SEBAL method of linear relationship between dT and LST provides a Ta map with information on the energy balance variations that point-based air temperature from meteorological stations fail to supply.

Considerations in extreme pixel selection
A fundamental advantage of SEBAL is the calibration of the result using pixels of extreme meteorological conditions: a very wet pixel with negligible H (H wet ~ 0) and a very dry pixel with negligible λE (λE dry ~ 0).Under the Mediterranean semi-arid summer conditions, the wet pixel during the irrigation season is usually a very well irrigated homogeneous patch (rice fields or other water demanding crop).Thus, the wet pixel was selected following the criteria: low temperature, low albedo, and high NDVI.For the dry pixel, it has been suggested to select a patch of desert, burned area, or the outskirts of a densely populated settlement [16].In the absence of desert or burned areas during the irrigated season in Greece, the dry pixel was searched for in poorly-vegetated areas or harvested winter cereals (low NDVI) presenting a high temperature value and a high albedo.

Integrating over the temporal dimension
The SEBAL and other suggested methods for ETa estimation on certain sample dates during the cropping season give a very good indication of its spatial distribution, but cannot be used directly in water balance studies.The reason behind this is the wide fluctuation of ETa from day to day, depending on meteorological conditions, availability of water, and plant growth (phenological stage).Therefore, in order to obtain an accurate estimation of the seasonal ETa, daily values have to be simulated.
In recent literature [38,39] a method for temporal integration was suggested, in which the fraction ETrF j = ETa j /ETo j (where j refers to the satellite image acquisition date and ETo j to the potential evapotranspiration for day j) was considered constant for the time period between two consecutive satellite images, and the seasonal ETa (termed ETs) could be simulated with the following equation: where t j and t j+1 delimit a short time period around the acquisition date j.
The above mentioned equation integrates the time component (ETo, estimated with the FAO 56 method [12]) and the spatial component (ETrF) to describe successfully the daily fluctuation of ETa across the study area.ETa also depends on the availability of water.Spatial variation of water availability can be described by the [t j , t j+1 ] time interval of the satellite images, which were assumed representative for each period.
This temporal integration method was applied to the MODIS data.Since input images were composites, the day of acquisition could differ from one pixel to another, thus the integration was performed at the pixel level.Despite the use of composite images, data could carry some null values due to cloudy conditions or low-quality data [40].The appearance of these null values even in a single composite product could influence the result in the whole irrigation season; therefore a dynamic selection of the integration period was designed.The start (resp.end) of the integration period for a certain pixel was designed to be in the middle of the period between this pixel and the temporally preceding (resp.following) pixel with non-null value in the same location.This dynamic selection is based on the assumption that the pixel acquired in the preceding or the following composite is the closest match for a missing pixel value.

Improving spatial resolution
The spatial resolution of ETs described above is limited by the resolution of the sensor (e.g., MODIS: 250-1000 m at nadir).Considering the small parcel sizes in Greece, the detailed crop pattern, and the interspersion of irrigated with non-irrigated fields on the fringe of the irrigated area, it is suggested that a high resolution ETs map is required to adequately describe the pattern of irrigation water use.In such cases, downscaling for higher resolution results can be achieved by merging high resolution data of similar content.
Merging high and low resolution images has been achieved at the pixel, feature or decision level [41].The resolution (or scale) of several hydrologic parameters has been improved in recent literature, such as precipitation from the TRMM sensor [42], soil moisture from MODIS data [43], and thermal infrared images [44].In similar studies [39,45] ETs estimated from NOAA AVHRR images (at 1.1 km resolution) were merged with a single day ETa estimated from a Landsat ETM+ image (30-60 m resolution) using the modified Brovey method.
The proposed methodology incorporated the modified Brovey method [39,45] to improve the resolution of the ETs map using higher resolution maps of similar content, such as ETa-Landsat at 30 m.According to this method, the detailed pattern of the ETa map was used as a weight to spatially redistribute the ETs calculated for each 1.21 km 2 pixel of the NOAA AVHRR satellite using the formula: Where is the average of the values of ET ETM+ (Landsat) that occupy a pixel of ETs (NOAA).Using this proportional distribution method, the total ETs value calculated for each area of 1.21 km 2 was preserved, but redistributed according to the detailed pattern of ETa derived from the ETM+ image.The advantage of the method is that it minimizes the influence of the high-resolution image induced in the methods found in literature (IHS and PC transformation).
The detailed spatial content of these satellite images is sufficient to cover the specific conditions of Greek irrigated river basins.However, care should be given to the selection of their acquisition date.A single date image may be selected, provided that it is representative of the spatial distribution of the seasonal water use, otherwise, if spatial distribution of the irrigation and crop patterns varies over time, several images scattered within the irrigation season should be used.

Development of the ITA-Water Tool
A stand alone tool, named ITA-Water tool (Integrated Thermodynamic Algorithms for Estimating Agricultural Water Use), has been developed in order to facilitate the implementation of the proposed methodology.The ITA-Water tool presents four separate parts: -Wizard 1 (resp.Wizard 1 L5): Generation of intermediate input maps (NDVI, albedo, Ts) for Terra/MODIS (resp.Landsat 5 TM) input data.
-Wizard 3: Generation of the ETs-MODIS low-resolution map through temporal integration of ETa-MODIS and daily ETo data.
-Wizard 4: Generation of ETs-Landsat map through resolution improvement of ETs-MODIS with ETa-Landsat maps.
Each wizard offers a user-friendly interface that is accessible through an initial launching wizard.The downloading and pre-processing of MODIS, Landsat 5, and meteorological data is a pre-requisite for the use of the ITA-Water tool.In order to generate a high-resolution ETs map from the various input data the operation scheme presented in Figure 4 should be followed.For wider access and increased flexibility, the same functionalities have been made available in GRASS GIS v6.3 or v6.4 as add-ons [46], in GRASS GIS v7 as generic modules [47] and in R statistical language through the "RemoteSensing" plug-in package [48].The availability of these free-of charge tools guarantees their wide use by interested researchers or water managers.

Experimental Results in Strimonas River Basin
The application of the integrated methodology in Strimonas river basin using the ITA-Water tool has provided the high resolution ETs-Landsat map (Figure 5).It has to be noted that results have not been provided for May 2007 because of missing meteorological data.Using standard overlay geographic analysis, mean ETs per irrigation district was estimated (labels on Figure 5).The total water consumption for irrigation over the period June to September 2007 in the Strimonas river basin is estimated at 365.49 × 10 6 m 3 .The spatial distribution reveals a pattern of higher water use closer to the lake and the river.The main branches of the conveyance network are fed by these water bodies, which being unlined, are responsible for significant water losses.Another reason for the high levels of ETs in this area is the cultivation of paddy rice, which not only requires high water input, but also synchronized irrigation applications, which is translated to low distribution efficiency of the irrigation network [49].Despite the improper agricultural water management in the area, water deficiencies have not been reported due to the low water utilization in the upstream Bulgarian territory.However, this could change as water development may be implemented there in the near future.Time series at three pixels located in rice, corn and wheat fields have been generated in order to visualize the temporal variation of ETa values in the dominant crop types in the study area (Figure 6).The curves' trend is coherent with the expected temporal evolution of evapotranspiration in the Mediterranean region, taking into account the seasonality of ETo and the rare rain events in the semiarid summer, which are minimal as compared to irrigation applications.Evapotranspiration is increasing while progressing into the hot and dry summer months, and peaks in mid July (doy 190-200).The existence of the two abrupt changes for rice and corn in mid June (doy 169) and mid September (doy 249) may be related to extreme meteorological conditions or low quality satellite data.For the rain-fed wheat pixel, the decrease of ETa in June (from doy 160) is in accordance with the crop's maturation and harvest season.During the following months, the wheat's low fluctuation of 1-2 mm/day represents ETa from weeds and other natural vegetation which nourish on soil moisture replenished by short summer storms.

Quality assessment breakpoints
Various Quality Assessment (QA) breakpoints were set throughout the processing chain.Two of the QA breakpoints concerned the albedo and the R n to G 0 ratio.For these two terms, a few sample pixels were selected for various cover types (e.g., corn field, clay, water, etc) and the value obtained at the pixel location was checked against the typical range proposed in Waters et al. [38].The values of albedo and R n to G 0 ratio at the selected pixels were within the reported range, or else the source of error was searched in the input parameters.The third QA breakpoint was related to the comparison of SEBAL ETo at the meteorological station location against the values obtained using the FAO 56 Penman-Monteith method [12].The error in ETo estimation was 14.7% (DF = 41, t = 6.252,P < 0.0001).

Validation of ETa-MODIS maps using ETc meteorological data
The validation of ETa-MODIS maps was performed using the crop evapotranspiration (ETc) values at the meteorological stations.The ETc values were calculated using the FAO 56 method [12] as the product of the ETo and the crop coefficients (K c ).The K c accounts for elements that characterize each crop's evapotranspiration from the reference crop's evapotranspiration (grass): the crop height, the albedo (reflectance) of the crop-soil surface, the canopy resistance, and the evaporation from soil.Therefore, the FAO 56 is a standard approach for the determination of K c according to the crop type, the climate, the soil evaporation, and the crop growth stage.The main advantage of this method is that it is considered as standard, thus it has been widely used leading to a subsequent availability of support material and results for comparison purposes.
Since corn (resp.cotton) fields surround the meteorological station(s) A9 (resp.A11 and A12), the K c values relative to these two crop types were calculated in order to derive the ETc values at the meteorological stations.The average difference between the ETc values and the ETa estimations over the whole irrigation season at the three stations was 19.3% (DF = 36, t = 0.949, P = 0.348).This rate is satisfactory considering that the K c coefficient method does not account for water stress factors, meaning that ETc represents an accurate estimation of ET only under good irrigation conditions.Also, ETa values were extracted from low spatial resolution maps (250 m), where the crop distribution over such large area may not be homogeneous.

Validation of ETs-Landsat map using reference water supply data
The final resulting map (ETs-Landsat) was compared against water supply data, reported by the Directorate of Water Resources of Serres, which oversees the Local Water Organizations in the floodplain of Strimonas basin [49].These have been collected daily by noting the gates opening at the head of each main canal, and calibrated regularly using a propeller current meter.The water released in each canal command area was converted to monthly evapotranspiration using the general irrigation efficiencies reported for each irrigation canal type.ETs' from reference datasets can be estimated using the formula: ETs' = a × WS, where a is the efficiency of the irrigation system and WS is the supplied water.This reference dataset was compared with the area averaged ETs per month and per canal command area (Figure 7).A slight tendency for underestimating ETs with the proposed methodology was noted in the high end values.However, the mean difference of -0.25Mm 3  was not statistically significant (paired T-test, t = -0.57,DF = 51, P = 0.5712).

Validation of the spatial resolution improvement method
The resolution improved ETs using a Landsat TM image (ETs-Landsat at 30m) was compared against the resolution improved ETs using an existing map of the Strimonas river basin termed ETs-ALOS (10m) [50].The ETs-ALOS map was derived from ALOS AVNIR-2 data by applying the K c method.According to this method, the ALOS AVNIR-2 image was subjected to a maximum likelihood spectral classification, using field surveyed training data in order to produce a crop distribution map with the required spatial distribution of K c .The classes of the crop map were assigned to ETs for each crop, which had been estimated using the FAO 56 method in a worksheet, taking into account meteorological data and crop type [50].
The resolution improvement of the ETs-MODIS map using the detailed spatial pattern provided by: (a) ETa-Landsat or (b) ETs-ALOS is displayed in Figure 8. Visual assessment did not reveal any differences, and the two products were compared (i) on a pixel level, using a random dataset of 1,000 points, and (ii) on an irrigation district area level, using a polygon vector dataset.Statistical comparison was performed with paired T-test and correlation analysis.Pixel level comparison revealed that the mean difference between the two sampled datasets was not statistically significant (mean difference of 2.34 mm, DF = 849, t = 0.8396, P = 0.4013), and displayed a strong positive correlation (r = 0.86).Area based comparison of the canal command area level revealed that the mean difference between the two datasets was not statistically significant (mean difference of 77 Mm 3 , DF = 10, t = 0.3268, P = 0.7505), and displayed a very strong positive correlation (r = 0.99).This proves that resolution improvement with a single Landsat 5 TM image was successful in describing the spatial pattern of water use during the irrigation season.

Advantages of the methodology
The advantages of the proposed methodology are: (i) the use of MODIS remotely sensed data that cover the study area every eight days guarantees that the factors which affect the spatial distribution of energy balance (crop phenology, plants stress, water availability) are updated frequently enough for the changing Mediterranean pattern [51], (ii) the use of daily meteorological data guarantees that the factors which affect the temporal variation of evapotranspiration (temperature, relative humidity, wind speed) are described frequently enough [12], (iii) the resolution improvement guarantees that the detailed crops pattern of Greek agricultural landscapes is adequately described to avoid the "oasis effect" which appears when irrigated fields are interspersed within non irrigated fields [52], (iv) using a single Landsat 5 TM image is sufficient to describe the spatial pattern of seasonal water use, thus minimizing the need to process multiple high resolution datasets which is time consuming, and (v) the cost of data acquisition is minimal by using free of charge (Terra MODIS, Landsat 5 TM) or low cost satellite images (ALOS AVNIR-2).
Based on the validation with reference datasets, the customized SEBAL approach proved to perform well with multi-spatial resolution input data (Terra MODIS at 250 m, 500 m and 1 km; Landsat 5 TM at 30 m and 120 m).The spatio-temporal integration using ETo-driven heuristics is a major asset for overcoming the unavailability of spatially distributed data on a daily basis.In addition, the use of multi-source data for ETs resolution improvement has been successfully validated.Consequently, the proposed methodology can be positively extended to use alternative combinations of data for water use mapping (TIR: AATSR, AVHRR, METEOSAT, ASTER; VNIR/SWIR: MERIS, ASTER, ALOS, SPOT, AVHRR, METEOSAT, and other upcoming missions).The fact that MODIS data are provided free of charge and that all Landsat imagery is freely available for download (since April 2009) were major motives for the selection of this specific satellite data combination.However, the use of alternative data may be dictated by specific needs such as higher-resolution, more frequent data availability or ad-hoc geographical coverage.

Regionalization and implications for water management
The Mediterranean region is characterized by small to medium sized river basins with conditions similar to the study area: arid summers, mild topographic relief, small agricultural field sizes or crop patterns with high spatial distribution, and lack of detailed datasets.Most of them face similar problems related to the unwise management of water resources: depleted aquifers [53], increased concentration of nitrogen and other pollutants [54,55], saltwater intrusion [56,57], and degradation of nearby natural environment [58][59][60].The changing climate is only expected to exacerbate the problems, adding the risk of desertification in these areas [61].Monitoring these areas is difficult because of the scattered dissimilar environment and the lack of adequate facilities, even a network of meteorological stations.
The proposed methodology can be applied in these Mediterranean regions providing considerable benefits.Through the user-friendly interface (ΙΤΑ-Water), and using low cost data, water managers and decision makers can produce actual water consumption thematic maps over the entire irrigated basin.Thus, an often unknown part of the agricultural water cycle may be accurately estimated and reporting to cover national and regional legislative requirements may be facilitated.Also, in having this information, several other secondary information may be calculated, such as irrigation system performance indicators, water valuation parameters, design of water savings measures, and desertification risk assessment, all of which allow for added value of water services.In addition, it can be applied to previous years to identify temporal changes using archived satellite images, even when field surveys had not been conducted.However, more sophisticated methodologies may be required for monitoring water use of large basins expanding in several agro-climatic zones, due to the inconsistent crop cycles [62,63].

Conclusions
An integrated methodology has been presented for improved estimation of agricultural water use in Greek basins using satellite earth observation.Several improvements on existing methods have been implemented to achieve a high performance under Greek conditions: higher spatial accuracy of surface energy fluxes estimation, consideration of variable Mediterranean agricultural pattern, and tailored to the low availability of datasets.All together, these methods have been embedded in an integrated methodology, which was rendered operational because it: accepts input from multiple data sources, offers options for alternative methods, and is implemented with a user friendly tool for further use of the methodology by the water managers.The proposed methodology has been tested in an irrigated river basin of 800 km 2 .The validation shows high accuracy, when compared to alternative methods and available reference datasets.

Figure 1 .
Figure 1.Location of the study area.

Figure 3 .
Figure 3. Three alternative equations relating surface temperature (Ts) to the air to surface temperature difference (dT).

Figure 4 .
Figure 4. Operation scheme of the ITA-Water tool.

Figure 5 .
Figure 5. Agricultural water use map (ETs-Landsat) from June to September 2007 in the agricultural area of Strimonas river basin.Labels indicate mean seasonal water use (mm) per irrigation district.

Figure 6 .
Figure 6.Temporal variation of ETa for irrigated rice and corn, and rain-fed wheat pixels from June to September 2007.

Figure 7 .
Figure 7.Comparison of ETs Landsat with reference data in Strimonas river basin, per month and per canal command area.