Article

# Assessing Hydrological Vulnerability to Future Droughts in a Mediterranean Watershed: Combined Indices-Based and Distributed Modeling Approaches

Youssef Brouziyne <fnref n="1" />,*, Aziz Abouabdillah <fnref n="2" />, Abdelghani Chehbouni <fnref n="1" />,<fnref n="3" />, Lahoucine Hanich <fnref n="4" />,<fnref n="5" /> ID,
Karim Bergaoui <fnref n="6" />,<fnref n="7" />,†, Rachael McDonnell <fnref n="7" />,† ID and Lahcen Benaabidate <fnref n="8" /> ID

Received: 15 May 2020; Accepted: 20 July 2020; Published: 19 August 2020

<figure><img src="image_1.png" /></figure>

check for
updates

**Abstract:** Understanding the spatiotemporal distribution of future droughts is essential for effective water resource management, especially in the Mediterranean region where water resources are expected to be scarcer in the future. In this study, we combined meteorological and hydrological drought indices with the Soil and Water Assessment Tool (SWAT) model to predict future dry years during two periods (2035–2050 and 2085–2100) in a typical Mediterranean watershed in Northern Morocco, namely, Bouregreg watershed. The developed methodology was then used to evaluate drought impact on annual water yields and to identify the most vulnerable sub-basins within the study watershed. Two emission scenarios (RCP4.5 and RCP8.5) of a downscaled global circulation model were used to force the calibrated SWAT model. Results indicated that Bouregreg watershed will experience several dry years with higher frequency especially at the end of current century. Significant decreases of annual water yields were simulated during dry years, ranging from –45.6% to –76.7% under RCP4.5, and from –66.7% to –95.6% under RCP8.5, compared to baseline. Overall, hydrologic systems in sub-basins under the ocean or high-altitude influence appear to be more resilient to drought. The combination of drought indices and the semi-distributed model offer a comprehensive tool to understand potential future droughts in Bouregreg watershed.

**Keywords:** drought indices; SWAT; spatiotemporal assessment; water yield; Bouregreg watershed

## 1. Introduction

The global climate has changed relative to the pre-industrial period, with many forms of evidence showcasing the impacts on organisms and ecosystems, as well as human systems and well-being [1–4]. Besides the increase in land and ocean temperatures, the manifestations of global climate change

<footnote n="1">¹ Mohammed VI Polytechnic University (UM6P), International Water Research Institute, Benguerir 43150, Morocco; Abdelghani.chehbouni@um6p.ma</footnote>

<footnote n="2">² Department of Agronomy and Plant Breeding, Ecole Nationale d'Agriculture de Meknès, Meknes 50001, Morocco; Aabouabdillah@enameknes.ac.ma</footnote>

<footnote n="3">³ Institut de Recherche Pour le Développement (IRD), Unité Mixte de Recherche (UMR), Centre D'études Spatiales de la Biosphère (Cesbio), 31401 Toulouse, France</footnote>

<footnote n="4">⁴ L3G Laboratory, Department of Earth Sciences, Faculty of Sciences & Techniques, University of Cadi Ayyad, Marrakech 40000, Morocco; L.hanich@uca.ma</footnote>

<footnote n="5">⁵ Mohammed VI Polytechnic University (UM6P), Center for Remote Sensing Applications, Benguerir 43150, Morocco</footnote>

<footnote n="6">⁶ Acquatec Solutions FZCO, Dubai-Digital Park, Dubai PO 6009, UAE; Bergaoui@gmail.com</footnote>

<footnote n="7">⁷ International Water Management Institute, 00054 Rome, Italy; R.Mcdonnell@cgiar.org</footnote>

<footnote n="8">⁸ Laboratory of Functional Ecology and Environmental Engineering, Faculty of sciences & Techniques, University of Sidi Mohamed Ben Abdellah, Fez 30000, Morocco; lahcen.benaabidate@usmba.ac.ma</footnote>

<footnote>* Correspondence: Youssef.brouziyne@gmail.com or Youssef.brouziyne@um6p.ma</footnote>

<footnote>† Formerly: International Centre for Biosaline Agriculture, Academic City, Al Ain Road, Dubai 14660, UAE.</footnote>
include the increase in frequency, intensity, and/or amount of extreme weather events (such as heavy precipitation and drought) [1,4].

Due to observational uncertainties and difficulties in distinguishing decadal-scale variability in drought from long-term trends, scientists are always cautious about the attribution of changes in drought to anthropogenic influence [5]. Nevertheless, most of the recent research assumes that climate change has led to an increase of frequency and intensity of droughts in some areas (Mediterranean, West Asia, South America, Northeastern Asia, etc.) [6,7].

Drought is generally defined as a "prolonged absence or marked deficiency of precipitation", a "deficiency of precipitation that results in water shortage for some activity or for some group", or a "period of abnormally dry weather sufficiently prolonged for the lack of precipitation to cause a serious hydrological imbalance" [8]. All these definitions of drought highlight the importance of precipitation, evapotranspiration, and runoff and their interactions in drought caused by climatic factors [9]. Water deficit during droughts generally affects the whole hydrological cycle, reducing the volume of water at many levels (water bodies, streamflow, and groundwater); water quality may decline too due to low-flow regimes and water temperature increase under drought conditions [10]. By affecting water availability and quality, drought can have significant effects on many other vital systems and services (drinking and sanitation water supply, food security, health and well-being, economic development, etc.) [11]. Effective water resource management under drought risks relies on forecast and the good understanding of characteristics and features of potential droughts [12].

For drought monitoring and understanding, many statistical probability-based indices have been created and have been continually developed until nowadays [13–15]. The usefulness of drought indices in forecasting, tracking, and measuring the recovery from droughts is confronted, in some situations, to the need for a large series of multi-variables over large areas [16,17]. The lack of consideration of interactions between many physical components [18] and the failure to include characteristics of some region-specific factors [10] might lead to weaknesses in investigating droughts. Under these situations, combining drought indices and distributed hydrological modeling offers a more complete methodology base to scrutinize hydrological drought [10,16,19–21]. Many recent studies used this combination to understand past droughts and to project future water resources under potential drought impacts in many regions worldwide. Hoerling and Eischeid [22] used large scale hydrological models and drought indices and have investigated future droughts and the impact of climate change on droughts in the southwestern and the central United States. Leng et al. [23] used the Variable Infiltration Capacity (VIC) model [24] in order to investigate future changes in China in terms of hydrology, and then they assessed climate change impact on drought patterns from many perspectives by calculating the SPI (Standardized Precipitation Index), SRI (Standardized Runoff Index) and SSWI (Standardized Soil Moisture Index) through VIC help. Kang and Sridhar [16] confirmed the ability of the Soil and Water Assessment Tool (SWAT) [25] and VIC models together with a set of drought indices to forecast and study future drought evolution in a multi-purpose watershed in USA.

In the Mediterranean region, droughts are considered one of the most severe and frequent natural hazards [26]. In terms of the future, most of the climatic models used in the last assessment report of the IPCC (Intergovernmental Panel on Climate Change) predict an increase of aridity in the Mediterranean region for the rest of the current century as fewer precipitations are expected throughout the whole year, with a clear deficit in summer season [5,6].

Because of its location in the Mediterranean region, Morocco has always suffered from extreme weather events such as drought. Although drought can occur in all climatic zones, arid and semi-arid regions (such as several parts of Morocco) are the most vulnerable because of rainfall variability and uneven temporal water distribution [27,28]. Many studies confirmed that the frequency of dry years has increased since the early 1980s in Morocco [29]. The frequency of dry years has changed from 1 year in 15 normal years in the 1930s, 1940s, 1950s, 1960s, and 1970s, 1 long dry year in 3 normal years in the last three decades [30,31].
Many severe drought episodes have been recorded in the recent history of the country, with clear impacts on many aspects of the environment, population well-being, and economic development [32]. For example, the drought that occurred in 1999–2000 was responsible for the decline in national economic growth and around USD 900 million in losses, with more than 275,000 persons being affected [29,31]. Other previous droughts in the 1990s led to a fall of national gross domestic product significantly (from –1.5% to –7.6%), pushing the government to import more cereals as national crop productions dropped down [31,33].

The Bouregreg watershed (BW) is located in the northwestern part of Morocco. With a total population exceeding 2.4 million and various natural development assets, this watershed is considered one of the most important watersheds in the country [34]. Medium and small urban centers are located within the BW territory. Farming and forestry are the major activities that are held over the watershed, with high economic as well as social importance for the local population [27]. The Sidi Mohammed Ben Abdellah dam is located a few kilometers from the outlet of BW, being considered one of the major five dams of Morocco as it secures most of drinking and sanitation water demand for all surrounding urban centers, including Rabat, the capital [34,35].

BW is a typical Mediterranean watershed, with large development potentialities, and having various challenges related to the growing water demand and uncertain climate in the future; however, no previous studies with future drought impacts on water as core scope have been performed except some isolate assessments of current water balance [36,37] and investigations of historic climate change trends [35]. Therefore, a comprehensive and integrated assessment of potential future drought impacts on water resources in BW is needed.

Thus, the primary objective of this study was to combine indices-based techniques and a large hydrological modeling approach, for the first time in Morocco, in order to forecast potential reactions of the hydrological system to future droughts in such a strategic watershed in Morocco. The specific objectives of the present study are (1) to combine drought indices and a semi-distributed model to forecast future drought episodes while taking into consideration the physical features of BW, (2) to assess potential changes in the hydrological systems under droughts, and (3) to identify the most vulnerable sub-basins of BW to drought impacts. This study aimed to highlight drought evolution in time and space in the watershed unit under increasing demand for water and provide a methodological basis for further water management systems under potential extreme weather events in the Mediterranean region.

## 2. Materials and Methods

### 2.1. *Study Area*

The study watershed is one of the main catchments in Morocco, being well characterized by topographic diversity, and hosting the main Moroccan climatic zones, as well as the main crops and vegetation of the country. The Bouregreg watershed is located in the northwestern part of Morocco (Figure 1), with a total surface area of 9656 km². Topography of the BW is divided into three different classes: mountainous area in the northeast (Oulmes area) with altitudes ranging from 1200 to 1726 m; relatively medium altitude area in the central and southeastern part of the study watershed (800 to 1200 m); and the western part, which is the lowest plateau of the watershed (0 to 800 m) [34].
<figure><img src="image_2.png" /><figcaption>**Figure 1.** topography, and streamflow network of Bouregreg watershed (BW) and the used gauging stations.</figcaption></figure>

BW is considered a Mediterranean watershed with different climatic zones and a clear ocean influence. The extreme western part of the watershed receives a generally cumulative rainfall of 500 mm with a mean annual temperature of 18 °C, which makes this area a sub-humid climate region. The central part of the basin is considered a semi-arid region and receives an annual rainfall of less than 400 mm with an annual mean temperature of 18 °C. The southeastern part is under a continental influence with a cumulative rainfall of 400 mm and an annual mean temperature of 15 °C, classifying this region as arid. The extreme northeastern part of the watershed, the highlands of Oulmes, is a mountainous sub-humid area that receives generally more than 600 mm rainfall, with a mean annual temperature of 15 °C [37]. The annual potential water resources of the Bouregreg watershed is 720 million m³, coming mostly from surface water (more than 90%). The Sidi Mohamed Ben Abdellah dam, located a few kilometers upstream of the watershed outlet, has a capacity of 1500 Million m³ and secures 250 to 450 million m³ of sanitation and drinking water for Rabat city and all the surrounding urban areas. Water inflow in the watershed is very linked to the rainfall regime, which is highly non-uniform due to variance between flood periods (December to April) and dry season (summer) [34,38].

The stream network of Bouregreg watershed is quite developed, being organized around three main rivers: Korifla, Grou, and Bouregreg rivers. Both of the two major rivers (Bouregreg and Grou) have origins from the same area, M'riret, which is located at the base of Middle Atlas Mountains [34,39] (Figure 1).

In terms of land use (Figure 2), bare lands and oak forests cover, respectively, 14% and 24% of the total area of the study watershed. Farming activity (28%) is very important, with the dominant crops being cereals, olive trees, and pulse crops (broad beans), in addition to minor parcels where some irrigated vegetables and grape crops are grown. Range and pasture lands (30%) are distributed across the watershed as husbandry is an important economic activity in the area [34,40].
<figure><img src="image_3.png" /><figcaption>Figure 2. Land use classes in BW.</figcaption></figure>

## 2.2. Hydrological Model

For this study, the SWAT model (Soil and Water Assessment Tool) was used to simulate the hydrological system and its reaction under baseline and future drought conditions in the study watershed. This model offers wide possibilities by integrating different modules (rainfall-runoff, plant growth, sediment transport, pesticides transport, etc.) in a physically based approach and with planning and management change possibilities [41]. This guarantees more reliable scenario simulation in the future and better confidence in the study conclusions in such a complex watershed [42].

A distributed rainfall-runoff model, such as SWAT, divides a watershed into smaller discrete units for which spatial variation of the major physical properties are limited and the hydrological processes can be treated as being homogeneous [43]. In the specific case of SWAT, HRUs (Hydrologic Response Units) are formed from the combination of soil types, land use classes, and slope ranges. All processes covered by SWAT (runoff, crops yields, sediment yield, etc.) are estimated at the HRU scale [44].

SWAT relies on the Soil Conservation Service curve number method [45] to capture the rainfall-runoff relationship. The daily water balance simulated by the SWAT model is based on the following equation [25]:

$$SW_t = SW_0 + (R_{day} - Q_{surf} - E_a - W_{seep} - Q_{gw}) \quad (1)$$

where

SWt: the daily final soil water content (mm H₂O);

SW0: the initial soil water content on day i (mm H₂O);

t: the time (days);

Rday: the amount of precipitation on day i (mm H₂O);

Qsurf: the amount of surface runoff on day i (mm H₂O);

Ea: the amount of evapotranspiration on day i (mm H₂O);

Wseep: the amount of water entering the vadose zone from the soil profile on day i (mm H₂O);

Qgw: the amount of return flow on day i (mm H₂O).

The SWAT version used in this study is the 2012 edition; calibration and sensitivity analysis processes were performed using the SWAT auxiliary program SWAT CUP [46].
SWAT CUP is a calibration, validation, and sensitivity analysis tool for SWAT and can be freely downloaded. This public program offers different choices in terms of optimization programs. In this study, the SUFI-2 (Sequential Uncertainty Fitting Version 2) was used. In SUFI-2, uncertainty in parameters, expressed as ranges (uniform distributions), accounts for all potential sources of uncertainties, such as variables (rainfall, temperatures, etc.), conceptual model, and parameters [47].

Before proceeding with SWAT calibration over BW, global sensitivity analysis tests were carried out using the SUFI-2 algorithm. The objective relies on identifying the most influential model parameters on both hydrology and plant growth in the study watershed. SWAT CUP uses a t-test to identify the relative significance of each parameter on objective function change [46]. For this study, parameters to be tested and their ranges were selected on the basis of professional knowledge on local processes and literature.

During streamflow (m³/s) calibration through SWAT CUP, we looked at three statistical indicators to confirm a maximal concordance between observed and predicted water budget components: the Nash–Sutcliffe coefficient, the determination coefficient, and the percent bias.

With

$$NSE = 1 - \frac{\sum_{i=1}^{n} (O_i - P_i)^2}{\sum_{i=1}^{n} (O_i - O_{\text{avg}})^2} \qquad (2)$$

and

$$R^2 = \left( \frac{\sum_{i=1}^{n} (O_i - O_{\text{avg}}) (P_i - P_{\text{avg}})}{\left[ \sum_{i=1}^{n} (O_i - O_{\text{avg}})^2 \sum_{i=1}^{n} (P_i - P_{\text{avg}})^2 \right]^{0.5}} \right)^2 \qquad (3)$$

and

$$PBIAS = \frac{\sum_{i=1}^{n} (O_i - P_i) 100}{\sum_{i=1}^{n} O_i} \qquad (4)$$

where

O_i: the observed streamflow for time period I;

P_i: the simulated value for the same period;

O_{avg}: the mean of observed streamflow per time period;

P_{avg}: the mean of simulated values for the same period;

n: the number of time periods.

The Nash–Sutcliffe coefficient (NSE) value can range from -∞ to 1, and a perfect fit between the simulated and observed data is reached when NSE value is close to 1 [47]. The determination coefficient, R², indicates how well the variance of observed values are replicated by the model simulations and can vary from 0 to 1, where 0 indicates no correlation and 1 represents the perfect correlation [48]. The percent bias (PBIAS) identifies average tendency of simulated data to be greater or lesser than their recorded counterparts. PBIAS starts indicating satisfactory simulation performance when its absolute value is less than 25% [48]. These goodness-of-fit indicators form the objective function for calibration tests and were supervised during the automatic modification of the sensitive parameters of SWAT modules in the Bouregreg watershed.

For optimal performances of the SWAT model over BW, and due to the importance of farming in the study watershed and its potential influence on water balance, we made specific efforts to calibrate the plant growth module of SWAT in addition to the rain-runoff module. The determination coefficient (R²) was adopted as a statistical indicator of goodness-of-fit between simulated and observed annual crop yields while amending the most influencing parameters on the plant growth process. SWAT calibration helper command was used to perform this task by modifying the influential crop parameters of the concerned crops in the HRUs where this crop is grown.
## 2.3. Drought Indices

According to Sun et al. [9], there are four different drought categories: meteorological, hydrological, agricultural, and socioeconomic droughts. Meteorological drought is generated by rainfall deficit over a region for a period of time [49]; hydrological drought is related to below-normal supply of water in rivers, water bodies, and groundwater [50]; agricultural drought is marked by a period with very low (to null) soil moisture levels leading to crops stress [50,51]; and socioeconomic drought is linked to the impact of drought conditions on the regular fluxes and mechanisms of a given economic commodity [49].

Depending on the scope, drought can be monitored through various statistical probability indices that have been developed for each scope perspective. The Standardized Precipitation Index (SPI) [52] and percent of normal precipitation [53] are some of the widely used meteorological drought indices. Palmer Hydrological Drought Index [54], Baseflow Index [55], and the recently developed Streamflow Drought Index (SDI) [56] are some of the common hydrological drought indices. Although precipitation deficiency is the root-cause of all types of droughts, hydrological drought is generally lagging behind the meteorological drought [57]. This is mainly because it takes some time before insufficient precipitation appears in the various subsurface component of the hydrological system [57,58]. In this study, the SPI was chosen for describing meteorological drought and the SDI was chosen for hydrological drought.

The SPI is probably the most popularly used meteorological drought index where the precipitation data is fitted to gamma distribution [51,59]. Further description of the SPI calculation description is reported by MacKee et al. [52].

SPI is a flexible index that can be calculated at multiple timescales: 3, 6, 9, and 12 months [51]; for example, the SPI-3 reflects the short and medium term moisture conditions, and the SPI-12 reflects the long term precipitation patterns used to monitor streamflow, reservoir, and groundwater [60]. For this study, and since the hydrological systems is the core focus, the SPI-12 was calculated for BW using monthly precipitation.

The SDI is a “simplistic” hydrological index proposed by Nalbantis [56], and it keeps the advantage of effectiveness found in meteorological drought indices such as SPI; the monthly streamflow data is the unique requirement to compute SDI.

This index assumes that a time series of monthly streamflow volumes $Q_{i,j}$ is available, where i denotes the hydrological year and j the month within that hydrological year. In the Mediterranean region, the typical hydrological year begins by October.

Then,

$$V_{i,j} = \sum_{j=1}^{3k} Q_{i,j}, \quad i = 1, 2, \dots, j = 1, 2, \dots, 12 \quad k = 1, 2, 3, 4 \quad (5)$$

where $V_{i,j}$ is the cumulative streamflow for the i-th hydrological year and k-th reference period, k = 1 for October–December, k = 2 for October–March, k = 3 for October–June, and k = 4 for October–September.

On the basis of cumulative streamflow volumes $V_{i,j}$, the Streamflow Drought Index is defined for each reference period k of the i-th hydrological year as follows:

$$SDI_{i,k} = \frac{V_{i,k} - \bar{V}_k}{S_k} \quad i = 1, 2, \dots, k = 1, 2, 3, 4 \quad (6)$$

where $\bar{V}_k$ and $S_k$ are, respectively, the mean and the standard deviation of cumulative streamflow volumes of reference period k, as these are estimated over a long period of time [56].

For both SPI and SDI, a drought starts in the month when their values fall below -1 and it ends when they return to positive values. In Table 1, we show the SPI and SDI classes as suggested [52,56].
<table><caption>**Table 1.** Drought classes definition and their corresponding SPI (Standardized Precipitation Index)/Streamflow Drought Index (SDI) values.</caption><thead><tr><td>SPI and SDI Value</td><td>Class</td></tr></thead><tbody><tr><td>Value ≥ 2.0</td><td>Extremely wet</td></tr><tr><td>1.5 ≤ value &lt; 2.0</td><td>Very wet</td></tr><tr><td>1.0 ≤ value &lt; 1.5</td><td>Moderately wet</td></tr><tr><td>-1.0 &lt; value &lt; 1.0</td><td>Normal</td></tr><tr><td>-1.5 &lt; value ≤ -1.0</td><td>Moderately dry</td></tr><tr><td>-2.0 ≤ value &lt; -1.5</td><td>severely dry</td></tr><tr><td>value ≤ -2.0</td><td>Extremely dry</td></tr></tbody></table>

## 2.4. Datasets

### 2.4.1. Projected Climate Variables

In order to force the SWAT model to simulate the future hydrology in BW, we used dynamically downscaled output data from the Global Climate Model CNRM-CM5 (used in the fifth version of the Coupled Model Intercomparison Project (CMIP5)) [61] in this study. Data from the CNRM-CM5 model were successfully used in climate change studies in the area, with low uncertainty risks [33,62,63].

The Bias-Corrected Spatial Disaggregation (BCSD) method [64] and Global Land Data Assimilation System (GLDAS) data at 25km resolution were used to bias correct and spatially downscale CNRM-CM5 simulations. The BCSD method is a two-step procedure that is based on a statistical bias correction and spatial disaggregation to transform bias-corrected future model simulations from coarse-scale to fine-scale [65]. In addition to the reduction of rainfall distribution biases at the GCM scale, the main advantage of this method is its ability to bring up more realistic variability at finer resolution than simpler spatial interpolation techniques [62]. On top of these advantages, the selection of this method was motivated by its previous use under similar climate and physical conditions [66–68].

GLDAS (Global Land Data Assimilation System) is a NASA (National Aeronautics and Space Administration) project to ingest satellite- and ground-based observational data products to generate optimal fields of land surface states and fluxes through advanced land surface modeling and data assimilation techniques [69]. Due to its high grid density over BW (guaranteeing optimal representation of climate variability) and its integration of ground observations, GLDAS data offered good compromise to reduce uncertainty related to future climate projection. Generally, satisfactory bias correction results are obtained when GLDAS data are used to downscale projected climate, including in complex and limited scale regions [70,71]

GLDAS-based baseline data were from January 1985 to December 2005, and for future datasets both the periods of 2030–2050 and 2080–2100 were considered. The first 5 years of both baseline and future periods of modeling process were exclusively dedicated to model warm up and were not considered in different comparisons and analysis.

Two Representative Concentration Pathways (RCPs) were selected: RCP4.5, described as an “intermediate” scenario, and RCP8.5, the worst-case climate change scenario. Emissions are expected to peak around the mid-century under RCP4.5, whereas they are supposed to keep rising throughout the 21st century under RCP8.5 [5]. These two scenarios have been selected for this study as they are believed to uncover climate change dynamics in both future modeling periods (2035–2050 and 2085–2100).

### 2.4.2. SWAT Inputs

Due to the physical diversity of BW and the importance of the performance of SWAT model in predicting reaction of hydrology to droughts in the study watershed, we dedicated large amounts of effort in collecting quality input datasets and calibrating the model.

Digital Elevation Model data of BW were obtained from the Shuttle Radar Topography Mission (SRTM) of NASA’s Space Endeavour flight; the resolution of the DEM data used was 30m. A soil
map of the study watershed has been obtained from the latest version of the Harmonized World Soil Database (HWSD V1.21) [72]. It is a 1km resolution database that was designed by IIASA (International Institute for Applied Systems Analysis) and FAO (Food and Agriculture Organization) in partnership with other institutions: ISRIC (International Soil Reference and Information Centre), ESN (European Soil Bureau Network), and ISC-CAS (Institute of Soil Science-Chinese Academy of Sciences).

Due to the major role played by land cover in hydrologic modeling, we developed a specific 30 m resolution' land cover map for the study area (Figure 2). The analysis and classification of different land cover classes were performed using three high-quality Landsat 8 OLI/TIRS satellite images. Landsat 8 is an American Earth observation satellite launched on 11 February 2013; OLI (Operational Land Imager) and TIRS (Thermal InfraRed Sensor) are the last improved sensors developed for the Landsat 8 satellite. The images were downloaded on 21 March 2019 from the Earth Explorer database of U.S. Geological Survey. ERDAS IMAGINE 2015 software was used for classification and all image processing. Supervised classification was performed in order to come up with the different classes. In addition to expertise on the ecology of the study area, field visits were carried out (during the third week of March 2019) to collect reference point coordinates of some spots where the dominating land cover classes are located. The aim was to develop spectral signatures of the major land cover classes to help the software in assigning the signatures to other equivalent pixels across the watershed map with high confidence. The land cover classification accuracy was assessed using a confusion matrix and the kappa coefficient of agreement. The kappa coefficient of BW's land cover classification was 0.86, which showed an acceptable precision according to Congalton [73].

Observed daily rainfall and temperature (max and min) series for 21 hydrologic years since 1990/1991 were provided by local authorities in charge of climate and hydrology monitoring from 11 different recording stations within the watershed territory. Since only air temperature and rainfall were available as observed climate data for BW, we used the Hargreaves method [74] to calculate the evapotranspiration in this study. The Hargreaves method is widely used and gives accurate results in evapotranspiration estimation over semi-arid conditions (including in Morocco) [75].

For calibration and validation, recorded monthly discharge data from the four different gauging stations (Figure 1) from 1990 to 2005 were confronted to the simulated discharge data at the same four streamgauges; similarly, average annual yields of major crops in many sub-basins (wheat, broad beans, and olive trees) recorded from 1992 to 2005 were compared to simulate annual crops yields.

### 2.4.3. Drought Indices Inputs

In this study, SWAT was used to compute the input variables for both adopted drought indices. The main objective behind using SWAT to generate input data for drought indices calculation is for the benefit of the spatial and physical interpolation capabilities of such a semi-distributed and physically based tool over such a large heterogeneous watershed.

For instance, the simulated monthly precipitation data over BW (from 2035 to 2050, and from 2085 to 2100) under both scenarios (RCP4.5 and RCP8.5) were used to derive 12-month SPI over the simulation periods, whereas the estimated runoff data under both scenarios and along both study periods were used to compute 12-month SDI over BW. Similarly, and for historic SDI-12 and SPI-12 calculations, simulated monthly runoff and rainfall along the baseline period (1990–2005) by SWAT were used.

## 2.5. Methodology

Figure 3 describes the different milestones to forecast future droughts and to evaluate their potential impacts on hydrologic processes in the study watershed.
<figure><img src="image_4.png" /><figcaption>Figure 3. Research methodology flowchart.</figcaption></figure>

Due to the good compromises between a good correlation with recorded climate variables, a higher grid density (compared to observation stations), and the similar statistical intensity distribution to bias-corrected CNRM-CM5 data, we used GLDAS 25 km reanalysis data to run SWAT throughout the baseline period. Then we used the downscaled CNRM-CM5 scenarios (RCP4.5 and RCP8.5) to simulate water balance during the periods of 2035–2050 and 2085–2100. After deriving the required data to compute historical SDI-12 and SPI-12, both hydrological and meteorological drought indices were calculated and tested for correlation with the baseline streamflow data to assess their reliability in predicting future drought. Future drought indices were then estimated for 2035–2050 and 2085–2100 periods under RCP4.5 and RCP8.5 scenarios. To assess the capabilities of both drought indices in capturing droughts, we used the method of simple linear regression with least squares

For hydrological process monitoring under identified dry years, we picked the Total Water Yield (TWYD) parameter, which is simulated by the SWAT model as the sum of surface runoff, lateral flow, and groundwater contribution net of transmission losses. The advantage of this parameter is that it offers the possibility to quantify most of the hydrologic components with high importance for water resource use systems in the study area. To evaluate the spatial extent of droughts over BW, we compared TWYD values (in each sub-basin) under the identified droughts to the average of baseline water yields.

# 3. Results and Discussion

## 3.1. *SWAT Performance*

The SWAT model was initially set up and run throughout January 1985 to December 2005. The first 5 years were allocated to the model warm-up process. Sensitivity analysis of the hydrological module at the four streamgauges revealed that the most influential parameters on BW's hydrology were CN2 (initial SCS runoff curve number for moisture condition II), SOL_AWC (available water capacity of the soil layer), ESCO (soil evaporation compensation factor), and GWQMN (threshold depth of
water return in shallow flow to occur). All the four parameters were adjusted accordingly at the four streamgauges to reach the best reproduction of recorded streamflows.

Taking into consideration that $R^2$ and NSE coefficient start to indicate acceptable goodness-of-fit at the value 0.5, where in the PBIAS value needs to be ±25% [48,76], the monthly streamflow calibration and validation results indicated a satisfactory model performance for hydrological processes over all three rivers in BW (Table 2).

<table><caption>**Table 2.** $R^2$ and Nash–Sutcliffe coefficient (NSE) values during calibration and validation at different gauging stations.</caption><thead><tr><th></th><th></th><th>Aguibate Ezziar</th><th>Ras Fathia</th><th>S.M.Chrif</th><th>Ain Loudah</th></tr></thead><tbody><tr><td rowspan="3">Calibration</td><td>$R^2$</td><td>0.69</td><td>0.73</td><td>0.77</td><td>0.71</td></tr><tr><td>NSE</td><td>0.7</td><td>0.63</td><td>0.87</td><td>0.69</td></tr><tr><td>PBIAS</td><td>12%</td><td>6%</td><td>-10%</td><td>-9%</td></tr><tr><td rowspan="3">Validation</td><td>$R^2$</td><td>0.72</td><td>0.69</td><td>0.52</td><td>0.67</td></tr><tr><td>NSE</td><td>0.62</td><td>0.51</td><td>0.46</td><td>0.5</td></tr><tr><td>PBIAS</td><td>8%</td><td>15%</td><td>-8%</td><td>-11%</td></tr></tbody></table>

The $R^2$ values in calibration and validation periods were all higher than 0.5 for all four gauges, ranging from 0.69 to 0.77 in the calibration period and from 0.52 to 0.72 in the validation period (Table 2). The NSE values were all ≥0.46 for all gauges, ranging from 0.63 to 0.87 in the calibration period and from 0.46 to 0.62 in the validation period. As per PBIAS, all absolute values ranged from 6% to 15%, which represents satisfactory goodness-of-fit (Table 2).

SWAT performance over hydrologic processes in BW is similar to other similar watersheds with equivalent conditions [77–79].

Hydrographs of simulated and observed discharges at the four gauging stations (over the three rivers) showed that the SWAT model succeeded in tracking the recorded flow overall (Figure 4). Most of peak and low flows were fitted well with the observed data along the simulation periods at all the four gauging stations.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Simulated and observed monthly discharges during calibration and validation at the stations (a) Aguibate Ezziar, (b) Ras Fathia, (c) S.M. Chrif, and (d) Ain Loudah.</figcaption></figure>
Similarly to the hydrologic module, results of the calibration and validation of plant growth module were satisfactory also. Four different parameters were influential on yields of the three major crops in BW: HVSTI (harvest index for optimal growing conditions), BLAI (maximum leaf area index), BIO_E (radiation use efficiency), and EXT_COEF (light extinction coefficient). All the four parameters were subject to adjustments in sub basins where each of the three crops are grown while monitoring the adopted goodness-of-fit indicator $R^2$. All $R^2$ values were higher than 0.5 when recorded annual yields of the three major crops were compared to the simulated ones by SWAT over the simulation period (Table 3).

<table><caption>**Table 3.** $R^2$ coefficient values during calibration and validation of plant growth module of SWAT.</caption><thead><tr><th></th><th>Wheat</th><th>Broad Beans</th><th>Olive Trees</th></tr></thead><tbody><tr><td>Calibration</td><td>0.67</td><td>0.64</td><td>0.67</td></tr><tr><td>Validation</td><td>0.86</td><td>0.68</td><td>0.8</td></tr></tbody></table>

The $R^2$ values ranged from 0.67 to 0.86, which confirms the capabilities of the SWAT model to simulate crop performance in BW under the baseline conditions and all related processes (such as soil moisture mechanisms, transpiration, percolation, etc.) [48].

## 3.2. Future Climate

To come up with future trends during 2035–2050 and 2085–2100 periods under both scenarios (RCP4.5 and RCP8.5), we assessed consolidated precipitation and temperature projections via the SWAT model over the whole study watershed, taking into consideration the 1990–2005 period as reference (baseline).

### 3.2.1. Temperature

Figure 5 represents the daily averages of maximum, minimum, and mean temperatures under different scenarios over the Bouregreg watershed, with their standard deviation values shown.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Averages of daily temperatures during simulation periods under baseline, RCP4.5, and RCP8.5, as well as values of year-to-year standard deviation.</figcaption></figure>

A clear increase of temperatures was expected during both study periods under both scenarios compared to baseline averages. The increase values of mean temperatures during the 2035–2050 period were +1.3 °C and +1.44 °C under RCP4.5 and RCP8.5, respectively, and +1.82 °C and 2.73 °C during the 2085–2100 period under RCP4.5 and RCP8.5, respectively.

Regarding the year-to-year variability, the interannual standard deviation values were higher in both periods under both scenarios in comparison to the baseline, indicating more fluctuations of mean temperatures on a year-to-year basis.
### 3.2.2. Precipitation

Figure 6 represents annual averages of total rainfall for baseline and all scenarios during both study periods, along with standard deviation values of interannual rainfall.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Annual rainfall averages under baseline and both Representative Concentration Pathways (RCPs), as well as standard deviation values.</figcaption></figure>

Overall, we expect an equivalent cumulative annual rainfall under RCP4.5 compared to baseline, but with different interannual variability. Standard deviation values were, respectively, 99.1 mm, 152.2 mm, and 161.2 mm for baseline, 2035–2050, and 2085–2100 periods.

Under RCP8.5, it is expected that there will be slightly equivalent cumulative yearly rainfall during the 2035–2050 period and less (327 mm) at the end of the century. The highest standard deviation values were found under this scenario: 186 mm during the 2035–2050 period and 171.3 mm for the 2085–2100 period.

The high standard deviation values under both scenarios (in comparison to the baseline) suggest high interannual rainfall variability, which may lead to a succession of dry (drought) and moist (flooding) periods with higher frequency.

This result is confirmed by the comparison between the evolution of annual TWYDs under baseline and future scenarios (Figure 7). High fluctuation can be clearly observed on TWYD graphics of both RCPs.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Annual distribution of mean precipitation averages under both RCPs during (a) 2035–2050 and (b) 2085–2100 periods in comparison to baseline (1990–2005).</figcaption></figure>

### 3.3. Droughts Characteristics

#### 3.3.1. Historic Drought Propagation

The historic time series of SPI-12 and SDI-12 and the historic annual TWYD in BW (Figure 8) show that both indices followed an identical trend along the baseline period. On the basis of these data, BW
experienced two severely dry years. In these specific hydrologic years, the lowest average annual TWYDs in BW have been recorded: 79 mm in 1992–1993 and 65 mm in 1999–2000.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Annual distribution of historical Total Water Yield (TWYD), SPI-12, and SDI-12.</figcaption></figure>

In 1992–1993, SPI-12 and SDI-12 values were –1.86 and –1.92, respectively, whereas in 1999/2000, their values were –1.96 and –2, respectively. According to observation authorities in charge of tracking extreme weather events in Morocco, the study region underwent severe drought in both years (similarly to most of other regions in Morocco). According to El Khatri and El Hairech [29] and Bouignane [80], Morocco was meant to manage a clear deficit in water resources during both years and to import more cereals to cope with internal demand.

To confirm the ability of both indices to capture droughts, we assessed a simple regression test between recorded historic flow (TWYD) and the computed drought indices. Figure 9 represents the correlation between SPI-12 and annual historic TWYD, and similarly the SDI-12 with historic TWYD.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Scatter plot for correlation tests between annual historical TWYD and historical (a) SPI-12 and (b) SDI-12.</figcaption></figure>

Correlation values for SPI-12 and SDI-12 when compared to historic annual TWYDs were 0.566 and 0.522, respectively, indicating that there was a (moderate) correlation between the both drought indices and the evolution of TWYD on an annual basis. As both indices were equivalently correlated to the selected hydrological indicator (TWYD) for this study, we considered both SPI-12 and SDI-12 to analyze future droughts under both scenarios (RCP4.5 and RCP8.5) in BW.

### 3.3.2. Future Droughts

Figure 10 shows future time series of SPI-12 and SDI-12 for 2035–2050 and 2085–2100 periods under both scenarios.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Annual distribution of SPI-12 and SDI-12 (a) under RCP4.5, (b) under RCP8.5 for the period 2035–2050, (c) under RCP4.5, and (d) under RCP8.5 within the window 2085–2100.</figcaption></figure>

Future drought simulation over BW revealed that the watershed will experience several droughts with different severity levels under both climate change scenarios. Most of dry years were detected by both indices, except in the moderately dry years 2044–2045 and 2048–2049 (Figure 10a) that were identified only by SDI-12 = -1.01 (while SPI-12 = -0.73) and the moderately dry year 2042–2043 (Figure 10b) that has been described by SPI-12 = -1.27 (while SDI-12 = -0.89).

Figure 11 describes the distribution of different severity classes in comparison under each scenario including the baseline.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Distribution of dry years by severity classes according to (a) SPI-12 and (b) SDI-12.</figcaption></figure>

On the basis of the SPI-12 index, and similarly to the baseline situation, there will be two dry years during 2035–2050 (under each RCPs). The difference to baseline situation is that one year should be moderately dry and the second one should be severely dry. During the last 15 years of the present century, it has been predicted that under a realistic scenario, BW should undergo two moderately dry years and one extremely dry year (2098–2099), whereas under the RCP8.5 scenario, three moderately dry years and one severely dry year were simulated.

According to SDI-12 index, and during the 2035–2050 period, BW will experience two moderately dry years and one extremely dry year (2049–2050) under RCP4.5, and one moderately dry year and
severely dry year under RCP8.5. During the 2085–2100 window, it was predicted that under RCP4.5, there will be one moderately dry year and two severely dry years. However, under the RCP8.5, there will be four severely dry years. In overall, around 60% of dry years were predicted for the last 15 years of the century.

All these results confirm the findings of Gudmundsson and Senvirane [26], as they reported that enhanced greenhouse forcing has contributed to increased drying in the Mediterranean region (including Northern Africa), and that this trend will keep increasing over the century, especially under high levels of global warming. In the technical summary of the fifth Assessment Report of the IPCC, many experts relied on various indicators (such as soil moisture drying that is of key importance to many hydrologic and agronomic processes) to conclude that the Mediterranean region (in addition to southwestern USA and southwestern Africa) should experience several droughts (and extreme heatwaves) with more frequency toward the end of the 21st century [4].

## 3.4. Droughts Impacts

### 3.4.1. Temporal Distribution

To understand water resource vulnerability in BW to drought, we focused on the evolution of annual TWYD along the study periods (2035–2050 and 2085–2100) with a particular focus on the identified dry years by SPI-12 and SDI-12.

Figure 12 shows future time series of annual TWYD averages along with the average decrease (in comparison to average TWYD of baseline situation) during dry years.

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Simulated annual TWYD under (a) RCP4.5 and (b) RCP8.5 for the period 2035–2050, and under (c) RCP4.5 and (d) under RCP8.5 for the period 2085–2100 compared to baseline TWYD average along with TWYD decrease during the future dry years.</figcaption></figure>
Hydrological modeling in BW watershed under future emission scenarios revealed that water resources in this watershed will show strong inter-annual variability under the previously described year-to-year variability of climate variables.

In the 2035–2050 period, a significant decrease of TWYD value was simulated in the dry years as all the identified dry years recorded very low TWYD values in comparison to baseline. The decrease ranged from –45.6% to –76.6% in comparison to the 1990–2005 period. The largest change under droughts should happen under the RCP8.5.

During 2085–2100 window, a dramatic TWYD decrease is expected too during the driest years. All the identified dry years using drought indices showed very low water yields. Annual TWYD drops during this period ranged from –65% to –95.6% in comparison to average TWYD in the baseline period.

These results are in line with the general forecasts that are formulated for this region and for the Mediterranean area in general. Despite the aforementioned "equivalent" averages of future precipitation to baseline, some years are expected to cumulate excessive water yields with flooding risks (such as 2041–2042 and 2088–2089 under RCP4.5, and 2035–2036 under RCP8.5). On the other hand, there will be several years with an alarming drop of water yield values in BW (such as 2042–2043 under RCP4.5 and 2091–2092 under RCP8.5). The impact is expected to be more severe under the pessimistic scenario (RCP8.5), as the average of TWYD decrease (over both study periods) under this scenario was found to be –82% (vs. –61% under the RCP4.5).

Regarding the possible difference between drought impact on water yield in the two simulation periods, the simulated drought around the end of the century had a significantly higher impact on water (average –81%) than the predicted impact for 2035–2050 period (average –59%). This difference between the two windows in this study conforms to findings of most specialists about the increasing impacts of climate change in the long term, mainly under pessimistic scenarios (no efficient actions being performed to mitigate greenhouse gas emissions).

In other similar climate change hotspots around the world, Koch and Cherie [81] found that streamflow of the Blue Nile River Basin (in Ethiopia) might drop by –61% during the 2090s under climate change impacts, as very dry years will be experienced by the study watershed. Abouabdillah et al. [82] predicted high occurrence rates of droughts with high intensities in a Tunisian watershed in 2070, leading to potential drop of TWYD by –69% in some years.

For further investigations of the contribution of both droughts (meteorological and hydrological) to such effects on hydrological regime in BW, we performed correlation tests between both drought indices (SPI-12 and SDI-12) and future annual TWYD. The monitoring of correlation results (Table 4) showed that both indices were consistent with respect to changes of the considered hydrological indicator. Correlation values ranged from to 0.81 to 0.91, indicating strong abilities of both indices to track annual water yield evolution.

<table><caption>**Table 4.** Values of correlation coefficients between annual TWYDs and the indices.</caption><thead><tr><th rowspan="2"></th><th colspan="2">2035–2050</th><th colspan="2">2085–2100</th></tr><tr><th>RCP4.5</th><th>RCP8.5</th><th>RCP4.5</th><th>RCP8.5</th></tr></thead><tbody><tr><th>SPI-12</th><td>0.81</td><td>0.89</td><td>0.85</td><td>0.81</td></tr><tr><th>SDI-12</th><td>0.78</td><td>0.92</td><td>0.87</td><td>0.86</td></tr></tbody></table>

The previous analysis of future climate in BW revealed that there will be temperature and precipitation increase (overall) in comparison to historic averages. Precipitation becomes runoff and soil moisture and other TWYD components through the interaction with various factors such as temperature, soil, and vegetation [9]. The increasing trend in precipitation has a direct impact on soil moisture and runoff, but its expression in the hydrological regime is generally linked to other factors. The use of both indices (meteorological and hydrological) was (initially) worthy; however, their equivalent correlations with water yields in historic and future projections can be explained by
low importance of underground components of the hydrological system in BW and the high ratio of land cover classes that promote superficial flow components.

Although there was a good correlation between both drought indices and future TWYD, their limitation in accurately considering the interactions between many mechanisms (such as potential evapotranspiration, plant water demand, etc.) is to be taken into consideration. Furthermore, water budget in some situations should be addressed physically, hence the use of the SWAT model and a comprehensive water budget indicator such as TWYD to understand drought impact on water resources. An example of the justification of this "cautious" approach are the inconsistent rates of TWYD reduction in 2089–2099 (RCP4.5) and in 2097–2098 (RCP8.5) with the computed dryness severities during both years. Although 2098–2099 was revealed as dryer than 2097–2098 (Figure 10c,d), the reduction of water yield was less than in 2097–2098 (Figure 12c,d).

### 3.4.2. Spatial Extent

The previous results help to describe the future occurrence of droughts and to evaluate the potential reaction of water resources during drought years in the whole BW. However, further understanding of the study watershed vulnerability at the individual sub-basin level is needed to put in place appropriate drought planning systems.

Four different dry years (one year by scenario and by time period) were selected to investigate the local impacts of drought on water yield by sub-basin during these specific years (in comparison to baseline). All the four dry years were identified jointly by both drought indices with different drought intensity levels. Table 5 shows the selected dry years along with dryness classes.

<table><caption>**Table 5.** The selected dry years for spatial drought analysis.</caption><thead><tr><th></th><th colspan="2">2035–2050</th><th colspan="2">2085–2100</th></tr><tr><th></th><th>RCP4.5</th><th>RCP8.5</th><th>RCP4.5</th><th>RCP8.5</th></tr></thead><tbody><tr><th>SPI-12</th><td>2049–2050</td><td>2045–2046</td><td>2098–2099</td><td>2097–2098</td></tr><tr><th>SDI-12</th><td>2049–2050</td><td>2045–2046</td><td>2098–2099</td><td>2097–2098</td></tr><tr><th></th><th>Extreme</th><th>Severe</th><th>Moderate</th><th></th></tr><tr><th></th><td>≤−2</td><td>−2 to −1.5</td><td>−1.5 to −1</td><td></td></tr></tbody></table>

The spatial distribution of TWYD reduction by sub-basin (in comparison to baseline) revealed an uneven response to droughts, suggesting different levels of vulnerability (Figure 13).

A significant decrease of water yield values was simulated at most of the BW sub-basins for the selected dry years. The decrease ranged from −10% to −33% in 2049–2050 (RCP4.5) and from −20% to −79% in 2045–2046 (RCP8.5). A more dramatic decrease of annual water yields is expected in the simulated drought years within the 2085–2100 period. The change to baseline ranged from −40% to −86% in 2098–2099 (RCP4.5) and from −60% to −89% in 2097–2098 (RCP8.5).

Besides the variability in terms of range reduction from one dry year to another, some sub-basins showed a consistent response to drought overall. In fact, sub-basins 1, 3, 5, 8, 9, and 13 were the least vulnerable to droughts across the four dry years; the other sub-basins recorded the highest decrease of TWYD in comparison to baseline.

In our perspective, the strong variability in sub-basin vulnerability to drought impact could be explained by the interaction of different biophysical factors related to regional differences in terms of climate, land use, topography, and proximity to the Atlantic Ocean.
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Water yield changes (to baseline average) by sub-basin during the dry years: (a) 2049–2050 under RCP4.5, (b) 2045–2046 under RCP8.5, (c) 2098–2099 under RCP4.5, and (d) 2097–2098 under RCP8.5.</figcaption></figure>

In addition to the potential impacts of climate on the land ecosystem functioning, land use has many ways to affect the atmosphere and thereby climate and weather [5]. For example, bare lands with dry soil conditions will lead to reduced evapotranspiration, which will strengthen the heat wave during dry season [6]. Moreover, evaporation exceeds rainfall in oceans, which leads to moisture transport from the ocean to exposed lands through the atmosphere. Such a phenomenon might moderate meteorological drought in the land phase [83]. In BW, the sub-basins 1, 3, 5, 8, and 9 host annual crops, forest, and pasture lands; the elevation from sea level is the lowest in BW, as most of the terrains are plains and are located in front of the ocean. All these elements means that these sub-basins area taking advantage from the ocean effect to mitigate the drought impact. The situation of sub-basin 13 is different; even that it is not directly exposed to the ocean. This sub-basin is the one with the highest altitudes in the study watershed with a sub-humid climate, which makes this specific sub-basin less vulnerable to drought. Deshmukh and Singh [84] and Stewart [85] argue that in such climatic contexts, watersheds with high altitude bands are generally less sensitive to warm climate extremes since their high portions are unaffected by warming. In BW, the sub-basin 13 still benefits from the moderating effect from the ocean from the West and is high enough to "escape" the warm air masses coming from the South and the East.
The vulnerable group of sub-basins (2, 4, 6, 7, 10, 11, 12, 14, 15, 16, and 17) is expected to be less resilient to droughts, making all processes and activities (such as farming, drinking, and sanitation water services) in these areas vulnerable also. These insights might help in the prioritization of any drought risk management and planning strategies with a region-specific approach.

## 4. Conclusions and Summary

Morocco, like many countries in the Mediterranean region, has undergone many drought episodes with clear impacts on water resources and all the related activities. Bouregreg watershed is a representative Mediterranean hydrological entity with a highly strategic "role", diversified development base, and various climate change-related risks (such as droughts). In order to manage future drought risks, reliable and context-adapted drought forecasting methodologies are needed. The 2035–2050 and 2085–2100 projections of climatic variables over BW revealed that the study watershed will experience temperature increase with equivalent precipitation averages to the baseline (overall). The interannual variability is expected to be higher in the future, suggesting the occurrence of future droughts. This has been confirmed by the 12-month SPI and SDI, which identified several dry years with different severity levels. Drought frequency has been shown to be higher in the last 16 years of this century.

Annual water yield averages are expected to be affected during the simulated dry years. BW sub-basins expressed different vulnerability levels to future droughts according to their properties (elevation, proximity to the ocean, etc.). These results support the necessity to take regional and local approaches in planning drought priorities.

Using multiple indices and a semi-distributed hydrological model provided a valid approach to forecast and to evaluate drought conditions in a Mediterranean context, and would contribute to enhancing drought risk management since multiple factors influencing the persistence and recovery from drought (land use, elevation, ocean effect, climatic zone, etc.) could be considered also.

Taking into consideration the complexity of drought phenomena and the various factors interfering in arid and semi-arid contexts, we recommend the exploration of more drought indices and further physical interpolation approaches to identify the best fit in terms of drought monitoring tools. There is also a need to measure the impact of future droughts on crops and other natural components of the Mediterranean ecosystem to provide inclusive and sustainable solutions to potential droughts.

It is evident that addressing drought phenomena is challenging, but the combination of multiple approaches at high resolutions is very promising to provide fundamental insights on the evolution of future droughts; thus, the drought modeling method suggested in this study can be useful for drought risk forecasting and drought mitigation strategies.

**Author Contributions:** Conceptualization, Y.B.; methodology, Y.B., A.A., K.B.; writing—original draft, Y.B., A.A., L.B.; writing—review and editing, A.C., L.H.; validation, A.C., L.B., R.M.; funding acquisition: K.B., R.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** As a part of multi-milestones project supervised by the International Center for Biosaline Agriculture (ICBA) in the UAE, this work was supported by the Office of Science and Technology, Bureau for the Middle East, U.S. Agency for International Development (grant number AID-ME-IO-15-003). The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the ICBA and U.S. Agency for International Development.

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

## References

1. IPCC. *Climate Change 2014: Impacts, Adaptation, and Vulnerability*. Part A: Global and Sectorial Aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, UK, 2014; p. 1132.

2. Gosling, S.N. The likelihood and potential impact of future change in the large-scale climate-earth system on ecosystem services. *Environ. Sci. Policy* **2013**, *27*, S15-S31. [CrossRef]
3. Flandroy, L.; Poutahidis, T.; Berg, G.; Clarke, G.; Dao, M.-C.; Decaestecker, E.; Furman, E.; Haahtela, T.; Massart, S.; Plovier, H.; et al. The impact of human activities and lifestyles on the interlinked microbiota and health of humans and of ecosystems. *Sci. Total Environ.* **2018**, *627*, 1018–1038. [CrossRef] [PubMed]

4. IPCC. *Climate Change 2013: The Physical Science Basis*. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, UK, 2013; p. 1535.

5. IPCC. *IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems*; IPCC: Geneva, Switzerland, 2019.

6. IPCC. *Impacts of 1.5 °C Global Warming on Natural and Human Systems*; United Nations Framework Convention on Climate Change (UNFCCC): New York, NY, USA; World Meteorological Organization (WMO): Geneva, Switzerland; United Nations Environment Programme (UNEP): Nairobi, Kenya, 2018.

7. Mearns, L.O.; Bukovsky, M.; Pryor, S.C.; Magaña, V. Downscaling of Climate Information. In *Climate Modelling: Philosophical and Conceptual Issues*; Lloyd, A.E., Winsberg, E., Eds.; Springer: Cham, Switzerland, 2018; pp. 199–269.

8. Trenberth, K.E.; Dai, A.; van der Schrier, G.; Jones, P.D.; Barichivich, J.; Briffa, K.R.; Sheffield, J. Global warming and changes in drought. *Nat. Clim. Chang.* **2013**, *4*, 17–22. [CrossRef]

9. Sun, F.; Mejia, A.; Zeng, P.; Che, Y. Projecting meteorological, hydrological and agricultural droughts for the Yangtze River basin. *Sci. Total Environ.* **2019**, *696*, 134076. [CrossRef]

10. Zhang, S.; Wu, Y.; Sivakumar, B.; Mu, X.; Zhao, F.; Sun, P.; Sun, Y.; Qiu, L.; Chen, J.; Meng, X.; et al. Climate change-induced drought evolution over the past 50 years in the southern Chinese Loess Plateau. *Environ. Model. Softw.* **2019**, *122*, 104519. [CrossRef]

11. Li, P.; Omani, N.; Chaubey, I.; Wei, X. Evaluation of Drought Implications on Ecosystem Services: Freshwater Provisioning and Food Provisioning in the Upper Mississippi River Basin. *Int. J. Environ. Res. Public Health* **2017**, *14*, 496. [CrossRef]

12. Jun, K.S.; Chung, E.-S.; Sung, J.-Y.; Lee, K.S. Development of spatial water resources vulnerability index considering climate change impacts. *Sci. Total Environ.* **2011**, *409*, 5228–5242. [CrossRef]

13. Guttman, N.B. Comparing the Palmer Drought Index and the Standardized Precipitation Index¹. *JAWRA J. Am. Water Resour. Assoc.* **1998**, *34*, 113–121. [CrossRef]

14. Vicente-Serrano, S.M.; Beguería, S.; López-Moreno, J.I. A Multiscalar Drought Index Sensitive to Global Warming: The Standardized Precipitation Evapotranspiration Index. *J. Clim.* **2010**, *23*, 1696–1718. [CrossRef]

15. García-Herrera, R.; Garrido-Perez, J.M.; Barriopedro, D.; Ordóñez, C.; Vicente-Serrano, S.M.; Nieto, R.; Gimeno, L.; Sorí, R.; Yiou, P. The European 2016/17 Drought. *J. Clim.* **2019**, *32*, 3169–3187. [CrossRef]

16. Kang, H.; Sridhar, V. Combined statistical and spatially distributed hydrological model for evaluating future drought indices in Virginia. *J. Hydrol. Reg. Stud.* **2017**, *12*, 253–272. [CrossRef]

17. Tsakiris, G.; Loukas, A.; Pangalou, D.; Vangelis, H.; Tigkas, D.; Rossi, G.; Cancelliere, A. Drought characterization [Part1. Components of drought planning. 1.3. Methodological component]. In *Drought Management Guidelines Technical Annex*; Iglesias, A., Moneo, M., López-Francos, A., Eds.; CIHEAM/EC MEDA Water: Zaragoza, Spain, 2007; Volume 58, pp. 85–102.

18. Pathak, A.A.; Dodamani, B.M. Comparison of two hydrological drought indices. *Perspect. Sci.* **2016**, *8*, 626–628. [CrossRef]

19. Chattopadhyay, S.; Edwards, D.R.; Yu, Y.; Hamidisepehr, A. An Assessment of Climate Change Impacts on Future Water Availability and Droughts in the Kentucky River Basin. *Environ. Process.* **2017**, *4*, 477–507. [CrossRef]

20. Lakshmi, V.; Piechota, T.; Narayan, U.; Tang, C. Soil moisture as an indicator of weather extremes. *Geophys. Res. Lett.* **2004**, *31*. [CrossRef]

21. Wang, D.; Hejazi, M.; Cai, X.; Valocchi, A.J. Climate change impact on meteorological, agricultural, and hydrological drought in central Illinois. *Water Resour. Res.* **2011**, *47*. [CrossRef]

22. Hoerling, M.; Eischeid, J. Past peak water in the Southwest. *Southwest Hydrol.* **2007**, *6*, 5.

23. Leng, G.; Tang, Q.; Rayburg, S. Climate change impacts on meteorological, agricultural and hydrological droughts in China. *Glob. Planet. Chang.* **2015**, *126*, 23–34. [CrossRef]

24. Liang, X.; Lettenmaier, D.P.; Wood, E.F.; Burges, S.J. A simple hydrologically based model of land surface water and energy fluxes for general circulation models. *J. Geophys. Res. Atmos.* **1994**, *99*, 14415–14428. [CrossRef]
25. Arnold, J.G.; Srinivasan, R.; Muttiah, R.S.; Williams, J.R. Large area hydrologic modeling and assessment-PartI. *J. Am. Water Resour. Assoc.* **1998**, *34*, 73–89. [CrossRef]

26. Gudmundsson, L.; Seneviratne, S.I. Anthropogenic climate change affects meteorological drought risk in Europe. *Environ. Res. Lett.* **2016**, *11*, 044005. [CrossRef]

27. WB. *Adaptation au Changement Climatique et aux Désastres Naturels des Villes d'Afrique du Nord. Phase 1: Evaluation des Risques en Situation Actuelle et À l'Horizon 2030 Pour la Vallée du Bouregreg*; World Bank: Washington, DC, USA, 2011.

28. IRES. *Anticipation et Gestion des Risques D'evenements Climatiques Extremes et de Catastrophes Naturelles au Maroc*; Royal Institute for Strategic Studies (IRES): Hay Riad Rabat, Morocco, 2014.

29. El Khatri, S.; El Hairech, T. *Drought Conditions and Management Strategies in Morocco*; WMO: Geneva, Switzerland, 2014.

30. Balaghi, R.; Badjeck, M.C.; Bakari, D.; De Pauw, E.; De Wit, A.; Defourny, P.; Donato, S.; Gommes, R.; Jlibene, M.; Ravelo, A.C.; et al. Managing Climatic Risks for Enhanced Food Security: Key Information Capabilities. *Procedia Environ. Sci.* **2010**, *1*, 313–323. [CrossRef]

31. Gommes, R.; ElHairech, T.; Rosillon, D.; Balaghi, R.; Kanamaru, H. Les Futurs Possibles. Available online: https://sites.google.com/site/aridoculture/home/climat/changement-climatique/ (accessed on 12 March 2020).

32. Schilling, J.; Hertig, E.; Tramblay, Y.; Scheffran, J. Climate change vulnerability, water resources and social implications in North Africa. *Reg. Environ. Chang.* **2020**, *20*, 15. [CrossRef]

33. Balaghi, R. *Le Changement Climatique dans la Region de Fes-Meknes-MAPM*; Agence pour le Développement Agricole: Rabat, Marocco, 2017.

34. ABHBC. *Le Plan Directeur d'Amenagement Integré des Ressources en eau du Bassin Hydraulique du Bouregreg et de la Chaouia*; Agence du Bassin Hydraulique du Bouregreg et de la Chaouia: Benslimane, Morocco, 2011.

35. Khattabi, A. The Proposal of IWRM in the Bouregreg Basin and How it Fits with the UNESCO/HELP Policy Program. In *Integrated Water Management*; Meire, P., Coenen, M., Lombardo, C., Robba, M., Sacile, R., Eds.; Springer: Dordrecht, The Netherlands, 2008; pp. 179–187.

36. Rahdou, M. Modelisation Hydrologique du Bassin Versant de Bouregreg (S.S) et Evaluation de Taux d'Envasement du Barrage Sidi Mohamed Ben Abdellah (Maroc) par la Methode Bathymetrique. Master's Thesis, Faculty of Sciences and Techniques Fes, Fez, Morocco, 2010.

37. Zellou, B.; Rahali, H. Evaluation de la Pertinence d'Un Modele d'Evolution du Paysage a Eveluer l'Hydraulique des Inondations; Cas d'Etude: La Vallée de Bouregreg. In Proceedings of the 3ème Edition du Colloque International des Utilisateurs du SIG, Oujda, Marocco, 22-23 November 2016.

38. Khomsi, K.; Mahe, G.; Sinan, M.; Snoussi, M. Evolution des Evenements Chauds Rares dans les Bassins Versants de Tensift et Bouregreg (Maroc) et Identification des Types de Temps Synoptiques Associés. In *Reunion Multi-Acteurs, Sur le Bassin du Bouregreg*; CERGéo: Rabat, Marocco, 2013.

39. Belmatrik, S. *Dispositif de Prevision des Risques de Crues*; ABHBC: Benslimane, Morocco, 2010.

40. Goussot, E.; BrouYao, T.; Laouina, A. Dynamique de l'occupation du sol et statitistiques agricoles sur le bassin versant du Bouregreg au Maroc. *Eur. J. Sci. Res.* **2014**, *126*, 191–205.

41. Arnold, J.G.; Moriasi, D.N.; Gassman, P.; Abbaspour, K.C.; White, M.J.; Srinivasan, R.; Santhi, C.; Harmel, R.D.; Van Griensven, A.; Van Liew, M.W.; et al. SWAT: Model use, Calibration, and Validation. *Am. Soc. Agric. Biol. Eng.* **2012**, *55*, 1491–1508.

42. Valdivieso, F.R.O. Modelamiento Semi Distribuido con Informacion Escasa, de Una Cuenca Hidrografica Para su Planificacion Territorial. Ph.D. Thesis, University of Alcala-Spain, Madrid, Spain, 2009.

43. Moriasi, D.N.; Wilson, B.N.; Douglas-Mankin, K.R.; Arnold, J.G.; Gowda, P.H. Hydrologic and water quality models: Use, calibration, and validation. *Trans. Asabe* **2012**, *55*, 1147–1241. [CrossRef]

44. Gassman, P.; Reyes, M.; Green, C.; Arnold, J. The Soil and Water Assessment Tool: Historical development, applications, and future research directions. *Am. Soc. Agric. Biol. Eng.* **2007**, *50*, 1211–1250. [CrossRef]

45. USDA-SCS. *National Engineering Handbook, V.4*; United States Department of Agriculture (USDA): Washington, DC, USA, 1972.

46. Abbaspour, K.C.; Vejdani, M.; Haghighat, S. SWAT CUP calibration and uncertainty programs for SWAT. In Proceedings of the MODSIM 2007 International Congress on Modelling and Simulation Land, Water and Environmental Management: Integrated Systems for Sustainability, Christchurch, New Zealand, 10–13 December 2007; pp. 1603–1609.
47. Krause, P.D.; Boyle, O.; Base, F. Comparison of different efficiency criteria for hydrological model assessment. *Adv. Geosci.* **2005**, *5*, 89–97. [CrossRef]

48. Moriasi, D.N.; Arnold, G.J.; VanLiew, W.M.; Bingner, L.R.; Harmel, D.R.; Veith, L.T. Model Evaluation Guidelines for Systematic Quantification of Accuracy in Watershed Simulations. *Trans. ASABE 2007*, *50*, 885–900. [CrossRef]

49. Mishra, A.K.; Singh, V.P. A review of drought concepts. *J. Hydrol.* **2010**, *391*, 202–216. [CrossRef]

50. Dracup, J.A.; Lee, K.S.; Paulson, E.G. On the definition of droughts. *Water Resour. Res.* **1980**, *16*, 297–302. [CrossRef]

51. WMO; GWP. *Handbook of Drought Indicators and Indices*; Integrated Drought Management Tools and Guidelines; Series 2; WMO: Geneva, Switzerland, 2016.

52. McKee, T.B.; Doesken, N.J.; Kleist, J. The Relationship of Drought Frequency and Duration to Time Scales. In *Proceedings of the 8th Conference on Applied Climatology, Anheim, CA, USA, 17–22 January 1993*; pp. 179–184.

53. Hayes, M.J. Drought Indices. In *Van Nostrand's Scientific Encyclopedia*; John Wiley & Sons: Hoboken, NJ, USA, 2006.

54. Palmer, W.C. Meteorological Drought. In *Research Paper No.45*; US Department of Commerce, Weather Bureau: Washington, DC, USA, 1965.

55. Institute of Hydrology. *Low Flow Studies Report*; Institute of Hydrology: Wallingford, UK, 1980.

56. Nalbantis, I. Evaluation of a Hydrological Drought Index. *Eur. Water* **2008**, *23*, 67–77.

57. Kamali, B.; Houshmand Kouchi, D.; Yang, H.; Abbaspour, K.C. Multilevel Drought Hazard Assessment under Climate Change Scenarios in Semi-Arid Regions—A Case Study of the Karkheh River Basin in Iran. *Water* **2017**, *9*, 241. [CrossRef]

58. Iglesias, L.; Garrote, L.; Cancelliere, A.; Cubillo, F.; Whilhite, D. (Eds.) *Coping with Drought Risk in Agriculture and Water Supply Systems*; Springer: Dordecht, The Netherlands, 2009; Volume 26, p. 320.

59. Heim, R.R. A review of tweentieth-century drought indices in the United States. *Bull. Am. Meteorol. Soc.* **2002**, *83*, 1149–1165. [CrossRef]

60. Svoboda, M.D.; Hayes, M.J.; Wood, D. *Standardized Precipitation Index User Guide*; Organization WMO, Ed.; WMO: Geneva, Switzerland, 2012.

61. Taylor, K.E.; Ronald, J.S.; Gerald, A.M. Overview of CMIP5 and the experiment design. *Bull. Am. Meteorol. Soc.* **2012**, *93*, 485–498. [CrossRef]

62. Brouziyne, Y.; Abouabdillah, A.; Hirich, A.; Zaaboul, R.; Benaabidate, L. Modeling sustainable adaptation strategies toward a climate-smart agriculture in a Mediterranean watershed under projected climate change. *Agric. Syst.* **2018**, *162*, 154–163. [CrossRef]

63. Hadour, A.; Mahé, G.; Meddi, M. Watershed based hydrological evolution under climate change effect: An example from North Western Algeria. *J. Hydrol. Reg. Stud.* **2020**, *28*, 100671. [CrossRef]

64. Wood, A.W.; Leung, L.R.; Srindhar, V.; Lettenmaier, D.P. Hydrologic implications of dynamical and statistical approaches to downscaling climate model outputs. *Clim. Chang.* **2004**, *62*, 189–216. [CrossRef]

65. Nahar, J.; Johnson, F.; Sharma, A. A rank-based approach for correcting systematic biases in spatial disaggregation of coarse-scale climate simulations. *J. Hydrol.* **2017**, *550*, 716–725. [CrossRef]

66. Fu, J.; Niu, J.; Kang, S.; Adeloye, A.J.; Du, T. Crop production in the Hexi Corridor challenged by future climate change. *J. Hydrol.* **2019**, *579*, 124197. [CrossRef]

67. Ullah, A.; Salehnia, N.; Kolsoumi, S.; Ahmad, A.; Khaliq, T. Prediction of effective climate change indicators using statistical downscaling approach and impact assessment on pearl millet (*Pennisetum glaucum* L.) yield through Genetic Algorithm in Punjab, Pakistan. *Ecol. Indic.* **2012**, *90*, 569–576. [CrossRef]

68. Wang, G.; Zhang, Q.; Yu, H.; Shen, Z.; Sun, P. Double increase in precipitation extremes across China in a 1.5 °C/2.0 °C warmer climate. *Sci. Total Environ.* **2020**, *2020*, 140807. [CrossRef]

69. Rodell, M.; Houser, P.; Jambor, U.; Gottschalck, J.; Mitchell, K.; Meng, J. The Global Land Data Assimilation System. *Bull. Am. Meteorol. Soc.* **2004**, *85*, 381–394. [CrossRef]

70. Du, E.; Vittorio, A.D.; Collins, W.D. Evaluation of hydrologic components of community land model and bias identification. *Int. J. Appl. Earth Obs. Geoinf.* **2016**, *48*, 5–16. [CrossRef]

71. Zhang, Q.; Fan, K.; Singh, V.P.; Song, C.; Xu, C.-Y.; Sun, P. Is Himalayan-Tibetan Plateau "dying"? Historical estimations and future trends of surface soil moisture. *Sci. Total Environ.* **2019**, *658*, 374–384. [CrossRef] [PubMed]
72. FAO. *Harmonized World Soil Database*; Version 1.21; FAO: Rome, Italy; IIASA: Luxenburg, Austria; ISRIC: Wageningen, The Netherlands; ISSCAS: Nanjing, China; JRC: Ispra, Italy, 2013.

73. Congalton, R.G. A review of assessing the accuracy of classifications of remotely sensed data. *Remote Sens. Environ.* **1991**, *37*, 35–46. [CrossRef]

74. Hargreaves, G.H.; Samani, Z.A. Reference crop evapotranspiration from temperature. *Appl. Eng. Agric.* **1985**, *1*, 96–99. [CrossRef]

75. Er-Raki, S.; Chehbouni, A.; Khabba, S.; Simonneaux, V.; Jarlan, L.; Ouldbba, A.; Rodriguez, J.C.; Allen, R. Assessment of reference evapotranspiration methods in semi-arid regions: Can weather forecast data be used as alternate of ground meteorological parameters. *J. Arid Environ.* **2010**, *74*, 1587–1596. [CrossRef]

76. Arnold, J.G.; Kiniry, J.R.; Srinivasan, R.; Williams, J.R.; Haney, E.B.; Neitsch, S.L. *Soil and Water Assessment Tool Input/Output File Documentation-Version 2009*. Texas Water Resources Institute Technical Report No.365; Texas A&M University System: College Station, TX, USA, 2011.

77. Brouziyne, Y.; Abouabdillah, A.; Bouabid, R.; Benaabidate, L.; Oueslati, O. SWAT manual calibration and parameters sensetivity analysis in semi-arid watershed in North-western Morocco. *Arab. J. Geosci.* **2017**, *10*, 427. [CrossRef]

78. Chaponnière, A.; Boulet, G.; Chehbouni, A.; Aresmouk, M. Understanding hydrological processes with scarce data in a mountain environment. *Hydrol. Process.* **2008**, *22*, 1908–1921. [CrossRef]

79. Fadil, A.; Rhinane, H.; Kaoukaya, A.; Kharchaf, Y.; Alami Bachir, O. Hydrologic Modeling of Bouregreg Watershed (Morocco) Using GIS and SWAT Model. *J. Geogr. Inf. Syst.* **2011**, *3*, 279–289. [CrossRef]

80. Bouignane, A. *Sechresse Dans le Bassin du Sebou Etat Actuel et Perspectives*; Agence du Bassin Hydraulique du Sebou: Fes, Marocco, 2010.

81. Koch, M.; Cherie, N. SWAT-Modeling of impact of future climate change on the hydrology and the water resources in the upper Nile River basin, Ethiopia. In *Proceedings of the 6th International Conference on Water Resources and Environment Research, Koblenz, Germany, 3–7 June 2013*.

82. Abouabdillah, A.; Oueslati, O.; De Girolamo, A.M.; Lo Porto, A. Modeling the impact of climate change in a Mediterranean catchment (Merguellil-Tunisia). *Fresenius Environ. Bull.* **2010**, *19*, 2334–2347.

83. Trenberth, K.E.; Fasullo, J.T.; Mackaro, J. Atmospheric Moisture Transports from Ocean to Land and Global Energy Flows in Reanalyses. *J. Clim.* **2011**, *24*, 4907–4924. [CrossRef]

84. Deshmukh, A.; Singh, R. Physio-climatic controls on vulnerability of watersheds to climate and land use change across the U.S. *Water Resour. Res.* **2016**, *52*, 8775–8793. [CrossRef]

85. Stewart, I. Connecting physical watershed characteristics to climate sensitivity for California mountain streams. *Clim. Chang.* **2013**, *116*, 133–148. [CrossRef]

<figure><img src="image_15.png" /></figure>