*2.2. Methods*

The main idea of the methodology applied in this study was to downscale the coarse resolution MODIS data based on a combination of other freely available optical remote sensing images, as well as ground measurements. An altitude-dependent SCD rate derived from MODIS data and DEM coupled with snow cover area (SCA) derived from the combination of Landsat and Sentinel images were used to produce annual SCD maps with 30 m spatial resolution from 2000–2017. We used the daily temporal resolution of MODIS Aqua and Terra satellites for deriving the temporal distribution of snow over the last 17 years, particularly the altitude dependent SCD rate, snow onset (So), and melt (S m) dates. The Landsat and Sentinel images, with high spatial resolutions, were used for constructing the spatial distribution of snow in the study area. The abbreviations used in this study are shown in Table A1, Appendix A, and the flowchart of the procedure can be found in Figure 2.

**Figure 2.** Flowchart of the procedure to map snow cover duration (SCD) and snow cover area (SCA) in the Sugnugur catchment. MODIS is Moderate Resolution Imaging Spectroradiometer. NDSI is normalized difference snow index. GF and CA are gap-filling and conditional adjustments, respectively. Ts and Ta denote daily mean ground surface temperature (GST) and daily mean air temperature, respectively. SD and SWE are snow depth and water equivalent, respectively.

For all satellite images, we adopted the normalized difference snow index (NDSI) with a threshold of ≥0.4 to denote snow cover [40,42–44].

$$\text{NDSI} = \frac{\text{Green} - \text{SWIR}}{\text{Green} + \text{SWIR}} \tag{1}$$

where *Green* represents reflectance in a visible band and *SWIR* is reflectance in a short-wave infrared band.

*Estimating SCD from ground temperature*: The timing characteristic of snow cover was studied at the hydro-climatic station to understand the general snow accumulation and ablation processes in the study area. SCD is commonly derived from SD, SWE, or satellite data. In this study, a daily mean air temperature (Ta) of <0 ◦C and a daily mean albedo of ≥0.25 (unitless) were used to distinguish snow from the existing vegetation (Figure 3). We assume a persistent snow cover (continuous) if these two conditions are consistent for at least 14 days, otherwise it will be assigned as intermittent (not continuous). Occasional autumn and spring (transition periods) snowfalls develop intermittent snow cover lasting for a couple of days, whereas decreased air temperature towards the winter months (DJF) allows the snow cover to become persistent. We classified snow duration for each year as the sum of both intermittent SCD and persistent SCD.

**Figure 3.** Snow depth (SD) collected from the SR50A-L sensor at the hydro-climatic station site by coupling with daily air temperature and albedo. The values are shown as 5-day averages.

We also used surface temperatures collected from the hydro-climatic station as proxies to classify persistent snow cover, such that if the daily mean surface temperature (Ts) remained below −1 ◦C for at least 14 days, the snow cover was determined as persistent. These conditions were then applied to our iButton measurements, which allowed us to validate our reconstructed SCD across the catchment. We examined the sign and strength of the relationship between the persistent SCD determined from the Ta and albedo and that approximated from Ts using Pearson's correlation. We assumed that if a positive and significant relationship existed, then persistent SCD may be reconstructed from Ts. Winter 2012/2013 was an exceptional year with relatively low temperatures in autumn, which led to a high number of below −1 ◦C days. Therefore, we did a manual correction by selecting the days between the date when first snowfall was observed after the surface temperature dropped below the threshold and the date when albedo dropped to below 0.25.

*MODIS derived snow metrics:* The timing characteristic of snow over the study area was analyzed using the daily MODIS Terra and Aqua data by determining snow onset (So) and snow melt (Sm) dates, as well as SCD for snow seasons 2000–2017. We first separated the cloudy pixels, designated with a value of 250, from the cloud free pixels and then detected the snow cover for the cloud free pixels using the NDSI to create daily 500 m snow maps.

To reduce cloud coverage, we applied a series of cloud removal approaches to Terra and the combination of Terra and Aqua data separately. First, we applied a 1-day backward and forward gap-filling (temporal GF) approach, with a limit of 3 days, to both Terra and Aqua such that when a pixel is cloudy, it goes one day backward and forward to test for a cloud free pixel [11,44]. Following this temporal gap-filling, a conditional adjustment (CA), which takes spatial and seasonality corrections into account, was applied. The CA assumes that if the given pixel after the temporal gap-filling is still cloudy, it can be reclassified as snow or no-snow under certain circumstances; if two of the 8-neighboring pixels at lower elevation are cloud free and the day of year is within the range between So and Sm dates of that certain pixel, we assigned snow. The So and Sm dates were determined for each MODIS pixel as the first and the last day of the classified snow cover days within the hydrologic year; it also defined the length of snow season for each pixel.

After adjusting for cloud cover, gap-filled images from both Terra and the combination data still underestimated SCD over heterogeneous topography and forested areas because the coverage area of a single 500 × 500 m MODIS image includes a high range of elevations as well as both forested north- and steep south-facing slopes. Thus, we generalized the annual number of SCD from the MODIS data for each year over the study area based on annual SCD rates. To do so, 53 random points were selected, either on mountain peaks or lowland areas where no forests or shrubs exist and the topography is homogenous, in order to employ linear regression analysis between the SCD from the improved MODIS data and DEM; later the regression parameters were used for deriving the SCD rate for each year.

Trend analysis using linear least squares regression were also conducted for snow onset (So), snowmelt (S m) dates, and SCD from the improved MODIS daily snow-cover. Since air temperature is one of the main parameters that drives the snow accumulation and ablation processes, we calculated a trend analysis for monthly mean air temperature during the snow season, using observations from the nearest long-term climate station at Chinggis Khan International Airport (60 km south from the hydro-climatic station). The results can be found in Section 3.2.2.

*Landsat and Sentinel images to correct for land-cover and topography:* The generalized SCD is not realistic because it includes the intermittent SCA, which mostly occurs on the steep south-facing slopes and wind-blown surfaces. To identify these areas, we created a sequence of snow cover maps from the combination of Landsat-7, Landsat-8, and Sentinel-2A images (40 images in total) taken on clear-sky days. These time sequence maps are helpful in showing the development of SCA with high spatial resolution in the study area. Although the spatial resolutions of the satellites are high, snow cover below the forest canopy and dense shrub is difficult to detect [45]. Therefore, a correction method, which considers altitude dependency and land-cover types, was applied to all Landsat and Sentinel images to correct snow cover in those areas. More specifically, if a pixel from either a Landsat or Sentinel image was classified as snow at the hydro-climatic station, the areas in higher elevations with forests or shrubs were also classified as snow. If not, this assumption was ignored. By repeating this correction for each image, the change detection analysis was completed for the snow season period 2016–2017, representing the development of SCA in the catchment. Based on the changes in SCA, we separated the study area as persistent SCA (pixels classified as snow in >70% of the data) and intermittent SCA (pixels classified as snow in <30% of the data). The result can be found in Section 3.3. We assumed that the general spatial distribution of snow would be similar for each year to obtain the mean SCD for winters 2000/2001 to 2016/2017 over the seasonal persistent SCA by using the empirical altitude-dependent SCD rates that were found from the improved MODIS data.

*Result assessment:* To validate the results, we first compared the estimated total SCD with the observed total SCD at the hydro-climatic station. Daily surface temperature measurements (Ts) from the iButtons, distributed in the catchment, were also used by comparing the linear relationship, which was found at the hydro-climatic station, to that between the estimated persistent SCD days and days with < −1 ◦C of Ts from the iButtons (Figure 2). Photographs from the time-lapse camera and snow field measurements were also used to verify the spatial distribution of snow visually.
