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

: Understanding the spatiotemporal distribution of future droughts is essential for e ﬀ ective 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–2050and 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. Signiﬁcant 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 inﬂuence appear to be more resilient to drought. The combination of drought indices and the semi-distributed model o ﬀ er a comprehensive tool to understand potential future droughts in Bouregreg watershed.


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][2][3][4]. Besides the increase in land and ocean temperatures, the manifestations of global climate change 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][14][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][20][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 1year in 15 normal years in the 1930s, 1940s, 1950s, 1960s, and 1970s, 1one 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 1990sled 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.

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 2 . 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]. 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 500mm 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 3 , 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 3 and secures 250 to 450 million m 3 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]. 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 3 , 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 3 and secures 250 to 450 million m 3 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].

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 homogenous [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 rainfallrunoff relationship. The daily water balance simulated by the SWAT model is based on the following equation [25]: where SWt: the daily final soil water content (mm H2O); SW0: the initial soil water content on day i (mm H2O); t: the time (days); Rday: the amount of precipitation on day i (mm H2O); Qsurf: the amount of surface runoff on day i (mm H2O); Ea: the amount of evapotranspiration on day i (mm H2O); Wseep: the amount of water entering the vadose zone from the soil profile on day i (mm H2O); Qgw: the amount of return flow on day i (mm H2O).
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].

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 homogenous [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]: 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 3 /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 and and 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 2 , 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 2 ) 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.

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, 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: where 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].

SPI and SDI Value Class
Extremely dry

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][67][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,whereasthey 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).

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), ESBN (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. 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.

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.  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.

SWAT Performance
The SWAT model was initially set up and run throughout January 1985 to December 2005. The first 5years 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. 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.

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). 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][78][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.
Water 2020, 12, x FOR PEER REVIEW 11 of 24 Taking into consideration that R 2 and NSE coefficient start to indicate acceptable goodness-offit 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). 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][78][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.   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). 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].

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). 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).

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).  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. 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. 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.

Precipitation
Water 2020, 12, x FOR PEER REVIEW 13 of 24 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.   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.

Precipitation
Water 2020, 12, x FOR PEER REVIEW 13 of 24 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.

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-1900.

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-1900. 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.  Correlation values for SPI-12 and SDI-12 when compared to historic annual TWYDs were0.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. 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.

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-1900. 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.  Correlation values for SPI-12 and SDI-12 when compared to historic annual TWYDs were0.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. Correlation values for SPI-12 and SDI-12 when compared to historic annual TWYDs were0.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.  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.   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. 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.

Future Droughts
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].

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.
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].

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.  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 2090sunder 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 watershedin2070, 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. 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-2099was revealed as dryer than 2097-2098 (Figure 10c,d), the reduction of water yield was less than in 2097-2098 (Figure 12c,d).

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. 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).
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. 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 subbasins 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 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.

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.