A method for robust estimation of snow seasonality metrics from Landsat and Sentinel-2 time series data over Atlas Mountains scale using Google Earth Engine

Abstract

Seasonal snow cover provides the majority of freshwater supplies for human society and natural ecosystems especially in semi-arid regions. For water resource managers, precise data regarding the spatiotemporal variability of snow cover and snow phenology is of paramount importance. Owing to the great spatial and temporal heterogeneity of the snowpack and the inaccessibility of high mountainous areas, gapless satellite remote sensing presents an unprecedented opportunity to monitor snow cover effectively and affordably on a fine scale, from different aspects, and with regular revisit time. This study derives the snow seasonality metrics (first day of snowfall, last day of snow melt; and snow cover duration) over a large semi-arid region in Morocco’s Atlas Mountains. We calculate these metrics by combining over 10,000 images from Landsat 8 and Sentinel 2 satellites for four hydrological years (2016–2021) to create a harmonized product with a time interval of about 3 days using the Google Earth Engine (GEE) platform. This dense and large time series facilitates a gap-filling method to minimize and overcome the effect of cloud cover, and its assessment shows a positive correlation between the masked pixels and the interpolated ones. These methods allowed us to realize a map of the snow cover area and extract a homogeneous Normalized Difference Snow Index (NDSI) profile over the four years whereby we were able to determine the optimal threshold to separate the presence of snow from its absence. The results showed that derived snow cover metrics provide considerable variation in both time and space, where an increase in snowpack measurement values at higher elevations can be noted. Overall, the snow duration ranges between November and April depending on the characteristics of each hydrological year. The retrieved Landsat 8 and Sentinel 2 snow dates had a high level of agreement with in-situ data observations with almost a day-and-a-half delay with an overall accuracy equal to 0.96. The analysis of snow cover dynamics via GEE has offered the ability to calculate the first day of snowfall, last day of snowmelt, and snow cover duration annually at a pixel level, providing the user with the ability to track the seasonal and interannual variability in the timing of snowmelt toward a better understanding of how the hydrological cycles of higher latitude and mountainous regions are responding to climate change.

Keywords:

• Snow metrics
• Landsat 8
• Sentinel 2
• Google Earth Engine
• cloud masking
• gap filling

1. Introduction

Snow cover is a highly sensitive indicator to climate variability due to its connection with the atmosphere and surface of the earth (Salminen et al. 2009), and fluctuates considerably from year to year, both in amount and duration due to environmental and climatic factors (Xuejin et al. 2019). The seasonal pattern of variation in the properties of snow cover on the regional or global scale is called snow phenology. Since snow is regarded as a freshwater source that fills reservoirs, recharges aquifers, and provides water supply for irrigation (Jackson et al. 2001), therefore, studying and monitoring the spatiotemporal variation in snow phenology parameters including the snow cover duration and snow melt-out dates has become paramount under current climate trends in which snow accumulation will diminish and aridity will be exacerbated (López‐Moreno et al. 2020).

In arid and semi-arid regions like central Morocco, this holds particularly true this, where the Atlas Mountains represent the most important water resource for the neighboring arid plains through liquid but also most importantly solid precipitation (Boudhar et al. 2007; 2009). However, its precipitation patterns are characterized by wide spatial and temporal variability and the dynamics of snow cover in mountainous areas change rapidly, with often snow falling and melting within a week (Boudhar et al. 2010; Marchane et al. 2015; 2021). Thus, it is important to quantify the existing snowpack and to better understand the snow water equivalent (SWE – the amount of available water in the snowpack) available during the melt process. In this context, snow cover can be monitored using in-situ measurements and weather station networks that constitute a long-term, reliable measurement of depth, density, snow water equivalent, and meteorological conditions. However, due to the limited coverage of ground observations and meteorological stations as they are almost always collected as discrete locations with uneven sparse distribution, the spatiotemporal characterization and monitoring of snow phenology in large-scale heterogeneous regions such as Morocco is challenging (Studer et al. 2007).

Bridging the gaps that exist between point-based measurements and models, these limitations can be addressed via satellite remote sensing technology which, via its global coverage and regularly repeatable observations, constitutes a powerful tool that since the 1960s, has provided the ability to quantitatively study the physical properties of snow distributed across different spatial scales, and also in otherwise inaccessible areas where measurements may be expensive, inaccessible, and/or dangerous (Nolin 2010).

In this context, the snow cover has been mapped and monitored globally by the National Oceanic and Atmospheric Administration-Advanced very high-resolution radiometer NOAA-AVHRR (since 1980), Moderate Resolution Imaging Spectroradiometer TERRA-MODIS, and SPOT-VEGETATION (since the end of 1990s). However, these sensors’ coarse spatial resolution (>250 m) results in a highly mixed composition and thus can obscure fine spatial patterning in snow variation (Bouamri et al. 2021; Bousbaa et al. 2022). Other satellites, like the Landsat 8/9 and Sentinel 2 provide finer spatial resolution measurements of snow at the expense of having a longer temporal return period when they have been used as a single optical sensor which may pose more challenges in the timely detection of snow phenology information (Baba et al. 2021).

To address this deficiency, the combined use of Sentinel-2 sensor images with Landsat images can result in additional time series implying a global average revisit observation interval of ∼ 2.9 days (Li and Roy 2017). However, one issue remains controversial: managing cloud cover is a major challenge in correctly estimating snow cover for optical satellite remote sensing applications (Christopher and Gupta 2010). Generally, for scientists trying to analyze objects on the surface of the earth, clouds are one of the primary obstacles. Their presence severely limits the ability and usefulness of optical satellite imagery, unless its purpose is specifically to analyze the cloud itself. Thus, cloud masking is a necessary pre-processing step for any information derivation particularly for multi-temporal satellite data analysis (Shu et al. 2022).

There is not many methods have been developed for estimating snow phenology parameters and events from satellite time series. A common approach is to calculate them based on the binary retreated snow cover data from Moderate-Resolution Imaging Spectroradiometer (Wu et al. 2022; Wang et al. 2023). However, new strategies could be developed using similar approaches reported by vegetation studies which are based on the use of multitemporal indices such as the Normalized Snow Difference Index (NDSI) to produce high-resolution snow phenological information (Mityók et al. 2018).

However, the use of high-resolution products over a large area leads to the creation of large repositories of satellite imagery; and the challenge here is its increasing volume, which is too large to be processed, and thus requires new ways for efficiently and consistently extracting relevant information (Ma et al. 2015; Htitiou et al. 2021). To overcome this issue, the use of cloud computing platforms such as Google Earth Engine (GEE) simplifies file management providing a full set of remote sensing data that can be readily accessed with a modern web browser and 1-2 lines of code. A multi-petabyte curated collection of widely used geographic datasets (Gorelick et al. 2017) may be found in the Earth Engine public data catalog, and are already processed on multiple levels. The bulk of the catalog consists of Earth Observation Remote Sensing imagery (including the entire Landsat archive and a complete archive of Sentinel), in addition to many other environmental, geophysical, and socioeconomic data sets. The average latency from scene capture is about 24 h, and the catalog is updated regularly at a frequency of nearly 6000 scenes per day from operational missions.

GEE provides two additional advantages (Sproles et al. 2018) all processing occurs on large cloud-based servers which greatly speeds up processing time and removes the requirement to download and archive large datasets. This is particularly relevant for remote sensing products that comprise a longer-term, ongoing analysis such as monitoring snow (Gorelick et al. 2017).

To address the challenges of snow cover monitoring using big data analytics, the general objective of this paper can be delineated as follows: (a) Establishing a systematic and regular snow cover monitoring system across the Moroccan Atlas Mountains; (b) Exploiting high spatio-temporal resolution data to enable accurate and detailed monitoring of snow cover dynamics; (c) Using the Google Earth Engine (GEE) cloud platform as a robust tool for efficient and scalable processing of large datasets involved in snow cover analysis; (d) The combined use of Landsat 8 (L8) and Sentinel-2 (S2) time series data for an integrated and comprehensive analysis of snow cover dynamics; and (e) The derivation of a spatially and temporally detailed set of seasonal snow cover metrics, including First Day of Snow (FDS), Last Day of Snow (LDS), and snow duration (SD).

2. Materials and methods

2.1. Study area

The scope of this research encompasses the Kingdom of Morocco, with a particular focus on snow-covered regions above 1,200 meters in the Atlas Mountains. The Moroccan Atlas Mountains (Atlas henceforth), are separated from the Rif by the Fez-Taza corridor, which runs diagonally down the backbone of the country from the Northeast to the Southwest, separating the plains and plateaus of central Atlantic Morocco from the vast pre-Saharan and Saharan areas. reaching elevations of 4000 m (Htitiou et al. 2021). Climatologically, the Atlas represents the confluence of two distinct air masses: the colder, wetter air masses that originate from the north and the hot and dry air masses that move up from the south. As a result, the climate of Morocco is characterized by a dry summer and wet winter, which starts in October and lasts until April, with the wettest months from December to February. Additionally, the whole region is characterized by high inter-annual precipitation variability with long-term mean precipitation reflective of many dry years and some relatively humid years (Schilling et al. 2012). Particular attention was given to a catchment within this region. This catchment, denoted in (Figure 1(c)), was strategically selected to enhance the resolution and focus of our analysis. It is situated in the Moroccan High Atlas, south of the city of Marrakech, and covers 225 square kilometers (Marchane et al. 2017). Rheraya catchment features a steep altitudinal gradient, with altitudes ranging from 1027 m to 4167 m, culminating at the summit of Mount Toubkal, the highest point in North Africa. The climate is characterized as semi-arid. The relative contribution of snowmelt to annual river flow varies between 30% and 48%, depending on the year. (Boudhar et al. 2009)

Figure 1. (a) Moroccan Kingdom (b) Moroccan Atlas Mountains (c) Rheraya catchment and Oukaimden_SM snow depth station.

2.2. Study environment and data acquisition

2.2.1. Google Earth Engine platform

The GEE platform was accessed and managed through a web-accessible application programming interface (API) with a web-based interactive development environment (IDE) for rapid prototyping and visualization of results. We implemented existing GEE client libraries and developed custom code in the IDE using JavaScript.

2.2.2. Satellite data

In mountainous regions, during the wet, winter months a higher percentage of satellite observations are cloudy, obscuring the land surface signal. To improve our ability to quantify the full range of snow change patterns we integrated L8 and Sentinel-2, imagery with similar spectral and spatial resolution. L8 images synthesized with Sentinel-2 fill data gaps obstructed by clouds, to create harmonized, higher resolution images that reduce error associated with long prediction periods, high spatial heterogeneity, and fine resolution changes. L8 was launched in 2013, and its Operational Land Imager (OLI) provides high-quality multispectral images at a resolution of 30 meters (15 for panchromatic) and a revisiting time of 16 days (Du et al. 2014). In this study we used L8 Collection 2 Tier 1 and real-time data calibrated Top-of-Atmosphere (TOA) Reflection Tier 1 includes Level-1 Precision Terrain (L1TP). All Tier 1 Landsat data can be assumed to be consistent and inter-calibrated (Chander et al. 2009). The Sentinel-2 mission provides (Baillarin et al. 2012) for a combination of two satellites Sentinel-2A and Sentinel-2B equipped with identical Multispectral Instruments (MSI) capable of acquiring data in 13 bands at different spatial resolutions (between 10 m and 60 m); Orthorectified Top-Of-Atmosphere (TOA) reflectance is provided by S2 Level-1C product, with sub-pixel multispectral registration. The item includes masks for the clouds and land/water.

These image collections were spatially filtered for our Region of Interest (ROI), the Moroccan Atlas, with a final data archive of 10260 scenes that extend from (2016 to 2021)

2.2.3. In situ data

In our study area, the Oukaimeden_SM automated meteorological station (AWS) provides a wide range of hydrometeorological variables as detailed by Boudhar et al. (2016). Located at the Oukaimeden ski resort, about 70 km south of Marrakech at an altitude of 3200 m a.s.l (Figure 1) (Hanich et al. 2022) on the northern side of the High Atlas. The station represents the progressive installation of its devices over six years from 2003 to 2010.

These in situ data provide the only field measurements to test the efficacy of the methods adopted for snow cover monitoring and also to examine the values of the combined derived snow metrics at a point scale (Table 1).

Table 1. Snow depth instrumentation at the Oukaimden_SM weather station and their specifications (2003–2010). (Boudhar et al. 2016).

Quantity	Instrument type	Measurement range	Installation Period
Snow depth	Sonic CR50	0.5–10 m	2003–2010

2.3. Methodology

We calculate the snow seasonality metrics across hydrological years from September 2016 to August 2021 to characterize the spatial and temporal distribution of snow cover using a harmonized L8 and Sentinel-2 product to quantify the snow seasonality metrics (FSD, LSD, SD). These satellite-derived metrics are compared to the AWS in-situ data at Oukaimeden_SM. We accomplished this through a four-step workflow described visually in Figure 2.

Figure 2. General workflow.

2.4. Data pre-processing

2.4.1. L8 and S2 band adjustment

Merging L8 & S2 collection offers the capacity to better detect ephemeral changes, nonetheless, different sensor characteristics (different spectral response functions, view, and solar angles) provide pre-processing challenges and must be corrected before analysis. Castain (Chastain et al. 2019) proposed cross-sensor transformation coefficients (Table 2), which we used to adjust the two common bands of those sensors (green, Short-Wave Infrared1 swir1) using Equation 1(1) $$\mathrm{\text{Sentinel}}\mathrm{ }2=\text{Landsat}8\text{* Slope}+\mathrm{\text{Intercept}}$$(1) for each Band. (1) $$\mathrm{\text{Sentinel}}\mathrm{ }2=\text{Landsat}8\text{* Slope}+\mathrm{\text{Intercept}}$$(1)

Table 2. Cross-sensor transformation coefficients for Landsat 8 (L8) according to Chastain et al. (2019).

Bands	L8 Intercept	L8 Slope
Green	0.0026	1.0043
SWIR1	0.0065	1.0049

2.4.2 Cloud masking

The L8 (30 m) cloud detection algorithm used for cloud masking incorporates the C implementation of the Function of Mask (CFMask) algorithm, producing a quality assessment band (QAB). CFMask, known for its superior accuracy, is recognized as one of the most advanced cloud detection algorithms (Mateo-García et al., 2018). This algorithm is not only optimized for L8 imagery but has also been extended to be compatible with Sentinel-2 data. While L8 uses the BQA band and thermal infrared bands, Sentinel-2 uses a set of three quality assessment bands, QA60 being a bit mask band providing valuable cloud mask information. The adaptability of the algorithm to the distinct spectral characteristics of the two satellites underlines its effectiveness in detecting clouds across a variety of satellite systems.

2.4.3. Normalized difference snow Index (NDSI)

To exploit the high reflectance of snow in the visible and the shortwave infrared wavelength bands, and to accentuate the presence of snow in our study area, we employed NDSI after masking the cloud process. Here, we derived our own NDSI from L8 and S2 products using the ‘normalizedDifference’ function in GEE using the normalized difference between two bands (Equation 2(2) $$\mathit{\text{NDSI}}=\frac{\mathit{\text{Green}}-\mathit{\text{SWIR}}}{\mathit{\text{Green}}+\mathit{\text{SWIR}}}\mathrm{ }$$(2) ). (2) $$\mathit{\text{NDSI}}=\frac{\mathit{\text{Green}}-\mathit{\text{SWIR}}}{\mathit{\text{Green}}+\mathit{\text{SWIR}}}\mathrm{ }$$(2)

NDSI values range from −1 to 1, with values greater than 0.4 generally considered to have snow present (Xiao et al. 2002; Kulkarni et al. 2006; Sibandze et al. 2014) and was used to quantify the spatiotemporal extent of snow cover (Riggs et al. 1994; Wang et al. 2018).

2.4.4. L8 and S2 combination and compositing

The combination of L8 and S2 collections is created based on the merging of the two resultant NDSI bands. The goal is to develop a data stream with a high revisit time (every 3 days), which increases the number of scenes in the collection for a better interpolation of the information lost with the presence of clouds. However, this combination leads to a dense time series (Figure 3(c)) that causes problems during processing. Plotting observations in GEE for the L8 and S2 NDSI imagery resulted in many computational timeouts, especially when charting time series data. As a result, we opted for the compositing method to address this problem. Conceptually and computationally compositing represents the integration of spatially overlapping images into a single image using an aggregation function. We extracted the maximum value of NDSI per pixel per 3 days, essentially reducing the number of images by a factor of 3.

Figure 3. Coverage footprints for Morocco of: (a) Landsat 8 (L8) Satellite, (b) Sentinel 2 (S2) Satellite, (c) Combination of L8 &S2 Satellites/Pixel frequency distribution Histogram of Satellites scenes.

2.4.5. Gap-filling method

Unlike built land or vegetation (such as cropland), snow is less affected by human activities and is primarily impacted by environmental factors that drive its formation and melt. The distribution of snow is driven by meteorological (i.e. precipitation, temperature, wind, and solar radiation) and topographic (altitude, slope) factors. Thus, when two pixels are remarkably similar, the two adjacent pixels’ NDSI snow state may also be very similar (Li et al. 2020). This similarity includes not only local neighborhood similar information (Tobler’s first law of geography (Gilboa and Osher 2009)) but also non-local similarity.

When a pixel is covered by clouds, its NDSI information is lost. Auspiciously, this information can be estimated from similar pixels in the remaining known areas at the same date and time.

To overcome the problem of cloud contamination and to fill spatial and temporal gaps created by the cloud masking process we adopted the Whittaker Smoother method (Weinert 2007), a spatio-temporal interpolation across images in the image collection to reduce noise due to clouds and atmospheric conditions.

The advantages of the Whittaker Smoother lie in that it has a single smoothing parameter (λ) and high computational efficiency. For example, a series of 100,000 observations can be smoothed in a few seconds with an ordinary computer (Weinert 2007) and has broad prospects for large-scale high-resolution remote sensing data applications on the GEE platform. The Whittaker Smoother smooths time series by minimizing a cost function (Zhou et al. 2016) that combines the precision and roughness of the smoothed time series. The cost function (C) can be expressed as: (3) $$C=A+\mathrm{ʎ}R\mathrm{ }$$(3)

Where A measures the accuracy of the fitting, R represents the roughness of the fitting result, and λ is a coefficient controlling the smoothness of the fitting result.

2.4.6. Validation of the gap-filling method based on a ‘cloud mask assumption’

When a pixel is covered by clouds, any NDSI information is null for that time image. However, to evaluate the precision of the estimated NDSI pixels we need the real NDSI pixels. To overcome this limit, an original validation method based on a ‘mask assumption’ was applied. The premise behind this approach is to create a synthetic cloud mask of known and no cloudy NDSI values and then compare the restored NDSI with the known NDSI. This quantifies the skill of the Whittaker Smoother to gap-fill pixels of areas that were synthetically masked compared to known NDSI values (Figure 4). The performance of the adopted gap-filling method in this study is evaluated using a standard correlation coefficient. The value of the correlation coefficient is between −1 and 1. ‘0’ means that there is no relationship between variables at all, and −1 or 1 means completely negative or positive correlation.

Figure 4. Gap filling’s evaluation strategy.

2.5. Snow seasonality metrics

The extraction and analysis of snow seasonality metrics involve several crucial steps. To initiate the process, a date range for analysis is defined. This analysis relies on a combination of L8 and S2 NDSI data. These bands are processed annually. They are then enriched with date-related information, including Calendar Day of Year, Relative Day of Year Milliseconds since the Unix epoch, and the corresponding year, which is vital for tracking temporal changes in snow cover. To ensure a focused analysis, an analysis mask is defined to filter out specific areas, such as water bodies, leaving the focus only on land areas. For each year, the dataset is filtered for the specific year of interest, date bands are added to images to facilitate temporal analysis, and images are sorted by date to identify temporal changes in snow cover. The sorting is reversed for snowmelt detection to pinpoint the disappearance of snow on the ground. a threshold is also applied: The first value of snow greater than this threshold is considered as FDS. Similarly, in the reverse sorting, the first value encountered after the determined threshold is considered the LDS. Once these key dates are determined, the SD for each year can be calculated by measuring the period between these two events.

2.5.1. Derived L8 and S2 snow onset and snow melt dates

To calculate the first day of snow (FDS), no snow (LDS), and snow duration (SD) (Figure 5) annually for all water years (September–August) at the pixel level, allowing us to track seasonal and inter-annual changes in snowmelt timing to better understand how the hydrological cycle in high altitudes and mountains responds to inter-annual climate variability.

Figure 5. Example of the extraction of the snow seasonality metrics.

We defined a date range where we specify the day-of-year (DOY) to start the search for the snow cover at the beginning of the water year (September 1) until the end (August 31). Each image has a metadata timestamp property, but since we created a yearly mosaic of images that combine pixels from many images, the date must be encoded as a value in each pixel for the final multidimensional mosaiced image. From here we defined a function in GEE that adds multiple date bands to the images containing a band for DOY (calculated from 1st January), Relative Day of Year calculated from 1 September (RDOY), Band for Time in Milliseconds, and band for the Year. Finally, we constrained our analysis geographically through a water body mask (i.e. lakes, reservoirs, and rivers), an elevational mask >1000 m), and individual pixels that do not meet a threshold for snow for at least 5 days a year.

2.5.2. NDSI threshold optimization

As mentioned above, the value of the NDSI index generally varies between −1 and 1, and a differentiation between snow and no snow is based on a standard threshold value which is commonly assumed to be 0.4 (Hall and Riggs, 2007), however, the spatiotemporal representativeness of this standard threshold is questionable at the local scale. Therefore, the optimization of the NDSI threshold is very important for obtaining the snow seasonality metrics and plays a major role in achieving a highly improved detection performance. In practice, we applied a binary (0,1) reclassification of each pixel within an image into a graduated series of thresholds at 0.05 increments based on its NDSI value (Brownlee 2020). For example, each NDSI image for a given date would be reclassified into a data cube of 20 sub-images (0.05–1.0). Simply put, for the 0.5 sub-image any pixel with a NDSI value ≥ 0.5 would be classified as a 1, referring to the presence of snow, or as a 0 if < 0.5, referring to its absence (Figure 6).

Figure 6. Converting probabilities to class labels.

The binary reclassification was based on the following logic for each NDSI threshold: $$\text{Prediction}<\text{Threshold}=\text{Class}0\to \text{lack of snow cover}$$ $$\text{Prediction}\ge \text{Threshold}=\text{Class 1}\to \text{indicates the presence of snow cover}$$

To validate our binary classification outputs, we created a classified image considered as a reference image for comparison with each predicted sub-image to identify which threshold can be an adequate solution for our case.

For its acquirement, we used all the bands of a cloud-free S2 image captured on 23 March 2018 with a dimension of 537,557 × 299,109 and its accompanying NDSI product to demonstrate the variation of the snow in the ROI.

We applied a Support Vector Machine (SVM) classification over the stack of these bands. The SVM classification is a supervised statistical learning technique used to solve classification problems (Mountrakis et al. 2011) and solves a quadratic optimization problem to determine the optimal separation limits (hyperplanes) between two classes in a multi-dimensional feature space. When searching for the optimal threshold, several functions based on the classification performance can be used (Zhang and Hu 2017), among which we can find Precision, Sensitivity, Accuracy, and F1 score (Zhang and Hu 2017), those parameters were included in the confusion matrix.

Precision is defined as the fraction of correct predictions for a certain class (Equation 4(4) $$\mathit{\text{Precision}}=\frac{\mathit{\text{True positives}}}{\mathit{\text{True positives}}+\mathit{\text{False positives}}}$$(4) ). (4) $$\mathit{\text{Precision}}=\frac{\mathit{\text{True positives}}}{\mathit{\text{True positives}}+\mathit{\text{False positives}}}$$(4)

Sensitivity is defined as the fraction of instances of a class that were correctly predicted (Equation 5(5) $$\mathit{\text{Sensitivity}}=\frac{\mathit{\text{True positives}}}{\mathit{\text{True positives}}+\mathit{\text{False negatives}}}\mathrm{ }$$(5) ). (5) $$\mathit{\text{Sensitivity}}=\frac{\mathit{\text{True positives}}}{\mathit{\text{True positives}}+\mathit{\text{False negatives}}}\mathrm{ }$$(5)

Accuracy is defined as the fraction of instances that are correctly classified.

F1 Score is defined as the harmonic mean between sensitivity and precision (Equation 6(6) $$\mathrm{F}1=2\mathrm{*}\frac{\mathit{\text{Precision}}\mathrm{*}\mathit{\text{Sensitivity}}}{\mathit{\text{Precision}}+\mathit{\text{Sensitivity}}}$$(6) ). (6) $$\mathrm{F}1=2\mathrm{*}\frac{\mathit{\text{Precision}}\mathrm{*}\mathit{\text{Sensitivity}}}{\mathit{\text{Precision}}+\mathit{\text{Sensitivity}}}$$(6)

Kappa: the difference between the overall accuracy and the expected accuracy divided by 1 minus the expected accuracy.

2.6. Accuracy assessments of L8 and S2 metrics against in situ snow depth data

To evaluate the satellite-derived data, we integrated the L8 time series from 2013 - 2016 by merging it with the combined data of both L8 and S2 from 2016 till 2021 in pursuance of getting more variables to adopt in the accuracy assessment. We obtained seasonally-paired dates of snow onset and snow melt between 1 September 2013 and August 2021 using 30-minute measurements of snow depth from the Oukaimeden_sm AWS (Marchane et al. 2015) down-sampled to daily data. Figure 7 visually describes the daily variation of snow depth at Oukaimeden_sm AWS for nine years from 2013 to 2021 as compared to the NDSI values.

Figure 7. Variation of snow Index against snow depth between 2013 and 2021 at the Oukaimeden station.

Using these datasets, in situ snow onset and melt dates from the AWS were converted to a day-of-snow-year for direct comparison to the L8 & S2 derived snow metrics (1 September to 31 August). We evaluated the skill of these time series of L8 & S2-derived metrics of paired snow onset and melt dates to data (first day of snow, last day of snow, and snow duration) from the in-situ AWS using R² and root mean square error (RMSE) over the entire time series (2013–2021).

3. Results

3.1. Combined use of Landsat 8 and Sentinel-2

Time-series analysis can be done at a frequency denser than when using a single data source by combining L8 and S2 imagery.

As shown in Figure 8, we mapped the NDSI dynamics of a snowy pixel over the Oukaimeden region by using only Sentinel-2 imagery (orange profile), L8 imagery (Red profile), and the combined imagery (blue profile) respectively. the twin Sentinel-2 satellites have a 5-day revisit cycle, enabling us to collect 72 images in 2018. then again L8 with 16 revisit times provided 20 more usable images. by contrast, the combined use of the two data sources led to a relatively dense time series and more efficient temporal variation characterization as shown in the figure (purple line) is composed of more than 120 values of NDSI with a time step of 3 days.

Figure 8. Illustration demonstrating the combined use of Landsat 8 (L8) and Sentinel 2 (S2) NDSI per pixel.

The fusion of L8 and S2 data improves the temporal resolution of the NDSI data and enhances the detail of the variability of snowpack.

3.2. Results of cloud removal and gap-filling method and generation of gap-filled NDSI time series

The spatial gaps in the NDSI images from the cloud masking process were filled with the Whittaker Smoothing method. Results show a new grid with a smoothed NDSI value displaying almost no difference between the smoothed values and the grid with the real NDSI values Figure 9.

Figure 9. Zoomed view of the gap filling result: (a) cloudy masked pixels; (b) smoothed NDSI value.

To validate the NDSI prediction and gap-filling results by the proposed method, we chose an image that does not have any cloud cover captured on the 9 March 2018 to obtain uncontaminated NDSI values. Masking a portion of the image (700 km²) to estimate temporally the lost information values of these pixels, we applied a set of multi-temporal reference images before and after the date of the image to interpolate. Using five images over a time interval of 3 days, to assess the efficacy of the Whittaker Smoothing (Figure 10). The correlation between the two variables indicates a strong relationship between the measured and estimated NDSI with an R² of 0.92 and an RMSE of 0.05 (Figure 10(d)).

Figure 10. Gap filling evaluation process (a) real NDSI value; (b) masked NDSI (c) smoothed NDSI. (d) Density profile between real & smoothed scene.

To expand the validation of the gap-filling we extended our validation to the entire time series. Figure 11 shows cloud effects over a 3-day time-averaged NDSI time series over the 2017-2018 hydrological year for a single location (1 pixel). The cloud masking feature caused some gaps that did not indicate any geographic information (green line), and we adopted as a solution the gap-filling techniques to overcome this problem in which a linear interpolation is applied on the time series to fill those cloud gaps; the blue line refers to the estimated values (replacing the cloudy masked ones) of NDSI from the before and after cloud-free pixels.

Figure 11. Cloud mask and gap filling impact on NDSI time series season 2017/2018 (the profile represents the NDSI values of one pixel).

3.3. Threshold optimization for snow seasonality metrics

Identifying a single ‘optimal’ threshold that produces the highest possible accuracy of snow extent mapping was done using sensitivity analysis to optimize threshold values for snow characteristics. Seventeen thresholds (T1–T 17) were tested by increasing the NDSI data by 0.05 from an NDSI of 0.1 to 0.9, and results were compared to the classifications produced by the SVM. Optimized skill scores of the largest Accuracy, Kappa, and F1 scores indicate an optimal NDSI of 0.45, with an accuracy of 0.86 (Figure 12).

Figure 12. Snow-based thresholding optimization.

3.4. Snow seasonality metrics

After summarizing the results with a series of visualizations that display spatial maps referencing the date of snowmelt/fall and SD respectively.

The figures below show the variation of the derived seasonal metrics (FSD, LSD, SD) that depict the beginning (Figure 13), the end (Figure 14), and the duration (Figure 15) of the time between the first snowfall and the last snowmelt. Calculated for the entire Atlas over the hydrological years from 2016 to 2021 and a zoomed view over the Rheraya basin as an example; where we can notice that the annual variance of the first day of snowfall and last day of snow melt was intense in high elevation areas and comparatively weak in low elevation areas of the Moroccan Atlas Mountain range over all the six hydrological years (Figures 13, 14), and SD is highly variable (Figure 16). The pixels with a high altitude represent a greater duration of snow cover compared to other pixels located at full height (low and medium altitude). We found variations in both time and space, with the high-altitude region often exhibiting the longest snow duration (250 ± 5 days), in contrast to the lowest locations, which have a significantly shorter snow duration (45 ± 3 days). As the value of snow duration varied between 42 and 250 days from 2016 to 2021, a significant interannual fluctuation was also noted. In-depth analysis of snow metrics over the hydrological years from September 2016 to August 2021 revealed nuanced temporal patterns for each specific metric. Looking at the first day of snow (Figure 16), distinct concentrations were observed for different years. In the 2016/2017 hydrological year, the majority of the first snow days occurred in November and December, in contrast to 2017/2018, when concentrations spanned October, November, and December. In 2018/2019, the focus was on October, and for 2019/2020, the first snowfalls were observed in September, October, November and December. The following year, 2020/2021, featured a concentration in November. Complementing this, the last snow day metric (Figure 17) showed a general trend for most years, where snow persisted into February, March, and May, except for the anomaly observed in 2019/2020, where snow lingered mainly in April. Furthermore, investigating snow duration (Figure 18) revealed a variability ranging from 100 to 200 days for most years.

Figure 13. Variation of snowfall dates over the six hydrological years at the Moroccan Atlas Mountain Range.

Figure 14. Variation of snow melt dates over the six hydrological years at the Moroccan Atlas Mountain range.

Figure 15. Snow duration in Moroccan Atlas Mountain between 2016 and 2021.

Figure 16. Bar graphs showing the frequency distribution of Snowfall date at the Rheraya catchment scale.

Figure 17. Bar graphs showing the frequency distribution of Snow Melt date at the Rheraya catchment scale.

Figure 18. Bar graphs showing the frequency distribution of snow duration days at the Rheraya catchment scale.

3.5. Validation of L8 and S2 metrics against in situ snow depth data

The difference between the in-situ and L8 & S2 derived first day snowfall data is positive, indicating that the L8 & S2 derived snow season starts earlier than the weather station data (positive bias). The overall RMSE for the first snow day was less than 6 days, with a correlation coefficient (R²) of 0.98.

Data from L8 & S2 are biased late (negative bias) relative to in situ observations for the last day of the snow season and snow duration. The last day of snow melting season data exported by L8 & S2 differs by about 3 days from the in-situ snow melting data; and under 5 days for snow duration. The RMSE of the snowmelt dates is below 15 days with a correlation coefficient equivalent to 0.7; as for the snow duration, the total RMSE is less than 5 days with an R²= 0.93.

Figure 19. L8 and S2 derived snow metrics for the first day of snow (FDS), last day of snow (LDS), and snow duration (SD) versus in situ derived dates Doy 2013–2021. axes are DOY of hydrological year.

4. Discussion

As like in most arid and semi-arid zones, snow cover provides an important proportion of the freshwater resources that Morocco relies on and has a direct influence on the economy and livelihoods of several provinces. Therefore ensuring accurate and timely monitoring of the seasonal snow cover and its phenology over the years is very crucial to provide valuable insights for water management and water security. In this study, we proposed a novel approach that automatically extracts key snow phenological parameters in the Moroccan Atlas from 2016 to 2021, and produces high resolutions maps of the snow seasonality metrics including the First Day of Snow (FDS), Last Day of Snow (LDS) and Snow Duration (SD). The proposed approach was developed as a proof-of-concept.

Many studies have proved that optical-sensor such as MODIS data are more advantageous to characterize snow phenology on different scales (Wang et al. 2023). However, these coarse-resolution sensors (>250 m) are less reliable when applied to heterogeneous and fragmented snow surfaces at the local scale because a large number of mixed pixels exist in these coarse images. By offering increased observation density and higher spatial resolution (Figure 3), the combined use of L8 and S2 satellite data in this study has enabled a deeper comprehension of changing snow cover patterns (Gascoin et al. 2019) and also apparently increasing the accuracy of snow monitoring than is possible with either satellite used alone (see Figure 8).

Although Landsat and Sentinel-2 have similar bands and bandwidth, there were some surface reflectance differences (Claverie et al. 2018). In this study, a spectral band adjustment step has been applied to reduce these biases and facilitate the synergetic use of both datasets as discussed in Section 2.4.1.

To produce a good product of snow data from L8/S2 combined data, we had to fill data gaps due to clouds and cloud shadows in the optical data. However, effective identification and differentiation between snow and clouds in snow-covered areas is a significant problem when processing remotely sensed cloud-cover data (Stillinger et al. 2019). As mentioned, in section 2.4.5, a gap-filling method was proposed to fill and restore cloud-contaminated images using the weighted Whittaker filter, as more robust for outliers even when the input NDSI time series is more contaminated which is in line with the study of Kong et al. (2019) who proposed this filter in GEE for vegetation analysis. The application of this approach on artificial cloud mask as shown in section 2.4.6 has shown a high agreement with the validation set, which proves the robustness and effectiveness of the proposed approach.

The very large amount of data, high computing power, and efficient acquisition techniques available in the GEE cloud platform are key to this research to extract the snow seasonality metrics over a large area using high-resolution satellite images. We could need several days to finish such a task if the conventional process scheme was used, which is unacceptable in terms of labor costs and time delays. However, this operation only took a few seconds with the GEE platform (Gorelick et al. 2017; Hird et al. 2017).

The findings showed that derived snow cover metrics revealed significant variation in both time and space, where an increase in snowpack measurement values at higher elevations can be noted (Figures 13, 15, and 17).

Our results show potential for the use of L8/S2 to retrieve snow metrics, as illustrated by the correlation between in-situ data observations and NDSI-derived metrics (Figure 19) with almost a day and a half delay with an overall accuracy equal to 0.96. On the whole the snow duration ranges between November and April depending on the climate of each hydrological year. Overall, this result is encouraging, although there is no study to compare these results to, so it is important to interpret this with caution.

In general, a pixel is regarded as snow if the NDSI is ≥ 0.4, however, numerous local studies have recently questioned the general applicability of a standard NDSI threshold in snow and glacier monitoring (Yin et al. 2013). Therefore, threshold optimization should be done.

Different thresholds were tested to derive snow cover metrics from the snow depth data by employing a range of thresholds. A value > 0.2 was used to extract the first metric which is the SDF and a threshold value of 0.3 to extract the LDS from the snow depth series. We relied on threshold optimization processing for each parameter to find the optimal threshold to adopt. The results indicate that the snow depth interval likely to have an optimal threshold is between 0.1 and 0.25 for the FDS metric (Appendix Figure A1), while the optimal interval is between 0.25 and 0.35 for the LDS metric represented by (Appendix Figure A2). The difference between these intervals is due to several factors, including landscape (topography) and hydroclimatic factors (temperature and precipitation), in addition to the influence of water storage and release cycles. Large amounts of water are stored in the snowpack and then released during spring and summer snowmelt, which affects soil moisture and therefore snow depth measurements; otherwise, snowfall in these regions can also be in the form of shallow snow (less than 20 cm deep). This is difficult to measure accurately due to their height spatio-temporal variability.

It is important to note that results are site-dependent, requiring in situ data to develop threshold values, whether it be for the NDSI or the snow depth.

Even if AWS can provide snow cover observations, they can still require some uncertainty analysis as data are often incomplete and represent only a small area around the measurement point. This could affect assessments of long-term trends and create problems when using snowpack data to estimate runoff and predict future changes in freshwater availability in river systems especially when we survey large areas with only one site.

The reason that snow seasonality was not changing at high altitudes is related to other factors in addition to the effect of altitude that influences the spatial variation of snow cover duration in Morocco. changes of snow in response to climate change and the associated effects on hydrological regimes are of great concern. The snow cover of the Moroccan mountains depends on the climate prevailing in the country. For a given winter, the snow cover in terms of quantity and duration depends on the succession of meteorological episodes (disturbances, anticyclones, hot or cold periods), it cannot be simply linked to average parameters such as precipitation or average temperature. Therefore, the correlation between climate change and these parameter variations needs to be examined more closely, to better understand the complex relationships between snowfall and changes in the environmental context.

5. Conclusions

We applied the Google Earth Engine platform to generate a cloud-free L8 and S2 NDSI product and used these derived data to create a spatially and temporally detailed set of seasonal snow cover metrics for the Moroccan Atlas Mountains from across snow seasons (September 2016 to August 2021). These metrics include snow cover duration, first day of snowfall, and last day of snowmelt dates. An accuracy assessment using in situ snow onset and melt dates from Oukaimeden_SM weather station where snow depth was monitored showed that dates measured at in situ weather station location generally precede Satellites derived dates.

The use of remote sensing and the GEE platform facilitated processing and analysis as our study area extended over 2000 km due to the pre-processed data library. This facilitated streamlined analysis of the resultant data. These L8 and S2 derived metrics provide a useful tool for documenting inter-annual variability and long-term trends in snow cover characteristics. We can also move towards the fusion to enhance the performance; the fusion result of the two satellites could be used to extend the obtained results at 10 m spatial resolution.

Moving forward with this project to treat and calculate other snow metrics based on the first and last day of snow cover values in this study obtained among which we can find the continuous snow season as the longest interval of the season for which a pixel was snow covered and more.

Acknowledgments

The authors would like to thank the Google Earth Engine team for providing cloud-computing resources; and gratefully acknowledge the reviewers for providing insightful comments and recommendations, which significantly enhanced the work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Data supporting the results will be available on request from the corresponding author.

Additional information

Funding

This study was funded by: GEANTech project with the assistance of OCP Foundation, ESRI, CNRST and UM6P, as well as by the MorSnow research program within Mohammed VI Polytechnic University (UM6P), Morocco, (Accord spécifique n° 39 entre OCP Innovation et UM6P).

Appendix A. 

Threshold optimization of NDSI and in situ snow depth data

Figure A1. RMSE and R square confusion matrix result between NDSI and snow depth different threshold, detection of threshold based on FDS (day of year) of both NDSI and snow depth.

Figure A2. RMSE and R square confusion matrix result between NDSI and snow depth different threshold, detection of threshold based on LDS (day of year) of both NDSI and snow depth.