Climate change impacts on surface water resources in the Rheraya catchment (High Atlas, Morocco)

ABSTRACT

This study aimed to quantify possible climate change impacts on runoff for the Rheraya catchment (225 km²) located in the High Atlas Mountains of Morocco, south of Marrakech city. Two monthly water balance models, including a snow module, were considered to reproduce the monthly surface runoff for the period 1989‒2009. Additionally, an ensemble of five regional climate models from the Med-CORDEX initiative was considered to evaluate future changes in precipitation and temperature, according to the two emissions scenarios RCP4.5 and RCP8.5. The future projections for the period 2049‒2065 under the two scenarios indicate higher temperatures (+1.4°C to +2.6°C) and a decrease in total precipitation (−22% to −31%). The hydrological projections under these climate scenarios indicate a significant decrease in surface runoff (−19% to −63%, depending on the scenario and hydrological model) mainly caused by a significant decline in snow amounts, related to reduced precipitation and increased temperature. Changes in potential evapotranspiration were not considered here, since its estimation over long periods remains a challenge in such data-sparse mountainous catchments. Further work is required to compare the results obtained with different downscaling methods and different hydrological model structures, to better reproduce the hydro-climatic behaviour of the catchment.

EDITOR M.C. Acreman

ASSOCIATE EDITOR R. Hirsch

KEYWORDS:

• water resources
• climate changes
• Med-CORDEX
• climate scenarios

1 Introduction

The mountainous chain of the High Atlas in Morocco constitutes a “water tower” providing most of the water supply for the neighbouring arid plains (Chehbouni et al. 2008, Jarlan et al. 2015). In spite of the importance of snow in the regional water balance, to our knowledge, no studies have been devoted so far to the hydrological modelling of snow-driven catchments in Morocco (or in North Africa) under changing climate conditions, while several such studies have been performed for mountainous ranges in Europe (Dobler et al. 2012, Gobiet et al. 2014, Ruelland et al. 2015, Szczypta et al. 2015), North America, (Furey et al. 2012), South America (Hublart et al. 2013) and Asia (White et al. 2014). However, this type of modelling is necessary to evaluate climate change impacts on water resources for this mountainous area in order to develop adaptation strategies.

Future climate projections indicate a possible temperature increase of 2–3°C for the Mediterranean basin by 2050, which might be associated with a decrease in precipitation of about 20–30% (Giorgi and Lionello 2008, Schilling et al. 2012). Most projections for North Africa agree on a future decrease of precipitation associated with a temperature increase, notably over Morocco (Schilling et al. 2012). These changes in the Mediterranean mountains may be associated with earlier declines in high flows from snowmelt in spring and increasing water resource uncertainty because of the reduced influence of snow accumulation and snowmelt processes (Garcia-Ruiz et al. 2011). The changes projected at the regional scale need to be analysed at the basin scale, which is relevant for water management. A common approach for hydrological impact studies of climate change at the basin scale is to use hydrological models with scenarios provided by climate model outputs downscaled to the catchment of interest (e.g. Etchevers et al. 2002, Chiew et al. 2009, Tramblay et al. 2013). Various downscaling methods exist to bridge the gap between the coarse outputs of global climate models (GCMs) and catchments (Teutschbein and Seibert 2013). These methods range from statistical approaches, requiring local calibration (Piras et al. 2015), to dynamical approaches, such as regional climate models (RCMs), that do not require such local calibration, thus providing added value in data-sparse catchments (Tramblay et al. 2013). For the Mediterranean region, a new generation of RCM simulations provided by the Med-CORDEX initiative (Ruti et al. 2016) is aimed at increasing the reliability of information on regional climate and understanding the processes responsible for the variability and evolution of the Mediterranean climate.

The hydrological impacts of climate change have been assessed across Mediterranean catchments using a variety of models, ranging from conceptual (Hreiche et al. 2007, Majone et al. 2012, Milano et al. 2013, Tramblay et al. 2013, Fiseha et al. 2014, Lespinas et al. 2014, Ruelland et al. 2015) to physically-based ones (Quintana Seguí et al. 2010, Senatore et al. 2011, Fatichi et al. 2014, Smiatek et al. 2014). However, the choice of hydrological model is strongly dependent on data availability and could also influence future projections (Bosshard et al. 2013). In order to provide future projections, hydrological models must be able to reproduce adequately various climatic conditions, which implies using data to calibrate and validate the models over long time periods (Vaze et al. 2010). As noted by Blöschl and Montanari (2010), a long observation period is required to better understand the relationships between climate and hydrology and the long-term variations in the data. In developing countries with data-sparse networks, this makes the use of complex and physically-based process models difficult, partly due to the lack of long-term measurements necessary to run such models and to validate intermediate processes taken into account by complex mechanistic models. That is why a parsimonious conceptual model has been favoured for this initial study as regards the available data. However, since conceptual models require calibration, which can also compensate for systematic errors in the data, they may fail under conditions different from those used for calibration (Klemeš 1990, Refsgaard et al. 2014, Thirel et al. 2015).

This study builds upon existing studies for the High Atlas catchments in Morocco, providing datasets available for the present work and a detailed description of the catchment function. Previous studies have focused on (a) the monitoring of snow by remote sensing (Schulz and De Jong 2004, Chaponnière et al. 2005, Boudhar et al. 2010, Marchane et al. 2015), (b) the modelling of snow cover and runoff using hydrological models such as GR4J and SWAT (Chaponnière et al. 2008, Simonneaux et al. 2008, Boudhar et al. 2009), and (c) the evaluation of the water and energy balance of the snowpack at the station scale (Schulz and De Jong 2004, Boudhar et al. 2016). These studies highlighted the significant spatiotemporal variability of snow cover and its sensitivity to climate variability, and the importance of sublimation in the ablation of snowpack at high altitudes. These characteristics are typical of the semi-arid climate encountered in the High Atlas, where elevations reach about 4000 m a.s.l. Indeed, the rugged topography and the steep slopes of the basins result in strong spatiotemporal variability of the weather conditions. The aridity and intense solar radiation (with low relative humidity of the air) are notably responsible for the non-negligible contribution of sublimation to snow ablation by contrast with temperate basins, of between 20% and 44% in the Atlas Mountains, according to Schulz and De Jong (2004) and Boudhar et al. (2016).

The objective of this paper is to quantify possible climate change impacts on runoff in the Rheraya catchment. To this end, two hydrological models, each with a snow module, are first run to reproduce runoff at the catchment outlet. Then, an ensemble of regional climate models is considered in order to produce precipitation and temperature scenarios that are then incorporated in the hydrological simulations via a perturbation approach. In Section 2, the study area and the available datasets are presented; in Section 3, we describe the hydrological models, together with the methods used to create future scenarios; Section 4 shows the results; and in Section 5 the results are discussed.

2 Study area and datasets

2.1 The Rheraya catchment and in situ measurements

The catchment area of Rheraya (225 km²) is located in the High Atlas Mountains, 40 km south of the city of Marrakech, between latitudes 30°10ʹN and 30°20ʹN, and longitudes 7°40ʹW and 8°W. The catchment presents a large altitudinal gradient with elevations ranging from 1027 m to 4167 m at the top of Mount Toubkal, the highest point of North Africa (Fig. 1). The characteristics of the Rheraya catchment are presented in Table 1. The climate is semi-arid, precipitation is very irregular in space and time (Chaponnière et al. 2008), and the annual peak runoff occurs in spring due to snowmelt (Fig. 2). The relative contribution of snowmelt to annual streamflow is between 30% and 48%, depending on the year (Boudhar et al. 2009).

Table 1. Catchment characteristics.

Catchment	Area (km^2)	Perimeter (km)	Average altitude (m)	Altitude min (m)	Altitude max (m)	Average slope (%)
Rheraya	225	78	2160	1027	4167	19.1

Figure 1. Location of the Rheraya catchment.

Figure 2. Mean inter-annual (a) daily precipitation, (b) temperature and (c) runoff for the Rheraya catchment.

Despite the low density of meteorological stations in the Rheraya catchment, it is the most densely instrumented sub-basin of the High Atlas. The data considered in the present study have been collected and processed in the work of Simonneaux et al. (2008), Chaponnière et al. (2008) and Boudhar et al. (2009). The station network in the Rheraya catchment was set up in the framework of the SUDMED (Chehbouni et al. 2008) and LMI TREMA (Khabba et al. 2013, Jarlan et al. 2015) programmes. For this study we used daily data for precipitation recorded at the Tahanaout (1027 m) and Oukaimeden (2650 m) stations for the period 1989‒2009. The Oukaimeden station also has measurements of snow heights carried out by the staff of the Club Alpin Français (CAF). The precipitation series over the Rheraya catchment, including both liquid and solid measurements, were computed by Simonneaux et al. (2008). In addition, average, maximum and minimum daily temperature data were collected at the Aremd, Imskerbour, Neltner, Oukaimeden and Tachedirt stations for the period 2007–2013. For the Aremd station, other available variables include relative humidity, wind and solar radiation for the same period (2007–2013). Figure 2 shows the seasonal variability of precipitation, temperature and runoff. A decrease in runoff is observed at Tahanaout, with a shift towards drier conditions around the year 1997 (Simonneaux et al. 2008). No clear trend could be found in precipitation; however, a slight decrease in snow amount at the Oukaimeden station might be a cause of the declining runoff. No significant trends were detected in the temperature records.

The runoff data used in this work were collected at the Tahanaout station for the period 1989‒2009, using a stage–discharge relationship designed by the water management agency of the Tensift basin (ABHT). Since the runoff is measured based on a scale in an unstable rocky bed, the section is sensitive to change as strong events occur, which requires frequent updating of the rating curves (Simonneaux et al. 2008). The hydrological simulations were thus only evaluated at the monthly time step in order to limit the errors that can be expected in the daily discharge values.

2.2 Climate model outputs from the Med-CORDEX initiative

In this study, five RCM simulations with a horizontal resolution of 50 km were used. These selected simulations belong to core simulations of the Med-CORDEX initiative (http://www.medcordex.eu), which aims to produce for the Mediterranean region a set of simulations at different resolutions taking into account the specificity of the Mediterranean climate (Table 2). Different GCM/RCM combinations were considered to assess the uncertainty in future projections.

Table 2. Med-CORDEX simulations available at 50 km resolution (MED 44 grid). GUF: Goethe University Frankfurt; IPSL:Institut Pierre Simon Laplace; CNRM: Centre National de Recherches Météorologiques; ICTP: International Centre for Theoretical Physics; CMCC: Centro Euro-Mediterraneo sui Cambiamenti Climatici.

Institute	RCM	GCM
GUF	CCLM	MPI-ESM
IPSL	LMDZ	IPSL
CNRM	ALADIN	CNRM
ICTP	RegCM4	HAD
CMCC	CCLM	CMCC

For each model, the following climate simulations were selected:

1. Evaluation runs (EVAL): simulations for the control period 1979‒2005, with the RCMs forced by ERA-Interim reanalysis (Dee et al. 2011). This allows the RCM outputs to be evaluated in comparison to the observed data.

2. Historical runs (HIST): simulations for the control period (1950‒2005, depending on the model) with RCMs forced by the different GCMs indicated in Table 2.

3. Scenario runs (RCP4.5 /RCP8.5): simulations for future projections (2005‒2100) with RCMs forced by different GCMs according to the emissions scenarios RCP4.5 and RCP8.5.

3 Method

3.1 Monthly water balance models

Two conceptual monthly water balance models were considered, the Génie Rural à 2 paramètres model (GR2M) proposed by Mouelhi et al. (2006) and the monthly water balance model (MWBM) of McCabe and Markstrom (2007). Both rely on a snow accumulation and melt model described in McCabe and Wolock (1999), whose inputs are monthly temperature and precipitation. The occurrence of snow is computed as:

(1)

where S is monthly snowfall (mm), P is monthly precipitation (mm), T_a is monthly air temperature (°C), T_rain is a threshold above which all monthly precipitation is rain, and T_snow is a threshold below which all monthly precipitation is snow (McCabe and Wolock 1999). Between T_rain and T_snow the proportion of precipitation that is snow or rain changes linearly. If snow occurs in a given month, it is added to the snowpack and is subject to melt if conditions are such that melting can occur. Snowmelt is computed as a fraction of snow storage by:

(2)

where MF is the fraction of snow storage that can be melted in a month; MF increases from 0 at T_snow to MaxMelt at and above T_rain. The MaxMelt parameter specifies the proportion of snowpack that becomes runoff (McCabe and Wolock 1999, McCabe and Markstrom 2007).

In the GR2M model, the production function relies on a soil moisture store, fed by net precipitation obtained from precipitation and potential evapotranspiration (Mouelhi et al. 2006), and snowmelt runoff. The parameter X1, the maximum capacity of the production store, is positive and is given in mm. The parameter X2 is the exchange coefficient, it is positive and has no dimension. Similarly, in the MWBM a soil moisture store receives net precipitation and snowmelt runoff; it has a soil moisture storage capacity denoted SMSC in mm. The Drofac parameter specifies the fraction of monthly precipitation that becomes direct runoff. The Rfactor parameter specifies how much surplus in a month becomes runoff. Hence the GR2M model has five parameters to calibrate and the MWBM six parameters, three of which are from the snow model. The inputs of the hydrological models are the total precipitation (P), produced by Simonneaux et al. (2008) from the liquid precipitation recorded at Tahanaout and solid precipitation at Oukaimeden, temperature (T) and potential evapotranspiration (PET).

3.2 Model calibration and assessment of parameter uncertainties

The calibration of model parameters requires the definition of a function quantifying the model error that is a measure of the difference between the observed and simulated runoff. The quantitative criterion used to measure the quality of adjustments made by the model is the multi-objective Kling-Gupta efficiency (KGE) criterion introduced by Gupta et al. (2009) and modified (KGE′) by Kling et al. (2012). As noted by Gupta et al. (2009), if solely the Nash-Sutcliffe criterion is used for calibration in basins with high runoff variability, the bias component will have a smaller contribution to the optimization of the Nash-Sutcliffe criterion, possibly leading to model simulations having large volume errors. The KGE′ combines linear correlation, bias and variability measures between observed and simulated runoff (Kling et al. 2012), and the higher the KGE′ value the better the agreement between the simulations and observations:

(3)

with:

(4)

(5)

where r is the Pearson product-moment correlation coefficient, μ is the mean runoff (m³/s), CV is the coefficient of variation and σ is the standard deviation of the runoff. The variables KGE′, r, β and γ are dimensionless and their optimum value is 1. Parameter calibration was carried out by maximizing the KGE′ using the Nelder-Mead simplex technique (Lagarias et al. 1998).

The performance of the hydrological models was estimated using different cross-calibration and validation periods to evaluate the transferability of the parameters under hydrologically contrasting conditions (see e.g. Andréassian et al. 2009). The first two methods were based on the split-sample test (SST) and the differential split-sample test (DSST) (Klemeš 1986). The SST is the most commonly used method in hydrology and involves dividing the data into two samples, one for calibration and the other for validation (here Period1 is 2000‒2009 and Period2 is 1989‒1999). The DSST consists of separating the period available into independent sub-periods and with different climatic characteristics between the periods, calibrating the model on the first period and validating it on the second, and vice versa. According to the median annual precipitation, the years were split into years wetter than the median (hereafter named Humid years) or drier than median (dry years). Therefore the wet and dry years are not necessarily continuous. In addition, a third method based on stratified sampling was considered. This method divides the time series into even years (Running1) and odd years (Running2), which allows one to account for possible trends in the climatic series (precipitation, temperature or flow). By using these three different calibration–validation schemes, six different parameter sets were obtained by calibration over the different periods considered (Period1, Period2, Humid years, Dry years, Running1, Running2) and used to assess the robustness of the model and the stability of its performance under different climatic conditions.

3.3 Estimation of potential evapotranspiration

Potential evapotranspiration (PET) is a key input variable for hydrological models, and its estimation in mountainous regions remains a major challenge (Herrnegger et al. 2012). To reduce the potential uncertainties related to the PET calculation, three formulas were compared over the period 2007‒2013 when many meteorological variables (temperature, relative humidity, wind and solar radiation) were available at the Aremd station. The three PET computations include, in order of complexity in their data requirements:

1. the Oudin (Oudin et al. 2005) formula, requiring only daily mean temperature and latitude;

2. the Hargreaves and Samani (1985) formula, requiring the average daily mean, maximum and minimum temperatures; and

3. the FAO–Penman-Monteith (Allen 1998) equation, requiring mean, minimum and maximum temperatures, relative humidity, wind speed and solar radiation.

As shown in Figure 3, the correlation between PET values computed using the three formulas is high, with correlation coefficients greater than 0.89. In particular, PET values computed using the FAO–Penman-Monteith and Hargreave-Samani equations are quite similar, as already noted by Er-Raki et al. (2010) in the plain downstream of the Rheraya catchment. The Oudin formula provides PET values slightly lower than the two other equations. When used as input for the hydrological models, no significant differences were found between monthly simulated runoff with the three different PET formulations. As previously noted by Oudin et al. (2005), conceptual hydrological models may not be very sensitive to the PET formulations as long as their calibration allows the water balance to be adjusted. As a consequence, the Oudin formula was retained for this study and served to compute a mean annual cycle of PET from the available daily temperature data (2007‒2013).

Figure 3. Correlation between PET computed using the Oudin, Hargreaves-Samani and Penman-Monteith formulas.

3.4 Construction of climate scenarios based on RCM outputs

Most often, bias-correction methods are considered with RCM outputs to produce future projections (Fowler et al. 2007). However, such methods assume that the bias in climate model outputs is stationary over time. Driouech et al. (2010), Maraun (2012) and Tramblay et al. (2013) have shown that this hypothesis could not be always validated in semi-arid regions such as Morocco. In addition, the bias correction of individual RCM variables may have an impact on the climate change signal, in particular in mountainous areas, since several variables are needed to reproduce snow dynamics (Dobler et al. 2012). Another reason not to use bias correction in the present study is the very low density of climatic gauges with long-term records (two stations for precipitation only), by comparison with the relatively coarse resolution of RCM outputs (50 km).

For the aforementioned reasons, perturbation methods were preferred for the present study. These methods rely on perturbing the observed series by a climate change factor, inferred from the comparison of RCM simulations between the historical and future periods. The so-called “delta-change” method is both simple and robust since it does not rely on any stationarity hypothesis and has been widely applied in the evaluation of the impacts of climate change (e.g. Wilby and Harris 2006, Minville et al. 2008, Dobler et al. 2012). The series of monthly averages of temperature and precipitation between the periods 1989‒2005 and 2049‒2065 were modified with a monthly change factor to obtain future climate scenarios (via an additional factor for temperature and a multiplicative factor for precipitation). The main limitation of such a method is that it does not account for the possible change in the distribution of extreme events, which can be important to consider in a semi-arid context with few precipitation days.

Therefore, in addition to the delta-change method, a quantile-perturbation method (Chiew et al. 2009), described in detail in Tramblay et al. (2013), was also used to produce future climate scenarios for precipitation. In comparison to the delta-change method, which considers only the monthly means, in the quantile-perturbation approach all the daily quantiles of the distributions are compared between the historical and the future periods, thus reflecting possible changes in extreme precipitation events. The computation was done for each month separately, between an observed period (CTRL), a simulation during the observation period (SIM) and a projected future period (PROJ). The approach followed these steps (Tramblay et al. 2013):

1. Identification of the frequency of dry days (when P < 0.1 mm) in CTRL, SIM and PROJ.

2. If the amount of dry days increases between SIM and PROJ, the corresponding amount of additional wet days in CTRL is randomly changed into dry days (100 randomly perturbed series (PERT) were generated).

3. If the amount of dry days decreases between SIM and PROJ, the corresponding amount of dry days is randomly changed into wet days in CTRL. The “new” wet days are randomly selected from wet days in CTRL (100 randomly perturbed series (PERT) were generated).

4. For each PERT series perturbed in steps (b) or (c), the empirical quantiles are adjusted by the change factor (CF) between the PROJ and SIM quantiles for the corresponding probabilities p:

(6)

4 Results

4.1 Assessment of RCM simulations in the reference period

This section assesses the quality of RCM simulations in the historical period 1989‒2005 by comparing the seasonal distributions of observed and simulated temperature and precipitation. Figure 4 shows the monthly variations in the observed temperature and those simulated by RCM available for the evaluation (EVAL) and historical periods (HIST). All RCMs reproduce the general shape of observed temperature well, but with an overestimation of 3–6°C. However, this comparison can be considered biased by the fact that RCM grid size is equal to 2500 km² while the target basin area is only 225 km². We do not have observed temperature and precipitation data spread over such an area (2500 km²) to compare the RCM simulations and observations at the same scale. Therefore this overestimation is expected as the complex orography of the Atlas is largely smoothed at the grid size of a RCM.

Figure 4. Observed and RCM temperatures during the reference period 1989–2005.

Regarding precipitation, Figure 5 shows the annual cycle simulated by the RCMs for the period 1989‒2005 compared to observed precipitation. Most RCMs underestimate precipitation, except the CNRM model, which strongly overestimates summer precipitation. For the evaluation period (EVAL), there is better reproduction of seasonality than in the historical period (HIST), with the maximum precipitation reproduced in March and the minimum in July. As observed in previous studies, precipitation simulated by RCMs is not precise enough to be used as input in hydrological models (Tramblay et al. 2013).

Figure 5. Observed and RCM precipitation during the reference period 1989–2005.

4.2 Projected changes in temperature and precipitation

In this part, we analyse the future changes in temperature and precipitation for the catchment according to the RCP4.5 and RCP8.5 scenarios between the historical period 1989‒2005 and the projection period 2049‒2065. The latter was chosen to evaluate changes at the medium term.

Figure 6 shows the projected monthly precipitation from the multi-model ensemble mean with the two RCP scenarios. Table 3 shows the climate change signal for each model, in relative values for precipitation (%) between the periods 1989‒2005 and 2049‒2065. The results show that all the climate simulations lead to a decrease in precipitation. This decrease appears particularly pronounced during winter and spring. The projected changes in precipitation are −22% in the RCP4.5 scenario and −31% for the RCP8.5 scenario. However, there are significant differences in the magnitude of these changes for different RCMs (Table 3). Note that the RCM IPSL projects the largest precipitation decrease, with −38% for the RCP4.5 scenario and −49% for the RCP8.5 scenario.

Table 3. The precipitation change signal (%) between 1989–2005 and 2049–2065 for each climate model and scenario.

	CMCC	CNRM	ICTP	GUF	IPSL
RCP4.5	−15.9	−7.8	−37.7	−11.2	−39.6
RCP8.5	−22.4	−8.0	−49.7	−31.0	−44.8

Figure 6. Projected precipitation for the future period 2049–2065.

For temperature, Figure 7 shows the projected temperature cycle according to the two RCP scenarios. Table 4 shows the climate change signal for each model, in absolute values of temperature, between the periods 1989‒2005 and 2049‒2065. It can be seen that both scenarios suggest a significant increase in temperature. On average, projected changes in temperature are +1.4°C for the RCP4.5 scenario and +2.6°C for the RCP8.5 scenario. In comparison to precipitation, these changes are relatively uniform throughout the year, with a larger increase in the summer months. During this season, the temperature changes are likely to have little impact on flows because there is almost no runoff or snow cover during the summer, while the increases in winter and spring temperatures could have a significant impact on the snow, reducing the period of snow cover and thus snow volumes available.

Table 4. The temperature change signal (°C) between 1989–2005 and 2049–2065 for each climate model and scenario.

	CMCC	CNRM	ICTP	GUF	IPSL
RCP4.5	+1.9	+1.2	+1.0	+1.4	+1.5
RCP8.5	+3.0	+2.0	+3.3	+2.3	+2.5

Figure 7. Projected temperature for the future period 2049–2065.

4.3 Hydrological modelling results

Tables 5 and 6 present the KGE′ values obtained for each calibration period (Period1, Period2, Humid years, Dry years, Running1, Running2; see Section 3.2) and the model parameter values. The KGE′ values range between 0.79 and 0.87 for calibration with the GR2M (Table 5), and between 0.64 and 0.8 for the MWBM (Table 6). In validation, the criterion values are only slightly degraded, with values between 0.44 and 0.68 for both models. The mean KGE′ over all calibration samples is 0.83 for the GR2M and 0.74 for the MWBM. On average, from the different calibration experiments, the relative bias on simulated runoff is 4% for the GR2M and 7% for the MWBM. It can be concluded from the analysis of KGE′ and model parameter values that there is relative stability of the results over the different calibration and validation periods. The drop in performance is most important when the model is validated using a standard SST on two periods (Period 1, Period 2). This is an indication that some changes affecting the rainfall–runoff relationship, such as land-use change, could have occurred in the catchment during the 1989‒2009 period.

Table 5. Calibration and validation results of the MWBM.

Calibration period	Calibration							Validation KGE′
	Rf (%)	DRf (%)	SMSC (mm)	Train (°C)	Tsnow (°C)	MaxMelt (%)	KGE′
Humid years	0.80	0.43	140.06	4.20	−0.90	0.53	0.74	0.56
Dry years	0.90	0.30	360.06	5.00	−15.00	0.79	0.64	0.50
Period1	0.54	0.45	124.91	4.41	−0.91	0.56	0.80	0.52
Period2	0.48	0.27	125.59	4.81	−14.98	0.33	0.75	0.49
Running1	0.68	0.37	100.51	4.68	−5.08	0.71	0.76	0.57
Running2	0.35	0.29	128.41	2.40	−3.51	0.51	0.75	0.61

Table 6. Calibration and validation results of GR2M.

Calibration period	Calibration						Validation KGE′
	X1 (mm)	X2 (–)	Train (°C)	Tsnow (°C)	MaxMelt (%)	KGE′
Humid years	822.31	0.99	3.95	−0.06	0.58	0.79	0.62
Dry years	285.60	1.05	4.70	−3.85	0.59	0.87	0.68
Period1	653.77	1.04	5.00	−4.52	0.50	0.81	0.49
Period2	1392.68	1.04	1.59	0.00	0.58	0.86	0.49
Running1	689.39	1.04	4.42	−0.01	0.57	0.83	0.55
Running2	1439.67	0.97	2.84	0.00	0.50	0.84	0.44

The analysis of hydrological simulations shows that both the GR2M and MWBM are able to correctly reproduce the shape of the observed monthly hydrographs (Fig. 8). The inter-annual mean monthly runoff dynamics are well reproduced; the six model runs with different parameter sets obtained for different calibration strategies encompass the observed runoff regime, except for the months of October and November (Fig. 8), when the GR2M model underestimates runoff while the MWBM overestimates it. In addition, the snow reservoir levels simulated by the model appeared realistic (as shown later, in Fig. 11), with snow accumulation between November and May, although no distributed snow data over the catchment are available to validate this result.

Figure 8. Observed and simulated mean runoff cycles for the reference period 1989–2009, according to different sets of calibrated parameters for the GR2M and MWBM.

Figure 9. Comparison of the projected changes in runoff with delta-change and quantile-perturbation methods.

Figure 10. Projected changes in runoff from five RCM simulations for the future period 2049–2065.

Figure 11. Projected changes in snow reservoir levels from five RCM simulations for the future period 2049–2065.

In order to produce runoff scenarios while keeping the whole uncertainty on the parameter estimates, no selection of the “best parameter set” was made for the projections: the hydrological model ensemble mean simulations obtained from the six parameter sets (the “ensemble mean simulation” as shown in Fig. 8) was considered as the reference simulation for the historical time period.

4.4 Future projections on streamflow and snow amounts

As mentioned above, the monthly runoff values simulated from the observed data with six sets of parameters for each hydrological model were averaged to produce a single simulation at the monthly time step, reflecting the uncertainty related to the transferability of the hydrological model parameters towards climate-contrasted conditions. These reference simulations (shown in Fig. 8 by the curves “GR2M ensemble mean” and “MWBM ensemble mean”) were then compared to future hydrological projections, obtained from the two hydrological models driven by perturbed observed climatic series, to assess the changes in monthly flows. The time series of observed precipitation and temperature were perturbed, according to a change signal simulated by the different RCMs between the 1989‒2005 and 2049‒2065 periods, using the delta-change or quantile-perturbation methods (see Section 3.4). The time series obtained were then used as input in the GR2M and MWBM, with the six sets of parameters described above, to obtain for each scenario RCP4.5 and RCP8.5 an average monthly runoff for the period 2049–2065.

Before analysing in detail the potential impact of climate change on catchment runoff, we first compare the impact on runoff of projected changes in precipitation only, based on the two perturbation methods tested: the delta-change by monthly factors and the quantile-perturbation method considering the full distribution of precipitation changes. Figure 9 shows the simulated runoff changes due to precipitation changes only, with the GR2M and MWBM according to the two methods for the two RCP scenarios at the mid-term horizon. The changes are almost identical according to the two methods, with for the RCP4.5 scenario −34% projected for the delta-change method and −30% for the quantile-perturbation method, and −68% and −66% respectively for the RCP8.5 scenario. This indicates that both methods give very similar results at the monthly time scale. The following results are presented using the delta-change approach with monthly change factors only.

Figure 10 shows the relative changes in monthly runoff simulated between the periods 1989‒2005 and 2049‒2065 for all RCMs and for both RCP scenarios. Overall, most climate scenarios lead to a significant decline in spring flows (March–May). Some RCMs suggest an increasing runoff for some seasons, such as GUF, CMCC and ICTP under the RCP4.5 scenario in late spring and summer. Considering the projections resulting from each of the RCMs (Fig. 10), high uncertainty is observed between the different models, although on average all agree on a decline. Under scenario RCP4.5, the projected decline in runoff is −38.5% for the GR2M and −19.2% for the MWBM, and with scenario RCP8.5, the projected decline in runoff is −62.9% and −37.5%, respectively.

The decline suggested by the RCP8.5 scenario is more marked than for the RCP4.5 scenario, especially in spring and autumn, while in summer and winter projected changes are of lesser amplitude. Figure 11 shows the changes in monthly levels of the snow reservoir simulated in the hydrological models between 1989‒2005 and 2049‒2065. Due to the increased temperature in the degree-day component of the snow module, the model simulates a significant decrease in snowfall in favour of liquid precipitation. This necessarily leads to a decrease in snow amounts, while the snow period is reduced. These results indicate that the impacts of climate change could be severe by 2050, with a major reduction in surface runoff during the spring, summer and autumn, and a considerable reduction in snow cover.

It should be noted that future changes in PET were not considered, since the formula of Oudin et al. (2005) is based only on temperature. Indeed, with such a formula, increasing temperature would increase PET, but several studies have shown that, despite increased temperature, a decrease in the evaporative demand could be observed related to a decrease in the radiation, humidity or wind (Johnson and Sharma 2010, Sheffield et al. 2012). Nevertheless, sensitivity tests have shown that it is mainly the joint decrease of precipitation with increase in temperature that is reducing the runoff in the Rheraya catchment. If in addition to changes in precipitation and temperature an increase in PET is considered, computed from the climate scenarios according to the Oudin formula, the impact on runoff would be minimal (with an additional change of −3% to −10% in runoff, according to the two scenarios, not shown).

5 Discussion

The strongest uncertainties are most often related to the choice of climate models and scenarios, as noted by many studies (Minville et al. 2008, Chiew et al. 2009, Quintana Seguí et al. 2010, Dobler et al. 2012, Majone et al. 2012). However, hydrological model uncertainty could also have an important role (Bosshard et al. 2013). In the present work, an ensemble of simulations was considered at a moderate spatial resolution (50 km). More simulations at higher resolution (12 km) are becoming available for the Euro-Mediterranean domain, but a recent study by Tramblay et al. (2013) in north Morocco showed that projected changes were similar at the 12 km and 50 km resolutions. However, this needs to be tested in further studies and in other catchments when a sufficient ensemble of high-resolution RCM simulations becomes available. Additionally, in the present study two perturbation approaches, based on monthly changes and on the daily quantiles of the distributions, were compared. Other downscaling approaches could also be tested to provide a broader range of projections according to different methods. In particular, there is a need to adapt the bias-correction methods in the semi-arid context, to overcome issues related to the validation of quantile-mapping bias-correction methods in semi-arid regions and to account correctly for extreme events (see e.g. Maraun 2012, Tramblay et al. 2013).

This study provides climate-driven quantitative estimates of the projected changes in runoff for the Rheraya basin; consequently, these results need to be complemented and completed by other studies with different hydrological or land-surface models, able to consider the possible feedback effects between the atmosphere and the land surface. In particular, the snow dynamics are modelled here via quite a simple degree-day approach, and the results of the present study could be compared to further work considering more complex energy-balance models (Boone and Etchevers 2001, Liston and Elder 2006, Magand et al. 2014). However, the application of such models in a data scarcity context remains challenging. As noted by Dobler et al. (2012), an improved understanding of the physical processes governing future river flows can help to better estimate the uncertainties related to this type of impact studies. In addition, no changes in land use have been considered in the present study. Since most parts of the catchment exhibit steep slopes and high elevation, the vegetation cover remains limited and the stability of the hydrological model results indicated no major land-use changes occurred to modify the rainfall–runoff processes. However, the lowest elevation areas near the stream bed are cultivated, mainly as orchards; therefore, changes in land use or agricultural practice may have an effect on runoff processes in the future (Cristea et al. 2014) and there is a need to test models able to account for this type of anthropogenic activity (Szczypta et al. 2015).

In the present study PET was not modified in the future scenarios. However, we tested the simple hypothesis of an increase in PET proportional to the temperature increase. The impact of PET changes on runoff was less notable than that of precipitation changes. As observed by Kingston et al. (2009), most climate change studies indicate an increase in PET associated with a warming climate, but very strong differences can be found between the PET models (Sheffield et al. 2012). As shown by Sheffield et al. (2012), a decreasing trend in the global average of the Palmer drought severity index (PDSI) is detected when the PDSI is computed with PET derived from temperature only, but no trend is found when PET is computed with the Penman-Monteith equation. In semi-arid areas, water losses via evapotranspiration for vegetation and via sublimation for snow are important and not properly taken into account by simple conceptual models (Boudhar et al. 2009). In particular, feedback effects are neglected using this type of model; indeed, with increasing surface temperature the sublimation could decrease and then produce more ablation by melting, and as a consequence reduce the projected decline in surface runoff. As noted by Adam et al. (2009), earlier snowmelt could increase water availability at a time of year when the evaporative demand is relatively low, which could have a negative feedback on annual evapotranspiration. In contrast, loss of snowpack could lead to an earlier removal of snow cover and an increase in albedo, leading to a positive feedback on evapotranspiration due to increased net radiation. Therefore, there is a need to take into account additional variables linked to the evapotranspiration processes to produce future scenarios (Johnson and Sharma 2010).

6 Summary and conclusions

This study is the first attempt, to our knowledge, to evaluate quantitatively the possible future climate change impacts in a mountainous catchment in North Africa. The study focused on the Rheraya catchment (225 km²) located in the High Atlas Mountains in southern Morocco, where sufficient measurements over a relatively long time period were available to perform such a study with simple yet robust methods. Two conceptual models including a snow module were considered to reproduce the monthly dynamics of runoff in the catchment, and the uncertainty of its calibrated parameters was assessed using six contrasted periods for model calibration and validation. Model performance was found adequate to reproduce monthly runoff dynamics, considering the available data. To provide future scenarios, an ensemble of five regional climate models from the Med-CORDEX initiative was considered together with a delta-change and a quantile-perturbation method. Using the future climate simulations with the hydrological model revealed that surface runoff is expected to decline at the mid-term horizon in the Rheraya catchment. The evaluation of simulated precipitation by RCMs showed a strong underestimation of about 50%, but a good reproduction of the temperature cycle. The future projections for 2049–2065 under the two scenarios RCP4.5 and RCP8.5 showed increasing temperature (+1.4°C to +2.6°C) in conjunction with a decrease in total precipitation (−19% to −31%). Projected changes indicated a decrease in surface runoff (−19% to −63%, depending on the scenario and hydrological model) and a significant decrease in snow amounts, related to reduced precipitation and increased temperature. Despite the potential limitations of this type of approach, this is a first step towards the assessment of future projected changes in surface water resources in Morocco using climate models, which needs to be complemented by other studies comparing different hydrological and climate models.

Acknowledgements

The authors are grateful to the Water Management Agency of the Tensift Basin (ABHT) in Marrakech for providing the necessary data for this study.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work has been carried out within the framework of the Joint International Laboratory TREMA (http://trema.ucam.ac.ma). In addition to the CNRST scholarship awarded to Marchane Ahmed, this study was also supported by the ENVI-Med CLIHMag (Changement cLimatique et Impacts Hydrologiques au Maghreb) project funded by the programme INSU-MISTRALS. This work is part of the Med-CORDEX initiative (www.medcordex.eu) supported by the HyMeX programme (www.hymex.org).