Prospective analysis of spatial heterogeneity influence on the concordance of remote sensing drought indices: a case of semi-arid agrosystems in Morocco (Moulouya and Tensift watersheds)

Abstract Several drought indices are commonly used to quantify the severity and magnitude of agricultural drought. A comprehensive assessment of temporal and spatial convergences is important. This study evaluates the spatio-temporal concordances of eight drought indices, between October and April, over the period 2003–2021 under the semi-arid conditions of two watersheds in Morocco. These include three monovariate, two bivariate and three composite indices. Taylor diagram, linear regression and Random Forest (RF) clustering were used to evaluate and compare the overall performance of the indices. The results indicate that spatio-temporal variability and spatial heterogeneity of the different agrosystems have little influence on the concordances of composite indices while TCI and VCI have a variable spatio-temporal concordance depending on the spatial extent considered. Nevertheless, depending on the spatial extent considered and the agroclimatic characteristics of each agrosystem, the seasonal correlations between the composite indices and precipitation anomalies, net primary productivity and cereal yield anomalies vary from 0.52 to 0.79, 0.23 to 0.78 and from 0.33 to 0.73, respectively. The correlations between the simple indices and precipitation anomalies, net primary productivity, and cereal yield anomalies range, respectively, from 0.33 to 68 0.13 to 0.72, and from 0 to 0.84. In general, the composite indices showed the best agreement between them and reference indicators but the VCI seems better suited for the assessment of the seasonal profitability of studied agrosystems.


Introduction
Agrosystems of semi-arid environments are spatial units that are distinguished from natural ecosystems by human-induced modifications for exploitation purposes. In recent years, under the combined and cumulative effect of climatic and anthropogenic factors, they have become particularly sensitive to extreme climate variability (Malone et al. 2016;Le Page and Zribi 2019;Fragaszy et al. 2020). In the particular context of Mediterranean arid and semi-arid agroecosystems, drought is the most doubtful climatic hazard (Cook et al. 2016;Ouatiki et al. 2019;Fragaszy et al. 2020;Meliho et al. 2020;Hadri et al. 2021). The unpredictable occurrence in time and space of drought sequences has serious socioeconomic and ecosystem consequences in regions heavily dependent on rainfall (Lui et al. 2018;Sinha et al. 2017;Cartwright et al. 2020;Liu et al. 2020;Liu et al. 2016;Wu et al. 2020;Jiao et al. 2021). Nevertheless, the negative effects of drought sequences can be mitigated by early warning and drought monitoring actions (Enenkel et al. 2015;Masinde, 2015;Didier et al. 2017;Raza et al. 2021;Jiao et al. 2021;Chandrasekara et al. 2021).
However, effective monitoring of drought conditions is often compromised by the occurrence and interactions of the many factors that greatly determine the magnitude and severity of the drought status (Jiao et al. 2021). In the same way, the naturally highly variable landscape heterogeneities, coupled with human-induced changes for food and ecosystem exploitation, challenge the spatial consistency of drought index concordance (Sun et al. 2013;Son et al. 2021;Naumann et al. 2014;Bushra et al. 2019;Bento et al. 2020;Mahyou 2020;Phiri et al. 2020). In our days, due to the diversity of indices and methods used for drought assessment, uncertainties exist as to the multi-scale consistency of biophysical index performance to adequately represent drought conditions taking into account the specific response of each type of agrosystem to drought-induced water stress (Hazaymeh and Hassan 2016;Jiao et al. 2021). These spatial heterogeneities are likely to induce highly differentiated spatial sensitivities Hoylman et al. 2019;Bento et al. 2020;Son et al. 2020;Wu et al. 2020;Bowell et al. 2021). The spatial differentiated sensitivities can soften and/or aggravate the impacts of drought. It should therefore be noted that the quantitative and qualitative assessment of the state of drought can vary significantly according to the specific characteristics of each type of agroecosystem, hence the interest of a prospective analysis of drought sensitivity per homogeneous unit.
In Morocco, in semi-arid agrosystems, such as the Tensift watershed as well as desert agrosystems, such as the eastern highlands, drought episodes have become more intense and more frequent on average 49% of dry years at Guenfouda station and 45% at Oujda station over the period from 1970 to 2016 (El Hafid et al. 2017). These drought episodes have had significant impacts on the productivity of agrosystems (Hanad e et al. 2020;Hadri et al. 2021). In these agroecosystems, the response of natural and anthropogenic units to drought is controlled by many processes. These processes include non-climatic factors namely, topographic factor, soil type, soil water conveyance and retention, interactions with underground water, the existence or lack of irrigation-related inputs Malone et al. 2016;Hoylman et al. 2019). Among these factors, the dominance of poorly evolved soils and the progressive degradation of the vegetation cover is the real factors of sensitivity and fragility of agrosystems in the face of the impacts of drought.
Thus, to cope with the multiple cyclical and structural challenges induced by drought sequences, considerable efforts have been made in terms of monitoring, surveillance and/ or early warning to drought (Fragaszy et al. 2020;Song et al. 2020). The most efficient and spatially adapted indices for monitoring and surveillance of agricultural drought are those based on geospatial and technologically contemporary techniques (Han et al. 2019;Jiao et al. 2019;Liu et al. 2020;Son et al. 2020;Vyas and Bhattacharya 2020;Chang et al. 2021;Tao et al. 2020). The spatial complexity of semi-arid agroecosystems and the need to monitor spatially almost in real-time the resurgence of extreme hazards including drought have led to the testing of several approaches ranging from mono-variable approaches to those integrating a synergy of variables and their spatio-temporal interactions. Recently, Liu et al. (2020) proposed a Drought Gravity Index (DSI) based on total water storage (TWS) and the elimination of the effect of non-climatic factors on the estimation of real drought conditions. Nevertheless, it should be noted that depending on the variables taken into account, the time scale and the approach used, there are very often contradictory conclusions as to the spatial concordance of biophysical indices and their statistical relationships with in situ measurements (Bushra et al. 2019;Li, He, et al. 2020;Bento et al. 2020;Han et al. 2020;Wu et al. 2020;Mahyou 2020;Hadri et al. 2021). Many studies worldwide showed that spatial coherence, as well as correlations between drought indices, vary greatly depending on the specific characteristics off each area (Naumann et al. 2014;Zhuo et al. 2016;Bento et al. 2018;Zhong et al. 2019;Huang et al. 2020;Han et al. 2020;Wu et al. 2020;Cartwright et al. 2020;Son et al. 2020). This suggests that the spatial complexity of semi-arid agroecosystems, including the non-climatic characteristics of the environment, play an important role in drought sensitivity. We can therefore hypothesize that surface and human practices are likely to exacerbate the drought sensitivity of natural and anthropogenic units of semi-arid agroecosystems. The latter coupled with the local variability of climatic and meteorological variables can lead to extreme situations of the state of drought hence the interest in assessing the differentiated sensitivity of drought conditions per homogeneous unit.
In light of previous research, on the same type of climate or sub-climate, which drought index has the best spatio-temporal coherence to describe drought stress? And which indices are the most appropriate and for what types of agrosystems? The answers to all these questions justify the approach developed in this study. This approach seeks to determine the convergences and spatio-temporal divergences of several drought indices in Mediterranean semi-arid agrosystems with particular emphasis on the latest composite models based on artificial intelligence algorithms.
The previous studies have addressed the issue of drought through the development of numerous indices and the assessment of their large-scale spatial coherence (Bijaber et al. 2018;Adnan and Ullah 2020 ;Liu et al. 2020; Vyas and Bhattacharya 2020) but the qualitative assessment of spatial sensitivity at the scale of homogeneous units has not been thoroughly studied in semi-arid agrosystems. The analyses and conclusions drawn remained very general (Bento et al. 2018;Huang et al. 2020;Bento et al. 2020;Cartwright et al. 2020) and do not provide a comprehensive assessment that takes into account the spatial heterogeneity at small temporal and spatial scales that are necessary to improve the accuracy of quantification of drought indices. However, the influence of spatial heterogeneity may be quite considerable on the relevance of drought indices at the scale of homogeneous units that make up an agroecosystem (Son et al. 2021).
In this investigation, spatial heterogeneity includes all surface characteristics that may affect the quantitative and qualitative estimation of the state of drought in the Mediterranean semi-arid agroecosystems. In the current state of knowledge, we still do not know the characteristics of surfaces that influence the sensitivity of agricultural drought indices. It is therefore necessary to study the concordance of drought indices on several spatial units that have different spatial characteristics. This can help in choosing the appropriate index for each type of agrosystem. In this context, this study analyses the influence of spatial heterogeneity on the spatio-temporal concordances of drought indices in several semi-arid agrosystems representative of the Moroccan Mediterranean context.

Study areas
The study area under investigation includes two major hydrological systems in Morocco: the Tensift watershed and the Moulouya watershed ( Figure 1). These two hydrological entities encompass the different Mediterranean sub-climates. They were chosen because they are well suitable for elucidating the question of the spatial concordance of drought indices by remote sensing because of their hydro-climatic and landscape characteristics. They are also quite representative of the arid and semi-arid agroecosystems of the Mediterranean fringe of Morocco. The specific characteristics of each of them are described as follows.

Moulouya watershed
Moulouya watershed is geographically located between latitudes 32 18 0 and 35 8 0 North and longitudes 1 11 0 and 5 37 0 West, on an area of 55,500 km 2 (Figure 1). It is characterized by a marked semi-arid climate of Mediterranean type with low average annual temperatures (12-14 C), low rainfall (<350 mm/year) and strong winds reaching a maximum speed of 50 m/s. The Moulouya river is the longest in northern Africa, with a length of 600 km from the upstream to the outlet, and flows into the Mediterranean Sea. The Moulouya watershed is of major socio-economic importance, particularly for irrigated crops. It is home to Morocco's third-largest citrus region on an area of 21,200 ha in 2018 ).

Tensift watershed
The Tensift watershed is a large continental domain located in west-central Morocco between latitudes 32 10 0 and 30 50 0 North and longitudes 9 25 0 and 7 12 0 West on an area of about 20,450 km 2 or 2.7% of the total area of the country (Meliho et al. 2020). It is one of the most important watersheds in Morocco by its extent, the heterogeneity of its relief and its strong mobilization of water resources. This area is particularly vulnerable to the increased frequency of extreme events (Fniguire et al. 2017). It is characterized by a water deficit due to the combined effect of increasing anthropogenic pressure and episodic droughts (Hajhouji et al. 2018). This watershed is the key sector of the region's economy and employs more than 50% of the working population on an area of about one and a half million hectares, more than 60% of which is reserved for cereal farming. Land use is characterized by poor vegetation that varies according to the altitude and nature of the land and the irrigated area of Haouz is divided into three main parts: Central Haouz, Tessaout upstream and Tessaout downstream (Meliho et al. 2020).

Remote sensing data
The image data used in this study and their characteristics are presented in Table 1. These are data from biophysical variables such as decadal NDVI with a spatial resolution of 500 m (Gidey et al. 2018;Bayissa et al. 2019;Benedict et al. 2021), 8-d land surface temperature (LST) of MODIS MOD11A2 of 1 km spatial resolution (Rhee et al. 2010;Huang et al. 2020;Chang et al. 2021), evapotranspiration abnormality version four (V4) (Bijaber et al. 2018;Sharma and Tare 2018) and the decadal cumulated CHIRPS precipitation with a spatial resolution of 5 km (Guo et al. 2017;Bayissa et al. 2019). These time series, covering the period between 2003 and 2021, were resampled at the spatial resolution of 250 m. In addition to these variables, the annual net primary productivity derived from MODIS data was used in the comparative evaluation the calculated drought indices. We note here that the annual NPP is a product derived from the sum of all products of net photosynthesis over 8 d (MOD17A0H) of each year. Net photosynthesis values are the difference between gross primary productivity and maintenance respiration. The DTM was uploaded to NASA'S earth explorer site data (Table 1) and the choice of agrosystems was based on auxiliary data (shp data, Google earth and the land cover product).  (Lui et al. 2016;Huang et al. 2020;Wei et al. 2020). It is therefore a key variable for the prospective assessment of agricultural drought indices. It is used in the formulation of several indices and is often used for the calibration and validation of drought characterization and prediction models (Zhuo et al. 2016;Huang et al. 2020). In the Moroccan context, a recent study (Bouras et al. 2020) showed the relative importance of the added value of using LDAS data concerning the moisture of the root zone compared to the moisture derived from a single remote-sensing product. In this study, soil moisture was proxied by the difference between day and night LST of the MOD11A2 product (Bijaber et al. 2018). The idea stems from the fact that the morning evolution of LST is strongly related to soil moisture conditions (Hain et al. 2011). The principle put forward is that a soil covered with water-stressed vegetation warms up rapidly than a wet surface (Anderson et al. 2007). This approach was chosen because it offers the best compromise to reconstruct a time series of seasonal soil moisture data for the different considered homogeneous units over the studied period (2003-2021).

Land use data
In this study, the identification and extraction of agroecosystems were facilitated by the use of land cover and land use data. To this end, we merged several land cover units from the ESA-CCI land cover map of Morocco to reconstitute the relatively homogeneous units ( Figure 2). This product has 22 land cover class including grasslands, irrigated land, cropland, savannah, deciduous forests and mixed forests ). The highresolution Google Earth visualization made it possible to extract the limits of the homogeneous units from which the overall evaluation of the spatio-temporal concordance of the various drought indices considered was conducted.

Reference data
The reference data series were used mainly for validation and graphical representation of the results. They include precipitation data from weather stations, data on the primary productivity of agroecosystems, cereal yield data for the period 2003-2020, reference maps. The World Bank's database of major cereal crops (Barley, durum wheat and soft wheat) was used in this study and rainfall data (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021) are the monthly cumulations of three stations per watershed obtained via the https://www.infoclimat.fr/climatologie/ site. We considered for each watershed the average precipitations for three weather stations taking into account the spatial distribution of the stations. For this purpose, the precipitation data series from the Safi, Marrakech and Essaouria stations were used for the Tensift watershed and the Taza, Oujda and Nador stations were considered for the Moulouya watershed.

Methodology
To study the multi-scale spatial coherence of drought indices by remote sensing, the approach used in this study combines several methods and techniques of treatments specific to each dataset. First, several remote sensing drought indices were constructed from several sources of the multiplier image data (Table 1). After this step, high-resolution visualization (Google Earth) and axillary data were used to delimit and extract the homogeneous entities (agrosystems) on which the concordance of the indices was evaluated on several time scales and the last part of the approach used is a validation of the different results obtained. Figure 3 presents the main steps in the methodology used in this study.

Multi-scale construction of biophysical indices
Three types of indices were calculated on a monthly and seasonal scale to study the spatio-temporal concordances of drought indices over the last 19 years. The first category of indices is presented by single monovarietal indices (TCI, VCI and PCI), the second category includes bivariate indices (VHI and NVSWI) and the last category gather three relatively recent composite models (CDMI, SDCI and CDMIa_RF). All these indices are constructed based on remotely detected biophysical parameters from climatological averages (monthly and seasonal) calculated over the period 2003-2021.

Monovariate indices
The three monovariable indices used in this study are the vegetation condition index (VCI), the temperature condition index (TCI) and the precipitation condition index (PCI). The VCI is calculated from the time series of vegetation index data from Equation (1) (Kogan 1990). Similarly, the PCI is calculated based on historical precipitation anomalies from Equation (2). The VCI monitors the moisture content of vegetation, while the TCI assesses heat stress in vegetation. The combined use of these two indices is increasingly recommended. The TCI index Equation (3) is based on the brightness temperature from the maximum and minimum conditions of the time series considered (Kogan 1995).
NDVI (i), NDVImax and NDVI min, are the pixel values of NDVI and maximum, minimum of it, respectively, in the time series considered.
P(i), Pmax and Pmin, are the pixel values of monthly and seasonal averages of CHIRPS data precipitation and maximum, minimum of it, respectively, in the time series considered.
where LST (i) is the monthly and seasonal averages of 8 d MODIS LST and LSTmax and LSTmin are, respectively, the max and min pixel values of the time series considered.

Bivariate indices
The bivariate indices used are the Vegetation Health Index (VHI) and the Vegetation Water Supply Index (VSWI). They were chosen because they would be better suited to agricultural drought in semi-arid areas where water is the main factor influencing vegetation growth (Karnieli et al. 2010;Rhee et al. 2010;Zhong et al. 2019). The VHI is an index that combines the information provided by both the VCI and the TCI. Similarly, the Normalized Water Stress Index (NVSWI) is considered a comprehensive drought monitoring method based on LST and NDVI (Equation (5)).
In this equation, VHI of month (i) or year (i) is equal to the sum of VCI (i) Ã 0.5 þ TCI (i) Ã 0.5 at the time scale considered.
NDVI and LST of the same month or year, VSWImin and max are the maximum and minimum values of the time series.

Composites index
Since the onset of drought is often linked to interactions of the many factors, three composite models were calculated to compare their spatio-temporal concordance with that of the mono and bivariate indices. The composite models under investigation are the CDMI, CDMIa_RF and SDCI. These composite models have been constructed (equations) from an objective aggregation of monovariable indices derived from monthly and seasonal anomalies of various factors involved in the development, evolution and/or worsening of the drought phenomenon. The CDMI model was developed in China by Han et al. (2019) by incorporating using a machine learning the anomalies of four variables (P, LST, NDVI and ET). The CDMIa_RF is also an additive model inspired by the CDMI model adding the anomaly component of soil moisture developed in the context of Tensift (Morocco). The SDCI is the result of the work of Rhee et al. (2010) in the arid/semi-arid regions of the Great Plains and the western United States. It combines three remote sensing variables (P, LST and NDVI) from a correlative weighting.
The classification used for the cartographic comparison of different indices is provided in the following Table 2.
The weighting of the composite models was inspired and comparable to the official composite model (drought composite index) adapted throughout the agricultural zone of Morocco (Bijaber et al. 2018).

Identification and extraction of homogeneous units
To establish an in-depth analysis of the influence of spatial heterogeneity on the concordance of remote sensing drought indices at the monthly and seasonal scales, the delimitation of natural and anthropogenic units of land cover and land use is important to relativize the applicability of each index and/or indicator. To achieve this, homogeneous entities were extracted from Morocco's land use mapping and agro-climatic zoning. In this study, a homogeneous unit corresponds to any type of natural and/or anthropogenic ecosystem composed of elements mainly of the same nature, characterized by a certain spatial continuity and which differ from the limiting surrounding context. Thus, given the landscape heterogeneity of our study area, we considered homogeneous units that best represent agro-climatic contrasts. For this purpose, the analysis was conducted on six spatial domains. There are a favourable, unfavourable, irrigated and mountainous agrosystems, to which are added analyses at the scale of the Tensift and Moulouya watersheds (Figure 1).

Construction of reference indicators
In this study, anomalies of cereal yields, precipitations and net primary productivity were used to assess the spatio-temporal concordances of the different indices. Crop yield data are widely used as an important factor in monitoring particularly agricultural drought (Zhuo et al. 2016;Benabdelouahab et al. 2019;Bouras et al. 2020;Wei et al. 2020). However, due to many factors that affect agricultural productivity, agricultural yield data need to be normalized and reconstructed into agricultural yield indicators so that they can serve as benchmarks of the impact of drought on agricultural productivity (Sun et al. 2013;Huang et al. 2020;Bushra et al. 2019). For this purpose, the three series of reference data used for the validation of the results have been normalized with the maximum and minimum values to approximate the order of scale of the different indices.

Evaluation and validation metrics
First, Pearson correlation and linear regression were used to compare the overall chronological coherence of biophysical index models over different seasons against the reference dataset Guo et al. 2019;Vyas and Bhattacharya 2020;Benedict et al. 2021). Second, analyses of standardized Taylor diagrams have been used to complement bias analyses because they allow for graphical visualization of the spatial Table 2. Classification of the different drought indices used (Kogan 1995;Rhee et al. 2010;Han et al. 2019 variability of indices relative to baseline data (Taylor 2001). Because of its synthetic representation in two-dimensional space of several statistical measurements simultaneously, the Taylor diagram is often preferred to compare the series of calculated or simulated data of several models concerning real observations (Samadianfard et al. 2018;Sharma and Tare 2018;Bisht et al. 2019;Malik et al. 2019Malik et al. , 2020He et al. 2021). The representations of the Taylor diagrams user were produced in the Rstudio environment with R code available online. For this purpose, the monthly and seasonal average time series of the values of each index on several homogeneous entities representative of the agrosystems of the area of interest was used in comparison with the reference indicators. The in-situ precipitation reference data series and historical anomalies in cereal yields were used to evaluate index concordances by agrosystem type.

Evaluation using a machine learning algorithm (clustering)
In addition to statistical validation using Taylor's diagram, and Pearson correlation, the Random Forest (RF) algorithm was used to evaluate index concordances by agrosystem type. For this purpose, the in-situ precipitation reference data series and historical anomalies in cereal yields were used to cluster the relevance of each index in six selected agrosystems. The aim is to make an analogue comparison with the results of the statistical indicators (Taylor diagram and correlation). In this study, the default parameter (ntree ¼ 500, mtry ¼ 1, et.seed (123)) of the RF algorithm in the R 4.1.0 environment was considered because they gave the most consistent results in the statistical analyses of linear regressions and the Taylor diagram.  (Figure 4), it can easily be observed that the minor normal distribution was not observed on the seasonality of the VCI index in the context of Tensift while it is well observed in the context of the Moulouya watershed. In the same direction and for the same indicator, the curve of the major normal distribution (blue colour, VCI index) shows the most marked phase shift of all index distributions for the four types of agrosystems. From this result, it appears that the three categories of indices used in this study have a very similar seasonal variance in the four agrosystems. As a result, only the statistical and cartographic analyses established in the following sections have made it possible to objectively assess the spatio-temporal concordances.

Identification of the most relevant indicators
In accordance to identify the most relevant indicators, we compare different matrices of Pearson correlations for different contexts and time scales ( Figure 5). For this purpose, the observed precipitation anomalies (reference indicator) of the three stations per watershed were used. The analysis of the result obtained shows that in the six spatial contexts considered, the correlations of the composite indices (CDMIa_RF, SDCI and CDMI) exceed those of the mono and bivariable indices (VHI and NVSWI). The inter-index composite correlative values reach the extreme value of 0.99 and their correlation with precipitation anomalies varies between 0.59 and 0.79. The monovariable indices are characterized by the lowest and highly variation correlations depending on the temporal, spatial and/or agro-climatic scale considered. Their correlations compared to monthly precipitation anomalies (reference data) converge towards negative correlations in the context of semi-irrigated agrosystems (À0.0037) as indicated by the black colour frame in Figure 5 for the case of the VCI index. This indicates that the vigour and development of vegetation in irrigated areas are very little dependent on monthly precipitation variability. On the contrary, in the unfavourable agrosystem, the VCI achieves significant correlations with precipitation anomalies. The indicator seems better suited to reproducing drought conditions in the arid zone where the development of vegetation cover is strongly conditioned by the monthly cumulative rainfall.  unfavourable agrosystem for the same index (Table 3). This clearly indicates the influence of spatial scale, climatic and surface characteristics (land cover) in the representativeness of correlations of biophysical indices. Indeed, the Tensift watershed consisting mainly of agricultural land has less noise related to spatial heterogeneities unlike the Moulouya watershed characterized by a dominance of more arid climatic conditions and an essentially steppe land use that would add additional noise in the seasonal average values of the indices.

Cartographic convergences of composite models in two watersheds
The prospective analysis of historical cartographic similarities of the state of agropastoral drought at the watershed scale (Tensift & Moulouya) using the CDMIa_RF model aims to determine the ability of composite indices to adequately reproduce the parameters of agricultural drought regardless of the agro-climatic context. To do this, the conditions of the drought state of the Tensift watershed with essentially semi-arid climatic characteristics were compared to the historical dynamics of the state of the drought of the Moulouya watershed, characterized by a more heterogeneous climatic context with a clear dominance of arid conditions in its pastoral zone. Thus, from the retrospective and comparative analysis of the spatio-temporal dynamics of drought conditions presented in Figure 6, it is worth highlighting the specific ability of the composite model CDMIa_RF to detect with absolute consistency the extreme conditions (wet and dry) in the two hydro-agro-climatic units. The irregular alternation of dry, wet and/or normal years was observed quite identically in both watersheds. In Tensift, as in the Moulouya, over the last 18 years of the time series used, at least 5 years have been characterized by a dominance of extreme drought conditions over almost half of the areas of the two watersheds. In Tensift, it was the years 2005,2008,2012,2016 and 2020 that were the driest and in Moulouya, it was the years 2005, 2014, 2016, 2019 and 2020 that were the driest. Normal and wet conditions were observed similarly in very variable proportions depending on the year. Exceptionally wet years were also detected identically in both watersheds. However, it should be noted that the spatial distribution and the occurrence of pockets of severe to moderate droughts in the two watersheds do not follow any spatial and/or temporal dependence. The thresholds of the intensity of moderate drought were observed with great spatial variability in Tensift as well as in Moulouya watersheds.
This comparative analysis of spatial similarities is very well illustrated by the result of the monthly index dynamics in the two watersheds (Figure 7). Thus, we deduce from this result that the overall monthly variability of the indices in the two watersheds is very similar to a few nuances. All maximum and minimum peaks from January 2003 to March 2021 were highlighted with a remarkable agreement in both hydro-climatic systems. It is also well observed in both watersheds that extreme wet conditions rarely reach 4 months out of the 7 months of the agricultural season. This remains a common feature in all seasons of the period (2003-2021) even for the cases of the wettest years of the period under review which record only 4 months where very favourable wet conditions were observed. However, the length of extreme drought sequences is less important, rarely reaching the three consecutive months. Moreover, in both contexts, the major potential risk of agricultural drought lies in the very random occurrence of the monthly sequences of extreme droughts. A normal or surplus year does not reflect the absence of pockets of prolonged monthly drought sequences. The negative impacts of drought are, therefore, a function of the seasonal distribution and length of these dry sequences. These two characteristics are the most difficult to assess and predict. They are one of the typical characteristics of each agricultural year and the very random nature of their intensity further aggravates the risk and impacts associated with this hazard.

Comparative mapping of the longest dry and wet sequences
Comparative mapping of short dry and wet sequences allows for a more detailed examination of the spatial concordances of indices. To do this, we reproduced the monthly conditions of the agricultural drought of the last agricultural season (2020-2021) using six  indices (two composite indices, two monovariable indices and two bivariate indices) on an agrosystem classified as favourable. The 2020-2021 agricultural season that recorded a surplus crop year was chosen for this purpose as a reference year to highlight the spatial concordances and the relevance of the indices at the scale of a favourable agrosystem to identify the divergences between and intra indices. The analysis of the result shown in Figure 8 makes it possible to note the merit of the very distinct spatial convergences of the SDCI composite index and the CDMIa_RF model, which result in a very coherent mapping of the monthly drought conditions of the 2020-2021 agricultural season produced by the two composite indices. Indeed, the signals of the unfavourable conditions that characterized the beginning of the agricultural season, particularly the month of November, which recorded a dominance of exceptional drought conditions, were well captured by the SDCI index and the CDMIa_RF.
On the other hand, for the other indices, although the unfavourable conditions at the beginning of the season are captured, the analysis of this result reveals a weak spatio-temporal convergence of the monovarietal indices (TCI and VCI). The weak convergences of the VCI and TCI indices are observed for certain months of the agricultural season, in particular, the month of January when the very favourable conditions were highlighted by the TCI. Conversely, the VCI showed a favourable dominance situation. Overall, except for the composite indices (SDCI and CDMIa_RF) which remained in line with the real situation that characterized the 2020-2021 agricultural season, the other mono and bivariable indices tend to overestimate drought conditions. This is well highlighted by the comparative mapping of the last 4 months of the 2020-2021 agricultural season which were considered favourable, according to the indicators of the Moroccan Ministry of Agriculture, and illustrated by the average monthly cumulative rainfall.

Random Forest clustering comparative of indices by type of agrosystem
The clustering of the relative importance of the indices according to the reference indicators and by types of agrosystems was carried out by machine learning with the RF algorithm. The outputs of the model (RF) clearly illustrate the results of the statistical and cartographic metrics established above. In this exercise, the default model parameters were used because the model outputs showed a satisfied consistency with the above analyses performed. The results obtained as shown in Figure 9 above reflect the rather distinctive merit of the relevance of composite indices in most agrosystems and agro-climatic context considered. In terms of consistency with seasonal anomalies in cereal yields, the enhanced composite model performs much better than the other two composite models (SDCI and CDMI). However, exceptionally at the scale of an irrigated agrosystem, the VCI slightly outperforms the CDMIa but its concordance with yield anomalies in other types of agrosystems remains very low, which considerably reduces its relevance for the need for large-scale heterogeneous analysis. Similarly, with reference to monthly precipitation anomalies, VCI and PCI outperform composite indices in the unfavourable agrozone, but their concordances with monthly precipitation anomalies fall considerably in other types of agrosystems. Overall, in all agrosystems and on a monthly and seasonal time scale, the classification of the output of the RF model distinguishes with a significant Figure 9. Random Forest clustering of the relative importance of indices with reference to monthly precipitation anomalies (red colour) and seasonal anomalies in cereal yield production (blue colour).
difference the relative importance of the indices by type of agrosystems and by reference indicator.

Results validation
The comparative validation of spatio-temporal concordance of the studied drought indices was carried out by comparing the average values of these indices with the anomalies of three reference indicators: observed rainfall, cereal yields and net primary productivity. For this purpose, comparative representations of the Taylor diagram of each agrosystem for all indices were preferred over baseline indicator anomalies. To do this, the result of the first validation comparing the monthly and seasonal variability of precipitation anomalies with the variability of indices by types of agrosystems and agro-climatic zones is provided in Figure 10. This result highlights indices whose spatio-temporal variance is closest to the monthly and seasonal variance of precipitation anomalies. This also made it possible to identify the indices that have the best spatial convergences with each other and with reference data. This result shows a better spatio-temporal concordance of composite models (SDCI, CDMI and CDMIa_RF) in the six agrosystems considered in comparison with the observed precipitation anomalies.
On the other hand, it is well observed (black and blue dots) representing respectively the VCI and the TCI, concordances with the anomalies of precipitation very variable depending on the context. The VCI index is well correlated with precipitation anomalies in the unfavourable zone unlike the scale of a semi-irrigated agrosystem, the VCI has monthly correlations with almost zero precipitation anomalies. In the same vein, the best concordances of the TCI index with the other indices and reference data were observed in the mountainous zone and the favourable agrosystem. The VCI and TCI are very often used as indicators for monitoring agricultural drought, but this study found some variable inconsistencies with reference data depending on the temporal and spatial scale considered. On the contrary, changes in temporal or spatial scale have very little effect on the concordance of composite models.
To further assess the comparative assessment of spatio-temporal concordances and divergences at the level of the six spatial domains with different agro-systemic characteristics, the seasonal average series of indices were compared with anomalies in cereal yields and net primary productivity. Figure 11 summarizes the statistical relationships of indices compared to benchmark indicators in the agrosystems considered. Unlike the comparison with rainfall, the comparison made concerning anomalies in cereal production reveals relatively high concordance of the VCI index in the Tensift watershed and the irrigated agrosystem of the Moulouya. In both contexts, the correlations of the VCI index with anomalies in cereal yields surpass those of the other indices. In the same sense, the VCI is distinguished by a variance quite identical to the variance of cereal yield anomalies in the irrigated agrosystem.
However, it should be noticed that the rather variable nature of the performance of the VCI index depending on the context or type of agrosystem significantly limits the relevance of this index for the assessment of the risk of agricultural drought on a large heterogeneous scale. Thus, the global analysis shows the merit allocated the more stable concordance of composite models, especially the CDMIa_RF. The latter is distinguished by significantly high correlations in all agrosystems considered and a variance very close to the variance of yield anomalies in the case of irrigated agrosystem. Moreover, except for the TCI index, the comparison of the indices to the anomalies of the NPP shows very little distinct concordance between the indices. However, the VCI slightly outperforms other indices in the favourable zone and composite indices outperform it slightly in other types of agrosystems.

Discussions
The choice of relevant indicators for the qualitative and quantitative assessment of real drought conditions is a major imperative in terms of early warning and mitigation of the impacts of climate hazards on a large spatial scale. This study examines the spatio-temporal concordances of seven agricultural drought indices by spatial remote sensing at the scale of several spatial domains. The spatio-temporal concordance of these composite indices was evaluated against several benchmark indicators (observed precipitation anomalies, yields and NPP) and other very popular indices (VCI, TCI, VHI and NVSWI). Thus, of all the statistical analyses carried out (Taylor diagram, correlation) the results obtained indicate differentiated applicability between simple, bi-variable and composite indices.
The concordance of composite indices (CDMIa_RF, SDCI and CDMI) far exceeds that of bivariate and monovarietal indices in most of the agrosystems considered, while monovariable indices (VCI and TCI) are particularly affected by the change in temporal and/or spatial scale. This result corroborates with the conclusions of the study of Li, He, et al. (2020) and  carried out in mainland China and which revealed that composite indices (VHI, SDI and SDCI) work better than simple drought indices and that all indices have been characterized by different applicability across mainland China. Another China-wide study (Wei et al. 2020) demonstrated that the correlations of composite indices with the in-situ drought index are higher than those of simple indices across all ecosystems considered. In the same vein, many other studies have highlighted the highly variable performance of the VCI depending on the climate context and scale. In the CONUS area, the study by Zhang et al. (2017) also revealed the highly variable characteristics of the consistency of the VCI index, which was more reliable in semi-arid areas of low latitudes than in the high-density wetland of vegetation cover. Indeed, under significant irrigation conditions, the VCI has significant inconsistencies compared to other indicators and its inconsistencies have been observed differently in all agrosystems and at the scale of the two watersheds studied in this article. Other factors also such as the impact of coarse resolution and the quality of image data can particularly affect the VCI index which is a direct derivative of NDVI. Other authors (Parviz 2016) suggest that the effect of topography limits the short-term application of vegetation-based indices. In this investigation, the much better concordance of composite indices compared to other indices and benchmarks was observed in all agro-systemic contexts.
Indeed, when the analyses are carried out on a monthly, there are some inconsistencies of the monovariable indices in some places and months of the agricultural year. These inconsistencies could stem from the influence of several factors, namely the stage of the phenolic cycle of crops, the production and agricultural systems used and the occurrence of extreme weather situations. For example, extreme heatwaves would explain some abnormal peaks in the TCI index and irrigation and intensive fertilization that could explain some peaks in the VCI index. Indeed, the first advantage of composite indices over variable united indices comes from the fact that the incorporation of several variables considerably reduces the influence of other factors on the biophysical variable. Therefore, simple indices can be reliable but with a high risk of overestimating or underestimating drought conditions in the event of extreme situations carried by the single variable used in the construction of a monovarietal index. In the same vein, as illustrated by the comparison of the indices for the 2021 crop year considered satisfactory, the VHI and NVSWI indices were characterized by an amplifying effect of the real state of agricultural drought.
However, regarding cereal yield anomalies, seasonal scale analyses (7 months), the VCI was very competitive but with limited performance in two spatial areas (irrigated agrosystems and at the Tensift watershed scale) where the correlations of the VCI index with yield anomalies slightly outperform those of the composite indices. This is consistent with the result of a recent study (Bouras et al. 2020) carried out at the level of the agricultural provinces of Morocco which revealed through analyses of the correlations shifted between the VCI index and the cereal yield correlations of 0.94 in the province of Khemisset. These two conclusions converge towards a separate relevance of the VCI index for forecasting cereal yields. Nevertheless, the specific capabilities of composite models (SDCI, CDMI and CDMIa_RF) to capture agricultural drought signals and yield forecasting have not been investigated compared to the indices used by Bouras et al. (2020) at the scale of Morocco's main agricultural provinces. While composite models have the advantage of better spatio-temporal concordances regardless of the agro-climatic context and the time scale considered. To this end, they offer very good prospects for seasonal monitoring of the overall profitability of agrosystems and forecasts of yields at the national scale, which are very heterogeneous. The major disadvantage of monovariable indices, in particular the VCI index, lies in their highly variable performance depending on the time scale considered and the agro-climatic context. Seasonal scale analyses have shown that VCI is strongly correlated with rainfall in the pre-Saharan arid zone than in semi-arid agrosystems. At the irrigated perimeter scale, the monthly average values of VCI showed correlations very close to zero with the monthly anomalies of the observed precipitations. In contrast, composite indices compared to other indices and benchmark indicators showed more general applicability and low variability regardless of agro-climatic context and time scale. The low variability in composite index performance has also been observed in the continental United States (Zhang et al. 2017) and recently in Central Asia (Chang et al. 2021). Composite indices, including SDCI, have distinguished themselves in other studies that are more effective in short-term operational forecasts of the state of drought (Park et al. 2020). In other words, and in a global way, it should be emphasized that the spatiotemporal concordances of the indices are very variable according to the temporal and/or spatial scale considered. This characteristic is particularly pronounced in monovarietal indices that are distinguished by extreme anomalies that limit their effectiveness.

Conclusion
This analysis shows that spatial heterogeneities as well as spatial scale change have little influence on the concordance of composite indices in the context of Morocco's semi-arid agrosystems. The spatio-temporal coherences of composite indices (SDCI, CDMI and CDMIa_RF) with reference indicators (agricultural yield anomalies, NPP, rainfall) are weakly affected by changes in the spatial, climatic, irrigated or non-irrigated context. On the contrary, the spatio-temporal concordance of the monovariable indices (CI, TCI) with the reference data is very variable according to the agroclimatic characteristics of each agrosystem. The seasonal VCI was strongly consistent with cereal yield anomalies in the Tensift watershed with a correlation of 0.84 slightly higher than the composite models of 0.73. But in other agrosystems, the concordance between the VCI and yield anomalies is much lower than that observed in the Tensift watershed. However, it should be noted that the comparative analyses carried out in this study are based on a series of data on cereal yield at the national level. Analysis based on local or parcel scale yield data is needed to refine the conclusions of this study. Similarly, the factors underlying the variability of the VCI and TCI concordances have not been studied in this study, therefore, future studies will need to focus on identifying and quantifying the sensitivity factors of each index using Sentinel 2 and Landsat images.