Linkages between Rainfed Cereal Production and Agricultural Drought through Remote Sensing Indices and a Land Data Assimilation System: A Case Study in Morocco

: In Morocco, cereal production shows high interannual variability due to uncertain rainfall and recurrent drought periods. Considering the socioeconomic importance of cereal for the country, there is a serious need to characterize the impact of drought on cereal yields. In this study, drought is assessed through (1) indices derived from remote sensing data (the vegetation condition index (VCI), temperature condition index (TCI), vegetation health ind ex (VHI), soil moisture condition index (SMCI) and soil water index for di ﬀ erent soil layers (SWI)) and (2) key land surface variables (Land Area Index (LAI), soil moisture (SM) at di ﬀ erent depths, soil evaporation and plant transpiration) from a Land Data Assimilation System (LDAS) over 2000–2017. A lagged correlation analysis was conducted to assess the relationships between the drought indices and cereal yield at monthly time scales. The VCI and LAI around the heading stage (March-April) are highly linked to yield for all provinces (R = 0.94 for the Khemisset province), while a high link for TCI occurs during the development stage in January-February (R = 0.83 for the Beni Mellal province). Interestingly, indices related to soil moisture in the superﬁcial soil layer are correlated with yield earlier in the season around the emergence stage (December). The results demonstrate the clear added value of using an LDAS compared with using a remote sensing product alone, particularly concerning the soil moisture in the root-zone, considered a key variable for yield production, that is not directly observable from space. The time scale of integration is also discussed. By integrating the indices on the main phenological stages of wheat using a dynamic threshold approach instead of the monthly time scale, the correlation between indices and yield increased by up to 14%. In addition, the contributions of VCI and TCI to VHI were optimized by using yield anomalies as proxies for drought. This study opens perspectives for the development of drought early warning systems in Morocco and over North Africa, as well as for seasonal crop yield


Introduction
Agronomical drought, which is mainly caused by a deficit in water supply, is a major natural hazard to rainfed agricultural systems [1][2][3], causing up to 80% of yield losses worldwide [4]. In the semiarid regions of the south Mediterranean, 80% of arable lands are rainfed, enhancing the vulnerability of agriculture to drought events [5][6][7]. As a consequence, agricultural production, and cereal yields in particular, which are critical for food security [8], shows a very high interannual variability in these regions owing to uncertain rainfall and drought periods [9][10][11]. Drought frequencies have already risen in the Mediterranean area. In Morocco, the frequency changed from one event every 10 years at the beginning of the 20th century to 5 or 6 events every 10 years at the beginning of the 21st century [12,13]. In addition, climate change is expected to increase the frequency and severity of drought in Europe and Mediterranean countries [14][15][16][17][18]. For some countries such as Morocco, national cereal production does not cover the country's needs due to low productivity, even during sufficiently rainy years with a coverage rate for the country's needs ranging between 30 and 80% [19]. The first step towards drought mitigation is drought monitoring; thus, it is highly important to develop an early warning system for drought monitoring and to evaluate its linkages with yields regarding the dramatic socioeconomic consequences of a failure of cereal production [20]. Wilhite and Glantz (1985) classified drought by characterizing the cascade of impacts of a rainfall deficit (i.e., the meteorological drought) on water resources (hydrological drought), on the agricultural sector (agricultural drought) and on the country's economy (socioeconomic drought) [21]. Agricultural drought is related to a persistent lack of moisture in the root zone [22], affecting crop health and production. Thus, assessments of agricultural droughts should rely on a complex combination of variables related to water inputs, climate factors impacting the atmospheric evaporative demand, such as temperature or incoming radiation, soil factors and factors related to crops such as its leaf area index [23]. In addition, the timing of drought occurrence during the crop season, potentially associated with temperature and biotic stresses, is of prime importance to determine the impact on yields [24], because even a short water deficit occurs during the most sensitive phenological stages of cereals can drastically decrease production [25,26]. Thus, an early warning system for drought in terms of anticipating wheat production should also be able to deliver timely and accurate information at specific growth stages throughout the crop season.
The first systems for drought monitoring were mainly based directly on rainfall amount [1] or on related indices, such as the Standardized Precipitation Index (SPI; McKee et al., 1993 [27]). Considering the complex features of agricultural drought already highlighted, several authors considered other variables for drought monitoring [10,[28][29][30]. In particular, air temperature is both a governing variable in the atmospheric evaporative demand and crop evapotranspiration and a limiting factor of wheat production when strong anomalies (negative during the early stage of winter or positive around the heading stages) negatively impact yield [31,32]. The same reason also motivated the development of an extended SPI called the Standardized Precipitation Evapotranspiration Index (SPEI; Vicente-Serrano et al., 2010 [33]) to consider the effect of atmospheric evaporative demand. Nevertheless, considering that the lack of station data is often an obstacle in semiarid areas, several drought indicators derived from remote sensing observations, which are now freely available on a global scale, were quickly considered to monitor the severity and intensity of drought conditions [34][35][36][37]. Drought might affect vegetation and soil properties, such as vegetation biomass and development, canopy and soil temperature and soil moisture [38]. Remote sensing can provide information on vegetation vigor thanks to vegetation indices acquired in the optical domain, such as the normalized difference vegetation index (NDVI), land surface temperature (LST) derived from thermal infrared data and soil moisture (SM) retrieved from active or passive microwave observations. Most of the remote sensing drought indices are based on normalized anomalies of these products. For instance,  developed the vegetation condition index (VCI) based on NDVI and the temperature conditions index (TCI) based on LST, which have been shown to be two useful tools for monitoring drought on regional or global scales [39]. VCI was also shown to be strongly correlated with crop yield [40][41][42][43]. In addition, some studies found that a joint consideration of VCI and TCI together through the vegetation health index (VHI; computed as a weighted sum of TCI and VCI) was better suited for drought monitoring than considering both indices separately [39,42,44]. As there is no a priori knowledge of the actual contribution of the chlorophyll tissues and of the temperature conditions for the vegetation health in a given region, the value of α as proposed by  was 0.7, as VCI already considers the temperature effect on vegetation, but, more recently, a value of 0.5 was used by several authors, including Kogan (1997). Several recent studies questioned the value of 0.5 and proposed some new approaches for optimizing the values of this index [45][46][47]. In addition to temperature and vegetation conditions, some indices related to the SM conditions have also been proposed, as several surface SM products, mainly derived from microwave sensors with several dedicated missions (Soil Moisture and Ocean Salinity (SMOS), Advanced Microwave Scanning Radiometer (AMSR), and Soil Moisture Active Passive (SMAP)), are now available. The soil moisture condition index (SMCI) was developed by Zhang and Jia (2013) to describe drought conditions from a soil moisture perspective, which does not consider meteorological factors [48]. Table 1 displays a non-exhaustive list of studies based on different drought indices.
Nonetheless, the shallow sensing depth and the uncertain accuracy of currently available satellite SM retrievals fostered the integration of land surface models and surface SM observations or the leaf area index from satellites using data assimilation techniques to obtain more accurate root-zone SM estimates. Data assimilation techniques allow for the integration of spatially and temporally observed information into land surface models (LSMs) in a coherent way [49,50]. The Land Data Assimilation System (LDAS) refers to the framework where LSMs are driven by and/or ingest such observations, generating enhanced estimates of the land surface variables (LSVs) [51]. Recently, several operational land data assimilation systems (LDAS) have emerged, such as the coupled land vegetation LDAS (CLVLDAS, [52]) and the famine early warning systems network (FEWSNET) LDAS (FLDAS, [53]). SM data assimilation has also been actively investigated as a tool for improving operational drought monitoring [51,54,55]. Previous studies have demonstrated that the assimilation of surface SM retrievals can improve the estimation of root-zone SM by a land surface model [50,54,[56][57][58][59][60][61]. In particular, for agricultural drought monitoring, Bolten and Crow (2012) described the benefit of assimilating surface SM retrievals from the Advanced Microwave Scanning Radiometer for Earth Observing System (EOS; AMSR-E) into the modified Palmer SM model developed by the U.S. Department of Agriculture (USDA) [62]. In addition, they demonstrated that the assimilation of AMSR-E surface SM retrievals substantially improves the performance of a global drought monitoring system, especially for instrumented areas of the world where high-quality rainfall observations are unavailable. Additionally, Albergel et al. (2018) noted that the LDAS-monde (the land data assimilation system developed in the research department of the French meteorological service) is better able to characterize agricultural droughts than an open-loop counterpart (i.e., a model without any assimilation of satellite-derived measurements) [63].
Within this context, the main objective of this work is to evaluate the impact of agricultural drought on cereal yields in Morocco. Drought is assessed with extensively used remote sensing indices as well as outputs of an LDAS.

Study Area
Morocco is located at the southern edge of the mid-latitude storm track with a semiarid climate [76] (Figure 1). The climate is influenced by the Atlantic Ocean, the Mediterranean Sea and the Sahara Desert, together with the very steep orography in the Atlas region [77]. Most of the precipitation falls during winter and spring, which is from the beginning of November until the end of April [13]. Winter cereals occupy more than 55% of the country's agricultural areas; common and durum wheat account for approximately 75% of these cereals (MAPMDREF, 2019). This wheat is cultivated both in rainfed and irrigated fields, depending on access to water supply and climate conditions. Wheat production is mainly rainfed and represents more than 80% of total cereal production [78]. Cereals can be sown as early as November 1st if significant rainfall occurs, while a persistent drought at this time can delay seeding until 15 January. Late seeding is often associated with production loss through a decrease in wheat crop areas. Indeed, many farmers are accustomed to waiting for regular rainfall events to seed at the beginning of the season. Harvests occur in June-July of the following calendar year, starting in the south and steadily progressing towards the north and mountainous areas.

Data
The dataset is composed of the remote sensing drought index and the output of an LDAS described below and of the cereal yield statistics on the agricultural provincial scale. In addition, cumulative rainfall amounts and the SPEI, which is considered to be a reference for drought monitoring in several studies [79][80][81][82], were extracted from ERA5 reanalysis surface variables.

ERA5 Data and the SPEI
Data used for computing SPEI and rainfall amounts were downloaded from the ERA5 database, which is the fifth generation of European reanalyses produced by the ECMWF and a key element of the EU-funded Copernicus Climate Change Service (C3S) (https://www.ecmwf.int/en/forecasts/d atasets/reanalysis-datasets/era5). Important changes in ERA5 relative to the ERA-interim's former ECMWF atmospheric reanalysis include (i) a higher spatial and temporal resolution, as well as (ii) a more recent version of the ECMWF Earth system model physics and data assimilation system [83]. In this study, the SPEI was computed using the SPEI R package version 1.7 [84], and reference evapotranspiration (ET0) was computed using the FAO Penman-Monteith (FAO-PM) equation [85]. However, there is some debate over the use of ET0 or actual evapotranspiration (ETa) in calculating the SPEI. Begueria et al. (2014) showed the advantages of using ET0 instead of ETa for calculating the SPEI: inclusion of ET0 in the SPEI formulation is valid for both humid and arid climates, producing reliable estimations of drought severity [79]. In addition, the SPEI was computed on several time scales (SPEI 1, 3, 6, and 12 months) to consider the effect of an accumulating precipitation deficit and high evapotranspiration, which are critical parameters for crop growth.
In addition, Figure 2 shows the cumulative rainfall during the crop season (from November to May) derived from ERA5 data. The monthly mean temperature is displayed in Figure A1 in Appendix A. There are high spatial variations in rainfall and temperature, as noticed in previous studies [24]. Higher amounts of rainfall are observed in the northern part of the country: the cumulative rainfall reached more than 900 mm in the Taounat province. In contrast, the center of the country is characterized by low rainfall: Kelaa des Sraghna province exhibited approximately 340 mm less rainfall than Taounat (approximately 63%). Concerning temperature, areas located at high elevation (high Atlas, middle Atlas and Rif Mountains) are characterized by low temperatures compared with the other regions. In February, the mean air temperature in Khenifra and Haouz is approximately 6 • C, while it reached 15 • C in the Safi province, which is located along the Atlantic Ocean.

Yield Data
The main types of cropped cereals in Morocco are bread wheat, barley and durum wheat. Rainfed wheat yield data were acquired from the Economic Services of the Ministry of Agriculture. The dataset contains the crop production and harvested areas for 14 crop seasons between 2000 and 2017 on the administrative province scale (2001-2002, 2004-2005 and 2005-2006 were not available). These data are compiled from subprovince sample surveys and released in official documents as provincial averages.
The annual values of rainfed cereal yields in tons per hectare (t.ha −1 ) were calculated as the ratio between crop production (tons) and harvested area (hectares) for each province. In this study, focus is placed on the main cereal crop regions in Morocco. For this objective, the 15 most productive provinces (more than 90% of national production) were selected for analysis. Table 2 provides the average cropped areas, ratio between cropped area and total province area in percentage, average total production and average yields during the study period, coefficient of variation (standard deviation divided by the average, in percentage), cumulative total rainfall, and average air temperature during the crop season. Cereal yields are low compared with potential yields in all provinces, which can reach 7.5 to 8.0 t.ha −1 in rainfed areas when adequate technologies and practices are used [86,87]. This low yield is attributed to water stress, heat stress and poor soil, as well as to traditional agricultural practices that may lead to inappropriate crop rotation, insufficient nitrogen and fertilization and late sowing [88,89].
To enable interpretation of the results, groups of provinces with similar cereal yield interannual variability are identified through a classification using the kmeans based on the correlative distance (d c = 1−r 2 with r being the Pearson correlation coefficient between two time series [90]). d c was preferred to classic Euclidian or Manhattan distances, as time series with similar interannual patterns are sought.
The classification results are mapped in Figure 1. Time series corresponding to the standardized anomalies of yields for each group of provinces are presented in Figure 3. Four groups of provinces are isolated: - The first group named "Zone 1" covers the southern province of the study area: Settat, El Jadida, Safi, Kelaa des Sraghnas, Haouz and Beni Mellal provinces. Because of the strong north-south rainfall gradient in Morocco, the group 1 area is characterized by low rainfall and high temperature and is consequently considered to be a less productive area. Please note that the Safi province receives a large amount of rainfall (355 mm from November to May on average over the study period using the ERA5 dataset) as it is located along the Atlantic Ocean. However, yields are low, which may be due to poor soils. In addition, a large rainfed area in the Beni Mellal province is located in the foothills of the Atlas Mountains, which has colder conditions than in the plains; yields are also quite low (1.5 t/ha on average over 2000-2017). - The second group named "Zone 2" covers the Ben Slimane, Khenifra, El Hajeb, Khemisset and Meknes provinces. They are located in the center of the study region with relatively high rainfall conditions (approximately 600 mm on average during the crop season), showing higher yields than those for group 1, with an average yield of approximately 1.9 t/ha. - The third group named "Zone 3" covers the provinces of Taza, Taounat, and Fes. Zone 3 is located north of Morocco, and it is characterized by higher rainfall compared with the other groups (approximately 680 mm). The average yield in this zone is approximately 1.5 t/ha. - The last group named "Zone 4" includes only the province of Kenitra. This province is located in the northeastern part of the country. Zone 4 is dominated by irrigated areas where cereals represent more than 60% of crops, including the Gharb plain. Although rainfed yields are only considered in this study, it is likely that any rainfed fields neighboring irrigated areas benefit from the supplementary water supply. The combined high rainfall amount and high coastal air moisture may justify the high yield observed in this province compared with the other zones.

Remote Sensing-Based Drought Indices
In this study, five satellite-based drought indices were selected, and their potentials to detect agricultural drought were assessed. These indices are based on the normalized anomalies of NDVI, LST and SM. The dataset characteristics are summarized in Table 3.

The Vegetation Condition Index (VCI)
Changes in the NDVI index related to weather conditions (precipitation, temperature and wind) are lower than those related to the ecosystem (climate, soil, topography and vegetation type) for agricultural regions. Thus, drought impacts on crops are difficult to detect directly from NDVI data [39,64,91]. Therefore,  developed the VCI to control local differences in ecosystem productivity [92]: NDVI min and NDVI max are the maximum and minimum values, respectively, at each desired time step over the study period, and NDVI t represents the NDVI value in month "t". Thus, the VCI is a pixel-based normalization of the NDVI in which the short-term climate signal of the NDVI was filtered by separating it from the long-term ecological signal. In this study, the temporal composite series of MODIS NDVI (MOD13A2 collection 6; available at (https://lpdaac.usgs.gov/) from 2000 to 2017 with a spatial resolution of 1 km are used. The NDVI value is based on the maximum value composite method (MVC) [93] using a sliding 16-day period.

The Temperature Condition Index (TCI)
The TCI is based on LST. For the region where incoming radiation is not limiting vegetation growth, LST is related to the plant water stress when water is limited, as a higher LST should be related to a lack of available water for evaporation and transpiration and to more energy being used for sensible heat flux. Thus, the index is related to changes in vegetation health due to thermal effects. The TCI is computed as follows : LST max and LST min are the maximum and minimum values, respectively, at each desired time step over the study period, and LST t represents the LST in month "t". In contrast with NDVI, high LST during the growing season indicates unfavorable or drought conditions, while low LST indicates generally favorable conditions [39]. Therefore, the low TCI value indicates an unfavorable climatic condition (high temperature), while higher values mainly reflect favorable conditions (low temperature). The daily 1 km resolution LST (version 6) MOD11A1 product available through the U.S. Land Processes Distributed Active Archive Center (LP DAAC, https://lpdaac.usgs.gov/) from 2000 to 2017 is used in this study. The LST is derived from the generalized split-window algorithm [94].

The Vegetation Health Index (VHI)
To combine the weather and thermal condition impacts on vegetation response, several authors proposed the VHI as a combination of VCI and TCI: Parameter α describes the weights of VCI and TCI in VHI, and the α value ranges between 0 and 1. A value of α = 0.5 is considered in this study following Kogan (1997), and part of the discussion is dedicated to the assessment of optimal α values.

The Soil Water Index (SWI)
The SWI is related to the SM content. It is not considered to be a drought index; however, when integrated over time, the SWI is useful for monitoring the variation in the SM in rainfed agricultural areas. In particular, [95] highlighted the importance of SM for drought monitoring. However, SWI from the surface layer that can be inferred from microwave active (radar) and passive sensors has been shown to be well correlated with vertically integrated SM estimates derived from both in situ observations and water balance modeling [96][97][98]. Lagged correlations have also been highlighted with vegetation growth conditions on a monthly time scale [99].
The SWI algorithm is based on a two-layer infiltration model describing the relationship between the surface SM and vertical SM profile as a function of time [100]. The SWI is formulated as follows: where for t i ≤ t n , t n is the observation time of the current measurement and t i represents the observation times of the previous measurements (both given in Julian days). This model assumes that the water content of the deeper layer is controlled by the past moisture conditions in the surface layer and thus the precipitation history [96]. The parameter T, called the characteristic time length, represents the time scale of SM variations in units of time T = L/C, where L is the depth of the reservoir layer and C is an area-representative pseudodiffusivity constant [101]. An increased T value is either due to an increase in reservoir depth or a decreased pseudodiffusivity coefficient, which means that, for a fixed pseudodiffusivity constant, an increased T value represents a deeper soil layer [102], meaning that a high (low) T describes a deeper (shallower) soil layer [97]. Giving a general rule on how to translate a given T value to a certain soil depth is currently not possible since it strongly depends on the application and soil composition of the area of interest [101]. Ceballos et al. (2005) found that, for their study region in Spain, the best T-values were 40, 50 and 60 for layer depths of 0-25, 0-100 and 50-100 cm, respectively [96]. The surface SM product is retrieved from the C-band ASCAT radar onboard the METOP satellites [103]. The product is distributed by the Copernicus Global Land Service (http://land.coperni cus.eu/global/products/swi) with a daily time step and a 12.5 km resolution.

The Soil Moisture Condition Index (SMCI)
Agricultural drought is characterized by low SM levels that negatively affect agricultural production [104]. The SMCI [48] is a normalization of soil moisture values relative to the absolute maximum (SMmax) and the absolute minimum (SMmin) of the time series to obtain normalized SM ranging from 0 (very dry, unfavorable conditions) to 100 (very wet, favorable conditions). SMCI can describe drought conditions from a soil moisture perspective without considering meteorological factors, such as precipitation.
where SM min and SM max are the historical minimum and maximum values, respectively, at each desired time step over the study period, and SM t is the SM in month t. The data used for computing SMCI are derived from the European Space Agency Climate Change Initiative (European Space Agency (ESA) CCI) SM COMBINED version 04.2 (https://www.esa-soilmois ture-cci.org/), which merges SM observations from seven microwave radiometers (SMMR, SSM/I, TMI, ASMR-E, WindSat, AMSR2, SMOS) and four scatterometers (ERS-1 and 2 AMI and MetOp-A and B ASCAT) into a harmonious dataset [105]. This product had a daily temporal resolution and a 25 km spatial resolution.

Land Data Assimilation System (LDAS) Outputs
Within the SURFEX modeling platform of Météo-France [106], the LDAS [63,107,108] developed in the research department of Météo-France, the CNRM (Centre National de Recherches Météorologiques) permits integration of satellite products (leaf area index (LAI) and soil moisture (SM)) into the ISBA (Interaction between Soil Biosphere and Atmosphere) LSM [109,110] using a data assimilation scheme. The satellite product used in LDAS is the GEOV1 LAI, which is produced by the European Copernicus Global Land Service project (http://land.copernicus.eu/ global/). The LAI observations are retrieved from the SPOTVGT and PROBA-V (from 1999 to present) with a 1 km resolution and a temporal frequency of 10 days. Surface soil moisture was obtained from ESA CCI SM COMBINED (v4.7) as described above.
The algorithm of this data assimilation system permits the trajectory of the ISBA model to be optimally corrected each time an observation is available. Indeed, unlike the approaches conventionally used in crop models or carbon flux models, phenology in the CO2-responsive version of ISBA, ISBA-A-gs [110,111] is not constrained by degree-day submodels. The dynamic evolution of the vegetation biomass and LAI variables are entirely controlled by photosynthesis in response to atmospheric and climate conditions. Photosynthesis enables vegetation growth resulting from CO2 uptake and responds directly or indirectly to all atmospheric variables, as well as SM and canopy density.
The outputs of LDAS used in this study are the LAI, plant transpiration (Tr), soil evaporation (E), evapotranspiration (ETR) and SM at different depths: WG2 (0 to 4 cm), WG4 (10 to 20 cm), WG6 (40 to 60 cm) and WG8 (80 to 100 cm). In addition, the same variables derived from the open-loop run of the ISBA-A-gs model (without data assimilation) are considered to analyze the added value of data assimilation with regard to the model predictions alone. For all variables, a standardized time series anomaly was computed based on Equation (6).
where Y t is the standardized anomaly, X t is the actual value of the variable (LAI, Tr...) for month t, and X and σ are the mean and standard deviation of the time series.

Identification of Rainfed Cereal Areas
The identification of areas where rainfed cereal fields are cropped is important to assess the linkage between remote sensing indices and cereal production because rainfed crops are more sensitive to drought [2,3]. Genovese et al. (2001) found that the application of a cropland mask to select the NDVI value corresponding to wheat in a crop yield model significantly improved the accuracy of crop yield forecasts [112]. For this reason, the identification of rainfed cereal areas is carried out using a two-step approach: (a) ECOCLIMAP-II land cover [113,114] developed by CNRM at a 1 km resolution (https://opensour ce.umr-cnrm.fr/projects/ecoclimap/wiki) was used to select the pixels corresponding to C3 crops, including both irrigated and rainfed fields. (b) The land cover map at a 300 m resolution provided by the Climate Change Initiative land cover project of the ESA (https://www.esa-landcover-cci.org/) was used to distinguish between irrigated and rainfed crops.
The rainfed cereal areas were identified through the combination of ESA land cover and ECOCLIMAP-II.

Identification of Major Phenological Stage
The linkages between drought variables, including remote sensing indices and LDAS outputs, are first evaluated at the seasonal and monthly time scales, but the use of a time period that is in close agreement with crop functioning is discussed. For this objective, detection of the phenological scale is described in this section. Indeed, phenology detection is of prime importance for crop growth and yield [115]. A dynamic threshold method was used to detect phenology in this study following  [116][117][118]. As a preliminary step, the NDVI data were smoothed using the Savitzky-Golay filter to reduce noise. Then, the three main phenological stages [115] that can be easily determined from NDVI time series were detected from the smoothed NDVI curves as follows (see Figure A2 for an example of the smoothed NDVI profile, Settat province, season 2016/2017):

•
Emergence stage: The emergence stage starts when the NDVI reaches 30% of the difference between NDVI max and NDVI min. • Development stage: A drastic increase in the NDVI is observed at the start of the season. The development stage was defined as starting when the NDVI value reached 30% of the difference between the maximum and minimum of the NDVI values. The stage ends at the heading stage (see below).

•
Heading stage: The heading stage starts when the NDVI reaches its maximum value and ends at the harvesting stage.

Correlation Analysis between Drought Indices and Rainfed Cereal Yield
As reported by previous studies [2,72,73,[119][120][121], a correlation analysis between the indices and cereal yield was performed to investigate the link between drought indices and crop production. Pearson's correlation coefficient (R) between all indices and rainfed cereal yields was calculated on a monthly scale during the growing season, from November to May for the most productive province in Morocco for 2000-2017 using the following equation (Equation (7)). Pearson's correlation coefficient (R) represents the degree and direction of the linear regression between two continuous variables that are measured at the equal interval.
where x i and y i represent drought indices and the value of rainfed cereal yields, respectively, n is the number of samples, and x and y are the average values of x i and y i , respectively.

Seasonal Scale
The time series of remote sensing drought indices (VCI, TCI, VHI, and SMCI) averaged during the cereal growing season from November to May are superimposed with the yield anomalies and the standardized anomalies of the cumulative rainfall amount for each agricultural zone in Figure 3. Correlations between rainfall, drought indices and yields are reported in Table 4. to the presence of the largest irrigated area of Morocco (the Gharb region is more than 105,000 hectares) and/or high coastal air moisture. Indeed, the crop mask used to identify rainfed cereal suffers from uncertainties, and irrigated cereal may also be considered in the analysis. This finding may explain the significantly lower correlation coefficient with yields than those for the first 3 zones, as reported in Table 4. In addition to the low spatial resolution of some products, the correlation of yields and SMCI is particularly low (R = 0.14), as the ASCAT data used to compute SMCI have a 25 km resolution that may favor observations of mixed pixels containing both rainfed and irrigated fields. The 2007-2008 crop season for the latter zone exhibits a production deficit, while an average to above-normal production is recorded for the other zones. This result is because the Kenitra region faced severe flooding in autumn 2008 that hampered the crops in the region. The above-normal rainfall that watered the country at this time favored cereal production in regions that did not experience any flooding, such as the agricultural provinces located in zone 1 in the center of the country. Winter 2009/2010 is well known for an exceptional and persistent negative phase of the North Atlantic Oscillation [122], but unfortunately, it did not trigger above-normal rainfall in North Africa, as could be expected [77], leading to average production for the four zones. Finally, another striking feature is that all drought indices match well with each other while they depict very different crop characteristics, such as surface temperature (TCI), surface soil moisture (SMCI) or vegetation vigor (VCI). It is probable that the integration over the crop season and on the scale of the agricultural zone masks some specific seasonal and geographical features, which will be analyzed in the next section. Table 4. Correlation between remote sensing drought indices (VCI, TCI, SMCI and VHI), rainfall amounts during the crop season and crop yield anomalies for the four zones (*, **, *** = significant at 5/1/0.1% probability levels).

Monthly Scale
To further understand the linkage between drought and cereal crop yield, lagged correlations between drought indices and yields are carried out separately on the monthly time scale and for each agricultural province. Figure 3 displays the lagged correlation results from November to May. Only correlation coefficients significant at the 99% level are displayed. In addition, Figure A3 displays the lagged correlation between yields and cumulative rainfall and average temperature (both average since the beginning of the season in November and monthly values are considered) for comparison purposes, as those variables that are measured by meteorological stations in the synoptic network have been used for a long time for drought assessments.
The correlation between yields and rainfall and temperature ( Figure A3) corresponds to well-known patterns [9,24]: (1) a significant positive correlation with monthly rainfall early in the crop season in November or December. It is obvious that seeding in a well-watered soil will speed up plant emergence. In addition, the farmers practice opportunistic agriculture by seeding when the soil has been well watered thanks to well-distributed rainfall events at the beginning of the season; (2) a significant positive correlation with cumulative rainfall during the entire season with a maximum correlation occurring in February or in March around the vegetation peak; (3) for temperature, the correlation is much more scattered, but the monthly temperature appears to be significantly negatively correlated with yields in January during the tillering phase for most of the provinces, meaning that mild temperatures during this time are favorable for crop production due to the prevention of heat stress. Wheat is well known to be highly sensitive to air temperature during the developmental stages [123] when the crop is photosynthetically active and when high biomass accumulation occurs. Mild temperatures at this time favor long culms, large flag leaves and more potentially fertile florets in each spikelet [124].
For the VCI (Figure 4a), the correlation with cereal yield is high during the end of winter and in early spring (from February to April), which corresponds to development stages from late tillering to grain filling. The higher correlation for each province ranges between 0.61 and 0.94 (significant at the 99% level). Higher correlation coefficient values are observed for the provinces of Ben Slimane, Khenifra, Haouz, Settat, Khemisset and Taounat (0.89, 0.90, 0.92, 0.94, 0.93 and 0.85, respectively). These provinces are distributed within the first three zones, meaning that a differentiated analysis by zone is not relevant for the linkages between the VCI and grain yields. This finding could be attributed to the mechanisms explaining the linkages between the VCI and yields that are not related to climate or soil types. Indeed, the LAI is highly linked to the NDVI, and, thus, the LAI is also linked to the VCI. The LAI for cereals is a proxy for biomass accumulation; high aboveground biomass is also often related to high grain yield apart from late water or heat stresses that could prevent good grain filling, particularly during the pollination stage [125]. A high LAI also reduces soil evaporation through shielding and may improve the water use efficiency of the field. The lowest correlations are observed for Beni Mellal, Safi and Taza. The Safi province is known to have poor soil, which explains this area having the lowest yield of the 15 provinces (Table 2) with regard to cumulative rainfall amount. The only difference observed apart from the strength of the correlations is related to the timing at which the correlation is at the maximum, which slightly differs among provinces. Eighty percent of the 15 provinces exhibit a maximum correlation in March, which is the typical timing of the maximum green LAI for wheat in Morocco, while, for the Khenifra and Haouz provinces, the highest correlation was found in April. This finding may be attributed to slightly lower temperatures during the growing season, meaning low cumulative growing degree days and longer crop seasons. Indeed, these two provinces present the two lowest average temperatures during the crop season (see Table 2). Our results are in line with the work of Salazar et al. (2007) that was carried out on wheat cropped in Kansas, and these researchers highlighted a maximum correlation between wheat yield and the VCI during the maximum development of wheat [42]. Likewise, Unganai and Kogan (1998) in south Africa also found a similar timing of maximum correlation but for corn yield [41]. In addition, Zhang et al. (2017b) found a high correlation between VCI and wheat yield during the anthesis stage (February and March) in the Indo-Gangetic Plain (IGP) region of India [119].
Concerning TCI, as presented in Figure 4b, high correlation coefficients significant at the 99% level are also observed, ranging from 0.61 to 0.82. The positive correlation values found with the yields means that lower than normal LSTs are favorable for grain yield production during the development stages, including tillering and stem extension. In contrast with the VCI, the higher correlation values are found earlier in the season by an average of approximately 1 month, around January or February for 14 out of the 15 provinces. The main rationale of the TCI is that when water is a limiting factor, transpiration decreases drastically or even stops in the more extreme cases, and vegetation temperature, which may represent the LST when the canopy is adequately developed to shield the soil, increases. Stated differently, water stress should be avoided during wheat growth. In line with our results, several authors have already highlighted a positive correlation between cereal yields and the TCI during the growth stages of wheat in the United States [42], Mongolia [126] and Austral [75], and, as well as in Unganai and Kogan (1998) in South Africa for corn and others. The lowest correlation is found for the Kenitra province. For this province, irrigated cereal extending over more than 60% of the agricultural areas can be considered in the analysis because of the land use map uncertainty, as already highlighted. Irrigated cereals are well known to be less sensitive to air temperature than rainfed fields [127] because irrigation water inputs limit the water stress period and introduce thermal inertia to the soil. This finding could explain the lower correlation found for the Kenitra province.
The correlation between yield and VHI (Figure 4c) is significantly positive during all stages of development (tillering and stem extension) until anthesis, apart from the Kenitra province.
The correlation values with the VHI are generally stronger than those with the VCI and TCI taken separately. The most striking feature is that the spatiotemporal patterns are very similar to the correlation obtained with the VCI (Figure 4a), as already expressed by Bachmair et al. (2018) [128]. Likewise, the correlations exhibit a peak during late spring, which is in line with the results obtained with the VCI. This result means that the contribution of the VCI to the VHI dominates over the effect of the thermal conditions seen with the TCI. Several authors already underlined high positive correlation values with yields during the development stages of cereals, such as Ribeiro et al. (2019) on the Iberian Peninsula (IP), who found that the correlation between the cereal yield and VHI was the strongest during the heading stage of cereal [3].
The correlations between the SMCI and rainfed cereal yield are displayed in Figure 4d, and the correlations between the SWIs at different depths and yields are displayed in Figure 5. Please note that the SMCI is highly linked to the SWI that is representative of the upper soil layer (approximately 5 cm). In contrast with the other indices, the highest correlation values are between the SMCI and the SWI10 on the one hand, and on the other hand, yields are observed at the beginning of the growing season around seeding and/or the emergence stage, occurring as early as November for three of the provinces. Positive correlations that are significant at 99% are obtained with a maximum either in December (11 provinces for the SMCI) or January (three provinces for the SMCI); the only exception is Kenitra, which has already been shown to behave differently than the other provinces. The highest correlations of R = 0.88 for the SMCI and R = 0.95 for the SWI10 were observed for the Settat province. This positive correlation means that above-normal moisture conditions are favorable for high yields at this time. This finding is consistent with that of Modanesi et al. (2020), who found a high correlation between the standardized SM index (SSI) derived from the satellite SM observations of the ESA Climate Change Initiative (CCI) and the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2) SM dataset and wheat production in November and December in India [120]. Additionally, Zhang et al. (2017b) found a significant correlation between the SSI and wheat during the emergence stage (October and November) in the IGP, which is the northern region that separates the Peninsula from the Himalayan chain [119]. SM is another indicator that is highly suitable for assessing agricultural drought, as plants start to wilt when sufficient soil water is not available to meet the evapotranspiration demand. A deficit of SM during the development stages will affect crop growth, and finally, a decrease in final crop yield will be expected. This finding is well reflected by the correlation obtained with the SWI, which represents the deeper layers (SWI40 and SWI60). Indeed, as the roots develop during the crop season, deeper layers of soil undergo water extraction by the plant. The observed differences among provinces in terms of correlation could be attributed to differences in rainfall distribution or/and soil type, as computation of the SWI of the deeper layers was carried out using the time constant T for a specific type of soil. In line with our study, Modanesi et al. (2020) already showed that SM is a suitable indicator for representing agricultural drought, as SM correlates more closely with reduced crop yield than precipitation [120]. Zribi et al. (2010) illustrated a tight link between rainfall and the SWI in a semiarid area of Tunisia [99]. In addition, the researchers found that the SWI is strongly related to the NDVI during the early season (R = 0.81 in December). Likewise, a high correlation is also observed between the NDVI and the SWI in the more productive region of Northwest Africa [36]. Additionally, the SWI anomaly showed a strong correlation with the SPI index over the central region of Tunisia [129].

LDAS Outputs and Yields Time Series
In this section, the results of the relationships between anomalies of the LDAS outputs and the cereal yields are investigated ( Figure 6). In contrast with the remote sensing drought indices, the results are summarized by group instead of considering each province separately because correlation patterns similar to those of the remote sensing drought indices are observed. Indeed, the LAI provides information close to that of the VCI; SM in the root zone is related to the SWI of the deeper layer, while upper profile SM mainly governs soil evaporation; the TCI is closely related to the transpiration and evapotranspiration processes. For the LAI, the correlation reaches its maximum around the heading stage, as is the case for the VCI, which is in March for groups 1, 3 and 4 and in February for group 2. The strength of the correlation reached approximately 0.91, 0.89 and 0.90 for groups 1, 2 and 3, respectively. This result is in line with that of Sawada et al. (2019), who highlighted a similar correlation between the LAI analyzed within an LDAS system and wheat production at the national level in Morocco [130]. Transpiration (Tr) and cereal yield are significantly correlated from at approximately the maximum development stage and the heading stage. Higher transpiration is expected when the root zone is well watered and when the leaf area index is high, explaining the same timing (March) of the high correlation that was already observed with the VCI. The correlation also reaches higher values than the correlations with the VCI and the TCI, at least for zones 1, 2 and 3 (0.87, 0.78, and 0.88, during March, respectively). In contrast, monthly soil evaporation is significantly correlated at the start of the season, from sowing to emergence (November-December) for zones 1, 2 and 3. Again, zone 4 is different from the others, with a correlation peak occurring in February. Higher evaporation means higher soil water content, and this result is in line with the correlation pattern found with the soil moisture indices for the upper soil layers (SMCI and SWI10).
For the deeper layers (40 to 60 and 80 to 100 cm), the correlation is significantly positive later in the season, during the development stage of cereal, because of the progressive root development, as already highlighted for the SWI.

Case Study: 2015/2016
The 2015/2016 season is characterized by a high negative production anomaly that strongly affected the country's economy. In particular, more than 90% of the national territory of Morocco experienced droughts during December 2015 and January 2016 [34]. The cereal yield decrease reached 55%, 58%, 67%, and 13% of the long-term mean for zone 1, zone 2, zone 3 and zone 4, respectively. Figure 7 shows the 2015/2016 time series of drought indices and anomalies of LDAS outputs for all zones from September 2015 to July 2016. Figure A4 in Appendix A is the same for 2011-2012, which is also characterized by below normal yields. The first striking feature is a close correspondence between the drought indices and LDAS outputs, although a contrasted behavior can be highlighted between the first three zones and zone 4, as already discussed, particularly for the VCI and the LAI, as well as for the superficial soil moisture-related variables (SMCI and WG2). More discrepancies can be observed for root-zone soil moisture-related variables (SWI60 and WG6). This result is probably related to the uncertain soil depth for the SWI (the choice of the time scale T-Equation (4)-is related to the soil type). The TCI and Tr also follow the same dynamic, at least for zones 1 and 2. Zone 3 is the most affected zone, while zone 4 exhibited nearly normal production, probably because of the large irrigated areas, as already mentioned. This season faced a severe rainfall deficit (meteorological drought) as early as December in the four zones, as reflected by the below normal conditions in terms of the SMCI and the SPEI3, while the VCI and the SWI60 were still indicative of the average vegetation and root-zone moisture conditions at this time. Likewise, the air temperature was 2.7 • above average during the growing stages (MAPMDREF, 2016), which is in line with the TCI dynamic whose negative anomalies peak in January for all four zones. In January and February, during the development stage of cereal, the persistence of the precipitation deficit is associated with the continuation of adverse conditions in terms of the surface soil moisture described here by the SMCI. In March, during the heading stage of cereal, rainfall was still in short supply, and then, the VCI dropped for zone 1, zone 2 and zone 3, indicating that crops were severely affected by the continuous deficit of rainfall and above-normal temperatures. In contrast, the VCI for zone 4 still indicated good conditions even if precipitation remained in short supply, which shows the resilience of irrigation systems facing a rainfall deficit. Interestingly, 2011-2012 ( Figure A4) exhibited a different seasonal pattern of drought. Indeed, while the first part of the season until January exhibited average (for the SPEI) to above-normal conditions (for the VCI and the SMCI), the drastic drop in rainfall during the core of wheat development (February and March) had a major impact on grain yields, and the recovery of good soil moisture conditions occurred too late, probably after the grain filling stage in April. This contrasting behavior reveals that the different remote sensing drought indices related to vegetation, temperature and SM conditions can be used to anticipate drought impacts on yields in a timely manner throughout the crop season.

Discussion
Linkages between monthly drought indices and the outputs of land data assimilation with cereal yields have been highlighted at different times during wheat growth throughout the crop season. Several questions arising from these results are discussed below: (1) What is the added value of using LDAS to monitor drought with regard to the much higher expertise and computing time needed for the implementation of such a complex system? (2) Remote sensing drought indices are usually computed at a monthly time scale in the literature, but would an averaging time that is closer to the functioning of the crops, such as the phenological stages, be more relevant? (3) Finally, could an empirical choice of the α coefficient for VHI computation equal to 0.5 be revisited as a function of climate?

The Added Value of an LDAS
In terms of providing a drought monitoring dashboard to stakeholders, among which include managers on the catchment scale or ministries on the country scale, the complexity of implementation is a key question. Indeed, the LDAS, combining physically based LSMs with remote sensing variables through complex techniques of filtering, parameter identification and data assimilation, are difficult for local managers to implement without intense and continuous training sessions. This question is assessed in two different ways. As a preliminary step, the added value of assimilating remote sensing variables in a land surface model is assessed by comparing the correlation with the cereal yields of the LDAS output with and without data assimilation (the "open-loop" runs). Then, the correlation values between drought indices and LDAS outputs and cereal yields are compared at the group scale. Table 5 shows the correlation between the outputs of the LSM without assimilation (open loop) and with assimilation (analysis) and cereal yields. The analysis improves the correlation compared with the open loop for all variables. For example, the assimilation improved the correlation between the LAI and yield by approximately 6% for zones 1 and 3 and by 7% for zone 2. Concerning Tr, the correlation improved by approximately 1 to 2%, depending on the study zone. Conversely, the improvement is higher for the surface SM (WG2) when compared with the root-zone SM (WG6), as the variable is more difficult to predict because (1) the vertical profile of soil hydraulic properties is uncertain and (2) the analysis increments (i.e., correction of the variables) for the deeper layers may be lower than those for the superficial layer. Our results are in line with the work of Albergel et al. (2019), who highlighted the ability of the LDAS-monde to better characterize agricultural droughts than its open-loop counterpart over the continental United States of America [55]. Sawada et al. (2019) used the Coupled Land and Vegetation Data Assimilation System (CLVDAS) based on the EcoHydro-SiB land surface model to show that the simulated LAI at the end of the growing season is well correlated with wheat production over North Africa, including Morocco. Data assimilation also improved the skill of an LSM to reproduce the satellite-derived phenology so that the LDAS can reproduce the nationwide crop production of this water-limited region. Table 5. Correlation between output of the LSM with and without data assimilation and cereal yields at the group scale. Only the maximum correlation time already identified in Section 3 is displayed. All correlation coefficients are significant at the 99% level. Finally, Table 6 reports the correlation of the LDAS outputs and the remote sensing indices with cereal yields at the correlation peak. The LDAS and remote sensing indices provide comparable correlations with yields for vegetation characteristic variables (VCI and LAI) and to a lesser extent, for the TCI and Tr. In contrast, one striking feature is the added value of the LDAS outputs for variables related to the SM vertical profile. The correlations are always higher with the LDAS outputs than the indices derived from remote sensing products. This result could be related to the shielding effect of the canopy for the microwave observations used to derive the soil moisture-related indices. Indeed, the C-band backscattering coefficient acquired by the ASCAT sensor is used to compute the SWIs at different depths. The penetration depth of the C-band signal is low. In addition, Sigma0 is highly sensitive to surface roughness, which may change significantly during the crop season in response to soil workings and rainfall events. Brightness temperatures measured by microwave radiometers used to derive the SMCI are also highly sensitive to the canopy optical depth. Within this context, a mechanistic-based model predicting the processes involved in the SM dynamic (water interception by the canopy, surface layer evaporation, root water extraction, water diffusion between the different layers and deep drainage, etc.) appears to be better suited than using the remote sensing-derived SM products that could be contaminated with errors, especially when the crop canopy is fully developed.

Phenological Stages Versus Monthly Scale
Drought influences crop production differently depending on the development stage at the time of its occurrence [131][132][133]. The correlation of yields with the drought indices averaged over the main phenological stages described above is assessed in this section. Figure 8 shows the correlations between the VCI, TCI, SMCI and SWI60 and rainfed cereal yield at the key phenological stages; all correlation coefficients are significant at the 99% level. Table 7 displays the improvement in the maximum correlation with regard to the correlation obtained by integrating the drought indices on the monthly time scale. In line with the results described above, the correlation between the VCI and yields is the strongest during the heading stage, corresponding to the development peak, while the TCI is highly correlated during the development stage (Figure 8). The correlation by integrating across the main phenological stages is slightly higher than that using the monthly time scale, but the correlation improvement is limited, remaining below 15% for the VCI and the TCI and can even decrease by up to 14% (Kenitra province for the VCI). This finding means that, in view of implementing an early warning system of drought or for choosing predictors for the development of a model for the seasonal prediction of wheat yields, the monthly time scale is adequate regarding the additional data processing needed to identify the phenological stages. The behavior is quite similar for SWI60, except for some provinces in which the maximum correlation increase can reach 31%. For the index related to the soil moisture of the superficial layer SMCI, the maximum correlation is obtained around emergence during the early stage of wheat growth. It is striking that (1) the correlation is improved for all the provinces and (2) the correlation increase is significantly higher than that for the other indices and can reach 49%. The seeding dates are very variable from one province to another and even from field to field depending on the farming practices and the rainfall distribution. Thus, it is better suited to choose an integrating time in accordance with the observed emergence of the crop detected in this study using NDVI. It is more substantial at this time than during the following stages when the crops are well developed. The main conclusion is that the monthly time scale is adequate for the drought indices related to yields during the development and heading stages, while integration during the observed emergence stage should be preferred for the superficial moisture conditions impacting yields the earliest in the crop season.

Alpha Value for the Computation of VHI
It is well known that the LST and NDVI contributions through the α value (Equation (3)) to the VHI are not equal during the growing season or for different climate conditions [121]. As pointed out by [46], who used the SPEI as a reference drought index, it is possible and recommended to estimate an optimal value of α to assess the relative contributions of VCI and TCI for different regions during each phenological stage. The objective of this section is twofold: (1) to propose optimal values of α for each province and (2) to discuss the factors impacting these values. For this objective, α is optimized to maximize the VHI correlation with yields (α Yield ), which is expected to reflect the agricultural drought conditions during the entire growing season. The optimization method is a simple brute force approach to sampling the space of the α value (from 0 to 1) with a step of 0.001. For ease of interpretation, optimization at the scale of the phenological stages is preferred to the monthly time scale, and the development and heading stages are chosen, as the VCI and the TCI have been shown to be correlated with yields at these times of the crop season. As yield values are not available on the province scale everywhere in the world, optimization is also performed based on the SPEI 1, 3 and 6 from January to April (α SPEI ), which can be computed from widely available observations of temperature and precipitation following . Then, the correlation of the VHI with the α SPEI values and cereal yields is computed. The retained results are those maximizing the latter for the 15 provinces on average; the correlation corresponds to the SPEI6 in March around the vegetation peak, which is in agreement with Bento et al. (2020). Table 8 displays the correlation of the VHI with yields before and after optimization for both strategies (optimization with the yields and with the SPEI) and the optimized values of α based on the cereal yields. Only the best results (with the SPEI6 in March) are provided for optimization on the SPEI.
The results show that more weight should obviously be given to TCI during the development stage (optimized α < 0.5) and to VCI during the heading stage (optimized α > 0.5) for all provinces, which is in line with the correlation obtained between the TCI/VHI and the yields described above. The correlation between the yields and the VHI is obviously significantly improved when considering the optimized values based on yields, but the correlation values also increase when optimization is performed on the SPEI 6 in March. This result means that (1) the optimized value of α should be preferred to the widely used value corresponding to an equal weight to properly reflect the drought conditions of a specific region; and (2) it is possible to find a better value than the equal weight based on ancillary data when the yields are not available. This result has already been shown by  for Europe [46]. Interestingly, a close link is found between the α values and the seasonal rainfall amount during the development stage (positive correlation of 0.65) and to a lesser extent (r = 0.56) during the heading stage (no relationship with the average temperature was found; Table 2). Stated differently, the VCI contribution increases with the annual rainfall amount, which is in apparent conflict with the results of Bento et al. (2020), who showed that, on the global scale, the optimized α values increase with aridity [47]. Nevertheless, all ecoclimatic regions of the globe are of concern in the study of Bento et al. (2020), while only the semiarid and subhumid parts of Morocco are studied here. In addition, non-wheat areas are masked in this study, while mixed pixels are analyzed using the 8-km AVHRR data in Bento et al. (2020) [47]. For our interpretation, more weight should be given to the TCI, as the heat stress drastically impacts yields during the development and heading stages in the semiarid areas that are likely to face more severe temperatures than in the more humid regions located north of the country.

Conclusions
The main objective of this study was to assess the linkages between rainfed cereal production and agricultural drought in Morocco through remote sensing indices and an LDAS over the most productive province. This assessment was carried out using a simple lagged correlation analysis between the different variables and remote sensing indices and the wheat yields on the agricultural province scale. The main results, which are schematically summarized in Table 9, are biophysically sound and in line with the rich literature existing on this subject: the SMCI, SWI10, WG2 and Ev (soil evaporation), as proxies of surface soil moisture conditions, are closely linked to the yields for the early stage of wheat development, while the yields are correlated with moisture in the deeper layers later in the season as the roots develop (WG4 to WG8, SWI40 and SWI60); temperature conditions provided by the TCI and indirectly by the Tr mainly impact the yields during the development stage,  Table 9. Summary of the linkages between the remote sensing drought indices, the LDAS output variables and the cereal yields during the crop season (++ significant correlation, +++ highly significant correlation). Then, several questions arising from the obtained results are discussed. First, the added value of using an LDAS requiring strong expertise and computational power regarding simple drought indices freely available from remote sensing is analyzed to provide a drought monitoring dashboard to stakeholders. The main result is that a clear added value can be evidenced for all variables and in particular for the variables related to the SM profile that are not directly observable from space for the deeper layers. This result is interesting for stakeholders in the view of an operational monitoring of drought as, with the development of LDAS systems, surface analysis of moisture, temperature and vegetation conditions will be available freely in a near future, thus limiting the need for expertise and computing facilities. Indeed, vegetation is an interfering element for the retrieval of SM from active and microwave sensors. Thus, surface SM products have the worst quality when the vegetation is fully developed. Second, the typical integration time scale for the drought indices chosen within the literature is the monthly time scale. Therefore, the correlation by integrating over the main phenological stages that should be closer to plant functioning is computed. Higher correlations are obtained, but the correlation improvement is limited, apart from the index related to the soil moisture of the superficial layer SMCI, meaning that the monthly time scale is adequate for the drought indices related to the yields during the development and heading stages, while integration during the observed emergence stage are preferred for the superficial moisture conditions that most impact the yields early in the crop season. Finally, the VHI is an extensively used drought index combining the temperature conditions through the TCI and the vegetation conditions through the VCI, and the relative weight of both indices for computation of the VHI is questioned following Bento et al. (2020) [47]. The optimal value of the weight (α) in terms of mapping drought in relation to wheat production is computed by maximizing the correlation between the VHI and the yields. As an alternative way to identify α values better than the equal weight, α is also optimized, but it is optimized on the SPEI that can be computed from freely available data. The main results are that (1) the optimized value of α should be preferred to the widely used value corresponding to an equal weight to properly reflect the drought conditions of a specific region; and (2) it is possible to find a better value than the equal weight based on ancillary data when the yields are not available.

Emergence Tillering Elongation Booting Anthesis
This study opens doors for the development of early warning systems for agronomic drought in North Africa as well as for the seasonal forecasting of yields. To this objective, our on-going work is twofold: (1) estimating threshold values triggering a warning as well as quantifying the potential drop of yields as a function of the warning level; (2) early forecasting of yields production. Finally, considering the small dimensions of Mediterranean fields, the Copernicus constellation providing high-resolution data with a high revisit time will significantly improve the drought monitoring capability from space and should be considered in future studies.