Mapping and Characterization of Phenological Changes over Various Farming Systems in an Arid and Semi-Arid Region Using Multitemporal Moderate Spatial Resolution Data

: Changing land use patterns is of great importance in environmental studies and critical for land use management decision making over farming systems in arid and semi-arid regions. Unfortunately, ground data scarcity or inadequacy in many regions can cause large uncertainties in the characterization of phenological changes in arid and semi-arid regions, which can hamper tailored decision making towards best agricultural management practices. Alternatively, state-of-the-art methods for phenological metrics’ extraction and long time-series analysis techniques of multispectral remote sensing imagery provide a viable solution. In this context, this study aims to characterize the changes over farming systems through trend analysis. To this end, four farming systems (fallow, rainfed, irrigated annual, and irrigated perennial) in arid areas of Morocco were studied based on four phenological metrics (PhM) (i.e., great integral, start, end, and length of the season). These were derived from large Moderate resolution Imaging Spectroradiometer (MODIS) normalized difference vegetation index (NDVI) time-series using both a machine learning algorithm and a pixel-based change analysis method. Results showed that during the last twenty-year period (i.e., 2000–2019), a signiﬁcant dynamism of the plant cover was linked to the behavior of farmers who tend to cultivate intensively and to invest in high-income crops. More speciﬁcally, a relevant variability in fallow and rainfed areas, closely linked to the weather conditions, was found. In addition, signiﬁcant lag trends of the start ( − 6 days) and end (+3 days) were found, which indicate that the length of the season was related to the spatiotemporal variability of rainfall. This study has also highlighted the potential of multitemporal moderate spatial resolution data to accurately monitor agriculture and better manage land resources. In the meantime, for operationally implementing the use of such work in the ﬁeld, we believe that it is essential consider the perceptions, opinions, and mutual beneﬁts of farmers and stakeholders to improve strategies and synergies whilst ensuring food, welfare, and sustainability.


Introduction
Accurate monitoring of land use changes in arid landscapes and understanding of these changes drivers are crucial in arid environment research, especially since changes in farming systems affect directly socioeconomic as well as environmental sectors. Improving global food security will need a good understanding of the behavior of farming systems Remote Sens. 2021, 13, 578 3 of 21 analysis. Investigating the changes in farming systems at a large scale in the OER basin will have a far-reaching impact on agricultural sector development and how the system is adapting or not to climate changes. In this context, the present research sought to (i) use phenological metrics derived from twenty years of NDVI MODIS datasets (i.e., [2000][2001][2002][2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013][2014][2015][2016][2017][2018][2019] to map and monitor changes in selected farming systems over a large arid-to-semi-arid region in Morocco (i.e., OER basin); (ii) investigate trends in selected farming systems at a large scale in the study area to evaluate how the systems are adapting or not to changes; and (iii) provide information on farming systems' changes for stakeholders to adopt more accurate and efficient strategies.

Study Area
This study was conducted within the Oum Er-Rbia (OER) basin in west-central Morocco, between 31 • 15 -33 • 22 N and 5 • 08 -8 • 23 W (Figure 1). The OER basin covers an area of 38,000 km 2 , while its administrative area is around 50,000 km 2 . The OER basin is made up of five geographical units from the Atlas Mountains, the foothill areas, the plain, the phosphate plateau, and the coastal area [4]. A combination of flat and mountain terrains generally characterize the topography of the basin. Elevation ranges between 100 (e.g., in the western and coastal zones) and 3890 m (e.g., in the eastern zone) above sea level [42] ( Figure 1). The OER River sources are located in the mountainous upstream zones, and the river covers a distance of 550 km, overpassing the Tadla irrigated perimeter (TIP), the coastal areas, and the northern zone of the Doukkala irrigated perimeter (DIP), and flows into the Atlantic Ocean at Azemmour city [4]. The climate is variable from humid in coastal and mountainous zones to semi-arid in the plains, with cold winters and dry summers [42]. Annual rainfall average varies from 230 to 1000 mm in the plains and the Atlas Mountains, respectively [43]. The agricultural season generally occurs between September and June, while the most important amount of rainfall is received between October and May (i.e., 70% to 80% of the annual rainfall). The annual temperature varies between a maximum of 46 • C in August and a minimum of −3 • C in January. The investigated areas are primarily agricultural, with irrigated crops (i.e., cereals, sugar beet, and alfalfa) and rainfed areas (wheat). The study area is also dominated by tree cultivations, especially pomegranate and citrus trees, which present a permanent activity during the season ( Figure 2).
In this study, four farming system types ( The phenology of the IPC farming system is characterized by permanent vegetative activity during the cropping season and high biomass production ( Figure 2). IAC represents high inter-annual variability with a high amplitude and rapid growth and senescence moments. The IAC farming system is described with high vegetative activity with a maximum NDVI generally in March, low dependency on climate conditions and soil type, and a maintained growth cycle from seed to harvest time. The RA and FA farming systems are typical on semi-arid lands, where the start of vegetative activity depends on the climate conditions, especially the first rainfall. Figure 2 gives deep informative illustrations of the four farming systems investigated in this study.

Multitemporal MODIS and CHIRPS Data Acquisition and Derived NDVI Processing
In this study, we used the MOD13Q1 NDVI product. A dataset of 437 images covering the OER basin was downloaded for the period between 2000 and 2019 via the U.S. Geological Survey (USGS) Land Processes Distributed Active Archive Center (LP-DAAC). The MOD13Q1 NDVI data are within a spatial resolution of 250 m and a temporal frequency of 16 days in sinusoidal projection [4]. Time series data were re-projected to UTM projection to be used further in this study. Rainfall data from Climate Hazards group Infrared Precipitation with Stations (CHIRPS) [44], available from 1981, were summed over the OER Remote Sens. 2021, 13, 578 4 of 21 basin for each year within the MOD13Q1 product period. The CHIRPS data are delivered at two spatial resolution levels, 0.05 and 0.25 • . In our study, CHIRPS data converted to a 16-day period with a spatial resolution of 0.05 • were used to extract rainfall data and plotted against NDVI profile.
We performed NDVI time-series filtering and PhM extraction with TIMESAT software [45], a commonly used package used for (1) data filtering and smoothing and (2) phenological metrics extraction ( Figure 3). In order to obtain high-quality time series and to fill the gaps in data, TIMESAT implements three different smoothing methods, which are Savitzky-Golay (SG) [46], asymmetric Gaussian (AG) [27,47], and double logistic (DL) [48]. The AG filtering method is less sensitive to noise than the other approaches [45,49]. The study performed by Fu et al. [50] shows the robustness of this algorithm in smoothing time series and preserving the characteristics of NDVI profile, contrary to the DL and SG methods. The AG method can generate a smoothed NDVI profile while describing minor changes in the NDVI sequence data [51]. Therefore, the AG method was chosen as the time series reconstruction method to be used in this study. Table 1 shows the setting used for NDVI time series processing in TIMESAT. More details on the filtering approach used and phenological metrics extraction for this study can be found in Lebrini et al. [4].    In this study, four farming system types ( The phenology of the IPC farming system is characterized by permanent vegetative activity during the cropping season and high biomass production ( Figure 2). IAC represents high inter-annual variability with a high amplitude and rapid growth and senescence moments. The IAC farming system is described with high vegetative activity with a maximum NDVI generally in March, low dependency on climate conditions and soil type, and a maintained growth cycle from seed to harvest time. The RA and FA farming systems are typical on semi-arid lands, where the start of vegetative activity depends on the climate conditions, especially the first rainfall. Figure 2 gives deep informative illustrations of the four farming systems investigated in this study.

Multitemporal MODIS and CHIRPS Data Acquisition and Derived NDVI Processing
In this study, we used the MOD13Q1 NDVI product. A dataset of 437 images covering the OER basin was downloaded for the period between 2000 and 2019 via the U.S. Geological Survey (USGS) Land Processes Distributed Active Archive Center (LP-DAAC). The MOD13Q1 NDVI data are within a spatial resolution of 250 m and a temporal frequency of 16 days in sinusoidal projection [4]. Time series data were re-projected to UTM projection to be used further in this study. Rainfall data from Climate Hazards group Infrared Precipitation with Stations (CHIRPS) [44], available from 1981, were summed over the OER basin for each year within the MOD13Q1 product period. The CHIRPS data are delivered at two spatial resolution levels, 0.05 and 0.25°. In our study, CHIRPS data converted to a 16-day period with a spatial resolution of 0.05° were used to extract rainfall data and plotted against NDVI profile. We performed NDVI time-series filtering and PhM extraction with TIMESAT software [45], a commonly used package used for (1) data filtering and smoothing and (2) phenological metrics extraction ( Figure 3). In order to obtain high-quality time series and to fill the gaps in data, TIMESAT implements three different smoothing methods, which are Savitzky-Golay (SG) [46], asymmetric Gaussian (AG) [27,47], and double logistic (DL) [48]. The AG filtering method is less sensitive to noise than the other approaches [45,49]. The study performed by Fu et al. [50] shows the robustness of this algorithm in smoothing time series and preserving the characteristics of NDVI profile, contrary to the DL and SG methods. The AG method can generate a smoothed NDVI profile while describing minor changes in the NDVI sequence data [51]. Therefore, the AG method was chosen as the time series reconstruction method to be used in this study. Table 1 shows the setting used for NDVI time series processing in TIMESAT. More details on the filtering approach used and phenological metrics extraction for this study can be found in Lebrini et al. [4].    The second step was PhM extraction; such measurements are usually calculated using a common method based on value thresholds of the seasonal vegetation amplitude, assuming that a particular phenomenon has started when the NDVI values surpass a given threshold ( Figure 3). Our study was, thus, carried out by setting the proportion of the seasonal amplitude to 10% measured from the left and right minimum values, respectively [45]. In general, four PhM were extracted using the TIMESAT program for trend analysis. Definitions of the PhM used in this study are explicitly defined in Table 2. Table 2. Definition of the computed phenological metrics (PhM) [52][53][54]. To characterize the main cropping systems, the ground data used in this research were collected through ground truthing exercises over farming systems during the 2018-2019 cropping season. The collected data were reported using a GPS receiver, with a positional error of less than 2 m to generate data for classification training and accuracy assessments for the 2018/2019 dataset. The reference points were collected in a way so as to represent the full variety of farming system elements in the study area. Accuracy assessment of other classification results was performed using similar or near-similar ancillary data from MODIS and Google Earth images and high-resolution aerial photographs. A land cover map originated from Glob Cover was used to mask the agricultural zones over the study area [4,55]. A summary of the ground data is given in Table 3. To classify farming systems over the OER basin, the supervised Random Forest (RF) classifier was used based on the CARET package within R [56]. RF is a non-parametric machine learning classifier that combines a random selection of training subsets of data with an ensemble of trees. Recent studies reveled the effectiveness of the RF classifier in the remote sensing field, including land phenology mapping [39,[57][58][59]. For measuring the accuracy of the classification results using the RF classifier, ground truth data were split into two sets of training (80%) and testing (20%) samples using a spatial cross-validation approach with 5 folds from the CAST package in R [60]. This spatial cross-validation helped to make sure that the ground truth sample of the same field will be either in the test or in the training data in order to avoid over-fitting. The accuracy assessment was performed by calculation of many accuracy metrics for the classifications results, including overall accuracy (OA), Kappa coefficient, producer's accuracy (PA), user's accuracy (UA), and F1-score. The classification of 2000/2001 was generated using the "predict" function from the CARET package in R [56]. In order to increase the reliability and validity of our RF classification model, and as an additional check of the resultant information of study area-specific classification accuracy, a second accuracy assessment was performed for the 2000/2001 classification map using the same methods as in the accuracy assessment for the 2018/2019 classification map. In order to update testing polygons based on the 2000/2001 situation, we used the testing polygons from the 2018/2019 data and plotted them against the smoothed NDVI profiles from MODIS data and Google Earth images. From the MODIS time series, we investigated the correspondence between the NDVI profile of each testing polygon and the farming system. From imagery, we added necessary polygons to perform the accuracy assessment. Furthermore, farming system (FS) maps obtained from the 2000/2001 and 2018/2019 classification results were used to map changes that occurred over the study area. The transition between FS classes revealed by comparing the classification results was used to extract unchanged FS during the study period. In order to have significant and robust results of the further trend analysis, we opted to compute trend on unchanged farming systems resulted from the change analysis step. Figure 4 shows a detailed flowchart of the adopted methodology in this study.

Variability and Trend Analysis
From the annual PhM previously retrieved, we calculated their per-pixel temporal mean and coefficient of variation (CV). In order to assess the spatial distribution of PhM showing improvement (positive change) or degradation (negative change), we employed the non-parametric Mann-Kendall (MK) trend test [61] to determine trends on PhM over the OER basin between 2000 and 2019 [34].
The Z statistic follows the standard normal distribution with zero mean and unit variance under the null hypothesis of no trend. A positive Z value indicates an upward trend whereas a negative value indicates a downward trend. Probability (p) represents the measure for evidence to reject the null hypothesis, and positive p values show a negative trend whereas positive p shows a positive trend. We calculated the MK trend test separately for every pixel for the start and end of season, the length of season, and the great integral over the OER basin. These parameters are crucial indicators of seasonal productivity. The Mann-Kendall statistical trend Z was determined as follows: where n is the length of time series data, X k and X j are the observations at k and j time, respectively, and polygons to perform the accuracy assessment. Furthermore, farming system (FS) maps obtained from the 2000/2001 and 2018/2019 classification results were used to map changes that occurred over the study area. The transition between FS classes revealed by comparing the classification results was used to extract unchanged FS during the study period. In order to have significant and robust results of the further trend analysis, we opted to compute trend on unchanged farming systems resulted from the change analysis step. Figure 4 shows a detailed flowchart of the adopted methodology in this study.  The probability linked to the Mann-Kendall statistic "S" and the selected n-data were determined to quantify the level of significance of the trend. Var(S) was calculated, and then the normalized test statistic Z was computed using the following equations: where t varies over the set of tied ranks and ft is the number of times that the rank t appears (i.e., frequency). The equation used to calculate the Mann-Kendall significance Z-score is as follows: where Var(S) is the variance of the dataset and n is the number of data points. The equation used to estimate the Theil-Sen (TS) median trend is: This robust non-parametric trend operator is highly recommended for assessing the rate of change in time-series data. It is calculated by determining the slope between every pairwise combination and then finding the median value. Using the Mann-Kendall test and Theil-Sen median trend analysis, the trends in PhM were described accordingly, and the significance level of the changes in NDVI trends was determined using the Z-score at p-value below 0.1 significance level.

NDVI and Rainfall Time Series Analysis
Given the importance of the NDVI in monitoring vegetation cover, the spatiotemporal variation of this index was assessed based on reference data collected from fieldwork over farming system types between the 2000/2001 and 2018/2019 cropping seasons and the average rainfall over the OER basin for the same period ( Figure 5).
A visual analysis of the temporal NDVI values shows different patterns associated with each farming system type. These patterns were characterized by a specific range of NDVI values and represent the seasonality and the growth cycle of each farming system. The NDVI values for IPC range between 0.45 and 0.9, which is an indicator of the permanent photosynthetic activity of perennial systems ( Figure 5). The studied irrigated tree crops are generally carried out in intensive systems. Owing to their received quantities of water and their long life cycle, this farming system shows a unique characteristic of persistence during the growing season from 2000 to 2019. The IAC class has some identical patterns to the IPC farming system such as the water supply and the high value of production. The IAC farming system shows NDVI values that range from 0.3 to 0.8. The NDVI values start increasing in September and decrease in May of each agricultural season. Their variability in length of the growing season is mainly related to sowing timing, the management of agricultural land, and other plant physiological disease problems.
For the RA farming system, the NDVI values varied from 0.2 to 0.7 for wet seasons and between 0. In general, their NDVI response to rainfall is systematic. Indeed, agricultural lands in the RA farming system are entirely dependent on rainfall, which varies in amount from one year to the next and directly affects agricultural productivity. Indeed, in semi-arid regions, this relationship between rainfall and NDVI has been demonstrated in numerous research studies [62][63][64]. Thereby, the RA farming system shows high variability in terms of the start and end of season and the length of season. This variability is linked to the strong dependency of this farming system on the climatic conditions, especially the rainfall amount. The vegetative activity begins earlier in the case of a wet season and later for a dry season, and the same pattern could be observed for the end of the vegetative activity. A similar pattern could be observed in the FA farming system, where the climatic conditions have more effects on the vegetative activity. These effects of rainfall could be well spotted in the 2009/2010 cropping season over FA and RA farming systems, and we can precisely observe the start of the vegetative activity after receiving an important amount of rainfall ( Figure 5). Through the coming sections, the paper explores more the behavior of farming systems using trend analysis techniques.

NDVI and Rainfall Time Series Analysis
Given the importance of the NDVI in monitoring vegetation cover, the spatiotemporal variation of this index was assessed based on reference data collected from fieldwork over farming system types between the 2000/2001 and 2018/2019 cropping seasons and the average rainfall over the OER basin for the same period ( Figure 5). A visual analysis of the temporal NDVI values shows different patterns associated with each farming system type. These patterns were characterized by a specific range of NDVI values and represent the seasonality and the growth cycle of each farming system. The NDVI values for IPC range between 0.45 and 0.9, which is an indicator of the permanent photosynthetic activity of perennial systems ( Figure 5). The studied irrigated tree crops are generally carried out in intensive systems. Owing to their received quantities of

Spatial Patterns of Phenological Metrics
To mutually characterize PhM behavior (i.e., TSOS, TEOS, LOS, and GINT) and their variability, we combined them with their coefficients of variation (CVs), which are used as a measure of reliability, into the bivariate maps shown in Figure 6. These bivariate maps simultaneously provide a spatial representation of (1) where high or low variability in PhM is expressed for each farming system and (2) the risk involved in the future management and seeding in these farming systems. This combination of factors can assist in developing knowledge about which farming systems can be cultivated in an average year while considering that in such regions, farming system choices are largely based on the potential for attaining good yields and the risk of season failure [65].  Figure 6D shows that areas located inside the irrigated perimeters with a start of the season (TSOS) that occurred between February and June generally have a CV below 0.1. Encouraged by the availability of irrigation water throughout the year, farmers in these areas are diversifying and intensifying their agriculture, thus putting a larger area under crops with high added value, whatever their water consumption [66]. Outside the irrigated perimeters, especially the center of the study area, strong variability is apparent (i.e., variability above 0.2). This high variability is more accentuated by moving towards the mountainous zones. Farmers prefer not to invest excessively in these zones in order to avoid the risks associated with drought, related mainly to the lower stability of TSOS and, therefore, the rainfall amount and distribution. Areas outside the irrigated perimeters with TSOS between August and September have CV values between 0.1 and 0.3 for practically all locations. Nonetheless, the declining availability of water during the studied period causes a decrease in productivity and leads farmers to opt for crops that require less water but offer good margins in order to maximize their profits and avoid the risk of season failure. Other regions outside the irrigated perimeters with a TSOS between days 242 and 365 of the year have a CV below 0.15. Since these areas may present a low-risk investment opportunity, they could be a good choice for the implementation of new farming systems. Towards the western part of the basin, low variability is apparent ( Figure  6D).
On the other hand, areas located outside the irrigated perimeters with an end of season (TEOS) that occurred between the day of year (DOY) 33 and 177 express two degrees of variability, the first with a CV below 0.15 and the second with a CV above 0.35. The majority of the area that has a TEOS between DOY 178 and 241 has low variability, with  Figure 6D shows that areas located inside the irrigated perimeters with a start of the season (TSOS) that occurred between February and June generally have a CV below 0.1. Encouraged by the availability of irrigation water throughout the year, farmers in these areas are diversifying and intensifying their agriculture, thus putting a larger area under crops with high added value, whatever their water consumption [66]. Outside the irrigated perimeters, especially the center of the study area, strong variability is apparent (i.e., variability above 0.2). This high variability is more accentuated by moving towards the mountainous zones. Farmers prefer not to invest excessively in these zones in order to avoid the risks associated with drought, related mainly to the lower stability of TSOS and, therefore, the rainfall amount and distribution. Areas outside the irrigated perimeters with TSOS between August and September have CV values between 0.1 and 0.3 for practically all locations. Nonetheless, the declining availability of water during the studied period causes a decrease in productivity and leads farmers to opt for crops that require less water but offer good margins in order to maximize their profits and avoid the risk of season failure. Other regions outside the irrigated perimeters with a TSOS between days 242 and 365 of the year have a CV below 0.15. Since these areas may present a low-risk investment opportunity, they could be a good choice for the implementation of new farming systems. Towards the western part of the basin, low variability is apparent ( Figure 6D).
On the other hand, areas located outside the irrigated perimeters with an end of season (TEOS) that occurred between the day of year (DOY) 33 and 177 express two degrees of variability, the first with a CV below 0.15 and the second with a CV above 0.35. The majority of the area that has a TEOS between DOY 178 and 241 has low variability, with a CV below 0.15. Inside the irrigated perimeters, low variability is expressed, with a CV below 0.15 for areas having the TEOS between DOY 178 and 241 and between 242 and 273. The high variability in terms of CV observed for the TSOS and TEOS could be related to the dependence of the start of a cropping season to climatic conditions-in rainfed farming systems, the season could not start without the first rainfall. Adversely, the end of the cropping season depends on the climatic conditions but also depends on the physiological properties of the plant and its persistence ( Figure 6C).
Considering the phenological metric LOS, the results of this research show that the variability is inversely proportional to the value of the length of the season, without regard to particulars or exceptions. Hence, the areas with an LOS value below 160 days generally have a CV value above 0.35, and the areas with an LOS value between 160 and 192 days have CV values between 0.15 and 0.35, especially in the center of the study area. However, these areas may express a deficit in agricultural production during years with low rainfall amounts. The low variability is mostly located inside the irrigated perimeters where the LOS is above 224 days. These zones shows CV values of below 0.15. Various areas in the western part of the basin also show high CV values which are above 0.28 ( Figure 6B).
Finally, by combining the fourth phenological metric explored during this study with its CV values (GINT), our results can help in the characterization of the biomass variability of farming systems encountered in this study area ( Figure 6A). The areas having GINT values below four have CV values above 0.35. These zones are generally rainfed areas and fallow, where climatic conditions have an effect. Generally, these areas are not recommended for the seeding of plants with high values of production. Instead, areas with GINT values above eight have CV values between 0.15 and 0.25. This moderate variability in biomass productivity could be considered as normal, considering the differences in crop types over the irrigated parameters (i.e., alfalfa, sugar beet, wheat, etc.). Other zones in the OER basin remain with low to moderate biomass productivity with high variability of the GINT metric.

Determination of Unchanged Farming Systems' Area
RF classification offers a powerful algorithm to classify the spatial patterns of farming systems in the study area. In this study, the classification of farming systems was produced from the implementation of the Random Forest (RF) classifier based on PhM. The classifier was applied and the accuracy assessment results are summarized and provided in the Supplementary Materials ( Figures S5-S7). In terms of the individual accuracy, the results present higher overall accuracy and much more balanced producer/user accuracy. The overall accuracy for RF was 97% and 93% with a kappa value of 0.95 and 0.91 for the 2018/2019 and 2000/2001 cropping seasons, respectively. In general, all farming system classes achieved over 90% user accuracy. The RF classifier also produced over 90% producer accuracy for most farming system classes. The change detection results from 2000 to 2019 reveal an important dynamic between the different FSs in the studied region. The transition between FSs from 2000/2001 to 2018/2019 shows that there has been a marked change in FS classes over the study area (Figure 7). To evaluate trends over the study area, we only considered the unchanged classes during the studied period.

Trend Analysis Results
This research further analyzed the inter-annual trends for the four PhM (i.e., SOST, EOST, LOS, and GINT). The Mann-Kendall test was used to identify significant trends at 90% confidence level and Sen's slope estimator was applied to compute the magnitude of trends for the period between 2000 and 2019, alternating between increasing and decreasing trends over the study area.

Trend Analysis Results
This research further analyzed the inter-annual trends for the four PhM (i.e., SOST, EOST, LOS, and GINT). The Mann-Kendall test was used to identify significant trends at 90% confidence level and Sen's slope estimator was applied to compute the magnitude of trends for the period between 2000 and 2019, alternating between increasing and decreasing trends over the study area. For the start of season (TSOS), Figure 8 displays its trends as obtained using Mann-Kendall and Sen's slope tests. Positive and negative values refer to delayed and advanced dates of TSOS, respectively. Most areas in the OER basin did not show a significant trend during the past twenty years for the TSOS, LOS, and GINT metrics (p > 0.1). On the contrary, the TEOS shows a significant trend over most parts of the study area (p < 0.1). In this section, only significant trend pixels are considered for description and analysis. A significant negative TSOS trend was discovered across small parts of the Doukkala irrigated perimeter, the Tadla irrigated perimeter, and in parts of the rainfed area. For the pixels with negative trends, over irrigated perimeters, Sen's slope shows an average decrease in TSOS of approximately −0.2 Most areas in the OER basin did not show a significant trend during the past twenty years for the TSOS, LOS, and GINT metrics (p > 0.1). On the contrary, the TEOS shows a significant trend over most parts of the study area (p < 0.1). In this section, only significant trend pixels are considered for description and analysis. A significant negative TSOS trend was discovered across small parts of the Doukkala irrigated perimeter, the Tadla irrigated perimeter, and in parts of the rainfed area. For the pixels with negative trends, over irrigated perimeters, Sen's slope shows an average decrease in TSOS of approximately −0.2 day/year; outside the irrigated perimeters, the TSOS shows a decrease between −0.2 and −0.4 day/year. The early start of the agricultural season in these zones could explain these results, especially the irrigated annual crops where we found different crop types and irrigation systems. Hence, this implies that over the studied period of 20 years, the TSOS has been delayed by more than 6 days (−6 day/20 years). Other regions of the OER basin showed a non-significant trend. These regions are not interpreted since the significance level is above 0.1.
On the other hand, we found a significant positive trend for the TEOS across large parts of the OER basin and in a part of the Tadla irrigated perimeter. However, we found a significant negative trend in the rainfed and fallow farming systems (Figure 9). For the pixels with positive trends, Sen's slope shows an average increase in TEOS above 0.2 day/year in the TIP and between 0 and 0.2 in the rainfed area. The TEOS shows a decrease between 0 and −0.2 days/year in numerous parts of the OER basin, especially in the southern part. These results could be explained by the satisfaction of plants requirements in the irrigated zones, whereas in other areas, the negative trend is essentially linked to climatic conditions and poor land management. This signifies that over the studied period of 20 years, the TEOS has changed by more than 3 days/20 years in irrigated areas and by −3 days/20 years in rainfed and fallow areas.
After the trend analysis of the TEOS metric, we directed our studies towards the analysis of trends in the length of season (LOS) (Figure 10) to investigate the link between TEOS and TSOS with the LOS, since the LOS is the result of the difference between the TEOS and the TSOS. The results show a significant positive trend in LOS over the south parts of the Doukkala irrigated perimeter and in the north parts of the Tadla irrigated perimeter. In addition, we found a significant negative trend in the rainfed and fallow farming systems. For the pixels with positive trends, Sen's slope shows an average increase in TEOS above 0.2 day/year in the TIP and between 0 and 0.2 in the rainfed area. The TEOS shows a decrease between 0 and −0.2 days/year in numerous parts of the OER basin, especially in the southern part. These results could be explained by the satisfaction of plants requirements in the irrigated zones, whereas in other areas, the negative trend is essentially linked to climatic conditions and poor land management. This signifies that over the studied period of 20 years, the TEOS has changed by more than 3 days/20 years in irrigated areas and by −3 days/20 years in rainfed and fallow areas.
After the trend analysis of the TEOS metric, we directed our studies towards the analysis of trends in the length of season (LOS) (Figure 10) to investigate the link between TEOS and TSOS with the LOS, since the LOS is the result of the difference between the For the pixels with positive trends, Sen's slope gives an average increase in LOS above 0.1 day/year in the Tadla irrigated perimeter and approximately between 0 and 0.1 days/year in the DIP. The LOS shows a decrease between 0 and −0.1 day/year in numerous parts of the OER basin. This signifies that over the studied period of 20 years, the LOS has changed by more than 1.6 days/20 years in irrigated areas and by −1.6 days/20 years in rainfed and fallow areas.
Comparing these results with trends in TSOS and TEOS indicates that increases in TEOS and decreases in TSOS dates are mostly responsible for the positive LOS trends in irrigated perimeters of Doukkala and Tadla. Negative LOS trends in the rainfed area and fallow farming systems are related to a delay in TSOS dates and an early TEOS.
Finally, considering the GINT metric, the trend of this index was obtained using Mann-Kendall and Theil-Sen median trend tests ( Figure 11). We found significant positive trends over large farming systems. These areas are specially located in the irrigated perimeters of Tadla and Doukkala, in addition to some zones with improved agricultural management practices. For the pixels with positive trends, Sen's slope gives an average increase in LOS above 0.1 day/year in the Tadla irrigated perimeter and approximately between 0 and 0.1 days/year in the DIP. The LOS shows a decrease between 0 and −0.1 day/year in numerous parts of the OER basin. This signifies that over the studied period of 20 years, the LOS has changed by more than 1.6 days/20 years in irrigated areas and by −1.6 days/20 years in rainfed and fallow areas.
Comparing these results with trends in TSOS and TEOS indicates that increases in TEOS and decreases in TSOS dates are mostly responsible for the positive LOS trends in irrigated perimeters of Doukkala and Tadla. Negative LOS trends in the rainfed area and fallow farming systems are related to a delay in TSOS dates and an early TEOS.
Finally, considering the GINT metric, the trend of this index was obtained using Mann-Kendall and Theil-Sen median trend tests ( Figure 11). We found significant positive trends over large farming systems. These areas are specially located in the irrigated perimeters of Tadla and Doukkala, in addition to some zones with improved agricultural management practices.
Generally, significant positive trends are located in the western and eastern parts of the OER basin. The central region is dominated by insignificant negative trends, while stable farming systems are dispersed over the study area. Irrigated perimeters show an increase in GINT by 0.2/year, while the rainfed area shows an increase between 0 and 0.1/year. Generally, the trend results reveal an important variation between the different FSs in the studied region.

Discussion
In this study, we have demonstrated that the phenological pattern of vegetation cover across different farming systems and across different regions over the last 20 years has an important implication in vegetation dynamics and climate change. It provides an insight of the vegetation cover status of a region affected by climate risks. To our knowledge, this study is the first to demonstrate the ability of annual phenological metrics (PhM) and RF modeling to map and characterize changes over the major farming systems in Morocco. The use of TIMESAT for NDVI time-series processing and analysis has provided a more comprehensive investigation of the NDVI behavior from a phenological perspective. As far as we know, this study is the only one that has explored phenological farming systems' variability over North Africa at 250-m spatial resolution.
Overall, the accuracies obtained in this study for the classification of farming systems were reliable and consistent with those revealed for other MODIS-based land cover classifications [30,31]. For comparison, the resulting overall accuracy of the RF classifications obtained over the OER basin scale (i.e., 90-96%) is similar to that obtained using MODIS time series for mapping cropping intensity in China at a regional scale (i.e., 92%; [31]). On the side of classification performance, the confusion between rainfed area and fallow, which was the main source of error, reflects the strong similarity between these two farming systems, especially when the climatic conditions are critical. Consequently, our findings have confirmed the advantage of using a combination of PhM from MODIS time series and an RF classifier to discriminate between farming systems [67]. In addition to the

Discussion
In this study, we have demonstrated that the phenological pattern of vegetation cover across different farming systems and across different regions over the last 20 years has an important implication in vegetation dynamics and climate change. It provides an insight of the vegetation cover status of a region affected by climate risks. To our knowledge, this study is the first to demonstrate the ability of annual phenological metrics (PhM) and RF modeling to map and characterize changes over the major farming systems in Morocco. The use of TIMESAT for NDVI time-series processing and analysis has provided a more comprehensive investigation of the NDVI behavior from a phenological perspective. As far as we know, this study is the only one that has explored phenological farming systems' variability over North Africa at 250-m spatial resolution.
Overall, the accuracies obtained in this study for the classification of farming systems were reliable and consistent with those revealed for other MODIS-based land cover classifications [30,31]. For comparison, the resulting overall accuracy of the RF classifications obtained over the OER basin scale (i.e., 90-96%) is similar to that obtained using MODIS time series for mapping cropping intensity in China at a regional scale (i.e., 92%; [31]). On the side of classification performance, the confusion between rainfed area and fallow, which was the main source of error, reflects the strong similarity between these two farming systems, especially when the climatic conditions are critical. Consequently, our findings have confirmed the advantage of using a combination of PhM from MODIS time series and an RF classifier to discriminate between farming systems [67]. In addition to the influence of climate conditions on phenological changes, changes in crop type were also an important determinant of phenological changes and trends in the OER basin. For this reason, in this study, we used a mask layer to mask only unchanged areas and ejected the effect of annual land use change over farming systems.
In terms of PhM variability and trends, our study revealed substantial variability over the studied PhM. Generally, TSOS and TEOS show high and low variability inside and outside the irrigated perimeters, respectively ( Figure 6). The TSOS inside irrigated perimeters displayed earlier onset trends than other zones, which could be due to the sowing date and irrigation (Figure 7). On the other hand, delayed TEOS could be explained by the variability of harvest time and biological factors and the climatic conditions that occurred during the cropping seasons. These delays and early onsets in time of occurrence between TSOS and TEOS were translated to a longer cropping season over the irrigated perimeters ( Figure 9). The remarkable shift expressed by GINT and LOS over the irrigated perimeters is the result of an intensification plan and a change of agricultural practices with the encouragement of farmers to plant perennial crops. The question raised here is how sustainable are the water and soil resources supporting this mode of production under the growing pressure applied?
Our study has also revealed that annual variability outside the irrigated perimeters is strong, especially in the center of the study area, where fallow and rainfed farming systems extend. This high variability can be explained by the changes in rainfall amount and distribution received in these regions during the cropping season and the changes to the agricultural production potential. As identified by Ouatiki et al. [68], climatic stations located in these regions experienced a significant decreasing trend and received a low rainfall amount during the season. Other factors leading to the high variability of PhM in these regions were the depletion of groundwater and the quality of agricultural management, since these regions are considered as rainfed and fallow areas where the accessibility to irrigation supplies is absent [4]. A similar situation was observed in the variability of the GINT metric ( Figure 6A) as the same regions experienced high variability and low values of GINT, which signify low biomass production and an unbalanced field for agricultural planning. Regarding the trend results, GINT increased during the studied period by 0.1/year in most areas of the basin, while for certain regions of FA and RA, it decreased by −0.1/year ( Figure 11). Therefore, we believe that all of these factors affect the productivity of these two farming systems, which are translated in the short cropping season and the low biomass production ( Figure 6A,B, respectively).
The wide spatial distribution of annual variability of PhM found in this study is reliable and comparable with the outcomes from previous studies [65,69,70]. For instance, Vrieling et al. [65] identified that a higher variability of Length of Growing Period (LGP)represented in our study by the length of season metric (LOS)-is generally found in arid and semi-arid areas, with coefficients of variation above 0.25 ( i.e., LOS in our study varied from 0.15 to 0.35). This slight difference between our studies could be due to the higher spatial resolution of MODIS data (i.e., 250 m). Contrarily, in the irrigated perimeters, where land receives significant amounts of water and fertilizer, the variability of PhM remained lower. Furthermore, the variability in crop types may have had some influence on the variability of PhM. However, this hypothesis is uncertain as our study focused only on unchanged farming systems of the OER basin.
Weather extremes and consecutive droughts in this region strongly affected the vegetation cover dynamics and resulted in adaptations of farming system management in response to climatic variation. Socio-economic policies and improving farm practices were also dominant drivers of farming system changes in addition to climate conditions. These factors could all have affected the variability of the start, the end, and the productivity of the cropping season, resulting in phenological changes over time. A good example of agricultural planning that made great changes in the agricultural sector in Morocco is the Green Morocco Plan strategy, which aimed to encourage tree-growing extensions that led to a high variability over some parts of the perennial farming system.
Notably, our findings showed that the variability of PhM is varied and related to farming system types. The results from this study are expected to constitute a background to other research about drought monitoring and desertification studies in arid regions. The PhM datasets and the trend results could be integrated with climate data to perform an estimation of crop water requirements and offer a tool for managers and stakeholders to make decisions for the extension of agricultural areas according to the available water resources in a context of water stress.

Conclusions
Finally, variation in phenological metrics over FSs was estimated at the OER basin level during 2000-2019 seasons using MODIS NDVI data and trend analysis tests. Our study findings are the following: (1) Over-irrigated perimeters' (TIP and DIP) mean LOS, GINT, and TEOS values showed low variability. On the other hand, moderate variability was observed for the mean TSOS values during the studied period. Contrary to the irrigated zones, PhM over the rainfed and fallow farming systems showed high variability. This variability over RA and FA is justified by the irregularity of rainfall amounts received over these farming system areas.
(2) Trends over farming systems are not uniform at the OER basin level. Most areas in the OER basin did not show a significant trend during the past 20 years for the TSOS, LOS, and GINT metrics (p > 0.1). Contrary to the TEOS, where a significant trend was observed (p < 0.1), TSOS shows early onset over the IPC and IAC (i.e., 0.2 days/year), while over RA and FA, it was delayed by −0.2 days/year, especially in the center of the basin. Other regions of the FA and RA showed extended TSOS by 0.2 days/year. TEOS shows early onset (i.e., −0.4 days/year) over the FA, RA, and part of the IPC. Other regions of the basin showed a TEOS extended by 0.2 days/year. LOS generally slightly increased over the farming systems, except in particular zones of the FA, and did not advance markedly during the study period. GINT increased during the studied period by 0.1/year in most areas of the basin, while for certain regions of FA and RA, it decreased by −0.1/year.
The investigation of changes in PhM over twenty years (19 cropping seasons) proved the ability of these metrics to characterize the spatiotemporal changes over the OER basin. Nevertheless, the need to take into account the perceptions and opinions of local populations is essential in order to reduce the process of vegetation cover degradation and to better manage natural resources.
Supplementary Materials: The following are available online at https://www.mdpi.com/2072-429 2/13/4/578/s1. Figure S1: Bivariate map showing, simultaneously, the start of the season time and its coefficient of variation. Figure S2: Bivariate map showing, simultaneously, the end of the season time and its coefficient of variation. Figure S3: Bivariate map showing, simultaneously, the length of season and its coefficient of variation. Figure S4: Bivariate map showing, simultaneously, the great integral metric and its coefficient of variation. Figure   satellite data" at the Center for Remote Sensing Applications (CRSA), Mohammed VI Polytechnic University (UM6P). We also thank the academic editor and three anonymous reviewers for their helpful comments which greatly improved earlier versions of the manuscript. We would like to acknowledge the Department of Strategy and Statistics (DSS) from the Ministry of Agriculture in Morocco for the financial support as part of the agreement between National Institute of Agronomic Research (INRA) and DSS.

Conflicts of Interest:
The authors declare no conflict of interest.