Future projections of High Atlas snowpack and runoff under climate change

The High Atlas, culminating at more than 4000 m, is the water tower of Morocco. While plains receive less than 400 mm of precipitation in an average year, the mountains can get twice as much, often in the form of snow between November and March. Snowmelt thus accounts for a large fraction of the river discharge in the region, and is particularly critical during spring, as the wet season ends but the need for irrigation increases. In the same region, future climate change projections point towards a significant decline in precipitation and enhanced warming of temperature. Understanding how the High Atlas 5 snowpack will evolve under such trends is therefore of paramount importance to make informed projections of future water availability in Morocco. Here, we build on previous research results on snow and climate modeling in the High Atlas to make detailed projections of snowpack and river flow response to climate change in this region. Using a distributed energy balance snow model based on SNOW-17, high-resolution climate simulations over Morocco, and a panel regression framework to relate runoff ratios to regional meteorological conditions, we quantify the severe declines in snowpack and river discharge that are 10 to be expected, even under a scenario of substantial mitigation of emissions. Our results have important implications for water resources planning and sustainability of agriculture in this already water-stressed region.

mountains to the south typically receive up to 800mm (Ouatiki et al., 2017). Most of that precipitation occurs between October and May, when the region is under the influence of North Atlantic westerlies (Knippertz et al., 2003;Tuel and Eltahir, 2018) ( Fig. 1-b). Mountainous areas also experience substantial precipitation during summer, due to small-scale convection (Born et al., 2008). As a consequence, vegetation outside the lower-elevation valleys is sparse, essentially limited to bare soil, grass and occasional shrubs (Baba et al. 2019). Snowfall is common between November and March above 1500m elevation, and a 70 somewhat persistent snowpack is not uncommon above 2500m (Marchane et al., 2015). Inter-annual variability is substantial, however, following that of precipitation (Boudhar et al., 2010). Melt is rapid, beginning in March and typically lasting 1-2 months at most (Tuel et al., 2020a).

Hydroclimatological data
We use for this study a mixture of model-, station-and satellite-based hydrometeorological data. Model-based data is described 75 in section 2.3. Daily precipitation data are available at seven stations in the study area, including three at more than 1200m elevation, over the 1980-2015 period ( Fig. 1-a). For each station, we discard the months for which more than 10% of the data is missing. Daily discharge measurements are available at seven locations as well, between 1978 and 2015. Each has at most 0.5% of missing data. These locations define seven sub-catchments for which runoff will be modeled, from north to south: Tarhat, Chacha, Ouchene, Tillouguite, Moulay Hassan, Segmine and Tamesmate ( Fig. 1-a). These sub-catchments include most of 80 the area within the Oum-Er-Rbia watershed that receives significant snowfall (Marchane et al., 2015). Their average elevation varies from 1460 to 2360m. We remove the contribution from base flow by subtracting the minimum monthly discharge value for each catchment and each hydrological year (September-August). This correction is minor for all catchments except Tarhat, the northernmost one, which includes the headwaters of the Oum-Er-Rbia river, and receives a substantial contribution of base flow to its annual discharge. In particular, the flow at Tarhat remains high during summer (≈35% of its wet-season peak), likely 85 due to groundwater discharge from deep mountain aquifers. Annual cycles of corrected monthly discharge are shown on Fig.

2-a.
Satellite-based data is used for basin-wide precipitation, temperature and snow cover. 3-hourly precipitation from the TRMM TMPA (TRMM Multi-Satellite Precipitation Analysis) 3B42 version 7 dataset is used as the reference precipitation dataset for the region. It consists in remotely-sensed data corrected with rain gauge data on a monthly basis (Huffman et al., 2007). The data cover the period 1998 to present. While satellite-based precipitation data suffers from numerous biases (Milewski et al., 2015;Derin et al., 2016;Hashemi et al., 2017), it is often the only option available in complex terrain where stations are scarce. Milewski et al. (2015) and Ouatiki et al. (2017) assessed the accuracy of the TRMM 3B42 V7 dataset in the Oum-Er-Rbia and found that, although unreliable at the daily timescale, it offered satisfactory estimates of precipitation if averaged in space or time. Annual cycles of TRMM precipitation for the seven catchments are shown on Fig.1 -b. 95 For comparison, we also consider the CHIRPS dataset, available from 1981 onwards at a 0.05 • resolution (Funk 2015). CHIRPS is produced by combining high-resolution satellite-based precipitation with station and fine-scale topography data. In our region of focus, TRMM and CHIRPS show some differences ( Fig. 3-a,b): CHIRPS is notably wetter, particularly near the Tizi N'isly station. A comparison of monthly values with the four available stations above 1000m suggests a rather dry bias in TRMM and inconsistent biases in CHIRPS d). Both datasets have a dry bias over the north-eastern corner of our domain (around 100 the Kenifra station). Across the four stations shown on Fig. 3, absolute biases range from 13.4mm to 25.6mm for TRMM and from 18.6mm to 25.2mm for CHIRPS. Our purpose here is not to perform an in-depth comparison of the performance of the two datasets, but to test the robustness of the runoff projection results to variability in reference precipitation. To bias-correct regional climate model output for snow modeling (Section 3.1), however, we follow T20a and consider only the TRMM dataset.
Reference surface air temperature is derived from MODIS Land Surface Temperature (LST) product MOD11A1 L3 version 6 105 at 1km resolution (Wan, Z., Hook, S., Hulley, 2015). We refer to T20a for details on data filling and correction. Observed snow cover area for the region is extracted from the MODIS Terra snow cover daily L3 product (MOD10A1) at 500m resolution (Hall and Riggs, 2016). Snow cover is detected using values of the Normalized-Difference Snow Index, based on reflectances in the visible/near infrared and middle infrared. We apply the correction methodology described in Marchane et al. (2015) which allows to substantially reduce the number of missing data points due, mainly, to cloud cover, and average the data at 110 the weekly timescale as in T20a. We refer to these two studies for a discussion of the accuracy of the MODIS dataset in this region. All MODIS data is available from February 2000 to present. Elevation data is taken from the Shuttle Radar Topography Mission 90-meter resolution dataset version 4.1 (STRM90) (Jarvis et al., 2008), and interpolated to the approximately 1km resolution of the MODIS land surface temperature data.
2.3 Regional climate simulations 115 We consider the regional climate downscaling data and projections developed by Tuel et a. (2020b) for the Western Mediterranean, at a 12km resolution, using the MIT Regional Climate Model (MRCM). MRCM is based on the Abdus Salam International Centre for Theoretical Physics Regional Climate Model Version 3 (RegCM3) (Pal et al., 2007), but with significant enhancements of model physics, and notably a coupling with the Integrated BIosphere Simulator land surface scheme (IBIS).
Dynamical downscaling is performed for ERA-Interim (1982-2011 (Dee et al., 2011) as well as three carefully-selected GCMs 6-hourly wind speed, specific humidity, air temperature, precipitation, and downward longwave and shortwave are extracted from the MRCM output over our domain. For all three GCM-driven simulations, as well as the ERA-Interim driven run (here-125 after referred to as ERA/MRCM), air temperature and precipitation data are statistically downscaled and bias-corrected at the 6-hourly timescale using MODIS LST-derived air temperature and TRMM precipitation at their native resolutions as respective targets, via the CDF-transform method (Michelangeli et al., 2009). Alone among the three GCMs, the IPSL-CM5A-LR model exhibits a negative bias in wet days that we correct at each grid cell by randomly generating wet days of magnitude drawn from the corresponding distribution of wet-day precipitation in the TRMM dataset. For bias correction, reference periods for All bias corrections are performed for the cold (November-April) and warm (May-October) seasons separately.
Additionally, we use wind speed, downward long-and shortwave radiation and specific humidity from the ERA/MRCM simulation over the 1982-2005 period as reference, since no observations are available. The corresponding variables in each GCM-135 driven simulation are therefore bias-corrected using the ERA/MRCM data as target. Specific humidity is further downscaled to the 1km MODIS LST resolution based on an empirical lapse-rate µ estimated at each time step: where q 12 is the specific humidity in a given 12-km resolution grid cell of elevation z 12 , and q the downscaled value at elevation z.

Methods
The methodological framework adopted in this study is summarized in Fig. 4. We start from GCM simulations from the CMIP5 archive, dynamically-downscaled with MRCM. We apply bias-correction to the MRCM output, which we then feed into a 1km-resolution snow model over our study region to reconstruct snowpack under the various emissions scenarios. Finally, a statistical model is developed for catchment runoff coefficients (RCs), in order to make projections of runoff under future 145 climate conditions.

Snow model
We apply the SNOW-17 model (Anderson, 2006) with a radiation-derived temperature index for melt (Follum et al., 2015) as We run the model at a 6-hourly time step and at the native MODIS LST 1km resolution over the same 13104 km 2 domain as T20a, that encompasses the seven catchments shown on Fig. 1-a. Elevation in this domain ranges from 621m to 3890m, with 155 an average of 1882m. SWE given by the model is then translated into snow cover fraction for each grid cell using the following relationship: with k = 100. The snow model requires optimizing three parameters: M f (melt factor), N M F max (maximum negative melt factor) and T IP M (coefficient used in updating snowpack temperature) (Anderson, 2006). In keeping with T20a, parameter 160 calibration is performed by maximizing the Nash-Sutcliffe coefficient (Nash and Sutcliffe, 1970)

Statistical modeling of runoff coefficients
We model catchment runoff coefficients (RCs), defined as total October-May discharge divided by total October-May precipitation, across time and space as functions of large-scale hydrological variables by adopting a panel regression framework.
Panel regression allows to enhance the effective size of a dataset and to obtain more robust estimates of the response to se-170 lected covariates compared to more traditional regression approaches (Steinschneider et al., 2013;Davenport et al., 2020). It also allows to account for static (space-dependent) and time-varying (time-dependent) factors, although here, with only seven catchments, we do not have enough data to make robust statements about static factors responsible for the disparity in average RC ( Fig. 2-b). Therefore, we focus on time-varying covariates, and consider a fixed-effects formulation: where j ∈ {1, ..., 7} is the catchment index, t is the time index, i is the covariate index, and RC j represents time-invariant, watershed-specific fixed effects (drainage area, land cover, mean climate), X i j,t are the covariates, β i are regression coefficients and j,t is random noise. For covariates, we consider catchment-averaged October-May precipitation (P), relative humidity (RH), temperature (T), snow water equivalent (SWE) and snow fraction of precipitation (SF). Enhanced precipitation or relative humidity lead to wetter soils and can be expected to be associated with higher RC values. Similarly, higher temperatures 180 increase evapotranspiration and tend to decrease runoff. Finally, increased snow cover favors losses by sublimation and shifts the distribution of runoff regimes towards slower runoff as opposed to rapid overland flow following rain storms. However, larger snow cover may also increase the risk of rain-over-ice events, which tend to have very high runoff coefficients (Davenport et al., 2020).
Therefore, we use the ERA/MRCM precipitation data, bias-corrected with TRMM, to estimate catchment-scale precipitation for the longer 1982-2011 period, during which runoff data are available for all catchments. This allows us to calculate 190 catchment runoff coefficients (RCs), defined as observed total October-June discharge divided by estimated total October-June precipitation. To assess the robustness of the results, we also calculate RCs based on CHIRPS precipitation (available from 1981-present).
Model selection is performed by stepwise regression: starting from a model with no covariates, covariates are added one at a time in the order of highest improvement to model skill, as determined by its Akaike Information Criterion (AIC). At each step, 195 we also test whether removing any of the currently selected variables and replacing it by one of the remaining, non-selected ones, brings any improvement. To estimate the sensitivity of RC to changes in climate conditions, we modify covariate values in the 1982-2011 ERA-Interim downscaled simulation by adding projected long-term changes in the GCM-driven simulations: where m ∈ {1, 2, 3} is model index, I is the set of optimal covariates, X i j,t are ERA-Interim downscaled covariate values and 200 X i m represent long-term covariate changes drawn at random according to: where µ i m,s (respectively σ i m,s is the average (respectively standard deviation) of covariate i in model m and scenario s. Results for the three models are then pooled together to yield a future distribution for RC j,t .  December-March average snow cover over the whole modeled area reaches 1460 km2 in MODIS observations but only 1275 km 2 in the ERA-Interim driven run and 1185 km 2 in the ensemble mean historical average -a bias largely concentrated at low elevations. Additionally, despite the statistical downscaling of the MRCM output to 1km, elevation-driven gradients in snow cover are also less sharp in the MRCM experiments compared to observations ( Fig. 6-a,b). All results tend to remain however in a narrow band around the MODIS values and, except for elevations below 2000m, inter-model spread generally 215 covers observed snow cover values. Inter-annual variability in basin-wide snow cover is lower in the simulations compared to MODIS (standard deviation of 220-440 km 2 compared to 480 km 2 in MODIS) but the discrepancy is mainly due to the negative bias at low elevations, where snow plays a much more limited role in the overall water balance (T20a). Unsurprisingly, future projections show a stark decline in snow cover across all the region (Figs. 6, 7). The greatest relative decline is at low elevations, as expected since they are already seldom above the 0 • C line (Boudhar et al., 2016). Above 2500m, projections still exhibit a Whatever the historical disparity in temperatures or precipitation between models, all models agree on the quasi-absence of snowpack under business-as-usual. This is consistent with observations that in mountain regions, above-freezing temperatures are common, even at high altitudes, where average winter temperatures are not far from freezing (López-Moreno et al., 2017).

Results and discussion
Therefore, under the RCP8.5 scenario, the projected 4 • C warming (Tuel et al. 2020b) would regularly bring all areas but the very highest peaks well above freezing, and prevent seasonal snowpack accumulation. Annual relative sublimation losses are strongly linked to annual-mean precipitation ( Fig. 11-a). Losses from latent heat fluxes are much smaller during wet years as compared to dry years, a relationship robust across all experiments. Wet years indeed bring higher RH over the region, due to enhanced moisture advection from the Atlantic which more than compensates for the 240 larger heat advection and increased air temperatures that also occur in parallel (Knippertz et al., 2003). In addition to higher RH limiting evaporation and increasing soil moisture, warmer temperatures in wetter years also tend to decrease the snowto-precipitation ratio, thus leading to reduced relative sublimation losses. Consequently, we expect the runoff coefficient to be higher in wet years. Under future climate conditions, average relative humidity will decline by 3-6% (Tuel et al., 2020b), and thus sublimation rates will be higher when snow is present. However, because the snow-to-precipitation ratio will also sharply 245 8 https://doi.org/10.5194/hess-2020-622 Preprint. Discussion started: 18 December 2020 c Author(s) 2020. CC BY 4.0 License. decline due to rising temperatures, the overall loss of annual precipitation by sublimation will tend to decrease by about a third ( Fig. 11-b).
For the areas as a whole, decreasing precipitation and warmer temperatures are the primary causes of the projected decline in snowpack. In addition, drier soils together with enhanced absorption of solar energy where snowpack disappears will also lead to enhanced warming locally, driving yet further snowpack melt. We do not explicitly take this into account in our model. In 250 particular, the snow albedo effect is largely absent from the MRCM simulations due to their 12km resolution, which is still too coarse to represent the complex topography. For areas at the highest elevations (near 4000m) which may still remain largely below freezing in future winters, melt may not increase very significantly in the middle of winter; however, a drier atmosphere will still be associated with reduced precipitation and increased sublimation losses, which will play a critical role in reducing the snowpack.

Runoff modeling and projections
The panel regression and model selection framework are applied to runoff coefficients and selected covariates over the 1982-2011 period. Stepwise regression yields as optimal covariates relative humidity and snow fraction. The adjusted r 2 is equal to 0.30, meaning that these two covariates explain a small third of inter-annual variability in RC. Fitted RCs against observations are shown on Fig. 12. Consistent with the model r 2 , fitted values have a much smaller variance. Still, for all the catchments, 260 except Segmine, we observe a significant positive relationship between fitted and observed values. The coefficients for RH and snow fraction are both significant, with β RH > 0 and β SF < 0 (Table 1). All else being equal, a larger RH yields more runoff -consistent with previous studies (Wang et al., 2016;Duan et al., 2017) -and more of the precipitation falling as snow yields less runoff. These results are robust to the choice of the precipitation dataset. When using the CHIRPS dataset, average RCs are generally lower, due to CHIRPS?s wet bias compared to TRMM (Fig. 3); the optimal RC model includes RH and SF as 265 well, but also precipitation as a third variable (Table 1). The values of β RH and β SF obtained with CHIRPS data are similar, although β SF is slightly less significant.
For both precipitation datasets, the effect of RH on runoff coefficient, as measured by the respective regression coefficient β RH , is 3-5 times that of snow. As discussed above, we can understand the influence of snow fraction by noting that a higher snow fraction means more opportunity for sublimation, particularly large at high elevations, and evaporation of melted snow, 270 consistent with our analysis of sublimation losses (Fig. 11). Precipitation in the area tend to occur in short and intense storms, and quickly saturate the dry soil, leading to rapid overland flow with limited opportunity for evaporation (El Khalki et al., 2018). By comparison, snowmelt is slow and leads to a more gradual surface flow with the potential for higher evaporative losses in this climate where evaporation tends to be water-limited. More winter snow also leads to a higher likelihood of rainon-snow episodes in spring, known to cause rapid flooding due to high runoff efficiencies (Davenport et al., 2020). More snow 275 may also mean more opportunity for infiltration and aquifer recharge (Hssaisoune et al., 2020); however, in six out of seven catchments considered here, base flow is negligible. Still, aquifer discharge may occur naturally further down the mountains, or artificially via direct groundwater pumping in the agricultural plains. Compared to surface runoff, groundwater remains a small fraction of water use in the Oum-Er-Rbia basin (<15%), and even less of available renewable water since aquifers are largely overdrawn (Hssaisoune et al., 2020). Groundwater data are quite scarce in this region; we make the choice of purely focusing on surface runoff, keeping in mind that a more complete picture of basin-wide water availability would also require taking aquifer fluxes into account.
Future runoff projections are characterized by consistent, steep declines in runoff coefficients of 5-17% under RCP4.5 and 15-30% under RCP8.5 ( Fig. 13-a). The impact of decreasing RH largely dominates over that of declining snow fraction. Tarhat and Chacha, the two watersheds which already receive almost no snow in the present climate, exhibit the greatest relative 285 RC decline, whereas in other watersheds, decreases in snow fraction help limit the decline in RC (Table 2). Combining now precipitation and runoff coefficient estimates, we find a 20-40% decrease in runoff in the RCP4.5 experiments, compared to a 50-65% decline under RCP8.5 ( Fig. 13-b, Table 2). Decreases in precipitation drive most of the runoff trends, especially in the RCP8.5 scenario. Projected RC declines are about the same when using CHIRPS data for the analysis, although tend to be slightly higher, on average by 3-6% on average (not shown). Our projected runoff trends are consistent with those of El Moçayd 290 et al. (2020) who focused on catchments part of the Sebou watershed, where snow plays a much smaller role than in our region (Marchane et al., 2015). This is not surprising given the weak impact of the declining snow fraction on runoff coefficients, compared to that of the relative humidity decline. As a final note, one can ask whether observed runoff trends are consistent with these projections. At the country scale, government data indicates that river discharge has indeed significantly declined over the last 60 to 70 years, at a rate of about 5% per decade (Fig. A1). The wet years of the 1960s tend however to bias the 295 result towards a steeper decline than expected. In this case, natural decadal variability linked to the North Atlantic Oscillation acted in the same direction as the expected climate change response. Runoff data from the seven sub-catchments studied here are less clear. They cover shorter time periods (by about 25 years) and, more importantly, begin in the late 1970s, just when the climate turned much drier in Morocco; drier-than-average conditions indeed prevailed from 1980 to the mid-1990s, associated with increased frequency in the positive phase of the North Atlantic Oscillation. Thus, runoff trends tend to be slightly positive 300 (about +1%/yr), though none are significant, even at the 10% level.

Conclusions
Based on the robust understanding of its snow water balance in the current climate, we quantified in this final chapter the response of the High Atlas snowpack to climate change using high-resolution downscaled climate projections. Unsurprisingly, given the warming and drying trends projected by climate models for this region, we find that the High Atlas snowpack will 305 significantly decline, even with substantial mitigation of emissions. By the end of the century, snow may become a rarity below 2000m, and even near the highest peaks, snowpack water equivalent could decline by 80%. Precipitation decreases of 40-60% are responsible for much of these trends, with background warming accounting for the rest.
The analysis of runoff coefficients for seven mountain catchments showed that a third of their inter-annual variability could be explained by large-scale meteorological factors like snow fraction of precipitation and relative humidity. Interestingly, in 310 this region, a larger snow fraction leads to less runoff. While the reverse is believed to be true at higher latitudes (Berghuijs et al., 2014), this finding is consistent with other analyses in warm, semi-arid regions that receive substantial amounts of snow