Monitoring Irrigation Consumption Using High Resolution NDVI Image Time Series: Calibration and Validation in the Kairouan Plain (Tunisia)

Water scarcity is one of the main factors limiting agricultural development in semi-arid areas. Remote sensing has long been used as an input for crop water balance monitoring. The increasing availability of high resolution high repetitivity remote sensing (forthcoming Sentinel-2 mission) offers an unprecedented opportunity to improve this monitoring. In this study, regional crop water consumption was estimated with the SAMIR software (SAtellite Monitoring of IRrigation) using the FAO-56 dual crop coefficient water balance model fed with high resolution NDVI image time series providing estimates of both the actual basal crop coefficient and the vegetation fraction cover. Three time series of SPOT5 images have been acquired over an irrigated area in central Tunisia along with a SPOT4 time series acquired in the frame of the SPOT4-Take5 experiment, which occurred during the first half of 2013. Using invariant objects located in the scene, normalization of the SPOT5 time series was realized based on the SPOT4-Take5 time series. Hence, a NDVI time profile was generated for each pixel. The operationality and accuracy of the SAMIR tool was assessed at both plot scale (calibration based on evapotranspiration ground measurements) and perimeter scale (irrigation volumes) when several land use types, OPEN ACCESS Remote Sens. 2015, 7 13006 irrigation and agricultural practices are intertwined in a given landscape. Results at plot scale gave after calibration an average Nash efficiency of 0.57 between observed and modeled evapotranspiration for two plots (barley and wheat). When aggregated for the whole season, modeled irrigation volumes at perimeter scale for all campaigns were close to observed ones (resp. 135 and 121 mm, overestimation of 11.5%). However, spatialized evapotranspiration and irrigation volumes need to be improved at finer timescales.


Introduction
In arid and semi-arid regions, water availability is a major limitation for crop production. In the Kairouan plain (Central Tunisia), the combined effect of drought spells and the increase of irrigated surfaces during the last decades have had a negative impact on the available water resources. Efficient agricultural water management is therefore a major issue, especially in irrigated areas. The design of tools that provide regional estimates of the water balance may help the sustainable management of water resources in these regions.
Evapotranspiration (ET) is one of the most important fluxes of the water balance in semi-arid areas; it is a key factor for optimizing irrigation water management [1]. Direct measurements of ET are only possible at local scale (single plot) using for example eddy-covariance devices. Scintillometers measure sensible heat flux along a given path and then latent heat flux (ET) is returned as a residual term of the surface energy budget. Furthermore, remote sensing (RS) capabilities for monitoring vegetation and its physical properties on large areas have been identified for years now [2]. It provides spatialized and periodic information about some major drivers of ET such as albedo, surface temperature and vegetation properties. Several methods for estimating ET using remotely-sensed data have been developed [3][4][5][6][7]. Most of them solve the surface energy budget for latent heat using thermal imagery. Instantaneous estimates at the time of satellite overpass have been successfully used to estimate ET at daily scale [8]. However, the main limitations of these methods are the difficulties in obtaining valid estimates of the aerodynamic surface temperature and the atmospheric resistance to heat transfer in single-pixel methods or the difficulty in identifying the wet and dry edges when using triangle-based methods [9]. A further limitation arises when trying to extrapolate ET beyond one day due to the limited availability of high resolution thermal imagery.
Another possible approach is to use Soil Vegetation Atmosphere Transfer models (SVAT) to simulate ET. These models can benefit from remote sensing since the latter provides periodic information about the vegetation development which is a primary factor driving evapotranspiration. The use of high-resolution image time series for monitoring irrigated crops was more recently discussed [10][11][12]. The low availability of such data, for financial as well as technical reasons, combined with the intermittent presence of cloud, has been a restraint to their use [13][14][15]. However, the forthcoming Sentinel-2 mission offers a unique opportunity to improve this monitoring thanks to high resolution (10 m) and high repetitivity (5 days) visible and near infrared (VIS-NIR) remote sensing.
For the operational monitoring of soil-plant water balance, the most common and practical approach used for estimating crop water requirement is the FAO-56 method [16]. The FAO 56 dual crop coefficient approach uses two coefficients to separate the respective contribution of plant transpiration (Kcb) and soil evaporation (Ke). However, standard basal crop coefficients (Kcb) profiles provided by FAO tables are average values not suited for specific growth conditions that can largely differ between plots. Remote sensing is a valuable asset to derive those temporal profiles of crop coefficients. It has been shown that the crop coefficients were linked to the spectral response of the cover, especially vegetation indices [12,[17][18][19][20]. Ke is linked to the bare soil fraction, complementary of the fractional vegetation cover (fc) which can also be related to visible RS data [21]. Although the relations proposed between Kcb, fc and vegetation indices are not theoretically fully linear, they can usually be approximated by linear relations [22][23]. Moreover, establishing a unique relationship between crop coefficient and spectral vegetation indices is an ongoing research topic [24] and many empirical linear relationships available in the literature have been derived experimentally.
The FAO-56 method has long been used to monitor water budget at plot scale with tools like CROPWAT [25]. The interest for coupling the FAO-56 method with remotely sensed crop coefficients is rising alongside the increasing availability of high resolution Normalized Difference Vegetation Index (NDVI) time series [26][27][28][29][30]. The SAMIR tool (SAtellite Monitoring of IRrigation) [31] used in this paper computes spatially distributed estimates of ET and crop water budget at regional scale. It is based on the coupling of the FAO-56 dual crop coefficient model with time series of high resolution NDVI imagery (Normalized Difference Vegetation Index) providing estimates of the actual basal crop coefficient (Kcb) and the vegetation fraction cover (fc).
In this study, regional evapotranspiration and crop water consumption were estimated over an irrigated area located in the Kairouan plain using the SAMIR model fed by SPOT high resolution time series. The model was calibrated on the basis of local ET measurements from flux towers and was validated at perimeter scale using known irrigation volumes. The objective of the work was to assess the operationality and accuracy of SAMIR outputs at plot and perimeter scales, in a context of high land cover complexity (i.e., trees, winter cereals, summer vegetables) and limited data available for parameterization.

Study Area
The experimental site is located in the Kairouan plain, a semi-arid region in central Tunisia (9°30′E to 10°15′E, 35°N to 35°45′N) (Figure 1), covering an area of more than 3000 km 2 , which is part of the Merguellil watershed. The rainfall patterns are highly variable in time and space with an average annual rainfall of approximately 300 mm (extreme values recorded in Kairouan city are 108 mm in 1950/51 and 703 mm in 1969/70). The mean daily temperature in the city of Kairouan is 19.2 °C (minimum of 10.7 °C in January and maximum of 28.6 °C in August). The relative humidity ranges between 70% and 55% in winter and 40% and 55% in summer. The mean annual reference evapotranspiration estimated by the Penman-Monteith method is close to 1600 mm. Dominant crops in this region are cereals, olive and fruit trees and market gardening [32].
Water management in the Merguellil basin is characteristic of semi-arid regions with an upstream sub-basin that collects surface and subsurface flows to a dam (the El Haouareb dam), and a downstream plain supporting irrigated agriculture ( Figure 1). Irrigation water comes exclusively from the groundwater, except for a very small part of the plain on the edge of the dam: the major part of dam water infiltrates to the downstream aquifer. The main user of the Kairouan aquifer is agriculture, which consumes more than 80% of the total amount extracted each year [1]. Most farmers in the Kairouan plain extract water for irrigation directly from private wells, while a few rely on public irrigation schemes based on collective networks of water distribution pipelines stemming from a main gauged borehole. Each borehole corresponds to one organizational unit named GDA ("Groupement de Développement Agricole"). Annual consumption exceeds the annual recharge of the water table resulting in a piezometric decrease of between 0.5 m and 1 m per year [33].

Experimental Setup and Data Pre-Processing
Half hourly meteorological measurements were recorded using an automated weather station installed in the study area. It includes measurements of solar radiation, air temperature and humidity, wind speed and rainfall. Cumulative precipitation and reference evapotranspiration (ETo) values between November and June were respectively 238 mm and 1008 mm, for the 2008/2009 season, 142 mm and 666 mm for the 2011/2012 season and 97 mm and 992 mm for the 2012/2013 season.
A flux station was installed in a plot located in the study area. It measures the various energy balance components using the eddy correlation method. Energy fluxes measurements were acquired over irrigated barley and irrigated wheat during the 2011-2012 and 2012-2013 seasons, respectively. These experiments allowed continuous monitoring of actual ET as well as soil moisture measurements. Energy balance closure of the EC measurements was checked and corrected using the residual method. The low uncorrected value obtained at both sites (around 60% of closure) led to discard fast response psychrometer measurements. Few isolated inconsistent peaks were also removed. Overall, the quality of eddy covariance measurements is mainly affected by instrumental errors, uncorrected sensor configurations or problems of heterogeneities in the area and atmospheric conditions [34].
Time series of SPOT5 image acquisitions were planned over the plain: Nine, six and ten SPOT5 images were acquired for the 2008-2009, 2011-2012 and 2012-2013 seasons, respectively ( Figure 2).The SPOT5 images for the three campaigns were georeferenced using orthorectification and then radiometrically corrected to obtain top of canopy (TOC) reflectance on the basis of physical modeling corrections using the Simplified Method for Atmospheric Corrections (SMAC) algorithm based on the 6S radiative transfer model [35]. The SMAC 6S model was applied for each image using values of atmospheric optical depth and water content taken from a photometer located in the area and part of the AERONET network [36]. Due to the uncertainties in the atmospheric parameters, and in order to eliminate time profile artifacts due to radiometric correction discrepancies within the time series, an additional inter-calibration between dates was achieved based on the identification of pseudo-invariant features for which a constant reflectance value is assumed over time. Indeed, the radiative transfer model shows that for a flat topography and an assumed spatially homogeneous atmosphere, the reflectance can be linearly related to the image DNs [37]. Thus, an additional normalization of the image time series was achieved by linear correction of the identified inconsistent dates [30,38]. For the last season, we also benefited from images acquired in the frame of the SPOT4-Take5 experiment which occurred during the first half of 2013 [39] and whose main purpose was to simulate the revisit frequency and resolution of Sentinel-2 images to help users set up and test their applications and methods before the mission is launched. In this frame, SPOT4 images at 20 meters resolution were acquired every 5th day from 3 rd February to 18 th June 2013 over the Kairouan plain, but only 14 dates were cloud free among the 28 images acquired ( Figure 2). The SPOT4 series was corrected using the Multi-sensor Atmospheric Correction and Cloud Screening (MACCS) algorithm taking into account both temporal and spectral approaches for retrieving the aerosol optical thickness (AOT) [40]. Monthly irrigation volumes used for validation were obtained at the scale of each GDA irrigated sector. It was assumed that these data were trustworthy since these entities manage a collective well equipped with a meter, providing the water to the plots inside the perimeter. However, some plots outside the official perimeter also benefit from this water and in the frame of the acquisition of our validation data, they were delineated with the help of the irrigation manager. Conversely, no private well is exploited inside the GDAs, so that the monthly volumes collected can be reliably linked to the declared cultivated area.

Model Description
The algorithmic basis of SAMIR is the FAO dual crop coefficient method under stress conditions, i.e., considering actual soil water status [16]. This approach has been largely used for irrigation scheduling and to compute crop evapotranspiration for operational purpose. It is based on the concept of reference evapotranspiration for a standard well-watered grass, modulated by crop coefficients to account for the specific development of any vegetation cover as well as its actual water status. We present here only the major equations of the model and the modifications implemented within SAMIR. The reader should refer to the FAO paper N°56 for other equations; we kept for clarity the same notations as FAO paper 56 [16]. The actual (adjusted) evapotranspiration ETa of a crop is defined as: where Kcb is the basal crop coefficient representing unstressed crop transpiration, Ke is the evaporation coefficient representing soil evaporation and Ks is a stress coefficient accounting for the reduction of transpiration due to water shortage in the root zone. ET0 was computed following the FAO paper 56 and is thus not presented here. Regarding Ke, in order to account for frequent overestimations of bare soil evaporation as observed by [41], we used the same formalism as the latter reference to modify the Kr coefficient, i.e., the evaporation reduction coefficient accounting for water availability in the evaporation layer which is used to compute Ke.
with TEW the total evaporable water, REW the easily evaporable water and De the depletion (water deficit) in the shallow surface layer used to compute soil evaporation. The m coefficient lies between 0 and 1 and allows to further reduce the maximum evaporation level when REW = 0 and is functionally equivalent to a minimum surface resistance to evaporation of the soil.
As an extension to the standard FAO-56 formalism, an additional layer of depth Zd was considered beyond the root depth Zr to account for capillary flow from below the root zone [42]. The total depth of soil involved in crop functioning Zsoil is defined as Zsoil = Zr + Zd. To allow the water stored in this deep layer to be used by the plant, a diffusion process was introduced in the SAMIR model to simulate capillary flow between deep and root layers (Difrd). For coherence, diffusion was also introduced between the root zone and the shallow evaporation layer (Difer), allowing especially evaporation to last long after a wetting event, and also to the deeper layers to sustain low evaporation fluxes observed after harvest. The depth of the root and deep layers evolves dynamically with root growth: when root depth increases, a portion of the deep layer is included in the root zone. The total available water in the deep compartment (TDW) is computed similarly to total available water in the root layer (TAW) using the following formula: • with θfc the water content at field capacity and θwp the water content at wilting point. Water diffusion is driven by the water gradient between layers as follows: . .
where De, Dr and Dd are the depletions in the evaporation , root and deep layers, respectively, while cdE and cdR are the diffusion coefficients for the transfers between the root zone and the surface layer, and the root zone and the deep layer, respectively (mm·day −1 ). Regarding the vegetation development, instead of using standard values provided by the FAO paper 56, the major specificity of SAMIR is to use remote sensing data to estimate the actual basal crop coefficient Kcb using a linear relationship: where NDVI is the Normalized Difference Vegetation Index, depending on near infra red (NIR) and red (R) reflectances: The NDVI time series is linearly interpolated for each day between the successive dates of image acquisition.
In the same manner, the vegetation fraction cover fc was derived from NDVI using the following linear relation: • Then the root depth is linked to the vegetation fraction cover using the following formula: where fcmax is the maximum fraction cover for which the maximum rooting depth Zrmax is reached and Zrmin is the minimum rooting depth when the vegetation is detected by the satellite (fc > 0). Finally, the model updates the water content of the three soil layers at a daily time step in order to compute the water budget using equations similar to those in the FAO paper 56 [16]. The main difference lies in the added Difer and Difrd terms for water diffusion, and the addition of the deep layer. The surface runoff is neglected.
The water content in the evaporation layer is updated as follows: where Dej is the cumulative depth of evaporation (depletion) at the end of day j (mm), Dej − 1 is the same variable at the end of day j − 1 (mm), Pj is the precipitation, Ij is irrigation depth, Ej is the soil evaporation, Tewj is the depth of transpiration from the exposed and wetted fraction of the soil surface layer, all on day j (mm), fw is the fraction of soil surface wetted by irrigation and few is the wetted soil fraction exposed to evaporation. The water content in the root layer is updated as follows: where Drj is the depletion in the root layer at the end of day j (mm), Drj − 1 is the same variable on day j − 1, ETcj is the soil evaporation on day j (mm). If the water content goes beyond field capacity after a heavy rain or irrigation (i.e., Drj < 0), it is assumed that the amount of water above field capacity is lost the same day by percolation to the deep layer (DPrj = −Drj) and then Dr is set to zero. The depletion in the deep layer is computed as follows: where Ddj is the depletion in the deep layer at the end of day j (mm), Ddj − 1 is the same variable on day j − 1. In the same manner as for the root zone, if Ddj < 0 then an amount DPdj of deep percolation is assumed to be lost for the crop (DPdj = −Ddj) and Ddj is set to zero. The water stress of vegetation in the FAO method is expressed by a coefficient Ks related to the actual root zone water content. A fraction p of TAW named Readily Available Water (RAW) is supposed to be available for the plant without stress (Ks = 1). The stress is presumed to start when Dr > RAW and is calculated using Equation (13) (Ks < 1). Conversely, when Dr ≤ RAW then Ks = 1.
1 (13) whereas rainfall inputs can be estimated using meteorological data, the irrigation inputs are very variable in space and time and cannot be known practically on large areas. Thus, the SAMIR tool simulates irrigations based on the daily soil water balance of the root zone. Irrigation is modeled specifically for each land cover class, so as to reproduce the various irrigation practices applied by the farmers on the ground. These practices are described by several parameters that the user has to set. A major parameter is the threshold water content in the root layer to trigger a water input, defining the management allowable depletion MAD [16], i.e., the amount of water which can be depleted between two irrigations. When the root zone depletion reaches this value, irrigation is automatically triggered to fill up the soil. Other parameters to set are the fraction of soil surface wetted by irrigation (fw) depending on the irrigation system, the minimum depth of each input (Min_ir) and the minimum number of days between two water inputs (Min_days). Finally, the Kcb threshold to stop irrigation during the senescence stage (Kcb stop) is defined as a percentage of the peak Kcb value (maximum development) below which irrigation is stopped. In simulating irrigations, SAMIR does not aim at detecting an actual vegetation stress at the time it occurs.
Instead, the hypothesis is that if a significant stress occurs, it will have an effect on the NDVI in the following days which will be accounted for by SAMIR.

Model Calibration and Validation at Plot Scale
The model was calibrated at plot scale using latent heat flux measurements for two seasons. Calibration consists in maximizing the Nash efficiency computed between observed and modeled ET. The Nash-Sutcliffe efficiency coefficient is a non-dimensional statistical performance index that determines the relative magnitude of the residual variance compared to the observed variance [43].
where ETi obs is the observation of evapotranspiration on day i, ETi sim is the modeled value of evapotranspiration on day i and is the observed mean over the entire growing season. The Nash efficiency ranges from −∞ to 1; an efficiency of 1 corresponds to a perfect match between model outputs and observations. Some of the model parameters were taken from the FAO paper 56 [16] or measured in situ (e.g., soil water content), while others were calibrated either because they were not available from the bibliography or because the model was particularly sensitive to these parameters. The procedure to prescribe each parameter value is detailed in the results section. For irrigation, a two-step approach was implemented. First, actual irrigations values resulting from soil moisture measurement analysis were used as inputs in the calibration of the model, then, known irrigations were removed from the model inputs and the automatic irrigation mode was switched on in order to calibrate the irrigation parameters.

Spatialization of ET and Irrigation.
SAMIR was run over the whole irrigated plain using the image time series for the three seasons to compute spatially distributed estimates of irrigation depths. Water balance components, including irrigation, are computed for each pixel. Land use information is required for each pixel since many parameters of the model are crop specific. This is the case for the root zone parameters, the irrigation rules and the fc(NDVI) and Kcb(NDVI) relations. Land cover maps were available for 2008-2009 [44], 2011-2012 [45] and 2012-2013 [46]; they include eight, six and six land use classes, respectively. These classifications were obtained by applying a multi temporal decision tree which allows the identification of crop types on the basis of NDVI thresholds derived from ground truth datasets. The SAMIR parameters required for each land use class were taken either from the previous calibration step (i.e., for cereals) or from bibliographic data since there was no calibration data for market gardening and fruit trees classes. The climatic forcing was considered homogeneous over the plain and meteorological data was taken from the sole station present in the area. Finally, the irrigation parameterization for crops other than cereals was defined from our knowledge of the farmers' practices, while considering the main differences in irrigation practices between classes (e.g., aspersion for cereals, drip irrigation for orchards). In order to validate the SAMIR estimates, irrigation depths were cumulated to monthly values at the scale of irrigation sectors (GDA) and compared to available official irrigation volumes gathered through our ground survey.

Remote Sensing Data Preprocessing
For the SPOT4-take5 time series acquired in 2012-2013, the longest gap was at the beginning of the period, as the first correct image was acquired on 10/03/2015, which means 40 days without image data. This is quite long regarding vegetation monitoring and emphasizes the limitation of a five-day revisit frequency even in semi-arid areas (frequent cirrus clouds can be observed over the study site). However, in our case, this gap was filled using the SPOT5 satellite which successfully acquired two images, thanks to the programming capabilities of this sensor and its oblique viewing agility allowing to observe areas on cloud free days. This is an interesting result showing that combining Sentinel-2 data with other sensors (Landsat 8, SPOT6, etc.) may still be necessary in many places to get consistent high resolution time series. Another way to bridge the gaps in the time series would be to use fusion methods using medium resolution images to estimate high resolution signatures [47].
We present hereafter the outcome of the images pre-processing for the 2012-2013 season; the same approach was applied to the 2008-2009 and 2011-2012 images. After application of an atmospheric correction, we identified manually 28 invariant objects in the scene by visually comparing pairs of distant dates (i.e., 5 November 2012 and 10 June 2013). Then, for these 28 objects, the reflectance of one band at each date is plotted against the reflectance of the average image of this band (Figures 3 and 4). The quality of the invariant objects is confirmed by the determination of the linear fit. However, whereas in some cases the regression fits the 1:1 line (Figures 3a and 4a), in other cases the regression line is significantly different from the 1:1 line, showing a problem in the quality of the atmospheric correction. These discrepancies are more frequent for the SPOT5 time series, which is not surprising as each date was corrected independently, whereas the SPOT4 series was corrected using the MACCS algorithm taking into account the temporal dimension of the series. When the deviation from the 1:1 line was important (Figures 3b and 4b), the linear correction was applied to the image to match it with the average image. For the SPOT4 series, seven dates were linearly corrected, and one date was discarded (20 March) because of the strong scattering of the reflectance due to haze (Figure 4c). For the SPOT5 time series, eight images were corrected.
Once the consistency of reflectance levels within the two time series had been checked, a similar analysis is performed between the two time series by plotting the average reflectance of the invariant objects ( Figure 5). A significant bias was observed which can be explained by (i) differences in the atmospheric correction algorithm used, (ii) difference in band definition between SPOT4 and 5 and (iii) the variations in viewing angle between both sets of images. Indeed the SPOT4 images were acquired at a fixed angle different from nadir, whereas SPOT5 images were acquired at any angle. The observed bias had a strong impact on maximum NDVI values observed in the images, which was 0.9 for the SPOT4 series, and only 0.7 for the SPOT5 series. Considering that fully covering vegetations were certainly present in the area (i.e., cereals or forage fields), a realistic maximum NDVI value of 0.9 was expected for all seasons. Therefore, the SPOT5 series was linearly normalized to match the SPOT4-Take5 radiometry on the basis of the linear regression established between the reflectance of the SPOT5 and SPOT4-Take5 average images ( Figure 5).   The SPOT4 series was delivered with cloud masks that were applied to avoid anomalies in the NDVI. For the SPOT5 series, two images included small cumulus clouds (5 November 2012 and 21 January 2013) which were masked. The clouds were identified using a simple threshold since they have a strong reflectance in the blue band. The cloud shadows were also easy to identify because they had the lowest reflectance in the near infra-red band. As small clouds were rarely at the same place, they have limited impact on the resulting NDVI profiles.
Finally, from the combination of these two time series, a NDVI profile was generated for each pixel.

Plot Scale Calibration of Evapotranspiration Parameters
In a subsequent step, the NDVI time profile for each flux site was extracted from the SPOT time series and was used in the SAMIR model for calibration. The NDVI-fc relationship was determined empirically considering that for a bare soil (NDVI = 0.1) the fraction cover was null (fc = 0) and that at full coverage (fc = 1) the NDVI was the maximum value observed in the image (0.9). The characteristic volumetric soil water contents were determined from soil texture analysis (clay loam, θfc = 0.29, θwp = 0.15). An initial water content of 10% was considered for the two plots since soils were mostly dry after the summer and before the first autumnal rainfalls. The soil fraction wetted by rain or irrigation (fw) was set to one because the irrigation technique for cereals was sprinkler irrigation. The depth of the evaporation layer (Ze) and the proportion of easily available water (p) were fixed following FAO paper 56 recommendations [16]. All other parameters were fixed by calibration (REW, m, Zrmax, Zsoil, Difer, Difrd, Kcb) and are summarized in Tables 1 and 2. In a second phase, irrigation parameters were calibrated by assuming an optimal management to avoid stress (i.e., water input when RAW is empty) and the depth of water inputs was selected in order to fill the depletion (Dr = 0). The calibration was applied simultaneously to both seasons to get a unique set of parameters for cereals. The results (Table 3) show that cereals are irrigated only a short time after the vegetation peak is reached (Kcb_stop = 99%) which is consistent with the conventional agricultural practice for cereals (no more water is needed after grain filling during the maturation stage. The results showed that although evapotranspiration simulations are on the whole correct ( Figure 6), they are better for the barley plot than for the wheat plot (Nash efficiency of 0.6 and 0.53 respectively and a root mean square deviation RMSD of 0.63 and 0.94 mm respectively). However, we also see that this difference might be due to problems in the observed data, as it is very clear that the daily variations of observed ET for wheat are much stronger than for barley and should be considered with care. However, they might be mainly noisy and we see no clue of any significant bias since they reproduce well the seasonal cycle and they are also coherent with the independent ET0 measurements.
The discrepancies were more frequent at the end and at the beginning of the growing season ( Figure 6), when vegetation cover is low and soil evaporation process dominates. The calibrated parameters are shown in Table 1.

Model Parameters Setting for Evapotranspiration and Irrigation Spatialization
Evapotranspiration and irrigation estimates were spatially distributed at perimeter scale for the seasons 2008-2009 (December to June), 2011-2012 (November to May) and 2012-2013 (November to June). Considering the difficulty to get individual parameter values for specific crop types, the land cover typology was grouped into three major classes: market gardening (about 35% of the GDAs' areas), cereals (about 17%) and orchards (about 34%, mainly olive trees). The linear relationship linking fc and Kcb with NDVI for the market gardening were estimated using the FAO paper 56 [16] and satellite data ( Table 2). For bare soil conditions, we assumed like [23] that fc and Kcb values were zero, with NDVI values for bare soils extracted from the images. For full vegetation cover, fc was assumed to be 1 and Kcb was taken as Kcb-mid in the FAO paper 56, while NDVI was also extracted from the images. For olive trees, the same method was applied to estimate the NDVI-fc relation. For Kcb, a relationship between fc and Kcb was obtained from a previous calibration achieved on experimental data for irrigated olive trees in the Haouz plain in Morocco (not shown here) which was considered more representative of the Merguellil area than using bibliographic data. In absence of a detailed soil map, soil properties (θfc, θwp, REW, Ze, Zsoil, Dif er, Dif rd ) were considered homogeneous in the study area and the parameters were taken from calibration at plot scale (Section 3.2). Crop specific parameters (Zrmax and p) were set to calibrated values for cereals and taken from FAO paper 56 [16] for the other land use classes ( Table 3). The initial soil water content (Init_RU) for annual crops was estimated considering the previous precipitations. For the 2008-2009 and 2012-2013, it was also set to 10% of the soil available water (between θfc and θwp) because no significant precipitations were observed since summer, while for the 2011-2012 season initial water content was set at a relatively larger value of 32% of the soil available water due to 85 mm of precipitations recorded ten days before the starting date of the simulation. Higher initial soil filling rate was used for trees (50%) since they are almost continuously irrigated. Regarding irrigation rules (Table 3), market gardening was irrigated using drip irrigation which means that less soil surface is wetted (fw = 25%) and irrigation lasts longer because vegetable require water until harvest (Kcb_stop = 75%). In addition, to reproduce the drip irrigation, MAD was set to a low value (MAD = 0.2 * TAW), which triggers frequent irrigation inputs. Trees are mainly irrigated by gravity all year round (fw = 100%, Kcb_stop = 0).

Comparison between Modeled and Observed Irrigation Volumes
After running the SAMIR model using the image time series over the plain for the three seasons (Figure 7), the monthly values of modeled irrigation were computed for three irrigated perimeters (GDAs) for which validation data were available (Ben Salem II, Mlelsa and Karma II, Figure 8). The cumulated ET values are also plotted on the figure in order to scale irrigation totals to the total water loss occurring during the month. The seasonal water budget was computed for all campaigns (Figure 9) and shows that on the whole the inputs (P and I) are close to evaporative consumption (ET) with little water remaining in the soil.   (Figure 8b), the estimated irrigations are also on the whole satisfying with two noticeable exceptions for the first and last months of the simulation. Indeed, November exhibits a strong overestimation of irrigation which can be due to an error in soil water content initialization but is also at least partially due to the fact that late vegetables (e.g., pepper or tomatoes) are not all removed although still green and not irrigated. The model is not able to manage such partial vegetation cycles and irrigates as if it was the crop that will follow in the crop rotation (i.e., the crop mentioned in the land use map). In May, the discrepancy is not clearly explained but can be also due to the growing importance of vegetables planted mainly in April. For the 2012-2013 season (Figure 8c), irrigations are also on the whole correct but with several exceptions. Overestimations in November and December for Karma II and November for Mlelssa are linked to previous crops as explained for 2011-2012. The problem of May for Mlessa and June for Karma II can be also related to market gardening which is poorly parameterized. When aggregated at seasonal scale, the irrigation estimates for the eight campaigns give a mean absolute percentage error (MAPE) of 25% ( Figure 10) and the overall difference for all campaigns is very low (121 mm irrigation observed for 135 mm modeled). This is an encouraging result considering the fact that (i) the calibration dataset is minimal, and (ii) the calibration protocol affects a limited number of parameters for a limited number of LU classes. Progress in parameterization would require ideally crop specific information, e.g., flux measurements on vegetables of trees, or at least irrigation volumes collected at plot scale so as to be crop specific. However, even with the current irrigation volumes at perimeter scale, it might be possible to achieve a global calibration of some parameters. Although there would be too many degrees of freedom to correctly calibrate all parameters with such aggregated data, we could focus on irrigation practices which are much uncertain though quite sensitive.

Conclusions
We show in this study that the five days acquisition frequency of Sentinel-2, as simulated in the SPOT4-Take5 experiment, will not completely solve the problem of cloudiness even in semi-arid areas like Tunisia. The combination with other VIS-NIR high resolution sensors like Landsat8 or SPOT should still be useful. We also showed that although the radiometric correction of images was performed with special care using the state-of-the-art MACCS algorithm, the invariant analysis proposed here can help improving the time series quality, especially in semi-arid areas where such objects can be easily found. Moreover, although a cloud masking is performed during the MACCS implementation, the subsequent invariant analysis helps identifying and discarding some remaining hazy images. Using these high resolution time series including clear images approximately every 20 days, we have shown that with limited local data and literature review it was possible to estimate irrigation volumes at perimeters scale. The seasonal volumes estimated by this method appear acceptable, even though results at finer timescales (monthly and below) need to be improved, in particular by translating our knowledge of the agricultural practices into algorithmic constraints in the model. Despite these shortcomings, we have demonstrated that combining HR data and simple water balance modeling offers an interesting method to monitor irrigation volumes.