Mapping snow cover variations using a MODIS daily cloud-free snow cover product in northeast China

Abstract Cloud contamination is one of the major barriers for wider applications of MODIS snow cover products. This study presents a cloud-removal approach, through multiday backward replacements based on Terra and Aqua daily MODIS snow cover products (MOD10A1 and MYD10A1), to generate a series of daily cloud-free snow cover products for advanced applications (MODMYD_MC). The products are evaluated using in situ snow depth data measured during 2000 to 2010 at 53 weather stations in the Heilongjiang Province, northeast China. The results show that the annual mean cloud covers of MOD10A1, MYD10A1, MODMYD_DC (the daily combination of MOD10A1 and MYD10A1), and MODMYD_MC are 50%, 54%, 35%, and 0%, mean snow covers are 6%, 6%, 10%, and 19%, and their mean agreements of snow cover mapping are 42%, 40%, 51%, and 91%, respectively. The snow-covered days (SCDs) derived from MODMYD_MC are also in good agreement (91%) with those obtained from in situ observations. The MODMYD_MC snow cover images are then used to investigate the detailed variation of snow cover in the XiaoXing’AnLing watershed. The snow-covered area in the watershed has an increasing trend in the recent decade, with the minimum present in the 2002 (hydrologic year) and the maximum present in 2010. The plains with lower elevation show shorter SCD but larger interannual variations than in the mountainous areas. This study indicates that MODMYD_MC can be applied to monitor the spatiotemporal variations of snow cover in northeast China and elsewhere in the world.


Introduction
Snow has a significant impact on surface energy balance, crop frost, and water cycles due to its high albedo, thermal insulation, and melting water. 1 In situ measurements and space observations are the two most common approaches for snow monitoring. Although in situ measurements provide the most accurate snow observations, they are subject to environmental accessibility and sparse distributions. 2 Passive microwave remote sensors, for example, scanning microwave multiband radiometer, special sensor microwave/imager, advanced microwave scanning radiometer-EOS (AMSR-E), are used mostly for continental and global scale studies because of their daily repeat cycle and ability to penetrate cloud. But their coarse spatial resolution (25 km or coarser) constrains further applications for mountainous areas in local and regional scales. In contrast, the active microwave remote sensors, such as advanced synthetic aperture radar and interferometric synthetic aperture radar, have higher spatial resolution but a longer repeat cycle and are not practical for daily operation of snow cover monitoring. [3][4][5] Optical remote sensors, such as Landsat and others with higher spatial resolutions, can be used to accurately map local and regional snow cover conditions under clear sky, but are limited for daily operational use due to their longer repeat cycle and inability to penetrate cloud.
The moderate resolution imaging spectral radiometer (MODIS) onboard both Terra and Aqua satellites scans the Earth twice per day during daytime in most parts of the world. MODIS snow cover products have been widely applied for snow cover mapping and climatic and hydrologic studies, although cloud seriously impairs the advantages of MODIS images. [6][7][8][9] For instances, in Northern Xinjiang, China, and on the Colorado Plateau, USA, the annual mean cloud cover is around 50% and can be up to 95% in several straight days. 2,10 Various approaches have been developed to mitigate the influence of cloud obscuration for the MODIS standard snow cover products. The first common step of most cloud-removal approaches [10][11][12][13][14][15][16][17] is to combine the MODIS Terra and Aqua daily snow cover products to generate a daily cloud-less product (∼10% less as compared with the MODIS Terra or Aqua alone). The difference of various approaches comes from the further processing steps to remove those cloudblocked pixels. The first type is to use passive microwave snow water equivalent product (AMSR-E, 25-km spatial resolution) to replace the remaining cloud-blocked pixels. [4][5][6] This method greatly mitigates the impact of cloud but also reduces the spatial resolution of MODIS snow cover images, especially for spatially extensive cloud-blocked areas (from 500 m to 25 km). The second type is to use the snowline or spatial filter to remove the remaining cloud-blocked pixels. [11][12][13][14] The snowline method does not replace all cloud-blocked pixels with snow or land. 13 The effect of spatial filter that interpolates from nearby cloud-free pixels is limited for spatially extensive cloud-blocked areas. 18 The third type applies a temporal filter to composite multiday snow cover images (both forward and backward) to minimize the impact of cloud. [10][11][12][13][14][15][16][17] The snowline and spatial filters methods are normally applied together with or after the temporal filter to further reduce the cloud-blocked pixels.
The aforementioned methods sacrifice either spatial resolution with no remaining cloud pixels (the first type) or temporal resolution with still remaining cloud-blocked pixels (the second and third types). The forward combination is a retrospective estimate which requires data past the date and is useful for monthly or seasonal reanalysis. Dozier et al. 18 adopted a continuous backward combination approach to assimilate cloud-free pixels in their daily fractional snow cover products for runoff forecasting, but did not generate a daily cloud-free snow cover product.
In this study, we present a new multiday backward cloud-removal approach to produce daily cloud-free snow cover product. This method includes two steps. The first step is basically the same as the aforementioned methods to combine daily MODIS Terra (MOD10A1) and Aqua (MYD10A1) snow cover images to generate a daily cloudless (∼10% less) snow cover image MODMYD_DC. The second step, however, is unique to all other methods and is a multiday backward replacement, i.e., the cloud-blocked pixels on the current day combination of MODMYD_DC (n) are replaced by the cloud-free pixels on the previous day combination of MODMYD_DC (n − 1), continuing the backward replacement until all cloud-blocked pixels are replaced. The resulting image is a daily cloud-free snow cover product (i.e., MODMYD_MC). This daily cloud-free snow cover product is then validated using the in situ snow depth observations in the Heilongjiang Province and is finally used to analyze the spatiotemporal variations of snow cover in the XiaoXing'AnLing watershed in northeast China.

Data and Test Area
Terra and Aqua satellites launched in December 1999 and May 2002 pass the equator at the local time of 10:30 am and 1:30 pm, respectively. MODIS instrument is one of the key sensors mounted on both satellites. Both Terra and Aqua MODIS snow cover products apply the same algorithm and have similar accuracy on snow cover mapping. 6 MODIS standard snow cover products include the daily MOD10A1 (Terra) and MYD10A1 (Aqua) and the 8-day maximum composite of MOD10A2 and MYD10A2 with a 500-m spatial resolution and 10 deg by 10 deg tiles, and the daily, 8-day maximum and monthly mean snow cover products of MOD10C1-3 and MYD10C1-3 with a 0.05 deg global climate model pixel. 19,20 MOD10A1 and MYD10A1 are used in this study.
The MODIS tile of h26v04 covering the central part of the Heilongjiang Province and neighboring provinces is selected for the algorithm development ( Fig. 1). A total of 3962 MOD10A1 and 3429 MYD10A1 images are used during the 11 hydrologic years from 2000 to 2011 (September 1 to August 31 of the following year).
The collected in situ snow depth data are daily measurements rounded to 1 cm from 2000 to 2010. There are 53 weather stations within the umbrella of MODIS tile h26v04 in the Heilongjiang Province. The elevations of those stations range from 53 to 550 m above sea level. In comparison, MODIS pixel values collocated with one weather station are extracted to compare with the corresponding in situ snow depth data. To minimize the impacts of the spatial representative of in situ measurements, the snow depth data ≥4 cm are used to validate the accuracy of snow cover mapping, and other data of 0 cm (no snow) and 1 to 3 cm (patchy snow) are also provided for assistance. 2,6,11

New Cloud-Removal Algorithm
The new cloud-removal algorithm consists of two automatic steps (Fig. 2). The first step is to combine the snow cover images of MOD10A1 and MYD10A1 to generate a daily maximum snow cover image called MODMYD_DC. 10 The second step is a multiday backward replacement, i.e., the cloud-blocked pixels on the current day combination of MODMYD_DC (n) are replaced by the cloud-free pixels on the previous day combination of MODMYD_DC (n − 1), continuing the backward replacement on day (n − k) until all cloud-blocked pixels are replaced. There is no limit for k, the number of cloud-persistent days (CPDs), but k is normally less than 5 days.
Along with the cloud-free snow cover image, the algorithm also generates an image to record the CPD on each pixel. CPD is a relative number (days) against the current day and indicates the clear sky date of each pixel on the new cloud-free snow cover image. In other words, CPD records the days when the cloud-blocked pixel on the current day image is replaced with a cloud-free pixel on a previous day image(s) and acts like a quality flag for the new cloud-free snow cover images. For instance, if CPD ¼ 0, it indicates that this pixel is cloud free in the current day image and does not require replacement. If CPD ¼ 1, it indicates that the cloud-covered pixels on the current day are replaced as cloud free one day before the current day; if CPD ¼ 2, it indicates that the cloud-covered pixels on the current day and the day before are replaced as cloud free two days before the current day, and so on. As the CPD gets larger, the new cloud-free snow cover image has a greater uncertainty.

Efficiency of Cloud-Removal
MODIS tile h26v04 covering northeastern China is applied to test the cloud-removal algorithm and the new daily cloud-free snow cover product. As shown in Fig. 3, in the first 60 days of 2004, the daily combination algorithm removes the cloud-blocked pixels by ∼10% of the total pixels after combining MOD10A1 and MYD10A1. After 3-day replacements, the cloud cover dramatically declines to below 10% and continues to decline as continuous replacements are made. Since the MOD10A1 was the first available at the end of February 2000, a new daily cloud-free snow cover image could be generated each day at the beginning of April 2000 after the first month's initial combination.

The first step: Daily combination
The second step: Multiday replacements Fig. 2 The flowchart of the cloud-removal algorithm, consisting of two steps: the daily combination and the multiday replacement. Day (n) is the current day, day (n − 1) is one day before the current day, and k is the number of cloud-persistent days (CPDs). The annual mean cloud cover derived from the daily MODIS snow cover product is around 50% in the 11 years within the MODIS tile h26v04, leading to a large underestimate of annual mean snow cover (6%) ( Table 1). The daily combination (MODMYD_DC) of MOD10A1 and MYD10A1 reduces the annual mean cloud cover to 36% (14% to 18% cloud reduction) and increases the annual mean snow cover to 10%. After multiday replacements, the new daily cloud-free snow cover product (MODMYD_MC) retrieves an annual mean snow cover of 19%, with 68 mean snow-covered days (SCDs), and a minimum (35days) in 2002 (hydrologic year) and maximum (86 days) in 2010. Figure 4 illustrates the mean percentage of the cloud-free pixels on MODMYD_MC images against the nine CPD categories of 0, 1, 2, 3, 4, 5, 6, 7 days, and ≥8 days during 2000 to 2011. The cloud-free pixels on the current day combination images were dominant (57%), and those for CPD ¼ 1 and 2 were next at 18% and 9%, respectively. The total cloud-free pixels were 88% after the 4 days' backward replacements. This is similar to the results of Xie et al. 10 who found that a snow cover image was generated with a cloud cover of less than 10% within 2.5 days on average in Northern Xinjiang, China, and on the Colorado Plateau, USA.

Assessment of Snow Classification
The in situ snow depth measurements at 53 meteorological stations within the MODIS tile h26v04 and in the Heilongjiang Province are used to validate the agreement of snow and land classification for the four snow cover products. In comparison, the in situ snow depth data are divided into three groups: 2,6,10 (1) snow depth ≥4 cm (20% of all data), which assumes that the corresponding entire MODIS pixel is covered by snow; (2) 1 cm ≤ snow depth < 4 cm (8%), which assumes that the corresponding MODIS pixel is partly covered by snow; and (3) snow depth ¼ 0 cm (72%), no snow. Other snow depth thresholds of 1 cm, 21 2 cm, 22 and 2.54 cm 23,24 were also used in literature.
The mean agreements of snow and land classification with in situ measurements for MOD10A1 and MYD10A1 are as low as ∼40% and ∼55%, respectively, mainly due to the blockage of cloud (Fig. 5). MOD10A1 in the morning has slightly less cloud and higher agreement than MYD10A1 in the afternoon. MODMYD_DC, the daily combination of MOD10A1 and MYD10A1, reduces the cloud cover by ∼10% and correspondingly increases the agreement of snow classification to 51%. MODMYD_MC, the new cloud-free snow cover product, increases the snow and land agreements to 91% and 93%, respectively. Moreover, the number of CPD for the MODMYD_MC negatively affects the agreements of land and snow classifications. The agreement of snow classification is 95% and 93% when pixels are cloud free on the current day (CPD ¼ 0) and one day before (CPD ¼ 1), respectively. It decreases to 88% when CPD ¼ 2 and is below 80% when CPD ≥ 5 days (Fig. 4). The agreements of snow classification for the fractional snow cover conditions (snow depth ¼ 1, 2, and 3 cm) in the above four products are 22%, 13%, 24%, and 62%, respectively (not shown). The low agreements are partially due to the limitations of spatial representatives of points in situ measurements, which cannot represent the patchy snow cover situations over a MODIS pixel of 500 m.
To illustrate the impacts of the widespread forest in northeast China on the agreement of snow classification, the land cover types around each station are extracted using the MODIS land cover product of MCD12Q1 in 2009. 25 MCD12Q1 was produced using ensemble supervised approaches with different classification systems. MCD21Q1 consists of five different land Fig. 4 Mean agreements of snow and land classification against in situ snow depth measurements and the mean percentage of cloud-free pixels during 2000 to 2011. On the horizontal X -axis, CPD ¼ 0 indicates that this pixel is cloud-free on the current day and does not need to be replaced; CPD ¼ 1 indicates that the cloud-blocked pixels on the current day are replaced as cloud-free one day before; CPD ¼ 2 means that the cloud-blocked pixels are replaced as cloudfree two days before; and so on. cover classifications in each calendar year. The 17-class International Geosphere -Biosphere Programme classification is used in this study. The land cover types around those stations include deciduous coniferous forest, deciduous forest, shrub, grass, and crops and are reclassified into forest (coniferous forest, deciduous forest, and shrub) and nonforest (grass and crops). The mean agreements of snow classification over forested areas are ∼15% lower than those for nonforested areas (Fig. 6).

Assessment of Snow-Covered Days
The annual SCDs in each hydrologic year (from September 1 to August 31) are computed using the new daily cloud-free snow cover images (MODMYD_MC) from 2000 to 2011. SCDs are also validated by in situ snow depth data at the same 53 stations, but with a loose constraint to the snow depth data (≥2 cm). 22 The agreement between the two SCDs is determined by Eq. (1) 7 where A is the agreement, and SCD mod and SCD grnd represent the SCDs derived from MODIS and calculated from in situ snow depth data. The mean agreement in SCD during the 11 hydrologic years is shown in Fig. 7. SCDs varied from 50 to 150 days at those stations and most were in a range from 80 to 120 days. The MODISderived SCD is in agreement with the in situ SCD [A ¼ 91%, mean difference (MODIS − in situ) ¼ −2 days], and the former is slightly shorter than the latter partially due to the patchy snow cover. We compare MODIS-derived SCDs to in situ SCDs using different snow depth thresholds of 1, 2, 3, and 4 cm. The best results we present here are from 2 cm.

Snow Cover Variations
The in situ measurements of snow depth data are limited to represent the spatial variations of snow cover. Satellite observations expand beyond in situ locations to monitor the detailed variation of land cover including snow cover at the regional and even global scales. The new daily cloud-free snow cover images generated in this study are used to illustrate the spatiotemporal variations of snow cover in the XiaoXing'AnLing area in northeast China as a case study.
MODIS SCD describes the total SCDs for each 500-m pixel in a hydrological year from September 1st to August 31st in the following year. 7 The MODIS SCD is calculated from MODMYD_MC images. If a pixel value is 200 (snow), then this pixel is covered by snow on the current day and its SCD equals to 1; otherwise its SCD is 0. The annual SCD is obtained

Snow-Covered Days
The spatial distribution of the 11-year mean SCD is illustrated in Fig. 8 in the XiaoXing'AnLing area. SCD is classified into seven categories: shorter than 30 days, 30 to 60 days, 61 to 90 days, 91 to 120 days, 121 to 150 days, 151 to 180 days, and longer than 180 days. Most pixels have an SCD in a range from 60 to 150 days. The mountainous areas have a much longer SCD than the plains.
The spatial distribution of SCD anomalies (deviation from the 11-year mean) of each hydrologic year are demonstrated in Fig. 9. The SCD anomalies are also divided into seven groups with different colors: shorter than −60 days, −60 to −31 days, −30 to −11 days, −10 to 10 days (black, no change), 11 to 30 days, 31 to 60 days, and longer than 60 days. Red colors are below the mean and the green colors are above the mean. SCDs in hydrologic years of 2000 to 2001 (not shown in Fig. 9   when SCDs were 30 days longer than the mean in almost the entire area. Overall, SCDs demonstrated an increasing pattern during the 11 hydrologic years in this area. This increase trend is also observed by Chen et al. 22 and seems negatively correlated with the air temperature decrease in the recent decade. Figures 10 and 11 display the histogram of SCD and SCD anomalies in mountainous areas (elevation < 500 m) and in plains (elevation ≤ 500 m). In plains, the dominant SCD (total 37%) ranges from 100 to 130 days (median ¼ 103), whereas in mountainous areas, the primary SCD (56%) is between 140 and 170 days (median ¼ 146). The SCD in mountainous areas is less variable than those in plains.
The mode values for SCD anomalies (Fig. 11) are 0 to 20 days (48%) in mountainous areas, while they are 5 to 25 days (34%) in plains. The frequencies of less than −20 days (23%) and more than 30 days (11%) in the plains are higher than those in mountainous areas with frequencies of 16% and 4%, respectively. This indicates that the SCD both in plains and mountainous areas shows an increasing trend, being consistent with the spatial distribution of SCD anomalies in Fig. 9, and the interannual variations of snow cover in plains are larger than those in mountainous areas.

Variations of Snow Cover Index
This study applies a snow cover index (SCI) in Eq. (2) to integrate the spatiotemporal snow cover variations in a certain region. 7 where A represents the area of each pixel (0.25 km 2 ), SCD i refers to the SCD in a hydrological year for pixel i, and N is the total pixels in the study area, i.e., the XiaoXing'AnLing watershed labeled by the white polygon in Fig. 8 in this study. SCI in Fig. 12 shows a similar pattern and an increasing trend with SCD anomalies in Fig. 9.

Discussion
Our improved cloud-removal approach generates a daily cloud-free MODIS snow cover product through a daily combination and multiday backward replacements. The daily combination of Terra and Aqua MODIS snow cover products is based on the snow priority algorithm developed by Xie Fig. 11 The histogram of SCD anomalies in mountainous areas (elevation > 500 m) and in plains (elevation ≤ 500 m) from 2000 to 2011. et al. 10 Other studies that applied this algorithm did not generate a completely cloud-free snow cover product although the cloud-blocked pixels were greatly reduced. [15][16][17] The cloud threshold set by Xie et al. 10 in their final cloud-less snow cover product was 10%. Thus, they generated a cloud-less snow cover image with cloud <10% per 2.5 days on average, and with mean cloud of ∼6%, which is similar to that of MODIS 8-day snow cover products. Hall et al. 15 made the combination by replacing the cloud pixels with the cloud-free pixels on the previous 1 to 3 days' images, but only the 0.05 deg (5 km) snow cover products (MOD10C1) were tested. Gao et al. 17 further improved Hall et al.'s 15 algorithm by extending the snow combination window to both the previous and succeeding 3 to 5 days to produce a cloud-less snow cover product, yet still not generating a daily completely cloud-free snow cover product. The total snow cover retrieved from Gao et al.'s approach 17 and from ours is similar by a difference of 0.5% to the 2003 to 2004 hydrologic year. Our new cloud-removal approach, through multiday backward replacements, not only produces a daily cloud-free snow cover product, but is also a predictive estimate that can be made in near real time, while the forward combination is a retrospective estimate that requires data past the date and is only useful for monthly or seasonal reanalysis. 18 The CPDs images are generated at the same time when the cloud-free snow cover images are produced. CPD acts as a quality flag to indicate the specific date when a cloud-free pixel is used to replace the cloud-blocked pixel on the current day. The longer the CPD, the larger uncertainty the replaced pixels have (Fig. 4). In the new cloud-free snow cover product, the cloud-free pixels are dominant (57%) on the current day image and on the previous 3-day images (31%). Therefore, compared with the former methods, the improved cloud-removal algorithm makes full use of those cloud-free pixels on both Terra and Aqua MODIS snow cover images and generates a daily completely cloud-free snow cover product. Although northeast China has the thickest snow depth accumulation among the three snow distribution regions in China, few studies investigate the detailed snow variation using the updated remote sensing techniques. Lei et al. 26 evaluated the snow classification accuracy of the MODIS 8-day snow cover products and AMSR-E SWE in the Heilongjiang Province during 2002 to 2007. Chen et al. 22 examined the spatiotemporal variations of snow cover in northeast China by using the multiday combinations of MODIS snow cover products, which have a mean cloud cover of 4% and mean combination period of 7 days. Lei et al. 26 directly applied the daily MODIS Terra snow cover product to compute the SCDs in northeast China by only counting the snow-covered pixels under a clear sky and omitted those blocked by cloud but which were covered with snow. Cloud is the major factor influencing the MODIS snow cover classification, leading to an underestimation of the snow-covered areas. In contrast, SCD computed from the daily cloud-free snow cover product in this study can better represent the real snow-covered duration (Table 1). This study illustrates the XiaoXing'AnLing watershed as an example to analyze the detailed snow cover variations using our daily cloud-free snow cover product. The algorithm developed in this study can also be applied to other areas in China and elsewhere in the world. The daily cloud-free snow cover product and other derived maps provide basic information for snow disaster prevention and mitigation, agriculture, and water resources management and are accurate inputs for hydrologic and climatic models as well. 22 The snow cover shows an overall upward trend in the XiaoXing'AnLing watershed. The increase trend of snow cover in northeast China negatively correlated with the air temperature decrease in the recent decade. 22 Mountainous areas with higher elevation had longer SCD, but shorter SCD anomalies than in plains. These demonstrate the possible applications of the daily cloud-free snow cover product. In addition, it can also be combined with other microwave snow products to downscale their spatial resolution. [3][4][5] Logically, a next step is to apply this cloud-free snow cover product to investigate the detailed variations of snow cover in the entire northeast China and in the Amur River Basin, which is beneficial to the local hydrology, water resources, climate, agriculture, and other fields.

Summary
This study develops a cloud-removal approach via multiday backward replacements based on Terra and Aqua MODIS daily snow cover products (MOD10A1 and MYD10A1), to generate a daily cloud-free snow cover product (MODMYD_MC). MODIS tile h26v04 was used to test the proposed algorithm, and the four snow cover products were evaluated with in situ snow depth measured from 2000 to 2010 at 53 meteorological stations in the Heilongjiang Province in northeast China. The annual mean cloud covers of MOD10A1, MYD10A1, MODMYD_DC, and MODMYD_MC are 50%, 54%, 35%, and 0%, the mean snow covers are 6%, 6%, 10%, and 19%, and their mean agreements of snow cover mapping are 42%, 40%, 51%, and 91%, respectively. The SCDs derived from MODMYD_MC are in agreement (91%) with and slightly shorter than those obtained from in situ observations. This study also illustrates the detailed variations of snow cover in the XiaoXing'AnLing watershed using SCD and SCI derived from MODMYD_MC and finds that the snow cover in this area had an increasing trend since 2000, with a minimum in the 2002 hydrologic year and a maximum in 2010. The results from this study demonstrate that the MODMYD_MC can be applied to monitor the spatiotemporal variations of snow cover in northeast China and elsewhere in the world. Moreover, the improved cloud-removal algorithm makes full use of cloudfree pixels on both Terra and Aqua MODIS snow cover images and can be applied to assimilate the latest MODIS daily snow cover images to realize the near real-time snow cover monitoring operationally.