Changes in Snow Cover Dynamics over the Indus Basin: Evidences from 2008 to 2018 MODIS NDSI Trends Analysis

The frozen water reserves on the Earth are not only very dynamic in their nature, but also have significant effects on hydrological response of complex and dynamic river basins. The Indus basin is one of the most complex river basins in the world and receives most of its share from the Asian Water Tower (Himalayas). In such a huge river basin with high-altitude mountains, the regular quantification of snow cover is a great challenge to researchers for the management of downstream ecosystems. In this study, Moderate Resolution Imaging Spectroradiometer (MODIS) daily (MOD09GA) and 8-day (MOD09A1) products were used for the spatiotemporal quantification of snow cover over the Indus basin and the western rivers’ catchments from 2008 to 2018. The high-resolution Landsat Enhanced Thematic Mapper Plus (ETM+) was used as a standard product with a minimum Normalized Difference Snow Index (NDSI) threshold (0.4) to delineate the snow cover for 120 scenes over the Indus basin on different days. All types of errors of commission/omission were masked out using water, sand, cloud, and forest masks at different spatiotemporal resolutions. The snow cover comparison of MODIS products with Landsat ETM+, in situ snow data and Google Earth imagery indicated that the minimum NDSI threshold of 0.34 fits well compared to the globally accepted threshold of 0.4 due to the coarser resolution of MODIS products. The intercomparison of the time series snow cover area of MODIS products indicated R2 values of 0.96, 0.95, 0.97, 0.96 and 0.98, for the Chenab, Jhelum, Indus and eastern rivers’ catchments and Indus basin, respectively. A linear least squares regression analysis of the snow cover area of the Indus basin indicated a declining trend of about 3358 and 2459 km2 per year for MOD09A1 and MOD09GA products, respectively. The results also revealed a decrease in snow cover area over all the parts of the Indus basin and its sub-catchments. Our results suggest that MODIS time series NDSI analysis is a useful technique to estimate snow cover over the mountainous areas of complex river basins.


Introduction
Snow is a crucial component of the available freshwater resources in the world. It does not only hold significant hydrological importance in mid-to high-latitude mountainous areas [1], but also dominates the climate over mountainous ranges across the globe [2]. Mountains are the source of freshwater and about one sixth of the world's population receives freshwater from mountains that hold snow as temporal water storage and may have dynamic effects on hydrology [3][4][5]. Snow plays a key role in the hydrological response of complex river basins [6] and snow cover and its contribution to runoff has been focused on by several researchers worldwide [7,8]. The Indus basin is also a complex river basin fed by the Asian Water Tower, which holds several mid-to high-latitude mountainous glaciers and has significant effects on downstream hydrometeorological conditions [9]. The agriculture-based economy of Pakistan is dependent on the waters of the Indus basin and 90% of its flow originates from mountainous ranges of the Karakoram, Hindu Kush and western Himalayas [10].
The Indus basin receives runoff from rainfall as well as snowmelt [11] from low-to high-latitude mountainous ranges which are highly susceptible to high rates of snow melting due to recent global warming [12,13]. The snowmelt in the western Himalayas (the Indus and Sutlej rivers) accounts for almost 50% of the total runoff budget, while this percentage is less than 20% in the eastern and central Himalayan catchments [14]. The report of the International Panel on Climate Change (IPCC 4AR) revealed the adverse impacts of global warming in the form of increase in ocean and air temperature, the global melting of ice and snow and a subsequent rise in the sea level [15]. Snow cover is sensitive to climate change [16,17] and influenced by different climatic parameters like temperature and precipitation [16,18]. The monitoring of frozen water reserves is very essential for the proper regulation of water distribution [19], the appraisal of climate change impacts [16,20], the management of freshwater resources and predicting subsequent runoff [21,22]. For all these reasons, the study of snow cover dynamics is of great importance and a basin-scale study will not only help researchers to estimate the snowpack reserves, but also summarize the runoff seasonality.
The remoteness of areas and extreme hydrometeorological conditions pose serious challenges in ground-based measurements of snow cover. Alternately, remote sensing could be an appropriate way to acquire spatiotemporal snow cover information [23]. The technological advancements and multiplicity of remote sensors have increased the reliability of remotely sensed data. However, a higher frequency of snow cover estimations is also needed in hydrological and regional water balance studies [24,25]. Different daily Moderate Resolution Imaging Spectroradiometer (MODIS) satellite levels are being utilized by researchers to delineate fractional and binary snow cover. The Normalized Difference Snow Index (NDSI) can be helpful to estimate snow cover [26] over large mountainous basins on daily basis. The MODIS and Landsat Enhanced Thematic Mapper Plus (ETM+)-based NDSI threshold (≥0.4) has been applied by several researchers to estimate the snow cover area worldwide [27][28][29][30][31], while global threshold of NDSI greater than 0.4 can eliminate a significant amount of snow from estimations [32], especially in coarse-resolution datasets.
Snow cover estimation from remote sensing is subjected to different types of errors due to water, dark forests, snow in forests, barren land and clouds. Water pixels may be misclassified as snow in NDSI-based snow cover mapping [33], thereby increasing the uncertainties in the estimations. Water is an absorber for Near Infrared (NIR) radiation and pixels with a reflection greater than 11% in NIR are not mapped as snow [34,35] even if they have an NDSI greater than 0.4 [34,36]. Normalized Difference Water Index (NDWI) is also a good indicator used globally to delineate water bodies [37][38][39]. Dark forest areas may also introduce errors by indicating higher NDSI values [40,41] and the error of commission, in this case, can be minimized by fixing a threshold for visible reflectance [42]. Dark forest pixels with less than 10% reflectance in the green band cannot be mapped as snow [35]. Another limitation is the differentiation of snow pixels from snow-forest areas where both the NDSI and Normalized Difference Vegetation Index (NDVI) are higher. This problem can be resolved by incorporating NDSI and NDVI at the same time to detect snow in forest areas [34,43,44]. If the NDSI is less than the threshold (0. 4) and NDVI is about 0.1, then the pixel is classified as snow [36,45]. Errors due to clouds can reduce the overall accuracy of the final product [46,47] and different cloud algorithms for high, low and mid-latitude regions have been developed for cloud masking [26,[48][49][50]. The cloud pixels misclassified as snow, aerosol effects and snow/sand mixing can also be rectified by introducing thermal masks [51]. The reflection in different MODIS bands can also be used to differentiate clouds from other Earth features [52] and the MODIS product (MOD021km) can also be used to detect clouds using thermal bands [53].
Several researchers have applied MODIS and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER)-based NDSI to estimate the snow cover changes over different parts and sub-catchments of the Indus basin [8,[54][55][56]. The Quantitative spatiotemporal knowledge on snow cover dynamics is lacking for the entire Indus basin, which influences the freshwater availability downstream. Keeping in view the mentioned reasons and all the possible uncertainties in the snow cover estimations, an innovative and comprehensive approach is discussed in this study to quantify the decadal snow cover changes over the entire transboundary Indus basin. The objective of this study is to establish a new NDSI threshold for snow cover estimations using coarser resolution MODIS products. Moreover, the study also aims to ascertain the MODIS-based decadal snow cover changes over the Indus basin and its sub-catchments from 2008 to 2018. The uncertainties due to dark forests, water pixels, barren land (sand), clouds and snow-forest confusion were minimized in the MODIS products-based snow cover estimations. The Landsat ETM+, in situ snow data and Google Earth imagery were used to fine tune the MODIS NDSI value for snow cover estimations. The MODIS daily (MOD09GA) and 8-day (MOD09A1) products were further used to observe the snow cover trends using regression analysis over the Indus basin and its sub-catchments. The MODIS-based snow cover estimations, especially the daily MOD09GA product, can help researchers and water managers to estimate the glacial reserves and hydrological response of snow-dominant rivers' basins across the globe more efficiently.

Study Area
The current study is focused on the Indus basin and its western rivers' catchments, covering the spatial domain of 23.6 • to 37.4 • N and 64.8 • to 83.6 • E, including Pakistan, India, China and Afghanistan, as shown in Figure 1. The study area is monsoon dominated and precipitation varies greatly over space and time for the Indus basin. All six major rivers of the Indus basin originate from snow dominant mountainous areas and flow in southwest directions. The flow of the three western rivers of the basin (the Indus, Jhelum and Chenab) are under the control of Pakistan, while the three eastern rivers (the Ravi, Sutlej and Bias) are under Indian administrative control. The ASTER Global Digital Elevation Model (GDEM) was used for the elevation mapping of the Indus basin, which ranges from sea level to 8611 m and comprises several major peaks of the world. The 13 meteorological stations of the Pakistan Meteorological Department (PMD), are also indicated for reference in situ snow data. The Indus basin receives water from both rainfall and snowmelt and flows downstream to the Arabian Sea. The Indus basin experiences higher temperatures in its lower and middle parts and very low temperatures in the upper rivers' catchments.  MODIS Aqua and Terra satellites are designed to observe the Earth's biosphere as a component of NASA's long-term Earth Observing System. The major land features being observed by these two sensors are vegetation, water, ice and snow cover, surface albedo, temperature, emissivity, etc., almost on daily basis. MODIS has the capability to observe the whole Earth within 1-2 days with 36 spectral bands ranging from 0.4 to 14.4 µm. MODIS has spatial resolutions of 250, 500 and 1000 m in 1-2, 3-7 and 8-36 bands, respectively, with 12-bit radiometric sensitivity. The MODIS MOD09GA daily Level 2G Global data (Table 1) are used in the current study to estimate snow cover area, as they provide daily surface reflectance at a 500-m spatial resolution in the visible, NIR and Shortwave Infrared (SWIR) regions of the electromagnetic spectrum [57]. MOD09A1 is an advanced level 3G 8day composite product, which is also used at 500 m, with the same characteristics of reflectance bands as MOD09GA. The Landsat ETM+ is a multispectral scanning radiometer which offers datasets of the Earth's surface in visible, near, shortwave and longwave infrared with a 30-m spatial resolution (Table 1).  MODIS Aqua and Terra satellites are designed to observe the Earth's biosphere as a component of NASA's long-term Earth Observing System. The major land features being observed by these two sensors are vegetation, water, ice and snow cover, surface albedo, temperature, emissivity, etc., almost on daily basis. MODIS has the capability to observe the whole Earth within 1-2 days with 36 spectral bands ranging from 0.4 to 14.4 µm. MODIS has spatial resolutions of 250, 500 and 1000 m in 1-2, 3-7 and 8-36 bands, respectively, with 12-bit radiometric sensitivity. The MODIS MOD09GA daily Level 2G Global data (Table 1) are used in the current study to estimate snow cover area, as they provide daily surface reflectance at a 500-m spatial resolution in the visible, NIR and Shortwave Infrared (SWIR) regions of the electromagnetic spectrum [57]. MOD09A1 is an advanced level 3G 8-day composite product, which is also used at 500 m, with the same characteristics of reflectance bands as MOD09GA. The Landsat ETM+ is a multispectral scanning radiometer which offers datasets of the Earth's surface in visible, near, shortwave and longwave infrared with a 30-m spatial resolution (Table 1).  (1)) [58], which is the normalization ratio between green and SWIR 1 bands, was applied on both MODIS products on coinciding and mostly cloud-free days to estimate snow cover. As the spectral signature of snow reveals a greater response to visible wavelengths (green band) and a very small response to SWIR 1 [36], the normalization can be applied for snow cover delineation.
For NDSI values ranging from −1 to +1, higher values represent snow and the MODIS algorithm fixes the NDSI values higher than 0.4 as snow [36], but this threshold can reduce the snow cover significantly [59]. The 120 scenes of finer resolution Landsat ETM+ were used to estimate the snow cover area using the globally recommended threshold (NDSI ≥ 0.4). The MODIS products were tested with different NDSI thresholds ranging from 0.3 to 0.4 to map the snow cover over low-altitude glaciers due to its coarser resolution.

Water and Forest Detection
Uncertainties may be involved in the analysis due to the presence of dense forests and water response relatively similar to that of snow. Therefore, it is necessary to check the other Earth features and atmospheric particles in several wavelengths to minimize the chances of errors in the analysis. In the current study, the water pixels were identified by two different approaches. In the first approach, the absorption of NIR 1 radiation in water was used as an algorithm. In NIR 1 , the threshold value of reflection was selected as 11% and higher reflectance regions could not be considered as water. The NDWI, initially proposed by McFeeters et al. [60], was also applied for water body separation from snow cover and other Earth features using green and NIR 1 bands [61].
The NDWI ranges from −1 to +1, where values greater than 0.0 represent water bodies [62]. The blended water mask was prepared using NIR 1 and NDWI threshold and a comparison was performed with Google Earth imagery.
In the case of black dense/spruce forests, the denominator in NDSI remains quite low due to the very low reflectance in SWIR 1 and any small increase in the visible (green) band may increase the NDSI value to misclassify that pixel as snow. If the green band shows a reflection of <10% then the pixel is not classified as snow even if all the other criteria are justified [36]. To account for the snow cover in forests, Klein et al. [45] integrated the canopy and snow reflectance models to map pixels as snow in forests.
For this purpose, NDVI and NDSI were used at the same time to delineate snow cover, even in dark forests. Dense green forests show a higher reflectance in the NIR 1 band and lower reflectance in the red band, which ultimately assigns higher NDVI values to dense green areas, while, if there is snow, then the reflectance in NIR 1 decreases and thereby causes a decrease in NDVI values. Thus, a pixel can be considered as snow even if its NDVI value ≈ 0.1 and its NDSI is <0.4 [63].

Clouds Detection and Snow/Sand Confusion in MODIS Products
The MOD021km was used along with MOD09GA and MOD09A1 to detect clouds in reflectance and thermal bands due to its advantage of a wide range of 36 spectral bands (visible, infrared and thermal bands). MOD021km contains reflectance and radiance at a spatial resolution of 1000 m for both reflective and emissive bands in Wm −2 µm −1 sr −1 . The spectral signature/contrast analysis of clouds and other Earth features in different bands can be helpful to differentiate clouds from all other Earth features [64]. For cloud detection, a combination of two reflectance bands (blue and SWIR 1 ) for MODIS products were used in this study. A reflectance value greater than 0.2 in the SWIR 1 band is suitable to detect clouds [65] and a range of 0.2-0.3 can also be suitable for this purpose [66]. A minimum threshold value of 0.07 was tested for reflectance in the SWIR 2 band to detect clouds and soil, but an error was involved due to the slightly higher reflection of sandy soil. The maximum and minimum reflectance used in the blue band were 0.3 and 0.37 to detect barren (sandy) soil and clouds, respectively, for MOD09GA and MOD09A1 products.
The MOD021km was geo-registered in the ENVI for ArcGIS software package for the conversion of data into reflectance, and atmospheric corrections were applied using the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercube (FLASH) method [67]. Cloud detection using MOD021km has been applied by different researchers over mountainous regions across the world [68,69]. MOD021km provides reflectance bands from 1 to 19 and 26 (NIR 2 ) and thermal bands range from 20 to 36, excluding band 26. The MOD021km reflectance in SWIR 2 and NIR 2 bands at a 1000-m resolution were used as another algorithm to resolve snow/sand confusion and cloud detection, respectively. The minimum reflectance thresholds of 0.07 and 0.1 were used for SWIR 2 and NIR 2 bands to detect sandy soil and clouds, respectively.
Another approach used to detect clouds is based on the thermal bands of MOD021km, as cloud detection has been performed by different researchers using thermal bands [70,71]. The radiance of all 16 emissive bands was converted to brightness temperature using the following equation: where L is the radiance in Wm −2 µm −1 sr −1 , h is Plank's constant (6.626 × 10 −34 Js), c is the speed of light (3 × 10 8 m/s), k is the Boltzmann gas constant (1.381 × 10 −23 J/K), λ is the central wavelength of the band (µm) and T is the brightness temperature in Kelvin (K). As the temperature of the cloud is less than the other Earth's features (the base of the algorithm), each thermal band was assigned a maximum brightness temperature threshold to detect clouds. Finally, a blended cloud mask was prepared using the MODIS blue, NIR 2 , SWIR 1 , SWIR 2 and thermal bands.

Snow Cover Mapping
The snow cover area of Landsat ETM+, in situ snow data and Google Earth imagery were used as reference data for MODIS products. Initially, Google Earth was based on Landsat data with the Shuttle Radar Topographic Mission digital elevation model [72]. Google Earth uses Landsat data, Orthophotos by the State Government, high-resolution datasets available from DigitalGlobe Thematic Mapper (TM) , GeoEye TM , Worldview TM (1 and 2), SPOT TM , FORMOSAT-2, KOMPSAT-2, Pleiades and IKONOS [73]. Google Earth high-resolution reference data were used for water and snow indices by Yan et al. [74] in the Tibetan Plateau. The rock glaciers were mapped using Google Earth imagery by Schmid et al. [75] over the Hindu Kush Himalayan region. The weekly snow data of 13 stations of PMD were taken as training sites which lie in the snow-covered regions of Pakistan. The Landsat ETM+ data were used as standard data to delineate the snow cover area of 120 cloud-free scenes for different time periods and different regions over the Indus basin from 2008 to 2018. The Landsat ETM+-based water and forest areas were delineated and removed from the snow cover area estimated by the Landsat ETM+-based NDSI (≥0.4). In situ snow data and Google Earth images were used to cross check the snow cover delineation by Landsat ETM+.
All of the possible uncertainties, like water, forests and clouds, were removed and snow-sand confusion was also resolved in the snow cover estimated by MODIS MOD09GA and MOD09A1 products. The NDSI threshold fixation of a coarser resolution product like MODIS is of significant importance for the snow cover lying over low-to mid-altitude glaciers. For this purpose, the NDSI values (0.3 to 0.4) below the globally accepted threshold (0.4) were tested to estimate the snow cover area using MODIS MOD09GA and MOD09A1 products for 120 scenes (coinciding with 120 Landsat ETM+ snow cover scenes). The output of MOD09GA and MOD09A1 was compared with the Landsat ETM+, in situ snow data and Google Earth imagery for accuracy assessments. These comparisons were also made to establish a new NDSI threshold for snow cover estimations using coarser resolution MODIS products. The Google Earth imagery along with the in situ data of 13 PMD stations were also used to ascertain the snow cover, estimated by MODIS MOD09GA and MOD09A1 products, over the low-to mid-altitude mountains of the Indus basin.
The MODIS-based snow cover analysis from 2008 to 2018 was performed mostly for cloud-free days. However, snow cover analysis for a few cloudy days was necessary due to the continuous cloud cover during the whole month. The MODIS cloud masks were developed using multiple algorithms and a blended cloud mask was prepared by combining all the cloud masks for an individual day over the Indus basin. The blended cloud mask was then removed from the NDSI-based snow cover area in order to minimize the overestimation of the snow cover area due to the mixing of clouds in the snow, but this removal also removed some of the snow cover areas. During the processing of MODIS cloud-free days, snow cover was monitored, using Google Earth imagery and in situ snow data, and minimum ASTER GDEM-based elevations were recorded against the presence or absence of snow in all parts of the Indus basin. Afterwards, the snow cover was declared above those minimum elevation thresholds under cloudy regions for the respective month of the same year.
The time series intercomparison of snow cover area of MODIS daily (MOD09GA and 8-day (MOD09A1) products was carried out using R 2 and Nash and Sutcliffe Efficiency (NSE) for the Indus basin and its sub-catchments from 2008 to 2018. The NSE indicates how well the plot of observed versus simulated data fits the 1:1 line. In this study, NSE was used to compare the snow cover area of both the MODIS products with the Landsat ETM+ based snow cover estimations for 120 selected scenes over the Indus basin. NSE was also used for the intercomparison of snow cover areas of both MOD09GA and MOD09A1 products in order to check the closeness of their estimations. A least squares regression analysis was performed with a 95% confidence interval using both the MODIS products for a time series snow cover trend analysis for the Indus basin and its sub-catchments.

Uncertainty Analysis in Snow Cover Estimates
The water, forest and cloud masks are presented in Figure 2a-d. A comparison of the blended water mask with Google Earth imagery is presented in Figure 2a,b. The MODIS and Google Earth-based Tarbela and Mangla dams of Pakistan are presented for comparison. The comparison indicates that the overall pattern of both the lakes is good, while there are some disturbances on the borders of lakes, owing to the coarser resolution of the MODIS products. The forest and vegetation mask is presented in Figure 2c, indicating forests mainly in the mid to high altitudes, and wheat crop mainly in the plain areas of India and Pakistan. One of the main challenges was to differentiate the snow from snow-forest area. The NDVI (≈0.1) and NDSI (<0.4) were used in this study to represent snow which was added to the NDSI analysis, and the remaining forest pixels were added to the forest and vegetation mask. The cloud mask is presented in Figure 2d for a specific day, indicating most of the clouds over the Indus river catchment. The details of snow-sand confusion and the cloud detection algorithm, using reflection and thermal bands, are presented in Table 2. Clouds and barren soil (sand) were effectively detected over the Indus basin using the reflection bands of MOD09GA, MOD09A1 and MOD021km. The clouds' brightness temperature helps us to differentiate the clouds from other features of the Earth. The thermal bands of MOD021km, from 20 to 36 (except 26), were used in the analysis to detect clouds over different wavelengths of the electromagnetic spectrum. The mean brightness temperature of clouds indicates that the brightness temperature increases from lower wavelengths to higher wavelengths. The new thermal thresholds of clouds' brightness temperatures were established, as provided in Table 2, indicating a non-linear trend, with maximum (258 K) and minimum (225 K) values for thermal bands 20 and 36, respectively. The clouds' brightness temperature was tested for different time periods over the Indus basin to fix one value for each thermal band. All the thermal masks were applied, especially the higher wavelength bands, due to the very low clouds' brightness temperatures , to make a final blended cloud mask using thermal and reflection bands. The removal of the cloud mask from the NDSI-based snow cover estimation introduced a significant amount of uncertainty into the analysis due to the underestimation by removing the snow, while the snow cover analysis without cloud removal caused an overestimation of snow due to the misclassification of clouds as snow pixels. Most of the rivers' catchments of the Indus basin are snow dominant and observe both precipitation and anvil clouds. This was the main reason why we used ASTER GDEMbased elevation information, the clouds were considered as snow above the minimum elevation for the presence of snow, fixed for the respective area during the analysis of cloud-free days of the same month. The elevation data of 13 in situ PMD stations were also used in combination with ASTER GDEM-based snow elevation thresholds to decide the snow or snow-free zone below the clouds.  The details of snow-sand confusion and the cloud detection algorithm, using reflection and thermal bands, are presented in Table 2. Clouds and barren soil (sand) were effectively detected over the Indus basin using the reflection bands of MOD09GA, MOD09A1 and MOD021km. The clouds' brightness temperature helps us to differentiate the clouds from other features of the Earth. The thermal bands of MOD021km, from 20 to 36 (except 26), were used in the analysis to detect clouds over different wavelengths of the electromagnetic spectrum. The mean brightness temperature of clouds indicates that the brightness temperature increases from lower wavelengths to higher wavelengths. The new thermal thresholds of clouds' brightness temperatures were established, as provided in Table 2, indicating a non-linear trend, with maximum (258 K) and minimum (225 K) values for thermal bands 20 and 36, respectively. The clouds' brightness temperature was tested for different time periods over the Indus basin to fix one value for each thermal band. All the thermal masks were applied, especially the higher wavelength bands, due to the very low clouds' brightness temperatures, to make a final blended cloud mask using thermal and reflection bands. The removal of the cloud mask from the NDSI-based snow cover estimation introduced a significant amount of uncertainty into the analysis due to the underestimation by removing the snow, while the snow cover analysis without cloud removal caused an overestimation of snow due to the misclassification of clouds as snow pixels. Most of the rivers' catchments of the Indus basin are snow dominant and observe both precipitation and anvil clouds. This was the main reason why we used ASTER GDEM-based elevation information, the clouds were considered as snow above the minimum elevation for the presence of snow, fixed for the respective area during the analysis of cloud-free days of the same month. The elevation data of 13 in situ PMD stations were also used in combination with ASTER GDEM-based snow elevation thresholds to decide the snow or snow-free zone below the clouds.

NDSI Threshold and Snow Cover Estimations
A comparison of the snow cover area of MOD09GA and MOD09A1 with two scenes of Landsat ETM+ coarsened to a 500-m resolution is presented in Figure 3. Scenes 1 and 2 on the left and right side, respectively, in Figure 3 represent the snow cover of Landsat ETM+ and both the MODIS (NDSI ≥ 0.34) products. In the first scene, the comparison is good and an underestimation in snow cover is observed by both the MODIS products, as highlighted by the red circles on the Landsat ETM+ scene. In the second scene, the underestimation is also clear in both the MODIS products compared to the Landsat ETM+ scene, as indicated by the red circles.

NDSI Threshold and Snow Cover Estimations
A comparison of the snow cover area of MOD09GA and MOD09A1 with two scenes of Landsat ETM+ coarsened to a 500-m resolution is presented in Figure 3. Scenes 1 and 2 on the left and right side, respectively, in Figure 3 represent the snow cover of Landsat ETM+ and both the MODIS (NDSI ≥ 0.34) products. In the first scene, the comparison is good and an underestimation in snow cover is observed by both the MODIS products, as highlighted by the red circles on the Landsat ETM+ scene. In the second scene, the underestimation is also clear in both the MODIS products compared to the Landsat ETM+ scene, as indicated by the red circles. A statistical comparison of the snow cover area of both the MODIS products with the Landsat ETM+, in situ snow data and Google Earth imagery for 120 scenes at different times is presented in Table 3. The comparison indicates that a minimum NDSI threshold of 0.34 has a good R 2 and a higher value of NSE. The comparison in Figure 3 also reveals that snow is removed in the MODIS MOD09GA and MOD09A1 products over some of the low-altitude glaciers at NDSI ≥ 0.34. Moreover, A statistical comparison of the snow cover area of both the MODIS products with the Landsat ETM+, in situ snow data and Google Earth imagery for 120 scenes at different times is presented in Table 3. The comparison indicates that a minimum NDSI threshold of 0.34 has a good R 2 and a higher value of NSE. The comparison in Figure 3 also reveals that snow is removed in the MODIS MOD09GA and MOD09A1 products over some of the low-altitude glaciers at NDSI ≥ 0.34. Moreover, at NDSI > 0.4, some of the mid-altitude glaciers were misclassified as non-snow pixels, so the threshold of 0.34 was preferred in this study for both the MODIS products, as a lower threshold value (NDSI > 0) has also been tested by Dong et al. [76]. The NDSI thresholds below 0.34 revealed that snow was also observed in plain and desert areas, which was leading to an overestimation of snow cover. The possible reasons for the underestimation of snow cover, at the newly established threshold (NDSI ≥ 0.34), are the coarser resolution of MODIS products and changes in the wavelengths of bands compared to Landsat ETM+. The snow cover delineation patterns derived from both MOD09GA and MOD09A1 products are presented in Figure 4a  The snow cover delineation patterns derived from both MOD09GA and MOD09A1 products are presented in Figure 4a   On 5 March 2008, both the MODIS products show a difference in their snow cover patterns, mainly on the western side (the Kabul river basin), northern side and the eastern rivers' catchments of the basin. On 26 February 2018, both the MODIS products indicate almost similar snow cover patterns over the Indus basin. The snow cover patterns of both the MODIS products over the Indus basin indicate that the snow cover on the western side of the Indus basin and the eastern rivers' catchments is decreased. The snow cover on the low-altitude mountains is also significantly decreased over the Indus basin. The low-to mid-altitude snow cover of the western rivers' catchments also revealed a significant decrease from 2008 to 2018.
The snow cover of daily MOD09GA product is presented in Figure 5a-f to account for the differences in delineated snow cover patterns over the western rivers' catchments of the Indus basin. The snow cover pattern of the Chenab river catchment indicates that the snow cover on the upper part of the catchment is preserved, while in the lower reaches the snow cover is reduced. The snow cover decline in the upper reaches is observed on the low-altitude mountains and on both the sides of the main Chenab river and its tributaries. The Jhelum river catchment showed a serious change in the snow cover pattern, the snow cover on low-to mid-altitude mountains and on the eastern side of the Jhelum river catchment revealed a significant decline. The western and northern parts of the Indus river catchment are affected badly and a change in snow cover pattern has been observed in the low-altitude parts of the catchment. Our MODIS products-based snow cover pattern analysis revealed that changes in snow cover patterns have been observed on the low-to mid-altitude mountains of the Indus basin and its sub-catchments.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 23 catchments is decreased. The snow cover on the low-altitude mountains is also significantly decreased over the Indus basin. The low-to mid-altitude snow cover of the western rivers' catchments also revealed a significant decrease from 2008 to 2018. The snow cover of daily MOD09GA product is presented in Figure 5a-f to account for the differences in delineated snow cover patterns over the western rivers' catchments of the Indus basin. The snow cover pattern of the Chenab river catchment indicates that the snow cover on the upper part of the catchment is preserved, while in the lower reaches the snow cover is reduced. The snow cover decline in the upper reaches is observed on the low-altitude mountains and on both the sides of the main Chenab river and its tributaries. The Jhelum river catchment showed a serious change in the snow cover pattern, the snow cover on low-to mid-altitude mountains and on the eastern side of the Jhelum river catchment revealed a significant decline. The western and northern parts of the Indus river catchment are affected badly and a change in snow cover pattern has been observed in the low-altitude parts of the catchment. Our MODIS products-based snow cover pattern analysis revealed that changes in snow cover patterns have been observed on the low-to mid-altitude mountains of the Indus basin and its sub-catchments.

Annual Snow Cover Trends
The annual snow cover trends of the Indus basin and its sub-catchments, on an annual basis, are presented in Figure 6 for both the MODIS MOD09GA and MOD09A1 products. It is observed that the snow cover area for both the products follows an almost similar trend throughout the individual year for the Indus basin and all of its rivers' catchments. It is obvious that the major contribution to the snow cover area of the Indus basin is from the Indus river catchment compared to the other rivers' catchments. Decreasing trends of snow cover are observed from the middle of March to the middle of June, followed by constant trends up to the middle of August, and then rising trends are observed for all the rivers' catchments and the Indus basin. Over and underestimations in the snow cover area

Annual Snow Cover Trends
The annual snow cover trends of the Indus basin and its sub-catchments, on an annual basis, are presented in Figure 6 for both the MODIS MOD09GA and MOD09A1 products. It is observed that the snow cover area for both the products follows an almost similar trend throughout the individual year for the Indus basin and all of its rivers' catchments. It is obvious that the major contribution to the snow cover area of the Indus basin is from the Indus river catchment compared to the other rivers' catchments. Decreasing trends of snow cover are observed from the middle of March to the middle of June, followed by constant trends up to the middle of August, and then rising trends are observed for all the rivers' catchments and the Indus basin. Over and underestimations in the snow cover area of both the MODIS products are observed for all the rivers' catchments and the Indus basin in every year.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 23 Figure 6. Annual time series trends of snow cover area estimated using MODIS MOD09GA and MOD09A1 products for the Indus basin and its rivers' catchments.

Time Series Snow Cover Dynamics from 2008 to 2018
The long-term snow cover trends are presented in Figure 7a-e from 2008 to 2018 for the Indus basin and its rivers' catchments. In the case of the Chenab river catchment, the maximum snow cover is observed for MOD09GA and MOD09A1 in 2011 and 2009, respectively, with very small differences in peaks. As most of the snow-dominant area of the Chenab river catchment lies higher than 4000 m Above Mean Sea Level (AMSL), the decrease in snow cover is not too much due to the presence of snow mostly at higher altitudes. A drastic decrease in the maximum snow cover over the Chenab river catchment is observed in 2015 and a slightly decreasing trend is observed up to 2018. Over the Jhelum river catchment, the maximum snow cover is recorded by MOD09GA and MOD09A1 in 2011 and 2010, respectively and a severe decrease is observed in 2014, which further

Time Series Snow Cover Dynamics from 2008 to 2018
The long-term snow cover trends are presented in Figure 7a-e from 2008 to 2018 for the Indus basin and its rivers' catchments. In the case of the Chenab river catchment, the maximum snow cover is observed for MOD09GA and MOD09A1 in 2011 and 2009, respectively, with very small differences in peaks. As most of the snow-dominant area of the Chenab river catchment lies higher than 4000 m Above Mean Sea Level (AMSL), the decrease in snow cover is not too much due to the presence of snow mostly at higher altitudes. A drastic decrease in the maximum snow cover over the Chenab river catchment is observed in 2015 and a slightly decreasing trend is observed up to 2018. A decreasing trend of snow cover area is observed from 2008 to 2018 over the eastern parts of the Indus basin, consisting of the catchments of all three major eastern rivers. The maximum snow cover area is recorded at the start of 2008 for both the MODIS products and a severe decline in the maximum snow cover area is observed at the start of 2015. Most of the snow cover area of the eastern rivers' catchments lies below 4500 m AMSL and some of it even lies below 4000 m, hence faced degradation of snow cover due to global warming and climate change. A drastic decrease in the maximum snow cover area is observed after 2012 and the decreasing trend continued, while a small increase in yearly maximum snow cover is observed in 2018.
The monitoring of snow cover dynamics is of great importance over the Indus river catchment as Indus river flow is mainly dependent on snowmelt compared to runoff produced by rainfall. Most of the snow-dominated area of the Indus river catchment has an elevation below 4500 m AMSL that experiences snowfall and snowmelt in each individual year. Due to the dynamic elevation range of the Indus river catchment, snow cover monitoring poses serious challenges to the researchers. The maximum snow cover area is observed in 2008 for Indus river catchment and a variable trend is observed up to 2018. A severe decrease in the maximum snow cover area is observed in the first quarters of 2011 and 2013. Over and underestimations of both the MODIS products are also observed from 2008 to 2018 for Indus river catchment. For the Indus basin, the maximum and minimum snow cover area is observed at the start of 2008 and 2013, respectively, and over and underestimations of A decreasing trend of snow cover area is observed from 2008 to 2018 over the eastern parts of the Indus basin, consisting of the catchments of all three major eastern rivers. The maximum snow cover area is recorded at the start of 2008 for both the MODIS products and a severe decline in the maximum snow cover area is observed at the start of 2015.
Most of the snow cover area of the eastern rivers' catchments lies below 4500 m AMSL and some of it even lies below 4000 m, hence faced degradation of snow cover due to global warming and climate change. A drastic decrease in the maximum snow cover area is observed after 2012 and the decreasing trend continued, while a small increase in yearly maximum snow cover is observed in 2018.
The monitoring of snow cover dynamics is of great importance over the Indus river catchment as Indus river flow is mainly dependent on snowmelt compared to runoff produced by rainfall. Most of the snow-dominated area of the Indus river catchment has an elevation below 4500 m AMSL that experiences snowfall and snowmelt in each individual year. Due to the dynamic elevation range of the Indus river catchment, snow cover monitoring poses serious challenges to the researchers. The maximum snow cover area is observed in 2008 for Indus river catchment and a variable trend is observed up to 2018. A severe decrease in the maximum snow cover area is observed in the first quarters of 2011 and 2013. Over and underestimations of both the MODIS products are also observed from 2008 to 2018 for Indus river catchment. For the Indus basin, the maximum and minimum snow cover area is observed at the start of 2008 and 2013, respectively, and over and underestimations of snow cover for both the MODIS products are observed from 2008 to 2018. All the rivers' catchments and the Indus basin showed a decline in the snow cover area from 2008 to 2018 causing a loss of reserved frozen water resources, which were present in form of snowpack.

MODIS Time Series Snow Cover Trends
A least squares regression analysis was performed with a 95% confidence interval to discuss the maximum snow cover trends of both the MODIS products over the Indus basin and its sub-catchments from 2008 to 2018, as shown in Figure 8a-j.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 23 snow cover for both the MODIS products are observed from 2008 to 2018. All the rivers' catchments and the Indus basin showed a decline in the snow cover area from 2008 to 2018 causing a loss of reserved frozen water resources, which were present in form of snowpack.

MODIS Time Series Snow Cover Trends
A least squares regression analysis was performed with a 95% confidence interval to discuss the maximum snow cover trends of both the MODIS products over the Indus basin and its subcatchments from 2008 to 2018, as shown in Figure 8a-j. The percentage difference between the maximum snow cover area of MOD09GA with respect to MOD09A1 is presented in Table 4. In most of the observations, the percentage difference is less than 5%, representing a reliable agreement, while greater differences are observed for the year 2009 and from 2013 to 2017 for the Chenab and Jhelum rivers' catchments. Low percentage differences are observed for the Indus river catchment, except for the year 2012, with a difference of −6.37%. The eastern rivers' catchments record a less than five percent difference, except for the years 2008, 2010, 2011, 2012 and 2018. The differences in yearly maximum snow cover are mostly recorded as less than 5% for the Indus basin, represent reliable agreements between MODIS products

Discussions
Snow cover estimation is very important for the planning, management and mitigation of biophysical processes. The MODIS products used in this study are preferred for snow cover mapping compared to other products like the National Operational Hydrologic Remote Sensing Center, as MODIS detects more snow over higher altitudes and its algorithm is not static but evolving [77]. The National Environmental Satellite Data and Information Service snow product of the National Oceanic and Atmospheric Administration overestimates the snow cover over higher latitudes [78]. MODIS also offers regular quantifications of snow cover with respect to data availability and spatiotemporal resolution. [79,80]. The accuracies of MODIS and the National Operational Hydrologic Remote Sensing Center snow products have been observed as 94 and 76%, respectively, over the Upper Rio Grande Basin [81]. The MODIS 8-day snow cover product (MOD10A2/MYD10A2) is the combination of 2 to 8 days of MODIS Aqua/Terra daily snow product (MOD10A1/MYD10A1). The MODIS MOD10A1 and MYD10A1 products are being used to produce multi-day snow cover products with fixed temporal windows [82][83][84][85]. Such approaches can reduce the cloud coverage, but products lack the ability to observe snowfall events to some extent due to the fixation of temporal windows [86]. The daily and 8-day composite products used in this study offer a reliable temporal resolution, but in poor weather conditions the composite products produce images of reliable quality after a gap of 8 days to several weeks. The MOD10_L2 can be analyzed for land or inland water pixels only, which are unobstructed by clouds in the daylight, for the presence or absence of snow, and coarser resolution MODIS products like MOD10C1 and MOD10C1C highly overestimate the snow systematically [87].

Discussions
Snow cover estimation is very important for the planning, management and mitigation of bio-physical processes. The MODIS products used in this study are preferred for snow cover mapping compared to other products like the National Operational Hydrologic Remote Sensing Center, as MODIS detects more snow over higher altitudes and its algorithm is not static but evolving [77]. The National Environmental Satellite Data and Information Service snow product of the National Oceanic and Atmospheric Administration overestimates the snow cover over higher latitudes [78]. MODIS also offers regular quantifications of snow cover with respect to data availability and spatiotemporal resolution. [79,80]. The accuracies of MODIS and the National Operational Hydrologic Remote Sensing Center snow products have been observed as 94 and 76%, respectively, over the Upper Rio Grande Basin [81]. The MODIS 8-day snow cover product (MOD10A2/MYD10A2) is the combination of 2 to 8 days of MODIS Aqua/Terra daily snow product (MOD10A1/MYD10A1). The MODIS MOD10A1 and MYD10A1 products are being used to produce multi-day snow cover products with fixed temporal windows [82][83][84][85]. Such approaches can reduce the cloud coverage, but products lack the ability to observe snowfall events to some extent due to the fixation of temporal windows [86]. The daily and 8-day composite products used in this study offer a reliable temporal resolution, but in poor weather conditions the composite products produce images of reliable quality after a gap of 8 days to several weeks. The MOD10_L2 can be analyzed for land or inland water pixels only, which are unobstructed by clouds in the daylight, for the presence or absence of snow, and coarser resolution MODIS products like MOD10C1 and MOD10C1C highly overestimate the snow systematically [87].
The MODIS binary snow products map the snow along with other major Earth features with some sets of rules and the codes do not contain valid land cover state information like clouds, missing data, etc. The MODIS NDSI-based fractional snow cover products represent snow pixels within a MODIS pixel using some empirical equations and in situ snow information, and these products also lack valid land cover state information. Fractional snow cover is no longer provided in version 6 products, and previous binary snow cover has been replaced by NDSI snow cover (equation 1) [88]. The most widely used indicator for snow cover estimations is NDSI, which requires deeper investigations for greater accuracy. The fixation of the NDSI threshold is very important for coarser resolution products, as a significant amount of snow can be misclassified as non-snow pixels at a minimum globally accepted NDSI threshold of 0.4 [32,89]. The results of our snow cover comparison (Table 3) indicate that higher values of NSE and R 2 are obtained at a minimum NDSI threshold of 0.34 for MODIS MOD09GA and MOD09A1 products when compared with finer resolution Landsat ETM+ snow cover, in situ snow data and Google Earth imagery. We also tested the NDSI values less than 0.34 in this study, as suggested by Huang et al. [88] and Shimamura et al. [90], which resulted in misclassification of many forests, water and barren lands (sand) as snow pixels.
The atmospheric disturbance, spatial resolution and the subsequent mixing of reflectance of other Earth features with that of snow, within the pixel, increase the probability of uncertainty. The water, forest, barren lands (sand) and clouds cause uncertainties in snow cover estimations and were rectified by a separate analysis in this study. The resemblance between water and snow poses a serious challenge for differentiation due to their similar properties and one threshold is not sufficient for the detection of water [91]. The blended water body mask can be developed using combined NIR 1 [92] and NDWI-based [60] thresholds. We prepared a blended water mask and the results are satisfactory over the Indus basin when compared with Google Earth imagery of Tarbela and Mangla dams of Pakistan. The dark forest areas can be misclassified as snow pixels [45,93], which may lead to poor snow cover estimations, and the visible (green) band threshold [36] played a key role in this study to differentiate the dark forests from the snow. The snow-forest differentiation is also a serious concern in snow cover estimations [27,94] as it can lead to underestimations of snow cover if snow-forest pixels are classified as forest. The snow-forest confusion was resolved using NDSI and NDVI thresholds [63,95] at the same time in this study, which proved to be very helpful for improving the accuracy of snow cover estimations. The clouds over snow-covered mountains also lead to the overestimation of snow cover [96] and cloud detection using reflection and thermal bands has been applied by several researchers worldwide [97][98][99][100]. We used both the thresholds of reflection and thermal bands for development of a blended cloud product for each individual cloudy day and ASTER GDEM-based month-wise snow elevation thresholds proved to be very helpful to classify snow and non-snow pixels for cloudy days. The analysis of thermal bands indicated that the mean clouds' brightness temperature showed an increasing trend towards longer wavelengths, which is also observed by Ahmad et al. [53]. The intercomparison of MODIS products in Figure 8a-e also indicates reliable agreements, as various researchers have also validated the overall accuracy, ranging from 85 to 90 percent for MODIS snow cover around the globe [46,81,87].
The study reveals a declining trend of about 458, 382, 840, 463 and 2459 km 2 per year for the Chenab, Jhelum, Indus and eastern rivers' catchments and Indus basin, respectively. Such high declining rates of snow cover may cause the smaller rivers of the Indus basin to become dry in the near future, which may be a serious threat to the freshwater supplies of the Indus basin. A decline in the snow cover has been observed by researchers over the Himalayas [8,101] and the upper Indus basin [102,103], mainly due to an increase in temperature [104,105]. The Global Land Ice Measurements from Space (GLIMS) project also revealed a consistent decline in the snow cover extent in the upper Indus basin [106]. Decreasing snow cover trends in the catchments of western rivers of the Indus basin during winter and autumn periods have been observed by Hasson et al. [56]. Snow cover decline due to climate change has been observed by researchers in the Sutlej river basin [107] and Kashmir Himalayas [108,109]. The IPCC report concluded that the warming of land is expected to be the highest in higher northern latitudes, resulting in a severe decline in the snow cover [110]. The results of this study draw attention to an astonishing fact: that the glaciers are melting at a very high rate under the prevailing climatic conditions and, if this rate continues, the Indus basin may face serious snow cover decline within a century. The IPCC report also unveiled that the snow cover on the Himalayas is receding at faster rates than the other parts of the world and the thinning of the Himalayan glaciers due to climate change may cause the Indus, Ganges and Brahmaputra rivers to become seasonal rivers in the near future [15]. The need of time is to take serious measures to preserve the frozen water reserves by involving all the stakeholders in joint watershed management and climate change mitigation programs in all the parts of the Indus basin and the Himalayas.

Conclusions
The main goal of this study was to estimate the high temporal MODIS-based snow cover dynamics over the Indus basin and its sub-catchments from 2008 to 2018 by fixation of new NDSI threshold. We compared the snow cover area of MODIS MOD09GA and MOD09A1 products with snow cover of 120 scenes of Landsat ETM+, in situ snow information and Google Earth imagery to establish a new NDSI threshold for snow cover estimations for coarser resolution MODIS products. The comparison indicated that the minimum NDSI threshold of 0.34 is more appropriate with minimum R 2 and NSE values of 0.97 and 0.85, respectively. These NDSI tests were conducted by focusing on the snow cover over the low-to mid-altitude mountains of the study area. The minimization of uncertainties, using water, forest and cloud masks, provided a better understanding of the spatiotemporal behavior of the snow cover over complex and data-scarce mountainous rivers' catchments of the Indus basin. The intercomparison of MOD09GA and MOD09A1-based maximum snow cover areas indicated that the percentage difference for most of the estimations was less than 5%, which indicates the reliability of the snow cover estimations. A linear least squares regression analysis of both the MODIS products indicated a minimum snow cover declining trend of about 458, 382, 840, 463 and 2459 km 2 per year for the Chenab, Jhelum, Indus and eastern rivers' catchments and Indus basin, respectively. The MODIS 8-day MOD09A1 product indicated greater declining snow cover trends compared to the daily MOD09GA product for the Indus basin and its sub-catchments. The Indus basin experienced greater variations in snow cover trends for both the MODIS products, mainly due to variations in the snow cover of the Indus river catchment. The low-to mid-altitude mountainous regions of the Indus basin and its sub-catchments have been affected badly, and snow cover decline has been observed over all the parts of the Indus basin. It is believed that localized NDSI threshold fixation for coarser resolution products will provide new opportunities to the researchers to estimate the frozen water reserves over complex mountainous regions across the world.