Snow Lapse Rate Changes in the Atlas Mountain in Morocco Based on MODIS Time Series during the Period 2000–2016

: The spatio-temporal distribution of snow cover metrics in a mountainous area is mainly related to the climatic conditions as well as to the prevailing morphological conditions. The present study aimed to investigate the altitudinal sensitivity of snow cover metrics using the MODIS Terra snow cover product (MOD10A1 v5). Annual snow metrics, including start of snow season (SOSS), end of snow season (EOSS), and snow cover duration (SCD) were extracted from snow-covered area (SCA) maps, which had been pre-processed using a cloud removal algorithm; the maps were of the Atlas Mountains, taken from the period of 2001–2016. In addition, a linear regression was applied to derive an annual altitudinal gradient for each snow metric in relation to various spatial scales in order to analyze the interdependency between snow and topography, and especially to assess the potential temporal trend of the snow gradient. Results indicated that elevation was the principal regulator of snow presence where snow was mostly accumulated above 2500 m. The annual altitudinal gradients for EOSS and SCD showed a marked negative trend beginning in 2007. However, the SOSS altitudinal gradient was marked by a positive trend. The mean SCD gradient for the entire Atlas Mountains decreased from 6 days/100 m to 3 days/100 m. This is a new and important ﬁnding since it may indicate the impact of climate change on the dynamics of snow metrics and provides guidance for water managers to better manage the snowmelt water with different terrain features.


Introduction
Meltwater issued from a snow-covered area is a significant component of the water balance in many of the world's catchments [1]. Snow accumulation during the cold season presents a major source of fresh water during the melt period, particularly in a mountainous region. Consequently, it also provides water for irrigation, aquifer recharges, and contributes to fill dams, which are important for agriculture and also hydropower generation [2][3][4].
In Morocco, a place dominated by an arid and semi-arid climate, the mountainous regions (e.g., High and Middle Atlas) receive a significant amount of precipitation as snowfall during winter. As a result, the runoff regime in the high-elevation catchments is highly dependent on snow cover, and snowmelt contribution is roughly estimated to be between 15% and 50% of the total annual discharge in the Tensift river sub-basins [5].
at large scale [23,29]; (3) monitoring snow cover variability at a large scale from remotely sensed data [6,7,30,31]; (4) assessing snow component contribution through hydrological models [5,32]; and (5) linkages between climate and snow-covered area (SCA) derived from remote sensing observations [33]. Despite all these investigations, there has not been any systematic study on the effects of topographic characteristics on snow cover seasonality over the Atlas Mountains.
The aims of this paper are: (1) to assess the inter-annual variations in snow cover parameters (i.e., SOSS, EOSS, and SCD) over the whole Atlas Mountain during the period of 2001-2016; (2) to analyze topographic factors, including elevation, aspect, and slope, in snow cover variability; (3) to estimate the altitudinal gradient using long-term trends of the selected snow metrics.

Study Site
This work covers the mountainous region of the northern part of the Moroccan Kingdom, stretching from 2 • W to 10 • W longitude, and from 30 • N to 35 • N latitude. In this study, the mountain domain is defined as the area above 1000 m ( Figure 1). The reason for this delimitation is to only account for regularly covered snow areas. The study domain covers an area of 121,000 km 2 , with 54.42% of the elevation range located between 1000 and 1500 m, 44.58% of the elevation band lying between 1500 and 3000 m, and with only 1% of the area located above 3000 m ( Table 1). The study area is composed of seven main catchments in Morocco (Moulouya, Oum Er rbii, Sebou, Ziz Rheris, Tensift, Souss, and Draa).     [34]. The climate in this region is governed by the complex influence of the Atlantic Ocean to the west, the Mediterranean Sea to the north, and the Sahara Desert to the south. The aspect of the study zone is predominantly northerly (29%), while the other exposures (West, East, and South) vary between 20% and 27% ( Table 2).

MODIS Daily Fractional Snow Product (MOD10A1)
We used data derived from version 5 of MODIS/Terra Snow Cover Daily (MOD10A1) at 500 m resolution [35]. The snow mapping technique used by the MOD10A1 data is based on an algorithm called SNOWMAP, which determines the occurrence of snow using the Normalized Differential Snow Index (NDSI) combined with additional criteria [35][36][37].
The MOD10A1 represents the fraction of snow cover (percentage of snow cover of the observed pixel). This snow cover fraction is calculated as a function of NDSI using the following equation [37]: The FSC dataset was pre-processed using a three-step cloud-correction method and evaluated over the Moroccan mountains by Marchane and others [6]. The post-processed dataset showed a drastic reduction in the number of cloud-covered pixels, from 22.6% (raw data) to 0.8% (post-filtered data), estimated for the period between 2001 and 2013.
The assessment of these adjusted data products by comparing them to the automatic snow depth measurements of 5 stations and to a unique satellite dataset from the formosat-2 sensor at 8 m resolution showed that the accuracy of this product was at 89% [6].
The daily MODIS product has shown good capabilities in the production of timing characteristics (onset, melt-out dates, and snow cover duration) despite the difficult detection conditions in the Moroccan mountains [6].

Extraction of Snow Cover Metrics
The snow parameters (SOSS, EOSS, SCD) were extracted from the snow cover maps for 16 snow seasons (from September 2001 to June 2016). They served to estimate the possible trends in the snow metrics gradient and to analyze the relationship with topographic variables (Figure 2).

Topographic Effects
The extraction of topographic characteristics was completed using the Digital Elevation Model (DEM), derived from NASA's Shuttle Radar Topography Mission (SRTM) images. The SRTM data were used for its proven accuracy in complex terrains, as reported by many authors [39,40]. This DEM was resampled at a 500 m spatial resolution to make it consistent with the MODIS data. Figure 3 shows the inter-annual variation of snow-covered areas across the entire Atlas Mountains for altitudes higher than 1000 m. The 10-day fractional snow-covered area (fSCA) showed a high seasonal variability, which differs greatly by elevation. The timing of the snow season can be considered hypsometrically by elevation bands, spatially by pixel, or globally (study area) [38], which gives multiple metrics. The start of snow season (SOSS) was calculated as the date when a pixel exceeded a threshold of fractional snow cover (FSC). The end of snow season (EOSS) was determined when data remained below the threshold for at least 20 days. In this case, the minimum threshold of FSC was set to 50%, consistent with the approach of [26]. The length of the SCD was calculated as the difference in days between the SOSS and the EOSS. A similar approach was adopted to detect timing indicators from the SCA at the study area level, except that the threshold was chosen to be 5% of the maximum snow-covered area. Calculating the snow metrics for each hydrological year from 2001 to 2016 also allowed for the calculation of the coefficient of variation (CV) as follows:

Spatial and Temporal Variation of Snow Cover
where σ is the standard deviation of metrics and µ is the mean over 16 years. CV is used here as an efficient way to calculate the relative dispersion of data on the mean [28].

Topographic Effects
The extraction of topographic characteristics was completed using the Digital Elevation Model (DEM), derived from NASA's Shuttle Radar Topography Mission (SRTM) images. The SRTM data were used for its proven accuracy in complex terrains, as reported by many authors [39,40]. This DEM was resampled at a 500 m spatial resolution to make it consistent with the MODIS data. Figure 3 shows the inter-annual variation of snow-covered areas across the entire Atlas Mountains for altitudes higher than 1000 m. The 10-day fractional snow-covered area (fSCA) showed a high seasonal variability, which differs greatly by elevation.

Spatial and Temporal Variation of Snow Cover
it consistent with the MODIS data. Figure 3 shows the inter-annual variation of snow-covered areas across the entire Atlas Mountains for altitudes higher than 1000 m. The 10-day fractional snow-covered area (fSCA) showed a high seasonal variability, which differs greatly by elevation. Generally, snow season starts effectively at the beginning of November and finishes in early April, except for altitudes superior to 3500 m, where snow can occur from October and persists until May ( Figure 4). As expected, the maximum snow cover is registered during the winter months (December, January, and February). The evolution of fSCA is elevation-dependent (Figure 4), while its behavior is different. For the 1000-1500 m altitude range, the snow is almost absent in winter and the fSCA median values are lower than 0.1%. In the altitude range of 1500-2000 m, the fSCA maximum median value is around 2%, with a delay in terms of snow season start (end of December). From 2000 m, the behavior of fSCA is similar to that described for the entire study area, with a maximum median value above 10% and a maximum value of fSCA at 40%. However, the median value in the altitudes above 3000 m is more than 50% during the winter season.  Generally, snow season starts effectively at the beginning of November and finishes in early April, except for altitudes superior to 3500 m, where snow can occur from October and persists until May ( Figure 4). As expected, the maximum snow cover is registered during the winter months (December, January, and February). The evolution of fSCA is elevation-dependent (Figure 4), while its behavior is different. For the 1000-1500 m altitude range, the snow is almost absent in winter and the fSCA median values are lower than 0.1%. In the altitude range of 1500-2000 m, the fSCA maximum median value is around 2%, with a delay in terms of snow season start (end of December). From 2000 m, the behavior of fSCA is similar to that described for the entire study area, with a maximum median value above 10% and a maximum value of fSCA at 40%. However, the median value in the altitudes above 3000 m is more than 50% during the winter season.

Spatial and Temporal Variation of Snow Cover
At altitudes between 2000 and 3000 m, the snow season starts at the beginning of November, while it starts a little early for altitudes above 3000 m (mid-October). The end of the snow season is also variable according to altitude range, with early dates for the 1000-1500 m altitude range (end of January) and later dates for altitudes above 3500 m (end of May). To characterize the spatial distributions of seasonal snow cover, the main snow metrics (SCD, SOSS, EOSS) and their coefficient of variation (CV) are averaged during the 16 years over the Atlas Mountain and its seven river basins ( Figure 5). The SCD map shows that the occurrence of snow increases with altitude and tends to follow the contours of the terrain, with the highest SCD concentrated on the mountain's top (Figure 5a). The highest At altitudes between 2000 and 3000 m, the snow season starts at the beginning of November, while it starts a little early for altitudes above 3000 m (mid-October). The end of the snow season is also variable according to altitude range, with early dates for the 1000-1500 m altitude range (end of January) and later dates for altitudes above 3500 m (end of May).
To characterize the spatial distributions of seasonal snow cover, the main snow metrics (SCD, SOSS, EOSS) and their coefficient of variation (CV) are averaged during the 16 years over the Atlas Mountain and its seven river basins ( Figure 5). The SCD map shows that the occurrence of snow increases with altitude and tends to follow the contours of the terrain, with the highest SCD concentrated on the mountain's top ( Figure 5a). The highest SCD values for more than 120 days were detected in the Moulouya, Tensift, and Oum'Erbia basins where the elevation exceeds 3500 m. variability associated with it. Ultimately, the control of elevation on the SCD is apparent in the Atlas Mountains. Snow seasons start earlier at high altitudes. Meanwhile, first snow at the foothills comes later during the snow season ( Figure 5b). The average of SOSS over the study area is spread over two months, with the earliest starting around the first of November. The latest is recorded at low altitudes of less than 2000 m around the first of January.
A distribution map of the CV-SOSS (Figure 5e) highlights the degree of spatio-temporal variability associated with SOSS, with values ranging from 18% to 40%. High CV-SOSS values are recorded at low altitudes, which are marked by high variability in SOSS dates. In contrast, at high altitudes, CV-SOSS values are relatively lower since these areas record less variable SOSS dates. The spatial variation of the average melt dates over the period of 2001-2016 is shown in Figure 5c; the earliest dates are recorded at lower elevations, while the latest dates occurred at the highest elevations. Generally, the melting dates vary from the beginning of March to the end of May. Concerning CV-EOSS, its relationship with elevation is the same as that between CV-SOSS and elevation.

Global Altitudinal Effect on Snow Cover Metrics
In previous studies, the snow extent in the Atlas Mountains was mainly associated with elevation [6,30]. The specific concern here is to quantify the rate of snow metrics variations against altitudes during the 16 years studied over the whole Atlas bands. We also focused on analyzing the change of the altitudinal gradient of the snow metrics in relation to both space and time scales. Figure 6a shows the average of snow cover metrics based on elevation zones, with an elevation step of 200 m. As we have visually noticed from the SCD map, the increase in the duration of the snow cover season (SCD) seems to be consistent with the increase in altitude. The control of the elevation on the SCD is very clear from altitudes above 2500 m. SCD is positively and linearly correlated with elevation. The main altitudinal gradient of SCD is about 4 days/100 m. However, SOSS is negatively correlated with elevation, but with strong variability at altitudes below 3000 m.
The EOSS is positively correlated with elevation. Moreover, a significant relationship occurs at elevations above 2500 m. The analysis of the effect of altitude on the snow cover However, the CV-SCD map (Figure 5d) highlights an opposite variation to SCD. At the high elevations, CV-SCD was generally less than 30%. The CV-SCD increases significantly at lower altitudes, indicating the relatively short SCD and the high inter-annual variability associated with it. Ultimately, the control of elevation on the SCD is apparent in the Atlas Mountains.
Snow seasons start earlier at high altitudes. Meanwhile, first snow at the foothills comes later during the snow season (Figure 5b). The average of SOSS over the study area is spread over two months, with the earliest starting around the first of November. The latest is recorded at low altitudes of less than 2000 m around the first of January.
A distribution map of the CV-SOSS (Figure 5e) highlights the degree of spatio-temporal variability associated with SOSS, with values ranging from 18% to 40%. High CV-SOSS values are recorded at low altitudes, which are marked by high variability in SOSS dates. In contrast, at high altitudes, CV-SOSS values are relatively lower since these areas record less variable SOSS dates. The spatial variation of the average melt dates over the period of 2001-2016 is shown in Figure 5c; the earliest dates are recorded at lower elevations, while the latest dates occurred at the highest elevations. Generally, the melting dates vary from the beginning of March to the end of May. Concerning CV-EOSS, its relationship with elevation is the same as that between CV-SOSS and elevation.

Global Altitudinal Effect on Snow Cover Metrics
In previous studies, the snow extent in the Atlas Mountains was mainly associated with elevation [6,30]. The specific concern here is to quantify the rate of snow metrics variations against altitudes during the 16 years studied over the whole Atlas bands. We also focused on analyzing the change of the altitudinal gradient of the snow metrics in relation to both space and time scales. Figure 6a shows the average of snow cover metrics based on elevation zones, with an elevation step of 200 m. As we have visually noticed from the SCD map, the increase in the duration of the snow cover season (SCD) seems to be consistent with the increase in altitude. The control of the elevation on the SCD is very clear from altitudes above 2500 m. SCD is positively and linearly correlated with elevation. The main altitudinal gradient of SCD is about 4 days/100 m. However, SOSS is negatively correlated with elevation, but with strong variability at altitudes below 3000 m.

Snow Cover Metrics in Relation to Aspect
The orientation aspect determines the duration of the insolation and the intensity of the received radiation [19]. In our study area, this parameter is more contrasted. The reduction of solar radiation in the shade areas helps to conserve soil moisture and reduces the evaporative power of the air. This parameter has a great influence on the snowmelt process and thus on snow distribution.
Generally, the south-facing slopes receive more solar radiation; therefore, snow ac-

Snow Cover Metrics in Relation to Aspect
The orientation aspect determines the duration of the insolation and the intensity of the received radiation [19]. In our study area, this parameter is more contrasted. The reduction of solar radiation in the shade areas helps to conserve soil moisture and reduces the evaporative power of the air. This parameter has a great influence on the snowmelt process and thus on snow distribution.
Generally, the south-facing slopes receive more solar radiation; therefore, snow accumulation and snow duration on south-facing slopes were more sensitive to temperature increment than those on north-facing slopes [19,21]. Figure 7 displays the evolution of the annual SCDs of four aspects in each hydrologic year, from 2001 to 2016.
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 16 a significant difference. In some years, we found a maximum SCD in the south-facing slopes. The variability of snow duration induced by the exposure effect is mainly lost and not captured with the MODIS resolution (500 m). This result is consistent with a recent study by [24], which indicated that MODIS is not able to capture the differences in snow cover due to the contrasting melt rates between north and south facing slopes.
The seasonality of the Atlas climate is also determined essentially by the atmospheric conditions [41]. During winter, this region is dominated by polar front perturbations whose cold air masses reach the Maghrib from continental and maritime areas, often in connection with the Azores anticyclone, or from the Eurasian continent [41]. The North Atlantic Oscillation (NAO) is the principal mode of atmospheric variability in the northern hemisphere [42], which includes Mediterranean mountainous areas. Marchane et al. [33] analyzed the linkages between snow cover and the NAO over the Atlas Mountains and suggested that the NAO is likely to strongly govern snowpack dynamics in the region. At the pixel scale, the effect of orientation on snow duration is not clear enough. From  Figure 8a, we can see that most of the pixels with SCD > 120 (days) have a southern aspect, between 120° and 180°, and an elevation above 3500 m. The role of the aspect disappears as the elevation increases. These results are also consistent with the study of [43], who explained that the control of topography is not important when the snow cover is homogeneous (at high altitudes).
The latest dates of SOSS are recorded from the aspect of 300° to 0°, and an altitude lower than 2500 m (Figure 8b). Regarding snowmelt dates, we also notice that the impact of the aspect on the distribution of EOSS dates is not very clear (Figure 8c).
Generally, based on the spatial resolution of MODIS (500 m), the aspect control on seasonal snow measurements is not clear, either at the Atlas level or even by basin. In contrast to previous studies highlighting that snow accumulates more on northfacing slopes [16,30], the control of the aspect on the SCD-extracted using the MODIS product-is not clear, either at the scale of the entire Atlas or at the scale of the basin. The results show that the snow duration on the north and south-facing slopes do not display a significant difference. In some years, we found a maximum SCD in the south-facing slopes. The variability of snow duration induced by the exposure effect is mainly lost and not captured with the MODIS resolution (500 m).
This result is consistent with a recent study by [24], which indicated that MODIS is not able to capture the differences in snow cover due to the contrasting melt rates between north and south facing slopes.
The seasonality of the Atlas climate is also determined essentially by the atmospheric conditions [41]. During winter, this region is dominated by polar front perturbations whose cold air masses reach the Maghrib from continental and maritime areas, often in connection with the Azores anticyclone, or from the Eurasian continent [41]. The North Atlantic Oscillation (NAO) is the principal mode of atmospheric variability in the northern hemisphere [42], which includes Mediterranean mountainous areas. Marchane et al. [33] analyzed the linkages between snow cover and the NAO over the Atlas Mountains and suggested that the NAO is likely to strongly govern snowpack dynamics in the region.
At the pixel scale, the effect of orientation on snow duration is not clear enough. From Figure 8a, we can see that most of the pixels with SCD >120 (days) have a southern aspect, between 120 • and 180 • , and an elevation above 3500 m. The role of the aspect disappears as the elevation increases. These results are also consistent with the study of [43], who explained that the control of topography is not important when the snow cover is homogeneous (at high altitudes). Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 16

Snow Cover Metrics in Relation to Slope
The results obtained for the duration of the snow cover confirm the data obtained for SCA, with a maximum duration at slopes > 35°. On the one hand, this is due to the fact that a steeper slope receives less radiation than a flatter surface; on the other hand, it is because a significant part of steep slopes is located at high altitudes.
The variation in the SCD based on the slope shows strong inter-annual variability ( Figure 9). The longest period is recorded in the hydrological year 2009, while the shortest period of snow cover is recorded in the year 2016. An analysis of the slope control on the snow metrics distribution in the Moroccan Atlas shows that its effect on SOSS and EOSS is not homogeneous and especially influenced by the altitude.

Snow Metrics Lapse Rate Change
The snow metrics variability in mountain regions can serve as a good indicator of climate change. In the region of the Atlas Mountains, an earlier melting could limit water availability in the neighboring plains in the summer. The very few investigations on the snow cover variability in the Atlas Mountains did not report any significant trend. In order to investigate the possible change in snow cover metrics more precisely, we analyzed the variation of the altitudinal gradient of the selected snow metrics over the study period (2001-2016). The latest dates of SOSS are recorded from the aspect of 300 • to 0 • , and an altitude lower than 2500 m (Figure 8b). Regarding snowmelt dates, we also notice that the impact of the aspect on the distribution of EOSS dates is not very clear (Figure 8c).
Generally, based on the spatial resolution of MODIS (500 m), the aspect control on seasonal snow measurements is not clear, either at the Atlas level or even by basin.

Snow Cover Metrics in Relation to Slope
The results obtained for the duration of the snow cover confirm the data obtained for SCA, with a maximum duration at slopes >35 • . On the one hand, this is due to the fact that a steeper slope receives less radiation than a flatter surface; on the other hand, it is because a significant part of steep slopes is located at high altitudes.
The variation in the SCD based on the slope shows strong inter-annual variability ( Figure 9). The longest period is recorded in the hydrological year 2009, while the shortest period of snow cover is recorded in the year 2016. An analysis of the slope control on the snow metrics distribution in the Moroccan Atlas shows that its effect on SOSS and EOSS is not homogeneous and especially influenced by the altitude.

Snow Cover Metrics in Relation to Slope
The results obtained for the duration of the snow cover confirm the data obtained for SCA, with a maximum duration at slopes > 35°. On the one hand, this is due to the fact that a steeper slope receives less radiation than a flatter surface; on the other hand, it is because a significant part of steep slopes is located at high altitudes.
The variation in the SCD based on the slope shows strong inter-annual variability ( Figure 9). The longest period is recorded in the hydrological year 2009, while the shortest period of snow cover is recorded in the year 2016. An analysis of the slope control on the snow metrics distribution in the Moroccan Atlas shows that its effect on SOSS and EOSS is not homogeneous and especially influenced by the altitude.

Snow Metrics Lapse Rate Change
The snow metrics variability in mountain regions can serve as a good indicator of climate change. In the region of the Atlas Mountains, an earlier melting could limit water availability in the neighboring plains in the summer. The very few investigations on the snow cover variability in the Atlas Mountains did not report any significant trend. In order to investigate the possible change in snow cover metrics more precisely, we analyzed the variation of the altitudinal gradient of the selected snow metrics over the study period (2001-2016).

Snow Metrics Lapse Rate Change
The snow metrics variability in mountain regions can serve as a good indicator of climate change. In the region of the Atlas Mountains, an earlier melting could limit water availability in the neighboring plains in the summer. The very few investigations on the snow cover variability in the Atlas Mountains did not report any significant trend. In order to investigate the possible change in snow cover metrics more precisely, we analyzed the variation of the altitudinal gradient of the selected snow metrics over the study period (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016).
A time series of altitudinal sensitivity of snow metrics in the Atlas Mountains are plotted in Figure 10. Generally, the snow metrics gradient exhibited large inter-annual variability. The results of linear regression analysis indicate that the altitudinal sensitivity of SCD and EOSS exhibit a decreasing trend, with slopes of −0.21 and −0.18, respectively. In 2007, the duration of snow cover increased by an average of 6 days for each 100 meters. This value was reduced to 3 days in 2016. However, the change of SOSS altitudinal gradient is marked by a positive slope of 0.14. The response of the SOSS to climate change at different altitudes is similar to that of the SCD.
A time series of altitudinal sensitivity of snow metrics in the Atlas Mountains are plotted in Figure 10. Generally, the snow metrics gradient exhibited large inter-annual variability. The results of linear regression analysis indicate that the altitudinal sensitivity of SCD and EOSS exhibit a decreasing trend, with slopes of −0.21 and −0.18, respectively. In 2007, the duration of snow cover increased by an average of 6 days for each 100 meters. This value was reduced to 3 days in 2016. However, the change of SOSS altitudinal gradient is marked by a positive slope of 0.14. The response of the SOSS to climate change at different altitudes is similar to that of the SCD.
A trend analysis shows that EOSS gradient increases below 2500 m, but decreases between 2500 and 4000 m, indicating that EOSS at high and low altitudes show a contrasting response to climate change. The analysis above further confirms that the effect of EOSS dates on the length of the SCD on the Atlas Mountains is greater than that of SOSS dates. This decrease in the SCD gradient was largely driven by an earlier snowmelt and partly by a later snow onset.
The direct reasons for this change in the snow metric gradient must be linked to the long-term change and variability of air temperature and precipitation. The proof of this connection can be provided by climate projections over Morocco that agree on robust warming and drying trends under greenhouse gas forcing. According to Marchane et al. [44], future projections in the High Atlas indicate an increase in temperature (+1.4 °C to +2.6 °C) and a decrease in total precipitation (from −22% to −31%), depending on the emissions scenario.  Table 3 shows the difference between the average of the snow metrics lapse rate over two periods of 8 years each. This comparison also shows that the SCD and EOSS gradients decreased by 0.67 day and 0.35 day, respectively, between the two periods. In contrast, the SOSS gradient shows a positive change signal of about 0.5 day.  A trend analysis shows that EOSS gradient increases below 2500 m, but decreases between 2500 and 4000 m, indicating that EOSS at high and low altitudes show a contrasting response to climate change. The analysis above further confirms that the effect of EOSS dates on the length of the SCD on the Atlas Mountains is greater than that of SOSS dates. This decrease in the SCD gradient was largely driven by an earlier snowmelt and partly by a later snow onset.
The direct reasons for this change in the snow metric gradient must be linked to the long-term change and variability of air temperature and precipitation. The proof of this connection can be provided by climate projections over Morocco that agree on robust warming and drying trends under greenhouse gas forcing. According to Marchane et al. [44], future projections in the High Atlas indicate an increase in temperature (+1.4 • C to +2.6 • C) and a decrease in total precipitation (from −22% to −31%), depending on the emissions scenario. Table 3 shows the difference between the average of the snow metrics lapse rate over two periods of 8 years each. This comparison also shows that the SCD and EOSS gradients decreased by 0.67 day and 0.35 day, respectively, between the two periods. In contrast, the SOSS gradient shows a positive change signal of about 0.5 day. The gradient variability of the snow metrics was analyzed by watershed for two reasons: firstly, because it is the scale of decision making for water managers, and secondly, in order to compare it with the results found at the Atlas scale. Figure 11  in order to compare it with the results found at the Atlas scale. Figure 11 presents the variability of the altitudinal gradient of the snow metrics for the seven watersheds covering the Moroccan Atlas, for the period of 2001-2016. The results highlight the strong interannual variability of the altitudinal gradient in relation to the strong spatio-temporal variability of snowfall for all the studied basins. From 2007 to 2016, we notice a decreasing trend in the altitudinal gradient of the SCD.
For the Oum-ER Bia watershed, which is the most important in terms of snow surface among the seven studied basins, the SCD gradient decreased from 7 days/100 m to 3 days/100 m. For the Tensift, Moulouya, Sebou, Souss, and Ziz basins, the SCD gradient varied from 6 days/100 m to less than 3 days/100 m. However, the Draa watershed has the maximum altitudinal gradient, at around 9 days/100 m, recorded in 2010-2011. In general, the variation of the SCD gradient at the watershed level is identical to that found at the Atlas scale from 2007 to 2016, with some differences in terms of the magnitude of the trend.
The variation of the EOSS gradient across watersheds is not consistent. For the period 2007-2016, there are only four basins (Oum-Er Bia, Draa, Souss, and slightly Tensift) that show a decreasing trend in the EOSS gradient (as seen in the Atlas). The results of a linear regression analysis for the basins Sebou, Ziz, and Moulouya show no clear trend during the period 2007-2016.
Regarding SOSS, the gradient variation results for all the basins examined are close to the Atlas findings, with an increasing trend. The greatest variation in the SOSS gradient is recorded in the southern watersheds, which are characterized by low snow surfaces and high variability in snow season starts.

Conclusions
The study highlights the topographic controls on the snow-covered area of the Atlas Mountains based on snow cover estimation during the 2001-2016 period. The following conclusions are drawn: For the Oum-ER Bia watershed, which is the most important in terms of snow surface among the seven studied basins, the SCD gradient decreased from 7 days/100 m to 3 days/100 m. For the Tensift, Moulouya, Sebou, Souss, and Ziz basins, the SCD gradient varied from 6 days/100 m to less than 3 days/100 m. However, the Draa watershed has the maximum altitudinal gradient, at around 9 days/100 m, recorded in 2010-2011. In general, the variation of the SCD gradient at the watershed level is identical to that found at the Atlas scale from 2007 to 2016, with some differences in terms of the magnitude of the trend.
The variation of the EOSS gradient across watersheds is not consistent. Regarding SOSS, the gradient variation results for all the basins examined are close to the Atlas findings, with an increasing trend. The greatest variation in the SOSS gradient is recorded in the southern watersheds, which are characterized by low snow surfaces and high variability in snow season starts.

Conclusions
The study highlights the topographic controls on the snow-covered area of the Atlas Mountains based on snow cover estimation during the 2001-2016 period. The following conclusions are drawn:

•
The snow cover is mainly affected by elevation, inducing a remarkable snow lapse rate. The quantity and timing (onset and melt-out dates and duration) of snow cover share the same pattern in their vertical distribution across the Atlas Mountains.

•
The altitudinal gradient of the SCD shows a decreasing trend. Since 2007, the SCD gradient in the Atlas Mountains has decreased from 6 days/100 m to 3 days/100 m. However, the change in the altitudinal gradient of SOSS is marked by a positive trend.

•
The altitudinal gradient of EOSS presents a potential contrasting response to climate change, with an increasing trend for elevations below 2500 m and a decreasing trend for altitudes between 2500 and 4000 m. The controls of orientation on snow cover metrics are not clear enough at the Atlas scale and even by basin, which can be related to the coarse spatial resolution of MODIS images that do not take into account the aspect control. The snowpack homogeneity can be considered as an important factor for spatial snow distribution associated with other topographic features. Both SCA and SCD increase with the slope. This can be explained by the fact that a significant part of the steep slopes is located at high altitudes, and steeper slopes receive less radiation than flatter surfaces.
The findings in this study represent an important advance towards the understanding of snow dynamics over the Atlas chain in Morocco and will contribute to improving our knowledge on the availability of water resources in relation to potential future climate change impacts. The analysis of snow lapse rate changes, as performed in this paper, constitute a novel approach that may be considered when producing a future scenario for water management at the river basin scale.