**3. Results**

#### *3.1. Seasonal Snow Cover from In Situ Measurements*

In general, seasonal snow cover lasts from the beginning of October until the end of March with strong variability during the transition period. This transition period refers to late autumn and early spring with occasional intermittent snowfalls that develop a non-continuous snow cover that is relatively thin and temporary. Observed snow depth derived with the combination of Ta < 0 ◦C and albedo ≥ 0.25 showed that the total SCD ranged from 145 to 157 days with a mean depth of 12–25 cm/year from 2012–2017. The intermittent and persistent snow cover days were distinguished for each year (Table 1). The duration of persistent snow cover was relatively stable while the intermittent snow cover days showed strong variability ranging from 6 to 21 days. The maximum snow depths

were observed in the second half of March before the main melt started. There was no clear relationship between the variability of snow depth and the duration of snow cover.

**Table 1.** Intermittent and persistent snow cover days (SCD) during the observation period. 1 and 2 denote the persistent SCD derived from the combination of albedo and air temperature and below −1 ◦C of daily mean surface temperature, respectively.


As snow cover influences the surface energy balance significantly [46], the temperature at the underlying surface can serve as a proxy of absence and presence of seasonal persistent snow cover. Therefore, we validated persistent SCD using Ts at several locations, where snow cover was underestimated from MODIS sensors, and which is demonstrated in Figure 4 by plotting snow depth together with surface temperature measurements for each year at the hydro-climatic station.

**Figure 4.** An example of the discrimination of intermittent and persistent snow cover, determined from the combination of albedo and air temperature measurements. The red and blue lines indicate the observed daily mean surface temperature and the threshold of −1 ◦C, respectively. The black line describes the typical snow accumulation during the snow season.

For the observation period at the hydro-climatic station, the distinguished persistent SCD showed a moderate and positive correlation (r = 0.64, p < 0.05) with the persistent SCD determined from the surface temperature record. These results sugges<sup>t</sup> further that persistent SCD in other parts of the catchment could be verified using the surface temperature iButton measurements. Snow field measurements were conducted at each iButton location, and south-facing slopes were snow free while both valley bottom and north-facing slopes had similar SD, indicating no changes with increasing altitude (see details of the snow field measurements in Table S2, Supplementary Materials).

#### *3.2. Temporal Distribution of Snow*

### 3.2.1. Cloud Reduction

To reduce cloud cover in the daily Terra and the combined Aqua and Terra images, the temporal gap-filling (GF) and conditional adjustment (CA) were applied to both datasets. Figure 5 shows the results as the mean cloud coverage. The most obvious decrease appeared after the temporal GF filtering [47]. The relative percentage of cloud cover for the entire study area decreased up to 100% from the initial MODIS data to the final CA filtering. The combination of the two initial data yielded only slight improvements in reducing cloud coverage because the time difference between the two overpasses was about 3 h; hence general atmospheric conditions remained similar. The difference between Terra and the combined version after the final CA was not significant.

**Figure 5.** Changes in percent of cloud coverage after a series of cloud reduction steps applied to daily MODIS snow products Aqua and Terra; Combined—The combination of both satellites; Terra\_GF—Gap-filling applied to only Terra; Comb\_GF—Gap-filling applied to the combination; Terra\_CA—Conditional adjustment applied to Terra\_GF; Comb\_CA—Conditional adjustment applied to Comb\_GF.

#### 3.2.2. SCD and Trend Analysis

Altitude, latitude, and solar radiation are known factors that govern SCD as well as their spatial distribution [48,49]. Due to the study area's relatively small size, stretched from east to west, the derived SCD from the improved MODIS data showed an insignificant relationship with potential solar radiation (p = 0.33); therefore, we ignored the influence of latitude and topography, and considered only the altitude dependency for deriving the SCD rate (days/m). The mean annual SCD ranged from 124–226 days from the lowland plain area to the mountain peak, showing clear underestimation over the forested and heterogeneous topography area (Figure 6). The mean increasing SCD rate was +6 days/100 m with mean coefficient of determination r2 = 0.92 (ranging from 0.71–0.96) with the probability value p < 0.001, and it agrees well with the result (+5.9 days/100 m) found in other Central Asian mountains [7].

**Figure 6.** Mean annual snow cover duration (SCD) for 2000–2017 from the improved MODIS-Terra and the derived mean SCD rate (inset image). Green and blue dots denote the locations of hydro-climatic station and selected random points, respectively. Areas with heterogeneous topography and forest cover mostly in the middle part of the catchment demonstrate significant underestimation, such as only 13 days, as shown in the scale.

The trend analyses of MODIS derived snow metrics (So, Sm, and SCD) and long-term air temperature of snow season (from October to April) are shown in Figure 7. The So and Sm dates exhibited large inter-annual variability with insignificant slightly increasing delays in both So (p = 0.53) and Sm (p = 0.77) dates, indicating a slight shift in the snow season. To implement trend analysis for the total SCD, we separated the study area into two parts (mountain and plain area) based on an elevation threshold (2300 m.a.s.l) of the tree line [25], excluding the area that is always underestimated from the MODIS data because of the forest canopies and topographic heterogeneity. Overall, the SCD over both the mountains and the plains showed insignificant decreasing trends (p = 0.96 for plains and p = 0.41 for mountains). The mean change rates were −0.5 days/year for mountains and −0.04 days/year for the plains, respectively. The mean SCD for mountain peaks dropped down to its minimum in winter 2014/2015 and increased to its maximum in winter 2012/2013. On the contrary, SCD in the plains showed high variability between 2000 and 2008, and then it stabilized. The hydrologic year 2001/2002 had the shortest SCD in the plains due to late snow cover onset and early snow melt, and is in good agreemen<sup>t</sup> with the other findings in Central Asia [50]. In addition, there was no decline in snow cover over our study period, which is also in an agreemen<sup>t</sup> with other studies in the region, e.g., [1]. Nevertheless, the time series of 17 years may be relatively short to produce a trend analysis of these snow metrics, and a longer time-series of data may be necessary for a more comprehensive analysis.

**Figure 7.** Trend analysis for the mean snow onset, snow melt, and SCD for 2001–2017, as well as monthly mean air temperature during snow season (ONDJFMA). Mountain peaks and plain areas were separated for SCD analysis. The temperature trend was analyzed using daily mean values from the nearest long-term observational site at the Chinggis Khan International Airport. There was no significant trend in each snow metric, while a slightly decreasing trend (p < 0.05) in air temperature was detected during the snow season, including intermittent and persistent snow periods.

Interestingly, the monthly mean air temperature during the snow seasons, including intermittent and persistent snow periods, showed a decreasing trend over the last three decades. This decrease was statistically significant (p < 0.05), but only for the mid-season months (DJF) (Table 2). Unlike in many regions in the world, a decrease or no increase in winter air temperature has been observed since the 1990s, replacing the general increasing trend since the 1960s over Central Asia under the influence of weakening SH [9,10,49]. The mean temperature of the snow season also showed an insignificant negative relationship (r2 = 0.19, p = 0.07) with the SCD on the mountain top, indicating the decreasing temperature may have played an important role in preventing SCD decrease.


**Table 2.** The significance of trend analysis in monthly mean air temperature over the last 31 years at the Chinggis Khan International Airport (60 km from the study area).

 Significant, when p < 0.05

\*

#### *3.3. Spatial Distribution of Snow*

For analyzing the spatial distribution of snow in the study area, we exploited the high spatial resolutions of Landsat-7, Landsat-8, and Sentinel-2A images to create a sequence of snow cover maps during the snow season 2016/2017 (Figure S1, Supplementary Materials). The first snow covered area was captured on 11 October 2016 over mountain peaks by Landsat-7, and then decreased by 16 October after 5 days as captured by Sentinel-2 during the autumn transition period. The percent of SCA in the study area increased following several snow events and eventually reached 99.6% by 21 November, including the hydro-climatic station site (Figure 8). The persistency of seasonal snow cover for the Sugnugur catchment began after the heavy snow events from 6–7 November 2016.

**Figure 8.** Precipitation events, distinguished as rain or snow, observed at the hydro-climatic station using the combination of albedo and air temperature measurements, and the inter-annual variability of SCA detected by different satellites over the Sugnugur catchment for winter 2016/2017.

We were able to use 28 images (70% of the total data) to identify the spatial distribution of seasonal persistent snow cover in 94.2% of the total area. The intermittent and snow free areas were 5.3% and 0.5% of the catchment area during the winter 2016/2017, respectively. The different spatial resolution of the two sensors may affect their combined use of time-series analysis [51] and result in minor uncertainties, mostly in the intermittent SCA. The intermittent snow areas appeared mostly on the steep south-facing slopes. Following several snow events, the south-facing slopes were covered by snow, but shortly became snow free due to higher sublimation rates or blowing-wind [52]. However, the intermittent SCD on the south-facing slopes still remain unclear because of the coarse spatial resolution of MODIS and non-daily coverage of the combined Landsat and Sentinel overpasses. The absence of snow cover on south-facing slopes was captured by our time-lapse camera, and was also visually inspected during the snow field campaigns (Figure S2 and Table S2, Supplementary Materials).

#### *3.4. Mapping Snow Cover*

Finally, we created a snow duration map (Figure 9) with 30 m resolution based on the findings that included the combination of both ground measurements and various optical remote sensing images. The map was produced only for seasonal persistent SCA because SCD for intermittent SCA had higher uncertainty. The simple correction, which was applied to Landsat and Sentinel images

using the supervised land-cover classification map, improved the spatial accuracy of SCA that was not detectable or often underestimated by MODIS satellites.

**Figure 9.** The total snow cover duration (SCD) days in Sugnugur catchment calculated only for persistent snow cover area. The intermittent and snow free areas were excluded due to their unknown snow accumulation periods.

The mean altitude-dependent SCD rate +6 days/100 m was consistent with that found in the same region [7]. The south-facing slopes were mostly snow free or with intermittent snow cover during winter due to high sublimation rate, wind blowing, and slope steepness. The mean annual SCD ranged from 124–226 days from the lowland plain area to the mountain peak (elevation difference of ~1700 m and horizontal distance of ~35 km), or showing high variability for such a small scale. Therefore, snow-melt water from the higher elevations with extended snow season is extremely important for balancing the Sugnugur River discharge during spring and autumn dry seasons.

To match our temporal record of snow depth observations at the hydro-climatic station, we compared the estimated SCD with the total observed SCD for the winters of 2012/2013 to 2016/2017. There was a clear underestimation of satellite derived SCD and resulted in an error of 12–13 days/year (Figure S3, Supplementary Materials). Previous studies have described MODIS snow products in clear-sky conditions to return high accuracies of >93% for Terra and >90% for Aqua [28,53]. Assuming these accuracies, our underestimation equates to 9% and is within the MODIS accuracy range. Nevertheless, errors are likely propagated through the relatively weak snow detection rate of MODIS sensors during the transition periods [13,54], the position of the hydro-climatic station relative to the receipt of incoming solar radiation, as well as the accuracy of annual altitude-dependent SCD rate.

The absence of snow cover on the south-facing slopes was captured by the automated time-lapse camera system, as well as was visually observed during fieldwork, which took place in the beginning of March 2017 and 2018, or the accumulated persistent snow cover period. The photographs from the camera showed a light snow event on 5 March 2017 together with the previous and following days, indicating the south-facing slopes were snow-free before the event (Figure S2, Supplementary Materials). The snow field measurements showed similar snow depths on both valley bottom and north-facing slope without any altitude-dependent increase. Therefore, we speculate that snow depth was homogeneous for the seasonal persistent SCA over the study area.

The days with below −1 ◦C of Ts at the hydro-climatic station showed a positive relationship with the observed persistent SCD days. This suggests the possibility of validating the estimated SCD against the days with Ts < −1 ◦C in the upper part of the catchment. For this analysis, surface temperature measurements at the valley bottom and north-facing slopes in different land-cover types were used

because the south-facing slopes were snow-free, and thus were excluded from the final SCD map. The estimated SCD determined from Ts revealed good agreemen<sup>t</sup> with an r2 of 0.85 (p < 0.001) and a mean bias of ±1 day compared to the hydro climatic station site. This would mean that the mean altitude-dependent SCD rate is consistent for persistent snow cover areas where MODIS snow data are contaminated by the topographic heterogeneity, forest, and shrubs.
