No evidence of widespread decline of snow cover on the Tibetan Plateau over 2000–2015

Understanding the changes in snow cover is essential for biological and hydrological processes in the Tibetan Plateau (TP) and its surrounding areas. However, the changes in snow cover phenology over the TP have not been well documented. Using Moderate Resolution Imaging Spectroradiometer (MODIS) daily snow products and the Interactive Multi-sensor Snow and Ice Mapping System (IMS) data, we reported daily cloud-free snow cover product over the Tibetan Plateau (TP) for 2000–2015. Snow cover start (SCS), melt (SCM) and duration (SCD) dates were calculated for each hydrological year, and their spatial and temporal variations were analyzed with elevation variations. Our results show no widespread decline in snow cover over the past fifteen years and the trends of snow cover phenology over the TP has high spatial heterogeneity. Later SCS, earlier SCM, and thus decreased SCD mainly occurred in the areas with elevation below 3500 m a.s.l., while regions in central and southwestern edges of the TP showed advanced SCS, delayed SCM and consequently longer SCD. The roles of temperature and precipitation on snow cover penology varied in different elevation zones, and the impact of both temperature and precipitation strengthened as elevation increases.

mountainous areas. However, Ngari zone, the south region of the Kunlun Mountains and the eastern part of TP had late SCS, where snow usually occurred after DOY 365. We found that sd of SCS increased from southeast to northwest. The highest sd of SCS were found at the eastern section of the Kunlun Mountains and the western section of the Gangdise Mountains.
The spatial distribution of SCM showed a contrasting shape to SCS ( Fig. S1e and f). Regions with high elevation had the latest SCM dates, and while SCM overall occurred earlier with decreasing elevation. Specifically, snow melting dates were later than DOY 165 in the western borders of the TP, e.g., the west Kunlun Mountains, Nyainqentanglha Mountains, and several high elevation areas in the central plateau. In the vast interior of the TP, SCM ranged between DOY of 60 to 90. The earliest SCM (before DOY 45) were found in the north of Qangtang and south of Qinghai. Most regions had small sd of SCM within 45 days and the regions with largest sd (above 45 days) were located in the south border of the TP.
We further investigated snow cover phenology along the altitude gradient. As shown in Fig.  S2, SCS substantially advanced, SCM delayed, and SCD decreased with an increasing altitude. In particular, these trends became much more obviously above the elevation of 5000 m a.s.l.
Text S2: Discussion of the algorithm and its effectiveness of cloud removal.
MODIS daily snow cover products have been widely used to monitor snow cover dynamics for various regions in recent years 1-4 . Unfortunately, cloud contamination greatly limits the application of these products. We firstly applied two previously reported steps: combining Terra and Aqua products and adjacent temporal combination successively 1,5 . However, not all cloud pixels can be removed after these two steps. Some researches chose to increase the days of temporal combination 6,7 , but in this way we would reduce the temporal resolution and might thus influence the accuracy of snow phenology date retrieval. Passive microwave sensors can penetrate cloud cover while they also provide the information of snow depth and snow water equivalent. Thus, several studies used passive microwave data to remove cloud pixels 8,9 . However, the spatial resolution of these data is 25 km, which is too coarse compared to MODIS data. This may also lead to rather low spatial resolution in the composite images if a large number of cloud pixels exist. Due to the aforementioned challenges, we proposed a new algorithm which employs the IMS data to remove the residual cloud pixels. This process effectively removed all the cloud pixels and kept higher spatial and temporal resolution as well. Yu et al 10 developed daily cloud-free snow products through combining MODIS Terra-Aqua and IMS, and the overall accuracy of their products is 94%. Our algorithm added adjacent temporal combination, which could further reduce the proportion of cloud pixels, before fusing MODIS and IMS. Consequently our products should have higher accuracy in snow identification.
We took the data of 1 February 2011 as an illustration example for the effectiveness of cloud removal approach. Both MOD10A1 and MYD10A1 had a large number of cloud pixels, accounting for 32% and 43.3% of the whole plateau respectively, and there was no cloud pixel in IMS data (Fig.  S3a, b and c). After the combination of MOD10A1 and MYD10A1, the percentages of cloud pixels were reduced to 22.6% (Fig. S3d). The adjacent temporal combination could further eliminate partial cloud coverage. However, there were still 16.3% of the pixels covered by cloud (Fig. S3e). Then, based on the IMS data, all the remaining cloud pixels were removed, and the resulting cloudfree data is shown in figure S3 f. This method effectively curtailed the problem of cloud interference as compared to the original MODIS data. Thus, it is more accurate to use the cloud-free data to analyze the spatial and temporal patterns of the snow cover phenology over the TP.

Text S3. Accuracy assessment of original MODIS snow products and our daily cloud-free snow map.
We evaluated the accuracy of the original MODIS snow products and our daily cloud-free snow map based on daily snow depth data from meteorological stations. We calculated the snow accuracy and overall accuracy in both clear-sky and all-sky conditions (Table S2). The overall accuracy of both MOD10A1 and MYD10A1 is more than 97% in clear-sky condition. However, it decreases to about 50% in all-sky conditions, because of the effect of cloud. After the process of cloud removal, the overall accuracy was significantly improved to 96.6% in all-sky condition.