2.1 Data

MODIS sensors onboard Terra and Aqua satellites provide the global snow cover datasets. The daily MOD10A1 and MYD10A1 datasets of C6 are available through the website of the National Aeronautics and Space Administration (NASA; https://search.earthdata.nasa.gov/, last access: 6 July 2022). The NDSI_Snow_Cover (hereafter referred to as NDSI) scientific dataset, with a range of 0, 10 to 100, was used in this study. As shown in Fig. 1, the NDSI of 19 tiles covering China (excluding sea area) from 1 August 2000 to 31 July 2020, was acquired to generate snow cover maps. The 90 m digital elevation model (DEM) dataset of the Shuttle Radar Topographic Mission (SRTM) was obtained from the United States Geological Survey (USGS). In addition, two existing cloud-free snow cover products, including the daily snow cover extent at a 5 km resolution (NIEER AVHRR SCE, https://doi.org/10.11888/Snow.tpdc.271381, Hao et al., 2021) and the daily cloud-gap-filled MODIS NDSI at a 500 m resolution (MODIS CGF NDSI/MOD10A1F, https://doi.org/10.5067/MODIS/MOD10A1F.061, Hall and Riggs, 2020), were used for comparison. For reference data, the snow-depth measurements respectively derived from 49 and 102 meteorological stations in NX and QTP (Tibet Meteorological Bureau and National Meteorological Information, 2018) were used for station-based validation. Since the snow depth data can only assess the classification performance of MODIS NDSI retrievals, the NDSI maps derived from Landsat OLI images were utilized for comprehensive validation.

Figure 1Topographic relief of China, meteorological stations in NX and QTP, and Landsat OLI scenes used for validation.

2.2 Algorithm description

MODIS NDSI datasets are unable to represent the daily conditions of snow accumulation and ablation accurately because the optical remote-sensed images are subject to severe cloud pollution. Therefore, STAR, which is derived from our two-stage spatiotemporal fusion method (Jing et al., 2019), is presented to produce a spatiotemporal continuous snow collection. As shown in Fig. 2, the generation procedure comprises the pre-process TAC and the key-process STAR. Then, a quality assessment (QA) approach is presented to provide a data reliability profile for users. On this basis, post-processing is used to further improve the data quality in individual abnormal areas.

Figure 2Schematic representation of the generation procedure of STAR NDSI collection.

Download

2.2.2 Spatio-Temporal Adaptive fusion with erroR correction (STAR)

Many regions with persistent clouds are out of the scope of TAC. To this end, an advanced STAR method, which comprehensively utilizes spatiotemporal contextual information, is proposed to remove the clouds thoroughly. As shown in Fig. 3, the method performs in two passes: spatio-temporal adaptive fusion (STAF) and error correction (EC).

Figure 3Detailed flowchart of the Spatio-Temporal Adaptive fusion with erroR correction (STAR).

The first pass involves the generation of new NDSI maps by adaptively merging the spatiotemporal contextual information, including space partition, adaptive space–time block determination, and Gaussian kernel function (GKF)-based fusion. The research area is first segmented into dozens of partitions considering the spatial heterogeneity of snow patterns. Thus, the subsequent processes can be performed on a partition basis. Moreover, the optimal query partitions (Q) to each target partition (T) are determined by a comprehensive consideration of temporal distance (t), regional correlation (r), and cloud-free fraction (f) concerning the temporal complexity of snow variations. The following optimal parameters are derived from the extensive experiments:

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{l}\text{Scheme 1:}r>\mathrm{0.7},\text{if}{f}^{\mathrm{C}\mathrm{\&}\mathrm{T}}>\mathrm{0.3}\\ \text{Scheme 2:}max(t^{-\mathrm{1}}+f^{\mathrm{C}}),\text{others}\end{array}\right.,\end{array}$$

where the regional correlation between the candidate partition (C) within an 8 d window and the target partition is considered representative if the fraction of the intersecting cloud-free areas (f^C&T) is higher than 0.3. The candidate partition is then determined as a query partition according to Scheme 1 when the regional correlation is larger than 0.7. Otherwise, Scheme 2 is activated. Two query partitions, with short distance and high cloud-free fraction, are identified within the preceding 8 d and the backward 8 d, respectively. Subsequently, the 3×3 neighborhoods for each pixel of the target partition in all the associated query partitions are determined as the space–time reference block. Lastly, each pixel is reassigned a fused value from the related space–time block, as expressed in Eq. (3):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\text{NDSI}}_i^{\mathrm{F}}={\sum }_{t=\mathrm{1}}^M{\sum }_{s=\mathrm{1}}^N{w}_{(i,\mathrm{st})}\times {\text{NDSI}}_{(i,\mathrm{st})}^{\mathrm{Q}},\\ & \text{where }{W}_{(i,\mathrm{st})}=r_{\mathrm{t}}^{\mathrm{2}}\times \exp \left({\displaystyle \frac{-\left(\right(\mathit{\epsilon }\times \mathrm{\Delta }{s}_{(i,\mathrm{s})}{)}^{\mathrm{2}}+\mathrm{\Delta }{t}_{(i,\mathrm{t})}^{\mathrm{2}}\left)\right)}{\mathrm{2}\times {\mathit{\sigma }}^{\mathrm{2}}}}\right),\end{array}\end{array}$$

where ${\text{NDSI}}_i^{\mathrm{F}}$ denotes the fused NDSI of pixel i in the target partition. The pre-processed NDSI is ${\text{NDSI}}_{(i,\mathrm{st})}^{\mathrm{Q}}$ in associated query partitions, and M is the number of query partitions, each of which contains N reference pixels. In addition, the weight W_(i,st) is assigned by a two-dimensional GKF involving the spatial distance (Δs_(i,s)) and the temporal distance (Δt_(i,t)), which is then normalized to w_(i,st). The standard deviation of GKF is σ. The dimensional difference, which is characterized by ε, is equal to ${\mathit{\sigma }}_{\mathrm{t}}/{\mathit{\sigma }}_{\mathrm{s}}$ with an expression of each single-dimensional GKF. The regional correlation between the query and target partitions is represented by r_t if Scheme 1 works; otherwise, it is ignored (i.e., r_t=1). The constant term ($\mathit{\epsilon }/\left(\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}\right)$) of GKF is ignored due to the normalization process. The important parameters in STAF are listed in Table 2.

Table 2Description and default values of STAF parameters.

Download Print Version | Download XLSX

The second pass corrects the fused NDSI maps considering the spatial correlation within a partition. Specifically, the residual errors of the intersecting cloud-free areas of the pre-processed and fused NDSI maps (NDSI^P after TAC and NDSI^F after STAF) are diffused to other cloud-free areas of the fused NDSI maps using the triangulation-based natural neighbor interpolation (NNI; Sibson, 1981). Then, the high-quality NDSI maps (NDSI^H) can be generated by removing all errors from the fused NDSI maps. The process is formulated as follows:

$$\begin{array}{}\text{(4)}& \left\{\begin{array}{l}{E}_{\mathrm{R}}={\text{NDSI}}_{\mathrm{R}}^{\mathrm{F}}-{\text{NDSI}}_{\mathrm{R}}^{\mathrm{P}}\\ E_{\mathrm{T}\left(i\right)}={\sum }_{n=\mathrm{1}}^{N^{\prime }}{\mathit{\varphi }}_{(i,n)}{E}_{\mathrm{R}(i,n)}\\ {\mathrm{NDSI}}_{\mathrm{T}}^{\mathrm{H}}={\mathrm{NDSI}}_{\mathrm{T}}^{\mathrm{F}}-E_{\mathrm{T}}\end{array}\right.,\end{array}$$

where “R” indicates the reference area which is the boundary of the intersecting cloud-free areas. The target area is indicated by “T”. The dynamic weights in the error diffusion from E_R to E_T are based on the Voronoi diagrams. As expressed in Fig. 3b (left), the original Voronoi cells (bounded by red and gray solid lines) of the reference pixels (gray dots) intersect with the new Voronoi cells (bounded by blue and gray solid lines) of the reference and target pixels. Taking the target pixel T₁ with the reference pixel R₁ as an example, the weight is assigned as the ratio of the area of the intersecting Voronoi cell (A_dabch) to that of the new Voronoi cell (A_defgh).

$$\begin{array}{}\text{(5)}& {\mathit{\varphi }}_{(\mathrm{1},\mathrm{1})}={\displaystyle \frac{A_{\mathrm{dabch}}}{A_{\mathrm{defgh}}}}.\end{array}$$

After all the partitions are processed in sequence, the next iteration of STAR begins until the clouds are completely removed.

2.2.3 Quality assessment (QA) approach

A revised QA approach for the gap-filled NDSI collection is proposed on the basis of the quality estimate of MODIS NDSI datasets (Riggs and Hall, 2015), and an example is presented in Fig. A1 (Appendix A). Users can examine the basic quality of the gap-filled NDSI collection considering cloud coverage and spatiotemporal consistency of the raw NDSI dataset by retrieving the bit flags from the integer stored in QA maps. The specific attributes are listed in Table 3.

Table 3Bit flags indicating the retrieval conditions according to the raw NDSI dataset.

Download Print Version | Download XLSX

The snow detection reversal of the pre-processed value in TAC is tracked in Bit 2, the post-processing (Sect. 2.2.4) is tracked in Bit 3, and the number of iterations primarily related to cloud coverage is indicated by Bit 4. If the range of reference values is larger than 30, then the Bit 5 flag is set; if the difference between pre-processed and fused values is larger than 20, the Bit 6 flag is set. Bit 7 reflects the cloud coverage of the space–time block. Furthermore, the combination of Bits 0 and 1 is a qualitative estimate of the cloud-removed NDSI collection based on the number of iterations (hereafter referred to as NI) and spatiotemporal consistency. The comprehensive estimate is determined as follows:

• if NI = 5 and Bit 5 or 6 is set to 1, then it is assigned “poor”;

• if NI = 4 and Bit 5 or 6 is set to 1, then it is assigned “OK”;

• if NI = 3 and Bit 5 or 6 is set to 1, then it is assigned “good”;

• otherwise, it is assigned “best”.

The QA maps are recommended for in-depth application of the cloud-removed NDSI collection.

2.3 Validation of the NDSI collection

The gap-filled NDSI collection is evaluated with the in situ snow-depth observations and Landsat NDSI maps considering classification and numerical accuracies according to the current mature verification methods. As shown in Table 4, the classification metrics based on the confusion matrix include overall accuracy (OA), commission error (CE), and omission error (OE) (Klein and Barnett, 2003). Moreover, 3 general metrics are introduced to measure numerical accuracy: correlation coefficient (CC), absolute error (AE), and root-mean-square error (RMSE).

Table 4Confusion matrix and validation metrics.

Download Print Version | Download XLSX

3.1 Validation against in situ snow depth measurements

As described above, the in situ snow-depth data in NX from 1 January 2001 to 31 August 2007 and in QTP from 1 August 2000 to 31 December 2013, were used as the ground truth to evaluate the classification accuracy of the TAC NDSI, NIEER AVHRR SCE, and STAR NDSI datasets. The nearest pixel was matched with each meteorological station, with a total of about 550 000 data pairs. Snow-covered pixels in NDSI datasets range from 10 to 100, whereas snow-free pixels are 0; thus, the classification threshold is set as 10 (Zhang et al., 2019). The discriminant threshold for in situ snow depth is set as 0 or 1 cm. In addition, the cloud-covered areas in the TAC NDSI dataset are considered to be snow-free.

Figure 4Monthly classification accuracy of TAC NDSI, NIEER AVHRR SCE, and STAR NDSI products in NX (group a) and QTP (group b). Note that the optimal values for OA, CE and OE are 100 %, 0 % and 0 %, respectively.

Download

Table 5Confusion matrices between TAC NDSI, NIEER AVHRR SCE, STAR NDSI datasets and in situ snow-depth (SD) data in NX from 1 January 2001 to 31 August 2007. Note that bold values indicate classification accuracy, including OE, CE, and OA.

Download Print Version | Download XLSX

Table 5 demonstrates that NIEER AVHRR SCE and STAR NDSI datasets preeminently capture the snow dynamics in NX referring to the in situ measurements, with OAs of more than 94 %. However, the TAC NDSI dataset is insufficient to accurately describe the snow cover variability. Although CEs perform well regardless of the snow-depth threshold, OEs of the TAC NDSI collection are extremely high, indicating that many cloud-covered areas are dominated by snow. The NIEER AVHRR SCE dataset partially retrieves snow pixels under cloud obstruction with an OE decreased by ∼43 %. The STAR NDSI collection completely removes clouds and accurately presents snow distribution, with an OE further decreased from ∼17 % to ∼7 %. The generation procedure in NX has two strengths. Firstly, the satellite-borne sensors can accurately capture the snow events on the ground due to the generally thick snow averaging approximately 20 cm. Secondly, the gap-filling approach with comprehensive consideration of spatial and temporal correlation has outstanding reliability due to the significant periodicity of snow variation. It can be inferred that the NDSI datasets in NC have high accuracy because of the similar snow conditions, despite the lack of in situ data in this region.

By contrast, despite the satisfactory performance of OAs and CEs, the OEs of three snow cover datasets over QTP are as remarkably high as 72 %, 42 %, and 39 %, even at the snow-depth threshold of 1 cm (Table 6). This finding indicates the omission of a large number of snow-covered pixels. The specific reasons are as follows: firstly, the original MODIS NDSI maps frequently underestimate the snow presence throughout the snow period; because discriminating the shallow snow pixels with an averaged snow depth of approximately 4 cm over QTP is challenging. Secondly, the credibility of the spatiotemporal contextual information is relatively low because the patchy snow rapidly and irregularly varies due to the extremely complex topographic and climatic conditions, leading to a further decrease in the accuracy of the gap-filled results. Lastly, the meteorological stations over QTP are unevenly distributed and are mostly located in low- and medium-altitude/latitude areas dominated by transient snow. Consequently, the evaluation results slightly exaggerate the real OEs.

For the in-depth analysis of the temporal characteristics, the monthly classification accuracies of TAC NDSI, NIEER AVHRR SCE and STAR NDSI products in NX and QTP are shown in Fig. 4 (the horizontal axis is the month in a hydrological year). In NX (group a), the monthly snow fraction in the in situ samples is greater than 85 % from December to next February. Therefore, the clouds in the TAC NDSI dataset seriously affect the snow cover estimation, while both cloud-free products exhibit superior OAs. Compared to the NIEER AVHRR SCE product, the STAR NDSI collection has slightly higher CEs but relatively lower OEs. The OEs of the STAR NDSI collection typically occur during snow accumulation and ablation periods, and almost disappear during stable snow-covered and snow-free periods. In QTP, the snow period is generally from October to next May, with the monthly snow fraction of less than 10 % in the in situ samples. Consequently, the underestimation of snow coverage caused by the clouds in the TAC NDSI dataset is slight. All three products achieve outstanding OAs and CEs, but exhibit relatively poor OEs. In NIEER AVHRR SCE and STAR NDSI datasets, these OEs are generally observed outside the snow period. As mentioned above, there are three reasons for this phenomenon. Nonetheless, the STAR NDSI collection presents superior classification accuracy to TAC NDSI and NIEER AVHRR SCE datasets.

Table 6Confusion matrices between TAC NDSI, NIEER AVHRR SCE, STAR NDSI datasets and in situ snow-depth (SD) data in QTP from 1 August 2000 to 31 December 2013. Note that bold values indicate classification accuracy, including OE, CE, and OA.

Download Print Version | Download XLSX

Overall, the STAR NDSI collection is capable of snow status estimation, eliminating cloud contamination in the TAC NDSI dataset, and capturing more snow events than the NIEER AVHRR SCE dataset. However, the accuracy of the STAR NDSI collection has regional and temporal heterogeneity. Firstly, the accuracy over QTP is lower than that of NX, which is consistent with the characteristic of the original MODIS NDSI maps. Furthermore, the permanent and periodic snow regime regions reconstructed by STAR have prominently high accuracy, while the transient snow-covered regions are easily omitted. Fortunately, the monitoring of permanent and periodic snow plays a key role in most snow-related investigations. Finally, the accuracy of stable snow-covered and snow-free periods is slightly better than that of snow accumulation and ablation periods.

3.2 Validation based on Landsat NDSI maps

Only the classification accuracy can be evaluated by in situ measurements due to the significant difference in the nature of the snow depth and NDSI data. Therefore, Landsat images with finer spatial resolution were commonly adopted for the numerical evaluation of NDSI datasets (Crawford, 2015). The NDSI values for the Landsat 8 dataset were calculated as follows: $\left(\text{Band3}-\text{Band6}\right)/\left(\text{Band3}+\text{Band6}\right)$. Subsequently, the average of the 17×17 neighborhoods closest to the center of the MODIS NDSI pixel in the Landsat NDSI map was considered to be the reference of this MODIS NDSI pixel to match the spatial resolution of Landsat with that of MODIS. Specifically, the cloud-contaminated pixels marked by the quality band in Landsat images were excluded, and the reference areas with cloud coverage larger than 30 % were discarded. A total of 19 Landsat NDSI maps with different snow coverages from January to April in 2018, which are distributed in NC (four scenes), the Central China region (CCR, two scenes), QTP (eight scenes), and NX (five scenes), were exploited as evaluation references for this validation experiment. Two evaluations including a cross-comparison of TAC NDSI, MODIS CGF NDSI, and STAR NDSI datasets, and an internal comparison of clear-sky and cloud-covered areas are described in detail below.

For the cross-comparison, the visual effects of three NDSI datasets on 8 January 2018 and 3 February 2018 are shown in Fig. 5. The TAC NDSI dataset is still heavily obscured by clouds. Although the MODIS CGF NDSI dataset completely removes clouds from the MOD10A1 product, it is difficult to accurately retrieve periodic and transient snow cover regions due to the simplicity of the cloud-gap-filled method. Specifically, the gaps are filled by retaining clear-sky observations from previous days. However, snow patterns under cloud cover are likely to change significantly during these days. Therefore, snow cover is significantly underestimated during accumulation (Fig. 5a2) and overestimated during ablation (Fig. 5b2). By contrast, the STAR NDSI collection pre-eminently captures the snow dynamics under temporally continuous clouds, attributing to the spatiotemporal adaptive fusion strategy. Furthermore, the three NDSI datasets are quantitatively assessed by Landsat NDSI maps.

Figure 5Comparison of TAC NDSI (column 1), MODIS CGF NDSI (column 2), and STAR NDSI (column 3) products on 8 January 2018 (group a) and 3 February 2018 (group b).

As the average cloud cover is as high as 40.7 %, the TAC NDSI dataset has a low correlation with Landsat NDSI maps, with an average CC of 0.46 (Table 7). The cloud-free MODIS CGF NDSI dataset enhances the correlation with Landsat NDSI maps, with an average CC of 0.73. By contrast, the snow dynamics presented by the spatiotemporally continuous STAR NDSI dataset are highly consistent with Landsat NDSI maps, with an average CC further increased to 0.84. Compared to TAC NDSI and MODIS CGF NDSI datasets, the average RMSE of STAR NDSI dataset is decreased by 9.06 and 5.58, respectively. The positive AEs reveal that NDSI values for snow pixels in MODIS CGF NDSI and STAR NDSI datasets are generally higher than those of Landsat NDSI maps. In terms of snow coverage, the STAR NDSI dataset notably improves the detection of snow events compared to the other two datasets, with an average absolute SRD decreased by 31.1 % and 2.4 % (SRD indicates the difference of snow rate between MODIS and Landsat NDSI datasets). Consequently, the STAR NDSI collection is a more promising snow cover product than TAC NDSI and MODIS CGF NDSI datasets, contributing to related hydrological and meteorological studies.

Table 7Performance statistics for TAC NDSI, MODIS CGF NDSI, and STAR NDSI datasets against Landsat NDSI maps.

Note that TAC, CGF, and STAR represent TAC NDSI, MODIS CGF NDSI, and STAR NDSI datasets, respectively.

Download Print Version | Download XLSX

In addition to the cross-comparison of TAC NDSI, MODIS CGF NDSI, and STAR NDSI datasets, an internal comparison of the STAR NDSI collection in clear-sky and cloud-covered areas was performed based on Landsat NDSI maps, to highlight the accuracy of the recovered pixels in the STAR NDSI collection. As described above, clear-sky and cloud-covered areas account for 59.3 % and 40.7 %, respectively. Table 8 indicates that the snow distribution of recovered areas in the STAR NDSI collection is relatively consistent with that of Landsat NDSI maps. Although the average CC decreases from 0.85 to 0.73 and the average RMSE increases from 13.48 to 16.30 compared to clear-sky areas, the accuracy of recovered areas is satisfactory. Since many recovered areas inherit errors from clear-sky areas because the cloud-removal procedure completely relies on the original dataset, a slight decrease in accuracy is reasonable. In addition, the average AEs of clear-sky and recovered areas are 7.81 and 6.83, respectively, revealing the systematic overestimation of NDSI values in snow cover regions (Landsat NDSI values are generally low). Except for a few areas, the snow conditions in most cloud-covered areas are well recovered, with an average SRD of −5.0 %. This finding highlights that the STAR NDSI collection can completely remove clouds with satisfactory accuracy.

Table 8Performance statistics for STAR NDSI collection against Landsat NDSI maps in clear-sky and cloud-covered areas according to the TAC dataset. Note that ** and * respectively indicate an improvement and degradation of cloud-covered areas compared to clear-sky areas (corresponding to four groups in Fig. 6).

Download Print Version | Download XLSX

For in-depth verification analysis, Fig. 6 shows the visual effects in four regions corresponding to four highlighted groups in Table 8. The accuracy of clear-sky areas in NC is prominently high with a CC of 0.92, while recovered areas notably reduce the performance with a CC of 0.42. However, Fig. 6a shows that clear-sky areas in the TAC NDSI dataset cannot reflect the snow events, whereas the STAR NDSI collection effectively retrieves these events. Inevitably, the snow edges are slightly inaccurate and blurry due to insufficient reference information. For group (b), the original accuracy of the NDSI dataset in CCR is relatively low, with high cloud coverage and a false acceptance rate, while the STAR NDSI collection presents a snow pattern consistent with Landsat NDSI. Nevertheless, CCR is a transient snow area with relatively low altitudes and latitudes. Therefore, the gap-filled result has visible uncertainty, in which commission (black box) and omission (red box) frequently occur. As mentioned above, MODIS NDSI datasets generally perform poorly over QTP. However, Fig. 6c demonstrates that TAC NDSI and STAR NDSI collections can accurately capture snow events despite a few omissions (red box). Similar to the NC region, CCs provide a positive indication of overall performance in NX. As shown in Fig. 6d, even the NX4_20180103 scene, with a remarkably low CC of 0.58, can effectively reflect the snow pattern. In addition, the NDSI dataset retrieved by STAR inevitably has a few extremely abnormal areas during 20 years due to the vast territory of China; an example is presented in Fig. A2 (Appendix A). These areas have severe cloud contamination and irregular snow dynamics, contributing to the challenges in reconstruction and evaluation. Therefore, these areas are corrected by post-processing as described in Sect. 2.2.4.

Figure 6Comparison of TAC NDSI (column 1), Landsat NDSI (column 2), and STAR NDSI (column 3) products, and classification consistency (column 4) corresponding to NC3_20180311, CCR2_20180203, QTP2_20180225, and NX4_20180103 (groups a to d).

Overall, the numerical accuracy of the STAR NDSI collection is rigorously evaluated based on fine-resolution Landsat NDSI maps. The cross-comparison indicates that the STAR NDSI collection is superior to both TAC NDSI and MODIS CGF NDSI datasets. In addition, the internal comparison reveals that the effectiveness of cloud removal is satisfactory, although recovered areas have slightly lower accuracy than clear-sky areas. Consequently, the STAR NDSI collection has considerable application potential.

3.3 Application of the STAR NDSI collection

In addition to quality evaluation, the exemplary analysis also contributes to understanding the potential of the STAR NDSI collection for hydrological and climatic applications. Therefore, the spatial distribution and temporal variability of snow cover in China are analyzed as two simple application demonstrations.

From a spatial perspective, the sequence shown in Fig. 7 indicates that the snow fraction first increases and then decreases in NC and NX regions from 5 December 2014 to 25 January 2015. However, due to the complex topographic and climatic conditions, the snow cover in QTP presents irregular distribution and considerable fluctuations. During this period, the snow cover pattern in China changed dramatically, with the overall snow fraction ranging from 19.25 % to 35.42 %. In addition, another sequence of cloud-free NDSI maps in early 2020 is shown in Fig. 8. During this period, the snow cover in China presents a single-wave depletion curve. The peaks of snow fractions are observed in three major snow areas between 11 January 2020 and 21 January 2020. Compared to the previous sequence, the snow cover variation in QTP in this sequence has significantly enhanced regularity.

Figure 7A sequence of the STAR NDSI collection from 5 December 2014 to 25 January 2015.

Figure 8A sequence of the STAR NDSI collection from 1 January 2020 to 21 February 2020.

Figure 9Daily average snow fraction of NX, NC, QTP, and China.

Download

Table 9The cloud-removal effectiveness of STAR compared to Landsat NDSI maps under different simulated cloud conditions. Note that bold values indicate a significant degradation of the accuracy under the current cloud cover compared to the previous one.

Download Print Version | Download XLSX

Figure 9 shows the daily average snow fraction in the three main subregions of China and the entire situation considering temporal analysis. In terms of intra-annual variability, the snow dynamics periodically evolve in NX and NC but substantially fluctuate in QTP. Similar snow depletion curves are present in NX and NC, demonstrating rapid accumulation and ablation in November and March, respectively. There is a relatively long snow period in QTP, with an average snow fraction varying from 20 % to 40 % from October to next May. Consequently, China is dominated by periodic snow. As for inter-annual variability, among the three major snow areas, the snow fraction in NC remarkably fluctuates with a standard deviation of 5.3 %. The snow coverage in QTP presented a slight decreasing trend from 2005 to 2017 but increased significantly in the past 2 years. In particular, rather than a significant rise in maximum snow coverage, the increase can be observed throughout the snow period with a slight reduction in intra-annual volatility. This finding implies that the regional climatic conditions tend to stabilize slightly. In addition, no significant trend has been detected in snow dynamics in China during the 20 years. Nevertheless, the significant fluctuation of maximum snow coverage in China indicates the presence of non-negligible large-scale transient snow cover regions. For example, Fig. A3 (Appendix A) shows the extreme snow event in southern China caused by the La Niña phenomenon, which resulted in heavy casualties and economic losses in the hydrological year 2007–2008.

Overall, the STAR NDSI collection can accurately reflect the spatial and temporal dynamics of snow cover in China. It is promising for hydrological and meteorological applications.