Development of a Daily Cloud-Free Snow-Cover Dataset Using MODIS-Based Snow-Cover Probability for High Mountain Asia during 2000–2020

Abstract

Investigating the changes in snow cover caused by climate change is extremely important and has attracted increasing attention in cryosphere and climate research. Optimal remote sensing-based snow datasets can provide long-term daily and global spatial-temporal snow-cover distribution at regional and global scales. However, the application of these snow-cover products is inevitably limited because of the space–time discontinuities caused by cloud obscuration, which poses a significant challenge in snowpack-related studies, especially in High Mountain Asia (HMA), an area that has high-elevation mountains, complex terrain, and harsh environments and has fewer observation stations. To address this issue, we developed an improved five-step hybrid cloud removal strategy by integrating the daily merged snow-cover probability (SCP) algorithm, eight-day merged SCP algorithm, decision tree algorithm, temporal downscaling algorithm, and optimal threshold segmentation algorithm to produce a 21-year, daily cloud-free snow-cover dataset using two daily MODIS snow-cover products over the HMA. The accuracy assessment demonstrated that the newly developed cloud-free snow-cover product achieved a mean overall accuracy of 93.80%, based on daily classified snow depth observations from 86 meteorological stations over 10 years. The time series of the daily percentage of binary snow-cover over HMA was analyzed during this period, indicating that the maximum snow cover tended to change more dramatically than the minimum snow cover. The annual snow-cover duration (SCD) experienced an insignificantly increasing trend over most of the northeastern and southwestern HMA (e.g., Qilian, eastern Kun Lun, the east of Inner Tibet, the western Himalayas, the central Himalayas, and the Hindu Kush) and an insignificant declining trend over most of the northwestern and southeastern HMA (e.g., the eastern Himalayas, Hengduan, the west of Inner Tibet, Pamir, Hissar Alay, and Tien). This new high-quality snow-cover dataset will promote studies on climate systems, hydrological modeling, and water resource management in this remote and cold region.

1. Introduction

Snow cover is an essential geophysical factor in the global climate system [1,2] and has a significant influence on the water cycle [3,4]. The snow-albedo feedback, along with the ice-albedo feedback, is invoked as a leading cause of amplified warming in polar and mountainous regions [5,6,7] and, thus, plays a vital role in regulating climate change at both regional and global scales. Snowpack takes on special significance in mountain regions where snow stores enormous quantities of water, altering the ecologic and economic balance of regions far downstream by delaying the release of water until months after precipitation events [8]. By providing more than one-sixth of the world’s population with available natural water resources, snowmelt serves several essential uses for human beings, such as irrigation, industry, and hydropower generation [9,10]. Accurately monitoring snow cover and its variations is, thus, crucial for the studies mentioned above.

Covering the whole area of the Tibetan Plateau and a series of mountain chains [11] (Figure 1), the High Mountain Asia (HMA) (64°–107°E and 22°–47°N) is a hotspot for investigating cryospheric change against the background of global warming [12]. Since the radiative forcing of the snow cover in this middle-latitude region is more significant than that in high-latitude areas [13], the snow cover in this particular region can profoundly impact the surface energy budget and regulate the thermal state of the plateau [14] and, thus, significantly influence weather and climate at both regional and global scales [15,16]. As reliable quantitative snow datasets on the HMA region are also scarce, the impact of climate change on snowpack is still under debate [17]. Thus, collecting high-quality snow-cover datasets from this cold and remote region is urgent.

Ground observations, passive microwave remote sensing, and optical remote sensing are widely utilized to collect snow data. In situ observations collected by meteorological stations, hydrological stations, and snow investigations represent the most direct and credible snow records [18]. However, ground observations are less capable of providing spatiotemporally continuous snow-cover data at regional and global scales and, thus, do not fully satisfy the needs of the modeling community [19]. Moreover, since traditional field snow surveying is time-consuming, expensive, and highly challenging [20], remote sensing satellites have served to detect snow properties for more than four decades [21,22], e.g., Landsat and the Système Probatoire d’Observation de la Terre (SPOT) [23], along with the advanced very high-resolution radiometer (AVHRR) [24], moderate-resolution imaging spectroradiometer (MODIS) [25], scanning multichannel microwave radiometer (SMMR), special sensor microwave/imager (SSM/I) [26,27], and advanced microwave scanning radiometer–earth observing system (AMSR-E) [28]. Passive microwave-based snow-cover extent data have misclassification errors, owing to several factors (e.g., signal obstructions, snow grain size, land cover influences, and algorithm limitations) [29,30,31], and, thus, exhibit less snow cover than optical remote sensing observations from the National Environmental Satellite Data and Information Service (NESDIS) [32]. Optical remote sensing technology is an effective way to quickly acquire reliable snow-cover data, due to its acceptable spatial resolution and long period of use.

Index-based methods provide an effective and time-saving approach for identifying snow cover using optical remote sensing datasets. The snow index methods originated in the 1980s when Dozier (1989) used the green and short-wave infrared 1 (SWIR1) bands of Landsat data for snow-cover identification research [33]. Relying on the snow spectral characteristic that the reflectance of snow cover in the green band is significantly greater than that in the middle infrared band, Hall et al. (1995) and their collaborators first proposed the concept of the normalized difference snow index (NDSI), which also normalizes the green and SWIR1 bands, and mapped global snow cover based on MODIS images for the first time [34]. The SWIR1 band was replaced by the SWIR2 band in the NDSI algorithm using the MODIS Aqua sensor, owing to its malfunctions in the SWIR1 band [35]. In addition, other similar snow indices (e.g., the snow grain size index (SGI) [36], normalized difference snow and ice index (NDSII) [37], normalized difference forest snow index (NDFSI) [38], and normalized difference snow index with no water information (NDSI_nw) [39]) were developed and are widely used for identifying snow-cover information. However, clouds block the electromagnetic waves recording land surface features from being received by optical remote sensing instruments [40,41,42], making it challenging to identify snow cover under clouds to a large extent. Thus, snow product developers have no choice but to sacrifice the spatial and temporal continuity of data from gridded snow-cover products that are usually masked by cloud cover (e.g., the MODIS snow product, Visible Infrared Imaging Radiometer/National Polar-orbiting Partnership (VIIRS/NPP) snow product, and JAXA Satellite Monitoring for Environmental Studies (JASMES) AVHRR snow product).

To enhance the spatiotemporal continuity of the original cloud-contaminated snow-cover product, many researchers have made progress in recent decades, especially in collecting multisource observations. The combining of multisource observations is a typical strategy for reducing cloud obscuration, which usually makes full use of the available snow observations and can be grouped into two categories: combinations of ground and remote sensing observations and combinations of multisource snow-cover products. It is a valuable approach combining remote sensing and ground observations to identify snow-cover information in areas that are regularly masked by clouds [43,44,45,46]. However, this approach relies on the distribution and density of the ground observation network to a large extent [47]. It is, thus, less capable of retrieving snow-cover information in remote and inaccessible regions. Additionally, collecting and combining as many snow-cover products as possible from various kinds of satellite sensors (e.g., MODIS, the interactive multisensor snow and ice mapping system (IMS), FengYun-3, JASMES, the long-term series of daily snow depth data in China (WESTDC) and AMSR-E) should be encouraged to develop cloud-free snow-cover products [20,28,41,48,49,50]. However, many factors (e.g., the system inconsistency of different sensors, differences in mathematical algorithms, and daily weather and cloud conditions) may influence the quality of new snow-cover products derived from various instruments.

Many mathematical algorithms based on spatiotemporal statistics have been developed to eliminate cloud cover from optical snow products. These algorithms can be grouped into three categories: spatial filters, temporal filters, and hybrid methods. First, the spatial cloud removal method relies on the snow-cover information from neighboring pixels of the cloud obscuration pixel or the environmental association information (e.g., DEM, snow line, and land use cover) and achieves acceptable cloud removal results for snow-cover products [51,52,53,54]. The spatial cloud removal method can fill in gaps in small, scattered clouds. However, the method fails to fill in the gaps in massive cloud-covered areas due to the unavailability of spatial neighboring pixels [55]. Second, the temporal cloud-removal method exploits multiple complementary observations of the same scene at adjacent times to mitigate the cloud obscuration inherent in primary cloud-obscured snow-cover products and performs better than the spatial cloud-removal method [42]. The method works well in snow-stable periods but has limitations in snow-transitional periods because it takes advantage of the time-series correlations of snow cover over a short period. Third, considering the respective advantages of the abovementioned two methods, hybrid cloud-clearing methods combining spatial filters, temporal filters, and other complementary strategies have been widely used to eliminate existing cloud-obscured snow-cover products [54,55,56]. Most of these hybrid cloud-clearing methods are associated with high computational and time costs when developing long-term and large-scale snow-cover products.

Numerous studies have been conducted to produce continuous space-time snow-cover products at regional or global scales. In previous studies, less attention has been given to utilizing the snow-cover probability (SCP) during a given period to generate cloud-free snow-cover products. Thus, the SCP is likely to provide a new opportunity for developing cloud-free snow-cover products in snow-dominated regions. Due to the abovementioned concerns, the primary purposes of this study are as follows: (1) to develop a gap-filling strategy based on snow-cover probability to remove cloud cover pixels from the original MODIS snow-cover product; (2) to develop a long-term gap-filled snow-cover product over the HMA region using data from the 2000 to 2020 snow year (from 1 August to 31 July of the following year); and (3) to analyze the characteristics of the daily snow-cover fraction of HMA and annual snow-cover duration (SCD) changes over the HMA region during this period.

2. Materials

2.1. Meteorological Station-Based Snow Depth Observations

In situ observations are essential for evaluating the accuracy of snow-cover products. The Chinese Meteorological Administration (CMA) stations have no snow-cover observations; however, there are enough snow-depth observations for the purpose available from these meteorological stations. The snow depth is measured when more than 50% of the area around the station is covered with snowpack [57]. The daily snow-depth observations from 86 meteorological stations from September to May of the following year during 2001–2010 were collected and classified as snow cover or non-snow cover, with a threshold of 1 cm [58,59]. The binary snow-cover and non-snow-cover datasets from these meteorological stations were used to evaluate the accuracy of the original MODIS snow-cover product and the cloud-free snow-cover dataset developed in this study.

2.2. MODIS Snow Product

The moderate-resolution imaging spectroradiometer (MODIS) snow-cover data from 2000 to 2020 were acquired from the National Snow and Ice Data Center (http://nsidc.org, accessed on March 2022). The MODIS Terra and Aqua satellites revisit global land cover twice daily and can provide land cover image information with daily temporal resolution. Moreover, since the overpass times of MODIS Terra and Aqua are at approximately 10:30 a.m. equatorial crossing time and 1:30 p.m. equatorial crossing time, this snow-cover product comprises two sub-datasets named the Terra (MOD10C1, Version 6) and Aqua (MYD10C1, Version 6) products, including snow-cover percentage, cloud obscuration percentage, and clear index (CI) data ranging from 0% to 100%. The snow-cover percentage and CI data were generated using a binning algorithm and were used to estimate the SCP in this study [60].

2.3. Downscaling Snow Depth Dataset

A daily snow-depth dataset for the Tibetan Plateau, developed by an improved spatial-temporal downscaling method from the National Tibetan Plateau Data Center, was utilized in this study (https://data.tpdc.ac.cn/, accessed on October 2021) [61]. The snow depth dataset outside the Tibetan Plateau area over the HMA region was produced using the aforementioned downscaling method presented in this study. This snow depth product is an improved version of the WESTDC product developed with Che’s algorithm [27]. The snow-cover probability generated from the MODIS fractional snow-cover product enhanced the original 0.25° passive microwave snow depth product. As a result, the spatial resolution of this improved snow-depth product improved to 0.05°.

3. Methods

This hybrid cloud removal method includes five steps (Figure 2): (1) daily merged snow-cover probability (SCP), (2) the 8-day merged SCP [61], (3) a decision tree algorithm, (4) the temporal downscaling algorithm of the 8-day SCP to generate daily cloud-free SCP, and (5) the threshold segmentation algorithm of daily cloud-free SCP. The final version of the cloud-free SCP was generated using a combination of the available SCP data without cloud contamination in step (1) and the gap-filling SCP product in step (3) and step (4) (Figure 2). The cloud-free snow-cover product was generated based on the final version of the cloud-free SCP and the optimal threshold value in step 5. This study used Python 2.7 and 3.9 to process the above-mentioned datasets in these steps.

3.1. Step 1: Daily Merged SCP

The ratio of the number of binary snow maps to the number of binary cloud-free maps during a half-month was employed to estimate the SCP [54]. However, the binary snow cover and cloud-free images that are separated by a fixed threshold provide little information, compared to the FSC and CI images ranging from 0% to 100%. Thus, replacing binary snow and cloud-free pixels with FSC and CI pixels is an improved strategy for estimating daily and 8-day SCP [61]. The daily available FSC and CI data from Terra and Aqua without cloud cover were thus used to estimate the daily SCP. If the daily clear index from the MOD10C1 or MYD10C1 was greater than 0 within a single day, the SCP was calculated using the ratio of the FSC_Terra to the CI_Terra when the Terra product was available from 2000 to 2002 and using the ratio of the sum of the FSC_Aqua and FSC_Terra to the sum of the CI_Aqua and CI_Terra when both the Terra and Aqua data were available for 2002–2020 (Equation (1)):

$$\mathrm{SCP}=\left\{\begin{array}{c}{\displaystyle \frac{{\mathrm{FSC}}_{\mathrm{Terra}}}{{\mathrm{CI}}_{\mathrm{Terra}}}}, \mathrm{if} \mathrm{only} \mathrm{Terra} \mathrm{is} \mathrm{available} \mathrm{from} 2000 \mathrm{to} 2002\\ {\displaystyle \frac{{\mathrm{FSC}}_{\mathrm{Terra}}+{\mathrm{FSC}}_{\mathrm{Aqua}}}{{\mathrm{CI}}_{\mathrm{Terra}}+{\mathrm{CI}}_{\mathrm{Aqua}}}}, \mathrm{if} \mathrm{both} \mathrm{Terra} \mathrm{and} \mathrm{Aqua} \mathrm{are} \mathrm{available} \mathrm{after} 2002\end{array}\right.$$

(1)

where SCP is the daily snow-cover probability, FSC_Terra and CI_Terra are the fractional snow cover and clear index from the original daily MOD10C1 snow-cover product, and FSC_Aqua and CI_Aqua are the fractional snow cover and clear index from the original daily MYD10C1 snow-cover product.

3.2. Step 2: 8-Day Merged SCP

The sum of the FSC and the sum of the clear indices over 8 days were calculated based on the MODIS Terra and Aqua products (Equations (2) and (3)). If the sum of the daily clear indices was greater than 0 within 8 days, the SCP during these 8 days was estimated using the ratio of the FSC_sum to the CI_sum (Equation (4)):

$${\mathrm{CI}}_{\mathrm{sum}}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{CI}}_{\mathrm{i}}$$

(2)

$${\mathrm{FSC}}_{\mathrm{sum}}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{FSC}}_{\mathrm{i}}$$

(3)

$${\mathrm{SCP}}_{8-\mathrm{day}}={\displaystyle \frac{{\mathrm{FSC}}_{\mathrm{sum}}}{{\mathrm{CI}}_{\mathrm{sum}}}}={\displaystyle \frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{FSC}}_{\mathrm{i}}}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{CI}}_{\mathrm{i}}}}$$

(4)

where CI_sum is the sum of daily clear index data over 8 days; FSC_sum is the sum of fractional snow cover over 8 days; SCP_8-day is the 8-day cloud-free snow-cover probability.

The above process can effectively remove most clouds from the original MODIS snow-cover product in 8 days. If a pixel was utterly (100%) covered by clouds all the time within an 8-day period, the SCP of these pixels was estimated using the cloud-free SCP for the preceding 8 days (PSCP_8-day) and that of the following 8 days (FSCP_8-day) as follows (Equation (5)). As a result, a cloud-free SCP_8-day product was produced for 2000–2020.

$${\mathrm{SCP}}_{8-\mathrm{day}}=\left\{\begin{array}{c}{\displaystyle \frac{{\mathrm{PSCP}}_{8-\mathrm{day}}+{\mathrm{FSCP}}_{8-\mathrm{day}}}{2}},\mathrm{both} \mathrm{PSCP} \mathrm{and} \mathrm{FSCP} \mathrm{are} \mathrm{available}\\ {\mathrm{PSCP}}_{8-\mathrm{day}},\mathrm{if} {\mathrm{FSCP}}_{8-\mathrm{day}} \mathrm{is} \mathrm{not} \mathrm{available}\\ {\mathrm{FSCP}}_{8-\mathrm{day}},\mathrm{if} {\mathrm{PSCP}}_{8-\mathrm{day}} \mathrm{is} \mathrm{not} \mathrm{available}\end{array}\right.$$

(5)

3.3. Step 3: A Decision Tree Algorithm

A decision tree algorithm was utilized to generate a cloud-free daily SCP product using an 8-day merged SCP and temporal downscaling algorithm. Specifically, if the pixels of the 8-day SCP were 1 and 0, the pixels influenced by clouds during this 8-day period were assigned 1 and 0. If the pixel values of the SCP_8-day ranged from 0% to 100%; and the daily snow depth was greater than 0, the daily SCP in the cloud-covered pixels was estimated by multiplying the SCP_8-day by 8 times and each daily temporal weight (Equation (6)):

$${\mathrm{SCP}}_{\mathrm{i}}=\left\{\begin{array}{c}1, {\mathrm{SCP}}_{8-\mathrm{day}}=100\%\\ 0, {\mathrm{SCP}}_{8-\mathrm{day}}=0\\ \begin{array}{c}{\mathrm{SCP}}_{8-\mathrm{day}}\times 8\times {\left({\mathrm{W}}_{\mathrm{tSD}}\right)}_i, 0<{\mathrm{SCP}}_{8-\mathrm{day}}<100\% \end{array}\end{array}\right.$$

(6)

where SCP_i is the daily snow-cover probability on the ith day over 8 days, (W_tSD)_i is the temporal weight of the ith day from the daily snow depth product over 8 days, and SCP_8-day is the 8-day cloud-free snow-cover probability.

3.4. Step 4: Temporal Downscaling Algorithm for the 8-Day SCP

Passive microwave instruments have exhibited a great ability to monitor Earth’s surface changes under clouds [62,63]. Moreover, previous studies have demonstrated a positive relationship between FSC and snow depth data over the TP [64,65], indicating the potential of the snow depth to reveal temporal variations in optical remote sensing-based FSC pixels. Moreover, passive microwave-based snow-depth products have great potential for the temporal downscaling of snow-cover probability for periods of several days in this snow-dominated region. It is, thus, believed that dividing the sum of the 8-day snow depth data by the daily snow depth data can be used to calculate the daily temporal weight of 8-day SCP information. In this case, the 8-day SCP grids must be multiplied by 8 before being multiplied by the temporal weight value. The equations of the temporal downscaling algorithm are as follows (Equations (7) and (8)):

$${\left({\mathrm{W}}_{\mathrm{tSD}}\right)}_{\mathrm{i}}={\displaystyle \frac{{\mathrm{SD}}_{\mathrm{i}}}{{{\displaystyle \sum }}_{\mathrm{i}=1}^8{\mathrm{SD}}_{\mathrm{i}}}}$$

(7)

$${\mathrm{SCP}}_{\mathrm{i}}={\mathrm{SCP}}_{8-\mathrm{day}}\times 8\times {\left({\mathrm{W}}_{\mathrm{tSD}}\right)}_{\mathrm{i}}$$

(8)

where SCP_i is the daily snow-cover probability on the ith day, W_tSD is the temporal weight from the daily snow depth product, and SD_i is the snow depth on the ith day over 8 days (1 ≤ i ≤ 8).

3.5. Step 5: Threshold Segmentation of Daily Cloud-Free SCP

The sensitivity of individual criteria tests can be studied by changing the threshold value incrementally and then analyzing the effect on the results [34]. Previous studies employed various thresholds to extract snow cover from the background in NDSI images over different regions [34,66,67]. The 0.29 was considered the optimal NDSI threshold for generating binary snow cover over the Tibetan Plateau [12,67]. Since the SCP was estimated based on fractional snow cover, the NDSI threshold (0.29) was converted to the SCP threshold (0.41), using the NDSI-based fractional snow-cover algorithm [35].

3.6. SCD Estimation

The SCD was calculated as the number of days in which a pixel is covered by snow in a snow year or a season, as follows [12,68]:

$$\mathrm{SCD}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{SC}}_{\mathrm{i}}$$

(9)

where SCD is the snow-cover duration during a given snow year or a season; SC_i is the daily snow-cover data on the ith day; n is the number of days in a snow year or a season.

3.7. Validation of the Daily Snow-Cover Dataset

The confusion matrix is utilized to evaluate the performance of the newly developed snow-cover data and the meteorological station-based snow depth (>1 cm) observations. The freshly developed daily snow-cover dataset was quantitatively evaluated against the meteorological station observations using several parameters, such as overall accuracy (OA), underestimation error (UE), and overestimation error (OE), based on a confusion matrix, as follows:

$$\mathrm{Overall} \mathrm{accuracy} =({\mathrm{N}}_{11}+{\mathrm{N}}_{22})/{\mathrm{N}}_{\mathrm{t}}$$

(10)

$$\mathrm{Underestimation} \mathrm{error}={\mathrm{N}}_{12}/{\mathrm{N}}_{\mathrm{t}}$$

(11)

$$\mathrm{Overestimation} \mathrm{error} ={\mathrm{N}}_{21}/{\mathrm{N}}_{\mathrm{t}}$$

(12)

where N₁₁, N₁₂, N₂₁, and N₂₂ are the sample counts of the snow or non-snow resulting from the comparison between the binary snow-cover product and station-based data, based on the confusion matrix (

Table 1. Confusion matrix of the sample points for the accuracy assessment of the cloud-free snow-cover product. The classified binary images are in rows, and the corresponding reference points are in columns. N₁₁, N₁₂, N₂₁, and N₂₂ are the sample counts between the classified binary image and the reference point.

Class	Snow from Product	Non-Snow from Product	Total
Snow from station	N11	N12	N1j
Non-snow from station	N21	N22	N2j
Total	Ni1	Ni2	Nt

); N_t is the total sample count.

4. Results

4.1. Cloud Removal Results

An example of a cloud removal result for each step on 1 January 2004 is shown in Figure 3e–h. Figure 3a,c shows the original available MODIS-terra and MODIS-aqua fractional snow-cover data. Figure 3b,d shows the original MODIS-terra and MODIS-aqua fractional cloud-cover data. The ratio of the sum of the FSC_Aqua and FSC_Terra to the sum of the CI_Aqua and CI_Terra during a day was used to estimate the daily merged SCP. Combining the two sensors somewhat improved the ability to detect snow-cover information. If the 8-day SCP was equal to 0 or 1, the 8-day SCP was used to fill the cloud gaps of the daily merged SCP in the corresponding pixels (Figure 3f,g). To further reduce the existing cloudy pixels, a temporal downscaling algorithm was used to estimate the daily SCP based on the 8-day SCP if the SCP was between 0 and 1 (0 < SCP < 1) (Figure 3h). The final classified binary snow-cover extent based on cloud-free SCP in HMA is illustrated in Figure 3i.

To further test the performance of the classified gap-filled snow cover, a comparison was made between the classified gap-filled daily SCP product and the classified original MODIS-terra and MODIS-aqua snow-cover products (Figure 4). The spatiotemporal continuity of the newly developed snow-cover product improved to a great extent. As a result, a cloud-free daily snow-cover product was successfully generated based on the abovementioned steps, and is shown on the right side of Figure 4.

4.2. Validations of the Binary Snow-Cover Products

To test the performance of the new snow-cover product at each station, the spatial distributions of the overall accuracy, underestimation error, and overestimation error of the new snow-cover product over these 86 stations were analyzed and are shown in Figure 5. More than 230,000 station-based binary snow-cover observations were collected from 86 stations for 2001–2010 and these were used to assess the accuracy of the new gap-filled snow-cover product. It is easy to see that the overall accuracy of the new snow-cover product is relatively low in the central regions of these stations at high altitudes, due to the underestimation and overestimation errors. In contrast, the newly developed snow-cover product maintains high overall accuracy and low underestimation and overestimation errors at low altitudes. The overall accuracy of this new snow-cover product decreases with increasing altitude (

Table 2. The overall accuracy, underestimation error, and overestimation error of the new cloud-free snow-cover product in the different elevation zones.

Altitude (m a.s.l)	Overall Accuracy	Underestimation Error	Overestimation Error
1000–2500	96.03%	1.03%	3.12%
2500–3000	94.82%	2.27%	3.02%
3000–3500	95.35%	4.03%	2.73%
3500–4000	92.71%	4.52%	2.93%
4000–4800	92.21%	4.88%	3.07%
All stations	93.80%	3.37%	2.97%

). The overall accuracy of this new cloud-free snow-cover product is highest in autumn (93.51%), followed by winter (93.02%) and spring (92.50%) (

Table 3. The overall accuracy, underestimation error, and overestimation error of the new cloud-free snow-cover product in the different seasons (autumn, winter, and spring).

Season	Overall Accuracy	Underestimation Error	Overestimation Error
Autumn	93.51%	2.79%	3.69%
Winter	93.02%	3.83%	3.16%
Spring	92.50%	4.97%	2.53%
All seasons	93.80%	3.37%	2.97%

). Notably, the mean value of the underestimation error (3.37%) is slightly greater than that of the overestimation error (2.97%) for these 86 stations, indicating that the new product maintains a relative balance between the underestimation error (0–14.88%) and overestimation error (0–12.31%).

Figure 6 compares the overall accuracy, underestimation error, and overestimation error of the new cloud-free snow-cover products under different weather conditions. The new snow-cover product achieves overall accuracies of 87.12% and 97.47%, respectively, under cloudy sky (CI < 50%) and clear-sky (CI > 50%) conditions. Notably, the overall accuracy of the snow-cover products is very sensitive to the fraction of cloud cover. This is a satisfactory result since the new cloud-free snow-cover product achieves an overall accuracy ranging from 82% to 99% under all sky conditions.

4.3. Spatial Pattern of Monthly Snow Cover over the HMA

Figure 7 illustrates the spatial distribution of monthly snow cover over the HMA. The seasonal maximum snow cover occurs frequently in western, southeastern, and eastern HMA during the cold seasons. The persistent snow cover is distributed in high-altitude areas of the western and southern TP in all seasons. The precipitation from the westerly and Indian monsoons continuously spills over the frontier ranges of the Pamir and Himalayan chains and maintains the snow-cover distribution in this region [15,69]. As a result, snow cover often occurs at the edge of HMA, which agrees with the identified distribution of mountain chains in this alpine region. Snow-cover distribution is, thus, profoundly influenced by precipitation and terrain.

4.4. Time Series of Daily Snow-Cover Fractions in HMA

Figure 8 illustrates the ratio of binary snow cover to the whole of the HMA region during 2000–2020. The multiyear mean daily snow-cover fraction over the HMA region experiences an increasing trend from July to December and a decreasing trend from December to June. Although the multiyear maximum and minimum daily snow-cover fractions also show similar trends, the former fluctuates more dramatically. At the subregion scale, the snow-cover fraction in the western HMA (e.g., Karakoram, Pamir, the Hindu Kush, Tien, Hissar Alay, and the western Himalayas) is greater than that in the inner and eastern HMA. Note that the snow-cover fraction in the inner and eastern HMA regions changes more dramatically than that in the western HMA region.

Additionally, the daily snow-cover fraction from 2000 to 2020 shown in Figure 9 gives more detailed information on snow cover changes from year to year. Although the minimum snow-cover fraction was approximately 5%, the peak of the snow-cover fraction fluctuated between 40.56% and 71.86%, suggesting that the annual maximum snow cover tended to change more dramatically than the annual minimum snow cover during this period.

4.5. Spatiotemporal Trend of the Annual SCD during 2000–2020

Since the annual SCD exhibits significant spatial heterogeneity due to its particular geographic location, analyzing the annual SCD trend at the pixel scale is, thus, encouraged over the HMA region during the study period (Figure 10). At the interannual scale, as shown in Figure 10a, an increasing trend was observed across a large proportion of the northeastern and southwestern HMA (e.g., Qilian, east Kun Lun, east of Inner Tibet, the western Himalayas, central Himalayas, and the Hindu Kush). A decreasing trend was observed over most of northwestern and southeastern HMA (e.g., the eastern Himalayas, Hengduan, west of Inner Tibet, Pamir, Hissar Alay, and Tien). Note that most of these pixels over the HMA showed an insignificant trend during 2000–2020 (p > 0.05). The trends of the annual SCD in autumn, winter, and spring (Figure 10c,d) were roughly similar to those of the annual SCD (Figure 10a). Note that the annual SCD in autumn over some parts of the southern HMA experienced a significant declining trend during this period (Figure 10c,d).

5. Discussion

It should be noted that the accuracy assessment was conducted mainly in the eastern and southern HMA regions due to the uneven distribution of station-based observations in the HMA region. The accuracy of this new snow-cover product in the western and inner HMA regions should, thus, be evaluated in further studies. High spatial resolution optical remote sensing images (e.g., the Landsat series) have great potential for assessing medium-resolution snow-cover products and often serve to evaluate the relative accuracy of MODIS snow-cover products. However, conducting reliable accuracy assessments of MODIS-based binary snow-cover products is complicated when based on high-resolution remote sensing images in the central Tibetan Plateau and the Qaidam Basin, where shallow and fragmented snowpacks are distributed randomly and widely. This is because shallow and fragmented snowpacks are omitted in MODIS-based binary snow-cover products. The observation problem of shallow and fragmented snow-cover information may be solved with the help of high-quality optical images from unmanned aerial vehicles and high-performance commercial satellites at the submeter scale.

Although the new snow-cover product developed in the present study achieved an overall accuracy of approximately 94% when the station-based snow depth of greater than 1 cm was considered snow cover, one should note that the station-based snow-depth threshold may affect the snow-cover accuracy to a certain extent. Previous studies have used various in situ snow depth thresholds to evaluate the accuracy of binary snow-cover products. For example, Hao et al. (2022) used 362 CMA station-based snow depth observations (> 1 cm) during snow seasons to calculate the overall accuracy of 3 datasets [58]. Pan et al. (2024) used snow-depth observations that were greater than 3 cm from 175 stations to assess binary snow-cover products and achieved an overall accuracy of 93.26% [70]. A hidden Markov random field model-based cloud-free snow-cover product over the Tibetan Plateau showed an overall accuracy of 98.29% in comparison with in situ snow-depth observations (> 3 cm) [71]. Yu et al. (2016) used multiple snow-depth thresholds to classify station-based snow depth observations for validating the accuracy of snow-cover products over the Tibetan Plateau, indicating that the snow-depth threshold influences the accuracy of snow-cover products to a certain extent [72]. Thus, it is not easy to compare the accuracy of different snow-cover products if a variety of in situ snow-depth thresholds are used. This is because the overall accuracies of the snow-cover products tended to increase overall with increasing snow-depth thresholds [72,73]. It is, thus, suggested that the accuracy of snow-cover products should be compared when the station-based snow-depth thresholds are the same.

6. Conclusions

This present study designed an enhanced 5-step hybrid cloud removal strategy and developed a 21-year, daily gap-filled snow-cover product over the HMA region. The accuracy assessment showed that the new cloud-free snow-cover product achieved a mean overall accuracy of 93.8% (ranging from 82–99%), suggesting the reliable and stable performance of this new snow-cover product. The annual SCD change over this remote and cold region showed significant spatial heterogeneity, with an increasing trend over a large proportion of the northeastern and southwestern HMA, and a decreasing trend over most of the northwestern and southeastern HMA. This snow-cover product, covering the whole HMA region, is freely available to the scientific community of the cryosphere. The development of this snow-cover product over the HMA will, thus, enrich snow data sources in this remote and inaccessible region and will further promote snowpack-related studies in the future.