Articles | Volume 14, issue 9
Earth Syst. Sci. Data, 14, 4445–4462, 2022
https://doi.org/10.5194/essd-14-4445-2022

Special issue: Extreme environment datasets for the three poles

Earth Syst. Sci. Data, 14, 4445–4462, 2022
https://doi.org/10.5194/essd-14-4445-2022
Data description paper
29 Sep 2022
Data description paper | 29 Sep 2022

HMRFS–TP: long-term daily gap-free snow cover products over the Tibetan Plateau from 2002 to 2021 based on hidden Markov random field model

HMRFS–TP: long-term daily gap-free snow cover products over the Tibetan Plateau from 2002 to 2021 based on hidden Markov random field model
Yan Huang1,2,, Jiahui Xu1,2,, Jingyi Xu1,2, Yelei Zhao1,2, Bailang Yu1,2, Hongxing Liu3, Shujie Wang4, Wanjia Xu1,2, Jianping Wu1,2, and Zhaojun Zheng5 Yan Huang et al.
  • 1Key Laboratory of Geographic Information Science, Ministry of Education, East China Normal University, Shanghai 200241, China
  • 2School of Geographic Sciences, East China Normal University, Shanghai 200241, China
  • 3Department of Geography, the University of Alabama, Tuscaloosa, AL 35487, USA
  • 4Department of Geography, Earth and Environmental Systems Institute, Pennsylvania State University, University Park, PA 16802, USA
  • 5National Satellite Meteorological Center, Beijing 100081, China
  • These authors contributed equally to this work.

Correspondence: Yan Huang (yhuang@geo.ecnu.edu.cn)

Abstract

Snow cover plays an essential role in climate change and the hydrological cycle of the Tibetan Plateau. The widely used Moderate Resolution Imaging Spectroradiometer (MODIS) snow products have two major issues: massive data gaps due to frequent clouds and relatively low estimate accuracy of snow cover due to complex terrain in this region. Here we generate long-term daily gap-free snow cover products over the Tibetan Plateau at 500 m resolution by applying a hidden Markov random field (HMRF) technique to the original MODIS snow products over the past two decades. The data gaps of the original MODIS snow products were fully filled by optimally integrating spectral, spatiotemporal, and environmental information within HMRF framework. The snow cover estimate accuracy was greatly increased by incorporating the spatiotemporal variations of solar radiation due to surface topography and sun elevation angle as the environmental contextual information in HMRF-based snow cover estimation. We evaluated our snow products, and the accuracy is 98.29 % in comparison with in situ observations, and 91.36 % in comparison with high-resolution snow maps derived from Landsat images. Our evaluation also suggests that the incorporation of spatiotemporal solar radiation as the environmental contextual information in HMRF modeling, instead of the simple use of surface elevation as the environmental contextual information, results in the accuracy of the snow products increases by 2.71 % and the omission error decreases by 3.59 %. The accuracy of our snow products is especially improved during snow transitional period, and over complex terrains with high elevation and sunny slopes. The new products can provide long-term and spatiotemporally continuous information of snow cover distribution, which is critical for understanding the processes of snow accumulation and melting, analyzing its impact on climate change, and facilitating water resource management in Tibetan Plateau. This dataset can be freely accessed from the National Tibetan Plateau Data Center at https://doi.org/10.11888/Cryos.tpdc.272204 (Huang and Xu, 2022).

1 Introduction

Snow cover has the characteristics of low thermal diffusivity, high reflectivity, and strong water storage capacity, which has a profound effect on climate change (Gao et al., 2012), radiation budget (Yang et al., 2014; Huang et al., 2019), hydrological cycle (Dong, 2018), and human activities (Cereceda-Balic et al., 2020). The Tibetan Plateau (TP) has abundant snow cover, with highest elevation and largest snow cover area in the middle latitudes of the Northern Hemisphere (Qiu, 2008; Yao et al., 2019). Snow is highly sensitive to climate change (Chen et al., 2018a), and snowmelt water quantity is closely connected to the supply of soil moisture on the plateau (Wang et al., 2018) and the runoff of numerous rivers such as the Yangtze and Yellow rivers (Immerzeel et al., 2010). Long-term and detailed snow cover information is fundamental to investigating climate change and hydrological cycle of the TP.

Remote sensing has allowed extraction of historical and near-real-time snow cover extent over inaccessible areas (Yang et al., 2015; Li et al., 2018). Snow products could be acquired via remote sensing satellites, such as the Landsat and Sentinel-2 series. Although these satellites have a high spatial resolution, they have a relatively coarse (10 or 16 d) temporal resolution, thereby rendering them insufficient to monitor the temporal variations in snow cover (Huang et al., 2022). Moderate Resolution Imaging Spectroradiometer (MODIS) has commenced producing snow products at 500 m resolution from 2000, which have been widely utilized as the primary datasets for monitoring snow cover (Muhammad and Thapa, 2020; Kilpys et al., 2020). The accuracy of MODIS snow products is greater than 85 % at the global scale under clear sky (e.g., Parajka et al., 2012; Yang et al., 2015). However, these products have many data gaps due to frequent clouds, causing discontinuity in the time and space of the products (Y. Liu et al., 2020​​​​​​​). In addition, the complex terrain of the TP makes it more challenging for accurate snow cover detection. Previous studies (Dong and Menzel, 2016; Dariane et al., 2017) have suggested the low accuracy of MODIS snow products for mountainous areas.

Various data-gap-filling techniques have been proposed to produce seamless MODIS snow cover products including multi-source combination, temporal, spatial, and spatiotemporal filters (Xiao et al., 2021; Hussainzada et al., 2021; Richiardi et al., 2021). The multi-source combination method combines MODIS with passive microwave sensor data (Y. Li et al., 2019; Li et al., 2020). Although this method can be used to fill most data gaps, the accuracy of the combined product depends more on the accuracy of the microwave sensor data, and the spatial resolution (>1 km) is not sufficient to meet the demands for accurately assessing the snow cover variability (Wang et al., 2009). Many attempts have been directed toward filling the gaps in MODIS optical remote sensing data based on spatial and temporal information. Temporal methods are used to reclassify the data-gap pixels by inferring the land cover types of the current pixels under clear sky within a few days before or after (e.g., the previous or next day) (X. Li et al., 2019; Tran et al., 2019). Spatial methods estimate data-gap pixels based on gap-free pixels in the spatial neighborhood (Hou et al., 2019). Relevant environmental information (e.g., snow lines and topography) has been introduced into spatial methods (Wang et al., 2019; Kilpys et al., 2020). Spatiotemporal methods have been integrated to fill as many data gaps as possible (Parajka and Blöschl, 2008; Kilpys et al., 2020). These studies integrated all information (i.e., spatiotemporal, and environmental information) according to heuristic rules instead of rigorous quantitative models. Huang et al. (2018) developed a hidden Markov random field (HMRF) framework to optimally combine all information. This technique not only fills data gaps but also provides fine improvement in accuracy compared to the original MODIS snow products.

The long-term series and high-precision products are the basis for snow cover research on the TP, which enable the monitoring and analysis of snow cover phenology and more comprehensive understanding of the snow cover trend. However, the existing daily products for the TP have an earlier end date (the latest version at 500 m resolution product ended in 2015), and some products still have a small number of data gaps (Yu et al., 2016; Xu et al., 2017; Zheng and Cao, 2019). Thus, long-term and high-precision daily snow cover products for the TP should be generated using reliable data-gap-filling methods.

Here we generate long-term daily snow cover products over the TP by applying an HMRF technique (Huang et al., 2018) to the original MODIS snow products from 2002 to 2021. In the previous HMRF modeling technique (Huang et al., 2018), surface elevation was considered as a surrogate for environmental information in mountainous regions. However, according to in situ photos from field survey in the TP, we found that the distribution of snow cover differed from the model results, even at the same elevation level (Fig. 1b and c). The snow on the TP is strongly infected by complex topography (e.g., slope, aspect, sunlight duration, and solar incidence angle). Thus, while generating our new snow cover product, we incorporated solar radiation as a comprehensive indicator of topographical effects to correct the snow identification errors caused by the complex terrain of the TP.

Our study is outlined as follows: first, we present our study area and datasets. Then, we present the HMRF framework and data-processing flowchart. By employing this modeling technique to the TP, we produced daily gap-free snow cover products for the TP from 2002 to 2021. We also validated our new snow products against in situ observations, snow cover mapped from Landsat series data, and snow cover data estimated from the initial MODIS and original HMRF modeling technique.

2 Study area and data

2.1 Study area

The TP is situated in western China, spanning 2600–3946 N and 7318–10446 E (Fig. 1). It encompasses approximately 2.56 ×106 km2, with an average altitude exceeding 4000 m, and is one of the most susceptible regions to climate change (Jing et al., 2019). The air temperature in the TP increased from the northwest to southeast, with more precipitation in the southeast and less in the northwest (Wang et al., 2018). Generally, the TP has comparatively low temperatures and has largest distribution of glaciers and snow in China (Liang et al., 2017). Seasonal snow cover on the TP has a great potential to influence the hydrological cycle and heat wave frequency in northern China (Wu et al., 2012). In addition, seasonal snow accumulation on the TP is an important part of surface water accumulation in southwestern China and surrounding countries. Several major rivers in China and surrounding Asian countries, such as the Yangtze, Yellow, Mekong, Salween, Brahmaputra, Ganges, and Indus rivers, all originate from the TP.

The status of the snow in the TP changes rapidly as a result of multiple factors, such as temperature, precipitation, synoptic forcing, and large-scale ocean–atmosphere oscillations (You et al., 2020), which may lead to sublimation and snowfall (Li et al., 2020). Affected by the Indian Ocean monsoon and East Asian summer monsoon, the southeast TP has plenty water vapor, which results in a large amount of cloud, particularly in spring and summer (Yang et al., 2015).

https://essd.copernicus.org/articles/14/4445/2022/essd-14-4445-2022-f01

Figure 1Topography, meteorological stations, and survey photos of the TP. (a) Surface elevation and distribution of meteorological stations in the TP. Landsat images utilized for validation and sample area are also shown; (b) and (c) in situ photos from field survey in the TP.

2.2 Datasets

2.2.1 Daily MODIS snow cover products

We used version six MODIS daily snow products on Terra (MOD10A1 v6) and Aqua (MYD10A1 v6) from 15 May 2002 to 31 December 2021 (Hall and Riggs, 2016a, b). Google Earth Engine (GEE) (https://earthengine.google.com, last access: 20 January 2022) platform was used for pre-processing the MODIS snow product. We used MODIS/006/MOD10A1 and MODIS/006/MYD10A1 datasets and NDSI_Snow_Cover and NDSI_Snow_Cover_Class bands. NDSI_Snow_Cover represents the value of the normalized difference snow index (NDSI), ranging 0–100. The values in the NDSI_Snow_Cover_Class band include 200, 201, 211, 237, 239, 250, 254, and 255, which represent “missing data”, “no decision”, “night time”, “inland water”, “ocean”, “clouds”, “detector saturated”, and “filled data”, respectively (Riggs et al., 2019).

The original NDSI data were reclassified as snow, non-snow, and data-gaps classes (Chen et al., 2020; Huang et al., 2018, 2022). The pixels with an NDSI value of 40–100 in the NDSI_Snow_Cover band were reclassified as snow (Riggs et al., 2017). The pixels with an NDSI value of 0–40 in the NDSI_Snow_Cover band and with corresponding values of 237 and 239 in the NDSI_Snow_Cover_Class band were reclassified as non-snow. The remaining pixels were reclassified as data gaps. Then, we combined the MOD10A1 and MYD10A1 reclassification results of the same day according to the following rules: when a pixel included gap-free data in both MOD10A1 and MYD10A1 products, value in the MOD10A1 product was used; when a pixel included gap free only in the one product, the gap-free data value was used. Finally, the combined Terra and Aqua results were re-projected onto Universal Transverse Mercator (UTM zone 45) at 500 m resolution, which was used as the initial snow cover products for filling data gaps in this research.

We identified the data gaps in original composite MODIS products, and found the annual average data gap proportion of the TP was 33.63 %–38.75 %, with an average of 36.12 % (Fig. 2). The monthly average data gaps are also shown in Fig. 2, and we found that the data gaps of the TP were the largest in summer (June to August), with an average of 45.20 %. The data gaps declined rapidly starting from September and were the smallest in November (with an average of 28.00 %), and gradually increased in winter and spring.

https://essd.copernicus.org/articles/14/4445/2022/essd-14-4445-2022-f02

Figure 2Annual and monthly average proportion of data gaps of original composite MODIS snow products.

Download

2.2.2 Digital Elevation Models (DEMs)

We used the Shuttle Radar Topography Mission (SRTM) DEMs to calculate the topographical parameters (e.g., elevation, slope, and aspect). The SRTM DEM at 90 m resolution in GeoTIFF format was from the United States Geological Survey (USGS). We then pre-processed the DEM data, including mosaicking, reprojection, resampling, and clipping.

2.2.3 In situ observations

We used snow depths from 137 in situ stations in the TP (Fig. 1a) obtained from 1 October 2002 to 31 March 2021, which were provided by the National Meteorological Centre of China. Figure 1a shows all stations over the TP, which record coordinates, observation time, snow depth, and snow pressure. The surface elevation of these stations ranges from 1187 to 5310 m, with an average of 3240 m. Overall, 72 % stations are situated in the east (>95 E), and only 22 % are situated in high-altitude areas (>4000 m). A 3 cm threshold was utilized to evaluate the snow classes using in situ snow depth (Huang et al., 2022). The snow depth was reclassified as no-snow if it was less than 3 cm; if not, it was considered as snow.

2.2.4 Landsat images

The snow cover mapped from Landsat images were used to validate and evaluate our snow cover product. Using the GEE cloud platform, we selected 81 Landsat series images with less than 2 % cloud coverage under clear sky for verification, including Landsat-5 TM, Landsat-7 ETM+, and Landsat-8 OLI images. The detailed Landsat images information is shown in the supplementary material (Table A1). The selected images were distributed throughout the study area from 22 November 2002 to 16 December 2021. We used the C topographic correction model to correct for the terrain effect on Landsat series data (Teilet et al., 1982). Then, using the SNOWMAP method proposed by Hall et al. (1995), we generated a snow cover map using Landsat images for the TP. The SNOWMAP algorithm classified the pixels with NDSI values >0.4, green band >0.10, and shortwave infrared band >0.11 as snow (Huang et al., 2011). To improve the accuracy of resampling 30 to 500 m resolution, we adopted the following steps for resampling: first, for each MODIS pixel, we calculated the amount of Landsat snow-covered pixels in the current MODIS pixel. Second, we divided the snow-covered Landsat pixels by the sum of the Landsat pixels contained in each MODIS pixel (a 500 m MODIS pixel is close to 277 Landsat pixels) (Crawford, 2015; C. Liu et al., 2020). Finally, we reclassified pixels whose results in the previous step were less than 0.5 as no-snow; otherwise, they were reclassified as snow. Hence, we obtained resampled 500 m resolution Landsat images for further verification.

3 HMRF modeling technique

The HMRF framework optimally combines MODIS spectral, spatiotemporal, and environmental information to fill data gaps, thereby increasing snow estimate accuracy (Huang et al., 2018). This framework is expressed as a linear energy function in which the total energy is modeled as the combination of all information (Eq. 1). Thus, the HMRF framework requires to specify energy function for each information and to determine the optimal parameters that minimize the total energy.

(1) U T β n , x i , N sp , N tp , I ev = λ x i U x i β n , x i + λ st U st β n , N sp , N tp + λ ev U ev β n , I ev ,

where λxi, λst, and λev are the spectral, spatiotemporal, and environmental parameters, respectively; and Uxi, Ust, and Uev are the spectral, spatiotemporal, and environmental energy functions, respectively.

In the original HMRF modeling technique, Huang et al. (2018) used surface elevation to represent the environmental association in mountainous areas (hereafter referred to as HMRFele). However, because of the spatial heterogeneity of the TP, using only relative elevation cannot reflect the influence of complex terrain on snow cover (Fig. 1b and c). Therefore, the influence of additional topographic factors (e.g., slope, aspect, sunlight duration, and solar incidence angle) on snow cover must be considered. Based on the above reasons, we incorporated solar radiation as a surrogate for environmental information in HMRF framework (hereafter referred to as HMRFsolar). Further details are provided in Sect. 3.3.

The overall flowchart of the HMRFsolar modeling technique is shown in Fig. 3. We used daily composite MODIS snow cover products as the initial dataset. First, we calculated the optimal parameters of each information using the Ho-Kashyap and DFBETAS algorithms (Bormann et al., 2012), based on randomly divided training samples on the initial snow cover product. Second, we determined the snow and non-snow classes of each pixel in the initial snow cover products using the optimal parameters and HMRFsolar algorithm. In the first round of the HMRFsolar algorithm, we employed a 3×3×3 spatiotemporal cubic neighborhood to model spatiotemporal energy. If data gaps remained, we further expanded the cubic neighborhood of these pixels. Finally, we validated and evaluated our HMRFsolar-based products against the in situ observations, snow mapped from Landsat images, and snow cover data from the initial MODIS and original HMRFele modeling technique.

https://essd.copernicus.org/articles/14/4445/2022/essd-14-4445-2022-f03

Figure 3Overall flowchart of the HMRFsolar-based framework. (SCP stands for snow cover products.)

Download

3.1 Spectral information

The spectral energy is considered as the possibility of the pixel pertaining to non-snow or snow according to its spectral information. Fractional snow cover (FSC) represents the proportion of snow in a pixel, which can be calculated using the NDSI value provided by the MODIS snow product (Salomonson and Appel, 2004). The general linear relationship between FSC and NDSI was derived over other regions, which has limited accuracy in the TP. Therefore, we re-fitted that empirical relationship of Terra and Aqua satellites in the TP using Landsat data over 20 years with 972 884 and 952 221 sample points, respectively. The spectral probability P(xi | β1) of snow class was re-fitted as follows:

(2)Pxi|β1Terra=(1.222×NDSI+0.038)/100(3)Pxi|β1Aqua=(1.164×NDSI+0.058)/100.

The correlation coefficients of the re-fitted equation of Terra and Aqua satellites were 0.86 and 0.89, respectively. The spectral energy Uxi(βn,xi) can be calculated as follows:

(4) U x i β n , x i - P x i | β 1 , if n = 1 - [ 1 - P x i | β 1 , if n = 2 ,

where if pixel i is classified as no-snow class, then n=2; if pixel i is the snow class, n=1.

3.2 Spatiotemporal information

To fully take advantage of spatiotemporal information, we combined the temporal and spatial information to form a spatiotemporal cubic neighborhood Nst (Huang et al., 2018). In modeling, initial snow products are updated iteratively using classification results from the previous iteration. Convergence ends when the proportion of pixels whose type changes between two subsequent iterations is less than 0.1 %, in which case the data gaps need to be calculated. If the data gaps still exist, the cubic neighborhood is temporally and spatially expanded in the next round.

We primarily set a 3×3×3 cubic neighborhood; that is, the current pixel and its neighboring pixels formed a cube on day t, the day before (t−1), and the day after (t+1). In this condition, the energy equation of the current pixel is expressed as the weight of the proportion of snow and non-snow in the spatiotemporal neighborhood of the current pixel. The weight is reversely proportional to the distance from the neighboring pixel to the center pixel:

(5) U st β n , N sp , N tp = - x = - 1 1 y = - 1 1 t = - 1 1 S ( x , y , t ) Dist ( x , y , t ) V ( x , y , t ) Dist ( x , y , t ) , if n = 1 - x = - 1 1 y = - 1 1 t = - 1 1 NS ( x , y , t ) Dist ( x , y , t ) V ( x , y , t ) Dist ( x , y , t ) , if n = 2 ,

where (x,y,t) is the relative coordinates of the pixel in the spatiotemporal cube; S(x,y,t), NS(x,y,t), V(x,y,t), and Dist(x,y,t) are defined respectively as follows:

(6)Sx,y,t=1,ifx,y,tpertains to snow0,ifx,y,tpertains to non-snow(7)NSx,y,t=1,ifx,y,tpertains to non-snow0,ifx,y,tpertains to snow(8)Vx,y,t=1,ifx,y,tpertains to valid class0,ifx,y,tpertains to data gaps(9)Distx,y,t=x2+y2+wt2,

where Distx,y,t is the distance from each neighborhood pixel to the center pixel, which determines the weight of each area pixel; w denotes the weight of the temporal distance related to the spatial distance. This value was calculated via semi-variogram analysis and was determined as 3 in this study.

If data gaps still remain in 3×3×3 spatiotemporal cubic neighborhood, we further expanded spatiotemporal neighborhood to 5×5×5. In this study, when we used a 5×5×5 cubic neighborhood, the data gaps were reduced to 0.0173 %, which had little effect on snow cover estimate accuracy. Hence, the temporal neighborhood remained at 5, whereas the spatial neighborhood continued to expand. Finally, all data gaps could all be filled by continuously expanding the spatial neighborhood, which remarkably improved the MODIS estimate accuracy.

3.3 Environmental contextual information

The geographic location and seasons can determine the total solar radiation received on the TP, and the complex topography of the TP determines the availability of solar radiation at specific locations in the region, which further affects accumulation and melting of snow cover and determines its distribution (You et al., 2020). Here we use solar radiation to include the influence of environmental factors on snow cover, such as slope, aspect, sunlight duration, and solar incidence angle.

Kumar et al. (1997) originally proposed the solar radiation model, SHORTWAVE, to calculate the direct, reflected, and diffuse solar energy received by the ground. Huang et al. (2015) used a graphics processing unit (GPU) methodology to improve the speed and effectiveness for modeling. In our study, we used the GPU-accelerated SHORTWAVE model as the environmental parameter in the HMRFsolar model to calculate the daily solar radiation of the TP.

A general-purpose desktop computer was used to test the parallel computational efficiency. The computer has an Intel Core™ i7-10700 CPU (16 cores and max clock rate is at 2.90 GHz), an NVIDIA GeForce RTX 2070 SUPER card with 2560 cores and 16 240 MB global memory, and Windows 11 Ultimate 64 bit operation system. The GPU-accelerated SHORTWAVE model uses latitude, slope, aspect, date, and interval time as inputs. To calculate the overall solar radiation in a day, we subdivided the day into subsequent time intervals, and set 15 min as the suitable interval time in this paper (Kumar et al., 1997; Antonic, 1998). The overall solar radiation Ia at time t can be calculated using Eq. (10),

(10) I a = I p + I d + I r ,

where Ip, Id, and Ir are the direct, diffuse, and reflected solar radiation, respectively, and can be estimated using the following equations:

(11)Ip=I0τbcosθ(12)Id=I0τdcos2β2sinα(13)Ir=rI0τrsin2β2sinα,

where I0 is the extra-atmospheric solar flux; τb denotes the atmospheric transmittance; τd denotes the diffuse skylight transmittance; τr is the reflectance transmissivity (Gul et al., 1998); r is set to 0.2 in the current paper (Kumar et al., 1997); α is the solar elevation angle, which varies in terms of the latitude, season, and time; and θ is the solar incidence angle, which is determined by solar elevation angle α, the solar azimuth angle, surface slope (β), and aspect (Huang et al., 2015). The effect of the geographic location and seasons on the total solar radiation are represented by solar elevation angle variable α, and the effect of the complex topography on the total solar radiation is represented by surface slope (β) and aspect variables.

We repeated the calculation for each time interval. Last, we accumulated the solar radiation from sunrise to sunset to acquire the daily values.

Areas that receive more solar radiation generally experience later snow accumulation and earlier snowmelt than areas that receive less solar radiation. Therefore, if pixels in a spatial neighborhood have higher solar radiation and the classification result is snow, the probability that the center pixel is also snow increases. The environmental energy equation is defined as follows:

(14) U ev β n , N tp , N sp = - x = - 1 1 y = - 1 1 SE ( x , y ) SE x , y + NSE ( x , y ) , if n = 1 - x = - 1 1 y = - 1 1 NSE ( x , y ) SE x , y + NSE ( x , y ) , if n = 2 ,

where (x,y) are coordinates of the pixel in the spatial neighborhood Nsp; and SE(xy) and NSE(x,y) are defined as follows:

(15)SEx,y=1,ifx,ybelongs to snowand solarx,ysolar(i)0,others(16)NSEx,y=1,ifx,ybelongs to non-snowand solarx,ysolar(i)0,others,

where solar(x,y) and solar(i) represent the solar radiation received by the neighboring pixels and the central pixel, respectively.

4 Results

4.1 Accuracy assessment based on in situ observations

Using HMRFsolar-based method, we produced gap-free snow cover products over the TP with a daily 500 m resolution from 15 May 2002 to 31 December 2021. We applied the validation indices based on the confusion matrix (Table 1), which is widely recommended by previous studies (Li et al., 2020).

https://essd.copernicus.org/articles/14/4445/2022/essd-14-4445-2022-f04

Figure 4Temporal variations in OA (overall accuracy) (a), OE (omission error) (b), and CE (commission error) (c) of HMRFele- and HMRFsolar-based snow products from 2002–2021.

Download

Table 1Confusion matrix comparing snow cover products with observation value.

Download Print Version | Download XLSX

The three validation indices are defined as follows:

(17)OA=a+da+b+c+d×100%(18)OE=ba+b×100%(19)CE=cc+d×100%.

Overall accuracy (OA) is the probability a pixel that is correctly classified; and omission error (OE) and commission error (CE) represent the underestimated and overestimated snow pixels, respectively. We first compared our daily HMRFsolar-based and original MODIS snow products with snow depth observed from 137 in situ stations from 2002–2021 (Table 2). The OA of all gap-free pixels in original MODIS products was 97.96 %. After employing the HMRFsolar-based method, the OA of these gap-free pixels was increased to 98.29 %. The OE was reduced by 0.57 %. Since the in situ stations in the TP are primarily distributed in low-altitude areas (Fig. 1a), only using in situ observations may not be representative and lead to biases in the accuracy assessment. Therefore, we also conducted evaluation in comparison with snow cover mapped from Landsat data obtained from 2002–2021, and regarded it as the true verification value in the following accuracy evaluation.

Table 2Confusion matrices between HMRFsolar-based, original MODIS snow products and in situ observation for gap-free pixels in original MODIS snow products during 2002–2021.

Download Print Version | Download XLSX

4.2 Accuracy assessment based on Landsat data

A total of 81 Landsat-8 images (Fig. 1a, Table A1) with less than 2 % cloud coverage were selected for this accuracy assessment. We constructed the confusion matrices for the HMRFsolar-based snow cover products, HMRFele-based snow products, and original MODIS snow cover products against snow cover mapped from Landsat series data products (Table 3). The OA for gap-free pixels in original MODIS products was 89.31 %. When we used surface elevation to represent the environmental association in HMRF modeling technique, the OA of these gap-free pixels was increased to 90.47 %. Since using surface elevation only cannot well-represent the influence of complex topography on snow cover distribution over the TP, this improvement is still limited. After incorporating the solar radiation to model the comprehensive influence of multiple topographic factors on snow cover, the OA of our HMRFsolar-based snow products was increased by 2.06 % compared with the OA of original MODIS products. C. Liu et al. (2020) indicated that the classification errors in original MODIS snow products were strongly affected by OE. Our snow cover product provided a considerable improvement in this respect. The OE of our HMRFsolar-based snow products was 3.24 % lower than that of original MODIS products (Table 3).

Particularly, our HMRF method is able to fill data gaps in original MODIS products in which spectral information is unavailable caused by cloud, by integrating spatiotemporal and environmental information together. Table 4 shows that the OA of the HMRFsolar- and HMRFele-based snow cover products for those data-gap pixels in original MODIS products was 83.86 % and 81.15 %, which is 7.50 % and 9.32 % lower than that of gap-free pixels. It is also showed that the OA of HMRFsolar-based snow cover product is 2.71 % higher than that of HMRFele-based snow cover product. Both Tables 3 and 4 confirmed that using solar radiation as the environmental information in the HMRF framework can develop more reliable snow cover products.

https://essd.copernicus.org/articles/14/4445/2022/essd-14-4445-2022-f05

Figure 5Effect of elevation on OA (a), OE (b), and CE (c) of HMRFele- and HMRFsolar-based snow products from 2002–2021.

Download

4.3 Accuracy improvement in different times of a snow season

To explore the details of snow cover variation, we define the snow season from September of the previous year to August of the following year, for example, the time range of the 2002 snow season is 1 September 2002–31 August 2003​​​​​​​ (Chen et al., 2018b). Previous studies (Huang et al., 2018; Klein, 2003; X. Li et al., 2019) have suggested that the snow cover estimation accuracy was generally low during snow transitional period (i.e., the beginning and end of a snow season). We plotted the temporal variations in the OA, OE, and CE of the HMRFsolar- and HMRFele-based snow cover products (Fig. 4). Except November, December, and February, the OAs of the new snow products were more than 90 % in all months. The accuracy was relatively low during snow transitional period (November, December, and February to April), whereas the accuracy was higher in snow stable period. The trend of the CE was almost the same as that of snow cover; that is, the months with more snow cover had a larger CE. However, the CE was generally low in all months. By comparison with the HMRFele-based products, the accuracy of our new products was higher in almost all months. In addition to the stable snow period, the improvement for transition period was also notable, which increased by an average of 2.18 %. In transition period, with rapid temperature changes, the snow status changed rapidly, and the snow mapped by MOD10A1 (crossing at 10:30) and MYD10A1 (crossing at 13:30) were different because of the different temperatures and solar radiation. The original MODIS snow cover products have relatively high error during snow transitional period (Y. Li et al., 2019​​​​​​​). After incorporating solar radiation as an environmental contextual information in HMRF, the accuracy of snow cover products has been improved effectively.

Table 3Confusion matrices between HMRFsolar-based snow products, HMRFele-based snow products, original MODIS snow products, and snow cover mapped from Landsat series data products for gap-free pixels during 2002–2021.

Download Print Version | Download XLSX

Table 4Confusion matrices between HMRFsolar-based snow cover products, HMRFele-based snow cover products, and snow cover mapped from Landsat series data products for data-gap pixels during 2002–2021.

Download Print Version | Download XLSX

4.4 Accuracy improvement for different surface topography

We divided the elevation of the study area into four categories: <3000, 3000–4000, 4001–5000, and >5000 m, to explore the effect of elevation on snow cover. We then calculated the snow cover products accuracy in each category. Figure 5 indicates that the OA of the products decreased with increasing elevation. The accuracy of our products was highest in the <3000 m category (94.26 %). With increasing elevation, the OE first decreased and then increased; that is, the OE value was the lowest in the 4001–5000 m category, whereas the CE was the highest in the 4001–5000 m category. Li et al. (2020) established the effect of elevation on snow cover and also demonstrated the accuracy of original MODIS snow products generally decreased with increasing elevation. Areas at higher elevation generally receive more solar radiation; hence, the snow status changes more rapidly at higher elevations, resulting in decreased snow cover estimation accuracy.

By comparison with the HMRFele-based products, the accuracy of our new products was higher in almost all elevation categories, and as the elevation increased, the accuracy improvement was more remarkable (reaching an improvement of 3.16 % in the >5000 m category). The improvements in OE first decreased and then increased; that is, the accuracy improvement of the OE value was the smallest in the 4001–5000 m category, with an improvement of 3.47 %. The CE exhibited the opposite trend; that is, the accuracy improvement in the CE value was the largest in the 4001–5000 m category, at 2.33 %. The accuracy improvement provided by our new snow cover product was substantially increased with the increase in elevation, indicating its effectiveness in high-altitude areas.

To explore the effect of slope on snow cover, we divided the slope of the study area into five categories: <10, 10–20, 21–30, 31–40, and >40. The OA of the snow cover products decreased with increasing slope, which is consistent with the trend of the changes in accuracy with increasing elevation (Fig. 6). The accuracy of our new products was highest in the 0–10 category, with a value of 93.78 %. The OE was high for relatively low-slope areas (<30), and the CE was high for high-slope areas (≥30).

By comparison with the HMRFele-based products, the accuracy of our new products was higher in almost all slope categories, and the accuracy improvement was remarkable in the 10–20 category (reaching an improvement of 2.90 %). This means topography of 10–20 slopes had the strongest impact on snow cover on the TP. The accuracy improvements in OE decreased with increasing slope, whereas the accuracy improvements in CE increased with increasing slope. In the <10 category, the improvement in CE was the smallest, at 0.56 %, whereas the improvement in OE was the largest, at 6.83 %. In the >40 category, the improvement in CE was the largest, at 2.20 %, whereas the improvement in OE was the smallest, at 0.81 %.

https://essd.copernicus.org/articles/14/4445/2022/essd-14-4445-2022-f06

Figure 6Effect of slope on OA (a), OE (b), and CE (c) of HMRFele- and HMRFsolar-based snow products from 2002–2021.

Download

Finally, we explored the effect of aspect on snow cover (Fig. 7). The OA of the snow cover products was higher on shaded slopes than on sunny slopes. The OAs of our new products on sunny and shaded slopes were 88.37 % and 90.07 %, respectively. By comparison with the HMRFele-based products, the accuracy improvement provided by our new products was remarkable on sunny slopes (improvement of 3.16 %) than on shaded slopes (improvement of 1.75 %), corresponding with the actual situation. The improvement in the OE on sunny slopes was also remarkable (improvement of 3.63 %).

https://essd.copernicus.org/articles/14/4445/2022/essd-14-4445-2022-f07

Figure 7Effect of aspect on OA (a), OE (b), and CE (c) of HMRFele- and HMRFsolar-based snow products from 2002–2021.

Download

5 Discussion

5.1 Effect of the HMRFsolar-based snow cover products over complex terrain

Long-term and high-precision snow cover products are key to the snow hydrology research in the TP. Compared to the heterogeneous land cover, complex terrain has more severe effects on MODIS snow products in the TP (C. Liu et al., 2020; Azizi and Akhtar, 2021). To enhance snow detection accuracy in mountainous areas, we used solar radiation as the environmental contextual information rather than surface elevation used in the original HMRFele modeling technique (Huang et al., 2018). The calculated optimal parameters used in Eq. (1) could reflect the importance of each energy to a certain extent. The optimal parameters of the HMRFsolar-based products for spectral, spatial–temporal, and environmental information were 0.117, 1.294, and 0.532, respectively, while that of the original HMRFele-based products were 0.338, 1.419, and 0.576, respectively. The proportion of environmental information of HMRFsolar-based products (27.37 %) were higher than that of HMRFele-based products (24.70 %), which shows the improvement of HMRFsolar-based products over complex terrain. Moreover, Fig. 8a shows a true-color Sentinel-2B image. Figure 8b, c, and d show the original MODIS, HMRFele- and HMRFsolar-based snow cover products, respectively, on 31 October 2018. The examples in Fig. 8 (where the topography has a dominant effect on snow cover) show that the original HMRFele-based snow cover products mapped all data gaps in the shaded and sunny slopes at the same elevation as snow-covered pixels; however, after incorporating solar radiation, our new products accurately identified them as non-snow-covered pixels. In general, since sunny slopes receive relatively more solar radiation, the area of snow cover on sunny slopes is smaller than that on shaded slopes. The new HMRFsolar-based snow products more effectively filled the data gaps on sunny slopes, producing more realistic results. Solar radiation was particularly effective as an environmental contextual factor to correct the snow detection errors in mountainous areas.

https://essd.copernicus.org/articles/14/4445/2022/essd-14-4445-2022-f08

Figure 8Comparison between true-color Sentinel-2B imagery (a), original MODIS snow products (b), HMRFele-based snow products (c), and HMRFsolar-based snow products.

5.2 Potential applications of the HMRFsolar-based snow cover products

Compared with in situ observations, the overall accuracy of snow cover products on the TP reported by other studies is in the range of 90.74 %–96.6 %, and the OE is greater than the CE (Yu et al., 2016; Xu et al., 2017; Zheng and Cao, 2019). The overall accuracy of our new snow products is 98.29 % in comparison with in situ observations, and the new product has a considerable improvement in OE. The new snow cover products provide spatiotemporally continuous snow cover distribution at 500 m resolution over the past two decades, thus having a great application potential for analyzing the processes of snow accumulation and melting. Figure 9 shows the snow cover trends over a sample area (location shown in Fig. 1) from our daily gap-free snow products and 8 d composite MODIS snow products on Terra (MOD10A2 v6.1) (Hall and Riggs, 2021​​​​​​​) during the 2017 snow season. We chose this sample region because it covers more mountainous areas including Hengduan Mountain and Bayankala Mountain, and has seasonal snow for multiple years. The MOD10A2 v6.1 products composite 8 d MODIS data to reduce the data gaps in the original daily MODIS snow products, and represent the maximum snow cover extent in 8 d, which means once snow is identified on any day in the 8 d time window, the pixel is mapped as snow. Although the overall variation trends of snow cover are similar between these two products, the snow cover percentage obtained from the MOD10A2 v6.1 data was 19.60 % higher than that of daily products during the 2017 snow season. The status of the snow in the TP changes rapidly. Compared with the 8 d composite MODIS snow products (MOD10A2 v6.1), the HMRFsolar-based snow cover products can get more detailed and accurate information on snow accumulation and melting processes, particularly during the snow transitional periods (February to April and November to December) when short-term snowfall events occur frequently (Fig. 9). In addition, owing to its high temporal resolution (daily) and long-term spatiotemporal continuity, the generated snow cover products also provide important baseline information for monitoring climate change, calibrating hydrological models, and simulating snowmelt runoff in the TP.

https://essd.copernicus.org/articles/14/4445/2022/essd-14-4445-2022-f09

Figure 9Snow cover percentage for daily HMRFsolar-based snow cover products and 8 d composite MODIS snow products in a sample area during the 2017 snow season of the TP. (SCP stands for snow cover products.)

Download

5.3 Limitations of the HMRFsolar-based snow cover products

The accuracy evaluation for the new snow cover products had some limitations. The majority of in situ observations in the TP are distributed in the low-altitude areas in the east, with only 22 % located in high-altitude areas (Fig. 1a). The number of in situ observations in the snow category was considerably smaller than that in the non-snow category, which created challenges in verifying snow category classification (Zhang et al., 2021). C. Liu et al. (2020) demonstrated that compared to in situ observations, snow cover mapped from Landsat series images is more related to MODIS snow cover products owing to the closer spatial match. Therefore, the conclusions about the new snow products were based on the accuracy evaluation against Landsat images, which may still contain some bias. Additionally, the solar radiation indicated used only imitates the solar radiation under clear-sky conditions. In a follow-up study, we will simulate the solar radiation throughout the year according to the weather conditions (e.g., the effect of cloud cover), and use parallel calculations to enhance calculation efficiency and further improve the snow cover product accuracy.

6 Data availability

The long-term daily snow cover products produced here are gap-free at 500 m resolution under the Universal Transverse Mercator (UTM zone 45) projection, and can be freely accessed from the National Tibetan Plateau Data Center at https://doi.org/10.11888/Cryos.tpdc.272204 (Huang and Xu, 2022), which is stored as a zip file ( 1.61 GB) for each year. By uncompressing the zip file, the daily snow cover data are provided in GeoTIFF format, and the values in the snow cover products are classified as snow (1) and non-snow (2). The name of each file is “HMRFSTP_yyyyddd.tif”, in which HMRFSTP is the abbreviation of “hidden Markov random field-based snow cover products for Tibetan Plateau”, yyyy stands for year, and DDD stands for Julian day, such as 2002135.tif describes the snow cover on Tibet Plateau on the 135th day of 2002.

7 Conclusions

In this study, we generated long-term daily gap-free snow cover products at 500 m resolution from original MODIS snow products in the TP over the past two decades. The snow cover estimate accuracy was greatly improved by incorporating solar radiation as a surrogate for environmental contextual information in HMRF framework in mountainous areas. We validated and evaluated our snow cover products through comparison with in situ observations and high-resolution snow cover mapped from Landsat images, with accuracy estimates of 98.29 % and 91.36 %, respectively. Our evaluation also suggests that incorporating solar radiation, instead of the simple use of surface elevation as the environmental contextual information in HMRF framework, results in the accuracy of the snow products increases by 2.71 % and the OE decreases by 3.59 %. Specifically, the accuracy of the new snow products is particularly improved during snow transitional period and over complex terrains with high elevation and sunny slopes.

We believe that the long-term and spatiotemporally continuous snow cover products generated in this study have great potential to analyze the processes of snow accumulation and melting, to monitor climate change, and to understand hydrological cycling in the TP.

Appendix A

Table A1Landsat series images used for assessment of the HMRF-based snow cover products in this study.

Download XLSX

Author contributions

YH​​​​​​​ and HL developed the HMRF methodology and designed the study. JiaX​​​​​​​ and YZ applied the methodology and performed the validation. JingX and WX performed the visualization. BY, SW, and JW analyzed the results. ZZ worked on data curation. Both the authors devoted to the writing of the paper and data quality control.​​​​​​​

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

The views and interpretations in this publication are those of the authors and are not necessarily attributable to the ICI-MOD.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Extreme environment datasets for the three poles”. It is not associated with a conference.

Acknowledgements

We would like to thank the National Meteorological Centre of China for providing the in situ observations over the Tibetan Plateau.

Financial support

This research has been supported by the National Natural Science Foundation of China (grant nos. 42125604, 42071306, and 41701502).

Review statement

This paper was edited by Tao Che and reviewed by two anonymous referees.

References

Antonic, O.: Modelling daily topographic solar radiation without site-specific hourly radiation data, Ecol. Model., 113, 31–40, https://doi.org/10.1016/S0304-3800(98)00132-X, 1998. 

Azizi, A. H. and Akhtar, F.: Analysis of spatiotemporal variation in the snow cover in Western Hindukush-Himalaya region, Geocarto Int., 1–23​​​​​​​, https://doi.org/10.1080/10106049.2021.1939442, 2021. 

Bormann, K. J., McCabe, M. F., and Evans, J. P.: Satellite based observations for seasonal snow cover detection and characterisation in Australia, Remote Sens. Environ., 123, 57–71, https://doi.org/10.1016/j.rse.2012.03.003, 2012. 

Cereceda-Balic, F., Vidal, V., Ruggeri, M. F., and Gonzalez, H. E.: Black carbon pollution in snow and its impact on albedo near the Chilean stations on the Antarctic peninsula: First results, Sci. Total Environ., 743, 140801, https://doi.org/10.1016/j.scitotenv.2020.140801, 2020. 

Chen, S., Wang, X., Guo, H., Xie, P., Wang, J., and Hao, X.: A conditional probability interpolation method based on a space-time cube for MODIS snow cover products gap filling, Remote Sensing, 12, 3577, https://doi.org/10.3390/rs12213577, 2020. 

Chen, X., Long, D., Liang, S., He, L., Zeng, C., Hao, X., and Hong, Y.: Developing a composite daily snow cover extent record over the Tibetan Plateau from 1981 to 2016 using multisource data, Remote Sens. Environ., 215, 284–299, https://doi.org/10.1016/j.rse.2018.06.021, 2018a. 

Chen, X., Long, D., Hong, Y., Hao, X., and Hou, A.: Climatology of snow phenology over the Tibetan plateau for the period 2001–2014 using multisource data, Int. J. Climatol., 38, 2718–2729, https://doi.org/10.1002/joc.5455, 2018b. 

Crawford, C. J.: MODIS Terra Collection 6 fractional snow cover validation in mountainous terrain during spring snowmelt using Landsat TM and ETM+, Hydrol. Process., 29, 128–138, https://doi.org/10.1002/hyp.10134, 2015. 

Dariane, A. B., Khoramian, A., and Santi, E.: Investigating spatiotemporal snow cover variability via cloud-free MODIS snow cover product in Central Alborz Region, Remote Sens. Environ., 202, 152–165, https://doi.org/10.1016/j.rse.2017.05.042, 2017. 

Dong, C.: Remote sensing, hydrological modeling and in situ observations in snow cover research: A review, J. Hydrol., 561, 573–583, https://doi.org/10.1016/j.jhydrol.2018.04.027, 2018. 

Dong, C. and Menzel, L.: Producing cloud-free MODIS snow cover products with conditional probability interpolation and meteorological data, Remote Sens. Environ., 186, 439–451, https://doi.org/10.1016/j.rse.2016.09.019, 2016. 

Gao, J., Williams, M. W., Fu, X., Wang, G., and Gong, T.: Spatiotemporal distribution of snow in eastern Tibet and the response to climate change, Remote Sens. Environ., 121, 1–9, https://doi.org/10.1016/j.rse.2012.01.006, 2012. 

Gul, M. S., Muneer, T., and Kambezidis, H. D.: Models for obtaining solar radiation from other meteorological data, Sol. Energy, 64, 99–108, https://doi.org/10.1016/S0038-092x(98)00048-6, 1998. 

Hall, D. K. and Riggs, G. A.: MODIS/Terra Snow Cover Daily L3 Global 500m SIN Grid, Version 6, NASA [data set], https://doi.org/10.5067/MODIS/MOD10A1.006, 2016a. 

Hall, D. K. and Riggs, G. A.: MODIS/Aqua Snow Cover Daily L3 Global 500m SIN Grid, Version 6, NASA [data set], https://doi.org/10.5067/MODIS/MYD10A1.006, 2016b. 

Hall, D. K. and Riggs, G. A.: MODIS/Terra Snow Cover 8-Day L3 Global 500m SIN Grid, Version 61, NASA [data set], https://doi.org/10.5067/MODIS/MOD10A2.061, 2021. 

Hall, D. K., Riggs, G., and Salomonson, V. V.: Development of methods for mapping global snow cover using Moderate Resolution Imaging Spectroradiometer (MODIS) data, Remote Sens. Environ., 54, 127–140, https://doi.org/10.1016/0034-4257(95)00137-P, 1995. 

Hou, J., Huang, C., Zhang, Y., Guo, J., and Gu, J.: Gap-filling of MODIS fractional snow cover products via non-local spatio-temporal filtering based on machine learning techniques, Remote Sensing, 11, 90, https://doi.org/10.3390/rs11010090, 2019. 

Huang, G., Li, Z., Li, X., Liang, S., Yang, K., Wang, D., and Zhang, Y.: Estimating surface solar irradiance from satellites: Past, present, and future perspectives, Remote Sens. Environ., 233, 111371, https://doi.org/10.1016/j.rse.2019.111371, 2019. 

Huang, K., Zhang, Y., Tagesson, T., Brandt, M., Wang, L., Chen, N., Zu, J., Jin, H., Cai, Z., Tong, X., Cong, N., and Fensholt, R.: The confounding effect of snow cover on assessing spring phenology from space: A new look at trends on the Tibetan Plateau, Sci. Total Environ., 756, 144011, https://doi.org/10.1016/j.scitotenv.2020.144011, 2021. 

Huang, X., Liang, T., Zhang, X., and Guo, Z.: Validation of MODIS snow cover products using Landsat and ground measurements during the 2001–2005 snow seasons over northern Xinjiang, China, Int. J. Remote Sens., 32, 133–152, https://doi.org/10.1080/01431160903439924, 2011. 

Huang, Y. and Xu, J.: Daily cloud-free snow cover products for Tibetan Plateau from 2002 to 2021, National Tibetan Plateau Data Center [data set], https://doi.org/10.11888/Cryos.tpdc.272204, 2022. 

Huang, Y., Chen, Z., Wu, B., Chen, L., Mao, W., Zhao, F., Wu, J., Wu, J., and Yu, B.: Estimating roof solar energy potential in the downtown area using a GPU-accelerated solar radiation model and airborne LiDAR data, Remote Sensing, 7, 17212–17233, https://doi.org/10.3390/rs71215877, 2015. 

Huang, Y., Liu, H., Yu, B., Wu, J., Kang, E. L., Xu, M., Wang, S., Klein, A., and Chen, Y.: Improving MODIS snow products with a HMRF-based spatio-temporal modeling technique in the Upper Rio Grande Basin, Remote Sens. Environ., 204, 568–582, https://doi.org/10.1016/j.rse.2017.10.001, 2018. 

Huang, Y., Song, Z. C., Yang, H. X., Yu, B. L., Liu, H. X., Che, T., Chen, J., Wu, J. P., Shu, S., Peng, X. B., Zheng, Z. J., and Xu, J. H.: Snow cover detection in mid-latitude mountainous and polar regions using nighttime light data, Remote Sens. Environ., 268, 112766, https://doi.org/10.1016/j.rse.2021.112766, 2022. 

Hussainzada, W., Lee, H. S., Vinayak, B., and Khpalwak, G. F.: Sensitivity of snowmelt runoff modelling to the level of cloud coverage for snow cover extent from daily MODIS product collection 6, J. Hydrol.: Reg Stud., 36, 100835, https://doi.org/10.1016/j.ejrh.2021.100835, 2021. 

Immerzeel, W. W., van Beek, L. P. H., and Bierkens, M. F. P.: Climate change will affect the asian water towers, Science, 328, 1382–1385, https://doi.org/10.1126/science.1183188, 2010. 

Jing, Y. H., Shen, H. F., Li, X. H., and Guan, X. B.: A two-stage fusion framework to generate a spatio–temporally continuous MODIS NDSI product over the Tibetan Plateau, Remote Sensing, 11, 2261, https://doi.org/10.3390/rs11192261, 2019. 

Kilpys, J., Pipiraitė-Januškienė, S., and Rimkus, E.: Snow climatology in Lithuania based on the cloud-free moderate resolution imaging spectroradiometer snow cover product, Int. J. Climatol., 40, 4690–4706, https://doi.org/10.1002/joc.6483, 2020. 

Klein, A.: Validation of daily MODIS snow cover maps of the Upper Rio Grande River Basin for the 2000–2001 snow year, Remote Sens. Environ., 86, 162–176, https://doi.org/10.1016/s0034-4257(03)00097-x, 2003. 

Kumar, L., Skidmore, A. K., and Knowles, E.: Modelling topographic variation in solar radiation in a GIS environment, Int. J. Geogr. Inf. Sci., 11, 475–497, https://doi.org/10.1080/136588197242266, 1997. 

Li, C., Su, F., Yang, D., Tong, K., Meng, F., and Kan, B.: Spatiotemporal variation of snow cover over the Tibetan Plateau based on MODIS snow product, 2001-2014, Int. J. Climatol., 38, 708–728, https://doi.org/10.1002/joc.5204, 2018. 

Li, M., Zhu, X., Li, N., and Pan, Y.: Gap-filling of a MODIS Normalized Difference Snow Index product based on the similar pixel selecting algorithm: a case study on the Qinghai–Tibetan Plateau, Remote Sensing, 12, 1077, https://doi.org/10.3390/rs12071077, 2020. 

Li, X., Jing, Y., Shen, H., and Zhang, L.: The recent developments in cloud removal approaches of MODIS snow cover product, Hydrol. Earth Syst. Sci., 23, 2401–2416, https://doi.org/10.5194/hess-23-2401-2019, 2019. 

Li, Y., Chen, Y., and Li, Z.: Developing daily cloud-free snow composite products from MODIS and IMS for the Tienshan mountains, Earth and Space Science, 6, 266–275, https://doi.org/10.1029/2018ea000460, 2019. 

Liang, H., Huang, X., Sun, Y., Wang, Y., and Liang, T.: Fractional snow-cover mapping based on MODIS and UAV data over the Tibetan Plateau, Remote Sensing, 9, 1332, https://doi.org/10.3390/rs9121332, 2017. 

Liang, T., Huang, X., Wu, C., Liu, X., Li, W., Guo, Z., and Ren, J.: An application of MODIS data to snow cover monitoring in a pastoral area: A case study in Northern Xinjiang, China, Remote Sens. Environ., 112, 1514–1526, https://doi.org/10.1016/j.rse.2007.06.001, 2008. 

Liu, C., Li, Z., Zhang, P., Zeng, J., Gao, S., and Zheng, Z.: An assessment and error analysis of MOD10A1 snow product using Landsat and ground observations over China during 2000–2016, IEEE J. Sel. Top. Appl., 13, 1467–1478, https://doi.org/10.1109/jstars.2020.2983550, 2020. 

Liu, Y., Chen, X., Hao, J.-S., and Li, L.-h.: Snow cover estimation from MODIS and Sentinel-1 SAR data using machine learning algorithms in the western part of the Tianshan Mountains, J. Mt. Sci., 17, 884–897, https://doi.org/10.1007/s11629-019-5723-1, 2020. 

Muhammad, S. and Thapa, A.: An improved Terra–Aqua MODIS snow cover and Randolph Glacier Inventory 6.0 combined product (MOYDGL06*) for high-mountain Asia between 2002 and 2018, Earth Syst. Sci. Data, 12, 345–356, https://doi.org/10.5194/essd-12-345-2020, 2020. 

Parajka, J. and Blöschl, G.: Spatio-temporal combination of MODIS images – potential for snow cover mapping, Water Resour. Res., 44, W03406, https://doi.org/10.1029/2007wr006204, 2008. 

Parajka, J., Holko, L., Kostka, Z., and Blöschl, G.: MODIS snow cover mapping accuracy in a small mountain catchment – comparison between open and forest sites, Hydrol. Earth Syst. Sci., 16, 2365–2377, https://doi.org/10.5194/hess-16-2365-2012, 2012. 

Qiu, J.: China: The third pole, Nature, 454, 393–396, https://doi.org/10.1038/454393a, 2008. 

Richiardi, C., Blonda, P., Rana, F. M., Santoro, M., Tarantino, C., Vicario, S., and Adamo, M.: A revised snow cover algorithm to improve discrimination between snow and clouds: a case study in Gran Paradiso National Park, Remote Sensing, 13, 1957, https://doi.org/10.3390/rs13101957, 2021. 

Riggs, G. A., Hall, D. K., and Román, M. O.: Overview of NASA's MODIS and Visible Infrared Imaging Radiometer Suite (VIIRS) snow-cover Earth System Data Records, Earth Syst. Sci. Data, 9, 765–777, https://doi.org/10.5194/essd-9-765-2017, 2017. 

Riggs, G. A., Hall, D. K., and Román, M. O.: MODIS Snow Products Collection 6.1 User Guide, https://modis-snow-ice.gsfc.nasa.gov/uploads/snow_user_guide_C6.1_final_revised_april.pdf (last access: 1 March 2022), 2019. 

Salomonson, V. V. and Appel, I.: Estimating fractional snow cover from MODIS using the normalized difference snow index, Remote Sens. Environ., 89, 351–360, https://doi.org/10.1016/j.rse.2003.10.016, 2004. 

Tang, Z., Wang, J., Li, H., and Yan, L.: Accuracy validation and cloud obscuration removal of MODIS fractional snow cover products over Tibetan Plateau, Remote Sensing Technology and Application, 28, 423–430, https://doi.org/10.16089/j.cnki.1008-2786.000237, 2013. 

Teilet, P. M., Guindon, B., and Goodenough, D. G.: On the slope-aspect correction of multispectral scanner data, Can. J. Remote Sens., 8, 1537–1540, https://doi.org/10.1080/07038992.1982.10855028, 1982. 

Tran, H., Nguyen, P., Ombadi, M., Hsu, K. L., Sorooshian, S., and Qing, X.: A cloud-free MODIS snow cover dataset for the contiguous United States from 2000 to 2017, Scientific Data, 6, 180300, https://doi.org/10.1038/sdata.2018.300, 2019. 

Wang, G., Jiang, L., Shi, J., Liu, X., Yang, J., and Cui, H.: Snow-covered area retrieval from Himawari–8 AHI imagery of the Tibetan Plateau, Remote Sensing, 11, 2391, https://doi.org/10.3390/rs11202391, 2019. 

Wang, X., Xie, H., Liang, T., and Huang, X.: Comparison and validation of MODIS standard and new combination of Terra and Aqua snow cover products in northern Xinjiang, China, Hydrol. Process., 23, 419–429, https://doi.org/10.1002/hyp.7151, 2009. 

Wang, X., Wu, C., Peng, D., Gonsamo, A., and Liu, Z.: Snow cover phenology affects alpine vegetation growth dynamics on the Tibetan Plateau: Satellite observed evidence, impacts of different biomes, and climate drivers, Agr. Forest Meteorol., 256–257, 61–74, https://doi.org/10.1016/j.agrformet.2018.03.004, 2018. 

Wu, Z., Jiang, Z., Li, J., Zhong, S., and Wang, L.: Possible association of the western Tibetan plateau snow cover with the decadal to interdecadal variations of northern China heatwave frequency, Clim. Dynam., 39, 2393–2402, 2012. 

Xiao, X., Liang, S., He, T., Wu, D., Pei, C., and Gong, J.: Estimating fractional snow cover from passive microwave brightness temperature data using MODIS snow cover product over North America, The Cryosphere, 15, 835–861, https://doi.org/10.5194/tc-15-835-2021, 2021. 

Xu, W. F., Ma, H. Q., Wu, D. H., and Yuan, W. P.: Assessment of the daily cloud-free MODIS snow-cover product for monitoring the snow-cover phenology over the Qinghai-Tibetan Plateau, Remote Sensing, 9, 585, https://doi.org/10.3390/rs9060585, 2017. 

Yang, J., Jiang, L., Ménard, C. B., Luojus, K., Lemmetyinen, J., and Pulliainen, J.: Evaluation of snow products over the Tibetan Plateau, Hydrol. Process., 29, 3247–3260, https://doi.org/10.1002/hyp.10427, 2015. 

Yang, K., Wu, H., Qin, J., Lin, C., Tang, W., and Chen, Y.: Recent climate changes over the Tibetan Plateau and their impacts on energy and water cycle: A review, Global Planet. Change, 112, 79–91, https://doi.org/10.1016/j.gloplacha.2013.12.001, 2014. 

Yao, T., Xue, Y., Chen, D., Chen, F., Thompson, L., Cui, P., Koike, T., Lau, W. K. M., Lettenmaier, D., Mosbrugger, V., Zhang, R., Xu, B., Dozier, J., Gillespie, T., Gu, Y., Kang, S., Piao, S., Sugimoto, S., Ueno, K., Wang, L., Wang, W., Zhang, F., Sheng, Y., Guo, W., Ailikun, Yang, X., Ma, Y., Shen, S. S. P., Su, Z., Chen, F., Liang, S., Liu, Y., Singh, V. P., Yang, K., Yang, D., Zhao, X., Qian, Y., Zhang, Y., and Li, Q.: Recent third pole's rapid warming accompanies cryospheric melt and water cycle intensification and interactions between monsoon and environment: multidisciplinary approach with observations, Modeling, and Analysis, B. Am. Meteorol. Soc., 100, 423–444, https://doi.org/10.1175/bams-d-17-0057.1, 2019.  

You, Q., Wu, T., Shen, L., Pepin, N., Zhang, L., Jiang, Z., Wu, Z., Kang, S., and AghaKouchak, A.: Review of snow cover variation over the Tibetan Plateau and its influence on the broad climate system, Earth-Sci. Rev., 201, 103043, https://doi.org/10.1016/j.earscirev.2019.103043, 2020. 

Yu, J., Zhang, G., Yao, T., Xie, H., Zhang, H., Ke, C., and Yao, R.: Developing daily cloud-free snow composite products from MODIS Terra–Aqua and IMS for the Tibetan Plateau, IEEE T. Geosci. Remote, 54, 2171–2180, https://doi.org/10.1109/tgrs.2015.2496950, 2016. 

Zhang, H., Zhang, F., Zhang, G., Yan, W., and Li, S.: Enhanced scaling effects significantly lower the ability of MODIS normalized difference snow index to estimate fractional and binary snow cover on the Tibetan Plateau, J. Hydrol., 592, 125795, https://doi.org/10.1016/j.jhydrol.2020.125795, 2021. 

Zheng, Z. and Cao, G.: Snow cover dataset based on multi-source remote sensing products blended with 1km spatial resolution on the Qinghai-Tibet Plateau (1995–2018), National Tibetan Plateau Data Center [data set], https://doi.org/10.11888/Snow.tpdc.270102, 2019. 

Download
Short summary
Reliable snow cover information is important for understating climate change and hydrological cycling. We generate long-term daily gap-free snow products over the Tibetan Plateau (TP) at 500 m resolution from 2002 to 2021 based on the hidden Markov random field model. The accuracy is 91.36 %, and is especially improved during snow transitional period and over complex terrains. This dataset has great potential to study climate change and to facilitate water resource management in the TP.