Contemporary Snow Changes in the Karakoram Region Attributed to Improved MODIS Data between 2003 and 2018

: Snowmelt signiﬁcantly contributes to meltwater in most parts of High Mountain Asia. The Karakoram region is one of these densely glacierized and snow-covered regions. Several studies have reported that glaciers in the Karakoram region remained stable or experience slight mass loss. This trend has called for further investigation to understand changes in other components of the cryosphere. This study estimates the comparative snow cover area (SCA) and snowline altitude (SLA) changes between 2003 and 2018 in the Karakoram region and its subbasins, including Hunza, Shigar, and Shyok. We used three di ﬀ erent 8-day composite snow products of the Moderate Resolution Imaging Spectroradiometer (MODIS) in this study including (1) Original Aqua (MYD10A2), (2) Original Terra (MOD10A2), and (3) Improved Terra-Aqua (MOYDGL06*) snow products from 2003 to 2018. We used Mann–Kendall and Sen Slope methods to assess trends in the SCA and SLA. Our results show that the original snow products are signiﬁcantly biased when investigating seasonal and annual trends. However, discarding a cloud cover of > 20% in the original products improves the results and makes them more comparable to our improved snow product. The original products (without cloud removal) overestimate the SCA during summer and underestimate the SCA during winter and year-round throughout the Karakoram region. The bias in the mean annual SCA between improved and Aqua and Terra cloud threshold products for the Karakoram region is found to be − 1.67% and 1.1%, respectively. The improved (MOYDGL06*) product reveals a statistically insigniﬁcant decreasing trend of the SCA on the annual scale between 2003 and 2018 in the Karakoram region and all three subbasins. The annual trends decreased at − 0.13%, − 0.1%, − 0.08%, and − 0.05% in the Karakoram, Hunza, Shigar, and Shyok, respectively. The monthly trends were slightly positive overall in December. The annual maximum SLA shows a statistically signiﬁcant upward trend of 13 m above sea level (m a.s.l.) per year for the entire Karakoram region. This ﬁnding suggests a signiﬁcant uncertainty in water resource planning based on the original snow data, and this study recommends the use of the improved snow product for a better understanding. A similar pattern of underestimation is demonstrated by the mean annual maximum SLAs. Original and cloud threshold products provide nearly equal maximum annual SLAs throughout the study area. The annual minimum SLA shows lower values in the original and cloud products compared to the improved product. This highlights the underestimation of the SLA in all the original products even at the annual scale. The three basins in the Karakoram exhibit analogous behaviour. The annual SLA ﬂuctuates in the improved snow product compared to the original snow products. The annual snowline exhibits immense spatial variability in the Karakoram region, with mean annual maximum values of 4655, 5259, and 5145 m a.s.l. in MYD10A2, MOYDGL06*, and MOD10A2, respectively. The mean maximum annual SLA in the Shyok is the highest among the basins and the lowest SLA occurs in the Hunza. Similar patterns were observed in the mean annual SLA estimates. The snowline elevation shows a gradually decreasing trend from the eastern to western Karakoram. Speciﬁcally, the eastern Shyok basin has a higher snowline altitude than the central and western regions.


Introduction
The cryosphere is an important water resource in High Mountain Asia (HMA), and the region is often recognized as the water tower or the third pole because HMA constitutes the largest volume of ice outside the polar region [1]. Understanding snow dynamics is of paramount importance for one-fourth (approximately two billion people) of the global population, who are dependent on meltwater originating from the mountains [2]. Realizing this importance, the Global Climate Observation System (GCOS) recognizes snow as one of the essential climate variables (ECVs). Most studies have focused on glacier change assessments, although the majority of meltwater in most of HMA catchments is from

Study Area
The Karakoram region extends between the longitudes 72.44 • to 79.64 • E and latitudes 33.14 • to 37.47 • N, as shown in Figure 1, with an area of 104,010 km 2 ( Table 1). The altitude of the region ranges from 1236 to 8611 m above sea level (m a.s.l.) with a mean altitude of 4709 m a.s.l. Westerlies are the dominant climatic system in the region. The Hunza, Shigar and Shyok are high-altitude basins of the upper Indus Basin (UIB). More than 60% of the glaciers of the Indus basin are found in the Karakoram region [27] and are mostly fed by avalanches, redistribution by wind, and seasonal snow. Approximately 22% of the Karakoram region is covered by glaciers, of which 27.8%, 24.9%, and 39.03% are found in the Hunza, Shigar and Shyok basins, respectively. The Hunza River Basin is a northern subcatchment of the UIB located in the western Karakoram region of Pakistan. The Hunza river is one of the major tributaries of the Indus River irrigation system, and it contributes approximately 12% of the upper Indus flow upstream of Tarbela Dam [7]. It has a catchment area of 13,745 km 2 with an elevation range of 1450 to 7404 m a.s.l. The SCA reaches up to 80% in winter and reduces to 30% in summer [18]. The Shigar Basin is located in central Karakoram. This catchment occupies an area of 7050 km 2 . The total glacierized area of the catchment is 5701 km 2 . Some of the well-known and important glaciers in the catchment are the Biafo and Baltoro. The SCA reaches 88% in January and decreases to 10% in August [28]. The Shigar Basin receives approximately 67% of its annual precipitation in winter [29][30][31][32]. The Shyok Basin has a catchment area of 33,419 km 2 and is located in the eastern Karakoram. Snow remains a dominant cryospheric component in the basin [33]. The seasonal variation in the SCA over the Shyok River Basin ranges from 27% to 59% with the maximum in February and minimum in July [24].
Water 2020, 12, x FOR PEER REVIEW 3 of 20 of glaciers. We also investigated the potential of the improved product for studying snow climatology and its applications for glaciohydrological modeling.

Study Area
The Karakoram region extends between the longitudes 72.44° to 79.64° E and latitudes 33.14° to 37.47° N, as shown in Figure 1, with an area of 104,010 km 2 ( Table 1). The altitude of the region ranges from 1236 to 8611 m above sea level (m a.s.l.) with a mean altitude of 4709 m a.s.l. Westerlies are the dominant climatic system in the region. The Hunza, Shigar and Shyok are high-altitude basins of the upper Indus Basin (UIB). More than 60% of the glaciers of the Indus basin are found in the Karakoram region [27] and are mostly fed by avalanches, redistribution by wind, and seasonal snow. Approximately 22% of the Karakoram region is covered by glaciers, of which 27.8%, 24.9%, and 39.03% are found in the Hunza, Shigar and Shyok basins, respectively. The Hunza River Basin is a northern subcatchment of the UIB located in the western Karakoram region of Pakistan. The Hunza river is one of the major tributaries of the Indus River irrigation system, and it contributes approximately 12% of the upper Indus flow upstream of Tarbela Dam [7]. It has a catchment area of 13,745 km 2 with an elevation range of 1450 to 7404 m a.s.l. The SCA reaches up to 80% in winter and reduces to 30% in summer [18]. The Shigar Basin is located in central Karakoram. This catchment occupies an area of 7050 km 2 . The total glacierized area of the catchment is 5701 km 2 . Some of the well-known and important glaciers in the catchment are the Biafo and Baltoro. The SCA reaches 88% in January and decreases to 10% in August [28]. The Shigar Basin receives approximately 67% of its annual precipitation in winter [29][30][31][32]. The Shyok Basin has a catchment area of 33,419 km 2 and is located in the eastern Karakoram. Snow remains a dominant cryospheric component in the basin [33]. The seasonal variation in the SCA over the Shyok River Basin ranges from 27% to 59% with the maximum in February and minimum in July [24].

MODIS 8-Day Maximum Snow Cover Product
Standard 8-day maximum snow cover products from MODIS Terra (MOD10A2, V006) and Aqua (MYD10A2, V006) with a spatial resolution of 500 m were used to extract the SCA and SLA. These data sets were downloaded from the National Snow and Ice Data Center Distributed Data Archive (https://nsidc.org/) for the period between 2003 and 2018. The MODIS snow data is available for the period from 18 December 1999 and 4 May 2002 from through Terra and Aqua satellites, respectively. A comprehensive description of the MODIS 8-day snow products version 6 is available at The National Snow and Ice Data Center (NSIDC) [26].

Improved MODIS 8-Day Maximum Snow Cover Product
The MOD10A2 and MYD10A2 products were processed to remove clouds and reduced the average snow underestimation of 3.66% and merged to reduce the overestimation of 46% in the improved 8-day maximum snow cover product MOYDGL06* [12]. All these products were used in this paper to analyze snow cover trends and estimate snowline dynamics. The improved cloud-free snow cover product is freely available from the Pangaea repository [17]. The detailed snow data improvement methodology and dataset are described in [12].

Elevation Data
The snowline elevation from the MODIS product was estimated using the United States Geological Survey (USGS) HydroSHEDS. This dataset is available at https://hydrosheds.cr.usgs.gov. The HydroSHEDS elevation data are hydrologically conditioned data derived from the elevation of the Shuttle Radar Topography Mission (SRTM) at a 3 arc-second spatial resolution. A sequence of the automatic algorithm has been applied to SRTM data, including void filling and filtering, sink identification, stream burning and upscaling techniques.

Snow Cover Area Calculation
The snow cover area was estimated using original and improved MODIS snow products. The original products were also used with the <20% clouds threshold (as commonly used by previous studies [18,21]) for comparison with the improved product. In the products with 20% cloud thresholds, the original Aqua and Terra images with a percentage cloud cover more than 20% were discarded for the comparative assessment among five different products. The snow cover area was calculated using the total land area instead of the total cloud free pixels to highlight the difference in snow cover area estimation from different products.

Snowline Altitude Estimation
The SLA can be defined as the elevation above which snow exists after the melt season. However, at the end of the melt season, snow is often patchy in mountainous catchments due to wind erosion, avalanches and other topographic factors. Therefore, it is challenging to represent a snowline at a specific altitude. To overcome this difficulty, the WMO [34] recommended using an elevation band rather than an elevation line, which represents a zone of approximately 50% snow coverage. In this study, the snowline altitude was defined as the elevation band above which there are negligible snow-free pixels.
The snowline altitude was derived from the MODIS snow product and elevation data using the regional snowline elevation (RSLE) method developed by [14]. The snowline has been widely derived using the MODIS snow product in different parts of the world [35][36][37]. Comparisons of MODIS-derived snowlines from high-resolution DEMs and field-based studies have shown good agreement [5,38]. Additionally, the RSLE method by [14] has been widely used to investigate snow dynamics [37,[39][40][41][42]. This method estimates the elevation at which the number of snow-free pixels above the elevation line and the number of snowy pixels below the elevation line are at their minimum above and below the elevation line. The snowline altitude was identified by iterating elevation bands from the lowest to highest altitudes in the basin, evaluating the number of snow-covered pixels below and snow-free pixels above the elevation bands. We used elevation data from the USGS HydroSHEDS DEM with a spatial resolution of 90 m. This spatial resolution fulfils the requirement of this study because the MODIS snow product resolution is 500 m. Both products were masked by the area of interest. For each area of interest, the elevation raster was divided into different elevation bins with 10 m intervals, and the corresponding snow-covered pixels below the elevation band and snow-free pixels above the elevation band were evaluated. Finally, the 8-day snowline altitude was selected at which there was lowest snow-covered pixels below and lowest snow-free pixels above the elevation. Detailed information about the snowline algorithm can be obtained from [14].

Snow Cover and Snowline Trend Analysis
A set of nonparametric tests, i.e., Sen's Slope and the Mann-Kendall (MK) test, were used to evaluate the trends. A nonparametric test is commonly used for trend analysis in climate research [43]. Sen's Slope was used to detect the magnitude and direction while the MK test was used to test the statistical significance of the trend. The Sen's Slope method uses a linear model to estimate the trend slope, and the variance in the residuals should be constant in time, as calculated in [44]: where x j and x k are the data values at times j and k (j > k), respectively. In general, Sen's Slope is defined as the median of all linear slopes between x j and x k . The MK test was originally proposed by [45] and rephrased by [46]. The approach is widely used to test the null hypothesis (H 0 ) of no trend. Each value in a time series is compared with subsequent data values, and depending on whether the data value is higher or lower, the MK test statistic "S" is increased by 1 or decreased by 1. The statistic S can be computed as follows: where n is the number of data points, x i and x j are the data values in time series i and j (j > i), respectively, and sgn(x j − x i ) is the sign function, which is calculated as follows: Water 2020, 12, 2681 6 of 20 The variance is computed as follows: where n is the number of data points, m is the number of tied groups, and t i denotes the number of ties of extent i. A tied group is a set of sample data that have the same value. In cases where the sample size n > 10, the standard normal test statistic Z s is used and is calculated using the following equation.
Trends are tested at the α = 0.1 level of significance (90% confidence level). When |Z s | > Z 1−α/2 , the null hypothesis is rejected, and a significant trend exists in the time series. Z 1−α/2 is obtained from the standard normal distribution table. In this study, the null hypothesis of no trend is rejected if |Z s | > 1.64 at the 90% confidence level.

Results
We assessed SCA and SLA dynamics over the Karakoram region and its subbasins: Hunza, Shigar, and Shyok using open-source software R (https://cran.r-project.org/) and RStudio (https://rstudio.com/). Ggplot2 package (https://ggplot2.tidyverse.org) was used for data visualization. The SCA and SLA were examined between 2003 and 2018 for two dimensions: (i) interannual variability and (ii) intra-annual variability. A comparative basin-wide assessment was carried out using the original snow cover products, MYD10A2, MOD10A2, original products with <20% cloud cover, and the improved snow product MOYDGL06*.

Interannual Variability in SCA
The mean annual SCA from 2003 to 2018 in the Karakoram region and subbasins is shown in Figure 2. A comparative analysis of the mean annual SCA between the NSIDC original products, original products with 20% cloud threshold and our improved snow products reveals a significant underestimation in the original product. The results are mainly attributed to a large number of misclassified snow pixels in the original products. The cloud cover statistics are shown as percentages (>5 to >50) against the number of images in the original products in Table S1. The original 8-day Aqua and Terra snow products contain 15% and 9% clouds on average, respectively, in the Karakoram region. The Terra and Aqua 295 and 417 images contain more than 5% clouds on average, respectively. The <20% cloud threshold of the original snow reduces the annual variability to some extent. However, the cloud threshold products are still biased as compared to the improved product as depicted in Figure 2. The bias between improved and Aqua and Terra cloud threshold products for the Karakoram region is found to be −1.67% and 1.1%, respectively. This value is even higher (3.9% and 3.6%) in the original Aqua and Terra product. Therefore, we analyzed the snow dynamics in the improved cloud-free MOYDGL06* product. The average annual SCA of the Karakoram region between 2003 and 2018 is 67,606.5 km 2 , which is 65% of the total extent. The winter annual SCA fluctuates between 82% and 93%, whereas the summer SCA varies between 32% and 43%. In addition, the SCAs in Hunza and Shigar are 77%, whereas 64% of the Shyok Basin is covered by snow. The interannual SCA varies significantly during the study period. The annual average area covered by snow for Hunza ranges from 74% to 82%, with the lowest and highest snow occurred in 2007 and 2009, respectively. In contrast, Shigar received its maximum annual snow in 2009 (81%) and minimum in 2007 (72%). Shyok experienced its minimum annual snow in 2016 and its maximum in 2015. The detailed annual SCA statistics are presented in Table 2 and Table S2. The time series of annual maximum and minimum SCA are also shown in Figures S1 and S2. The original MODIS snow products underestimate the average annual SCA in the subbasins due to significant cloud cover variability (Figure 2), indicating high spatial variability across the region. The annual snow cover area is negatively biased in the original M*D10A2 snow products. On the contrary, the Aqua threshold product has random underand overestimation bias compared to the improved product. The Terra threshold product shows modest underestimation throughout the study period. In addition, the difference for the year 2012 is significantly large (>30% less snow cover area in MYD10A2 compared with the improved snow MOYDGL06* in the Shigar Basin). Original and threshold products show a comparable mean annual maximum SCA variability in both Aqua and Terra because of fewer clouds (<20%) during the snow accumulation season.   The MK test was applied to determine the annual SCA trend at a 90% confidence level in the Karakoram region and subbasins. The Sen's Slope estimator was used to determine the magnitude of the trends. Statistics analysis of the trends is given in Table 3 and Table S3. The annual minimum and mean snow exhibit contrasting trends in the original (MOD10A2 and MYD10A2) and improved products. This difference is reduced in the cloud threshold products except in a few cases with remaining bias. The magnitude of the annual maximum SCA trend using original and threshold products is equal across all spatial domains with a statistically significant negative trend in the Shigar and Shyok subbasins. The improved product shows that the annual minimum SCA is significantly (p = 0.1) decreasing in Shyok and insignificant in the other two basins. The annual average snow cover decreased by −0.13% in the entire Karakoram and −0.1% in the Hunza Basin. Similarly, the annual maximum SCA shows a decreasing trend, except in Shigar, where it is increasing at a rate of 0.22% per year. The most negative maximum snow trend is 0.23% per year in the Shyok Basin.

MOD10A2
Original The MK test was applied to determine the annual SCA trend at a 90% confidence level in the Karakoram region and subbasins. The Sen's Slope estimator was used to determine the magnitude of the trends. Statistics analysis of the trends is given in Table 3 and Table S3. The annual minimum and mean snow exhibit contrasting trends in the original (MOD10A2 and MYD10A2) and improved products. This difference is reduced in the cloud threshold products except in a few cases with remaining bias. The magnitude of the annual maximum SCA trend using original and threshold products is equal across all spatial domains with a statistically significant negative trend in the Shigar and Shyok subbasins. The improved product shows that the annual minimum SCA is significantly (p = 0.1) decreasing in Shyok and insignificant in the other two basins. The annual average snow cover decreased by −0.13% in the entire Karakoram and −0.1% in the Hunza Basin. Similarly, the annual maximum SCA shows a decreasing trend, except in Shigar, where it is increasing at a rate of 0.22% per year. The most negative maximum snow trend is 0.23% per year in the Shyok Basin. Table 3. Annual average snow cover trends (% per year) in the Karakoram region and subbasins. Italic bold numbers are statistically significant at the 90% confidence level.
We compared the annual average snow cover trend for two equally divided temporal windows and the entire study period, 2003-2010, 2011-2018, and 2003-2018, to understand the short-term and overall SCA variability. There are no statistically significant trends in any of the observed temporal windows in the observed period. Overall, the Karakoram region and Shyok Basin exhibit a positive trend in the annual average SCA between 2003-2010 and 2011-2018, whereas the long-term trend from 2003 to 2018 is negative. This result is due to the large interannual fluctuation in the SCA in the respective period. In contrast, the annual average SCA in Shigar indicates a negative trend in all three temporal windows. Similarly, the snow cover decreases during the first half and the overall study period in Hunza, although a slightly positive SCA is shown during 2011-2018. Despite contrasting characteristics in the short-term analysis, a decreasing regional signal was observed in the long-term SCA trend.

Intra-Annual SCA Variability
The seasonal variability in the SCA (%) derived from the original M*D10A2, <20% cloud threshold, and improved snow MOYDGL06* products in the Karakoram region is shown in Figure 3. The box and whisker plots represent quartiles of 8-day SCA variability at a monthly scale. The results indicate that the mean SCA is lowest in August in the original Terra (45%), cloud threshold Terra (44%) and improved product (41%) coverage. The lowest snow from Aqua is 48% in February, which is unexpected and attributed to persistent clouds (40%). However, this signal was not observed in the Aqua < 20% cloud threshold product. The mean SCA is the highest by 81% (median SCA = 82%) during late winter and early spring in the improved snow product. In contrast, the highest mean snow cover in the M*D10A2 products are comparatively low with 72% in Aqua and 78% in Terra during midspring. The monthly SCA variation shows conspicuous differences between different products. The maximum median snow cover extent was observed during March (81%) and April (79%) in the Terra and Aqua <20% cloud threshold products, respectively. The original MODIS product produces a bimodal distribution of intra-annual variation in the SCA, whereas the cloud threshold and improved products show a unimodal distribution. Generally, a bimodal distribution is common in the monsoon-dominated catchments and unimodal in the westerly-dominated regions [47]. However, the original product shows a bimodal distribution in the westerly-dominated Karakoram region. The bimodal pattern is attributed to the inherent uncertainty in the M*D10A2 products causing higher SCA variability. However, the cloud threshold products have a similar pattern to the improved data, with a noticeable difference during the winter months. Thus, the 8-day improved snow cover provides accurate snow cover information compared to the 8-day MODIS original and cloud threshold snow cover products.
analysis of the improved MODIS data shows a smooth increase in the SCA from September to February and a decrease from March to August. The duration from September to February is the snow accumulation period and March to August is the ablation period. For example, Figure 4 shows that 80% of the Shyok Basin was covered by snow between February and March, which reduces to <40% in August.  The monthly mean SCA variability is shown in Figures S3-S6 and the trends are presented in Figure 5 and Table S4 for all spatial domains. A distinct overestimation during summer and underestimation during winter were observed across all the regions in the original and cloud threshold products. The comparison of all the five snow products shows a significant difference. The magnitude of the trend ranges from 0.007% to −0.47% per year in the improved product and from −0.84% to 1.64% per year in the cloud threshold products. The trend reaches between −0.98% and 1.51% per year in the original products. This finding indicates a significant reduction and increment in the monthly SCA trend in the original and cloud threshold products. The M*D10A2 and cloud threshold products (except Aqua in the Shyok Basin) show a positive trend throughout January. In contrast, the improved snow (MOYDGL06*) product signifies a negative trend, although all of these trends are statistically insignificant. All products show a similar decreasing tendency for the months from March to August and October in the Karakoram and Shyok. Conversely, September and December show increasing SCA trends in the original Aqua and improved products in the The intra-annual (mean monthly) SCA variability using the original M*D10A2, 20% cloud threshold and MOYDGL06* snow products for individual basins is shown in Figure 4. The individual basins and the entire Karakoram region show a similar pattern on a monthly scale. However, a significant difference is observed in the winter snow between the M*D10A2 and improved MOYDGL06* products ( Figure 4). The original M*D10A2 products incorrectly show lower SCAs during winter. However, the improved product shows that the monthly mean SCA reaches its peak in February and is lowest in August. The <20% cloud threshold products indicate a maximum snow extent during April, which is somehow shifted from than the peak season. This may cause a dilemma in assessing the climate change impact using the product. The data-limited maximum snow shift may cause confusion between global climate change and data limitation. Thus, the original and cloud threshold M*D10A2 products are not recommended for snow climatology and hydrology. The analysis of the improved MODIS data shows a smooth increase in the SCA from September to February and a decrease from March to August. The duration from September to February is the snow accumulation period and March to August is the ablation period. For example, Figure 4 shows that 80% of the Shyok Basin was covered by snow between February and March, which reduces to <40% in August.
Water 2020, 12, x FOR PEER REVIEW 9 of 20 cause confusion between global climate change and data limitation. Thus, the original and cloud threshold M*D10A2 products are not recommended for snow climatology and hydrology. The analysis of the improved MODIS data shows a smooth increase in the SCA from September to February and a decrease from March to August. The duration from September to February is the snow accumulation period and March to August is the ablation period. For example, Figure 4 shows that 80% of the Shyok Basin was covered by snow between February and March, which reduces to <40% in August.  The monthly mean SCA variability is shown in Figures S3-S6 and the trends are presented in Figure 5 and Table S4 for all spatial domains. A distinct overestimation during summer and underestimation during winter were observed across all the regions in the original and cloud threshold products. The comparison of all the five snow products shows a significant difference. The magnitude of the trend ranges from 0.007% to −0.47% per year in the improved product and from −0.84% to 1.64% per year in the cloud threshold products. The trend reaches between −0.98% and 1.51% per year in the original products. This finding indicates a significant reduction and increment in the monthly SCA trend in the original and cloud threshold products. The M*D10A2 and cloud threshold products (except Aqua in the Shyok Basin) show a positive trend throughout January. In contrast, the improved snow (MOYDGL06*) product signifies a negative trend, although all of these The monthly mean SCA variability is shown in Figures S3-S6 and the trends are presented in Figure 5 and Table S4 for all spatial domains. A distinct overestimation during summer and underestimation during winter were observed across all the regions in the original and cloud threshold products. The comparison of all the five snow products shows a significant difference. The magnitude of the trend ranges from 0.007% to −0.47% per year in the improved product and from −0.84% to 1.64% per year in the cloud threshold products. The trend reaches between −0.98% and 1.51% per year in the original products. This finding indicates a significant reduction and increment in the monthly SCA trend in the original and cloud threshold products. The M*D10A2 and cloud threshold products (except Aqua in the Shyok Basin) show a positive trend throughout January. In contrast, the improved snow (MOYDGL06*) product signifies a negative trend, although all of these trends are statistically insignificant. All products show a similar decreasing tendency for the months from March to August and October in the Karakoram and Shyok. Conversely, September and December show increasing SCA trends in the original Aqua and improved products in the Karakoram and all the products except the Terra cloud threshold in Shyok. In general, the original products display growing SCAs in December and January throughout the study region. Though, the cloud threshold products also exhibit opposite trends during December, which stresses the need for improvement during winter months. The improved product found slightly increased snow cover trends during March and April in Hunza. The SCA in November increased in the Karakoram, Shigar and Shyok, and decreased in Hunza. This analysis shows a heterogeneous snow cover trend within the Karakoram region. Karakoram and all the products except the Terra cloud threshold in Shyok. In general, the original products display growing SCAs in December and January throughout the study region. Though, the cloud threshold products also exhibit opposite trends during December, which stresses the need for improvement during winter months. The improved product found slightly increased snow cover trends during March and April in Hunza. The SCA in November increased in the Karakoram, Shigar and Shyok, and decreased in Hunza. This analysis shows a heterogeneous snow cover trend within the Karakoram region.

Interannual Variability in Snowline Altitude
Interannual snowline altitude variability is estimated using all the snow products and HydroSHEDS Elevation data are presented as an annual boxplot in Figure 6. Compared to the original and cloud threshold snow, improved MOYDGL06* product provides enhanced snowline altitude information. The median SLAs throughout the study period (2003-2018) derived from the improved snow product is 4019 ± 736, 3507 ± 594, 3783 ± 696 and 4475 ± 812 m a.s.l. for the Karakoram, Hunza, Shigar, and Shyok, respectively. All of these values are higher than the SLA values derived from the MODIS original snow products. An underestimation was observed in the SLA derived from the Aqua cloud threshold product. On an annual scale, the median SLA from Terra cloud threshold product closely resembles the improved product. The original snow products (M*D10A2) and Aqua cloud threshold product were not successful in representing the annual median SLA dynamics represented by the thick black solid line in Figure 6 due to the underestimation of the SCA. The size of the box whiskers indicates that the first and third quartiles of the original products are lower than the improved product, indicating the underestimation of the SLA. Among all four areas, the 8-day snowline variability in Shyok is the highest of all three products. A regional signal of lower SLA fluctuation was observed in 2015 in all the products. The SLA values in the improved product are highest in Shyok and lowest in Hunza and Shigar. This noticeable difference between Shyok and Hunza-Shigar is due to the steep topography in Hunza-Shigar that enhances snow avalanches from

Interannual Variability in Snowline Altitude
Interannual snowline altitude variability is estimated using all the snow products and HydroSHEDS Elevation data are presented as an annual boxplot in Figure 6. Compared to the original and cloud threshold snow, improved MOYDGL06* product provides enhanced snowline altitude information. The median SLAs throughout the study period (2003-2018) derived from the improved snow product is 4019 ± 736, 3507 ± 594, 3783 ± 696 and 4475 ± 812 m a.s.l. for the Karakoram, Hunza, Shigar, and Shyok, respectively. All of these values are higher than the SLA values derived from the MODIS original snow products. An underestimation was observed in the SLA derived from the Aqua cloud threshold product. On an annual scale, the median SLA from Terra cloud threshold product closely resembles the improved product. The original snow products (M*D10A2) and Aqua cloud threshold product were not successful in representing the annual median SLA dynamics represented by the thick black solid line in Figure 6 due to the underestimation of the SCA. The size of the box whiskers indicates that the first and third quartiles of the original products are lower than the improved product, indicating the underestimation of the SLA. Among all four areas, the 8-day snowline variability in Shyok is the highest of all three products. A regional signal of lower SLA fluctuation was observed in 2015 in all the products. The SLA values in the improved product are highest in Shyok and lowest in Hunza and Shigar. This noticeable difference between Shyok and Hunza-Shigar is due to the steep topography in Hunza-Shigar that enhances snow avalanches from steep snow accumulation areas to the bottom of a hill slope. Approximately 56% of the Hunza and 47% of the area in the Shigar Basin has a slope greater than 25 degrees, which favours snow avalanches and redistribution. In contrast, the SLA of the Shyok Basin in the eastern Karakoram is the highest, where the intensity of the westerlies is weakest (compared with the other two basins). The annual SLA summary for the Karakoram region and its three basins is shown in Table 4 and Table S5. The annual average SLAs for the Karakoram in MYD10A2, 20% MYD10A2 cloud threshold, MOYDGL06* MOD10A2 and 20% MOD10A2 cloud threshold are 3652, 3807, 3991, 3817 and 3968 m a.s.l., respectively. A similar pattern of underestimation is demonstrated by the mean annual maximum SLAs. Original and cloud threshold products provide nearly equal maximum annual SLAs throughout the study area. The annual minimum SLA shows lower values in the original and cloud products compared to the improved product. This highlights the underestimation of the SLA in all the original products even at the annual scale. The three basins in the Karakoram exhibit analogous behaviour. The annual SLA fluctuates in the improved snow product compared to the original snow products. The annual snowline exhibits immense spatial variability in the Karakoram region, with mean annual maximum values of 4655, 5259, and 5145 m a.s.l. in MYD10A2, MOYDGL06*, and MOD10A2, respectively. The mean maximum annual SLA in the Shyok is the highest among the basins and the lowest SLA occurs in the Hunza. Similar patterns were observed in the mean annual SLA estimates. The snowline elevation shows a gradually decreasing trend from the eastern to western Karakoram. Specifically, the eastern Shyok basin has a higher snowline altitude than the central and western regions.  The annual SLA summary for the Karakoram region and its three basins is shown in Table 4 and Table S5. The annual average SLAs for the Karakoram in MYD10A2, 20% MYD10A2 cloud threshold, MOYDGL06* MOD10A2 and 20% MOD10A2 cloud threshold are 3652, 3807, 3991, 3817 and 3968 m a.s.l., respectively. A similar pattern of underestimation is demonstrated by the mean annual maximum SLAs. Original and cloud threshold products provide nearly equal maximum annual SLAs throughout the study area. The annual minimum SLA shows lower values in the original and cloud products compared to the improved product. This highlights the underestimation of the SLA in all the original products even at the annual scale. The three basins in the Karakoram exhibit analogous behaviour. The annual SLA fluctuates in the improved snow product compared to the original snow products. The annual snowline exhibits immense spatial variability in the Karakoram region, with mean annual maximum values of 4655, 5259, and 5145 m a.s.l. in MYD10A2, MOYDGL06*, and MOD10A2, respectively. The mean maximum annual SLA in the Shyok is the highest among the basins and the lowest SLA occurs in the Hunza. Similar patterns were observed in the mean annual SLA estimates. The snowline elevation shows a gradually decreasing trend from the eastern to western Karakoram. Specifically, the eastern Shyok basin has a higher snowline altitude than the central and western regions. The results of the annual SLA trend analysis for the Karakoram region and its subbasins between 2003 and 2018 are presented in Table 5 and Table S6. The annual SLA shows a similar pattern in the improved and original Terra products throughout the region, whereas the Aqua estimates an opposite pattern in some cases in the annual minimum and average series. This signal is more prominent in Aqua cloud threshold product. For example, negative trends in the annual minimum SLA were identified using the Aqua products in the Karakoram and Shyok, whereas improved products show a constant and positive trend. Likewise, the original Aqua shows a negative trend in the annual maximum SLA for Hunza and Shigar, whereas the opposite trends are noticed in the Terra and improved products. The annual minimum SLA trend in the Karakoram region is strongest using the improved snow data, with a rate of 14.33 m a.s.l. per year. The annual average SLA in the improved product exhibits an increasing trend with a statistically significant trend of 8.03 m a.s.l. per year in the Karakoram region. The snowline is often taken as a proxy of the ELA and used for understanding the glacier mass balance [22]. The maximum SLAs are extracted using the improved snow product for each year ( Figure S7) and used to estimate the trend (Table S6). The increasing trend of the maximum annual SLA during 2003-2018 is observed in all the subbasins. The rate of increase in the annual maximum snowline is highest (13 m a.s.l. per year) and statistically significant at the 90% confidence level on a regional scale and lowest (13 m a.s.l. per year) in the Shigar Basin. The annual average SLA trends in the three defined temporal windows are estimated similarly to the SCA (Table S7) [27,48]. Our findings do not support a positive glacier mass balance in the Karakoram region and its subbasins throughout the study period.

Intra-Annual Variability in Snowline Altitude
The spatiotemporal variations of the SLA are assessed in the Karakoram region and its subbasins. The seasonal variability in the 8-day SLA is presented in Figure 7. The box whiskers represent quartiles of the monthly variability in the snowline estimated during 2003-2018. A clear underestimation of the SLA throughout the study period was observed by the original and the cloud threshold products across all the study sites as indicated by the thick line drawn inside each box in Figure 7. All products show a significant fluctuation in SLA during autumn and winter months. The snowlines are lowest between January and February throughout the region using all products. The lowest monthly median SLA using the improved MOYDGL06* product is 3149, 2957, 2623, and 3385 m a.s.l. in the Karakoram, Hunza, Shigar, and Shyok, respectively. The time of the minimum snowline coincides with the observed maximum SCA. In contrast, the minimum median SLA estimates in the original snow products are dissimilar to the SCA statistics. The highest monthly median SLA was observed in August using all products throughout the region. This observation indicates that the SLA estimates represent the snow dynamics better than the SCA in the original products at the seasonal scale. Original and cloud threshold products provide nearly identical monthly median values indicating the robustness of the snowline algorithm to deal with the cloudy image to some extent. However, we recommend the use of the improved product to understand snow variability more confidently. The median monthly SLA in Shyok using the improved product is the highest (5545 m a.s.l.). For the Karakoram region, the highest monthly median SLA is 4999 m a.s.l. The snowline varies significantly in autumn and winter due to snowfall and seasonal snow onset. The snowline starts to shift downward after August due to the combined effect of a lower temperature and fewer precipitation events and reaches its minimum value in January and February.
Water 2020, 12, x FOR PEER REVIEW 13 of 20 indicating the robustness of the snowline algorithm to deal with the cloudy image to some extent. However, we recommend the use of the improved product to understand snow variability more confidently. The median monthly SLA in Shyok using the improved product is the highest (5545 m a.s.l.). For the Karakoram region, the highest monthly median SLA is 4999 m a.s.l. The snowline varies significantly in autumn and winter due to snowfall and seasonal snow onset. The snowline starts to shift downward after August due to the combined effect of a lower temperature and fewer precipitation events and reaches its minimum value in January and February. The annual SLA trends suggest that finer temporal resolution data capture the snow dynamics. The monthly trend analysis is presented in Figure 8 and Table S8. A corresponding time series of monthly SLA are plotted in Figures S8-S11. An underestimation of the SLA (lower SLA values) is observed in both original and cloud threshold products. There is a large discrepancy in the trend calculated from the original and cloud threshold products. In particular, the Aqua snow product mostly exhibits a contrasting trend compared with the other products. The improved product shows an increasing tendency in the entire monthly SLA series in the Karakoram region, although the trends are not statistically significant, owing to the relatively short period of analysis. Nonetheless, the trends are stronger in the winter months. Similarly, most of the monthly SLA trends in the individual basin are also negative except March and April for Hunza, November for Shigar and September, November and December for Shyok. The spatial and temporal monthly SLA trend indicates an inconsistent variation. The annual SLA trends suggest that finer temporal resolution data capture the snow dynamics. The monthly trend analysis is presented in Figure 8 and Table S8. A corresponding time series of monthly SLA are plotted in Figures S8-S11. An underestimation of the SLA (lower SLA values) is observed in both original and cloud threshold products. There is a large discrepancy in the trend calculated from the original and cloud threshold products. In particular, the Aqua snow product mostly exhibits a contrasting trend compared with the other products. The improved product shows an increasing tendency in the entire monthly SLA series in the Karakoram region, although the trends are not statistically significant, owing to the relatively short period of analysis. Nonetheless, the trends are stronger in the winter months. Similarly, most of the monthly SLA trends in the individual basin are also negative except March and April for Hunza, November for Shigar and September, November and December for Shyok. The spatial and temporal monthly SLA trend indicates an inconsistent variation.
Water 2020, 12, x FOR PEER REVIEW 14 of 20 Figure 8. Bar plot illustrating monthly SLA trend. Bars with plus signs represent a statistically significant trend at the 90% confidence level.

Discussion
We used an improved MODIS snow product (MOYDGL06*) and compared the results with the original M*D10A2, and cloud threshold products to characterize spatiotemporal SCA and SLA variations in the Karakoram region and its major subbasin. Most of the SCA and SLA trends are statistically insignificant in the relatively short observation period which limits the climate-driven physical and statistical explanation of snow dynamics which are often important aspects of climate research. We admit that the derived trend should not be understood as the signal of long term climate trends. However, the analysis of this study does provide the pattern of snow cover area change and highlights the importance of improved products over original and cloud threshold snow products. The annual SCA indicates an insignificant decreasing trend in the Karakoram region and its basins except for Shigar, where the annual maximum SCA increases by 0.22% per year. We also analyzed the monthly mean SCA, indicating a slightly positive trend, especially in the months of November and December throughout the region. However, the trends decrease during spring, particularly in March and April and during summer months. The overall spatiotemporal SCA analysis shows significant heterogeneity in the Karakoram region. In addition, the annual and monthly SLAs show higher snowline altitudes in Shyok compared with Shigar and Hunza. The higher SLA of Shyok is mainly due to weak westerly disturbances compared with the Hunza and Shigar basins. The annual average and maximum SLA of the entire Karakoram region show a statistically significant positive trend in contrast to the monthly insignificant SLA trend.
The quantitative comparison of SCAs in the original products (MYD10A2 and MOD10A2), cloud threshold and improved products (MOYDGL06*) show a large underestimation in the original products on the annual scale ( Figure 2). A noticeable difference is observed in the original products with a <20% cloud threshold. This difference is mainly due to the presence of clouds at higher altitudes in the original products. The underestimation is high in the Aqua snow product due to its afternoon revisit time when clouds become dense. In contrast, the seasonal SCA comprises both overestimation and underestimation. The overestimation is mostly found during the summer and early autumn, and the underestimation is found during winter and spring. The summer overestimation is mainly due to the misclassification of clouds as snow at lower altitudes and snow as clouds during winter at higher altitudes. Clouds are a primary constraint in optical remote sensing

Discussion
We used an improved MODIS snow product (MOYDGL06*) and compared the results with the original M*D10A2, and cloud threshold products to characterize spatiotemporal SCA and SLA variations in the Karakoram region and its major subbasin. Most of the SCA and SLA trends are statistically insignificant in the relatively short observation period which limits the climate-driven physical and statistical explanation of snow dynamics which are often important aspects of climate research. We admit that the derived trend should not be understood as the signal of long term climate trends. However, the analysis of this study does provide the pattern of snow cover area change and highlights the importance of improved products over original and cloud threshold snow products. The annual SCA indicates an insignificant decreasing trend in the Karakoram region and its basins except for Shigar, where the annual maximum SCA increases by 0.22% per year. We also analyzed the monthly mean SCA, indicating a slightly positive trend, especially in the months of November and December throughout the region. However, the trends decrease during spring, particularly in March and April and during summer months. The overall spatiotemporal SCA analysis shows significant heterogeneity in the Karakoram region. In addition, the annual and monthly SLAs show higher snowline altitudes in Shyok compared with Shigar and Hunza. The higher SLA of Shyok is mainly due to weak westerly disturbances compared with the Hunza and Shigar basins. The annual average and maximum SLA of the entire Karakoram region show a statistically significant positive trend in contrast to the monthly insignificant SLA trend.
The quantitative comparison of SCAs in the original products (MYD10A2 and MOD10A2), cloud threshold and improved products (MOYDGL06*) show a large underestimation in the original products on the annual scale ( Figure 2). A noticeable difference is observed in the original products with a <20% cloud threshold. This difference is mainly due to the presence of clouds at higher altitudes in the original products. The underestimation is high in the Aqua snow product due to its afternoon revisit time when clouds become dense. In contrast, the seasonal SCA comprises both overestimation and underestimation. The overestimation is mostly found during the summer and early autumn, and the underestimation is found during winter and spring. The summer overestimation is mainly due to the misclassification of clouds as snow at lower altitudes and snow as clouds during winter at higher altitudes. Clouds are a primary constraint in optical remote sensing when estimating snow cover, particularly in mountainous catchments. Additionally, a larger SZA [49] and a wide swath of MODIS cause significant overestimation of snow [50,51]. The presence of clouds in the original products could even alter the trend pattern (e.g., min and mean annual SCA trends, as shown in Table 3). At the monthly scale, the original products show an inverse trend compared with the improved snow product, mostly during winter months.
We compared the annual minimum SCA calculated in this study with the Randolph Glacier Inventory version 6.0 (RGI6.0) glacier boundaries and found that the annual minimum snow cover estimates are higher than the total glacier cover. The snow cover is usually more than the glacier cover, as there are firn and perennial snow at high altitudes. Although the old glacier ice is not necessarily mapped by MODIS, the uncaptured ice is less than the firn and perennial snow at high altitudes. Hasson et al. [52] removed clouds in the snow product and found a decreasing trend in Hunza, Shigar and Shyok's annual average snow cover analysis. Despite the increasing trend in winter precipitation in the region [18,29], the annual SCA is decreasing according to the results of this study. This finding could be explained by enhanced melting from the observed warming during winter [53,54]. The Khunjerab climate station in Hunza, which is located at 4710 m a.s.l., also shows an increasing trend in both annual and seasonal temperatures [21,55]. On the other hand, earlier studies have shown an increasing trend in the SCA [18,24,33]. The difference could be associated to the method used for removing clouds and the observation period. We used different filters to remove clouds and convert cloudy pixels to snow or not snow, whereas many studies adopted the threshold method (cloud <10%, 20%) to avoid cloudy images. The 20% threshold method provides a similar trend as the improved product in annual SCA series in most cases but the magnitude of the trend differs. The threshold method is susceptible to missing important snow information because the original products have significant uncertainity. Our analysis shows that 122 Terra images and 189 Aqua images (out of 735) have more than 20% clouds persisting on eight consecutive days in the Karakoram region during our study period. The cloudy images are even more abundant in the individual basins (Table S1). A comparison of the annual average SCA and SLA trends between the two different temporal windows shows heterogeneous snow characteristics (Table S7) The improved product indicates a slight rise in the annual SLA across the Karakoram region from 2003 to 2018. The mean annual maximum SLA is 5259, 4636, 4637, and 5603 m a.s.l. in the Karakoram, Hunza, Shigar, and Shyok, respectively. These values are close to the range identified by an earlier study [23]. However, our improved MODIS-based mean annual maximum SLA for Hunza is slightly lower than those of [22] (average ELA = 5176 m a.s.l.), estimated from Landsat. The slight difference is due to the spatial resolution and methodology used for snowline estimation. The annual maximum SLA shows a rising pattern over the entire Karakoram region. Huss and Hock [56] also projected an upward trend in the regionally average glacier ELAs between 100 and 500 m during 2010 to 2100 for southwest Asia. The magnitude of the SLA trend is high in Hunza and low in Shigar and Shyok. These findings suggest significant spatial variability in the SLA. In addition, the ELA in Hunza and Shigar are comparatively lower, mainly due to the influence of strong westerly disturbances and steep terrain. Our results agree with earlier studies [22,57], who also identified a strong east-west gradient in precipitation patterns and the ELA in this region. No trend is statistically significant in the individual basin, whereas the Karakoram region shows a significant positive trend. Our snow cover estimates from 2003 to 2018 in the Karakoram region are in line with the recent slight mass loss of glaciers [27,48,58,59]. The short-term estimates during 2003-2010 and 2011-2018 mostly indicate an increase in the SCA. Our areal expansion of snow cover may not necessarily represent glacier mass changes because thickness and surface elevation change is the most reliable estimate of glacier mass balance [60,61]. Although the annual snowline is useful to approximate the ELA and characterize the glacier mass balance, for an improved understanding of the accumulation and melt processes of snow, more field-based observations of the variables, e.g., snow depth and snow water equivalent, should be carried out.
Our analysis shows that the original NSIDC snow cover product (M*D10A2) contains significant uncertainty, which is mainly due to the cloud coverage and overestimation in the original products because of the larger SZA. The twenty percent threshold method also shows some uncertainty at all monthly, seasonal and annual scales. The summer overestimation is due to the larger SZA and the winter underestimation is the result of the misclassification of snow as cloud. The underestimation of snow during winter will result in a reduced snowmelt during the spring season when it is the major discharge source. Similarly, this underestimation may result in undercatching the enhanced snowmelt. Additionally, summer snow overestimation leads to more river runoff and meltwater generation. The improved product reduces these uncertainties and will eventually improve the snowmelt estimation when used in glaciohydrological modeling. The improved product uses spatial and temporal filters to recover the snow pixels from the adjacent before and after images. This criterion may lead to some uncertainty in snowmelt simulations from the improved product, particularly in the peak ablation season. However, the uncertainty in glaciohydrological modeling is expected to be negligible because the snow underestimation (recovered in our improved product) is insignificant on average by 3.66% in the region [12]. In addition, only part of the 3.66% snow recovered in lower altitudes may cause some uncertainty as the high altitude areas are usually snow-covered until the late ablation season. Therefore, we suspect a significant uncertainty in water resource planning based on the original snow data. Although the MOYDGL06* product has significantly improved the snow dynamics in the region, for better understanding, we suggest observing the snow depth, snowmelt, and snow water equivalent. It is important to mention that the SCA does not provide information about the water content in the snowpack. Thus, in addition to SCA and SLA assessments, we also recommend the assessment of the snow water equivalent to investigate the change in the water content of snowpacks and understand the change in water availability in the region. The MOYDGL06* product may be useful to improve the quality of the existing coarser snow water equivalent products by reducing overestimation in the extent and underestimation due to cloud cover.

Conclusions
This study assesses the variability in SCA and snowline altitude in the Karakoram region and its major river basins. In comparison to the previous snow studies, we examined the monthly and annual snow cover patterns using the improved MODIS snow product MOYDGL06* in the Karakoram region. Our results show that there is a significant difference in SCA and SLA statistics derived from original, <20% cloud threshold and improved products. The annual SCA in the improved product exhibits an overall decreasing pattern except for the annual minimum SCA series in Shigar, whereas the original product shows mixed trends. However, none of the trends are statistically significant in the improved product except for the annual maximum SCA of 0.22% per year in Shyok. The negative trend of the snow cover area is also verified by the rising pattern of the snowline altitude. Thus, it is found that the simultaneous assessment of the SCA and SLA provides a more robust estimate of the snow cover change. The cloud threshold Terra product successfully captured the annual SCA trend in most cases. The improved product reveals that the highest snow coverage was during February-March and the lowest snow coverage was in August. In contrast, the original and cloud threshold products overestimated the SCA during summer and underestimated the SCA during winter, which may constrain seasonal snow cover changes and associated glaciohydrological studies. The annual SLA analysis shows a decreasing trend throughout all of the study regions. In particular, the annual maximum SLA for the entire Karakoram region shows a statistically significant upward trend of 13 m a.s.l. per year, indicating a negative glacier balance. The seasonal SLA exhibits a minimum median value in January and February and a maximum median in August. The evolution of the SLA in the original product, cloud threshold and improved products is more similar than the SCA analysis.
From the overall assessment, we conclude that the SCA and SLA in the Karakoram are highly heterogeneous in space and time. A comparative analysis reveals a significant uncertainty in the original snow products. Additionally, the 20% cloud threshold product remained uncertain at the monthly, seasonal and annual scale. The SCA trend during 2003-2018 does not exhibit any significant variations; however, a noticeable difference is found in the subdivided temporal windows of 2003-2010 and 2011-2018. Despite the global spatial coverage of the MODIS, its relatively short observation period limits the comprehensive assessment of snow in response to climate change. We believe that our SCA and SLA estimates and analysis using the improved product (MOYDGL06*) are more robust than using the original and the cloud threshold snow products. We recommend improving the NSIDC MODIS snow cover product to remove clouds and overestimation before using it for any snow hydrology and climatology applications.  Table S1: Cloud cover in the original M*D10A2 products. Number in the table represents the images. A total of 735 images are analyzed, Table S2: Mean and standard deviation of annual minimum and maximum SCA (%) of the Karakoram region and subbasins using all the snow products, Table S3: Annual minimum and maximum snow cover trends (% per year) in the Karakoram region and subbasins. Italic bold numbers are statistically significant at the 90% confidence level, Table S4: Trend analysis of monthly mean SCA (% per year) from 2003 to 2018. Numbers in italicized bold are statistically significant at 90% confidence level, Table S5: Mean and standard deviation of annual minimum and maximum SLA (m a.s.l.) for the Karakoram region Table S6: Annual average SLA trend analysis (m a.s.l. per year) using all the products. Italic bold numbers are statistically significant at the 90% confidence level, Table S7: Trend analysis of annual average SLA (m a.s.l. per year) for the different temporal windows, Table S8