Optimizing the Sowing Date to Improve Water Management and Wheat Yield in a Large Irrigation Scheme, through a Remote Sensing and an Evolution Strategy-Based Approach

: This study aims to investigate the effects of an optimized sowing calendar for wheat over a surface irrigation scheme in the semi-arid region of Haouz (Morocco) on irrigation water requirements, crop growth and development and on yield. For that, a scenario-based simulation approach based on the covariance matrix adaptation–evolution strategy (CMA-ES) was proposed to optimize both the spatiotemporal distribution of sowing dates and the irrigation schedules, and then evaluate wheat crop using the 2011–2012 growing season dataset. Six sowing scenarios were simulated and compared to identify the most optimal spatiotemporal sowing calendar. The obtained results showed that with reference to the existing sowing patterns, early sowing of wheat leads to higher yields compared to late sowing (from 7.40 to 5.32 t/ha). Compared with actual conditions in the study area, the spatial heterogeneity is highly reduced, which increased equity between farmers. The results also showed that the proportion of plots irrigated in time can be increased (from 40% to 82%) compared to both the actual irrigation schedules and to previous results of irrigation optimization, which did not take into consideration sowing dates optimization. Furthermore, considerable reduction of more than 40% of applied irrigation water can be achieved by optimizing sowing dates. Thus, the proposed approach in this study is relevant for irrigation managers and farmers since it provides an insight on the consequences of their agricultural practices regarding the wheat sowing calendar and irrigation scheduling and can be implemented to recommend the best practices to adopt.


Background and Context
In Morocco, cereals are the main cultivated crop, both in terms of utilized agricultural area (UAA), and in terms of value in agricultural production, domestic food security and livestock survival [1]. This agricultural crop covers approximately 5.3 million ha and accounts for 65% of the total cropped area [1]. The production of wheat (durum wheat and common wheat) practiced by almost all Moroccan farmers concerns 3 million ha [1]. However, production remains highly variable from year to year depending on climatic variability, mainly annual precipitation levels, and does not follow national demand, which increases dependence on imports particularly during dry years [2]. Recent studies expected that global warming will threaten most of the world's wheat production and may cause major droughts in 60% of wheat-growing areas [3]. For Morocco, drought is increasing in frequency and intensity [4]. For instance, the observed annual average temperatures increased every decade since 1970 and the Intergovernmental Panel on Climate Change (IPCC) projected a continued warming trend with an increase of about 1 to 1.5 degrees Celsius by 2050 [5]. Additionally, a recent World Bank study on climate change in Morocco shows that precipitation will decrease by 10 to 35 percent [2].
Several studies showed that high temperature and water deficit during the key reproductive phases have significant effects on the length of the wheat development periods [6] and decrease of wheat yield and grain quality [7][8][9][10][11][12]. Thus, to mitigate the effect of climate change on wheat productivity, conventional agricultural practices such as sowing date must be adjusted to increase yield [13][14][15][16]. Sowing time impacts the duration and timing of developmental phases by avoiding climatic constraints in wheat sensitive stages [15,[17][18][19]. Different approaches have been developed to estimate crop sowing dates by defining a set of rules based on crop and climate characteristics [20,21]. In other common approaches, crop models were used to identify the best agricultural practices in a specific climate [22,23], or to study the effect of shifting the sowing date on crop yield [24][25][26][27][28] or to assess crops performance depending on climatic conditions [25,29,30]. Other studies have investigated the effect of spatiotemporal variability in sowing dates, at the regional level, on yield [31][32][33]. But among these studies, few focused on the optimization of the sowing window in arid and semi-arid regions such as Morocco under present climatic conditions [18].
Conversely, water scarcity is one of the key factors limiting agricultural development and thus food security, especially in arid and semi-arid regions where irrigation consumes more than 85% of the available water resources [34]. This is the case of the Haouz plain (central of Morocco), classified as a semi-arid region [35], where water availability is a major limitation for cereal production [36].
In this study, we focused on surface irrigation systems, which occupy about 80% of irrigated areas worldwide [37] and about 76% in Morocco [1], particularly the large-scale irrigation schemes under rotational water supply. It has been demonstrated that existing rules for water distribution among farmers, often based on a fixed-area proportionate depth of water applied in every irrigation round irrespective of the crops and their growth stages, induce an inadequate water supply and have negative effects on wheat productivity [38,39]. The irrigation scheduling approach adopted by the managers results in either an excess of water or a lack of water leading to a more complex management: (1) at the beginning of the agricultural season, the sowing dates depend on the first rainfall and the number of irrigation rounds and water amounts to be applied, depending mainly on the available water resources (the water stored in the dams), (2) for each irrigation round, the irrigation distribution is mainly based on hydraulic and technical constraints (location of the fields, canals discharge, staff operating the water intakes, etc. . . . ) and the starting date is based on visual observations, and (3) water is applied based on a fixed-area proportionate water rate without considering the crop water demand, the crop type and the level of water stress of each crop, which is the most influential factor on crop growth and yield and the heterogeneity of farmers' behaviors and objectives.
Consequently, efficient and optimized spatiotemporal allocation of irrigation water in a surface irrigation system to improve crop water use efficiency and to ensure adequate productivity in irrigated agriculture during water deficit years is crucial for the sustainability of agricultural development in the arid and semi-arid regions.

State of the Art
Optimized allocation of limited available irrigation water resources can be achieved by adopting best agricultural practices related mainly to sowing dates and irrigation schedules, while taking into account both the heterogeneity of environmental and agronomic conditions (water stress levels, type of crop, soil characteristics, rainfall, etc.) and the technical and human constraints (behavior of the actors, the different interactions between them, their decision-making strategies, the constraints of the irrigation network, etc.). This is a complex problem that requires a large amount of precise information. However, with help of optimization approaches using mathematical calculation, objective function and existing constraints, it has been possible to handle the complexities of surface irrigation systems.
Many studies have been conducted to address different irrigation water management issues using initial optimization techniques, such as linear programming (LP), nonlinear programming (NLP), integer programming (IP), dynamic programming (DP), and quadratic programming (QP) [40][41][42][43][44]. To cover the drawbacks of these traditional methods, which are time-consuming, intelligent evolutionary algorithms, such as genetic algorithm (GA), evolutionary strategies (ES) and evolutionary programming (EP), has also been widely used in an array of agricultural problems, and their results shown to be more efficient and more promising for finding global optima [45]. For instance [46], comparing five evolutionary algorithms (real valued genetic algorithm, particle swarm optimization, differential evolution, and two evolution strategy-based algorithms) on the problem of irrigation optimization demonstrated that they are able to optimize irrigation schedules, achieving results that are extremely close to the theoretical optimum. The efficiency and the strength of genetic algorithm approach to provide optimal irrigation scheduling has been also demonstrated in several studies [47][48][49][50][51][52].
Among this evolutionary algorithms, the covariance matrix adaptation-evolution strategy (CMA-ES) algorithm is recognized as one of the most powerful methods in the black box continuous optimization framework [53], and it is particularly suited to problems containing many local minimums [54], as presented by the results of the competition at IEEE-CEC 2005 [55] and BBOB workshops at ACM-conferences GECCO 2010, 2011 and 2013 [53]. Comparative analysis with several algorithms has also shown that the CMA-ES algorithm is one of the most effective meta-heuristics for dealing with difficult optimization problems [53]. Unlike genetic algorithms based on biological analogies, CMA-ES uses a space sampling of the parameters by a population of points drawn according to a Gaussian law. The number of individuals depends on the dimensionality of the problem. The CMA-ES algorithm is generally used to search for the solution of complex problems in order to obtain a reasonable calculation time.
In a previous study [56], we have proposed an approach based on the CMA-ES evolutionary strategy algorithm for optimizing irrigation scheduling. The obtained results demonstrate that it can be considered as an efficient tool for planning irrigation schedules by considering crops water needs and timing. This research was conducted in our study area and proposed an optimization during only one irrigation round in a surface irrigation network, without optimizing the spatiotemporal distribution of sowing dates and irrigation rounds along wheat season.
In this context, the present study aims to propose a decision support approach for rational and efficient management of irrigation in large surface irrigation systems, based on an optimization of the spatiotemporal distribution of sowing dates and irrigation, while evaluating the effects on the irrigation rounds scheduling, the consumed irrigation water and the wheat yield. This approach can contribute to the adaptation of agricultural practices to water scarcity conditions and meet food production needs in response to climate change.

Study Area
Our study area is located in the Haouz plain (Figure 1a) in the center of the Tensift basin (Central Morocco). It is a vast plain 6000 km 2 in area (including approximately 3100 km 2 of irrigated area), characterized by a semi-arid continental climate, with low and irregular rainfall of 250 mm/year on average, a high temperature in summer (38 • C on average, in July) and low temperature in winter (5 • C on average, in February) [57]. Among the irrigated sectors of the Haouz plain, the R3 district (Figure 1b) was chosen as the main area to validate our research results. It is located about 40 km east of Marrakech city and covers an area of 2800 ha composed of around 745 individual fields, used mainly for cereal crops (45% in 2011-2012 and the remainder surface was occupied by annual crops, trees and fallow). R3 is flat and considered by many previous studies as a suitable area for satellite applications (e.g., [34,[58][59][60]). At field level, water is applied by flooding method. A fixed water rate is received by each plot, at each irrigation round, without taking into account the type of crops, water requirement, and the cultivated area. Water is provided to the R3 irrigated area from the "Hassan 1st dam" located in the neighboring catchment of "Oum Er-rbia" via the "Rocade" canal, and through a classical gravity network made of lined canals (Figure 1b). The latter is composed of two primary canals: P1 on the right bank of the sector and P2 on the left bank, and a network of secondary and tertiary open canals convey water to the cropped fields. The entire area is divided into irrigation blocks composed of 6 plots of 4-5 ha each, generally supplied by a tertiary canal feeding 1 to 10 blocks.
Irrigation and land use in this study area is managed jointly by the Regional Office of Agricultural Development of Haouz (ORMVAH) and the water users associations (WUA). Depending on the water allocated to irrigation, managers and WUAs decide the number of irrigation rounds and water rates to apply at each one. Then, the irrigation schedule is implemented throughout the season in collaboration between the ORMVAH and the farmers (i.e., the starting date and duration of each irrigation round), which is usually based on simple observations of weather conditions and/or on a visual assessment of the plants and soil moisture status.
This approach generates a high spatiotemporal heterogeneity of the applied irrigation volumes at the plot scale, whose complexity is usually not considered in the irrigation scheduling process. Thus, an optimized irrigation management for an efficient and productive water use in such complex system should be adopted.

The Proposed Approach
To help farmers and irrigation scheme managers make good decisions about the winter wheat planting calendar and irrigation schedules that allow optimized water consumption and crop productivity, we propose an innovative approach (Figure 2), which takes place in three phases of the agricultural season:

1.
At the beginning of the agricultural season, the spatiotemporal distribution of the sowing dates over the irrigation scheme is optimized, considering the irrigation network constraints and the adopted sowing scenario. The number of plots to be sown in a same day is defined based mainly on the capacity of the irrigation network, i.e., the maximum volume of water that can be delivered to the plots irrigated by the same canal, according to its maximum discharge and location (tertiary, secondary or primary). Thus, the total area that will be sown in a given day is determined such that it can be completely irrigated in one single day during an irrigation round; 2.
During the agricultural season, before starting each irrigation round, we elaborate an optimal irrigation planning (irrigation date and irrigation water requirement for each plot) over the irrigation scheme according to the sowing calendar scenario and the crop water stress levels in each field. The starting date of each irrigation round is calculated based on the water stress coefficient (Ks), which is used to measure the soil water deficit levels for the crop, as defined in the FAO56 model [61], based on the Normalized Difference Vegetation Index (NDVI), which is a common and widely used remote sensing index to assess the green biomass and the state and health of crops by measuring the difference between near-infrared and red light. Thus, the Ks is equal to 1 at the day of sowing (because the farmers sow after a heavy rainfall event or a pre-irrigation at the beginning of the season). Then, Ks is calculated using the daily soil water balance at the root zone (see Section 2.3). An irrigation round is decided when the Ks reaches a value less than 0.7 at most for 5% of the total cultivated area. This threshold of Ks has been identified as the value above which the water stress does not affect significantly wheat yield [36,62,63]. Then, for each irrigation round, the spatiotemporal irrigation planning is optimized based on the Irrigation Priority Index (IPI) [56] (see Sections 2.4 and 2.6). When a plot is irrigated, the Ks coefficient becomes equal to 1 and the process continues to define the starting of the next irrigation round. Additionally, no more irrigation is decided for wheat crop at the end of the season when the canopy cover reaches the half of the maximum canopy cover during the senescence stage (see Section 2.2); 3.
At the end of the agricultural season, an evaluation of the irrigation water consumed and the final yield expected is conducted in order to decide on the most optimal spatiotemporal sowing date calendar to be adopted, i.e., the best recommendations concerning the cultivation schedule which will allow to increase incomes and reduce the use of water resources during the agricultural season.
The initial distribution of sowing dates is simulated randomly (taking into consideration the window of sowing scenario) to avoid local solutions and to ensure that the best results are independent from these initial values. We reiterate the algorithm 10 times. The solution corresponding to the higher value of the total yield and the lower value of the consumed water resources is considered as the optimal solution.
We applied this approach for wheat plots during 2011-2012 agricultural season over the left side of R3 irrigation scheme irrigated by the primary canal R3P2 (composed of one primary canal, two secondary canals and 30 tertiary canals). The simulations concerned 116 plots, which constitute 45% of the global area. Six scenarios of sowing calendar were considered using a gap of 15 days between scenarios on a sowing window from 1 November 2011 to 31 January 2012 (1st sowing scenario 1 November 2011; 2nd sowing scenario 16 November 2011; 3rd sowing scenario 1 December 2011; 4th sowing scenario 16 December 2011; 5th sowing scenario 1 January 2012; 6th sowing scenario 16 January 2012). These predefined windows are identified based on physically-based crop models and validated by farmers. In our study area, the cereals were sown between November and January, reaching their maximum development between end of March and mid-April, and they are harvested in June. For each scenario, a simulation of the irrigation schedules for the entire agricultural season and an evaluation of the global irrigation water consumed and the crop yield expected are carried out. This approach enables the analysis of a set of scenarios-consequences in order to define the best sowing calendar that allow for optimal irrigation water needs and wheat yields. The proposed approach allows to identify the consequences of each sowing scenario, according to different criteria (environmental, economic, social, etc.) on water resources, crop growth and final yield at plot level over the irrigation scheme.

Dataset 2.3.1. Meteorological Data
To implement and evaluate the proposed approach, we used the meteorological and irrigation data for estimating water requirements, calculating the Ks coefficient, and for estimating the NDVI profile for each spatiotemporal sowing. We also used the Landsat-TM5 image to calculate the observed NDVI profile of each plot. The required data was collected during the 2011-2012 agricultural season in the study area.
The reference evapotranspiration (ET 0 ) is an indicator of the evaporative demand of the atmosphere and it is calculated using the FAO Penman-Monteith equation [61]. The climatic variables (air temperature and humidity, solar radiation and wind speed) needed to estimate ET 0 and rainfall were measured by an automatic weather station that was installed within the R3 district. More description about this weather station as well as the different sensors used for measuring the different climatic parameters are previously described by [59,64]. Figure 3 shows the daily ET 0 , temperature and rainfall calculated during the 2011-2012 agricultural season. The evolution of ET 0 is characteristic of semiarid continental climate: it is moderate during the autumn (2 to 3 mm/day) and higher in summer (6 to 7 mm/day, in May-June). The annual precipitation is low (equal to 220.7 mm) and irregular and temperature is high in summer (33.5 • C, in June) and low in winter (4.33 • C, in February).

Irrigation Water Supply
The detailed data about effective irrigation schedules (Table 1) at the field scale was time-consuming, and it was surveyed locally in collaboration with the ORMVAH service. The simulations concerned 116 plots which constitute an important global area of about 1260 ha (45% R3 area). Five irrigation rounds for the wheat growing season in 2011-2012 were applied, which durations lasted between 12 and 22 days as presented in Table 1.

Wheat Fields Identification
Land cover mapping was established using a classification method based on NDVI time series obtained from satellite images and was validated by field surveys and data provided by the ORMVAH. In our case, among 18 Landsat images collected during the 2011-2012 wheat season, only 2 images were affected by clouds and were not used in this study. We used 16 Landsat-TM5 images well distributed throughout the agricultural season 2011-2012. Thus, the daily values of the NDVI are easily obtained by linear interpolation [65]. Radiometric calibration and atmospheric corrections were performed on these images based on the reflectance of invariant objects (bare soil, fallows and houses) to provide NDVI time series over the irrigation scheme according to the reported method by [66]. The NDVI profiles obtained for each pixel were used to provide land cover map of the area, identifying the wheat fields concerned by our study, using a decision-tree algorithm based on NDVI thresholds characterizing particular crops phenological aspects [65].

Applied Sowing Dates Identification
To estimate the applied sowing date of each field at the agricultural season 2011-2012, we used the NDVI time series calculated from the collected Landsat-TM5 images. The sowing date is estimated from the emergence date, which is identified when the canopy cover (CC) reaches a value of 0.1 [61]. The values of the canopy cover are estimated from the observed NDVI at the plot level (see Section 2.4). The sowing date is then obtained by subtracting 13 days starting from the emergence date for the case of early sowing and 20 days for the case of late sowing. These values (13 and 20 days) were estimated based on observations monitored in the study area [62,64]. The date of December 20 is commonly considered in our region to distinguish between early and late sowing in accordance with the impact of agro-environmental conditions on wheat growth and development and the farmer's agricultural practices, particularly the wheat cultivars used [57,62,67].

Irrigation Network Identification
The irrigation network serving the farmer's fields with water was identified and mapped in a geographic information system (GIS) composed mainly of: • Canals: line type, characterized by an id number, maximum discharge, the flow direction, the geographic coordinates, the category (primary, secondary or tertiary) and plots that are fed; • Water intakes: point type, representing the devices that supply a plot or a set of plots. Each intake is characterized by an id number, the geographical coordinates, the maximum discharge, the number of plots it supply, the opening and closing times and the id number of the canal to which it is linked; • Plots: polygon type, characterized by an id number, the farmer's name, the type of crop, the id number of the canal and the water intake, the date of sowing, the overall area, the irrigated area and the effective irrigation water applied in 2011-2012 season (in terms of time and quantity).

NDVI Profiles Simulation
For the simulation of the NDVI profiles for each irrigated plot according to the adopted sowing date, we first applied the equations defined by the FAO crop growth simulation model "AquaCrop" [68,69] to generate the green canopy cover profile along the wheat season.
The decline in green crop canopy is described by: where CC is the canopy cover at time t (%), CC 0 is the initial canopy cover at t = 0 (%), CC x is the maximum canopy cover at the start of senescence (t = 0), CGC is the canopy growth coefficient (% or cumulative growing degree days), CDC is the canopy decline coefficient (% or cumulative growing degree days), and t is time (days or cumulative growing degree days).
The canopy cover profile is then described by the four parameters: CC 0 , CGC, CC x and CDC. CC 0 , CGC and CC x determine the time required to reach maximum canopy cover, whereas the canopy decline coefficient CDC determines the rate of the green canopy decline in the late season. For our proposed approach, we used the calibrated values of CC 0 , CGC, CC x and CDC obtained by [36] for wheat crop grown in the study area (R3) ( Table 2). Once the canopy cover profiles generated, the NDVI profiles were obtained using the relationship between NDVI and the canopy cover (CC) previously established by [70] for the wheat in R3 area (Equation (4)): NDVI = (CC + 1.18 × NDVI min )/1.18 (4) with NDVI min = 0.14 as the minimum value of NDVI, associated with bare soil and dense vegetation in our study area [70].
In the proposed approach, the generated NDVI profiles are used in the optimization process, considering that the available water is the only factor limiting the optimal crop growth.
The values of the canopy cover (CC) are also used to calculate the stopping irrigation window at each sowing scenario (see Section 3.2). We consider that the irrigation of a plot must be stopped once the CC reaches a value equal to the half of the maximum canopy cover at the start of senescence (i.e., CC x /2). In the proposed approach, the optimization process (see Section 2.7) has been defined such that the irrigation window never overlaps the stopping irrigation window.

Spatialized Estimates of Crop Water Needs and Water Stress Level
In order to predict the amount of irrigation water needed by each plot, we used the soil water balance approach based on the FAO-56 dual crop coefficient model, which is widely used and allows simple and robust modeling [39,61,70,71]. This approach quantifies the crop evapotranspiration under standard conditions (ET c ) as the product of the reference evapotranspiration (ET 0 ) and a crop coefficient that account separately for the crop transpiration K cb (basal crop coefficient) and for the soil evaporation K e (soil evaporation coefficient). Additionally, in order to account for the level of the crop water stress to decide when to trigger irrigation, we evaluated the water stress coefficient (K s ). Thus, the actual evapotranspiration (ET) is given by Equation (5): where K cb and K e are derived from NDVI observations using Equations (5) and (6) where CC is the canopy cover calculated with Equation (4). The water stress coefficient is given by Equation (8) [61]: where Dr is the root zone depletion in (mm) and RAW is the readily available soil water in the root zone (mm). When Dr ≤ RAW, Ks is equal to 1. TAW is the total available soil water in the root zone [mm], expressed by the difference between the water content at field capacity and wilting point (Equation (9)): where θ FC is the water content at field capacity in (m 3 m −3 ), θ WP is the water content at wilting point in [m 3 /m 3 ] and Zr the rooting depth in (m) which varies according to the plant growth stages. RAW is expressed by Equation (10): where ρ is the average fraction of TAW that a crop can extract from the root zone without suffering water stress (0-1). The daily soil water balance, expressed in terms of depletion at the end of the day by Equation (11): where D r,i−1 and D r,i are respectively the root zone depletion at the end of day i − 1 and i (mm), and Pi is the precipitation, RO i is the runoff, Ii is the net irrigation depth, CR i is the capillary rise from the groundwater table, ET c,i is the crop evapotranspiration, and DP i is the water loss out of the root zone by deep percolation on day i (mm).
2.6. Grain Yield Estimation 2.6.1. NDVI-Based Approach Ref. [72] is one of the first to examine the relationships between developed NDVI and the biophysical parameters of vegetation cover. Many subsequent studies have been based on the ability of vegetation indices to predict agricultural production and yield [73][74][75]. Among these indices, NDVI has been used successfully to predict wheat yields in Morocco both at regional and national levels [76]. Thus, in this study, we considered the work carried out by [74] who established an operational method for forecasting crop yields based on NDVI data. It consists in defining a linear relationship between wheat yield and the 10-day spatial accumulation of NDVI values during critical phases of crop development. This approach was tested and validated using data of real grain yield of winter wheat obtained on ten fields of the study area during the 2002-2003 agricultural season [62,77,78]. Thus, for the estimation of wheat yield, we used a relationship developed and validated by [78] which considers the climatic and phenological characteristics of winter wheat in the Haouz plain and showed that the last 10 days of March present the best correlation between yield and ten-day NDVI (Equation (12)). This is consistent with the phenological calendar of winter wheat in the Haouz region. During the second half of March, the wheat is fully developed, photo-synthetically active and the crop reaches the grain filling stage [67,77,79].
The AquaCrop is a water-driven model developed by the Land and Water Division of the Food and Agriculture Organization of the United Nations (FAO). In this model, transpiration (Tr) is calculated first based on a water balance model and then converted into biomass (B) using a conservative crop specific parameter called water productivity (WP) (Equation (13)). The biomass is then translated into grain yield (Y) using harvest index (HI) (Equation (14)) [69,80]. To run simulations, AquaCrop requires inputs of climate data, crop characteristics, soil characteristics and description of management practices [69]. A detailed description of the AquaCrop model is available in [69,80].
Ref. [36] calibrated and validated the AquaCrop model for irrigated wheat in our study area, on the basis of a comparison between measurements and the model estimations.
Ref. [15] used AquaCrop model to evaluate the impact of climate change on the wheat grain yield and water requirement in our study area. In this study, to simulate wheat grain yield, AquaCrop (v.6.1) was used, with calibrated parameters reported in [36].
where B is the above ground biomass (kg/m 2 ), WP is the normalized water productivity (kg·m −2 ) and Tr is the crop transpiration (mm).
where Y is the grain yield (kg/m 2 ), B is the total ground biomass (kg/m 2 ) and HI is the harvest index (%).

The Optimization Methodology
The optimization problem presented in this study is considered as a nesting of three sub-problems: sequential scheduling, linear optimization and traveling salesman type. The problem of scheduling displacements from one water intake to another to serve all the plots is equivalent to the travelling salesman problem, where we seek the best path to move from a starting point to a destination point. In our case, we have other problems to solve, such as the management of flow variations and the management of crop water requirements and water stress. This justifies that the studied scheduling problem is particularly difficult to formulate in a way to be able to solve it using an exact method or by an algorithm of the graph theory. Thus, the approach adopted in this work is based on a CMA-ES algorithm belonging to the family of meta-heuristic algorithms inspired by the theory of evolution.

Irrigation Priority Index (IPI)
The optimization of the irrigation schedules over the irrigation scheme in our approach is based on the Irrigation Priority Index (IPI) (Equation (15)), previously developed by [56], which characterizes the irrigation priority of a plot during each irrigation round. The IPI index is calculated for each plot i according to two main terms: a first term that measures the level of water stress of the field j relatively to all the plots concerned by the irrigation round, which allows to quantify the water needs just before the start of the irrigation round and a second term that measures the duration (in days) between the starting date of the irrigation round and the date the plot i receives water. The proposed expression for the IPI is given as follows: where K si is the water stress at the beginning of the irrigation round, and K smin and K smax are the spatial minimum and maximum of the K s map that represent the most and less stressed plots, respectively. δt i the time delay (in days) between the beginning time of the irrigation round and the irrigation time of plot i and T the duration of the irrigation round. Values of δt i range between 0 and T and are chosen as the decision variables of the optimization process. The proposed IPI equation limits its values between −1 and +1. An IPI = −1 corresponds to an extreme late irrigation (the most stressed plots irrigated the last day of the irrigation round); IPI = +1 corresponds to an extreme early irrigation (the less stressed plot irrigated in advance). Contrarily, an IPI = 0 corresponds to the plots irrigated at the appropriate time, which means a reasonable tradeoff between the water stress level and the irrigation time. Thus, values of IPI close to 0, on both sides, are considered as an indicator of optimal irrigation distribution during an irrigation round.

Sowing Date Optimization
The spatiotemporal sowing date optimization approach is defined in two steps: (1) the choice of a sowing scenario and (2) the development of a first optimal spatiotemporal sowing date distribution according to the chosen scenario. The optimization of this distribution is based on the constraints of the gravity irrigation network as well as on the available water resources at the beginning of the agricultural season. For the resolution of this optimization problem, we propose an approach based on the CMA-ES evolutionary strategy algorithm such as: (a) Inputs: available water resources, geospatial data and meteorological data; (b) Decision variables: the sowing date of each plot. We have 116 plots to be sown, which means 116 decision variables are considered; (c) Optimization constraints: we identified five constraints which we present by a decreasing level of priority.

1.
The capacity constraint: supplies must never exceed the total capacity of the canal. This constraint is expressed by Equations (16) and (17): where θ, Q and q are respectively the discharge of the primary, secondary and tertiary canals (in m 3 /s), i and j are respectively the tertiary and secondary canals' number, opened at the same time; 2.
The interval constraint: all the irrigation tasks must be scheduled within the specified irrigation round interval, expressed by Equations (18) and (19): where δt start is the starting date of the irrigation round; δt i the opening time scheduled for canal i; E the allowed duration of the irrigation round, and Di the irrigation duration scheduled for the canal I; 3.
The overlap constraint: all the practical actions must be applied dependably (not simultaneously), considering the geographical distance between water intakes and the irrigation time span required for each plot supplied by the same canal; 4.
The traveling time: the time required for the operator to travel from one canal gate to another to start/stop an irrigation which depends on the distance be-tween canals that can be significant. This time is calculated linearly based on the spatial distance between two intakes and assuming a moving speed of 30 km/h; 5.
The daily working time of operators: each scheduled task (start/stop an irrigation) must be scheduled within the specified working time that is between 8:00 h and 18:00 h; (d) The objective function: the first objective function aims to propose the best spatiotemporal sowing distribution which minimize the gravity irrigation network constraints by exploring the space of decision variables in an efficient manner. In case the constraints are not met, we adopt the penalty method by adding penalty function. Thus, the first objective function is defined by Equation (20): where r j is the weight associated to each constraint j, and P j the corresponding penalty expressed by Equation (21): where N is the total number of plots to be irrigated and ϑ j,k = 0 if the task is satisfying the constraints, ϑ j,k = 1 otherwise.

Irrigation Scheduling Optimization
As previously mentioned, the simulation of the optimal spatiotemporal distribution of irrigation dates at the start of each irrigation round is based on the water stress coefficient level (K s ), while considering the irrigation network constraints. For that, we consider that the best distribution is the one that minimizes the IPI index (in terms of absolute value). Thus, for the resolution of this second optimization problem, we applied the CMA-ES algorithm, such as: (a) Inputs: available water resources, geospatial data, meteorological data, soil proprieties, NDVI profiles, and K s ; (b) Decision variables: During each irrigation round, two types of decisions must be made: when to irrigate the plots and in what order (i.e., the sequence and scheduling of the irrigations). The answer to these two questions is integrated into a single decision variable that consists of defining the start date to irrigate each plot (the stopping date is calculated based on the water amount to be supplied to a plot given the calculated irrigation duration and the constant flow rate of the tertiary canal). The stating date to irrigate the plots can therefore be considered as the only decision variables in the optimization process. In this work, 116 plots are scheduled, which means that 116 decision variables are considered during the optimization process while considering the constraints related to the irrigation network; (c) Optimization constraints: the same as the first optimization (i.e., irrigation network constraints); (d) The objective function: the second objective function aims to propose the best spatiotemporal irrigation distribution, which minimize the IPI index and the gravity irrigation network constraints. Thus, the second objective function is defined by Equation (22): where k is the canal index and N is the total number of plots to be irrigated.
This optimization allows us to define the optimal opening dates of each canal. Then, to calculate the date to irrigate plots fed by a same canal having limited flow, we applied the shortest job first (SJF) non-preemptive scheduling algorithm [81], which is usually used for scheduling processes. In this algorithm, the process having the smallest value of K s coefficient is chosen for the next execution, which means that we first begin with the irrigation of the most stressed plots (depending on the canal flow rate), which imply to minimize the level of water stress within the block of plots irrigated by the same canal.
The irrigation duration of each plot is calculated according to Equation (23): where flow is the maximum canal discharge that should not be exceeded. It depends on the canal capacity, a design characteristic which is equal to 90 L, 60 L, or 30 L/s in our irrigation network. Thus:

Optimization of the Spatiotemporal Sowing Date Distribution
The first step of our proposed approach is the optimization of the spatiotemporal distribution of the sowing dates. Six sowing scenarios with different starting date were simulated. The results of this optimization for each sowing scenario are presented in Figure 4. This figure shows practically a similar optimal distribution of the sowing dates from one scenario to another attributed to the fact that this optimization is based on the same inputs data and the same irrigation network constraints, only the initial individuals corresponding to the decision variables (distribution of the sowing dates) in our CMA-ES optimization algorithm change. These latter does not affect the optimal sowing distribution from one scenario to another because we start with a random initial individuals and then during the first iterations of the optimization process, the exploration of the control parameters (decision variables) is weak ( Figure 5), which implies the adaptation of the covariance by the CMA-ES algorithm through enlarging the range of these parameters and then refining search of parameters to obtain a rapid reduction of the objective function. The convergence of this function for each sowing date scenario is presented in Figure 5.  We also notice in Figure 4, that for the different scenarios, a high number of plots are sown between the 5th and the 10th day (counting from the starting date of each sowing scenario). This is attributed to the fact that our optimization approach maximizes the number of plots sown at the beginning, in such way that we can irrigate the same number of plots at the same time when starting an irrigation round. Thus, the sowing date of the remaining plots (which are no longer numerous) is defined based on the irrigation duration of the first group of sown plots which will be irrigated by the same canal (to not exceed the maximum flow rate). This method was applied to lead to an irrigation schedule that follows the same sowing schedule in order to ensure a good balance between the water stress level and the irrigation date. An example of the sowing date distribution (first sowing scenario) is presented in Figure 6. It can be seen that a high number of plots (54%) are sown between the 5th and the 10th day.

Effect of Sowing Date Optimization on Wheat Growth and Irrigation Schedules
Once the distribution of sowing dates was optimized, the NDVI profiles were generated for each plot according to its sowing date using Equations (1)-(4). Figure 7 presents an example of the difference in NDVI profiles for a single plot (with area = 9.72 ha), by varying sowing date window from one scenario to another. This figure shows clearly that the duration of wheat season decreases from early to late sowing (from 208 to 137 days), which affects the crop development and growth [15,36]. The NDVI profiles of all plots in the study area were then used to calculate the water stress coefficient Ks, initially between the sowing date and the starting date of the first irrigation round, and between the different irrigation rounds, according to the approach presented in the Section 2.2. An example of KS maps at the start of the third irrigation round for the first sowing scenario (early sowing) and the sixth sowing scenario (late sowing) is displayed in Figure 8. The spatial heterogeneity of the Ks values between plots irrigated by the same canal and the effect of the spatiotemporal distribution of sowing dates on water stress are clearly observed. This heterogeneity is also illustrated in Figure 9, which presents the Ks values for each plot at the starting of the third irrigation round for each sowing scenario. We also notice in this figure that the optimal irrigation date follows the evolution of the Ks level, i.e., the most stressed plots are irrigated first and the least stressed are irrigated last.   The next step is the irrigation scheduling optimization process based on the IPI index derived from Ks stress coefficient maps elaborated at the beginning of each irrigation round. An irrigation round is decided when the K s reaches a value less than 0.7 at most for 5% of the total cultivated plots.
The results of IPI optimization at plot scale for each sowing scenario are presented in Figure 10. It can be seen that for almost all sowing scenarios, the majority of plots is irrigated at the appropriate time considering their water stress level, which illustrates the performance and the efficiency of the optimization approach. The optimized sowing distribution and irrigation schedules at the plot scale for each scenario leads to a high range of plots irrigated on a favorable time for most irrigation rounds (an average of 76.17% of plots with −0.2 =< IPI <= +0.2) and a low range of plots irrigated early (an average of 14.19% of plots with IPI > +0.2) or lately (an average of 9.64% of plots with IPI < −0.2). Based on the obtained results, we can conclude that the 1st and the 2nd sowing scenarios (starting dates 1 November 2011 and 16 November 2011, respectively) are the best for a better optimization of the IPI indices and then a good irrigation scheduling optimization [56]. By comparing with the other four scenarios, these two early sowing scenarios have a high percentage of plots irrigated at the appropriate time (between 42.2% and 97.4%) and a low percentage of late irrigated plots (between 0.86% and 11.2%). The number of irrigation rounds can be reduced from five (for first to fourth scenario) to four (for both fifth and sixth scenarios). Conversely, Table 3 presents for the six scenarios the optimized irrigation and sowing windows and the simulated windows to stop irrigation rounds. The results showed a clear impact of the sowing window on the irrigation schedules during the agricultural season: (1) irrigation window is delayed depending on the sowing window, (2) the gap between successive irrigation rounds ranges from 13 to 16 days, the duration of the rounds ranges between 18 and 21 days and the number of irrigation rounds can be minimized when sowing dates are delayed (from five to four). The obtained values of the gap between two successive rounds, the duration and the number of irrigation rounds are consistent with the observed ones at the agricultural season 2011-2012.
Additionally, to evaluate the impact of the sowing date's optimization on irrigation scheduling, we compared the results of the IPI index distribution obtained by our proposed approach, with the results of a previous study conducted on the same irrigation district [56], but without optimizing the sowing date (i.e., using the applied sowing date at the 2011-2012 agricultural season). The comparison concerned only the third irrigation round of 2011-2012 agricultural season. Three classes (late irrigation, early irrigation and optimal irrigation) of IPI index distribution are presented in Figure 11, for each sowing scenario and compared with the IPI index distribution without sowing dates optimization (from the previous study). We notice that the proposed optimization approach greatly improves the number of plots irrigated at the appropriate time (optimal IPI increase from 40% without sowing date optimization to 82% with sowing optimization according to the second sowing scenario). Table 3. Optimized irrigation windows, sowing windows, stopping irrigation window, number and mean duration of irrigation round (IR) and the mean gap between two successive irrigation rounds.

Effect of Sowing Date Optimization on Water Requirements
To evaluate the effect of sowing optimization on crop water requirement, we compared the simulated water needs at each irrigation round and for each sowing scenario with the applied (observed) irrigation amounts at the agricultural season 2011-2012. The results of this comparison are presented in Figure 12 and summarized in Table 4. For all of the scenarios, the overall demand for irrigation water is greatly reduced compared to the observed case. This reduction varied from 48.54% to 35.04% and afterward, to 44.53% for the first, third and sixth scenarios, respectively. The statistical comparison of the performance of these six scenarios, using the Student's t-test [82], shows that there are two groups of identical scenarios at p > 0.05 (scenarios S1, S5 and S6 as a first group and S2, S3 and S4 as a second group).  Additionally, for all the scenarios, the variation in water demand from one plot to another is largely reduced; it is about ±3 mm for the 1st scenario against ±182 mm for the observed case. This result is important for irrigation equity between fields and between farmers [39]. Table 4 presents the amounts of water applied or simulated at each irrigation round. We notice that they increase from one irrigation round to another depending on the development of the root zone from one stage to another. The low values of the total amount of consumed water resources estimated for the six scenarios vary from 141 to 178 mm, but by adding the rainfall which was taken into account when optimizing the water requirements, we obtain a total volume varying from 362 to 399 mm, and a total amount of 495 mm for the observed irrigation, which is close to the requirements of wheat as defined in other studies [36,38,83].

Effect of Sowing Date Optimization on Wheat Yield
The last objective of our study is to evaluate the impact of the optimized sowing dates and consequently the optimized irrigations on wheat yields (biomass and grain). For this purpose, we compared the simulated wheat yields for each sowing scenario with the estimated ones at the real case. As presented in Section 2.6, two approaches were applied. The first one is based on the NDVI profiles using Equation (12). For each sowing scenario, the NDVI index is calculated based on the sowing dates. We also used the NDVI profiles derived from satellite images which reflect the real crop growth conditions to estimate the real obtained grain yield at the 2011-2012 agricultural season [78]. We remind here that this first approach is used to identify the best spatiotemporal repartition of sowing date corresponding the optimal grain yield. The second approach is based on the AquaCrop model considering crop growth, the sowing dates and the irrigation schedules (times and amounts). In this case, the evaluation of the proposed approach concerns the produced biomass and grain. Figures 13 and 14 present the results of the first and second approach respectively.
The results showed that the wheat yields (grain and biomass) are maximal for the early sowing dates and decrease for the late sowing scenarios. Except for scenario 6, the obtained yields (biomass and grain) are clearly high compared to the observed ones.
Otherwise, for scenario 1, the grain yield obtained by the NDVI-based approach ( Figure 13) and by Aquacrop ( Figure 14) reach 7.4 and 7.8 t/ha, respectively, and the biomass estimated by Aquacrop is 15.5 t/ha ( Figure 15). These yields are comparable to the potential yields that can be obtained under the agro-environmental conditions of our study site [38,62,67]. This excellent result shows the validity of the proposed optimization approach.
The statistical comparison of yields (grain and biomass) estimated for the different scenarios, using the Student's t-test, shows that overall, the scenarios results, presented in Figures 13-15, are significantly different at p ≤ 0.05, except for scenarios S1 and S2 ( Figure 13) that are evaluated identical at p > 0.05.

Discussion
The main objective of this study is to evaluate the effects of optimized wheat sowing date, based on our proposed approach, on irrigation schedules, water requirements and yield.
Thus, we first illustrated the performance and the efficiency of our optimization approach, by optimizing irrigation schedules based on the IPI index derived from spatialized Ks coefficient, at the beginning of each irrigation round. This optimization leads to a high range of plots irrigated on time for each irrigation round (Figure 10), particularly for the first and the second sowing scenarios (starting dates 1 November 2011 and 16 November 2011, respectively), while promoting early irrigation over late irrigation because it avoids crop water stress and leads to increased yields [62,78]. Then, to evaluate the impact of sowing date's optimization on irrigation scheduling, we compared the obtained results with another approach based only on IPI index optimization without optimizing the sowing date. Results show that with our proposed approach including sowing date optimization we can reach a high percentage of plots irrigated at the appropriate time ( Figure 11). This revealed that the optimization of the spatiotemporal distribution of sowing date over the irrigation scheme has a significant impact on the irrigation schedules and their optimization process, thereby increasing the adequacy of the supply of irrigation supply to meet crop water needs.
Conversely, to evaluate the effect of sowing optimization on crop water requirement, we compared the simulated water needs with the observed irrigation amounts. Results show that with an early sowing dates we can reach a high reduction (up to −40%) of water resources (Table 4), and that with the scenario 3 (stating date 1 December 2011), we have a relatively high demand, compared to the other scenarios, which is justified by the occurrence of the grain maturity stage in April, which is associated with high air temperatures and consequently high evapotranspiration [62]. This situation is much more critical for the three scenarios (4th to 6th), but the relative high air temperatures during all the phenological stages of wheat greatly reduces the crop cycle [15] and therefore a decrease in water demand. However, this climate condition for these three scenarios impacted significantly the wheat yield in biomass and grain [60,63].
The evaluation of the effect of sowing date optimization on wheat yields (biomass and grain) showed that with the early sowing scenarios we can increase the wheat yields compared to the late sowing scenarios ( Figure 13). This result is in agreement with previous research results conducted for the wheat crop in the same study area [36,63]. In this semiarid zone, early sowing allows wheat to have a relatively long cycle allowing crops to absorb a high amount of solar radiation, which consequently leads to a high yield [63,84]. Moreover, late sowing dates correspond to pollination and grain filling stages in April month characterized by high air temperatures. This increases the risk of heat stress, which significantly affects the grain weight and the number of grains per ear [85].
Thus, the promising outputs of this study highlighted the interest of the adoption of the suggested approach by the managers of gravity irrigation schemes as an efficient tool for optimal sowing date shifting and optimal irrigation scheduling, which can be considered as key adaptation strategies to ensure agricultural water sustainability.

Conclusions
In order to contribute to the adaptation of agricultural practices to water scarcity conditions and to meet wheat production needs in the context of climate change, the effect of shifting sowing calendars was considered in this study. An approach based on remote sensing and on the CMA-ES evolutionary algorithm was developed to optimize the spatiotemporal sowing date distribution of wheat crop in a large irrigation scheme located in Haouz plain (Morocco), considering the different constraints related to the irrigation network as well as the climatic conditions and available water resources. Different sowing scenarios were simulated to evaluate the performance of our approach and to evaluate the effect of the sowing date optimization on irrigation management and crops yield, using weather data, field data and irrigation practices data collected in the study area during the 2011-2012 agricultural season. The results showed a high efficiency of the proposed approach by increasing yield, reducing irrigation water demand and rationalizing irrigation scheduling, compared to observed data.
Conversely, the simulation results showed a significant impact of the sowing window optimization on the irrigation schedules during the agricultural season, mainly the obtained irrigation windows, the gaps between successive irrigation rounds, and the duration and the number of irrigation rounds.
The effect of sowing optimization on crop water requirement was also evaluated by comparing the estimated water needs derived from crop evapotranspiration (ET c ) with the observed irrigation amounts applied by farmers. Results showed that a reduction of up to −40% of water needs can be achieved. The simulated wheat yields after sowing dates optimization were compared with the estimated ones at the real case by applying two different approaches: NDVI-based approach and AquaCrop model. The results revealed that the early sowing scenarios lead to higher wheat yields compared to the late sowing scenarios for both approaches with a homogeneous distribution over the plots compared to the observed one, which implies a better equity among farmers.
Further simulations and validation in different agricultural seasons and different climate conditions and geographic regions, while considering more complicated cropping patterns, are needed to improve the proposed approach.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.