Hydrological Modeling and Impact of Climate Change on Water Resources in the Ziz Valley, Central High Atlas, Morocco

The Upper Ziz basin located in the southeast of Morocco, has a total area of 4,351 km 2 . The surface water feeds El Hassan Addakhil dam, which insures water supply for the downstream cities of Errachidia, Rissani, Erfoud and others along the Ziz valley. This study aimed to evaluate the availability of water resources in this basin known by its arid climate and strong climatic changes. Several global hydrological models at different times were used to simulate the discharge at the outlet. The Statistical Downscaling Model (SDSM) method has been used to reduce the average rainfall and the temperature to predict future climate change related to various Representative Concentration Pathway (RCP) scenarios such as RCP 4.5 and RCP 8.5. The results of the hydrologic models are available, with an NSE of 0.8 for the monthly model during calibration and 0.77 at validation. Future precipitation shows an increasing trend in both scenarios. As for future mean temperature, it will recognize great seasonal variability, such as warming winter and spring and cooling summer and autumn. As a result, simulated future discharge will decrease by 26% under RCP 4.5 and by 24% under RCP 8.5 in the near future.


INTRODUCTION
Human activities such as land use change, land use, and dependence on fossil fuels have increased the anthropogenic greenhouse gas (GHG) concentrations (IPCC, 2013; Hiraishi et al., 2014). This has led to global warming and a global energy imbalance (Wentz et al., 2007;Huang and Pike, 2011). An increase of 0.74 °C in global average temperature has been reported in the past 100 years, particularly between 1906 and 2005 (Van Gameren et al., 2014). In the last 50 years, a significant increase of 0.13 °C has been reported every 10 years. Thus, an increase of 1.1 to 6.4 °C in global average surface temperature is projected in the 21st century (IPCC, 2013). The IPCC in its fifth assessment report (IPCC, 2014) revealed an increase in atmospheric temperature in Africa of 2 °C to 3 °C for the period (2046-2065), and 3 °C to 6 °C for the long-term period (2081-2100) compared to the 1986-2005 reference period. Therefore, the effects of global warming will be felt most strongly through increases in extreme air temperatures that will encourage the occurrence of extreme events and heat stress conditions. This global warming can have a major impact on society, particularly on the hydrological cycle, and thus affect water resources. Indeed, climate change and anthropogenic impacts are widely recognized as the two main drivers of streamflow Hydrological Modeling and Impact of Climate Change on Water Resources in the Ziz Valley, Central High Atlas, Morocco change (Piao et al., 2007). Agriculture, tourism, public health, water demands of hydroelectric operations, and the ecosystem may be affected by this climate change (Chu et al., 2010;Huang and Pike, 2011). However, few assessments of the impact and vulnerability of climate change have focused on the Mediterranean region. Limited examples include the socioeconomic prospects of climate change (Giupponi and Shechter, 2003), showing a 2 °C increase in global temperature (Giannakopoulos et al., 2009). The Saharan areas of North Africa are among the region most exposed to the adverse effects of climate change impacts worldwide (IPCC, 2013). These areas classified as climate change hotspots, experience heterogeneity in the spatial and temporal distribution of precipitation and temperature (Diffenbaugh and Giorgi, 2012). These heterogeneities may lead to socio-economic and environmental crises depending on the area.
In Morocco, this issue is of major importance and is the subject of several national debates regarding the future of water resources, especially since the Moroccan economy depends on the agricultural sector as the mainstay. This sector will be the first to suffer from global warming, as it mainly depends on the water from the El Hassan Addakhil dam (El Ouali et al., 2022). Moreover, it not only provides irrigation, but also a drinking water supply for several downstream cities.
This problem of the impact of climate change on water availability has led to think about the future water resources in the Upper Ziz basin. This basin has an arid climate that limits its water resources due successive periods of drought (Riad, 2003). Hence, multiple hydrologic models were calibrated and validated at different time steps to better characterize and quantify these resources. In parallel with the hydrologic modeling, future climate data was downscaled from the global circulation model under different climate change scenarios. The future climate data (precipitation and temperature) was then fed to the best performing model to generate the future runoff from the downscaled future precipitation and temperature.

GEOGRAPHICAL AND GEOLOGICAL CONTEXT
The Upper Ziz basin lies within The Draa-Tafilalet economic region close to the High Atlas chain (Fig. 1). The study area corresponds to the upstream Ziz valley, part of the Central High Atlas. The Upper Ziz collects rainfall water from an area estimated at 4,351 km 2 .
The study area is located in the Central High Atlas, the main segment of the High Atlas chain which extends from the Atlantic Ocean to the West to the Morocco-Algeria borders to the East. This intracontinental chain structured during the Cenozoic is the result of the tectonic inversion of a subsident basin resulting from a Triassic-Jurassic rifting contemporaneous to the break-up of Pangea and the opening of the Central Atlantic. The filling sequences of this basin, are characterized by thick marine Triassic and Jurassic deposits overlain by continental red-beds, attributed to Upper Jurassic-Cretaceous stages. The Lower Jurassic layers consist of massive dolomites and limestone, which gives the Central High Atlas its unique structure corresponding to a succession of narrow anticline ridges and large synclines (e.g., Jbel Bouhmid, Jbel Tillicht) (Combe, 1977 , 2023). This is followed by the well-bedded limestone and more important marl interlayer and the much more diverse Sinemurian layers (Piqué et al., 2006).
From the Lower Toarcian to the end of the Lower Bajocian, the sedimentation changes abruptly and becomes noticeably marlier. The end of this stage is marked by thick limestone layers forming the ledges that may be Upper Bajocian age (Sadki, 1996). The Plio-Quaternary deposits are represented by puddings, gravels, pebbles, chalky and sandstone fluvial-lacustrine elements, sometimes impermeable marl, and finally, generally soil-forming silts (El Kayssi et al., 2019).

Climatology
There are four monitoring stations in the study area, all of which measure rain, whereas only two (Foum Zaabel & Foum Tillicht) measure temperature. However, for other climatic parameters, only one station (Foum Zaabel) is equipped with such instruments. The latter is considered as the reference station for all the other measuring stations. Figure 3 below shows that the Upper Ziz is characterized by two wet periods, during spring and autumn, separated by short dry winter and summer. From the north (the Zaouia Sidi Hamza station) towards the south (the Foum Zaabel station), a slight decreasing trend in precipitations can clearly be seen. Annual precipitation varies greatly from one year to the next other ( Figure 4). The study area went through some periods of drought, such as the late 90s followed by wet years, which indicates an alternation of dry and wet years.

Temperature
The average monthly temperature from the two available stations shows that the highest records were registered in July and August. The coldest months are both January and February (Fig. 5).
Conversely, it is noticeable that an increasing tendency of precipitation takes place from north to south. This is because the study area is located in the mountains and the downstream of the basin is in a Saharan region.

Watersheds characterization
The Upper Ziz basin is drained by the Ziz river and its tributaries such as the Zaouia Sidi  . Therefore, it can be divided into four watersheds using the monitoring stations as outlets ( Figure 6a). This dissection of the study area was made to compare the contribution of each watershed.
The study area, Upper Ziz, includes all previous watersheds. Its total area is about 4355 km 2 , and a perimeter of 447 km. The outlet of the Upper Ziz is marked by the El Hassan Addakhil dam, which was put into service in 1971 (El Ouali, 1992).
The following table summarizes the principal physical characteristics of the Upper Ziz basin and its tributaries. A compactness coefficient, greater than 1 shows that the watersheds in the study area are elongated. The mean altitude for the Foum Tillicht and Mzizel basins is high, because they represent the upstream part of the study area, and the altitude decreases as the lower basins move towards Foum Zaabel and the Upper Ziz. The two watersheds Foum Tillicht and Foum Zaabel are physically close to each other with equivalent rectangles that are of almost the same size. Regarding the longest flow path, Upper Ziz has the longest drainage with a total length of 185 km.

Depth of the annual rainfall
The annual rainfall depth for the study area and its sub-basins was calculated based on the isohyet map generated using the average of 30-year rainfall data ( Figure 6b). It is seen on the map that the Upper Ziz basin is located between 300 and 150 mm isohyets. It is also observable that precipitation decreases from north to south in accordance with the altitude gradient. As a result of the calculation of the annual mean depth of rainfall using the isohyets method (Table 2), one can see that the Foum Tillicht and Mzizel watersheds are the most irrigated ones.

GR4J
GR4J is a well-known and thoroughly tested conceptual model with only four parameters    (Perrin et al., 2003). It was built with a bottom-up approach by increasing the complexity of simpler modeling step-by-step from the previous version GR3J ( Edijatno et al., 1999).

MORDOR
The French MORDOR (Modèle à Réservoirs de Détermination Objective du Ruissellement) is a model that simulates the functioning of a watershed. It was developed by the Water Resources Service for hydrological monitoring and forecasting in 1990 to extend the scope of methods used. It radically schematizes the behavior of the watershed in a system of five reservoirs, each representing a different form of water storage.

IHACRES
The IHACRES model is an Australian rainfall-runoff model, which stands for (Identification of unit Hydrographs and Component flows from rainfall, evaporation, and stream flow data). Jakeman developed the first version at the Centre for Resource and Environmental Studies at the Australian National University (Jakeman et al., 1990). The model operates at a catchment scale and typically operates in a daily time step. The IHACRES version used in this study, is a modified version called IHAC and has six parameters as described in Oudin et al., (2005).

GR2M
The GR2M model is a global two-parameter rainfall-runoff model. It is a monthly time-step model. It has two reservoirs, a production (or ground reservoir) and a routing reservoir, where adjustment and interception are given differently on the inputs.
The model uses the monthly rainfall and evapotranspiration as inputs and outputs the flow rate. Its development was initiated at Cemagref in the late 1980s with application objectives in the field of water resources and low water levels.

Data
The climatic data used came from the Guir-Ziz-Rheris Hydraulic Basins Agency. They provided two datasets: onein daily time step and the other in monthly time step. The monthly time series for the four stations ranges from 1990/1991 to 2019/2020. In turn, the daily time series only covers the period from 2000/2001 to 2019/2020.
For the climate change scenarios, the outputs of the second-generation Canadian Earth System Model (CanESM2) were downscaled. These Figure 7. Flowchart of the adopted methodology outputs can be downloaded from (http://climatemodelling.canada.ca/climatemodeldata/data. shtml). The downloaded file contains five subfolders; each contains 26 predictors. The NCEP folder holds the normalized predictors used to calibrate the SDSM; the other four folders represent the historical and projected scenarios of the CanESM2 model.
The methodology used in this study was carefully designed to ensure reliable and valid results. The following flowchart (Fig. 7) provides a summary of the steps taken to collect and analyze data.

Data processing
The raw data received contained some missing values estimated using the Inverse Distance Weight (IDW) method. This method has proven to estimate good results when it comes to climatic data influenced by the distance factor (Campozano et al., 2014).
The models require as inputs precipitation (P) and potential evapotranspiration (PET) averaged over the watershed at the relevant time step. In turn, P values were estimated by calculating the mean of the considered stations, while PE values were calculated by taking the mean temperature using Odin's formula (Oudin, 2004;Oudin et al., 2005) for the two stations that have temperature measurements. This formula has proven to be the most efficient for hydrological modeling, using Oudin et al., (2005) after testing on 27 formulas. It only needs the station's mean temperature and the latitude as inputs. The formula is defined by the following equation: where: λ -the latent heat of vaporization of water (2.25MJ.KJ -1 ); ρ -the density of water (1000 Kg.m -3 ); Re -the terrestrial radiation (MJ.m -2 .j -1 ); T m -the mean temperature of the air in °C for a given Julian day. The observed discharge is used for calibration and evaluation of the model.
Two statistical tests were applied for the description of future flow time series. The first one is Pettit's test (Pettit, 1979), a non-parametric test that not need any assumption about the distribution of the data. Checks for an unknown break moment in the series from a formulation derived from the Mann-Whitney test. This test is more sensitive to any average variation and offers an estimation of the date of breakout date if the H0 hypothesis of homogeneity of the series is rejected. This test is based on the calculation of the U tT variable defined by (Eq. 2): where: The second one is the Man-Kendall test (Mann, 1945;Kendall, 1975), which is a nonparametric test. It is used to detect an upward or downward linear trend in a time series. In this test, the null hypothesis H0 (absence of trend) is accepted if the p-value is greater than the alpha significance level. Mann-Kendall's statistical coefficient defines the trend direction. It is upward if UMK > 0 and downward if UMK < 0.

Data criticizing
The data may contain possible errors. These uncertainties are mainly due to the fact that stream flows are measured in an unstable rocky bed based on a scale. The bed is changed from time to time as strong events occur. Thus, despite the monthly gauging frequency, rating curves are difficult to compute and may soon lose their validity, especially for the Foum Tillicht station, which is said to have more than 44 rating curves (AHBGZR, 2018). Furthermore, the location of the measuring point and the low spatial resolution can also be a possible source of error because precipitation in arid and semiarid zones is largely due to convective cloud mechanisms that typically produce short duration, relatively high intensity and limited-areal extent (Unland et al., 1996). However, the major source of errors is undoubtedly the protocol for collecting such data in terms of quality and not quantity.

Hydrological modeling
Daily hydrologic models were applied only to the Foum Tillicht watershed because testing the average daily precipitation over the whole study area always led to errors.

Foum Tillicht watershed
The inputs of the models for this watershed were calculated from two stations ZSH and Foum Tillicht. P was obtained by averaging the two stations while PET was estimated using the mean temperature from ZSH only. The data for three years from 01/09/2006 to 31/08/2009 are missing in the daily precipitation time series for the Foum Tillicht station. As a result, it was necessary to exclude this period from time series. Therefore, the dates from 01/09/2000 to 31/08/2006 are taken as the calibration period for all the daily models. The verification period started from 01/09/2009 to 31/08/2006.
In order not to waste a year of data as a warm-up period, the initial conditions were defined by running the model for one year and then feeding it the with the end-state conditions. This method has been proven to give good results for the warm-up period without wasting any data (Hajhouji et al., 2018).

GR4J
For the calibration period, the GR4J performs well. Its NSE value is 0.67 and KGE value is 0.6. However, the model seems to underestimate the discharge as evidenced by the bias term in the KGE formula with a value of 0.67 (Fig. 8).
For the validation period from September 2009 to August 2016, the performance of the model is considered satisfying. The criteria values are 0.45   (Fig. 9). These results are similar to those found in the literature for basins with the same characteristics as the one adopted in this study (Simonneaux et

IHAC
The IHAC model performed satisfactorily during the calibration period (Fig. 10). Its NSE value is 0.43 and the KGE value is 0.46. During the validation period, the performance of model did not change much (Fig. 11). The benchmark values are 0.4 and 0.43 for both NSE and KGE respectively.
The third and final model applied to the Foum Tillicht watershed is MORDOR7. This model showed a good performance during the calibration period (Fig. 12) with values of 0.61 and 0.63 for both the NSE and the KGE, respectively. During the validation period, the model gave a satisfying performance (Fig. 13). The criteria values are 0.47 for the NSE and 0.51 for the KGE.

Monthly model "GR2M"
The GR2M model showed very good performance. The overlaying of the observed and simulated graphs shows a very good fit (Fig. 14). The benchmarks also show good values, 0.7 for the NSE and 0.74 for the KGE in the calibration period.
During the validation period, the model dropped to NSE to 0.5 and KGE 0.65 (Fig. 15). This deterioration in performance quality may result from diversion dams because their diverted flows are not taken into account by the model causing an overestimation of the outflows at the outlet.

Foum Zaabel watershed
Due to its huge area, the application of the daily models has always led to errors and the NSE has been negative. Therefore, the authors decided to adopt the monthly model for this watershed.
The period of 09/2000 to 08/2007 was taken as a calibration period for the GR2M model. In order not to waste a year's data as a warm-up period, the initial conditions are defined by running the model with ending state conditions after it has been run for one year. The result of the simulation   showed a perfect fit between the simulated and observed graphs (Fig. 16). This perfect overlap of the graphs is also confirmed by the very good values of the numerical criteria with an NSE of 0.80 and KGE of 0.85.
The validation period has been changed from 09/2007 to 08/2014. The results show a very satisfactory fit between the observed and the simulated graphs (Fig. 17). This good fit is confirmed by the good NSE of 77% and a KGE of 75%. These results demonstrate the robustness of the GR2M model in adapting to different climate conditions. This quality is what makes it capable of giving reasonable and more realistic simulations for future projections.

Climate change
The outputs of the CanESM2 model were downscaled using statistical downscaling to assess climate change under different scenarios in the study area. First, the model was calibrated based on the observed data by scanning the variables to determine which accounts for the weather in the study area (Table 3). Second, the weather generator option was used to generate the weather for the same period of data observed in the model. After the model was finally calibrated, a scenario generator was used. This option enables to downscale the different scenarios simulated by the CanESM2. For this study, RCP 4.5 and RCP 8.5 used in IPCC 5th assessment report, were chosen.

Model calibration
The visual presentation shows that the simulated mean monthly precipitation can reproduce the observed data but not the rainfall characteristics. However, the results can be considered satisfactory for simulating the future climate, as other studies raise the same issues but still progress by downscaling future data (Ouhamdouch and Bahir, 2017; El Hafyani et al., 2023).   When the graphs representing three periods for four stations are examined, an increase in precipitations is observed for four stations within the scope of the RCP 4.5. The ZSH station shows an increase of 18% and 17.2% for the near and far future periods respectively, and the increase is mainly noticed in the spring season. The increase rates for the Foum Tillicht station are 22.4% and 17.4% over the historical period. The Foum Zaabel station shows a steady increase throughout the year, with rates of 25% for both the near and far future periods. Finally, the Mzizel station has the lowest increase rates for the next two periods with 8.5% and 7%, respectively (Fig. 19).
Examination of the results obtained for the projection under RCP 8.5 for the four stations shows an increase in precipitations. The ZSH station shows an increase of 19% and 16.2% for the near and far future periods, respectively, and the increase is mainly seen in the spring season. The increase rates for the Foum Tillicht station are 24.1% and 19.4% compared to the historical period, and contrary to the historical period, an increase is noticed in the winter and spring seasons, where a tendency to turn into a rainy season in winter is seen. The Foum Zaabel station shows a steady increase throughout the year, with rates of 27.7% and 24.4 for the near and far periods. Finally, the Mzizel station has the lowest increase rates for the next two periods, with 9.1% and 7%, respectively (Fig. 20). Figure 21 represents the variation of the observed mean monthly temperature and the simulated temperature by the SDSM. The simulated mean temperature is in good agreement with the observed temperature. This perfect fit shows that the model is well-calibrated, providing greater reliability to future projections.

Zaouia Sidi Hamza (ZSH)
In both scenarios, the mean temperature in ZSH will decrease during fall, stabilize from December to January, and increase from February to July (Fig. 22).

Foum Zaabel
The mean temperature at the Foum Zaabel station will increase from January to June and decrease from August to December in both scenarios (Fig. 23).

Impact on future runoff
To assess the future impact of climate change on runoff in the study area, the mean precipitation the study area would receive under two scenarios was calculated. The mean temperature from the two stations was used to calculate the potential evapotranspiration using Oudin's formula. The potential evapotranspiration was then averaged for each RCP scenario. Average precipitation and potential evapotranspiration were then fed into the calibrated model to obtain future runoff.

Monthly runoff
In both scenarios, the mean monthly runoff at the Foum Zaabel station shows a great decrease compared to the historical period, consistent with other studies in the literature (Driouech et   Under RCP 4.5, the period 2040/2070, which represents the near future, shows a mean decrease of 26% compared to the historical period. This rate jumps up to 29% for the far future period. Regarding RCP 8.5, the near future period will experience a mean decrease of 24%, which will increase to 28% for the far future 2070/2100 period (Fig. 24).

Annual runoff
To assess the annual variation of the runoff, monthly means were collected and then multiplied by the surface of the catchment to obtain variations in volume. The graphical representation shows a clear downward trend for both scenarios confirmed by the Man-Kendal test ( Table 4).
The Pettit test was also applied to identify change points, if any. The annual runoff under RCP 4.5 will recognize a change point in 2044/2045. The mean will drop from 1318 Mm 3

CONCLUSIONS
The Upper Ziz basin is located in the southeast of Morocco. Its total area is 4351 km 2 . This area is drained by the Upper Ziz River and its tributaries. The study area is dominated by the Jurassic outcrops. Upper Ziz is characterized by two wet seasons, spring and fall, separated by a short and dry winter, with high in temperature recorded in July and August and low values in January. The study area, which is the source of the Ziz basin, is the most irrigated part with an average of 225.95 mm per year. This relative water wealth is one of the reasons behind the construction of El Hassan Addakhil Dam. The latter is constructed to store the floods of the Upper Ziz and provided water supply to downstream cities.
Application of hydrologic models, especially diurnal ones, has revealed that the flow measurements relating to rating curves are questionable, as deposits vary from one event to the next. However, the monthly measurements appear reliable, as evidenced by the good performance of the monthly model in reproducing the observed discharge, particularly for the Foum Zaabel station, considered as the reference station. The monthly model was the best performing alternative with an NSE of 80% and 77% for the calibration and validation period, respectively.
The precipitation projection under RCPs 4.5 and 8.5 shows an increasing tendency for all the stations. The mean increase rates under RCP 4.5 are as follows; 7% for Mzizel, 18% for ZSH, 19% for Foum Tillicht, and 25% for Foum Zaabel compared to the annual total for the same length historical period. Mean increase rates under RCP 8.5 are 8% for Mzizel, 17% for ZSH, 20% for Tillicht, and 26% for the Foum Zaabel station. These increases are mainly in the spring season and will likely be in the form of storms given the aridity of the study area.
Future temperature predictions indicate that under the RCP 4.5 the study area will become warmer up to 2.2 °C from December to July and colder by -2.27 °C from August to November compared to the historical period 2000/2020. As for  Figure 25. Annual evolution of the future discharge at the Foum Zaabel station RCP 8.5, the temperature will increase by 2.1 °C from December to July and decrease by -2.33 °C from August to November, similarly to the RCP 4.5. Future runoff was calculated by feeding future projections of both precipitation and PET from future temperature. The result shows that the mean runoff will decrease by about 26%under RCP 4.5 and 24% under RCP 8.5 in the near future. As for the far future, the decrease will be 29% for RCP 4.5 and 28% for RCP 8.5. This decrease is confirmed by the Man-Kendall test, which indicates a downward trend. Similarly, the Pettit test marks a turning point in the annual runoff time series scenarios. The change point for the RCP 4.5 scenario is expected in 2044/2045, while for the RCP 8.5 scenario the change point is expected to be at 2064/2065. This decrease in runoff will reduce the amount of water stored in the dam, and increase the water stress that the region is already experiencing.
In light of these findings, it is strongly recommended to consider good policies along with integrated water resources management strategies to reduce and mitigate the impact of climate change on water resources. For example, investing in rainwater harvesting techniques can be a great solution, given that PET represents the greatest loss of the hydrological budget. Another aspect to consider is to improve irrigation methods and water uses as well as friendly agriculture by changing water-consuming crops to adapted alternatives convenient for the arid climate of the study area and other similar regions.