A Simple Light-Use-Efficiency Model to Estimate Wheat Yield in the Semi-Arid Areas

In this study, a simple model, based on a light-use-efficiency model, was developed in order to estimate growth and yield of the irrigated winter wheat under semi-arid conditions. The originality of the proposed method consists in (1) the modifying of the expression of the conversion coefficient (εconv) by integrating an appropriate stress threshold (ksconv) for triggering irrigation, (2) the substitution of the product of the two maximum coefficients of interception (εimax) and conversion (εconv_max) by a single parameter εmax, (3) the modeling of εmax as a function of the Cumulative Growing Degree Days (CGDD) since sowing date, and (4) the dynamic expression of the harvest index (HI) as a function of the CGDD and the final harvest index (HI0) depending on the maximum value of the Normalized Difference Vegetation Index (NDVI). The calibration and validation of the proposed model were performed based on the observations of wheat dry matter (DM) and grain yield (GY) which were collected on the R3 irrigated district of the Haouz plain (center of Morocco), during three agricultural seasons. Further, the outputs of the simple model were also evaluated against the AquaCrop model estimates. The model calibration allowed the parameterization of εmax in four periods according to the wheat phenological stages. By contrast, a linear evolution was sufficient to represent the relationship between HI and CGDD. For the model validation, the obtained results showed a good agreement between the estimated and observed values with a Root Mean Square Error (RMSE) of about 1.07 and 0.57 t/ha for DM and GY, respectively. These correspond to a relative RMSE of about 19% for DM and 20% for GY. Likewise, although of its simplicity, the accuracy of the proposed model seems to be comparable to that of the AquaCrop model. For GY, R2, and RMSE values were respectively 0.71 and 0.44 t/ha for the developed approach and 0.88 and 0.37 t/ha for AquaCrop. Thus, the proposed simple light-use-efficiency model can be considered as a useful tool to correctly reproduce DM and GY values.


Introduction
In the southern Mediterranean area, drought is increasing in frequency and intensity [1]. Associated with global climate change, this trend will likely be more evident in the future [2]. Drought damage to the agricultural sector affects both rural livelihoods and the national economy as a whole.
Morocco is an arid and semi-arid regions where water resources are very limited [3,4], and characterized by a high sensitivity to climatic conditions [5]. In this region, cereals are the dominant crop at the national scale [6]. With a high evaporative demand (about 1600 mm/year), frequent irrigation is crucially needed to reach potential growth and yield [7]. Indeed, about 85% of total annual water consumption is used for agriculture irrigation in the region [4]. In this context, accurate and real-time estimation of crop yield at different scales (e.g., local, regional, or national) is becoming increasingly important [8]. Additionally, rapid and precise acquisition of field information (e.g., land use, water status, and crop phenological stages) and crop production is important for formulating agricultural development planning and agricultural policy [9]. More specifically, to preserve water resources, the rational management of irrigation water is necessary [10,11].
Modelling can be a useful tool to assess and develop promising irrigation scheduling strategies under limited available water for increasing crop water productivity [12]. The challenge is to link meteorological conditions, soil properties, and field management controlling crop production, in order to better describe the processes leading to crop biomass and yield. In this context, crop modeling is an effective way to simulate the mechanisms and semi-empirical procedures associated to crop growth [13][14][15]. Many process-based crop models have been developed [16], such as AquaCrop [17], STICS [18], DSSAT [19], APSIM [20], SUCROS87 [21], WOFOST [22], and CERES [23]. Nevertheless, a common challenge for using these models is the necessity of knowing several input parameters describing the agro-environmental conditions, which are not always available. The diversities of wheat cultivars, water and fertilizer management, pests and diseases control, and field application (e.g., mulching) in most farmlands have greatly limited the application of such crop models at local and large-scale regions [13,14,24]. Thus, these models are used principally by scientists or engineers for academic research and cannot be used on a daily basis by actual field managers [25,26]. To address this issue, two types of models have been used. The first are the empirical models linking a vegetation index observed by remote sensing (like Normalised Difference Vegetation Index, NDVI), at a given time, to the final biomass and grain yield [27][28][29][30][31]. This kind of relationship is simple to develop but remains very sensitive to agro-environmental conditions and varies from one agricultural season to the other [29]. In addition, this approach integrates the effects of all types of stress into the vegetation index, which limits its use to drive the crop management and/or to interpret the crop yields [30,32]. The second alternative, which is based on light-use-efficiency [33], provides great potential for plant biomass estimation. This approach was implemented in several models such as CASA [24], LINTUL [34], and GLO-PEM [35]. In these models, biomass accumulation is transformed from the effective radiation intercepted by the crop canopy. Most of these models were widely used to estimate net primary productivity on either grassland or forest by using remote sensing data [36][37][38], but not as much for crops, because of less interpretation of processes of crop growth and field management [32]. Then, these models need to couple the light-use-efficiency theory with biomass partitioning [39], to simulate growth process and estimate yield for herbaceous crops. In real conditions, the climatic and soil conditions can modify the temporal evolution of the conversion efficiency [32]. Any stress having an influence on the stomata opening leads to a decrease in crop transpiration, thus a decrease in photosynthetic activity which results to a reduction in the conversion efficiency [32]. The most important external factors that can cause this closure of the stomata are water content, nitrogen content, and temperature.
This study proposes an enhancement of the Monteith light-use-efficiency model [33], to estimate wheat production in semi-arid areas. More specifically, stress factors introduced by the shortage of water and temperature will be expressed in a conversion coefficient of light intercepted to the biomass. The manuscript is organized as follows: Section 2 provides a description of the study area, in-situ measurements and satellite imagery, and the proposed approach to estimate wheat production.
Agronomy 2020, 10, 1524 3 of 23 Section 3 presents the results and discussions about model calibration and validation. Finally, Section 4 highlights the main conclusions and suggests some perspectives.

Study Area
The experiments were carried out in the R3 irrigated district located 40 km east of Marrakech city ( Figure 1) during three growing seasons of winter wheat (2002-2003, 2008-2009, and 2012-2013). R3 is an irrigated area of about 2800 ha. The dominant crops are cereals, more wheat than barley [40]. Flood irrigation is widely practiced by the majority farmers in this area. The wheat is generally sown between mid-November and mid-January, depending on climatic conditions and the start of the rainfall season. The climate of this region is semi-arid, typically South-Mediterranean with high temperature in summer (38 • C, in July) and mild temperature in winter (5 • C, in February) with significant daily and monthly differences [41]. The average annual precipitation is about 250 mm, whereas the potential evapotranspiration (ET 0 ) is about 1600 mm/year [41].
Agronomy 2020, 10, x 3 of 23 area, in-situ measurements and satellite imagery, and the proposed approach to estimate wheat production. Section 3 presents the results and discussions about model calibration and validation. Finally, Section 4 highlights the main conclusions and suggests some perspectives.

Study Area
The experiments were carried out in the R3 irrigated district located 40 km east of Marrakech city ( Figure 1) during three growing seasons of winter wheat (2002-2003, 2008-2009, and 2012-2013). R3 is an irrigated area of about 2800 ha. The dominant crops are cereals, more wheat than barley [40]. Flood irrigation is widely practiced by the majority farmers in this area. The wheat is generally sown between mid-November and mid-January, depending on climatic conditions and the start of the rainfall season. The climate of this region is semi-arid, typically South-Mediterranean with high temperature in summer (38 °C, in July) and mild temperature in winter (5 °C, in February) with significant daily and monthly differences [41]. The average annual precipitation is about 250 mm, whereas the potential evapotranspiration (ET0) is about 1600 mm/year [41]. The R3 district is practically flat. The total area is divided into different blocks, themselves separated into six rectangular plots with an area ranging from 3.5 to 5 ha. The irrigation in this area is managed jointly by three farmer associations in collaboration with the local center of the Regional Office of Agricultural Development of Haouz (ORMVAH). They play a key role on irrigation rounds management and on preparation of the irrigation scheduling according to sowing dates and the availability of water in the dam. At the beginning of agricultural season, a global amount of water is assigned to the irrigation sector. Managers and users decide the number of irrigations rounds and quantities assigned for the irrigation. At each round, the farmers receive an amount of water according to the owned area, without taking into account the type of crops and their water requirements. Additionally, ground water may be used for irrigation in limited cases, with priority for orchards, forages (alfalfa), and vegetables, but can also be used for cereal in case of high shortage The R3 district is practically flat. The total area is divided into different blocks, themselves separated into six rectangular plots with an area ranging from 3.5 to 5 ha. The irrigation in this area is managed jointly by three farmer associations in collaboration with the local center of the Regional Office of Agricultural Development of Haouz (ORMVAH). They play a key role on irrigation rounds management and on preparation of the irrigation scheduling according to sowing dates and the availability of water in the dam. At the beginning of agricultural season, a global amount of water is assigned to the irrigation sector. Managers and users decide the number of irrigations rounds and quantities assigned for the irrigation. At each round, the farmers receive an amount of water according to the owned area, without taking into account the type of crops and their water requirements. Additionally, ground water may be used for irrigation in limited cases, with priority for orchards, forages (alfalfa), and vegetables, but can also be used for cereal in case of high shortage of dam water. This complexity is compounded by the spatial heterogeneity of the sowing dates, the size of tilled plots, and the quantity of nitrogen used [42]. The result is a high spatial and temporal heterogeneity of irrigation water requirement that is not considered in water attribution.

Field Data
The meteorological data has been controlled by an automatic weather station installed near the study area. This station was equipped with instruments to measure global solar radiation (R g ), wind speed, precipitation, air temperature (T a ), and relative humidity.
For monitoring vegetation growth, measurements of dry matter (DM) and grain yield (GY) were performed on ten fields of wheat, during the 2002/2003, 2008/2009, and 2012/2013 growing seasons ( Table 1). For more details about these experiments see Belaqziz et al. [30], Duchemin et al. [41], and Le Page et al. [43]. Along each season, vegetation samples were cut at ground level within squares of 0.5 × 0.5 m. On each sampling date, five samples were chosen randomly within each field for spatial representativeness issues. The fresh biomass was first weighed in the field by using a portable electric balance. Next, the various samples were placed in an oven at 85 • C for 72 h. Afterward, dry biomass was weighed again to determine the DM. The GY was also weighted at the start of the grain filling. The fields named C i (i = 1 to 4) were used for model calibration (Equations (12)-(15)). These fields were selected for their representativeness of the growth conditions for the three studied seasons (Table 1). However, the fields V i (i = 1 to 6) were used for model validation (see, Section 2.2.3). Additionally, for additional validation of the model against AquaCrop model estimates, we used the data (climate, irrigation, and nutrient), relating to 112 wheat fields conducted in the study site (R3) during the 2008-2009 season, to run the two models [30].

Satellite Data
In this study, a set of surface reflectance images in the infrared (0.73-1.11 µm) and red (0.58-0.68 µm) bands have been used to determine the Normalized Difference Vegetation Index (NDVI) and the vegetation fraction cover, with high spatial resolution (10 to 30 m), for three agricultural seasons ( Table 2). For those images, a radiometric calibration and atmospheric corrections were made based on the reflectance of invariant objects (bare soil, fallow, and houses) and were then transformed into NDVI maps.  [30]. Atmospheric correction was performed by the SMAC correction algorithm [44] using the aerosol optical depth measured at the Saada station near Marrakech city [45]. -During the 2012/2013 agricultural season, 18 images obtained from the SPOT4and Landsat 8 Operational Land Imager (OLI) sensors were used. These images were acquired during 5 months between 31 January and 15 June 2013 [43]. Radiometric correction was performed by the Multi-Sensor Atmospheric Correction Software-Prototype [45,46]. For both SPOT and Landsat satellite images, the geometric correction was performed using the Ground Control Points collected during the experiment [40], followed by the radiometric correction which was achieved in three steps. Firstly, we corrected a "reference" image from atmospheric effects using the SMAC correction algorithm and standard values of atmospheric components. Secondly, the radiometry of each image was homogenized against this reference image thanks to a set of reliable invariant features. This normalization was performed by applying linear relationships, between the digital numbers of raw images and the reflectance values of the reference image [40]. Thirdly, an additional linear correction has been applied between the satellite NDVI values and the NDVI ground measurements collected using the hand-held Cropscan Multispectral Radiometer. This inter-calibration ensures a maximal agreement between satellite and in situ NDVI values.

Proposed Model
The proposed approach is based on the light-use-efficiency model developed by Monteith which uses the three efficiencies method as a basis [33,47]. This model exploits an energy approach relating the production of dry matter to the amount of solar radiation received by the plant, assuming that the accumulation of dry matter is proportional to the accumulation of photosynthetically active radiation (PAR) absorbed by the plant.

Dry Matter
The light-use-efficiency model is an agro-meteorological crop model, which is used to estimate the dry matter by introducing the radiation intercepted by the crop during a growing period. A semi-empirical equation based on three efficiencies is then used to calculate the daily production of dry matter as follows: with ∆DM (g m −2 ) is the dry matter produced during a period ∆t. ε s , ε i and ε conv represent the climatic efficiency, the light interception efficiency and the conversion efficiency of the radiation absorbed by the vegetation into biomass (g MJ −1 ), respectively. R g is the cumulative daily value of incident global radiation (MJ m −2 ) received by a horizontal plane above the vegetation.
-Climatic efficiency, ε s ε s is the ratio between PAR and R g . PAR is the spectral band of solar radiation between 380 nm and 710 nm [48]. The value of ε s varies between 0.42 (for direct radiation) and 0.65 (for diffuse radiation) with an average of 0.5 [49]. Additionally, Monteith [33] used the value of 0.5 as an average of ε s for tropical and temperate regions. In the same range, Varlet-Grancher et al. [50], obtained an average value of ε s equal to 0.48 for the Mediterranean climate. They also demonstrated that ε s varies slightly depending on the solar elevation, cloud cover, and site latitude. In this study the value ε s = 0.48 was used.

-
Interception efficiency, ε i ε i is defined as the ratio between the PAR intercepted by the canopy (PARc) and the downward PAR (PAR i ). ε i varies between 0 in the absence of green vegetation, and 0.95 for highly developed vegetation with full photosynthetic activity [51]. The calculation of this parameter depends essentially on the vegetation index, such as NDVI [52]. Often ε i is linked to NDVI by a linear relationship [53]. In this study we used: NDVI min and NDVI max are the minimum and maximum NDVI values observed during the wheat season in the study site. They correspond to the bare soil (NDVI min = 0.14) and the totally covering vegetation for the entire time series of remote sensing data (NDVI max = 0.92), respectively [41]. NDVIn = NDVI−NDVI min NDVI max − NDVI min is the normalized NDVI, and ε imax is the maximum interception efficiency corresponding to the value of ε i when NDVI=NDVI max . ε imax is then considered as a local characteristic of wheat growth. It depends on the canopy high, leaf area density, and the plant geometry at full development stage [54,55]. These canopy characteristics affect the scattering radiation within the canopy which is not measured by NDVI. Then, the coefficient ε imax is considered as an empirical constant that must be determined under local conditions of wheat growth.
-Conversion efficiency, ε conv ε conv (g MJ −1 ) characterizes the ability of a canopy to convert intercepted radiation into biomass. Indeed, ε conv is defined as the ratio between the quantity of dry matter (DM) produced during a given period ∆t and the PARc absorbed during the same period. This efficiency can be obtained experimentally, however there are great variations between studies [56]. Gosse et al. [57], recommended average ε conv values of 1.93 g MJ −1 and 2.51 g MJ −1 for C3 and C4 crops, respectively. However, ε conv can reach around 3.5 g MJ −1 in optimal conditions. For wheat crops, the values of ε conv are set between 1. 80 [61], and 1.9 g MJ −1 [62]. In addition, ε conv declines under extreme agro-environmental conditions (water, nitrogen, temperature, and salinity stresses) [32,59].
The following equation is proposed to estimate ε conv : with ε conv max (t) is the maximum conversion efficiency (g MJ −1 ), depending on the agro-environmental conditions, variety, and crop phenological stages [32,57,[63][64][65]. K s conv is the water stress coefficient that affects vegetative production [42,66]. K T is the temperature stress coefficient. •

Water stress coefficient K s conv
To estimate K s conv , a simple formulation based on the results of Toumi et al. [66], reinforced by Hadria et al. [42], and Jackson et al. [67], has been proposed. These studies demonstrated that the effect of water stress coefficient K s (defined by the FAO as the ratio between actual and maximal evapotranspirations) on yield production remains insignificant (less than 4%) when its value is situated between 0.7 and 1. Furthermore, keeping K s above 0.7 requires a significant increase in irrigation water [42,66]. In order to minimize the wasting of irrigation water so as not to affect the yield, we proposed the following formulation for K s conv estimating: Agronomy 2020, 10, 1524 7 of 23 The water stress coefficient K s depends on the soil water availability. It is calculated from the water balance in the root zone following exactly the FAO-56 model detailed by Allen et al. [68]: Dr is the depletion of water in the root zone (mm) calculated by using daily rainfall (P) and irrigation (I). TAW is the total available water in the root zone (mm). RAW is the readily available water used by the plant (mm). Equation (5) is valid if Dr > RAW. When Dr < RAW, Ks is equal to 1. TAW and RAW are calculated as a function of the soil moistures at field capacity (θ f c ) and at wilting point (θ wp ).
where Z r is the thickness of the root zone (mm) and p is a crop parameter which depends on the evapotranspiration under standard conditions (ET c ). For winter wheat and most cereals, the recommended value is p = 0.55 when ET c is around 5 mm/day (FAO-56, Table 2). When ET c differs to 5 mm/day, p can be adjusted using the following approximation: ET c is estimated from the NDVI by using the relationship developed by Duchemin et al. [41], for wheat crop grown in the studied site (R3).

•
Temperature stress coefficient K T The thermal stress factor K T measures the restrictive effect of temperature on the conversion of intercepted solar radiation into plant biomass. Indeed, the low and high air temperatures, compared to the optimal temperature for growth and yield, decrease the rate of phytomass production [69]. The effect of temperature on the rate of production of phytomass is taken into account. For that reason, the daily average of air temperature (T a ) is integrated to an equation based on the optimal temperature (T opt ) (necessary for the plants growth) as well as on the two extreme values T min and T max (below and beyond which plant growth is assumed insignificant) as follows [69]: The values of T min , T opt and T max used in this study were 5, 26, and 33 • C, respectively [42,66]. Following the above development, the Equation (1) of dry matter becomes: Finally, in order to reduce the number of unknown parameters, we replaced the product (ε imax .ε conv max (t)) by a single parameter ε max (t). Then, Equation (10) is rewritten as: In conclusion, ∆DM is expressed by using a single parameter (ε max ), three indices (K s conv , K T and NDVI) and one meteorological variable (R g ). ε conv max (t) can then integrate the effects of nitrogen stress, the carbon assimilation rate of the leaves (which depends on the phenological stages) and the possible interactions between these two factors and water and temperatures stresses.
From the appearance of the ears, the simple calculation of the dry matter dynamic remains insufficient to monitor the grain filling. The latter is controlled by the process of partitioning the produced DM between grains and straw. The Harvest Index (HI) evaluates this partitioning [70].

Harvest Index, HI
HI is the ratio between the grain's weight and the total dry matter produced at the time t. In this study, the dynamic of HI is modeled in two periods: -Linear increase of HI which begins at flowering (time t start , after sowing date) and ends when HI reaches its final value HI 0 (time t end ). This type of correlation has been used by other studies [71][72][73]. Thus, HI is modeled by the following equation: -Constant value of HI (HI = HI 0 ). t start and t end and HI 0 are local parameters, which will be derived from the field measurements. The advantage of this approach lies in its simplicity which intrinsically combines: (i) The storage of a part of DM assimilation in the ear, and (ii) avoid the requirement of the number and weight of the grain/ear to predict the grain yield [23,[74][75][76]. In fact, HI integrates the effect of water stress over different periods of the crop development phases. This has been confirmed for different types of crops by several studies carried out on cereals in Africa [77][78][79], in North America [80], in Europe [81], and globally [82].
For the determination of the final harvest index HI 0 , we were inspired by the work carried out by Jianqiang et al. [28], Becker-Reshef et al. [29], and Belaqziz et al. [30], who implemented operational methods for forecasting grain yields based on NDVI data. HI 0 is calculated as follows: where HI 0max is the maximum value of HI 0 observed in the study area and ∆HI 0 its variation range. NDVI max max and NDVI max min are the maximum and the minimum values of the NDVI max over the study area, respectively. Then following this formulation (Equation (13)), HI 0 varies from (HI 0max − ∆HI 0 ) to HI 0max , depending on the values of NDVI max . Finally, the grain yield GY, at the time t, is deduced by: where DM(t) is the sum of ∆DM values from plant emergence until the time t. As a summary, the Table 3 groups together the various parameters of the presented simple model.

Model Calibration and Validation
The model calibration was performed in two stages: -Inverting the Equation (11) to determine the local values of the parameter ε max as follows: Initially the values of ε max will be calculated by Equation (15) fed by the observed data of the evolution of DM, R g , and NDVI for the fields named C i (i = 1 to 4, Table 1). In the second step, ε max is correlated to the cumulative growing degree days (CGDD) calculated from the sowing date. The CGDD values were daily calculated by subtracting the wheat base temperature (T min ) from the daily average of air temperature (T a ) [76]: For analyzing the obtained relationship ε max (CGDD), the Brower phenological scale is used. This scale, described in detail by Laguette et al. [56], presents a finer agronomic description of the cereal development cycle. It allows classifying visible phenological stages in 21 steps easy to interpret: from 1 to 5 for the vegetative phase, from 6 to 17 for reproduction period, and from 18 to 21 for grain ripening.
-Adjustment of Equation (12), by using data for the fields C i (i = 1 to 4), to determine the local values of t start , t end , and HI 0 .
After, the calibrated model is validated by using the observed data collected for six fields V i (i = 1 to 6, Table 1), especially the temporal evolution of biomass and grain weight. In addition, further validation has been done through the comparison with the final biomass and grain yield calculated by AquaCrop model [17,83], previously calibrated for the study area [66]. For that purpose, the concerned fields are V i (i = 1 to 6) and the 112 fields mentioned in Section 2.1.2.

Model Evaluation Metrics
The estimated DM and GY by the proposed model were compared to the observed values and those estimated by AquaCrop model using the determination coefficient (R 2 ) of the regression line and the Root Mean Square Error (RMSE): where Y i is the predicted value by the proposed model, X i is the observed value or that estimated by AquaCrop model, Y is the mean of Y i values, X is the mean of X i values, and n is the number of observation. Additionally, the performance of proposed model was compared to AquaCrop model by using the statistics of linear models presented in Coursol (1983) [84], which evaluate the statistical significance of the similarity between two regression lines (see Appendix A).

Calibration of ε max
As mentioned above, the calibration of ε max consists of modeling its temporal dynamic as a function of CGDD. This evolution of ε max is attributed to the seasonal change of conversion efficiency [32,63]. Indeed, Equation (15) is used to calculate the decadal values of ε max by using data (R g , T a , NDVI, DM) observed over four fields of wheat (C1-C4, Table 1). The obtained results, presented in Figure 2, show that four stages can be distinguished in the evolution of ε max (CGDD): Agronomy 2020, 10, x 10 of 23

Calibration of
As mentioned above, the calibration of ε consists of modeling its temporal dynamic as a function of CGDD. This evolution of ε is attributed to the seasonal change of conversion efficiency [32,63]. Indeed, Equation (15) is used to calculate the decadal values of ε by using data (Rg, Ta, NDVI, DM) observed over four fields of wheat (C1-C4, Table 1). The obtained results, presented in Figure 2, show that four stages can be distinguished in the evolution of ε (CGDD): Stage 1: ε is practically zero, from emergence to thermal time t1, Stage 2: Rapid increase in ε to its maximum value 1.15 g MJ −1 . This stage occurs between times t1 and t2, Stage 3: ε remains practically constant at its maximum value 1.15 g MJ −1 (the average value between t2 and t3), Stage 4: Almost linear decrease of ε , from its maximum value (in t3) to a low value (in t4).   Table 1). The stars represent the calculated values of ε max (Equation (15)) and the solid black lines are their linear regressions. The green dotted line shows phenologic development along the Brower scale shown on the right y-axis. t 1 -t 4 indicate cumulative growing degree days (CGDD) at the growth stages t 1 = 230, t 2 = 260, t 3 = 1186, and t 4 = 1500 • C-day.
Stage 1: ε max is practically zero, from emergence to thermal time t 1 , Stage 2: Rapid increase in ε max to its maximum value 1.15 g MJ −1 . This stage occurs between times t 1 and t 2 , Stage 3: ε max remains practically constant at its maximum value 1.15 g MJ −1 (the average value between t 2 and t 3 ), Stage 4: Almost linear decrease of ε max , from its maximum value (in t 3 ) to a low value (in t 4 ). As mentioned above, to determine the meaning of articulation times (t 1 -t 4 ), we used the Brower scale of wheat development. Thus: I t 1 = 230 • C-day: This time corresponds to the stage 2 of the Brower scale. This shows that t 1 coincides perfectly with the start of tillering. This result is reinforced by Hadria et al. [42], who obtained tillering at a CGDD equivalent to 260 • C-day. This work was calibrated and validated by the STICS model on wheat crops over our study area R3. I t 2 = 620 • C-day: Is the thermal time corresponding to the start of the plateau phase. This time coincides with the stage 4 of the Brower scale wich corresponds to the end of tillering and the start of upstream stages. This time is characterized by the transition from leaf production to that of spikelets [85]. The obtained t 2 value is slightly lower than the CGDD = 696 • C-day estimated by Toumi et al. [66], for the start of the maximum of cover fraction.  According to this description, we note that phase 1 (from emergence to t 1 ) is characterized by very low values of ε max . Thus, for this period ε max is assumed to be zero. Additionally, for phase 3 we have assigned to ε max the average of the values calculated between t 2 and t 3 . However, for the two phases 2 and 4, the equations of ε max were derived based on the linear regression. The corresponding coefficients of determination (R 2 ) are 0.96 (n = 15) and 0.91 (n = 7) for phases 2 and 4, respectively. Thus, according to the results of this calibration, ε max is parameterized as follow: This expression reflects clearly the impact of the phenological stages on the coefficient ε max . In this context, Arkebauer et al. [32], and Asrar et al. [63] have shown that the values of ε conv vary significantly with the wheat phenological stages. Generally, the values of ε i and ε conv are low during the initial phase [63]. However, during the growing stage (between t 1 and t 2 ) the increase in ε max may be explained by the gradual increase of the carbon assimilation rate of the leaves required for the canopy expansion and for DM storage within the leaves fully emerged [17,18]. During the plateau phase (from t 2 to t 3 ), there is no canopy expansion, but the conversion coefficient remains maximum, because of grain filling which requires a high photosynthetic activity [17,18]. During the last phase (from t 3 to t 4 ), which coincides with the grain maturity stage and progressive senescence of the vegetation, the carbon assimilation rate of the leaves decreases [17,63]. This can explain the linear decrease in ε max .
Moreover, the conversion coefficient ε conv is also affected by vegetation water status. In this context, Muchow and Davis [86], have shown that ε conv is more affected by crop water stress than ε i . Additionally, for cereals, Muchow [87], showed that the plant water status may impact the interception coefficient (ε i ) during the first six weeks. After this period, ε conv is the most affected. In the same context, Asrar et al. [63] showed that just after irrigation, in the middle of the wheat season, ε conv increased from 3 to 3.22 g MJ −1 . Figure 3 shows the evolution of HI as a function of CGDD obtained for the calibration fields C i (i = 1 to 4, Table 1). This trend is segmented into two parts. The first one is characterized by a linear increase between t start and t end which fits well with the grains filling and maturity stages. The linear adjustment of Equation (12) led to the values of t start = 750 • C-day and t end = 1313 • C-day. For the second part (after t end ), HI keeps constant at its final HI 0 , reached at the end of the wheat growing season. Following the data used in this calibration HI 0 = 0.50. This value is comparable to those obtained by other studies on wheat, such as Toumi et al. [66] in Morocco (HI 0 = 0.46), Moriondo et al. [88] in South America (HI 0 = 0.48), Jin et al. [89] in the northern plain of China (HI 0 = 0.46) and Hammer and Muchow [90] in northern Europe (HI 0 = 0.55).

Calibration of HI
Agronomy 2020, 10, x 12 of 23 Figure 3. Variation of the harvest index (HI) obtained for the calibration fields Ci ( Table 1). The solid line represents the adjustment of Equation (12).
Consequently, for the study region, the Equation (12) becomes: The value 563 °C-day is the difference between tend and tstart. Furthermore, HI0 can exhibit significant annual variability [28,29], which essentially results in the combination of the agro-environmental and climatic effects [6]. Hence the interest of the equation 13, which calculates HI0 according to its maximum value observed in the study region (HI0max) and its variation interval ΔHI0. The determination of these two characteristics (HI0max and ΔHI0) is made by using data from 10 Figure 4). It can be seen, that the HI0max value observed in the study region is about 0.59 and the interval ΔHI0 is 0.15. The obtained value of HI0max is in the range (0.51-0.64) reported by Balaghi et al. [6], observed in other Moroccan regions. This study shows that HI0max depends on the climate, the sowing date, and the rainfall/irrigation of the agricultural season. These conditions affect the value of the field's NDVImax. In our case, depending on the NDVImax value, HI0 of the field can vary from 0.44 (=0.59-0.15) to 0.59.  Table 1). The solid line represents the adjustment of Equation (12).
Consequently, for the study region, the Equation (12) becomes: The value 563 • C-day is the difference between t end and t start . Furthermore, HI 0 can exhibit significant annual variability [28,29], which essentially results in the combination of the agro-environmental and climatic effects [6]. Hence the interest of the equation 13, which calculates HI 0 according to its maximum value observed in the study region (HI 0max ) and its variation interval ∆HI 0 . The determination of these two characteristics (HI 0max and ∆HI 0 ) is made by using data from 10 Figure 4). It can be seen, that the HI 0max value observed in the study region is about 0.59 and the interval ∆HI 0 is 0.15. The obtained value of HI 0max is in the range (0.51-0.64) reported by Balaghi et al. [6], observed in other Moroccan regions. This study shows that HI 0max depends on the climate, the sowing date, and the rainfall/irrigation of the agricultural season. These conditions affect the value of the field's NDVI max . In our case, depending on the NDVI max value, HI 0 of the field can vary from 0.44 (=0.59-0.15) to 0.59.
Finally, for this proposed model, the two relationships ε max (CGDD) and HI(CGDD) (equations 19 and 20, respectively) represent the local agro-environmental characteristics of wheat production. They implicitly integrate the effects of different factors affecting crop development and yield which were not considered explicitly in the Equation (11), especially nutrient stress, soil quality, wheat variety, etc. Thus, these two parameters represent the calibration keys of the proposed method for other sites having other agro-environmental characteristics.
(Table 1, Figure 4). It can be seen, that the HI0max value observed in the study region is about 0.59 and the interval ΔHI0 is 0.15. The obtained value of HI0max is in the range (0.51-0.64) reported by Balaghi et al. [6], observed in other Moroccan regions. This study shows that HI0max depends on the climate, the sowing date, and the rainfall/irrigation of the agricultural season. These conditions affect the value of the field's NDVImax. In our case, depending on the NDVImax value, HI0 of the field can vary from 0.44 (=0.59-0.15) to 0.59.   Regarding the final DM and GY yields ( Figure 5, red circles), the performance of the proposed approach is very promising. The RMSE values are 1.07 and 0.54 t/ha for DM and GY, respectively. These values are acceptable by comparison to the observed average final yields: 5.65 t/ha for DM and 2.70 t/ha for GY. The accuracy of the simple proposed approach is practically similar to the one obtained by more complex models of crops development and yield. Indeed, for the study area (R3), the STICS model [75] reproduced the final biomass and grain yields with RMSE equal to 1.35 and 0.87 t/ha [42]. Additionally, using the AquaCrop model, Mkhabela and Bullock [91] estimated the wheat grain yield for five sites in Western Canada from 2003 to 2006, with an error of about 24%. Similarly, the WOFOST model estimated the final yield, for six European regions, with an average error of 13 and 19% for the DM and the GY, respectively [56].

Model Validation by Using the Observed Data
Furthermore, for the field V1 invaded by oats during the 2002-2003 season [41], the observed DM and GY (indicated by the blue oval marking in Figure 5) are clearly overestimated by the proposed approach, which simulates the absorbed solar radiation using remotely sensed NDVI. The same result was underlined by Hadria et al. [42], and Toumi et al. [66], who monitored the wheat yield in the same study area (R3) by using STICS and AquaCrop models, respectively. Otherwise, for the DM between 1.5 and 5 t/ha, corresponding mainly to the stages before maturity, the model overestimates the observations with an almost systematic difference of about 1 t/ha (Figure 5a). Fortunately, this overestimation of DM did not affect the estimations of grain yield (Figure 5b).
Overall, the yield predictions (DM and GY) by the proposed approach are very encouraging for the six validation fields. The observed differences between observed and simulated production can be explained by the difference between the calibration and validation fields in terms of soil texture [92] and the nutrient applied. Currently, the proposed model does not take explicitly into account the nutrient stress. The latter is considered implicitly through the calibrated agro-environmental characteristics ε max (Equation (19)) and HI (Equation (20)) overall for the studied area. Figure 5 displays the comparison between the simulated and the observed dry matter (DM) and grain yield (GY) during the wheat development, for the six validation fields (Vi), conducted during the 2002/2003, 2008/2009, and 2012/2013 wheat growing seasons. From this figure, it can be seen that there is a good agreement between the simulated and observed yields (DM and GY). Indeed, for DM ( Figure 5a, black circles), the values of R 2 and the slope of the linear regression line are 0.81 and 0.82, respectively. Additionally, the RMSE value is around 1.18 t/ha, which represents about 17% of the average value of DM observed. For the evolution of the grain filling (Figure 5b, black circles), the statistical values of the comparison between simulation and observation results are also satisfying. The values of R 2 , the slope, and the RMSE are 0.77, 0.94, and 0.53 t/ha, respectively.

Model Evalution against AquaCrop Model
In order to go further in the evaluation of the proposed approach, the estimated DM and GY by this approach were also compared with the AquaCrop model simulations on six fields V i (Table 1) monitored during three wheat seasons ( Figure 6). These results showed that, despite its simplicity, the statistical performance of the proposed approach is close to that of the AquaCrop model, which is much more complex and very greedy in input parameters [76]. Indeed, for the DM estimations, the slope, the intercept of the regression line, the determination coefficient (R 2 ) and the RMSE are respectively 0.   The differences obtained between the two model estimates may be due to their schemes used to represent the effect of water stress on the growth, development, and yield [93,94]. Indeed, the AquaCrop model describes finely the effects of the plant water status on DM development by using a dynamic function of water stress taking into account the expansion of the leaves, the stomata closing and the senescence for the crop [17,76]. Thus, in AquaCrop, under limited water conditions, the plant transpiration is directly affected. On the other hand, the developed model used a simple relation (Equation (3)) to monitor the water stress effects on the conversion coefficient ε conv . This effect has been alleviated by the proposal of a new formulation of the stress coefficient (K s_conv , Equation (4)) to optimize the irrigation water supply [42,66]. This may explain the overestimation obtained in some fields, by the proposed approach, of the yields calculated by AquaCrop, especially for DM lower than 3 t/ha (Figure 7a) and GY lower than 1.5 t/ha (Figure 7b).
Additionally, for AquaCrop, carbon sequestration and nitrogen absorption are mainly linked to radiation and nitrogen limitations [17,76]. However, for the proposed model, the effect of nitrogen is incorporated into local ε max calibration (Equation (19)). This correlation represents the integrating effect of the regional agro-environmental conditions (e.g., soil fertility and salinity) on crop yield. The differences obtained between the two model estimates may be due to their schemes used to represent the effect of water stress on the growth, development, and yield [93,94]. Indeed, the AquaCrop model describes finely the effects of the plant water status on DM development by using a dynamic function of water stress taking into account the expansion of the leaves, the stomata closing and the senescence for the crop [17,76]. Thus, in AquaCrop, under limited water conditions, the plant transpiration is directly affected. On the other hand, the developed model used a simple relation (Equation (3)) to monitor the water stress effects on the conversion coefficient ɛconv. This effect has been alleviated by the proposal of a new formulation of the stress coefficient (Ks_conv, Equation (4)) to optimize the irrigation water supply [42,66]. This may explain the overestimation obtained in some fields, by the proposed approach, of the yields calculated by AquaCrop, especially for DM lower than 3 t/ha (Figure 7a) and GY lower than 1.5 t/ha (Figure 7b).
Additionally, for AquaCrop, carbon sequestration and nitrogen absorption are mainly linked to radiation and nitrogen limitations [17,76]. However, for the proposed model, the effect of nitrogen is incorporated into local εmax calibration (Equation (19)). This correlation represents the integrating effect of the regional agro-environmental conditions (e.g., soil fertility and salinity) on crop yield.

Conclusions
The main objective of this work was to present a simple and enhanced light-use-efficiency model for monitoring the wheat crop yield in the semi-arid area. The model links the biomass production to the solar radiation received by the plant, assuming that the increase of the biomass produced is proportional to the accumulation of photosynthetically active radiation absorbed. The originality of the method consists in: (1) The expression of the conversion coefficient (ε conv ) by considering an appropriate stress threshold (k sconv ) for triggering irrigation, (2) the substitution of the product of the two maximum coefficients of interception (ε imax ) and conversion (ε conv_max ) by a single parameter ε max , (3) the modeling of ε max as a function of the Cumulative Growing Degree Days (CGDD) since sowing date, and (4) the dynamic expression of the harvest index (HI) as a function of the CGDD and the final harvest index HI 0 depending on the maximum values of Normalized Difference Vegetation Index (NDVI) in the studied area.
The model calibration reflects the impact of the phenological stages and water stress on the ε max coefficient. Thus, the expression of ε max as a function of CGDD has been segmented into four parts. However, a linear evolution was sufficient to represent the evolution of HI according to the CGDD up to its final value (HI 0 = 0.50) obtained in 1313 • C-day. Thus, the ε max and HI equations are considered representative of the agro-environmental conditions of the Haouz plain. Therefore, they are considered the calibrating keys of the proposed approach in other different agro-environmental situations.
The model validation was performed against dry matter (DM) and grain yield (GY) observed during three wheat growing seasons and calculated by AquaCrop (more complex model). The obtained results reflect a satisfactory ability of the proposed approach to reproduce DM and GY. The values of R 2 and RMSE are of about 0.81 and 1.18 t/ha and 0.77 and 0.53 t/ha for monitoring the observed dynamics of DM and GY, respectively. For the final yield, the R 2 and RMSE are 0.8 and 1.07 t/ha and 0.81 and 0.57 t/ha for DM and GY, respectively. Against Aquacrop, the performance of the proposed model is very encouraging. Indeed, for DM, the R 2 and RMSE values are respectively 0.8 and 0.67 t/ha for the proposed model and 0.84 and 0.44 t/ha for AquaCrop. For the GY estimate, the values of these statistical metrics are respectively 0.71 and 0.37 t/ha for the developed approach and 0.93 and 0.26 t/ha for AquaCrop.
The integration of satellite data to the proposed model seems to be one of the most appropriate quantitative analysis methodologies that may be adopted for yield estimation at large spatial scales [95]. Indeed, the recent advent of new satellite sensors with both high spatial (less than 10 m), temporal resolution (less than 1 week), and multi-bands (optic, thermal, and microwave) as well as the developed computer processing tools have opened up new perspectives for mapping crop type [95] and field-level crop productivity [96]. This will allow the daily monitoring of the stress coefficient and the determination of the irrigations carried out necessary for the establishment of yield maps at large scales. Funding: This research was conducted within the Joint International Laboratory-TREMA (https://www.lmitrema.ma/). It was supported by the projects: SAGESSE PPR/2015/48 "Système d'Aide à la décision pour la GEstion des reSSources en Eau", ERANETMED3-062 CHAAMS "global CHange: Assessment and Adaptation to Mediterranean region water Scarcity", ACCWA-823965/H2020-MSCA-RISE-2018 "Accounting for Climate Change in Water and Agriculture management" and PRIMA-S2-ALTOS-2018 "Managing water resources within Mediterranean agrosystems by Accounting for spatiaL sTructures and cOnnectivitieS".
where Y i1 and Y i2 are the predicted values, X i is the observed value, Y 1 and Y 2 are the predicted means, X is the measured mean, and n is the number of observation.
Assuming that the couples (â 1 ,b 1 ) and (â 2 ,b 2 ) are independent, the variance of the estimators is: Comparison of the two slopes (a 1 and a 2 ): The null hypothesis (often denoted H 0 ) is: a 1 = a 2 . The observable value of the Fisher-Snedecor distribution (F obs ) is calculated as: Comparison of the two intercepts (b 1 and b 2 ): The null hypothesis (H 0 ) is: b 1 = b 2 . The observable value of the Fisher-Snedecor distribution (F obs ) is calculated as: Finally, for each comparison (of the two slopes or the two intercepts), the null hypothesis H 0 is accepted if F obs is less than the theoretical value of the Fisher-Snedecor distribution F th (1, n − 2) that exists in tabular form for different values of the accepted error (α). For example, if α = 1% , 1%, or 5% the hypothesis H 0 is accepted with an accuracy often denoted p < 0.001, p < 0.01 or p < 0.05, respectively.