Delft University of Technology Detecting the response of irrigation water management to climate by remote sensing monitoring of evapotranspiration

This study on a large irrigation scheme in Morocco addressed a two-fold question: (a) does the local water management authority adapt water releases to atmospheric water demand ET0-P? and (b) does crop actual evapotranspiration respond to interand intra-annual variability in water releases? We have evaluated the inter-annual variability of ET0-P during the period 1992–2017 and compared its anomalies (i.e., deviations from average) with anomalies in annual water release. Overall, it appeared that anomalies in water release were consistent with anomalies in ET0-P. The actual evapotranspiration (ETa) was estimated using a time series of multi-spectral satellite image data by applying the Surface Energy Balance (SEBAL) algorithm in a dry, wet, and reference year. We have determined the quartiles of the distribution of the ET0-P time series to identify these three years. The dry year was 2015–2016, the wet year was 2014–2015, and the reference (median of ET0-P) was 2013–2014. Finally, we compared seasonal and annual anomalies in ET0-P, ETa and release, Wd of irrigation water. In the period 1992–2017, the relative anomalies in ET0-P and Wd were similar only in 2000–2001 and 2015–2016. The analysis of anomalies in water inflow and stocks confirmed the response in increased Wd following wet years with higher inflow and replenishment of the reservoir. The response of crop water use to the irrigation water supply was evaluated by comparing anomalies in the ratio of actual to maximum ET, i.e., ETa/ETc with anomalies in Wd. As regards the Sidi Bennour, Faregh, and Gharbia districts, the relative anomalies in ETa/ETc and Wd were consistent, i.e., they had the same sign and comparable magnitude. Overall, the study shows that water delivery Wd responds to inter-annual variability in the pre-irrigation season water inflows into the reservoirs, rather than in ET0-P. On the other hand, actual crop water use, i.e., ETa/ETc, does respond to interand intra-annual variability in Wd. This evidence suggests that there is scope for adaptive water management based on a flexible adaptation of water release to interand intra-annual variability in water demand.


Introduction
At the present time, water resources are severely stressed and particularly scarce in arid regions of the world. In North Africa, one of the most water-scarce regions of the world, water use problems are already apparent, aquifers are over-pumped, water quality is deteriorating, and water supply and irrigation services are often rationed with consequences for human health, agricultural productivity, and the environment. The challenge is escalating as the region's population continues to grow, as the water availability is set to fall by 50% by 2050, and if climate change affects weather and precipitation patterns as predicted, the North African region may see more frequent and severe droughts and floods [1]. Therefore, good and sustainable management of this resource is even more important.
In recent years, agricultural regions in North Africa have been subject to extensive and increasing water constraints. Climate change is projected to increase the fluctuations in precipitation and surface water supplies, which will affect the crop water requirements. These water challenges are expected to strongly impact agriculture, undermining the productivity of rain-fed and irrigated crops, and livestock activities. These changes could, in turn, further impact markets, trade, and broader food security. Irrigated agriculture remains the largest user of water accounting about 83% of annual freshwater withdrawal in North Africa [2] will be important for policymakers to focus on efforts that increase the overall efficiency of water used by the agricultural sector, reduce the sector impact on freshwater resources, and improve its resilience to water risks.
Adaptive water management is an emerging concept [3], where more emphasis is placed on variability and seasonal to annual forecast in contrast to the traditional approach of designing and operating water management systems on the basis of long-term averages and climatology. Adaptive management is a decision process that "promotes flexible decision-making that can be adjusted in the face of uncertainties; e.g., water deficit and water supply, as outcomes from management actions and other events become better understood" [4]. Adaptive irrigation water management requires seasonal to annual forecasting of, e.g., precipitation, evapotranspiration, and water demand in addition to current hydrological state variables, e.g., soil water content [5]. Although the practical implementation of adaptive irrigation water management is a very active area of research and development, the literature provides evidence of ways and means to advance in that direction. Tanaka et al. (2006) [6] analyzed the potential of the California water supply system to adapt to long term changes in water availability and demand in relation with climate. Drieschova et al. (2008) [7] analyzed treaties governing shared international river waters to assess how the response to inter-annual variability in available water was regulated.
Interventions aiming at improving the efficiency of water use in large irrigation schemes, such as the Hetao Irrigation District in China, require accurate planning and operation of water distribution in response to the inter-annual variability in water supply and demand [8]. The rehabilitation of this irrigated area was meant to reduce water diversion from the Yellow River without reducing agricultural production. Water diversion was reduced from 768 mm in 2000 to 576 mm in 2010, while both the atmospheric water demand ET 0 -P and the actual evapotranspiration (ET a ) remained rather stable and lower than ET 0 , notwithstanding large inter-annual variability in water diversion.
This shows that the water management system could adapt to the reduction in water diversion by keeping constant agricultural water use (i.e., ET a ). It should be noted that water diversion was comparable with ET 0 -P in 2000, and much smaller in 2010.
In Spain, river basin organizations (RBO) are responsible to develop river basin management plans (RBMP) by taking into account the expected evolution of water availability and water demand in a long period [9]. The RBMP-s are updated every five years. Within each RBO the Withdrawal Commission is responsible for the yearly adjustment of the water allocation plan in three stages from September to May. The adjustment is mainly based on current water storage in the reservoirs. In Western Morocco, our study area, the procedures to adapt water allocations to the inter-annual variability of water availability and demand are similar to the procedures applied in Spain [10]. In particular, the initial water allocations plan is adapted during the irrigation season to the evolution of climate conditions, although a combination of precipitation and storage in water reservoirs, rather than ET 0 -P only, is used as an indicator of climate conditions and water availability.
From the point of view of irrigation water management, drier years are associated with larger irrigation water demand, determined by the different ET 0 -P between atmospheric water demand ET 0 and precipitation (P). Alternate metrics of water deficit might be applied and are briefly reviewed below. Zargar et al. (2011) [11] reviewed indicators of water deficit aiming at capturing the difference between precipitation and evapotranspiration. These authors highlighted the wider relevance of indicators based on precipitation and potential evapotranspiration. Tsakiris and Vangelis (2005) [12] proposed the Reconnaissance Drought Index (RDI), specifically designed to capture water deficit and defined as the ratio of the cumulative values of precipitation and potential evapotranspiration over a period of time.
The Palmer Drought Severity Index (PDSI) is very widely applied in the US and regards drought as an imbalance between water supply and water demand [13]. Zargar et al. (2011) and Vicente-Serrano et al. (2010) [11,14] proposed a new drought index combining precipitation and potential evapotranspiration into a water balance-based metric of water deficit. Zargar et al. (2011) [11] grouped indicators according to the nature of drought; particularly, they reviewed different indicators of meteorological and agricultural drought. The implementation of adaptive irrigation water management requires an indicator of crop water use to detect timely agricultural drought and an indicator of climatological water deficit to detect meteorological drought.
Accurate characterization of spatial and temporal variability of actual evapotranspiration is vital information for accurate and dynamic irrigation management [8,15]. These authors presented a detailed analysis of the seasonal and inter-annual variability of actual evapotranspiration based on Eddy-covariance (EC) observations and remote sensing retrievals, in response to precipitation and irrigation water.
The study by Zeng et al. (2010) [16] provides a good illustration of the spatial variability in ET a in relation with hydrological processes: the highest ET a was observed for water-bodies, wetland, and woodland where water is concentrated by surface and subsurface flow, while the lowest ET a was observed for grasslands where precipitation was the only water source. Accordingly, the actual water deficit ET a -P was largest for wetlands and waterbodies, thus providing a measure of water supplied by sources other than precipitation. We can, therefore, evaluate the response of irrigated lands to the inter-annual variability in water supply by analyzing the inter-and intra-annual variability in ET a -P.
Monitoring of ET a can be done by, e.g., EC devices, but capturing both, spatial and temporal variabilities is clearly needed for precision and adaptive irrigation water management. Very few experiments have been reported in the literature, (e.g., Li et al. (2013) [17]) where a sufficiently large number of EC devices were put in place to capture the spatial variability of actual crop water use. On the other hand, literature on studies where ET a has been estimated, mapped, and monitored using remote sensing methods, is abundant (see e.g., Gowda et al. 2008;Li et al. 2009 [18,19]). In addition, several reviews and evaluations against in situ measurement on the accuracy of remote sensing retrievals of actual ET have been published (e.g., Hu et al. 2014 [20]).
Various models have been proposed for estimating ET with remotely sensed data over the past decades, such as SEBI [21], TSEB [22], SEBAL [23,24], S-SEBI [25], SEBS [26], LSA-SAF MSG ET [27], STSEB [28], GLEAM [29], and MODIS-ET [30]. Ruhoff et al. (2013) [31] evaluated the daily MODIS ET product against EC measurement in the Rio Grande basin: RMSE was 0.78 and 0.46 mm d −1 at two different EC sites, the basin-wide mean annual ET was 21% lower than the estimate from a hydrological model. Hu et al. (2014) [20] evaluated MODIS and LSAF-SAF MSG ET data product against EC measurement at 50 locations in Europe. RMSE varied from 0.33 to 1.57 mm/d for MODIS ET and from 0.27 to 1.27 mm −d for LSAF-SAF MSG ET a . In summary, these models generally performed well, with the relative errors of 5-30% in comparison with ground-based flux measurements on different ecosystems around the world [8]. A detailed discussion of remote sensing-based ET a estimation models can be can be found in Gowda et al. (2008) and Li et al. (2009) [18,19].
The evidence we have reviewed above suggest an overall accuracy of monthly to seasonal ET a data products between 5% and 10%. In particular the accuracy of SEBAL estimates of actual ET is well documented and the relative error has been reported to be around 10%. We have applied SEBAL to estimate ET a , but we have used only anomalies in ET a , rather than its absolute value, in our analysis. Accordingly, we have concluded that the documented accuracy of SEBAL was appropriate for our purpose and did not carry out additional experiments to validate our estimates of ET a . More generally, we have also concluded that remote sensing retrievals of ET a can be applied to assess the effectiveness of adaptive irrigation water management to cope with inter-and intra-annual variability of water availability and atmospheric water demand.
Previous studies, e.g., [10], in the general area of information services to support irrigation management in the Doukkala irrigation scheme in Morocco have been designed and carried out by assuming that sound irrigation water management should be based on observable indicators of water demand. Moreover, in the Doukkala irrigation water is delivered to farmers if they request it. This lea to defining the broader issue to be addressed by our study: "Does the irrigation water management system responds to irrigation water demand? If this is the case, is this response observable?" This two-fold issue can be dealt with by addressing the two research questions (a) is water diversion in the Doukkala irrigation scheme coherent with irrigation water demand and are water diversions adapted yearly to the expected irrigation water demand? (b) Do changes in water diversion have a negative impact on water availability to crops? In principle, the inter-annual climate variability may be compensated by water management by varying water releases on the basis of the expected evaporative demand. Accordingly, we have first evaluated the intra and inter-annual variability of water release in relation with evaporative demand, taking into account (P). Then we analyzed, in detail, a dry and a wet year to assess the variability in ET a and whether this variability was consistent with changes in water releases.

Description of the Study Area
The current study was conducted in Doukkala region which is located in Western Morocco (32 • 15 N and 33 • 15 N, 7 • 55 W and 9 • 20 W), in the downstream portion of the Oum Er-Rbia basin ( Figure 1), This area includes two geomorphological units: the Doukkala and the Sahel. The plain of Doukkala covers approximatively 35,000 km 2 . The elevation is about 120-130 m above sea level [34].
The irrigated area of Doukkala is one of the largest irrigated areas in Morocco (currently 96,000 ha), remarkable for its size and national relevance. It contributes 38% of sugar beet and 20% of commercialized milk [35]. Sugar beet, spring wheat, alfalfa (perennial crop) and summer maize are the main crops. The area is characterized by a semi-arid climate, being temperate and mild in winter and warm and dry in summer. In the period (1992-2017) mean annual rainfall was about 308 mm/year and mean annual reference evapotranspiration ET 0 about 1358 mm/year. The Doukkala irrigation scheme includes two sub-areas: the higher section with an irrigated area of 35,000 ha and lower section with irrigated area of 61,000 ha ( Figure 1). The latter includes four main districts: Faregh, Sidi Bennour, Zemamra and Gharbia, from the East to the West. Each district is divided into a number of Centers of Irrigation Management (CGRs) irrigated with different irrigation systems. The Doukkala irrigation scheme is managed by the Regional Office of Agricultural Development of Doukkala (ORMVAD) [10].
The dominant irrigation technique is surface irrigation, sprinkler irrigation is applied in a relatively large area, while drip irrigation is applied in a rather small area. According to Brouwer et al. (1989) [36], the efficiency of water application with these methods is approximately 60%, 75%, and 90%, respectively. the ORMVAD technical staff in charge of each CGR carries a final check on whether the requested water volume can, and should be, delivered. This may include corrections to take into account precipitation preceding the rotational interval or short-term precipitation forecasts. In case of heavy precipitation, no irrigation water is delivered. It should be noted that this is a short-term risk-edging process, that does not evaluate a longer temporal horizon, e.g., saving irrigation water to get prepared for a drought expected later in the irrigation season.

Meteorological Data
We have used the observations of precipitation (P) and the estimates of reference crop evapotranspiration ET0 collected by ORMVAD at the three stations indicated in the (Table 1) for the period 1992-2017.
The observations at the three stations have been merged to create one continuous time series. More precisely, the observations at Gharbia and Mettouh stations were used to fill gaps in the Sandy soils are dominant with fine and coarse texture. These soils consist of approximatively 62% of sand, 23% of clay, and 15% of silt [37].
The water resource mobilized to irrigate the Doukkala comes mainly from the Al Massira dam, one of the main water storage structures in the Oum Er-Rbia basin, with a capacity of about 2.7 × 10 9 m 3 . Al Massira reservoir is one of the largest reservoirs in Morocco, with multiple uses: drinking water, agriculture, industry, and energy.
Water delivery in each rotational interval (15 days) is calculated by ORMVAD staff on the basis of current irrigated area, i.e., determined on the basis of requests by farmers for each interval. The gross irrigation water depth delivered to farmers is 864 m 3 /ha, equivalent to 172.8 mm/month for each hectare registered yearly with ORMVAD. This implies that the total water delivery to a farm or to a CGR, is variable in response to the requests for water by farmers. Prior to actual delivery of water, the ORMVAD technical staff in charge of each CGR carries a final check on whether the requested water volume can, and should be, delivered. This may include corrections to take into account precipitation preceding the rotational interval or short-term precipitation forecasts. In case of heavy precipitation, no irrigation water is delivered. It should be noted that this is a short-term risk-edging process, that does not evaluate a longer temporal horizon, e.g., saving irrigation water to get prepared for a drought expected later in the irrigation season.

Meteorological Data
We have used the observations of precipitation (P) and the estimates of reference crop evapotranspiration ET 0 collected by ORMVAD at the three stations indicated in the (Table 1) for the period 1992-2017.
The observations at the three stations have been merged to create one continuous time series. More precisely, the observations at Gharbia and Mettouh stations were used to fill gaps in the observations at Zemamra (Table 1). In this study we took this time series as a reference for the four districts of Sidi Bennour, Zemamra, Gharbia, and Faregh ( Figure 1). In addition, the hourly data collected at the same stations for the irrigation seasons (2013-2014, 2014-2015 and 2015-2016) have been used to estimate and map (ET a ) using SEBAL (see Section 3.5). We have used hourly observations of air temperature (T a ), relative humidity (RH), wind speed (U), precipitation (R) and solar radiation (Rs) on the same dates as the satellite overpasses.

Data on Water Delivery
The water management authority (ORMVAD) collects daily data on water delivered to each CGR in the Doukkala. ORMVAD provided data for the period 1992-2017 aggregated to the entire irrigation season for each CGR. In addition, data on the water delivered to each CGR in each quarter ( Figure 1) by the district water management office (Arrondissement ORMVAD) were collected for the irrigation seasons 2013-2014, 2014-2015 and 2015-2016. The irrigation season is from September to August in the following year. The selection of the three irrigation seasons is explained in Section 3.2. In this study we use the wording water delivery and the symbol (W d ) to indicate the water deliveries described in this paragraph. The data on water delivery provided by ORMVAD include a correction to account for conveyance and the distribution losses, estimated to be 35% of water allocated to irrigation at the Al Massira reservoir.

Satellite Data
All available LANDSAT-8-OLI and TIRS images for the irrigation seasons 2013-2014, 2014-2015, and 2015-2016 (see the previous Section 2.3 for additional explanations) with a nominal temporal resolution of 16 days and spatial resolution of 30m at nadir, were downloaded from https://earthexplorer. usgs.gov/. The technical specifications of the Landsat-8 data products can be found in USGS (2016) [38]. A few images were not usable due to cloud cover (see Table 2). All satellite images were geometrically corrected to the system of coordinates: UTM, WGS-84, zone 29 and atmospherically corrected using DOS (dark object subtraction) method [39]. For the VNIR and SWIR spectral ranges, this procedure provides at-surface spectral reflectance. The procedure applied to the TIR spectral range is described in Section 3.3.

Work Flow
The approach briefly described in the introduction leads to the detailed workflow in Figure 2, which shows how the evaluation of the net atmospheric water demand relates to the detailed analysis on actual crop evapotranspiration and how the different streams of data analysis lead to the evaluation (see red arrow, Figure 2) of the response of irrigation water management to irrigation water demand and of crops to water supply. The indicators of water demand are evaluated by merging meteorological and remote sensing data (see blue arrow, Figure 2). The aerodynamic resistance, once corrected for stability is applied to determine the relationship between dT and T s (see brown arrow, Figure 2). The elements of the workflow are described in some detail in the following sections.

Work Flow
The approach briefly described in the introduction leads to the detailed workflow in Figure 2, which shows how the evaluation of the net atmospheric water demand relates to the detailed analysis on actual crop evapotranspiration and how the different streams of data analysis lead to the evaluation (see red arrow, Figure 2) of the response of irrigation water management to irrigation water demand and of crops to water supply. The indicators of water demand are evaluated by merging meteorological and remote sensing data (see blue arrow, Figure 2). The aerodynamic resistance, once corrected for stability is applied to determine the relationship between dT and Ts (see brown arrow, Figure 2). The elements of the workflow are described in some detail in the following sections.

Analysis of Climate Data
The first step in the work-flow described in Section 3.1 is the evaluation of the net atmospheric water demand ET 0 -P. We have calculated first daily ET 0 -P for each year in the period 1992-2016. The precipitation P is the time series obtained by merging available observations as described in Sect, 2 (meteorological data).
The reference crop evapotranspiration ET 0 has been calculated by applying the Penman-Monteith equation (Equation (1)) adapted from the original equation of Monteith, (1965) [40] applied to a reference crop conditions according to FAO-56 [41]: where ET 0 = reference evapotranspiration (mm d −1 ); R n = net radiation at crop surface (MJ m −2 d −1 ); G = soil heat flux density (MJ m −2 d −1 ); R n is determined by the albedo of a standard crop and G is estimated as explained in Section 3.5; γ = psychrometric constant (kPa • C −1 ); T = mean daily air temperature at 2 m height ( • C), u 2 = wind speed at 2 m height (m s −1 ), e s = saturation vapor pressure (kPa), e a = Actual vapor pressure (kPa), e s − e a = Saturation vapor pressure deficit (kPa), ∆ = slope vapor pressure curve (kPa • C −1 ). The reference crop is assumed to have albedo = 0.23, crop height h c = 0.12 m, and minimum canopy resistance r c,min = 70 s m −1 . The daily P and ET 0 values were aggregated to the entire irrigation season from September to August in the following year to obtain the annual net atmospheric water demand. Note that the annual values were attributed to the same year as the months January to August; see Section 3.6.1 for further details about this analysis. Next, a wet, dry and reference years were determined in this period by analyzing the frequency distribution of the annual ET 0 -P. Quarterly values of ET 0 -P for the wet, dry and reference years were calculated using the daily meteorological observations described in Section 2.2 over the entire irrigation season.

Preprocessing of Satellite Data
Landsat 8 LT1T products have been downloaded and processed using t 'Preprocessing' module of the 'Semi-automatic classification' plugin of QGIS [42]. Preprocessing of Landsat8 images included the following steps: Conversion from digital numbers (DNs) to Top of Atmosphere (TOA) radiance (L λ ) is done for each VNIR and SWIR band using the band gain (ML) and offset (AL) (both present in the product metadata file).
Atmospheric correction in the VNIR and SWIR range to obtain bottom of atmosphere (BOA) radiance was performed using the revised dark object subtraction method (DOS) described by Chavez, (1996) [39]. This approach estimates the atmospheric path radiance L p by observing dark targets in each image, i.e., by assuming that the at-satellite radiance is due to the atmosphere only. This method estimates atmospheric transmittance using a very simple parameterization evaluated by Chavez (1996) against field measurements and accurate radiative transfer models. The algorithm present in the QGIS plugin calculates L p for each band as the difference between the maximum radiance L min recorded among the 0.01% darkest pixels and the presumed radiance LD01% of dark objects.
In the TIR spectral range we retrieved LST by applying a four-steps procedure. In Step 1 we converted the DN-values in at-satellite radiance as done for the VNIR-SWIR spectral range. In Step 2 we converted the at-satellite Radiance into at-satellite brightness temperature using Equation (2): where T = at-satellite brightness temperature (K); L λ = TOA spectral radiance (Watts/(m 2 × rad × µm)); K 1 and K 2 band-specific thermal conversion constants; values included in the metadata for each image. In Step 3, the Land Surface Emissivity is estimated from NDVI and in Step 4 LST is retrieved according to Jimenez-Munoz et al. (2014) [43] by applying the Split Window algorithm provided by Latif (2014) [44].

General
The maximum crop evapotranspiration ET c was estimated by the single crop coefficient approach where crop transpiration and soil evaporation are combined into a single K c coefficient [41]: where ET 0 was estimated as described in Section 3.2.

Estimation of K c Using NDVI Observations
Several authors found a good correlation between the Kc and Normalized Difference Vegetation Index [45,46], where the Normalized Difference Vegetation Index (NDVI) is defined as:  [46][47][48][49] compiled many linear relationships between K c and NDVI to conclude that a generic, crop independent relationship could be used:

Retrieval of Actual Evapotranspiration ET a : The SEBAL Model
The Surface Energy Balance Algorithm for Land (SEBAL) model, developed by Bastiaanssen (1995), Bastiaanssen et al. (1998a) and Bastiaanssen et al. (1998b) [23,24,50] accounts for bio-physical processes at the land surface and can be used to retrieve the actual evapotranspiration ET a from multispectral imagery.
SEBAL is a single -source resistance transfer scheme based on remotely sensed data and limited weather data for deriving components of the surface energy balance equation. SEBAL computes the complete radiation and energy balance components along with the resistances for momentum, heat, and water vapor transport for each pixel [23,24]. The input data for SEBAL consists of spectral radiance in the visible, near-infrared, and thermal infrared regions of the spectrum. In addition to satellite images, the SEBAL model requires weather observations, i.e., wind speed, air humidity, solar radiation, and air temperature. Since the satellite image provides information for the overpass time only, SEBAL computes an instantaneous ET a at the time of image acquisition. More precisely, the latent heat flux λE associated with ET a , i.e., λET a is calculated for each image pixel as a "residual" of the surface energy budget equation: where λE is the latent heat flux (W/m 2 ), R n is the net radiation flux at the surface (W/m 2 ), H is the sensible heat flux to the air (W/m 2 ), G is the soil heat flux (W/m 2 ). The daily value of λE is calculated by assuming that the evaporative fraction f = (λE/ET 0 ), and remains constant during the day so that the daily value of λEd = EF × ET 0 , d.

Net Radiation R n
The net radiation R n is obtained from the surface radiation balance, i.e., as the difference of outgoing and incoming radiant fluxes. The surface radiation balance equation was applied as described by [51].

Soil Heat Flux (G)
Soil heat flux is the rate of heat energy transferred from the earth's surface to the subsurface. It is calculated as an empirical fraction of the net radiation using surface temperature (T s ), surface Albedo (α) and NDVI [52]: If the NDVI < 0, we assume that the observed target is water and for this condition we applied G/R n = 0.5

Sensible Heat Flux (H)
Sensible heat flux is the heat loss to the air by conduction and convection due to a temperature difference. The sensible heat flux is computed as: where ρ is the air density (kg/m 3 ), cp is the specific heat of air at constant pressure (J/kg K), dT (K) is the vertical near-surface temperature difference, and r ah is the aerodynamic resistance to heat transport (s/m). The sensible heat flux (H) is, therefore, a function of the temperature gradient, surface roughness, and wind speed. There are two unknowns in Equation (8), r ah and dT, to be estimated. dT represents a temperature difference (T 1 − T 2 ) between two heights (z 1 and z 2 ). This temperature difference is used because satellites measure radiometric temperature which can differ from aerodynamic temperature by several degrees [53]. It is assumed that a linear relationship exists between dT and T s : This is a key element of the entire SEBAL approach, since this linear relationship is applied in a self-calibration procedure as explained below.
The coefficients a and b in Equation (9) are determined by fitting this linear relationship to two reference values of dT assumed to apply at the surface temperature of reference hot/dry and cold/wet targets, to be selected within the observed scene. Ideally H cold = dT cold = 0 at the surface temperature of the Cold/Wet target and Hhot = R n − G, dT = dT hot at the surface temperature of the hot/dry target. The selection of the cold and hot reference targets is critical to the accuracy and reliability of the SEBAL estimates of ET a , especially when time series of ET a are being generated.
The selection of the cold/wet and hot/dry reference targets has been done in different ways in SEBAL applications reported in literature. In this study we have selected "cold" pixels by sampling water bodies where H cold = dT cold = 0. The water bodies are identified by a water mask generated beforehand by using a multi-spectral satellite image at high spatial resolution. On these pixels, the surface energy balance equation is of the form: R n = G + λET. The "hot" reference targets are identified by selecting the pixels in the upper 20% tail of the T s distribution in each image. The "hot" reference temperature is estimated by calculating the average of pixels selected in this way. In the "hot" pixels we assume that there is no latent heat flux, i.e., λET = 0 and H = R n − G.
Temporal integration of instantaneous ET a values: Our objective is to estimate actual crop water use over a long period of time, i.e., we need to estimate quarterly and seasonal ET a . This requires a two-step temporal integration of the instantaneous ET a estimated as described above.
The integration of instantaneous to daily ET a is done by assuming that the evaporative fraction f = ET a /ET 0 remains constant during daytime or at least during the central hours of the day when the contribution to daily ET a is largest. With regard to the period in-between two sub-sequent image acquisitions, the same assumption f = const is applied to obtain the daily ET a values from daily ET 0 and finally the daily ET a is added up to obtain the required quarterly and seasonal ET a . For the period 1992-2017 and for the entire Doukkala irrigation scheme, we calculated the anomalies, i.e., the deviations from average, of ET 0 -P and of the water deliveries (see Section 2.3) to assess whether the water management system responds to water demand by adjusting water deliveries. To compare with ET 0 -P, expressed in mm for the entire irrigation season, the total water delivery to the entire Doukkala was divided by the total irrigated area, as determined by ORMVAD, to obtain water delivery in mm for the entire irrigation season. It should be noted that the total irrigated area was 61,000 ha from 1992-1993 to 2000-2001, when the higher Section (see Section 2.1) became fully operational, and 96,000 ha after that, i.e., from 2001-2002 to 2016-2017. The data on water delivery is normally collected by CGR (see Section 2.3), while ORMVAD provided the aggregated values described above for this study.

Evaluation of the
For the years 2013-2014, 2014-2015, and 2015-2016, ORMVAD provided the data on water deliveries disaggregated by quarter and by CGR, along with the corresponding irrigated area. The quarterly data on volumetric water delivery was divided by the irrigated area to obtain the quarterly water deliveries in mm for each CGR.
The response of the water management system to water shortages was evaluated with both data sets by comparing both absolute and relative anomalies in ET 0 -P and water delivery. The comparison was done by evaluating both the timing and the correlation of such anomalies.

Evaluation of the Response of Crop Water Use to Water Delivery
Given the observed inter-annual and intra-annual (inter-quarter) variability in water deliveries and ET 0 -P, we have evaluated the impact of this variability on actual crop water use, as captured by the ET a maps, generated by applying the SEBAL procedure (see Section 3.5).
On the basis of the analysis described in Section 3.6.1, we have identified (see Section 4.1) in the period 1992-2017 a dry year, i.e., the largest ET 0 -P, a wet year, i.e., the smallest ET 0 -P, and a reference year, i.e., the median ET 0 -P. The availability of satellite data restricted this analysis to the years 2013-2014 (reference), 2014-2015 (wet), and 2015-2016 (dry) (see Section 2.4).
We have applied two indicators of crop water use, i.e., ET a and the ratio ET a /ET c , and compared the absolute and relative anomalies in both indicators with the anomalies in water deliveries.

Inter-Annual Variability 1992-2017
The region of Doukkala is characterized within a year by clearly separated dry, i.e., with ET 0 P, and wet, i.e., with ET 0 P periods (Figure 3). In addition, the inter-annual variability of precipitation is rather large, while is relatively low in ET 0 (Figure 3). Extreme hydrometeorological events are not rare in the autumn and winter, e.g., 279 mm in December 1996 and 196 mm in November 2014. The high values of summer ET 0 , however, lead to a very large intra-annual variability of ET 0 -P. This led us to investigate in detail the possible response of W d to both the inter and intra-annual variability of ET 0 -P. events are not rare in the autumn and winter, e.g., 279 mm in December 1996 and 196 mm in November 2014. The high values of summer ET0, however, lead to a very large intra-annual variability of ET0-P. This led us to investigate in detail the possible response of Wd to both the inter and intra-annual variability of ET0-P.
This equation is approximate since it does not include other water users besides the Doukkala irrigation scheme.
The exceptionally high precipitation in the year 1995-1996 led to a sharp increase in the inflow into the water reservoir of Al Massira and to its replenishment with a delay of approximately one year. Wd increased gradually as the reservoir filled up. Following the dry period of 2001-2009, during which water inflow, stock and Wd were quite consistent, a wet period followed, and the reservoir filled up at the same time as the inflow increased. In the same period the Wd remained rather constant, i.e., the improved water availability did not lead to higher Wd, although the latter remained much lower than ET0-P. It should be recalled that that the Higher Section (35,000 ha) became fully operational in the irrigation season 2001-2002, thus raising new challenges in the equitable allocation of available irrigation water. We have evaluated ET 0 -P as an indicator of atmospheric water demand (see Section 3.2). It should be noted that W d was rather stable and around 482 mm year −1 in the period 2013-2016. In earlier years the inter-annual variability is rather large in both ET 0 -P and W d (Figure 4). The former varies between a minimum of 767.5 mm year −1 and a maximum of 1305.8 mm year −1 . The latter, between a minimum of 191.6 mm year −1 and a maximum of 708.6 mm year −1 . Overall, W d is much smaller than ET 0 -P and there is no evident correlation between ET 0 -P and W d . The reservoir water balance followed approximately the rule:  The analysis of anomalies in water inflow and stocks confirms the response in increased Wd following wet years with higher inflow and replenishment of the reservoir. The time lag between changes in water stock and water inflow was approximately one year from 1995-1996 (peak inflow) to 1996-1997 (maximum water stock) while it took two years to reach the maximum Wd in 1999-2000. This delay seems determined by the reservoir operation, but it was beneficial since ET0-P was much larger than average in 2000-2001. In the extended dry period 2001-2009 the reduction in Wd was mitigated by exploiting the water stored in the preceding wet years, but anomalies in Wd remained  This equation is approximate since it does not include other water users besides the Doukkala irrigation scheme.

Water deliveries (mm) Inflow (mm) Stock (mm) ET0-P(mm)
The exceptionally high precipitation in the year 1995-1996 led to a sharp increase in the inflow into the water reservoir of Al Massira and to its replenishment with a delay of approximately one year. W d increased gradually as the reservoir filled up. Following the dry period of 2001-2009, during which water inflow, stock and W d were quite consistent, a wet period followed, and the reservoir filled up at the same time as the inflow increased. In the same period the W d remained rather constant, i.e., the improved water availability did not lead to higher W d , although the latter remained much lower than ET 0 -P. It should be recalled that that the Higher Section (35,000 ha) became fully operational in the irrigation season 2001-2002, thus raising new challenges in the equitable allocation of available irrigation water.
Likewise, the annual anomalies in W d were not well correlated with ET 0 -P ( Figure 5). The absolute anomalies (Figure 5a) do not suggest a clear response of W d to evaporative demand. In the period 1992-2017, the relative anomalies in ET 0 -P and W d (Figure 5b)  The analysis of anomalies in water inflow and stocks confirms the response in increased Wd following wet years with higher inflow and replenishment of the reservoir. The time lag between changes in water stock and water inflow was approximately one year from 1995-1996 (peak inflow) to 1996-1997 (maximum water stock) while it took two years to reach the maximum Wd in 1999-2000. This delay seems determined by the reservoir operation, but it was beneficial since ET0-P was much larger than average in 2000-2001. In the extended dry period 2001-2009 the reduction in Wd was mitigated by exploiting the water stored in the preceding wet years, but anomalies in Wd remained negative for the entire period. The increase in water inflow in 2009-2010 allowed to replenish the water stock depleted during the dry years, but Wd was kept rather constant and close to the average. The latter was consistent with ET0-P slightly fluctuating around the average except in 2015-2016 when ET0-P was significantly larger than average and a large fraction of the water stock was depleted to increase Wd (Figure 5a, b).

Selection of Dry, Wet, and Reference Years
To identify the dry, wet, and reference years (see Section 3.2) we have used the 25 yearly values (Figure 3) of ET0 and P provided by ORMVAD, where ET0 was estimated according to FAO 56 [41], as detailed in Section 3.2. Next, we have determined the quartiles of the distribution of ET0-P (Table  3). By crossing the years for which the high-resolution satellite data are available with this distribution, we see that the year 2015-2016 falls between Q3 and the maximum, i.e., the range of highest atmospheric water demand, the year 2014-2015 between the minimum and Q1, i.e., the range of lowest water deficit. The ET0-P in the year 2013-2014, i.e., 1005 mm year −1 , is slightly lower than the median ET0-P, i.e., 1042 mm year −1 , thus we took 2013-2014 as a reference year. ET0-P was 914 mm year −1 in 2014-2015 and 1296 mm year −1 in 2015-2016.
We can then conclude that the years for which satellite data are available to fit well our objective, notwithstanding the limitations due to the available satellite datasets.  The analysis of anomalies in water inflow and stocks confirms the response in increased W d following wet years with higher inflow and replenishment of the reservoir. The time lag between changes in water stock and water inflow was approximately one year from 1995-1996 (peak inflow) to 1996-1997 (maximum water stock) while it took two years to reach the maximum W d in 1999-2000. This delay seems determined by the reservoir operation, but it was beneficial since ET 0 -P was much larger than average in 2000-2001. In the extended dry period 2001-2009 the reduction in W d was mitigated by exploiting the water stored in the preceding wet years, but anomalies in W d remained negative for the entire period. The increase in water inflow in 2009-2010 allowed to replenish the water stock depleted during the dry years, but W d was kept rather constant and close to the average. The latter was consistent with ET 0 -P slightly fluctuating around the average except in 2015-2016 when ET 0 -P was significantly larger than average and a large fraction of the water stock was depleted to increase W d (Figure 5a,b).

Selection of Dry, Wet, and Reference Years
To identify the dry, wet, and reference years (see Section 3.2) we have used the 25 yearly values (Figure 3) of ET 0 and P provided by ORMVAD, where ET 0 was estimated according to FAO 56 [41], as detailed in Section 3.2. Next, we have determined the quartiles of the distribution of ET 0 -P (Table 3). By crossing the years for which the high-resolution satellite data are available with this distribution, we see that the year 2015-2016 falls between Q3 and the maximum, i.e., the range of highest atmospheric water demand, the year 2014-2015 between the minimum and Q1, i.e., the range of lowest water deficit. The ET 0 -P in the year 2013-2014, i.e., 1005 mm year −1 , is slightly lower than the median ET 0 -P, i.e.,  We can then conclude that the years for which satellite data are available to fit well our objective, notwithstanding the limitations due to the available satellite datasets.

Intra-Annual Variability
Assuming that any difference in W d in relation with ET 0 -P would be more evident when comparing hydrologically different years, we focused our analysis of intra-annual variability on the dry, reference and wet years described above. To evaluate whether quarterly W d is determined to guarantee adequate water supply through the yearly sequence of low and high-water deficit periods, we have compared quarterly values of anomalies in ET 0 -P and W d for each district of the Doukkala irrigation scheme ( Figure 6). The anomalies are the deviations from the mean water delivery over all quarters and years, calculated separately for each district. Overall, the anomalies in W d are much smaller than the anomalies in ET 0 -P, although higher than average W d coincides with higher than average ET 0 -P. The latter suggests that the water management authority aims to increase W d when evaporative demand is larger, although the increase in water delivery is much smaller than the increase in evaporative demand. This feature applies to all districts although it is less evident in the case of Zemamra. To interpret the results in Figure 6 it should be taken into account that the annual W d in the years 2013-2016 was rather stable (Section 4.1.1 and Figure 4). It should also be taken into account that the quarterly values shown in Figure 6 relate to the sequence of the reference, wet and dry years. It should be noted that water deliveries in the 3rd quarter are slightly higher than average, if at all, although the atmospheric water demand in the same quarter is always much larger than average (i.e., large positive absolute anomalies).
Wd in the years 2013-2016 was rather stable (Section 4.1.1 and Figure 4). It should also be taken into account that the quarterly values shown in Figure 6 relate to the sequence of the reference, wet and dry years. It should be noted that water deliveries in the 3rd quarter are slightly higher than average, if at all, although the atmospheric water demand in the same quarter is always much larger than average (i.e., large positive absolute anomalies).

Mapping and Monitoring of Crop Water Use
As explained above we have calculated and mapped ETa for the three selected years. The results obtained for late-spring and summer (Figure 7) are consistent with the expected impact of the dry, wet and average conditions in the three years. The highest ETa values were observed in the year with the highest atmospheric water demand, i.e., the dry year 2015-2016. The lowest ETa values in the year with the lowest atmospheric water demand, i.e., the wet year 2014-2015, while ETa values in-between these two extremes were observed in the reference year 2013-2014. In addition, the spatial variability of ETa was largest in May 2015-2016. This leads to a hypothesis, i.e., that water supply meets actual water requirements in all conditions and ETa increases in response to the ET0-P.
The evolution through the year 2015-2016 ( Figure 8) reflect the expected seasonality, but the higher ET0, which characterized this year lead to a larger winter-summer difference in ETa. Likewise, the drier conditions in 2015-2016 led to a larger spatial variability in ETa in spring and summer. Differences up to 6 mm d −1 can be noticed in ETa on DoY 157, which suggest large differences in the distribution and use of irrigation water, these results clearly require detailed data on water delivery to and within the areas served by tertiary canals for a better understanding of irrigation water management in the Doukkala.

Mapping and Monitoring of Crop Water Use
As explained above we have calculated and mapped ET a for the three selected years. The results obtained for late-spring and summer (Figure 7) are consistent with the expected impact of the dry, wet and average conditions in the three years. The highest ET a values were observed in the year with the highest atmospheric water demand, i.e., the dry year 2015-2016. The lowest ET a values in the year with the lowest atmospheric water demand, i.e., the wet year 2014-2015, while ET a values in-between these two extremes were observed in the reference year 2013-2014. In addition, the spatial variability of ET a was largest in May 2015-2016. This leads to a hypothesis, i.e., that water supply meets actual water requirements in all conditions and ET a increases in response to the ET 0 -P.
The evolution through the year 2015-2016 ( Figure 8) reflect the expected seasonality, but the higher ET 0 , which characterized this year lead to a larger winter-summer difference in ET a . Likewise, the drier conditions in 2015-2016 led to a larger spatial variability in ET a in spring and summer. Differences up to 6 mm d −1 can be noticed in ET a on DoY 157, which suggest large differences in the distribution and use of irrigation water, these results clearly require detailed data on water delivery to and within the areas served by tertiary canals for a better understanding of irrigation water management in the Doukkala.  Although no detailed data are available on the irrigation system applied in each parcel, some CGRs are dominated by surface irrigation, others by sprinkler and very few by drip irrigation. It is, therefore, possible to extract samples of the ETa time series of maps, e.g., as in Figure 8, and average the ETa values over the CGRs stratified by irrigation system. This was done separately for the reference, wet, and dry years (Figure 9). ETa remains rather similar in areas dominated by surface  Although no detailed data are available on the irrigation system applied in each parcel, some CGRs are dominated by surface irrigation, others by sprinkler and very few by drip irrigation. It is, therefore, possible to extract samples of the ETa time series of maps, e.g., as in Figure 8, and average the ETa values over the CGRs stratified by irrigation system. This was done separately for the reference, wet, and dry years (Figure 9). ETa remains rather similar in areas dominated by surface Although no detailed data are available on the irrigation system applied in each parcel, some CGRs are dominated by surface irrigation, others by sprinkler and very few by drip irrigation. It is, therefore, possible to extract samples of the ET a time series of maps, e.g., as in Figure 8, and average the ET a values over the CGRs stratified by irrigation system. This was done separately for the reference, wet, and dry years (Figure 9). ET a remains rather similar in areas dominated by surface and sprinkler irrigation, regardless of the season and the year. On the contrary, ET a was clearly lower in areas dominated by drip irrigation in both the reference and the dry year (Figure 9), while no difference was observed in the wet year. Moreover, this difference appears earlier in the dry year, which clearly suggests a better control of water supply by drip irrigation, where water is supplied to match crop transpiration, while limiting soil evaporation. In the wet season, the soil surface remains wet because of precipitation and the evaporation remains high under both surface/sprinkler and drip irrigation. In the dry period in the plots with drip irrigation the soil surface dries out and evaporation becomes rather small. Under these conditions the sum of evaporation and transpiration tends to be lower under drip irrigation than under surface or sprinkler irrigation. and sprinkler irrigation, regardless of the season and the year. On the contrary, ETa was clearly lower in areas dominated by drip irrigation in both the reference and the dry year (Figure 9), while no difference was observed in the wet year. Moreover, this difference appears earlier in the dry year, which clearly suggests a better control of water supply by drip irrigation, where water is supplied to match crop transpiration, while limiting soil evaporation. In the wet season, the soil surface remains wet because of precipitation and the evaporation remains high under both surface/sprinkler and drip irrigation. In the dry period in the plots with drip irrigation the soil surface dries out and evaporation becomes rather small. Under these conditions the sum of evaporation and transpiration tends to be lower under drip irrigation than under surface or sprinkler irrigation.

Response of Relative Evapotranspiration to Irrigation Water
As explained in Section 3.6.2 we have evaluated crop response to irrigation under varying climate forcing as measured by ET0-P and we have applied changes in the ratio ETa/ETc as a measure of response. In a similar way as done earlier (see Figure 5) to evaluate anomalies in Wd, we have evaluated relative anomalies in ETa/ETc. This evaluation has been done separately for each irrigation district.

•
Sidi Bennour: Overall, the relative anomalies in ETa/ETc and Wd ( Figure 10) are consistent, i.e., they have the same sign and comparable magnitude. In the third quarter (April-June) higher

Response of Relative Evapotranspiration to Irrigation Water
As explained in Section 3.6.2 we have evaluated crop response to irrigation under varying climate forcing as measured by ET 0 -P and we have applied changes in the ratio ET a /ET c as a measure of response. In a similar way as done earlier (see Figure 5) to evaluate anomalies in W d , we have evaluated relative anomalies in ET a /ET c . This evaluation has been done separately for each irrigation district.

•
Sidi Bennour: Overall, the relative anomalies in ET a /ET c and W d ( Figure 11) are consistent, i.e., they have the same sign and comparable magnitude. In the third quarter (April-June) higher than average W d apparently lead to higher than average ET a /ET c . To the contrary, in the second quarter (January-March) lower than average W d seem to induce lower than average (ET a /ET c ) In the first quarter (October-December) the response of W d to ET 0 -P and of ET a /ET c to W d is more variable across the reference, wet and dry year, possibly due to high precipitation in this period of the year. • Zemamra: The pattern in relative anomalies is more complex than in the case of Sidi Bennour. W d in the first and second quarter ( Figure 10) is lower than average in the wet year while it is higher than average in the third quarter. On the other hand, W d is higher than average in the first and second quarter while it is lower in the third quarter of both dry and reference year. The response of ET a /ET c to changes in W d is more variable than in the case of Sidi Bennour. • Fregh: In most cases, relative anomalies in ET a /ET c have the same sign and similar magnitude as anomalies in W d (see Figure 12). As in the case of Sidi Bennour W d was higher than average in the third (generally drier) quarter of all the three years. To the contrary, W d was lower than average in the second quarter, also in all the years. Anomalies in W d were more variable in the first quarter. Higher than average ET a /ET c was generally associated with higher than average W d , with the notable exception of the first quarter of the wet year. Lower than average ET a /ET c was clearly associated with lower than average W d in the second quarter of all years. • Gharbia: Additionally, in this district relative anomalies in ET a /ET c and W d were consistent in most cases ( Figure 13). Large negative anomalies in ET a /ET c in the first and the second quarter of the wet year and in the second quarter of the reference year were associated with lower than average W d , with the negative anomaly being rather large in the first quarter in the wet year. Similar responses were observed in the other districts particularly in Zemamra. Relative anomalies in ET a /ET c were much larger than relative anomalies in W d in the second and third quarter of the wet year, although of opposite sign in both quarters.
Water 2019, 11, x FOR PEER REVIEW 18 of 24 than average Wd apparently lead to higher than average ETa/ETc. To the contrary, in the second quarter (January-March) lower than average Wd seem to induce lower than average (ETa/ETc) In the first quarter (October-December) the response of Wd to ET0-P and of ETa/ETc to Wd is more variable across the reference, wet and dry year, possibly due to high precipitation in this period of the year. • Zemamra: The pattern in relative anomalies is more complex than in the case of Sidi Bennour. Wd in the first and second quarter ( Figure 11) is lower than average in the wet year while it is higher than average in the third quarter. On the other hand, Wd is higher than average in the first and second quarter while it is lower in the third quarter of both dry and reference year. The response of ETa/ETc to changes in Wd is more variable than in the case of Sidi Bennour. • Fregh: In most cases, relative anomalies in ETa/ETc have the same sign and similar magnitude as anomalies in Wd (see Figure 12). As in the case of Sidi Bennour Wd was higher than average in the third (generally drier) quarter of all the three years. To the contrary, Wd was lower than average in the second quarter, also in all the years. Anomalies in Wd were more variable in the first quarter. Higher than average ETa/ETc was generally associated with higher than average Wd, with the notable exception of the first quarter of the wet year. Lower than average ETa/ETc was clearly associated with lower than average Wd in the second quarter of all years.  Figure 10. Zemamra district; relative anomalies in quarterly water delivery W d and relative evapotranspiration (ET a /ET c ); 1,2,3 and R, W, D: see Figure 6.
Water 2019, 11, x FOR PEER REVIEW 18 of 24 than average Wd apparently lead to higher than average ETa/ETc. To the contrary, in the second quarter (January-March) lower than average Wd seem to induce lower than average (ETa/ETc) In the first quarter (October-December) the response of Wd to ET0-P and of ETa/ETc to Wd is more variable across the reference, wet and dry year, possibly due to high precipitation in this period of the year. • Zemamra: The pattern in relative anomalies is more complex than in the case of Sidi Bennour. Wd in the first and second quarter ( Figure 11) is lower than average in the wet year while it is higher than average in the third quarter. On the other hand, Wd is higher than average in the first and second quarter while it is lower in the third quarter of both dry and reference year. The response of ETa/ETc to changes in Wd is more variable than in the case of Sidi Bennour.  • Gharbia: Additionally, in this district relative anomalies in ETa/ETc and Wd were consistent in most cases ( Figure 13). Large negative anomalies in ETa/ETc in the first and the second quarter of the wet year and in the second quarter of the reference year were associated with lower than average Wd, with the negative anomaly being rather large in the first quarter in the wet year. Similar responses were observed in the other districts particularly in Zemamra. Relative anomalies in ETa/ETc were much larger than relative anomalies in Wd in the second and third quarter of the wet year, although of opposite sign in both quarters. Overall it should be noted that Wd is higher than average in the 3rd quarter in all districts and for all years, likewise ETa/ETc, thus suggesting a beneficial adaptation in irrigation management to the higher water requirements in late spring/early summer. As noted earlier (see also Figure 6), it appears that a relatively small increment in Wd, although much smaller than the increment in ET0-P, is sufficient to achieve a lower evapotranspiration deficit, i.e., an increment in ETa/ETc.

Inter-and Intra-Annual Climatic Variability
The analysis presented in Section 4.1 documented the large inter and intra-annual variability of ET0-P, but did not provide evidence on a real adaptive response of the water management system. Water delivery did change in wet and dry years but mainly as consequence of the reservoir operational rules in response to changes in water inflow and stock. In 1995-1996 an exceptionally large inflow led to completely fill up the reservoir with a delay of approximately one year and to slowly increasing water deliveries over a period of a few years. The subsequent dry spell led to sharply reduced water deliveries and residual water storage. These statements are supported by the observed poor correlation of Wd with ET0-P, although an analysis of actual rather than atmospheric    • Gharbia: Additionally, in this district relative anomalies in ETa/ETc and Wd were consistent in most cases ( Figure 13). Large negative anomalies in ETa/ETc in the first and the second quarter of the wet year and in the second quarter of the reference year were associated with lower than average Wd, with the negative anomaly being rather large in the first quarter in the wet year. Similar responses were observed in the other districts particularly in Zemamra. Relative anomalies in ETa/ETc were much larger than relative anomalies in Wd in the second and third quarter of the wet year, although of opposite sign in both quarters. Overall it should be noted that Wd is higher than average in the 3rd quarter in all districts and for all years, likewise ETa/ETc, thus suggesting a beneficial adaptation in irrigation management to the higher water requirements in late spring/early summer. As noted earlier (see also Figure 6), it appears that a relatively small increment in Wd, although much smaller than the increment in ET0-P, is sufficient to achieve a lower evapotranspiration deficit, i.e., an increment in ETa/ETc.

Inter-and Intra-Annual Climatic Variability
The analysis presented in Section 4.1 documented the large inter and intra-annual variability of ET0-P, but did not provide evidence on a real adaptive response of the water management system. Water delivery did change in wet and dry years but mainly as consequence of the reservoir operational rules in response to changes in water inflow and stock. In 1995-1996 an exceptionally large inflow led to completely fill up the reservoir with a delay of approximately one year and to slowly increasing water deliveries over a period of a few years. The subsequent dry spell led to sharply reduced water deliveries and residual water storage. These statements are supported by the observed poor correlation of Wd with ET0-P, although an analysis of actual rather than atmospheric  Figure 13. Gharbia district; relative anomalies in quarterly water delivery W d and relative evapotranspiration ET a /ET c ; 1,2,3 and R, W, D: see Figure 6.
Overall it should be noted that W d is higher than average in the 3rd quarter in all districts and for all years, likewise ET a /ET c , thus suggesting a beneficial adaptation in irrigation management to the higher water requirements in late spring/early summer. As noted earlier (see also Figure 6), it appears that a relatively small increment in W d , although much smaller than the increment in ET 0 -P, is sufficient to achieve a lower evapotranspiration deficit, i.e., an increment in ET a /ET c .

Inter-and Intra-Annual Climatic Variability
The analysis presented in Section 4.1 documented the large inter and intra-annual variability of ET 0 -P, but did not provide evidence on a real adaptive response of the water management system. Water delivery did change in wet and dry years but mainly as consequence of the reservoir operational rules in response to changes in water inflow and stock. In 1995-1996 an exceptionally large inflow led to completely fill up the reservoir with a delay of approximately one year and to slowly increasing water deliveries over a period of a few years. The subsequent dry spell led to sharply reduced water deliveries and residual water storage. These statements are supported by the observed poor correlation of W d with ET 0 -P, although an analysis of actual rather than atmospheric water demand would be needed for a robust assessment. In the period 1992-2016 there was just one dry spell, although rather long from 2000 to 2009. From the point of view of adaptive water management, it would be rather relevant to analyze a much longer climate record to estimate the frequency of such events, in order to develop an adaptation strategy either by increasing water storage capacity or by controlling land use to balance water supply and actual water demand. These findings are in line with the review on the evolution of the irrigation management system in Morocco presented by [54], which emphasized the urgent need for renewed efforts by the ORMVA-s towards better implementation and performance of participatory irrigation management. Particularly, this author emphasized that irrigation water conservation implies the adoption of measures to modify water demands and maximize efficiency in water use.
Likewise, the detailed analysis of actual water demand and water use versus water supply in extreme wet and dry years did not provide evidence of seasonal adaptation of water delivery. Again, the operational rules of ORMVAD promote short term adaptation, e.g., on a weekly basis, to water availability (precipitation), rather than to the seasonal evolution of land use and crop water use. To further support this statement, we note that no systematic relation was observed between ET a /ET c and W d under the extreme wet (2014-2015) and the extreme dry (2015-2016) conditions with very different ET 0 -P.
Our analysis of the response of ET a /ET c to W d for all the CGRs did provide evidence of a larger variability in such response in the CGRs where drip irrigation is dominant. This suggests larger differences in the efficiency of water application notwithstanding the higher water application efficiency of drip irrigation [55].
It should be taken into account that we have analyzed the response of W d to ET 0 -P using seasonal data and the response of ET a /ET c to W d using quarterly data, this approach may have hidden the subtle aspects of the response of crop water use to short term adaptation in water delivery in relation with precipitation events.

Response to Climatic Variability
Sidi Bennour ( Figure 11): Relative changes in water delivery and relative evapotranspiration ET a /ET c were comparable for the first, second, and third quarters in the reference year, in the second and third quarters in the wet year and in the third quarter in the dry year. In the first quarter of the wet year the relative anomalies in ET a /ET c were much smaller than relative anomalies in water delivery. In this case water delivery was reduced taking into account the very large precipitation (P) in same quarter. Specifically, the precipitation in this quarter ( Table 4) was 226% of the average over all quarters and all years. The large precipitation was sufficient to maintain high ET a /ET c . In the 2nd quarter of the dry year, the water delivery was not increased to compensate the larger than average ET 0 -P and was even slightly lower than the average. As a consequence, ET a /ET c was equal to 0.66, i.e., much lower than the average over all quarters and all years. In the first quarter of the dry year, ET 0 -P was much larger than in the same quarter of the other years. Water deliveries were increased by 18% above the average, which apparently was sufficient to maintain the ET a /ET c roughly constant. This is what adaptive irrigation water management should achieve. Zemamra ( Figure 10): In the 3rd quarter of the wet year ET a /ET c was much larger than the average. This is explained by the water deliveries which were 30% above the average in the same quarter. In addition, the precipitation in the first quarter of the same year was 226% higher than the average, thus fully replenishing the soil water reservoir to maintain high ET a /ET c . In the 1st quarter of the dry year water deliveries were increased by 30% over the average, which was sufficient to maintain ET a /ET c roughly constant. In the 2nd quarter of the dry year, water deliveries were 16% higher than the average while ET a /ET c was 27% lower than the average. This can be explained by the precipitation in the first quarter being 29% lower than average. A similar response in ET a /ET c was observed in the 1st quarter of the reference year when deliveries were 10% higher than average, but the precipitation was 65% lower than the average. Faregh ( Figure 12): The relative changes in deliveries and ET a /ET c were roughly similar in the 1st, 2nd and 3rd quarter of the reference year, in the 2nd and 3rd quarter of the wet year and the 3rd quarter of the dry year. The 1st quarter of the wet year and the 3rd quarter of the reference year, water deliveries were adapted successfully to maintain ET a /ET c roughly constant notwithstanding the large variations in ET 0 -P, which was 134% lower than average in the 1st case and 83% higher than average in the second case. In the 2nd quarter of the dry year water deliveries were correctly reduced in response to ET 0 -P being 30% lower than average with precipitation being 16% higher than average but concentrated in few days in February. This led to a 32% decrease in ET a /ET c . Gharbia ( Figure 13): The relative changes in ET a /ET c and water deliveries had a similar magnitude in the 1st, 2nd and the third quarter of the wet year, in the 1st and 3rd quarter of the dry year and in the 2nd quarter of the reference. In the 2nd quarter of the dry year water deliveries were 16% higher than average, notwithstanding ET 0 -P being 30% lower than average with precipitation being 16% higher than average but concentrated in few days in February. This led to a 29% decrease in ET a /ET c .

Conclusions
This study evaluated whether the rigid water deliveries rules applied by ORMVAD, the water management authority of the Doukkala irrigation scheme, combined with the autonomous decision by farmers at each rotational interval, would lead to actual delivery of irrigation water consistent with an assessment of demand based on established agro-hydrological practices. Water delivery did respond to inter-annual variability in hydrological conditions, but rather to water inflows to reservoirs than to atmospheric water demand ET 0 -P. This notwithstanding, we have shown several instances where anomalies in water delivery do correlate with anomalies in ET 0 -P, which is surprising to some extent, since no such quantitative evaluation of water requirements plays a role in the decision by farmers about requesting water delivery.
Our evaluation of the crop response to variable water delivery showed that in 70%-90% of the cases the ratio ET a /ET c changes with water delivery as expected, i.e., it increases with larger water deliveries and vice versa. There are a few cases, however, where anomalies in ET a /ET c and W d have the opposite sign. Most likely this is due to the distribution in the time of precipitation and water delivery, but we do not have sufficiently detailed data to evaluate this hypothesis.
If we assume that the decision by farmers about whether or not to request water at a given time is based on a perception of ET 0 -P, this conclusion raises a new research question, i.e., whether farmers perceive correctly an impending water shortage and take the right decision about requesting delivery of irrigation water in a given rotational interval.