Monitoring of Irrigation Schemes by Remote Sensing: Phenology versus Retrieval of Biophysical Variables

Abstract

The appraisal of crop water requirements (CWR) is crucial for the management of water resources, especially in arid and semi-arid regions where irrigation represents the largest consumer of water, such as the Doukkala area, western Morocco. Simple and (semi) empirical approaches have been applied to estimate CWR: the first one is called K_c-NDVI method, based on the correlation between the Normalized Difference Vegetation Index (NDVI) and the crop coefficient (K_c); the second one is the analytical approach based on the direct application of the Penman-Monteith equation with reflectance-based estimates of canopy biophysical variables, such as surface albedo (r), leaf area index (LAI) and crop height (h_c). A time series of high spatial resolution RapidEye (REIS), SPOT4 (HRVIR1) and Landsat 8 (OLI) images acquired during the 2012/2013 agricultural season has been used to assess the spatial and temporal variability of crop evapotranspiration ET_c and biophysical variables. The validation using the dual crop coefficient approach (K_cb) showed that the satellite-based estimates of daily ET_c were in good agreement with ground-based ET_c, i.e., R² = 0.75 and RMSE = 0.79 versus R² = 0.73 and RMSE = 0.89 for the K_c-NDVI, respective of the analytical approach. The assessment of irrigation performance in terms of adequacy between water requirements and allocations showed that CWR were much larger than allocated surface water for the entire area, with this difference being small at the beginning of the growing season. Even smaller differences were observed between surface water allocations and Irrigation Water Requirements (IWR) throughout the irrigation season. Finally, surface water allocations were rather close to Net Irrigation Water Requirements (NIWR).

1. Introduction

The interest for the assessment of irrigation performance using satellite data developed in the late 1980s due to growing consensus on the difficulty of collecting the required ground data continuously, reliably and in a consistent way across all major irrigation schemes worldwide [1–4]. Since in a large part of irrigated lands water is allocated proportionally to irrigated area, work was initially focused on the relation between allocated water and irrigated area, observable with multispectral satellite data [5]. Later on, other aspects of irrigation water management were evaluated, like crop water requirements [6–12], actual consumptive water use [13–18], water productivity [17,19–22] and water and salinity stress [23–29]. On a higher level of abstraction, irrigation performance may be evaluated for different objectives such as equity, adequacy, or effectiveness [30–32]. In this study, we focus on the evaluation of adequacy by relating water allocation to water requirement.

The most common and practical approach used for estimating crop water requirements (CWR) is the FAO-56 method (the Food and agricultural Organization of United States (FAO) Irrigation and Drainage Paper No. 56) [33], based on the combination of reference evapotranspiration ET₀ and crop coefficients (K_c) to determine crop evapotranspiration (ET_c) under unrestricted water availability. In the majority of the studies, the K_c values are obtained by the single crop coefficient approach, where crop transpiration and soil evaporation are combined into a single K_c coefficient. Infrequently, the dual crop coefficient (k_cb) approach is used, where the effects of crop transpiration and soil evaporation are determined independently [34–38]. The FAO-56 method is based on the use of crop specific parameters. While this is not an issue for on-farm evaluations of CWR, it becomes rather challenging when the objective is to monitor CWR of large irrigation schemes. In this paper, we have adopted the definitions CWR = ET_c and irrigation water requirement IWR = ET_c – P, where P is precipitation.

Frequent mapping of crop types is, in principle, feasible using multispectral and multi-temporal satellite data [39–47], but accurate classification requires, in most cases, ground reference data and analysis of images acquired at multiple dates [48]. This makes timely availability of crop maps rather unlikely, thus reducing the timeliness and the reliability of a CWR monitoring service based on satellite data.

Remote sensing methods which do not require knowledge of crop type to determine CWR have been developed taking advantage of the strong physical relationship between the spectral response of cropped surfaces and the corresponding values of CWR and K_c. Examples of these approaches can be found in [11,48–51]. To such end, empirical relationships have been found to retrieve the value of K_c or K_cb from simple vegetation indexes, i.e., NDVI (Normalized difference Vegetation Index) [52].

Several indices have been proposed as alternatives to NDVI to estimate Kc, such as the Perpendicular Vegetation Index PVI [53], the Soil-Adjusted Vegetation Index SAVI [54–57], the Weighted Difference Vegetation Index WDVI [58] and the Global Environment Monitoring Index GEMI [59]. These indices have been formulated in order to reduce the influence of perturbing effects such as the soil background or the atmospheric influence, which may alter significantly the reflectance of vegetated surfaces.

A different approach is by using directly the Penman-Monteith equation either to determine a generic K_c (crop identification not required) or ET_c [6,10,60,61]. This approach is based on retrieving from satellite data the crop properties which determine ET_c, i.e., crop height (h_c), surface albedo (α), and leaf area index (LAI). The crop height influences the aerodynamic resistance (r_a) term of the FAO Penman-Monteith equation and the turbulent transfer of vapor from the crop into the atmosphere. The r_a term appears twice in the full form of the FAO Penman-Monteith equation. The surface albedo of the crop-soil surface influences the net radiation at the surface, Rn, which is the primary forcing factor of the transpiration and evaporation processes. The surface albedo (α) is affected by the fraction of ground covered by vegetation and by the soil surface wetness. The canopy resistance of the crop to vapor transfer is affected by leaf area (number of stomata), leaf age and condition, and leaf-level stomatal control. The fraction of exposed soil also affects K_c.

In this study, we have applied and evaluated two methods to determine CWR, i.e., the one based on the correlation between K_c and NDVI (K_c-NDVI method) and the one using directly the Penman-Monteith equation (analytical method). We have applied both methods to evaluate the adequacy of water allocation in the Doukkala irrigation scheme.

In the Doukkala area (Western Morocco), water demand has significantly increased over the last decades while fresh water resources are becoming increasingly scarce. This is mainly due to the combined effect of climate change, persistent drought and the increase of water demand related to increase in irrigated area, urbanization and recreational projects. This shows the necessity to use available water resources as effectively as possible in order to avoid or at least mitigate the consequences of recurring droughts. This is particularly important for water management in agricultural areas, where irrigation represents the biggest water consumer.

This paper is organized as follows. After the description of the study area and the data collected (Section 3), we describe in detail our implementation of the two methods (Section 4), followed by the presentation of results, including the evaluation against in-situ observations, and concluding with global interpretation of our results.

2. Study Area

The Doukkala region lies in western Morocco, between 32°15′N and 33°15′N latitude and 7° 55′W and 9°20′W longitude in the downstream portion of the hydraulic basin of Oum Er_Rabia (Figure 1). Geo-morphologically speaking, Doukkala is divided into three parts: the coastal area, the Sahel and the plain. The latter extends over an area of 3500 km² and it is located at about 120–130 m above sea level [62], with favorable conditions for agricultural development as regards arable land and soil fertility.

The irrigated area of Doukkala is among the largest (96,000 ha) and earliest developed areas in Morocco, remarkable for its size and strategic importance for national production, specially sugar beet (38%) and commercialized milk (20%). The important crops grown in the study area include wheat, corn, sugar beet, and alfalfa [63].

The climate of Doukkala is typically semi-arid to temperate and mild in winter whereas in summer it is generally warm and dry, with a large inter and intra-annual variability of rainfall, which amounts to 316 mm/year on average (1964–2009). Reference evapotranspiration (ET₀) is 1434 mm/year on average 2000–2008.

The irrigation is performed by means of different techniques with different efficiency (surface, sprinkler and drip irrigation) with a rotational interval of 15 days. Surface irrigation is the dominant system. Water allocation is calculated on the basis of current irrigated area, as determined on the basis of requests by farmers for each rotational interval. The gross irrigation water depth allocated to farmers is 864 m³/ha by rotational interval, equivalent to 1728 mm/month, which is largely sufficient to meet the water needs of the dominant crops. The water resources mobilized for irrigation come mainly from the dam Al Massira, a major water storage structure in the watershed of Oum Er-Rbia with a capacity of approximately 2760 × 10⁶ m³.

The Doukkala irrigation scheme is divided into the High and Low Sections. We focused on the Low Section, which contains three main districts: Faregh, Sidi Bennour and Zemamra, respectively, from the East to the West. Each district is divided into a number of Centers of Irrigation Management (CGR) irrigated with different irrigation systems (

Table 1. Irrigation Systems in different Centers of Irrigation Management (CGRs) in three districts of the Doukkala Irrigation Scheme (ORMVAD-personal communication).

District	CGR	Area (Ha)	Irrigation System
SidiBennour	330	5305.25	Surface Irrigation
	331	3520.06
	333	4293.39
	335	3202.49
	336	4197.1
	337	3112.55
	338	1738.92
	338 sprinkler	1905.3	Sprinkler Irrigation
Zemamra	320	2995.18	Sprinkler Irrigation
	321	5327.4
	322	3243.68
	324	4565.71
	325	3122.27
Faregh	312	4840.24	Drip Irrigation
	332	1490.46	Surface Irrigation
	310	4468.60
	311	5021.25

). The Doukkala irrigation scheme is managed by the Regional Office of Agricultural Development in Doukkala (ORMVAD).

3. Data Collection

3.1. Satellite Data

To meet the combined requirements of high spatial resolution and frequent imaging to monitor crop development and water requirements, we have combined image data acquired by multiples sensors at high and very high spatial resolution: (Landsat 8 (30 m), Spot 4 HRVIR1 (20 m) and Rapideye (5 m)) during the irrigation season 2012/2013 (

Table 2. Overview of image data characteristics and coverage of the study area.

SENSOR	DATE	Area	SPECTRAL RESOLUTION (μm)	Spatial Resolution	ORBIT
SPOT4-HRVIR1	From January to June 2013	FAREGH	XS1: 0.500–0.590	20 m	Altitude:832 km revisit: 5 days
			XS2: 0.610–0.680
			XS3: 0.790–0.890
			SWIR (HRVIR): 1.530–1.750
RapidEye-REIS	10 December 2012	ZEMAMRA SIDIBENNOUR	Blue: 0.440–0.510	5 m	Altitude:630 km revisit: Daily (off-nadir); 5.5 days (at nadir)
			Green: 0.520–0.590
	8 February 2013		Red: 0.630–0.685
			Red-Edge: 0.690–0.730
			NIR: 0.760–0.850
Landsat 8-OLI	19 April 2013	SIDIBENNOUR ZEMAMRA	Coastal/Aerosol: 0.433–0.453	30 m	Altitude:705 km revisit: 16 days
			Blue: 0.450–0.515
	26 April 2013		Green: 0.525–0.600
			Red: 0.630–0.680
	13 June 2013		NearInfrared: 0.845–0.885
			SWIR: 1.560–1.660
	29 June 2013		SWIR: 2.100–2.300
			Cirrus: 1.360–1.390
	15 July 2013		Panchromatic: 0.500–0.680	15 m

). All satellite images acquired were corrected geometrically with the following system of coordinate: UTM, WGS-84, zone 29.

An unprecedented time series consisting in Spot4 data (level 2A) with a five-day revisit interval has been acquired and analyzed. The Level 2A product is in-band surface reflectance corrected from atmospheric effects using the SMAC model [64].

3.2. Meteorological and Water Flow Data

Daily meteorological data (i.e., Air temperature (T_a), relative humidity (RH), solar radiation (R_s), wind speed (U), precipitation (R) and reference evapotranspiration ET₀) on the same dates of the satellite observations were collected at the Zemamra and Khmiss Mettouh meteorological stations and provided by ORMVAD. The observations at Zemamra were used for the district of Zemamra and the district of Sidi Bennour and the observations of Khmiss Mettouh were used for the district of Faregh.

Data on monthly water allocation by CGR were provided by ORMVAD. The water allocation is calculated monthly for each tertiary and secondary irrigation unit on the basis of requests submitted by farmers for each rotational interval. The irrigated area of each farm is multiplied by the irrigation module (i.e., 864 m³/ha) to determine the duration of water delivery to each farm and the total volume to be delivered. These farm-level irrigation water depths were added up to obtain the monthly water allocation by the CGRs we have used in our assessments on the adequacy of water allocations.

3.3. Data In-Situ

Periodical field campaigns were carried out with multiple objectives: to explore the structure of the irrigation system over the area as well as the irrigation systems (surface irrigation, sprinkler irrigation and drip irrigation); to collect phenological data on the dominant crops and to collect data on fractional cover and crop height.

The fractional cover (f_c) and crop height (h_c) of the dominant crops (wheat, sugar beet, alfalfa, berseem and corn) were measured in 22 pilot plots in the Zemamra and Sidi Bennour districts on 17–19 December, 27–28 February, 3–4 April, 23–24 May, and 20 July in the growing season 2012/2013. During the May and July surveys, 14 additional plots in the Faregh district were sampled.

For each plot (max. 1 ha), we measured h_c in several points (5–6 points) with a graduated stick, and used the mean value of h_c as representative of the plot. In the same plot, we estimated f_c at 5–6 locations using an approximate target of 1 m², and took the mean value of f_c as representative of the plot. Then, we calculated in-situ crop evapotranspiration ET_c using the ground measurement of h_c and f_c combined with meteorological data (see Section 4.3.2 Dual crop coefficient approach). For each plot where we had such in-situ ET_c, we obtained the corresponding satellite ET_c by sampling image data as described in Section 5.3.3.

In some cases, the difference between the dates of satellite data acquisition and field work date was significant. In these cases, we interpolated the temporal observations of h_c and f_c linearly to obtain estimates on the dates of acquisition of the satellite data. This gave 80 pairs of ground and satellite observations throughout the growing season.

4. Methods

4.1. Work Flow

Crop evapotranspiration ET_c is the basic information for the evaluation of crop water requirements and irrigation management. In this work, we estimated ET_c with two different methods: FAO-56 (k_c—NDVI) and analytical approach. After the pre-processing of satellite images, ET_c was estimated for the entire study area and the 2012–2013 growing season and validated using in-situ observations (Figure 2). The FAO-56 method is the most widely used method to compute ET_c and is based on the estimation of the so-called “crop coefficient” K_c, defined as the ratio of total evapotranspiration ET_c by reference evapotranspiration ET₀. We evaluated our RS estimates of ET_c using the FAO-56 dual approach to calculate the basal crop coefficient k_cb and the evaporation coefficient k_e. with the f_c and h_c observed in the field during the above mentioned surveys.

The analytical approach is based on the direct application of the Penman-Monteith equation. The required vegetation properties are surface albedo (r), Leaf Area Index (LAI) and the crop height(h_c) that are obtained from the processing of E.O data. LAI is used to compute canopy resistance (r_c,min),_,crop height to calculate the aerodynamic resistance (r_a) and surface albedo is used to calculate the net shortwave radiation (R_ns).

4.2. Pre-Processing

The study area is rather flat at an elevation of 150 m ± 15 m, thus we did not carry out an additional geometric correction taking into account the topography.

The radiometric calibration of satellite images was achieved in two steps:

Firstly, by conversion of DN values into radiance; the following Equation (1) was used for conversion:

$$L_{\lambda }=(X/A)+B$$

(1)

where, L is Top Of Atmosphere (TOA) band spectral radiance observed by the satellite. X is a Digital Number. A is the absolute calibration gain for relevant spectral band and B is the absolute calibration bias for relevant spectral band.

Secondly, the atmospheric correction was performed on Landsat 8 and RapidEye data using the FLAASH model that incorporates the MODTRAN 4 model, and the input parameters used in this study are presented in

Table 3. MODTRAN input parameters used in this study.

PARAMETER	VALUE	Date
Model Atmosphere	Mid Latitude Summer	10 December 2012
		8 February 2013
		19 April 2013
		26 April 2013
	Tropical	13 June 2013
		29 June 2013
		15 July 2013
Aerosol Model	Rural
Aerosol Retrieval	None
Visibility	40 km (Default)
Ground Altitude Above Sea Level	150 m

. The model calculates the Top Of Canopy (TOC) in-band Lambertian reflectance ρ_p by means of a radiance-to-reflectance conversion, using the Equation (2):

$${\rho }_p=\frac{\pi ·L_{\lambda }·d^2}{E_{\lambda }^0·\text{cos\hspace{0.17em}}{\theta }_s}$$

(2)

where:

L_λ: TOA in-band spectral radiance observed by the satellite; θ_s: Solar zenith angle; $E_{\lambda }^0$ is the mean in-band solar exo-atmospheric constant; d is the Earth-Sun distance at sensor’s aperture, in astronomical units; d = 1 – 0.01674 cos(0.9856 (JD-4)), where JD is Julian Day.

For each band (VIS, NIR and SWIR), a time series of the pre-processed multispectral reflectance data was constructed using the nearest neighbor resampling to match the spatial resolution of the master image (UTM-WGS84). This layer stack was then exported as a multilayer GeoTIFF file, which is easily read and analyzed by Matlab.

To evaluate the consistency of the multi-sensor reflectance data (Table 2), we have simulated [65] the Top of Canopy (TOC) spectral reflectance (400–2400 nm with 1 nm spectral resolution) for a very heterogeneous soil vegetation scene with widespread irrigation. We constructed a data set including 60,000 samples. For each sample, the 2101 spectral bands were convolved with the spectral response functions of the three sensors OLI (Landsat8), HRVIR1 (Spot4) and REIS (RapidEye) to simulate the TOC reflectance observed by each sensor. The red and NIR reflectances were used to calculate the NDVI and WDVI for each sensor and for each sample. The OLI sensor was taken as the reference and the other two sensors (REIS and HRVIR1) were compared with OLI.

For NDVI, the RMSE values were respectively 0.0469, 0.0328 and 0.02 for the pairs OLI-HRVIR1, OLI-REIS and HRVIR1-REIS, while for WDVI the RMSE values were 0.0288, 0.0343 and 0.006 for the same sensor combinations. Accordingly, we neglected these differences and concluded that OLI, HRVIR1 and REIS gave consistent observations of WDVI and NDVI.

4.3. Application of FAO-56 Model

In 1998, FAO proposed the FAO-56 Penman-Monteith reference evapotranspiration (ET₀) for irrigation scheduling [33]. This method has been widely used because it gives satisfactory results under various climate conditions across the world [66–69].

According to this model, two parameters are required to estimate the ET_c: the crop coefficient kc and ET₀. Crop coefficient curves provide simple, reproducible means to estimate ET_c from weather-based ET₀ values. ET₀ is defined as the evapotranspiration of a reference grass, completely covering the soil, well-watered and actively growing under optimal agronomic conditions.

In FAO-56, two approaches to determine k_c are presented: the single crop coefficient approach (k_c), which we have applied with remote sensing data, and the dual crop coefficient approach (k_cb), which we applied only with ground measurement.

4.3.1. Single Crop Coefficient Approach: k_c-NDVI Method

In the single crop coefficient (K_c), the effect of crop transpiration and soil evaporation is combined into a single K_c. The remotely sensed spectral reflectance data can be used to estimate K_c, because both K_c and spectral vegetation indices are correlated to leaf area index and fractional ground cover [70]. The simplest approach to derive K_c from remotely sensed data uses a linear relationship between K_c and NDVI (Normalized Difference Vegetation Index). NDVI is obtained from red (R) and near infrared (NIR) reflectance (0.6–0.7 μm and 0.7–1.3 μm, respectively), which are present in most imaging radiometers (Equation (3)).

$$NDVI=\frac{NIR-R}{NIR+R}$$

(3)

This approach was introduced by [71] and used and validated in further case studies by [72–75]. The theoretical basis has been established by [70]. This approach is one of the most promising ones for operational applications [8].

We have further simplified this concept, following [49,50] by using their K_c-NDVI relationship (Equation (4)). This is a relationship between the maximum NDVI (set as 0.8) and the maximum K_c (1.2 at effective full cover) and the minimum (bare soil) NDVI (0.16) and bare soil K_c (0.4), respectively. These values are valid for NDVI calculated from in-band surface reflectance and they are not crop-dependent.

$$K_c=1.25\times NDVI+0.2$$

(4)

4.3.2. Dual Crop Coefficient Approach

The dual crop coefficient approach of FAO-56 is intended to improve daily estimates of ET_c by considering separately the contribution of soil evaporation (k_e) and crop transpiration (k_cb). The dual method utilizes “basal” crop coefficients (k_cb) representing ET of a crop interspersed with dry soil, where:

$$k_{\text{c}}=k_{\text{cb}}+k_{\text{e}}$$

(5)

As crops grow, the crop height and the leaf area change, and due to the differences in evapotranspiration during the various growth stages, the k_c for a given crop will vary over each period. Following the FAO-56 approach (page 187–189) [33], growth season of the crop is divided into four distinct growth stages: initial, crop development, mid-season stage and late season. The k_cb mid-season can be estimated from simple field observations and measurement of fractional cover (f_c) and crop height (h_c):

$${\text{K}}_{\text{cb}-\text{mid}}={\text{K}}_{\text{c}-\text{min}}+\left({\text{K}}_{\text{cb}-\text{full}}-{\text{K}}_{\text{c}-\text{min}}\right)\left(\text{min}\left(1,2{\text{f}}_{\text{c}},{\left({\text{f}}_{\text{c}-\text{eff}}\right)}^{\frac{1}{1+hc}}\right)\right)$$

(6)

where:

• K_cb-mid is the estimated basal K_cb during the mid-season when plant density and/or leaf area are lower than for full cover conditions;
• K_cb-full is the estimated basal K_cb during the mid-season (at peak plant size or height);
• K_c-min is the minimum K_c for bare soil (in the presence of vegetation) (K_cmin ≈ 0.15–0.20),
• f_c is the observed fraction of soil surface that is covered by vegetation as observed from nadir [0.01–1], f_c-eff is the effective fraction of soil surface covered or shaded by vegetation [0.01–1].
• h_c is the plant height (m).

In Equation (6), we have applied the K_c-min and K_c-full values given by FAO-56 approach [33] (Table 17, page 137).

The soil evaporation coefficient (K_e) can be estimated from (1 – f_c) using an empirical relationship given by [76] and applied in [15] for irrigated wheat field in Morocco:

$$K_e=K_{e\hspace{0.17em}\text{max}}\times (1-f_c)$$

(7)

It gives the maximum soil evaporation coefficient (K_{e max}) when the soil is bare (f_c ≈ 0) and = 0 when the vegetation attains full cover (f_c ≈ 1). We adopted K_emax = 0.25 according to [15].

The choice of this value is not random but is based on the frequency and quantity of water supply and the rate of the reference evapotranspiration ET₀. We have estimated K_{e max} = K_{c ini} at low f_c, taking into account the frequency of irrigation and ET₀. In our case, the frequency of irrigation is 15 days during the growing season and the mean value of ET₀ is 3.5–4 mm per day.

Finally, we have applied the procedure described above to calculate ET_c for all reference plots and dates using our ground measurements of fractional cover (f_c) and crop height (h_c) and FAO-56 guide [33] (Table 17, page 137). We have used this ground based ET_c data set to evaluate the remote sensing estimates of ET_c obtained with the K_c-NDVI and the analytical methods described below.

4.4. Analytical Approach

Analytical approach is based on the direct application of the Penman-Monteith Equation (8) using crop characteristics estimated from satellite images, in analogy to the direct calculation proposed by FAO [10]:

$$E{T}_c=\frac{1}{\lambda }\frac{\mathrm{\Delta }(R_{ns}-R_{nl}-G)+1.013\rho De/r_a}{\mathrm{\Delta }+\gamma (1+r_{c,\text{min}}/r_a)}$$

(8)

where:

λ is the latent heat of vaporization [MT/kg]; R_ns is the net short wave radiation (MJ/m²·d); R_nl is the net long wave radiation (MJ/m²·d); G is the soil heat flux (KJ/m²·s), D_e is the vapor pressure deficit of the air (KPa); ρ is the mean air density at constant pressure (kg/m³); γ is the psychometric constant (KPa/°C); Δ is the slope of the saturation vapor pressure temperature relationship (KPa/°C); r_c,min and r_a are the minimum surface (in the absence of water stress) respective of the aerodynamic resistance.

Note that Equation (8) gives an estimate of K_c by dividing ET_c by ET₀.

In Equation (8), we have a radiative and an aerodynamic term. The former is the net shortwave radiation, R_ns, while the latter accounts for turbulent transport of heat and vapor.

The crop resistances r_c,min and aerodynamic resistance r_a require the knowledge of canopy geometrical characteristics.

The surface albedo (r), the Leaf Area Index (LAI) and crop height (h_c) can be integrated in the Equation (8) as follows:

$$R_{ns}=(1-r)\times R_s$$

(9)

$$r_{c,\text{min}}=\frac{100}{0.5\times LAI}$$

(10)

$$r_a=\frac{\text{ln}(\frac{Z_U-(2/3)h_c}{0.123h_c})\hspace{0.17em}\text{ln}(\frac{Z_h-(2/3)h_c}{0.0123h_c})}{0.168U_z}$$

(11)

where: R_s is the total incoming solar radiation (MJ/m²·d); Z_U and Z_h are the measurement heights for wind and humidity, respectively (m); h_c is the crop height (m) and U is wind speed at height z (m/s).

To determine the aerodynamic resistance (r_a) with Equation (11), we have used the Equation (4) page 20 in FAO-56 [33] where we have set the zero plane displacement height (d = 2/3 h_c), the roughness length for momentum Z_0m = 0.123 h_c and the roughness length for heat Z_0h = 0.0123 h_c, since we are dealing with full homogeneous vegetation canopies.

The surface albedo (r) is the spectrally integrated hemispherical solar reflectance and is the driving variable of the radiation budget of a surface. The estimation of (r) can be done using measurements of the reflected solar radiance K^↑(ϑ, Φ, λ) (Wm⁻²·sr⁻¹) at a wavelength λ (nm) and can be expressed as a function of viewing zenith, ϑ, and azimuth, Φ angles, respectively (Equation (12)). However, the current sensor capabilities impose several simplifications. In the first instance, the observed surface is considered as Lambertian. In this case, the dependence of K^↑ on ϑ and Φ will be neglected and r can be estimated from any direction of observation, by means of the Equation (13), using the reflectance corrected values for atmospheric effects, ρ_λ and the weighting coefficient W_λ [77]

The weighting coefficients calculated from the extraterrestrial solar irradiance E°_λ for each band for RapidEye (REIS), Landsat8 (OLI) and SPOT4 (HRVIR1) used in our study area are summarized in

Table 4. Weighting coefficients for the calculation of albedo α by using Equation (13) for different sensors.

Sensor	Spectral Band (μm)	Weighting Coefficient Wλ
RapidEye (REIS)	Blue: 0.440–0.510	0.2455
	Green: 0.520–0.590	0.2989
	Red: 0.630–0.685	0.1973
	NIR: 0.760–0.850	0.2583
Landsat 8 (OLI)	Blue: 0.450–0.515	0.2935
	Green: 0.525–0.600	0.2738
	Red: 0.630–0.680	0.233
	NIR: 0.845–0.885	0.1554
	SWIR: 1.560–1.660	0.0322
	SWIR: 2.100–2.300	0.0121
Spot 4 (HRVIR1)	XS1: 0.500–0.590	0.3925
	XS2: 0.610–0.680	0.3339
	XS3: 0.790–0.890	0.224
	SWIR:1.530–1.750	0.0496

.

$$r=\underset{300}{\overset{2500}{\int }}\frac{\left[\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\pi /2}{\int }}{K}^{\uparrow }(\vartheta ,\mathrm{\Phi },\lambda )\hspace{0.17em}\text{cos\hspace{0.17em}}\vartheta \hspace{0.17em}\text{sin\hspace{0.17em}}\vartheta d\vartheta d\mathrm{\Phi }\right]}{K^{\downarrow }(\lambda )}d\lambda$$

(12)

$$r={\sum }_{\lambda }{w}_{\lambda }{\rho }_{\lambda }$$

(13)

The Leaf Area Index (LAI) quantifies the amount of foliage area per unit ground surface area, and is an important structural property of vegetation canopies [58]. It is a crucial variable controlling many biological and physical processes associated with vegetation on the earth’s surface, such as photosynthesis, respiration, transpiration, carbon and nutrient cycle, and rainfall interception. In an operational context, the estimation of LAI from measurements of spectral reflectance has been mostly based on the (semi) empirical relationships between this parameter and vegetation indices.

The Weighted Difference Vegetation Index (WDVI) (Equation (14)) has the advantage to reduce to a great extent the influence of soil background on the spectral signal [58], by means of the factor C (Equation (15)). The soil line slope (C) represents a linear relationship between red and NIR reflectance of bare soil, and accounts for the effects of the soil background on the vegetation index, and depends on soil type, soil texture and soil moisture.

$$WDVI=r_{ir}-C·r_r$$

(14)

$$C=r_{s,ir}/r_{s,r}$$

(15)

where r_ir is the total measured NIR reflectance, r_r is the total measured red reflectance; The r_s,ir is the NIR reflectance of the bare soil, and r_s,r is the red reflectance of bare soil.

Once WDVI and WDVI_∞ (representing the asymptotically limiting value for WDVI when LAI tends to infinity) are determined, a light extinction coefficient α* has to be estimated in order to determine the LAI through the Equation (16):

$$\text{LAI}=(-1/{\alpha }^{*})\text{ln}(1-\text{WDVI}/{\text{WDVI}}_{\infty })$$

(16)

where α* represents the light extinction through the vegetation canopy, while it is dependent on crop geometry and solar zenith angle. We have used the average value (α* = 0.37) established by [78] from field measurements of LAI and WDVI for different crops.

An accurate estimation of crop height (h_c) using spectral reflectance data is quite difficult. Several studies were conducted in this framework using airborne laser altimeter [79–81], and other studies, i.e., [82] using logarithmic relationships between h_c and different vegetation indices (SAVI, WDVI, NDVI, TVI, etc.). The same author evaluated these relationships against field measurements of alfalfa and grass. Since h_c can be estimated indirectly from the canopy roughness length (Z_0m), [83] Brutsaert proposed a formula using (Z_0m) (Equation (17)). Several relationships between NDVI-Z_0m have been proposed by [84–86]. We have tested some of these equations and finally we have chosen the Equation (18) [84]. The values of the C2 and C3 coefficients have been determined by comparing estimates with Equation (18) and the C2, C3 values given by [84]. The values C2 = (−5.2) and C3 = 5.3 gave the best agreement with our observations.

$$h_c=Z_{0m}/0.123$$

(17)

$$Z_{0m}=\text{exp}(C2+C3\times NDVI)$$

(18)

4.5. Irrigation Performance Indicators

Several researchers have demonstrated the use of satellite remote sensing derived information in conjunction with canal flow data for the evaluation of irrigation command [87–89]. A considerable amount of work has been undertaken in the past 30 years to develop a framework for irrigation performance assessment, related to equity, adequacy, reliability, productivity and sustainability [5,20,30–32,90–93]. A list of irrigation performance indicators that can be quantified by use of remote sensing has been proposed by [20].

In our study on the Doukkala irrigation scheme, irrigation performance was assessed on the basis of adequacy of irrigation water allocations by CGR. For a real assessment of irrigation performance, the precipitation should be taken into account by using Irrigation Water Requirement (IWR), i.e., CWR—Precipitation. Subsequently the value of performance indicator 2 (IP2) is determined for all CGRs [5]:

$$IP{2}^k=\frac{{\sum }_{x,y}E{T}_c^k(x,y)-P^k}{V^k}$$

(19)

where: ET_c is the total crop evapotranspiration (m³); k is the reference unit (in our case is the CGR), P^k is the total precipitation in the reference unit k (m³), and V^k is the volume received at reference unit k (m³).

In this study, we calculated directly ET_c (x, y) using the k_c-NDVI and the analytical methods, which can then be integrated over the area of interest and compared with irrigation volumes to determine IP2 with Equation (19). In Equation (19), ET_c (x, y) in (m³) is obtained by multiplying ETc in (mm) by the area of the pixel (x, y).

5. Results and Discussions

5.1. Retrieval of Crop Bio-Physical Variables

The bio-physical variables required for ET_c estimation (the surface albedo (r), Leaf area index (LAI) and crop height (h_c)) are derived by applying the equations described in Section 4.3, to the time series constructed with the images listed in Table 2.

5.1.1. The Surface Albedo r

The surface albedo r depends on the sun elevation and zenith and azimuth view angles. The effect of sun elevation on the surface reflectance has been quantified by [77]. Since we used multi-sensor satellite data to construct albedo time series, the differences in terms of viewing angle and sun elevation were taken into account in the atmospheric correction and by considering the surface observed as lambertian. However, the differences in the spectral (number and width of bands) and spatial properties of the sensors used affected the estimated r.

To assess the impact on our retrievals of surface albedo (r) due to the differences in spectral and spatial properties of the sensors, we compared (r) estimated by means of different sensors in the same date and area. This was possible within the Faregh district on 26 April 2013 by comparing HRVIR1 (SPOT4) and OLI (Landsat8) r-estimates.

We chose 50 random pixels and then performed a linear correlation of OLI vs. HRVIR1 r-estimates (Figure 3). The RMSE = 0.0135 and the correlation coefficient R² = 0.768 indicate a good agreement of the estimated albedo with two sensors. The values of coefficient (0.77) and the offset (0.036) suggest that the residual error is not negligible. To evaluate whether this might be due to collocation errors of our data points, we repeated the comparison using larger samples.

We selected 10 independent and heterogeneous samples (20 pixels × 20 lines) of (r) estimated with each sensor (HRVIR1 and OLI). The analysis of the two populations of samples is presented in Figure 4. OLI data gave higher mean values of albedo than HRVIR1, with a lower spatial variability. The average of the 10 values of albedo estimated with HRVIR1 is 0.12 against 0.13 estimated with OLI: this slight difference is mainly due to the contribution of the blue band (not sampled by HRVIR1) to the OLI albedo.

We have evaluated the significance of the differences in the mean values of (r) estimated with the two sensors by calculating the mean and standard deviation value of (r) over all samples and applied a t-test (Student test) for the case of two variables with different variances. Taking into account the sample size and the standard deviations, the t-test (t) confirmed that the difference between these two estimates of albedo is not significant with 0.18 bilateral where t is the value of the t statistic, and p is the threshold of significance). The t-test gave very similar results for all samples except sample 5. The difference between the mean OLI-albedo and the mean HRVIR1-albedo was 0.01: this gives a difference in ET_c with 0.01 mm/day in winter and 0.05 mm/day in summer, which are both negligible.

Generally, in the entire study area, the surface albedo is low in winter, with a high frequency of values arranged between 0.1 and 0.15, and larger in summer with a high frequency of values between 0.14 and 0.23 (Figure 5).

5.1.2. The Leaf Area Index (LAI)

For the estimation of LAI using Equation (16), we have to determine three parameters α*, the coefficient C (Equation (15)) and WDVI_∞. In our application, we applied α* = 0.37 estimated by [78].

The C values were estimated by fitting a soil line to the scatter plot (Figure 6) of red versus NIR reflectance for all images of each study area (see Table 1). To determine the soil line, at first the dataset of each time series (study areas) was divided into multiple 0.002 intervals of red reflectance. Thus, within each interval the minimum value of infrared was selected [94]. To determine the soil line slope, a linear regression model was applied to the resulting r_ir vs. r_ir (minimum) subset taking into account only r_ir values less than 0.4 (bare soil). The slope of the soil line gives the value of C (Equation (15)) for each area: C = 1.20 (Sidi Bennour), C = 1.02 (Zemamra) and C = 1.25 (Faregh).

The values of WDVI_∞ were calculated from the WDVI time series for each study area. In each image, the mean WDVI (WDVIm) and the standard deviation (σ) were calculated. The WDVI_∞ of each image was estimated as WDVIm + 3σ to filter out outliers. Finally, the WDVI_∞ for each zone was calculated as the mean value (over all images) of WDVIm + 3σ. These values were equal to 0.46, 0.4 and 0.51 in Sidi Bennour, Zemamra and Faregh, respectively.

The time series of the maps of the variables, i.e., LAI in Figure 7, clearly show the pattern of land cover and its temporal evolution. In the Sidi Bennour district, LAI values were small and quite variable in space on 15 December (mean = 1.78; σ = 1.97) when some crops are at the beginning of their development stage (sugar beet and alfalfa) and others are going to be sown (wheat) while sugar beet, wheat and alfalfa reach the maximum vegetative development in February (mean = 2.98; σ =2.11). In July (mean = 1.28; σ =1.36), LAI values are lower because of the smaller cultivated area.

LAI maps present sparse outliers that are generally isolated except for some plots where the high LAI values are further extended in space (see arrow in Figure 7c). Outliers in LAI maps are mostly due to saturation effect [95], i.e., the received radiance at the satellite exceeds the maximum value that can be measured by the sensor. This occurs in general for the NIR band. In the Faregh area, saturated pixels (filtered out in the generation of the SPOT 4 data product) were surrounded by pixels with very high, although not saturated, NIR reflectance (Blooming effect [96,97]). This case gives very high WDVI values and therefore very high LAI values (Figure 8).

5.1.3. Soil Moisture and Radiation Control: LAI vs. Albedo

The evolution over time of LAI and albedo is correlated (Figure 9). During spring and summer, the soil surface becomes drier and crops reach maturation and are harvested, leading to an increase in albedo and a decrease in LAI. However, we notice that both albedo and LAI increase in autumn and winter because of initial wetting by precipitation and crop development at the beginning of the growing season.

5.1.4. The Crop Height h_c

To estimate the aerodynamic properties of a vegetation canopy, we need the zero plane displacement height (d), the roughness length for momentum, z_0m and heat, z_oh, transport. When the surface is uniformly covered by vegetation, these properties are simply related to the crop height (h_c):

We analyzed the spatial variability of the estimated h_c in winter (December 2012) and summer (July 2013) (Figure 10). Our own field observations (see Section 3.3) and local knowledge (ORMVAD personal communication) indicate that dominant crops in our study area (sugar beet, wheat, maize and fourage) have a maximum h_c of approximately 1.2 m for wheat in winter, 2 m for maize in summer and 0.6 m for sugar beet. For the perennial crops, such as trees, we expect some high values in a few plots (max h_c = 3 m).

In winter, we notice a significant dominance of h_c values between 0.1 and 0.5 m, with a lower frequency of values ranging between 2 and 3 m. In summer, a high frequency of small values of h_c was noted since winter crops have just been harvested and summer crops are at the beginning of the development stage (maize, fruits and vegetables). The spatial variability of h_c is larger in winter than in summer.

The high values of h_c in winter could be due to an over estimation depending on the relationship used to derive Z_0m (see Equation (18)). When crops are completing the vegetative development, their chlorophyll content starts to increase even though their h_c remains constant. During this phase, the high values of NDVI give high values of h_c estimated by Equation (18). The effect of over-estimated h_c values on ET_c estimation was evaluated by a sensitivity analysis described in the following Section 5.2.

5.2. Sensitivity Analysis of ET_c to Bio-Physical Variables

The ET_c calculated with the analytical method (Equation (8)) depends explicitly on r, LAI and h_c, while the k_c-NDVI does not depend directly on any of these variables, since ET₀ is calculated using constant values, i.e., r = 0.23, LAI = 2.88 and h_c = 0.12. The overall dependence of ET_c on crop conditions is accounted for by the value of k_c and its evolution over time. The spatial variability of ET_c calculated with the analytical method can, therefore, be different from ET_c calculated with the K_c-NDVI method, because of this different sensitivity to land surface conditions.

5.2.1. ET_c (Analytical) versus (r-LAI)

Figure 11 shows ET_c vs. r for different LAI values. Here, ET_c was calculated using meteorological data observed on 10 December 2012 and 15 July 2013, i.e., for low and high values of R_n. The value of h_c was set to 0.40 m. These variables affect directly the values of ET_c calculated by means of the analytical approach (see the points in blue color). The relationship between ET_c and r can be very well approximated by a linear function, i.e., ET_c (analytical approach) decreases with increasing value r, and increases with increasing value of LAI. The sensitivity of ET_c to LAI is higher than to r. The impact of r on ET_c (analytical approach) is slightly more pronounced in summer with higher solar irradiance.

To assess the sensitivity of ET_c calculated from K_c-NDVI method to LAI, we have derived NDVI using the following formula [98]:

$$\text{NDVI}=0.0653\hspace{0.17em}\text{ln\hspace{0.17em}}(\text{LAI})+0.5872$$

(20)

We notice that ET_c (K_c-NDVI) increases slightly with increasing value of LAI in both winter and summer, showing a small deviation and sensitivity to this geometric and structural variable. ET_c (K_c-NDVI) does not depend on r. The difference between ET_c estimated with the two methods, however, becomes greater with increasing r, and even more pronounced in summer with increasing solar radiation. In this context, we analyzed the temporal and spatial variability of ET_c estimated with the two methods at different scales in the following section.

5.2.2. ET_c (Analytical) versus (r-h_c)

We have assessed the sensitivity of ET_c to h_c and r under winter conditions (10 December 2012) and summer conditions (15 July 2013). In both cases, the value of LAI was set to 2, the surface albedo was varying from 0.05–0.25 and h_c from 0.05–3 (Figure 12).

In general, we notice that h_c is not a critical variable for the estimation of ET_c under both winter and summer conditions. As illustrated in Figure 12a, ET_c hardly depends on h_c, and increases very little with increasing h_c. The increase is higher in summer when high values of the vapor pressure deficit (D_e) occur (Figure 12b). According to Equation (11), the direct effect of h_c is negligible because it appears in both the numerator and denominator of the argument of the logarithm.

The variability of ET_c as function of albedo r is lower. Assuming a value of h_c equal to 0.4 m, this assumption determines an error of ET_c not larger than 0.2 mm during winter, and not more than 1 mm under summer conditions. The increase of ET_c that corresponds to the decrease in r is slightly more pronounced in summer conditions than in winter because of the higher solar irradiance in summer.

5.3. Estimation of Crop Evapotranspiration ET_c: kc-NDVI vs. Analytical Approach

5.3.1. Temporal Variability

Figure 13 shows the temporal evolution of daily ET_c estimated by the k_c-NDVI method and the analytical approach in the Sidi Bennour, Zemamra and Faregh irrigation districts during the growing season of 2012/2013. ET_c estimated by the two methods has rather similar evolution while the difference increases significantly from winter to summer (Figure 13), thereby supporting the hypothesis raised previously in Section 5.2, about the combined effect of LAI and r in ET_c values in summer. The values of ET_c (analytical approach) are slightly higher than ET_c values calculated by K_c-NDVI method. The radiative term of Equation (9) increases with decreasing albedo. The albedo used in the K_c-NDVI method is constant (r = 0.23) and higher than the mean value of albedo observed in our study, which explains the higher ET_c values obtained with the analytical method. In another study [11], it has been shown that the k_c-NDVI approach without a local calibration produces an average over-estimation of ET_c of 17% in the case of corn and 19% for alfalfa.

In the Faregh district, both methods show a high value of ET_c at the beginning of the growing season. This is mainly due to the values of solar radiation Rs, relative humidity RH, air temperature Ta and especially the vapor pressure deficit D_e observed on the day of acquisition of the satellite data.

In general, the temporal evolution of ET_c reflects the dominant crop development in the study area. In Sidi Bennour and Zemamra, the dip in June represents the transition between winter crops (wheat and sugar beet) and summer crops (Maize).

5.3.2. Spatial Variability

Since the mean value of ET_c estimated by the two methods is similar, we evaluated the spatial distribution of ET_c at different spatial scales by calculating statistics (

Table 5. Statistical variability of the daily ET_c values at different spatial scales.

		Kc-NDVI Approach		Analytical Approach
		Mean	Standard Deviation	Mean	Standard Deviation
20 × 20	December	1.22	0.25	1.47	0.55
	April	4.06	0.74	4.08	0.90
	June	4.22	0.76	4.24	0.92
200 × 200	December	1.26	0.36	1.55	0.74
	April	4.58	0.72	4.78	1.08
	June	4.43	1.08	4.61	1.31

) of ET_c with samples of 20 pixels × 20 pixels with respect to 200 pixels × 200 pixels.

The mean values of ET_c are quite similar for the two samples while the standard deviation (i.e., the spatial variability) of ET_c estimated by means of the analytical approach is significantly larger than ET_c estimated with the k_c-NDVI method. This applies to all months and depends on the combined effect of surface albedo, LAI and crop height in the analytical method.

The spatial variability of ET_c must be taken into account when evaluating remote sensing estimates against point observations (see next section).

5.3.3. Validation

The primary variable of interest to compute the performance indicator (Equation (19)) is the maximum evapotranspiration ET_c. Accordingly, we compared our satellite-based estimates of ET_c, by the k_c-NDVI and analytical methods, with values calculated by means of the dual crop coefficient approach (k_cb), using our ground observation of f_c and h_c. We used the k_cb method as reference, because it is the most accurate for partial canopies and it takes explicitly into account vegetation fractional cover and crop height. It should be noted that the three methods are completely independent except for the use of the solar irradiance and the vapor pressure deficit.

We noted in some cases anomalous values in the satellite-based estimates of ET_c. We identified outliers in two ways: by filtering out estimates deviating more than 2.5σ and more than 2σ from the mean value of the difference. The latter filter leaves out about 18% of data points for either method, while the former about 13%.

The K_c-NDVI method gave a better agreement with the reference ground based ET_c with RMSE=0.86 mm/d and RMSE=0.79 mm/d, when applying the 2.5σ–2σ filter (Figure 14(Left)). Contrariwise, the analytical method gave a RMSE = 0.99 mm/d and RMSE = 0.89 mm/d when applying the same filters (2.5σ and 2σ) (Figure 14(Right)).

The analytical approach gave slightly higher RMSE than the K_c-NDVI method, namely it was 13% higher when filtering values greater than 2.5σ and 11% when applying the 2σ threshold to identify outliers. In a previous statistical analysis, the two methods provided rather similar mean values both at the regional scale (Figure 13) and for the different sample sizes (20 × 20) and (200 × 200) (Table 5). However, since in the study area the typical plot size is 1 ha, to perform the validation of satellites based estimates using ground observations, we used 3 pixel × 3 pixel samples to extract the ET_c values from the maps obtained with Landsat 8 data. We carried out a statistical analysis of ET_c estimated with the K_c-NDVI method and the analytical approach at the plot scale using this sampling scheme. The ΔET_c samples were analyzed separately for each acquisition date (Figure 15). The standard deviation of ΔET_c gives an indication of how likely the agreement is between the ground measurements and satellites estimates. The time series of σETc was evaluated for each sample (from A–M) (Figure 16).

As expected, we noticed that ΔET_c is significantly larger at the plot scale (0.28 mm·d⁻¹ in December and 1.16 mm·d⁻¹ in June) than at the regional scale (0.25 mm in December and 0.02 mm in July (Table 5)), especially under summer conditions with high values of deficit vapor pressure (D_e) and solar radiation (R_s). The spatial variability of ET_c within (1 ha) was very large especially for the samples B, C, F, H, K, L and M. At such locations, a large difference should be expected between in-situ and remote sensing observations. The analytical approach captures the spatial variability of ET_c better than the Kc-NDVI method.

5.4. Irrigation Performance Indicator

5.4.1. CWR and IWR versus Water Allocation

We collected the monthly surface irrigation water volumes allocated in the Low Section of the Doukkala irrigation scheme, in each district and each CGR (Table 1) (ORMVAD-personal communication). Monthly values of the performance indicator IP2 were calculated for each CGR and for each district, since data on water allocation were provided at this temporal and spatial aggregation. To obtain the monthly IP2 values, the remote sensing pixel-wise estimates of K_c available on specific days have been interpolated linearly to obtain monthly ET_c values by multiplying the interpolated daily K_c value by the daily ET₀.

Less irrigation water is allocated in some cases when precipitations are significant. Irrigation Water Requirement (IWR) has been obtained by subtracting precipitation from the mean value of CWR (mm/month) estimated using the analytical approach. The comparison between CWR, IWR and water allocation (mm/month) for different districts and CGR, irrigated with different irrigation systems, is shown in Figure 17.

CWRs were larger than water allocation for both the entire irrigation scheme and the irrigation units (CGR). The irrigation water deficit was low at the beginning of the growing season (December–February), and larger at the end of the season (June and July). Except for Faregh district, the mismatch between CWR and allocations is significant for all months. Water allocation was roughly constant throughout the year, irrespective of the increasing water requirements during summer (for Maize). In the district of Sidi Bennour and Zemamra, it should be taken into account that precipitation (IWR) slightly reduces the mismatch between water allocation and demand, especially in winter when in some cases, the irrigation water deficit is converted into excess, i.e., February for a number of CGR.

Unfortunately, the gaps in satellite data time series (January, March, May) in Sidi Bennour and Zemamra did not allow capturing fully the temporal evolution of water allocation and requirements. The total values of pixel-wise interpolated ET_c in January, February and March were respectively 194 mm and 204 mm against 80.5 mm and 117.19 mm of irrigation water depth for Sidi Bennour and Zemamra district, and a total of 130.5 mm of rainfall for each district. This implies that IWR = 64 mm (Sidi Bennour) and IWR = 74 mm (Zemamra), i.e., water allocation was adequate. We concluded that from January to March, which represents the critical stage and the months of maximum development of the dominant crops i.e., sugar beet and wheat, the contribution of precipitation to meet CWR was significant.

In the Faregh district, the precipitation does not reduce the mismatch between requirement and allocations in all cases, at the district and CGR scale.

The ratio between IWR and water allocation (IP2) allows us to assess and understand the adequacy of water allocation in the entire area, in different locations of the primary canal of irrigation.

5.4.2. Spatial Distribution of IP2

The ratio between IWR volumes and surface irrigation water volumes gives the value of IP2. We have calculated IP2 for each CGR and district. An example of the spatial distribution of IP2 in February 2013, for the Sidi Bennour and Zemamra districts and the individual CGRs is shown in Figure 18.

We noticed that the irrigation performance is not uniform over the whole study area. In some cases, the IP2 was lower than 1, which means that water allocation exceeded irrigation water requirements, i.e., CGR 338 in Sidi Bennour and CGR 322, 325 in Zemamra. To some extent, this excess is necessary to compensate for water losses and it remains to be evaluated whether a fraction of it could be used for supplemental amount of irrigation water volumes in other CGR where crops suffer from a significant water deficit, such as CGR 330, 333, 336, 337, 320, and 321.

The irrigation performance indicator IP2 was used to assess the spatial pattern of adequacy between water consumption and allocation in the head–tail reaches of the primary irrigation canal (Figure 19). The water allocations show for some dates that no irrigation was provided to farmers in a given CGR. Accordingly, we have evaluated the IP2 head–tail end pattern only for February, April and June of the 2013 growing season.

Large differences were found in irrigation performance between the head, the middle and the end of the system.

In February 2013, the mean IP2 was approximately 1.06 in the beginning and the middle of the primary irrigation system which means a properly performing irrigation, but was around 1.33 at the end of the system. Although water is reported to be sufficiently available during the main season, it can be concluded that there are significant differences in adequacy towards the tail-end of the system (CGR 322, 321, 320 and 325).

In April and July 2013, the adequacy was lower at the beginning and the end of the system and better in the middle of the system, i.e., no clear head–tail end pattern in irrigation performance was observed.

These two results can be explained by the continuous control of water flows in the irrigation system, which apparently was less effective in offsetting head–tail end patterns in February 2013. On the other hand, it should be noted that the irrigation performance was lower in April 2013, when water allocations were about half the water requirement.

6. Application in Irrigation Water Management

In this study, we have demonstrated the potential of using satellite remote sensing as a practical tool for CWR estimation for improved understanding of water use in major irrigation schemes such as the Doukkala. Repetitive multispectral and high resolution imaging of this agricultural area was used to provide a precise and quantitative evaluation of the crop water needs during different irrigation periods during the growing season of 2012/2013.

In practice, the information provided by remote sensing could be used for irrigation water management in two ways: pixel-wise CWR data or aggregated CWR data by CGR or district.

The pixel-wise CWR data provide a reference for better precision in quasi-real time scheduling of irrigation water depth. The primary users of this information are farmers and the operators of the tertiary canals. As shown in Figure 14b, the difference between reference ground based and estimated CWR (using the analytical approach) was RMSE = 0.86 when applying the 2.5σ filter. In this case, the pixel-wise CWR data will present useful information for precise irrigation scheduling, when the spatial variability of CWR in the plot scale is higher than RMSE.

The CWR data aggregated by CGR and district provide a reference to adjust water allocation. The primary user of this information is the water management body, in our case ORMVAD, at the different management levels involved in planning and operation of water distribution. In general, it is necessary to take into account the difference between CWR, irrigation water requirement (IWR), and net irrigation water requirement (NIWR) in order to determine water allocation.

The assessment of the irrigation performance can only be done by simultaneously assessing the CWR, IWR, and more precisely using NIWR, which is the quantity of water necessary for crop growth, taking into account the rainfall. Information on irrigation efficiency is necessary to be able to estimate IWR given NIWR. The water balance in the soil–plant–atmosphere continuum can be described by models such as the Soil Water Atmosphere Plant (SWAP) model to estimate the NIWR by parameterizing root water uptake as a function of soil pressure head and soil water deficit. We have estimated NIWR by adding the soil water deficit on all dates of irrigations, where the latter are determined by maintaining crop transpiration at the potential rate. An estimation of the monthly NIWR (mm) for the dominant crops in the study area (wheat, sugar beet and alfalfa) using the SWAP model [99,100] for the growing season of 2000/2001 gave a mean value of 54.7 mm·m⁻¹ [101].

As illustrated in Figure 17, the CWR is significantly larger than water allocation for the entire study area with 20–30 mm·m⁻¹ of mismatch in winter for both the Sidi Bennour and the Zemamra districts. In summer, the CWR becomes much larger than water supply by 90–145 mm·m⁻¹ in June and July, respectively. In the Faregh district, CWR is much higher than irrigation water depth, i.e., around 50 mm·m⁻¹ in winter and 70 mm·m⁻¹ in summer.

This mismatch between requirement and allocation is improved for the entire study area when taking into consideration rainfall by means of IWR. For example, in February 2013, water allocation was almost equal to IWR in Sidi Bennour, while being just 10 mm lower than IWR in Zemamra. In the same month, water allocation exceeded IWR by 38 mm in CGR 325 and by 14 mm in CGR 335. In this case, water allocations were adequate.

In general, under summer conditions with an absence of rainfall, the mismatch between requirement and allocation remains high. However, the NIWR (54.7 mm·m⁻¹) is rather close to water allocation in winter and adequacy is reasonable in summer with the mismatch decreasing from 90–145 mm·m⁻¹ to 30–50 mm·m⁻¹. We cannot conclude, however, that water allocations can meet NIWR since we should take into account conveyance and operational irrigation water losses from the secondary canal to the plot. Bos, M.G. et al (1974) [102] evaluated over 250 irrigation schemes worldwide and estimated irrigation water losses at 50%. Taking into account water losses, net (on farm) water allocation would still be lower than NIWR in winter and much lower in summer.

Spatially speaking, and as shown in Figure 18, the adequacy of water allocation could be improved by reducing the water excess in some CGR and by using it in others where a deficit has been assessed. Likewise, the temporal distribution of water allocation could be improved by reducing water allocation at the beginning of the growing season and increasing it in summer.

7. Conclusions

The study confirmed that crop water requirement (CWR) can be estimated with satisfactory accuracy using a generic algorithm, which does not require prior classification of crops. The appraisal of irrigation performance in terms of adequacy between requirements and water allocation at both the district and CGR (Centers of Irrigation Management) level documented a significant mismatch of requirements and allocations. Taking rainfall into account, the difference between requirements and water supply becomes acceptable in winter, but the irrigation water deficit increases in summer (90–145 mm·m⁻¹). The mismatch in both winter and summer becomes even lower when the net irrigation water requirement (NIWR) is taken as a reference.

This was achieved by constructing a time series of multi-spectral satellite image data with different spatial, temporal and spectral resolutions (SPOT4 HRVIR1, Landsat8 (OLI) and RapidEye (REIS)), and implementing (semi-) empirical algorithms to assess phenology (K_c-NDVI method) and retrieve canopy biophysical parameters (analytical approach), such as surface albedo (r), leaf area index (LAI), and crop height (h_c). The spatial distribution of K_c, r, LAI and h_c was used in conjunction with ground-based meteorological data for mapping maximum crop evapo-transpiration (ET_c). These methods are fast, robust and easily applicable to large data sets and thus suitable for operational purposes.

Calera, A.B. et al (2005) [49] recommended the use of K_c-NDVI method to estimate CWR, and mentioned that the methods based on the retrieval of canopy biophysical variables are very complicated and give similar final results to the K_c-NDVI method. In our case, we have assessed the difference between the two methods, and we have concluded that they give rather similar mean values of ET_c but the analytical approach captures much larger spatial variability which is very useful for precision irrigation scheduling using pixel-wise CWR data.

In general, both spatially and temporally, the adequacy of water allocation to requirements could be improved by judicious management of irrigation, i.e., by reducing the water excess in some CGR/date and using it in others in deficit.