Present and Future High ‐ Resolution Climate Forcings over Semiarid Catchments: Case of the Tensift (Morocco)

: In semiarid areas, the climate will be accompanied by a marked decrease in precipitation ( − 16% for RCP8.5). These present and future spatialized data sets should be useful for impact studies, in particular those focusing on water resources.


Introduction
In semiarid areas, most of the available water resources for the production of hydroelectric energy, for drinking water, as well as for irrigated agriculture located downstream in the surrounding plains [1][2][3][4], comes from the mountainous regions [5]. This role of the water tower is nowadays threatened by the increase in water needs due to population growth and higher standards of living, as well as by the intensification of irrigated agriculture [6] and by climate change. The south Mediterranean region is now well known as a "hot-spot" for the latter [7][8][9] and there is reasonable evidence indicating that mountainous regions should be more impacted than their surrounding plains [10]. Within this context, the development of monitoring systems of climate variability and change in the higher elevation area of the southern Mediterranean regions are not only crucial for studying the climate variability and trends but also for proper water resource management.
The climate of the south Mediterranean mountain ranges faces a high spatiotemporal variability owing to steep precipitation and temperature gradients associated with a marked orography, to a strong seasonality of climate with most of the precipitation concentrated during the winter months [11,12] and, for the western regions, to the joint and complex influence of the synoptic flow originating from the Atlantic Ocean, from the Mediterranean sea and from the Sahara [13]. This high spatiotemporal variability, in addition to the inherent difficulty of access and the hard measuring conditions in mountainous regions, leads to strong observational issues: gauge network and radar observation are poor or nonexistent. This has led the scientific community to new developments, including new observation concepts based on remote sensing data and reanalysis systems merging atmospheric modelling with the available observations [14].
There are various data sets of precipitation estimates originating from different sources of remote sensing observation, both with and without gauge data, which have recently become available (see Sun et al. 2018 [15] and Levizzani and Cattani, 2019 [16] among others for a review). They are even sometimes preferred to the gauge network data thanks to their global coverage [17]. The methods of precipitation estimates are based on Visible-Infrared (IR) observations sensing cold cloud with a high vertical development associated with convection and/or on passive and active (radar) microwave (MW) data directly sensitive to precipitation particles [18]. The most accurate products merging IR and MW data [19] includes TRMM Multi-satellite Precipitation Analysis [17], PERSIANN [20] and CMORPH [21] among others. Despite the continuous improvement of the retrieval algorithms, several studies indicate that although these products perform reasonably well at monthly and larger time scales, they also display large uncertainties in complex terrain [22][23][24][25][26]. A few studies dedicated to the evaluation of these satellite products over the Mediterranean area draw very close conclusions [27][28][29].
Another appealing tool providing high-resolution gridded precipitation data together with other meteorological variables that cannot be observed directly using spaceborne sensors is the reanalysis system merging atmospheric model prediction with available ground observations (see Table S1 for a non-exhaustive review of the different existing systems in the Supplementary Materials). Despite a continuous increase in spatial resolution, the reanalysis data sets freely available for the whole world (NCEP, [30] or ERA-5, [31]) or at the regional level like MESAN [32,33] still do not meet the resolution requirement for mapping the high spatial variability of the climate in these complex terrains. By contrast, several tools were proposed for the spatialization of climate variables based on the merging of a "background", usually the forecast of an atmospheric model, and observations from the ground network. The MESCAN system [34] provides gridded precipitation only over the French territory. By contrast, the SAFRAN [35,36] system, a French acronym standing for "Système d'Analyse Fournissant des Renseignements Adaptés à la Nivologie", analyses temperature, specific humidity, wind speed, precipitation and cloud cover. In addition, it was initially developed for providing gridded meteorological variables at a high-resolution of 8 km [37] for the operational forecast of avalanche risks in mountainous areas. SAFRAN was also successfully implemented and evaluated in different geographical areas characterized by a dense network of ground observations: in France [38][39][40][41], in Spain [42] and in Tunisia [43]. A first question arises about the performance of such a system when the ground network is scarce, a very likely situation for numerous catchments in the south Mediterranean region such as in Morocco and Algeria. Many watersheds in these regions are also characterized by a very steep orography while its representation in the atmospheric models is often a key issue [44]. Indeed, a large part of the discrepancies between reanalysis and observations in complex terrain can be attributed to the large differences of elevation between a grid cell and the ground station [40,45]. Within this context, which spatial resolution should be chosen to properly predict the rainfall/snowfall partition in the higher elevation areas? MicroMet is also a quasiphysical system designed to produce high-resolution atmospheric forcing up to 30 m in mountainous areas [1]. Nevertheless, by contrast with SAFRAN, it is based on extrapolations from meteorological station data only. It was also positively assessed in regions where observational data do not exist, and where precipitation quantities are generally very important: Arctic Alaska, central Norway, Svalbard, Greenland and Antarctica, Wyoming, Colorado, and Idaho [46][47][48][49][50][51][52][53].
The objective of this study is twofold: (1) to evaluate and adapt the SAFRAN analysis system to a complex catchment typical of the south Mediterranean region, including an area of water production in the mountain with a marked orography and a downstream plain-the comparison to the MicroMet data-driven approach is also carried out; (2) to disaggregate the future climate scenarios to provide spatially explicit fields of precipitation and temperature.

Study Area
The study was carried out on the Tensift watershed (in the region of Marrakech, Morocco). Due to its location and geographical characteristics as well as its semiarid to arid climate, the Tensift catchment area is considered to be representative of many basins in the Southern Mediterranean. It has been the subject of meteorological, micrometeorological and hydrological monitoring since the early 2000s within the framework of the Tensift observatory co-led by the Cadi Ayyad University of Marrakech and the CESBIO (Toulouse).
The Tensift catchment is located in central Morocco ( Figure 1) and extends over 19,800 km 2 . The climate is governed by the complex orography, the Western flows bringing humid and cold air from the Atlantic Ocean, and is locally modulated by convective situations. The basin can be separated into two zones with contrasted hydrological functioning:  The southern part includes the northern slope of the High Atlas; its highest point (Jbel Toubkal) is located at 4167 meters above sea level. Average annual precipitation varies between 600 and 700 mm/year [54,55] and snowfall is localized in the highest elevation areas above 1800m [3,56]. The melting of the snow stored in winter contributes to the shift of flows towards spring and the support of base flows during summer.


The second part is the Haouz plain, characterized by average annual precipitation varying between 200 and 400 mm/year, located towards the west of the reliefs. The climate of this zone is semiarid. This part supports intensive irrigated agriculture, whose main crops are olives, citrus fruits, and wheat.

Meteorological Observations
The ground network in the Tensift catchment is sparse. Within this study, several data sets originating from different institutions have been gathered. The location of the stations is provided in Figure 1. Table 1 displays its characteristics. The stations provide either daily data including cumulative rainfall RR, minimum Tmin and maximum Tmax temperatures or 6-hourly observations of rainfall, pressure, relative humidity, wind (strength, direction), temperature and cloud cover. These data cover a hydrological period of 10 years from 1 August 2004 to 31 July 2014. The period was chosen to optimize for data availability.
The data was quality-controlled [57][58][59][60] to detect outliers. Note that most of the stations in the mountains are installed in a subcatchment of the Tensift named the Rheraya because it has been the support of several hydrological experiments since the beginning of the 2000s [12,54]. Finally, daily snowfall (mm/day)acquired at the Oukaimeden station from 1 August 2004 to 8 June 2011 are also used for the evaluation of the rainfall/snowfall partition by SAFRAN and MicroMet. In this study, Euro-CORDEX [61] climate change scenarios are used. They were designed to provide high-resolution scenarios (12 km) over a large Mediterranean window. A comprehensive evaluation of the Euro-CORDEX scenarios confirmed the ability of the Euro-CORDEX to reproduce the basic characteristics of the Mediterranean climate [61]. We selected four regional climate model (RCM) simulations from the Euro-CORDEX initiative with a 12 km horizontal resolution based on two contrasted representative concentration pathways (RCPs) (4.5 and 8.5). RCP4.5 is characterized by an atmospheric CO2 concentration that peaks in 2040 and decreases thereafter, while the CO2 concentration in the RCP8.5 scenario continues to increase throughout the 21st century. These simulations are provided by different institutes: the CNRM (ALADIN5.3 regional climate model driven by CNRM-CM5 general circulation model), KNMI (RACMO2.2 driven by EC-EARTH), IPSL (WRF3.3.1 driven by IPSL-CM5) and MPI (REMO2009 driven by MPI-ESM). Each model has two types of climate simulations: (1) historical runs (HIST) or control period available from 1960 to 2005; (2) scenario runs: forced simulations by General Circulation Models (GCMs) according to the selected emission scenario (RCP4.5, RCP8.5) available from 2006 to 2100 with RCMs. For our study, the chosen control period is 1988-2005 which will then be compared with the 2041-2060 horizon.
The variables processed in this study are maximum temperature (Tmax), minimum temperature (Tmin) and precipitation (pr). They are extracted with a daily time step at 12 km resolution for the Tensift catchment. The reanalysis of SAFRAN was also aggregated to a daily time step in line with Euro-CORDEX by a simple average or accumulation for precipitation in order to have a futuristic projection of SAFRAN forcing at the 2041-2060 horizon for the two scenarios.
For illustration purposes, the historical runs are compared to the observations from the stations of Marrakech and of Oukaimeden_CAF (see Figure S1 of the Supplementary Materials) during the historical period through average yearly cycles. The seasonal variability of temperature is well represented by all RCMs apart from some biases that can reach several degrees (see Oukaimeden during the winter month for instance). For precipitation, by contrast, stronger discrepancies are observed between the different RCM and between the RCM prediction and the observations. Some of the RCMs (CNRM) are not even able to reproduce the winter peak of precipitation. Those discrepancies have already been observed and have motivated several studies [62,63] to apply a bias correction approach.

Daily Snow-Covered Area MODIS
In order to evaluate the ability of MicroMet and SAFRAN to predict snowfalls, the daily snow-covered fraction (SCF) product MOD10A1 at 500 m resolution version 6 derived from the TERRA satellite data is used. The daily products provided by the National Snow and Ice Data Center (NSIDC) are extracted from August 2004 to July 2014 on the Tensift and on the Rheraya sub-basin. Version 6 of MOD10A1 does not include the snowcovered fraction anymore. In this study, the Normalized Difference Snow Index (NDSI (1); [64]) is used following (Hall et al. [65]; Baba et al. [66]) to compute SCF as follows: In a second step, a filter that searches for clear neighbors in space and in time is applied to decrease the number of cloudy images and to minimize the problem of cloud/snow misclassification following Marchane et al. [56]. The total Snow Cover Area (SCA) at the scale of the basin is then obtained by adding up the SCF and multiplying by the area of the MOD10A1 pixel (500 × 500 m²).

SAFRAN
SAFRAN analyses the low level of 2 m temperature, 2 m specific humidity, 10 m wind speed, precipitation, and cloud cover. The analysis between the observations and the first estimate field is carried out based on the Optimal Interpolation (OI) of a forecast and observations. The forecast was extracted from the weather forecast system used operationally at the Direction Générale de la Météorologie (DGM) named ALADIN-MO-ROCCO METEO model [67]. The analysis is carried out at 6-hourly time steps and then interpolated at hourly time steps according to different extrapolation methods. The precipitation is analyzed at the daily time steps and then redistributed throughout the day using the diurnal specific humidity cycle. The partition between rain and snow is carried out using the 0.5 °C isotherm. The analysis is performed on climatologically homogeneous areas, which are determined based on the elevation. SAFRAN assumes that within these climatic zones, the meteorological variables vary according to topography only. Stated differently, only one vertical profile of the meteorological variables for each climatic zone is analyzed and the variables are distributed on a regular grid of 8 km according to topography [38]. The main limitation of this zonal approach is the creation of artificial discontinuities at the borders of the zones [42]; these discontinuities are particularly strong for the precipitation variable. Incoming radiation in the visible and infrared range is predicted following the model of Ritter and Geleyn (2002) [68] because observations are often lacking for these two variables. For more details, please see Quintana-Segui et al. (2008) [39].

MicroMet
MicroMet [1] has been developed to produce high-resolution weather data of air temperature, surface pressure, relative humidity, incoming solar radiation, incoming longwave radiation, precipitation, wind speed and direction based on ground observations. It is part of the distributed snow model package developed by Liston et al. (2018) [69]. Unlike SAFRAN, MicroMet is "data-driven": it relies on the spatial interpolation of observations using the Barnes objective analysis scheme [1,[70][71][72] and the spatial distribution of variables is then corrected using values of temperature-elevation, wind-topography, humidity-cloudiness, and radiation-cloud-topography relationships. The system assumes that there is at least one observation for each of the meteorological variables (air temperature, relative humidity, precipitation, wind direction and strength) somewhere in the study area, for each time step, except for the incoming short-wave and long-wave radiation and for pressure. For these latter variables, MICROMET uses simple models described in Liston et al. (2006) [1] for the spatialization.

Implementation and Assessment
SAFRAN and MicroMet Implementation SAFRAN's 23 climatologically homogeneous zones, ranging from 306 km 2 to 1457 km 2 , were identified based on elevation using the Digital Elevation Model (DEM) 30m ( Figure 1). Figure 1 clearly illustrates the scarcity of the network with only 20 stations, mostly located in the southwest of the area: 17 zones among 23 have no observation available, though a minimum of 2 stations per zone is recommended for the implementation of SAFRAN [39]. Two implementations of SAFRAN are performed: one based on the projection on a regular grid of 8 km and another one with a variable grid to better represent the orography and assess the impact on the rain/snow partition and the temperature fields. Indeed, with a regular grid of 8 km, the elevation range within one grid cell can reach up to 3726 m. The variable grid is built up as follows: the grid size is divided by two as long as the elevation range is higher than 300 m down to a 1 km grid size.
MicroMet is implemented for the purpose of evaluating SAFRAN on the Rheraya mountainous subcatchment only with a 1 km resolution. Regarding the dense network on the catchment, including measurements of shortwave and longwave incoming radiation, it can be considered as a reference field for SAFRAN. The default value of the lapse rate and the precipitation gradient were replaced by those computed from the observation network.

SAFRAN Assessment
The SAFRAN assessment is carried out in three ways: (1) by a station-level validation based on the leave-one-out (LOO) approach [73] from August 2004 to July 2014. It consists of removing a station from the database, running SAFRAN without this station and evaluating the SAFRAN reanalysis at it. This validation method is used for a mountain station (Aremd) and the plain station (Chichaoua). Both stations were chosen due to data availability; (2) via a comparison with MicroMet monthly averages over the Rheraya catchment; (3) by comparing the rain/snow partition to the snow surface area derived from the daily MODIS product and to the in situ snow depth measurements.

Scenario Downscaling
The main spatiotemporal patterns of precipitation and temperature of the Mediterranean region are reasonably well reproduced by the Euro-CORDEX runs as shown by [74,75]. Nevertheless, some systematic discrepancies between RCM runs and local observations are also highlighted because of the inherent error in the model parameterization or in the representation of orography. The reduction of biases through the implementation of downscaling methods is thus a prerequisite for impact studies. Our objective is to provide high-resolution future climate scenarios spatially explicit for the 2041-2060 period. Only precipitation and temperature are considered for downscaling in this study. The quantile-quantile (Q-Q) approach was chosen [76]. The Q-Q is a nonlinear bias correction method aiming to correct the biases between the quantiles of the distribution of the RCM outputs and the observations. Note that this method has been applied with success, including for our study region [62,77,78].
As the period of availability of SAFRAN is too short (10 years [79]), the Q-Q approach was first applied to the Marrakech (elevation: 464 m) and the Oukaimeden_CAF (elevation: 2687 m) stations that provide daily data from 1988 to 2005 (control period of Euro-CORDEX). The correction was applied to each season (3 month periods) separately and each model was corrected independently, in order to evaluate the scattering between the different RCMs. Based on an elevation threshold equal to 1000 m, the grid points of Euro-CORDEX scenarios are separated into plains and mountains. The correction function obtained at the Marrakech station was then applied to the plain grid points while the one obtained at the Oukaimeden_CAF was used to correct the mountain grid points. Finally, the "Delta change" approach [80][81][82] is used to project the SAFRAN forcing at the 2041-2060 horizon for the two scenarios. The delta change simply consists in computing the monthly differences between the disaggregated Euro-CORDEX scenarios using the Q-Q approach and the history runs of Euro-CORDEX and applying these changes (delta) to SAFRAN considering the nearest Euro-CORDEX grid point.

At the Station Level
The results of the LOO approach are shown for the 8 km regular grid and for the new irregular grid.  At the plain level, the variable grid mesh based on orography has almost no impact on the SAFRAN reanalysis. Indeed, the results of the 8 km grid and of the irregular grid superimposed almost perfectly for the three variables ( Figure S2). Both are also in very good agreement with the observations of the station, even considering rainfall. This probably means that the ALADIN forecast used as a first guess is good, at least with regards to monthly cumulative rainfall, as there is no station in the surrounding area that could have been used by the SAFRAN system to supplement the removed station of Chichaoua. Some rainfalls predicted by SAFRAN during summer in 2009, 2010 and 2013 are not observed by the Chichaoua station. They probably originates from localized convection events that were well caught by the ALADIN model. The incoming solar radiation (ISR) is also perfectly reproduced with a correlation coefficient of R = 0.96 for both implementations of SAFRAN apart from a small bias in summer. This is remarkable as it is not analyzed in SAFRAN but calculated by a radiation scheme. In line with the results obtained by [42] over Spain, close conclusions can be drawn here for the air temperature with a bias lower than 0.5 °C and a reasonable RMSE of 2.8 °C on average over the 2004/2014 period over this station surrounded by a flat terrain, although a significant positive bias can be observed during winter 2010-2011 and 2011-2012. By contrast, the statistical metrics for wind speed are significantly worse, with a strong negative bias above 2 m/s and RMSE reaching 3.5 m/s. The wind is known to be highly variable in space and the metrics are close to Quintana-Seguí et al. (2008) [39] who also found an overall trend to under-estimation attributed to the linear interpolation at the hourly time scale that could mask some important atmospheric phenomenon such as breeze effects.

Aremd Station (2058 m)
A different behavior can be observed at the Aremd station. Although SAFRAN for both implementations reproduced the seasonal cycle of ISR, air temperature and rainfall rather well, the amplification of the error on this mountainous station with regards to Chichaoua appeared clearly for all variables (see Figure 2 and Table 2). First, a strong positive bias reaching almost 60 W/m 2 is observed for ISR with higher biases in summer. This stands for both implementations of SAFRAN with no clear added-value of the irregular grid on this variable. MicroMet overestimation of ISR is even higher than SAFRAN, with a bias and an RMSE of 86.2 W/m 2 and 351.7 W/m 2 , respectively. Interestingly enough, a positive bias was also highlighted on SAFRAN reanalysis for a mountainous station located in the French Alps by Quintana-Seguí et al. (2008) [39] while plain stations were characterized by slight negative biases. It is well known that shortwave radiation is strongly influenced by topography including slope [83] and this effect is not taken into account in SAFRAN. The large bias of ISR could also suffer from poor prediction of cloudiness by SAFRAN, as highlighted by several authors [39,40,42,84]. By contrast, the irregular grid outperformed the 8 km grid for air temperature, precipitation, wind speed and, to a lesser extent, specific humidity ( Table 2). The lower bias and RMSE for precipitation and temperature when using the irregular grid was expected as the grid cell elevation is closer to the station elevation. MicroMet provided slightly better statistical metrics in terms of the RMSE and the correlation coefficient, but the results are close to the SAFRAN irregular grid results.

Catchment Scale (Rheraya)
As the irregular grid outperformed the regular grid, only the results obtained with the variable mesh grid will be compared to MicroMet. The comparison over the Rheraya catchment is shown in Figure 3 (cumulative values for precipitation and the annual average for the other variables) for the period 2004-2012. Table 3 displays the annual average of the meteorological variables for SAFRAN and MicroMet. Table 3. Annual average (or cumulative for rainfall and snowfall) of the meteorological variables predicted by SAFRAN and MicroMet over the Rheraya sub-basin.
The most striking feature is that SAFRAN provided significantly smoother fields than MicroMet. This has been already observed and attributed by several authors to the use of zones that are assumed to be climatologically homogeneous [39,42,84]. For temperature, the spatial patterns of SAFRAN and MicroMet differed slightly, with a stronger altitudinal gradient for MicroMet, while the average air temperatures are close at the scale of the Rheraya (11.3 °C versus 11.7 °C for SAFRAN and MicroMet, respectively). The larger discrepancies are observed for the ISR with significantly higher values of about 30 W/m 2 on average over the catchment (Table 3) predicted by MicroMet and very variable spatial patterns. This is probably due to the representation of the slope and aspect effects in MicroMet, while these effects are not considered by SAFRAN. This difference could have significant impacts for snow prediction studies including the partition between snowmelt and sublimation [2,84].  The average yearly amounts of rainfall and snowfall by elevation range are displayed in Figure 4. On average over the catchment, the differences in the precipitation amounts are highly variable between the two models; MicroMet systematically predicts higher pre-cipitation than SAFRAN, except for in the lowest elevation ranges. There is also a significant difference in the partition between solid and liquid precipitation between the two models. MicroMet predicted significantly more rainfall than SAFRAN on the highest part of the catchment, while the opposite is true for the lowest elevation area to the north. By contrast, SAFRAN recorded a significant amount of snow spread over the southern half of the basin, while the snow prediction of MicroMet followed a very steep elevation gradient: a very high yearly snowfall of around 450 mm/year is observed for the extreme highest part above 3500 m, but the average amount decreased significantly when going northwards towards lower elevations.

Snow Assessment
In this section, we will focus on assessing the snowfall/rainfall partition prediction by comparing the predicted snowfall to the snow covered area derived from MODIS and to the in situ snow depth. Note that this is a very indirect validation, as the snow covered area (SCA) and snow depth result from complex processes of accumulation, melt and possibly sublimation in this semiarid area. The Snow Covered Area derived from SAFRAN and MicroMet obtained by summing up the area of all the pixels where snowfall occur is compared to the daily MODIS products on the Rheraya sub-basin (Figure 5a). Figure 5b is the same but over the Tensift catchment and considering SAFRAN only. The snowfall amount from SAFRAN and MicroMet is then compared with in situ data recorded at the Oukaimeden_CAF station from 2004 to 2011 (Figure 5c). The way SCA is computed (i.e., with no representation of the accumulation and melting processes) probably explains why both systems failed on average to properly reproduce the observed snow cover fraction during the core of winter. For instance, the high observed SCA in 2006-2007 or 2007-2008 probably originated from continuous accumulation of snow thanks to regular events and low melt since the beginning of the snow season. By contrast, the SCA peak predicted by both systems in [2004][2005][2008][2009] and, to a lesser extent, in 2013-2014 (SAFRAN only), may result from several abundant snowfalls during the month, while the mild temperature may have favored melting, explaining the low observed monthly SCA. Nevertheless, for the remaining years, the results remain very good, with low discrepancies between the models (SAFRAN-irregular and MicroMet) and the observations. In particular, the beginning and the end of the snowy season are properly reproduced by both systems. The biases are around +0.3% for SAFRAN and −0.2% for MicroMet, while the correlation coefficient is equal to R = 0.74 for SAFRAN compared to R = 0.72 for MicroMet (both are significant at the 99% level). At the scale of Tensift (Figure 5b), similar conclusions can be drawn, with a slight overestimate in some hydrological years (2004-2005, 2005-2006 and 2008-2009) as already discussed. The snow fraction estimated by SAFRAN is still highly correlated with MODIS observations (R = 0.73) at this scale with a very low bias (+0.4% over the entire study period).
Finally, snowfall in mm/day from MicroMet and SAFRAN is compared to daily observations gathered at the Oukaimeden station in Figure 5c. The overall monthly dynamic is reasonably reproduced by both systems (R = 0.54 and 0.75 for SAFRAN and MicroMet, respectively) considering the difference of scale [56] between the 1 km grid point of SAF-RAN and MicroMet and the point scale observations. In particular, the start and the end of the snowy season predicted by both systems are in line with the observations. Nevertheless, large negative biases are observed during the winter months, particularly in 2005-2006 and 2008-2009. This is attributed to the location of the grid points of both systems exhibiting an average elevation of 2676 m and 2373 m for SAFRAN and MicroMet, respectively, while the elevation of the station is of 2687 m. This is in agreement with the results obtained at the scale of the Rheraya catchment showing that MicroMet predicted a significantly higher amount of snowfall than SAFRAN. Indeed, by considering a grid point with an elevation closer to the station, it is likely that the strong bias observed for MicroMet would be significantly reduced.
In conclusion, despite the identified limitations mainly attributed to the complex orography, SAFRAN and MicroMet predict the snowfall reasonably well.

Downscaling of Climate Change Scenarios
The projection of SAFRAN forcing is conducted in three steps corresponding to the three subsections below: (1) disaggregation of the climate scenarios using the Q-Q approach based on stations data; (2) computation of the spatialized delta-change between historical and future Euro-CORDEX runs after Q-Q correction; (3) futurization of SAF-RAN using the spatialized delta change values.

Projection of Euro-CORDEX Based on the Q-Q Approach
As a first step, the Q-Q debiasing functions are determined on the historical period for the two stations of Marrakech in the plains and Oukaimeden_CAF in the mountains. Future projections are then derived for the two stations using the Q-Q functions and the nearest grid point of Euro-CORDEX for RCP4.5 and RCP8.5 at the horizon 2041-2060. The Q-Q functions are applied to the Euro-CORDEX predictions by separating plains and mountainous grid points. Figure 6 displays the Euro-CORDEX run during the historical period after debiasing using the Q-Q approach superimposed with the in situ observations gathered at the plain station of Marrakech and at the mountain station of Oukaimeden_CAF. The nearest grid point of the two stations is used for Euro-CORDEX. The discrepancies between Euro-CORDEX temperatures ( Figure S1) are almost perfectly corrected by the Q-Q approach even if some limited biases lower than 1 °C remain for some winter months (January and February), in particular for the mountainous station (Figure 6d,e). Interestingly enough, scattering disappears between the different RCMs. These results are in line with other studies that also highlighted the good performance of the Q-Q approach for the temperature variable [62,77,78]. By contrast, for precipitation, there is still a significant scattering between the different RCM after Q-Q debiasing and between the RCM and the observations (Figure 6c,f). For instance, monthly amounts ranged between 26 mm/month and 43 mm/month for January at the Marrakech station while observations are slightly lower with 25 mm/month recorded on average during the period from 1988-2005. Likewise, the precipitation peak was observed in March while the RCM runs after the Q-Q correction exhibit peaks occurring earlier in the season, either in February (ECEARTH, MPI) or in January (CNRM, IPSL). Nevertheless, the Q-Q approach is efficient in providing RCM runs with monthly precipitation amounts reasonably close to the observations compared to the raw predictions (see, for instance, IPSL in Figure S1). Likewise, the seasonal dynamic is in much better agreement with the observations (see, for instance, CNRM before debiasing in Figure  S1).  Table 4. A high degree of warming is observed for both stations, as evidenced by the negligible scattering of the four RCMs, and also because the observations are outside of the hatched section showing the standard deviation of the RCMs. This was already observed by several authors in the South Mediterranean area [62,77]. Interestingly enough, stronger warming may be expected during the winter and spring months than for the rest of the year. This is even more prominent for the mountainous station, which could face warming of up to 4.3 °C and 3.2 °C in March for RCP8.5 and May for RCP4.5 scenarios, respectively. As for precipitation, even if a drought is predicted, in particular for the more pessimistic scenario, it is much more uncertain, as the observations fall within the hatched section for most of the year-apart from November for the Marrakech station and the months of November and December for the Oukaimeden station. An increase in precipitation reaching +8.3 mm/month at the Marrakech station and +15.1 mm/month at the Oukaimeden station can even be observed in January and February. The start (November-December) and the end of the rainy season The Q-Q functions are then applied to the Euro-CORDEX historical runs by separating plain grid points (below 1000 m) from the mountainous grid points (above 1000 m). The (strong) assumption underlying this approach is that while not being able to reproduce the seasonal dynamic and the yearly average of the three considered variables as already observed ( Figure S1), the RCMs can predict the main pattern of the spatial variability. As orography plays a significant role in the dynamic predicted by the RCMs together with the large synoptic flows originating from the Atlantic, it may be a reasonable assumption. As it is an intermediate step before the projection of the SAFRAN reanalysis, the figures showing the projected maps of minimum and maximum temperature and precipitation (yearly average on the 2041-2060 horizon) and the deviation from the historical runs are displayed as Supplementary Materials (Figures S3 and S4, respectively). The main conclusion drawn from these two figures can be summarized as follows: (1) the whole catchment will face a temperature rise affecting both the minimum and maximum temperatures, reaching +3.5 °C in the mountainous region of the highest elevation for RCP4.5, and +4.5 °C for RCP8.5 in the same areas; (2) precipitation is also expected to decline throughout the basin. This decrease is more pronounced on the reliefs than on the plains (Table 5). This drop will reach −69.33 mm/year (−31%) on the plain against −80.3 mm/year (−32%) in the Atlas mountains for RCP8.5.

Projection of SAFRAN Reanalysis
Finally, the delta change determined between the historical and the future projections of Euro-CORDEX both corrected using the Q-Q functions are applied to the SAFRAN reanalysis using the nearest Euro-CORDEX grid points. The yearly maps for the 2041-2060 horizon and the two scenarios are displayed in Figure 8. Table 5 shows the annual changes. The SAFRAN projection is consistent, on average, with the Euro-CORDEX runs, with changes of the same order of value and higher warming projected for the mountains than for the plains. The only noticeable exception is the decrease in precipitation in the mountainous areas (elevation > 1000 m) for the RCP8.5 scenario, which is much lower for the SAFRAN projection than for the Euro-CORDEX prediction. Indeed, the delta change method computes the relative change between the historical data and the future runs. Thanks to the higher grid cell resolution, SAFRAN is better than Euro-CORDEX in mapping the strong contrast precipitation amount between the plain and the mountain: the average annual precipitation estimated by Euro-CORDEX in the mountains and the plain is equal to 275.34 mm/year and 223.48 mm/year, respectively, while it is equal to 353.0 mm/year and 221.54 mm/year for SAFRAN. Finer details of the change can be observed with the high-resolution SAFRAN projection by comparison to the 12 km Euro-CORDEX run.

Summary and Conclusions
The SAFRAN system of Météo-France was first assessed on the Tensift catchment thanks to station data using a leave-one-out approach, observations on snow cover both in situ and remotely sensed, and by comparison to a more data-driven approach for climate data spatialization, namely the MICROMET model [1] on a mountainous subcatchment. The main results can be summarized as follows: (1) due to complex orography, a variable mesh grid ranging from 8 km to 1 km based on terrain elevation outperformed the regular grid usually implemented in SAFRAN for temperature, wind speed, precipitation and, to a lesser extent, incoming solar radiation; (2) the leave-one-out approach is applied to a plain and a mountain station with significantly better performance for all variables obtained in the plain; (3) the timing (start and end) of the snow season is properly reproduced both by MICROMET and SAFRAN in comparison to the MODIS daily snow cover area product, although some large discrepancies are observed around the snow-covered area peak; in addition, the comparison of the predicted snow amount by SAFRAN and MICROMET to data of the Oukaimeden station supports a good timing prediction but also exhibits a large underestimation by both models with slightly better results obtained with SAFRAN.
Finally, the projection of SAFRAN reanalysis was carried out using a three-step approach. It consists in applying the quantile-quantile disaggregation to a station in the plain and a station in the mountain. The obtained transfer functions are then used to downscale the Euro-CORDEX scenario by separating the plain and mountainous grid cells using a 1000 m threshold on elevation. Finally, the delta change coefficients between the historical runs and the future scenarios of the disaggregated Euro-CORDEX are computed and applied to SAFRAN using the nearest Euro-CORDEX grid point for RCP4.5 and RCP8.5 at the 2041-2060 horizons. The obtained results are consistent with the disaggregated Euro-CORDEX scenarios and showed that the mountainous areas are expected to face a higher increase in air temperatures than the plains (up to +2.5 °C for RCP8.5 2041-2060 and +1.71 °C for RCP4.5 -same horizon-for SAFRAN), in agreement with the literature. This high-resolution historical and future climate forcing, available for the first time in this region, is to be used for impact studies, in particular those concerning water resources. It will allow for more localized impact studies at the subcatchment scale as required by water managers to anticipate their future catchment planning.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/, Figure S1: Euro-CORDEX historical data from the 4 RCM models on two stations: Marrakech (left) and Oukaimeden_CAF (right), before the application of the bias correction method (the Quantile-Quantile method), Figure S2: Figure S2. Comparison of the evolution of the reanalysis of two types of grids (Regular of 8 km and Irregular) from the SAFRAN model with the observation at the CHICHAOUA station, for the period 2004-2014., Table S1: Some analysis systems used in the world., Figure S3. Spatialization map of Euro-CORDEX 12Km futuristic projections on the Tensift Basin, according to the two scenarios RCP4.5 and 8.5 for 2041-2060., Table S2. Validation result of the two grids (regular 8 km and irregular) of the SAFRAN model, averaged over the 10-year period (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014), at the CHICHAOUA station., Figure S4. Deviation map of the Euro-CORDEX futuristic projections relative to the two scenarios (RCP4.5 and RCP8.5) compared to the control period used.