Estimating fractional snow cover from passive microwave brightness temperature data using MODIS snow cover product over North America

The dynamic characteristics of seasonal snow cover are critical for hydrology management, the climate system, and the ecosystem functions. Optical satellite remote sensing has proven to be an effective tool for monitoring global and regional variations in snow cover. However, accurately capturing the characteristics of snow dynamics at a finer spatiotemporal resolution continues to be problematic as observations from optical satellite sensors are greatly impacted by clouds and solar illumination. Traditional methods of mapping snow cover from passive microwave data only provide binary information at a spatial resolution of 25 km. This innovative study applies the random forest regression technique to enhanced-resolution passive microwave brightness temperature data (6.25 km) to estimate fractional snow cover over North America in winter months (January and February). Many influential factors, including land cover, topography, and location information, were incorporated into the retrieval models. Moderate Resolution Imaging Spectroradiometer (MODIS) snow cover products between 2008 and 2017 were used to create the reference fractional snow cover data as the “true” observations in this study. Although overestimating and underestimating around two extreme values of fractional snow cover, the proposed retrieval algorithm outperformed the other three approaches (linear regression, artificial neural networks, and multivariate adaptive regression splines) using independent test data for all land cover classes with higher accuracy and no out-of-range estimated values. The method enabled the evaluation of the estimated fractional snow cover using independent datasets, in which the root mean square error of evaluation results ranged from 0.189 to 0.221. The snow cover detection capability of the proposed algorithm was validated using meteorological station observations with more than 310 000 records. We found that binary snow cover obtained from the estimated fractional snow cover was in good agreement with ground measurements (kappa: 0.67). There was significant improvement in the accuracy of snow cover identification using our algorithm; the overall accuracy increased by 18 % (from 0.71 to 0.84), and the omission error was reduced by 71 % (from 0.48 to 0.14) when the threshold of fractional snow cover was 0.3. The experimental results show that passive microwave brightness temperature data may potentially be used to estimate fractional snow cover directly in that this retrieval strategy offers a competitive advantage in snow cover detection.

1 Introduction

Snow cover is a critical indicator of climate change, playing a vital role in the global energy budget (Flanner et al., 2011), water cycle (Gao et al., 2019), and atmospheric circulation (Henderson et al., 2018). Snow cover directly modulates the release of carbon and methane from underlying soil (Zhang, 2005; Zona et al., 2016) and influences permafrost conditions and active layer dynamics (Zona et al., 2016). Snowpack also stores a huge amount of water, providing for both domestic and industrial water needs (Sturm, 2015; Cheng et al., 2019). Accurate and timely monitoring of the spatiotemporal variation in snow cover is beneficial for hydrologic forecasting, climate predictions, and water resources management (Barnett et al., 2005; Bormann et al., 2018).

Snow cover data are typically obtained from meteorological stations or in situ manual measurements, which is spatially discontinuous and labor intensive. Remote sensing has become an attractive alternative tool to ground-based measurements as it is able to cover a wide area and is capable of high-frequency observations. Numerous studies have focused on snow cover detection, and snow cover products used optical and microwave satellite data (Tsai et al., 2019; Liu et al., 2018; Hori et al., 2017). Most of these snow cover products provide binary information at the pixel level: snow-covered or snow-free. However, snow cover often varies within a limited scale area and is characterized by high spatial heterogeneity, especially in alpine terrain areas. Dobreva and Klein (2011) demonstrated that the use of binary snow cover classification in snow cover area estimation may produce considerable uncertainties. Binary snow cover lacking fractional features hinders an accurate characterization of the spatial distribution of snow cover and cannot accurately capture variations in seasonal snow cover dynamics. In terms of the energy budget perspective, binary snow cover introduces significant uncertainties into the global energy budget estimation because of the large surface albedo differences between snow-covered and snow-free surfaces (He et al., 2014). Thus, there is an urgent need to acquire the snow cover area at the sub-pixel level to provide accurate snow cover information. Fractional snow cover allows for the derivation of snow cover area at the sub-pixel level; this is a better option compared to binary snow cover (Salomonson and Appel, 2004).

Fractional snow cover maps derived from optical imagery have been produced for over 40 years. Optical satellite observations have been recognized for their suitability in estimating fractional snow cover because of their high spatial resolution. Moderate- to high-resolution optical observations have been popular in previous snow cover studies, including, for example, Fengyun (FY) series sensors (0.5–4 km) (Wang et al., 2017), Moderate Resolution Imaging Spectroradiometer (MODIS) (500 m) (Kuter et al., 2018), and Landsat (30 m) (Berman et al., 2018). There are also many predictive methods for fractional snow cover, such as linear regression (Salomonson and Appel, 2004, 2006), spectral mixture analysis (Wang et al., 2017; Rosenthal and Dozier, 1996), machine learning (e.g., artificial neural network, ANN) (Liang et al., 2017; Moosavi et al., 2014), and multivariate adaptive regression splines (MARS) (Kuter et al., 2018). A simple linear regression cannot fully describe the complex relationship between satellite observations and fractional snow cover. As such, nonlinear approaches have recently been developed to replace this traditional method (Berman et al., 2018). Kuter et al. (2018) estimated fractional snow cover from MODIS data using the MARS technique, for which the Landsat 8 binary snow cover data served as the reference fractional snow cover data. They found that the estimated fractional snow cover using the MARS method was in good agreement with the reference fractional snow cover, with the average correlation coefficient being R=0.93 (Kuter et al., 2018). However, polar regions contend with clouds and limited solar illumination which are the greatest challenges for snow cover detection using optical satellite data. This has resulted in snow cover maps with incomplete spatial coverage, at times with gaps of up to 70 % (Parajka and Blöschl, 2008). Although there have been constant efforts to fill the gaps mainly caused by cloud contamination by fusing multi-source data (Chen et al., 2018) such as passive microwave snow cover products (Hao et al., 2018; Huang et al., 2016) and different spatiotemporal information on snow cover (Dong and Menzel, 2016; Gafurov and Bárdossy, 2009), most studies have focused on binary snow cover.

When there are consecutive cloudy days, the use of data fusion technology introduces significant uncertainties in detecting snow cover from optical imagery. Passive microwave sensors are largely advantageous because they have the capacity to measure microwave radiation emitted from the ground under the clouds and in darkness. Compared to active microwave sensors, passive sensors have a large swath width and generate a large amount of daily observations that extend for several decades (Cohen et al., 2015). To date, passive microwave brightness temperature data have been widely applied in monitoring soil moisture (Qu et al., 2019), sea/lake ice (Peng et al., 2013), frozen soil (Han et al., 2015), and snow cover. Previous studies on snow cover have typically focused on snow depth (Xiao et al., 2018; Che et al., 2008), snow water equivalent (SWE) (Takala et al., 2011; Lemmetyinen et al., 2018), and snow cover area (Liu et al., 2018; Xu et al., 2016). All studies on snow cover area were limited to binary information. Specifically, they involved the application of common passive microwave snow cover mapping algorithms, such as Grody's algorithm (Grody and Basist, 1996), National Aeronautics and Space Administration (NASA) Advanced Microwave Scanning Radiometer–Earth Observing System (AMSR-E) SWE algorithm (Kelly, 2009), Singh's algorithm (Singh and Gan, 2000), Neale's algorithm (Neale et al., 1990), the FY3 algorithm (Li et al., 2007), and the South China algorithm (Pan et al., 2012). All these algorithms utilize different thresholds for brightness temperature to identify binary snow cover. Recently, Xu et al. (2016) applied the brightness temperatures of different channels and their linear combinations into the presence and background learning (PBL) algorithm to identify global binary snow cover.

As the effect of environmental factors (e.g., vegetation, topography, and wind) on snow cover distribution produces great heterogeneity, snow cover monitoring still bears larger uncertainties when only using passive microwave data. These large uncertainties may result from “patchy” (shallow/discontinuous) snow cover and the use of coarse resolution (25 km) (Xiao et al., 2018). Despite the coarse resolution of passive microwave sensors, their ability to detect snow cover in the presence of clouds demonstrates their effectiveness as a snow cover monitoring tool. There is an urgent need for daily time series and full space-covered sub-pixel snow cover area data for climate and reanalysis studies. Thus, it is necessary to derive high resolution fractional snow cover that can describe snow cover distribution patterns and capture its rapid evolution processes. Brodzik et al. (2018) recently published the Calibrated Enhanced-Resolution Passive Microwave Daily Equal-Area Scalable Earth Grid (EASE-Grid) 2.0 Brightness Temperature data (see Sect. 2.1 below), which have high spatial resolution (3.125 and 6.25 km) depending on frequency (Brodzik et al., 2018; Long and Brodzik, 2016). This passive microwave data with enhanced resolution provide an opportunity for fractional snow cover estimation.

The main objective of this study is to develop a feasible method utilizing the enhanced-resolution passive microwave brightness temperature data to retrieve daily fractional snow cover at a 6.25 km resolution. The datasets used in this study are described in Sect. 2, including the enhanced-resolution passive microwave data, ground-based measurements, MODIS snow cover and land cover products, and topographic data. Section 3 details the proposed retrieval algorithm using the random forest method as a retrieval function. Section 4 presents the results from the comparison of methods, evaluation, and validation experiments. Finally, Sect. 5 discusses the possible factors that impact on the accuracy of the fractional snow cover estimates derived from passive microwave data.

2 Datasets

2.1 The enhanced-resolution passive microwave data

The NASA Making Earth System Data Records for Use in Research Environments (MEaSUREs) program provides a new version of passive microwave brightness temperature data known as the Calibrated Enhanced-Resolution Passive Microwave Daily (EASE-Grid) 2.0 Brightness Temperature. This passive microwave gridded data span from 1978 to mid-2017 using the level 2 satellite records from multiple passive microwave sensors (Brodzik et al., 2018). This enhanced-resolution data may be downloaded from the National Snow and Ice Data Center (NSIDC; https://nsidc.org/data/NSIDC-0630/versions/1, last access: 10 February 2021). To explore the feasibility of estimating fractional snow cover using passive microwave data, this study mainly selected January and February of 2008–2017 as the study period. The Special Sensor Microwave Imager/Sounder (SSMIS) sensor (F-16) used in this present study offers three channels (19, 37, and 91 GHz) in both horizontal (H) and vertical (V) polarization and 22 GHz with vertical polarization. These datasets were gridded onto EASE-Grid 2.0 projections at two spatial resolutions (19 and 22 GHz with 6.25 km and 37 and 91 GHz with 3.125 km). Only observations from descending orbit (morning) were used to avoid the effects of wet snow as much as possible (Derksen et al., 2000). To achieve a common resolution, we aggregated the 3.125 km spatial resolution data to 6.25 km by averaging the surrounding four pixels.

2.2 Ground measurements

Although ground measurements of snow cover have limited spatial representation in passive microwave coarse spatial resolution, in situ measurements continue to be the most authentic and reliable data source for snow depth estimation or snow cover detection (Chen et al., 2018; Sturm et al., 2010). Ground measurements from the Global Historical Climatology Network-Daily (GHCN-Daily) data were used to assess the snow cover detection capability (Menne et al., 2012a). The GHCN-Daily dataset was provided by the National Climatic Data Center (available in https://doi.org/10.7289/V5D21VHZ, Menne et al., 2012a), and it integrates daily observations from approximately 30 different data sources. The new version of the dataset was updated on 13 June 2018 and contained measurements from over 100 000 stations worldwide. These stations record various aspects of meteorological observations, including snow depth and snowfall (Menne et al., 2012b). More than 50 000 measurement sites across Canada and the United States were collected, and all available records applied in the validation stage are from approximately 18 000 sites.

2.3 MODIS land surface products

2.3.1 Snow cover product

MODIS snow cover products were considered the most suitable reference data because of their wide application, high accuracy (Hall and Riggs, 2007; Zhang et al., 2019; Coll and Li, 2018), and high spatiotemporal resolution (1 d; 500 m). The accuracy of the version 6 MODIS snow cover products has improved compared to that of version 5 (Dong et al., 2014; Huang et al., 2018). The most noticeable change for version 6 is that the Normalized Difference Snow Index (NDSI) snow cover has replaced fractional snow cover, while binary snow-covered area (SCA) datasets are no longer available (Riggs and Hall, 2016). A snow cover detection method using NDSI was applied in version 6 to alleviate commission errors (Riggs et al., 2017). The NDSI index helps to distinguish snow from other surface features and to describe the presence of snow (Hall et al., 1998, 2001). These products are available from NSIDC website (MOD10A1: https://nsidc.org/data/MOD10A1, last access: 10 February 2021; MYD10A1: https://nsidc.org/data/MYD10A1, last access: 10 February 2021) (Hall and Riggs, 2016a, b). The local equatorial crossing times of MODIS on board the Terra and Aqua satellites are approximately 10:30 and 13:30 local time, respectively. This study used both MOD10A1 and MYD10A1 NDSI snow cover products to generate reference fractional snow cover for North America. The NDSI snow cover data were initially converted to binary snow cover for aggregation into fractional snow cover data at 6.25 km spatial resolution (see Sect. 3.2).

2.3.2 Land cover product

Generally, the retrieval accuracy of snow cover parameters is strongly dependent on the land cover type (Xiao et al., 2018; Kuter et al., 2018; Dobreva and Klein, 2011; Huang et al., 2018). We indirectly considered the land cover effect when estimating fractional snow cover by establishing retrieval models on different land cover classes derived from MODIS land cover data (2008–2017). MODIS Land Cover Type Yearly Product (MCD12Q1, version 6) incorporates five different classification schemes and is globally available at a 500 m spatial resolution spanning 2001 to the present (https://search.earthdata.nasa.gov/, last access: 10 February 2021). The International Geosphere–Biosphere Program (IGBP) classification scheme categorizes land cover in 17 classes (Sulla-Menashe and Friedl, 2018). In this study, MCD12Q1 data were resampled onto the 6.25 km grid using a simple majority method, then integrated into five classes: forest, shrub, prairie, bare land, and water (see Xiao et al., 2018). Fractional snow cover retrieval models were established for four of these land cover types, excluding water.

2.4 Topographic data

Previous studies have demonstrated that topography plays an important role in snowpack distribution (Dai et al., 2017) and snow evolution (Savoie et al., 2009). The ETOPO1 data were used as the topographic auxiliary data; these data have a 1 arcmin spatial resolution and were developed by the National Geophysical Data Center of the National Oceanic and Atmospheric Administration (NOAA) (Amante and Eakins, 2009). This study also considered elevation, slope, and aspect factors. Elevation was directly acquired from ETOPO1, which was re-projected and resampled onto the grid at 6.25 km spatial resolution. The slope and aspect data were obtained from ETOPO1 data by ArcGIS 10. Figure 1 shows the elevation pattern for North America, limited to Canada and the United States in this study.

Figure 1 Topographic map of North America.

3 Methodology

Microwave radiation constantly emitted from the substratum can be measured by passive microwave sensors. However, the overlying snow pack attenuates the upward microwave radiation (Chang et al., 1987). This microwave radiation attenuation was mainly dominated by volume scatter relying on properties of the snow cover. Previous studies have demonstrated that there is great heterogeneity in snow properties and the distribution of snow cover, both of which may be influenced by many factors (Xiao et al., 2020), including the most prevalent land-cover (Che et al., 2016; Kim et al., 2019), topography (e.g., elevation, topographic relief) (Smith and Bookhagen, 2016; Revuelto et al., 2014), time (Sturm et al., 2010; Dai et al., 2012), and climatic conditions (e.g., wind speed, near-surface soil temperature, and air temperature) (Dong et al., 2014; Grippa et al., 2004; Josberger and Mognard, 2002). Satellite sensors receive reduced upwelling microwave radiation in proportion to a greater snow cover area or a larger mass of snowpack (Chang et al., 1987; Dietz et al., 2011; Saberi, 2019). A number of published works have demonstrated the potential to derive snow depth and SWE using passive microwave radiation data (Kim et al., 2019; Wang et al., 2019). Despite the high uncertainties associated with snow depth and SWE estimations, using passive microwave data can provide useful snow cover extent information (Brown et al., 2010; Foster et al., 2011).

3.1 Overview

To develop a fractional snow cover prototype retrieval method combined with optical and passive microwave data, we only used the January and February datasets as snow cover areas are at a maximum and snowpack properties are relatively stable during this period (Xiao et al., 2018). The influential factors on snow cover, including topography factor, land cover, location, and time, were indirectly or directly considered during retrieval of the fractional snow cover. To date, many researchers have applied machine learning techniques for the retrieval of snow cover parameters to explore the relationship between passive microwave signals and snow properties (Xiao et al., 2018; Tedesco et al., 2004). Random forest is an ensemble learning method gaining the attention of many researchers because it is more efficient and robust than the single method (Breiman, 2001). As a classifier, random forest has been successfully employed to detect snow cover (Tsai et al., 2019), land cover (Rodriguez-Galiano et al., 2012), and woody invasive species (Kattenborn et al., 2019). The random forest regression method can also successfully estimate land surface temperature (Zhao et al., 2019), biomass (Mutanga et al., 2012), and soil moisture (Qu et al., 2019). In this study, random forest regression (described in Sect. 3.4.4) was selected as the retrieval method to mine the relationship between passive microwave brightness temperature and fractional snow cover. We also compared random forest with three other methods (linear regression, MARS, and ANN) widely used in fractional snow cover retrieval from optical remote sensing data in model performance. Figure 2 provides an overview of the workflow that consists of four parts.

Figure 2 Workflow diagram illustrating the processing of fractional snow cover retrieval.

First, a ground “truth” observation was necessary to produce snow cover areas at sub-pixel levels. Under clear-sky conditions, the reference fraction of snow cover was generated within a 6.25 km pixel cell by applying the aggregation method to the MODIS binary snow map (see Sect. 3.2). To make the experiment fully independent, the reference fractional snow cover data were divided into three parts: the data from 2011 to 2016 used in the training stage; the data from 2010 used in the testing stage; and the independent datasets from January and February (2008, 2009, 2017) and from December (2007, 2008, 2016) used in the evaluation stage.

Second, to the best of our knowledge, there are few attempts to directly develop fractional snow cover from passive microwave brightness temperature data. This meant a series of sensitivity experiments of input variable selection were required. Input parameters were selected based on a series of tests described in Sect. 3.3.1. Moreover, we conducted several sensitivity experiments to determine the optimal training sample size for the retrieval method used in this study (Sect. 3.3.2).

Third, many studies found that the separate estimation of fractional snow cover (Dobreva and Klein, 2011) and snow depth (Xiao et al., 2018) on different land cover types produced better results than those obtained from the combined retrieval model. As such, the random forest models were developed separately for the four land cover types.

Fourth, the last stage consisted of the evaluation and validation of the established model. Data from 2010 were used to assess the performance of four different approaches (linear regression, MARS, ANN, and random forest) for estimating fractional snow cover. Additionally, the independent datasets were used to evaluate the performance of the random-forest-based retrieval algorithm for the four land cover types. Independent validations of snow cover detection capability were conducted using the fractional snow cover retrieval results and station snow depth measurements across North America. They were compared with the results of Grody's snow cover mapping algorithm.

3.2 Preprocessing of MODIS snow cover products

The base data for this study were the reference fractional snow cover data obtained from the interpretation of MODIS snow cover products. The highest priority was to produce daily binary snow cover area from NDSI snow cover. Previous snow cover detection studies recommend a 0.4 NDSI threshold for global- and regional-scale snow cover investigations (Parajka et al., 2012; Hall et al., 1995); however, for the new version of MODIS snow cover products, several studies employed a threshold of NDSI > 0 to identify snow cover (Dong et al., 2014; Riggs et al., 2017; Huang et al., 2018). The NDSI of other features (e.g., cloud-contaminated pixels at the edges of clouds, salt pans, and pixels with very low visible reflectance) may also be greater than 0 (Riggs et al., 2017). For this reason, Zhang et al. (2019) demonstrated that a 0.1 NDSI threshold was more reasonable than 0.4 for snow cover identification in non-forest regions, whereas forest-covered regions had insufficient station measurements for a reliable and complete evaluation. MODIS snow cover performance is better for non-forest landscapes than their forest-covered counterparts, where it is less accurate for snow cover identification (Hall and Riggs, 2007; Parajka et al., 2012).

This study selected conservative NDSI thresholds of 0.1 and 0.4 for non-forest-covered areas and for forest-covered areas, respectively (Riggs and Hall, 2016) to determine snow-covered and snow-free areas. The original NDSI snow cover layer classes were reclassified into five types: snow-covered, snow-free, water, cloud, and fill value (refer to Table S1 in the Supplement). In addition, MCD12Q1 datasets (500 m) were used as auxiliary data to mask water bodies (Fig. 3) in order to alleviate the uncertainty caused by frozen water bodies when using passive microwave data to detect snow cover (Tedesco and Jeyaratnam, 2016). The MODIS binary snow cover data were generated based on the NDSI snow cover basic quality assessment (QA) with values of 0 (best), 1 (good), and 2 (OK) (Liang et al., 2017).

Figure 3 The generation of MODIS fractional snow cover.

Despite the high spatiotemporal resolution and overall accuracy of snow cover detection (85 %–99 %) using MODIS snow cover products (Parajka et al., 2012; Tran et al., 2019; Zhang et al., 2019), the cloud effect hinders its widespread applicability. Previous studies have reported that clouds may cover more than 40 % of MODIS snow cover data, in some cases exceeding 60 % (Dong and Menzel, 2016; Yu et al., 2016; Parajka and Blöschl, 2006). As such, cloud removal processing is essential to mitigate the cloud obstruction of MODIS products. This study adopted the cloud removal method combining the MOD10A1 and MYD10A1 snow cover products, as proposed by Gafurov and Bárdossy (2009). This method consists of two main filters, as shown in Fig. 3.

1. Combining snow cover images from two sensors on a given day. The first simple filter was applied under the assumption that snowmelt and snowfall did not occur within the two sensor observations. Whether a pixel in the Terra ($S_t^{\text{Terra}}$) or Aqua ($S_t^{\text{Aqua}}$) snow cover image on a given day (t) was observed as snow-covered or snow-free, the pixel in the output image (MCD10A1) was assigned the same ground status (shown in Eq. 1). The results showed that about 3 % of cloud cover was removed compared to MOD10A1 (Gafurov and Bárdossy, 2009).

2. Short-term temporal filter. If the status of a pixel in the input image (MCD10A1) on a given day (t) was cloud-covered and both the preceding (t−1) and succeeding (t+1) days were snow-covered (or snow-free), the pixel in the output image (MCTD10A1) on the given day (t) was assigned as snow-covered (or snow-free) (summarized by Eq. 2). Compared to the first filter, this short-term temporal filter may markedly reduce the number of days (10 %∼40 %) for cloud coverage and increase the overall accuracy of snow cover detection (Gafurov and Bárdossy, 2009; Tran et al., 2019).

   $$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{S}_{(\text{output},t)}=max\left(S_t^{\text{Aqua}},S_t^{\text{Terra}}\right),\text{(2)}& {\displaystyle }& {\displaystyle }{S}_{(\text{output},t)}=\mathrm{1}\text{ if }\left(S_{(t-\mathrm{1})}=\mathrm{1}\text{ and }{S}_{(t+\mathrm{1})}=\mathrm{1}\right),\end{array}$$

   where t is the time, S represents the ground status observed in the image (0 or 1), 0 denotes cloud presence, and 1 indicates snow-covered or snow-free.

Theoretically, the MODIS fractional snow cover map should calculate the percentage of snow cover in a strictly delimited area of the passive microwave pixel. Calculated areas should have a larger footprint area than the pixel resolution to avoid MODIS geolocation uncertainties (Wolfe et al., 2002; Dobreva and Klein, 2011). In this study, a window of 15×15 pixel of MODIS binary snow cover data (MCTD10A1; 500 m) was used to calculate the fraction of snow cover in a 6.25 km pixel. We adopted the most rigorous pixel filtering rule, by which one clouded pixel cannot be allowed within a 15×15 pixel window. This is slightly different from a previous study that allowed 10 % of clouds (Dai et al., 2017).

3.3 Sensitivity study

3.3.1 Selecting input variables

After determining the retrieval function, selecting the fewest number of variables to establish an efficient estimation model is a major challenge (Mutanga et al., 2012). Many factors influence snowpack distribution, and the consideration of all factors in the properties of snow cover estimation is unrealistic. Therefore, we conducted six scenarios to evaluate and finally screen input variables. The topographic factors (digital elevation model, DEM, slope, aspect) (Revuelto et al., 2014) and location information (longitude and latitude) (Xiao et al., 2018; Sturm et al., 2010) were directly taken as the basic input variables. Additionally, consideration was also given to the passive microwave brightness temperature (19, 37, and 91 GHz; both H and V polarization) (Xiao et al., 2018; Xu et al., 2016) and the difference in brightness temperature between different channels (Xu et al., 2016; Liu et al., 2018) (listed in Table 1). The 22 GHz channel was excluded because it is sensitive to water vapor.

Table 1 The input variables list. The line means this variable is not selected; the asterisk indicates the variable is selected. The numbers in brackets denote the number of variables. T19H is the brightness temperature (T) of the 19 GHz channel with H polarization. For example, T_19V_19H denotes the difference in brightness temperature between 19V and 19H channel; the others are similarly defined.

A decision tree was established using all variables shown in Scenario 1 (Table 1) and was utilized to be compared with five scenarios in terms of prediction performance and efficiency. Note that these 19 input variables were determined by using the correlation attribute evaluation method in the Waikato Environment for Knowledge Analysis 3.8.3 (WEKA) data mining software. This method evaluates the importance of the attribute by measuring the correlation between the attribute and the target (Frank et al., 2004; Witten et al., 2016). The brightness temperature and its linear combination can also directly be used to detect snow cover based on the study of Xu et al. (2016); thereby, Scenario 2 only contained brightness temperature and its linear combination without consideration of the effects of location and topographic factors. Wiesmann and Mätzler (1999) reported that V and H polarizations were dominated by scattering and snow stratigraphy, respectively. Thus, Kim et al. (2019) only assimilated V polarization with an ensemble snowpack model to estimate snow depth. Therefore, in Scenario 3, we attempted to evaluate the performance of the established retrieval model by only using the brightness temperatures for 19, 37, and 91 GHz (V polarization) based on Wiesmann and Mätzler (1999) and Kim et al. (2019). In Scenario 4, we used similar input variables to those used for snow depth estimation in Xiao et al. (2018) and examined whether these same parameters can or cannot estimate the fractional snow cover. In Scenario 5, unlike the variables used in Scenario 4, we attempted to use the basic input variables coupled with the brightness temperature linear combination for fractional snow cover retrieval.

There are other variable selection strategies based on the importance rank when using the random forest method. For example, Mutanga et al. (2012) implemented a backward feature elimination method to progressively eliminate less important variables, whilst Nguyen et al. (2018) summarized the grade of the variable and selected the top eight important variables as the input variables in the training model. Similarly, this study assessed the importance of input variables on four land cover types using the same size as the training sample (15 000) (Xiao et al., 2018). We then counted the number of times that each variable was ranked in the top nine important variables (summarized in Table S2 in the Supplement), which were then used as the input variables for Scenario 6 (listed in Table 1). By assessing the performance of models established by these six scenarios, an optimal combination of input variables for the fractional snow cover retrieval model may be selected (see Sect. 4.1.1). All input variables were normalized to [0, 1].

3.3.2 Determining sample size

Although the random forest method can avoid overfitting (Breiman, 2001), it is important to evaluate its sensitivity to sample selection types and the size of the training sample (Belgiu and Drăgu, 2016; Millard and Richardson, 2015; Nguyen et al., 2018; Colditz, 2015). The performance of predicted models trained by machine learning methods is strongly dependent on the quality of the training sample (Dobreva and Klein, 2011). Good quality training samples indicate that the sample data are not biased towards a certain value. The distribution of the fractional snow cover value from our dataset shows that more than 70 % of values were near 0 and 1. As such, the use of the random selection or equal proportional selection method (Millard and Richardson, 2015; Lyons et al., 2018; Nguyen et al., 2018) would hinder the interpretation of the final fractional snow cover estimation model by reducing the accuracy of the estimation. To address this, we adopted stratified random sampling as a sample selection strategy (Xiao et al., 2018; Dobreva and Klein, 2011), in which stratification was performed on the fractional snow cover at 0.01 increments.

From previous studies, the sample size, approximately 0.25 % of the total study area, was adopted by Colditz (2015) when using the random forest method. This value has also been evaluated in optical and active remote sensing studies (Nguyen et al., 2018; Du et al., 2015). In this study, we generated the training sample datasets separately from 0.15 % to 0.35 % of the total covered area for each land cover class (in 0.05 % increments). Then, sensitivity tests were carried out for the four land cover types. This means the training dataset would represent the values of fractional snow cover categories for each land cover type (see Sect. 4.1.2). All selection operations were completely random.

3.4 Description of different estimation methods

In this study, we compared the random forest method with the other three methods for retrieving fractional snow cover, including linear regression, ANN, and MARS. Note that the four methods input the same variables selected by the sensitivity test, including 12 characteristic variables and 1 target variable (see Sects. 3.3.1 and 4.1.1).

3.4.1 Linear regression

For optical remote sensing studies, there is a classical and general linear regression method used to estimate the sub-pixel snow cover area in a medium- to high-spatial-resolution image. This only involves the relationship between NDSI and fractional snow cover derived from high-resolution snow cover maps (Salomonson and Appel, 2004, 2006). This type of regression method has been applied in generating the standard MODIS fractional snow cover product Collection 5. Similarly, the multiple linear regression method was used as a reference method in this study to estimate fractional snow cover based on passive microwave data. The inputs were the same as the other three methods in this study. This method was undertaken in WEKA 3.8.3 and did not use any attribute selection method. In the Supplement, we present the linear regression formulas of fractional snow cover estimation for the four land cover types (Eq. S1 and Table S6).

3.4.2 ANN

ANN is a popular machine learning technique widely applied in remote sensing studies. Tedesco et al. (2004) developed an SWE and snow depth retrieval algorithm based on an ANN technique using passive microwave brightness temperature. Xiao et al. (2018) also used ANN to derive snow depth, and Kuter et al. (2018) and Czyzowska-Wisniewski et al. (2015) used ANN to retrieve fractional snow cover from MODIS data.

ANN consists of multiple layers: an input layer, one or more hidden layers, and an output layer (Hecht-Nielsen, 1992). The network with a multilayer perceptron can easily handle the nonlinear relationship between the input and output without any prior knowledge (Haykin, 2009). The inputs of each neuron were multiplied and summed by the connection weight. The output results were subsequently computed using a nonlinear logistic sigmoid transfer function. For numerical data, the transfer function in WEKA substitutes the pure linear unit function for the logistic sigmoid.

Aside from data preprocessing, a crucial step in this process is to design and optimize the ANN structure for improved estimation performance and good generalization capability (Kuter et al., 2018). Kuter et al. (2018) demonstrated that multidimensional function modeling can be successfully achieved with one hidden layer network. All parameters were set to the default with the exception of the learning rate which was optimized through a simple trial-and-error method. Based on the accuracy index and the modeling speed, Table S3 (Supplement) shows that a learning rate of 0.2 generated the best performance for the ANN retrieval model.

3.4.3 MARS

The MARS technique has been applied in a number of studies and in many fields, such as classification and mapping (Quirós et al., 2009), atmosphere correction (Kuter et al., 2015), pile drivability prediction (Zhang and Goh, 2016), and fractional snow cover estimation (Kuter et al., 2018). Unlike ANN, the modeling process of MARS is flexible and straightforward. Friedman (1991) first proposed a MARS technique that organizes a simple model for the complex and high-dimensional relationship between input variables and the target by smoothly connecting piecewise linear polynomials (known as basis functions, BFs). The ranges of the input variables were cut into a series of sub-ranges by knots; these were the connection points for two pieces of BFs. A simple BF format of MARS is expressed in Eq. (3), in which $max\left(\cdot \right)$ indicates that only positive parts were taken (otherwise, it was assigned a zero), and τ is a univariate knot.

$$\begin{array}{}\text{(3)}& {\displaystyle }max(\mathrm{0},x-\mathit{\tau })=\left\{\begin{array}{ll}x-\mathit{\tau },& \text{if }x>\mathit{\tau }\\ \mathrm{0},& \text{otherwise}\end{array}\right.\end{array}$$

The MARS method involves two phrases to establish a regression model: forward phase and backward phase. In the forward phase, BFs were generated using a stepwise search of all univariate candidate knots and all variable interactions. These adopted knots and their corresponding pair of BFs should produce the greatest decrease in residual error. The BFs were successively added to the model until it reached the maximum number of BFs, resulting in an over-fitted and complicated model. In the backward phase, the redundant BFs that contribute least to the model prediction are completely excluded from the regression model. These two phases are an iterative process (Kuter et al., 2018; Zhang and Goh, 2016).

Two important parameters of MARS determine the model “growing” and “pruning” processes: the maximum number of basis functions (max_BFs) and the maximum degree of interactions among input variables (max_INT) (Kuter et al., 2018). Kuter et al. (2018) reported that the increase in the structural complexity of the model does not significantly contribute to improving the performance of the MARS model. We conducted several tests to optimize the structure of MARS and found that more complex structures had a longer modeling time; however, it did not significantly improve model performance. Specifically, the modeling time of the complex structure (max_BFs = 100, max_INT = 2) was 4 times greater than the simple structure (max_BFs = 40, max_INT = 2) based on our analysis experiments. As such, the simple structure was chosen, as per Kuter et al. (2018). We implemented an open MARS MATLAB source code available from http://www.cs.rtu.lv/jekabsons/ (last access: 17 February 2021) for fractional snow cover estimation. These codes were compiled on a 2.40 GHz Intel Xeon central processing unit server.

3.4.4 Random forest

Random forest builds a large series of decision trees by applying the bootstrap sampling method. During the training stage, each tree grows by randomly selecting several variables and samples from input datasets (Mutanga et al., 2012). Input data were repeatedly split into training and test data using the bootstrapping method. Each randomly selected bootstrap sample in each iteration contained approximately two-thirds of the input elements. The remaining data, referred as out-of-bag (OOB) data, were used for validation. The predicted value of OOB data was produced from all the tree results that were generated, and the OOB error was subsequently calculated. For classification, the output was determined by voting for the results from all decision trees, whereas for regression, the output results were determined by averaging. The random forest was conducted in WEKA 3.8.3. As several attempts to optimize the parameters of random forest structure had failed, all parameters used were the default values.

3.5 Snow cover identification

The microwave radiation characteristics of precipitation, cold deserts, and frozen ground are similar to those of snow cover (Grody and Basist, 1996), and as such, snow cover area is likely to be overestimated. Grody and Basist (1996) proposed a snow cover identification algorithm distinguishing snow cover from precipitation, cold desert, and frozen ground. Many researchers have since used Grody's algorithm and its derivative algorithm to detect snow cover (Che et al., 2008; Xiao et al., 2018; Wang et al., 2019). Liu et al. (2018) reported on the assessment of different passive microwave snow cover detection algorithms and demonstrated that Grody's algorithm had a higher precision (positive predictive value) than that of other algorithms. We adopted the revised snow cover decision tree of Grody's algorithm (Table 2) as the highest frequency in this study was 91 GHz instead of the 81 GHz in the Special Sensor Microwave Images/Sounder (SSMIS) sensors (Che et al., 2008).

Table 2 The description of the revised Grody's algorithm. The unit is Kelvin (K).

There were two main objectives for using the revised Grody's algorithm (hereafter referred to as Grody's algorithm) in this study. The first was to compare the snow cover identification capability of the proposed fractional snow cover estimation algorithm with respect to ground snow depth measurements (see Sect. 4.4). The second purpose was to assess the effect of non-snow scatterers in estimating fractional snow cover due to this algorithm's special ability to distinguish the non-snow scatterer (i.e., precipitation, cold desert, and frozen ground). In both optical and microwave remote sensing research, the capability assessment of snow-free detection has been regularly neglected in most snow cover detection studies.

3.6 Validation of snow cover identification

When using the in situ snow depth (or SWE) measurements to quantitatively validate the accuracy of snow cover area data, converting snow depth into binary snow cover using an appropriate threshold is the first challenge. Many different depth thresholds have been suggested in previous studies, for instance 2 cm for 20 m spatial resolution (Gascoin et al., 2019); 0 cm (Parajka et al., 2012), 1 cm (Zhang et al., 2019), 3 cm (Hao et al., 2018), 4 cm (Huang et al., 2018; Wang et al., 2008), and 15 cm (Gascoin et al., 2015) for 500 m spatial resolution; 2.5 cm for 5 km spatial resolution (Hori et al., 2017); 3 cm (Xu et al., 2016) and 5 cm for 25 km spatial resolution (Liu et al., 2018); and 2 cm for 0.75^∘ grid resolution (Brown and Derksen, 2013). Due to these significant disagreements in the depth thresholds, Gascoin et al. (2019) conducted a sensitivity experiment that tested the agreement between in situ measurements and optical snow cover area products. The sensitivity of passive microwave snow cover identification results to snow depth at 6.25 km spatial resolution was also tested by computing the accuracy metrics with snow depth increasing from 0 to 10 cm.

Table 3 Confusion matrix defining the accuracy of the predicted snow cover map reference to the in situ snow cover observation. The characters (TP, FP, FN, and TN) represent the number of snow-covered or snow-free records in a particular condition.

Then, we needed to determine the threshold for converting fractional snow cover to binary snow cover. To date, few studies exist on fractional snow cover from the passive microwave pixel level. Dai et al. (2017) considered snow cover on the grid if fractional snow cover (25 km) was larger than 10 %. If the fraction of snow cover was less than 0.25, the SWE was set to 0 mm to correct the snow cover area in the daily SWE product based on Luojus et al. (2018). However, optical remote sensing studies often adopted 0.5 as the threshold of fractional snow cover (Hall and Riggs, 2007). Sensitivity experiments of fractional snow cover similar to snow depth were conducted to obtain the optimum conversion threshold. Both sensitivity experiments were carried out using 2017 bare-land-type datasets in Sect. 4.4.

3.7 Performance accuracy assessment

When evaluating the estimation performance of fractional snow cover in Sect. 4.1–4.3, we used conventional accuracy metrics: correlation coefficient (R; Eq. 4), mean absolute error (MAE; Eq. 5), and root mean squire error (RMSE; Eq. 6). Where $\stackrel{\mathrm{‾}}{x}$ is the mean value of all predicted values x_i, $\stackrel{\mathrm{‾}}{y}$ is the mean value of all target values, y_i, and n denotes the amount of data used.

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }R={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left(x_i-\stackrel{\mathrm{‾}}{x}\right)\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}{\sqrt{{\sum }_{i=\mathrm{1}}^n{\left(x_i-\stackrel{\mathrm{‾}}{x}\right)}^{\mathrm{2}}{\sum }_{i=\mathrm{1}}^n{\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}}}}\text{(5)}& {\displaystyle }& {\displaystyle }\text{MAE}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n\left|x_i-y_i\right|\text{(6)}& {\displaystyle }& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n{\left(x_i-y_i\right)}^{\mathrm{2}}}\end{array}$$

We evaluated the predicted accuracy of fractional snow cover and assessed snow cover identification performance (see Sect. 4.4). Six accuracy assessment indices were used for the analysis of snow cover detection capability (Liu et al., 2018; Gascoin et al., 2019): overall accuracy (OA), precision (that is, a positive prediction value), recall, specificity (that is, the true negative rate), F1 score (Zhong et al., 2019), and Cohen's kappa coefficient (Foody, 2020). OA refers to the proportion of correctly classified pixels as snow-covered and snow-free. The F1 score is a weighted average measurement of precision and recall ranging from 0 to 1 to measure the accuracy of binary classification. Cohen's kappa coefficient measures the agreement between snow cover products and ground measurements. All these indices were calculated from the confusion matrix (Table 3). OE is the omission error, and CE is the commission error.

4 Results analysis

4.1 Sensitivity in the training sample

4.1.1 Influence of input variables on model performance

We evaluated the results from 24 random forest fractional snow cover retrieval models (four types times six scenarios) to better understand which input variables have a good relationship with fractional snow cover and the combination of the variables that can improve retrieval model performance. The data used for variable sensitivity tests in this part spanned only 2 years (2014–2015) as the 91 GHz H polarization data were missing over the area south of 50^∘ N for 2016–2017. The OOB error and 10-fold cross-validation error were used to measure the performance of fractional snow cover retrieval models in each scenario (Mutanga et al., 2012). Table 4 shows the results of the six scenarios for the bare-land-type datasets.

Table 4 Variable selection tests in six scenarios on bare-land-type data for the random forest method. The accuracy indexes of the estimation are calculated using OOB error estimates and 10-fold cross-validation (CV).

The variable selection tests were used to identify a better combination of different variables (Table 4). Scenario 3, which only involves V polarization data, yielded the smallest R (0.590) and the largest MAE (0.197) and RMSE (0.248) of OOB error, as well as for the 10-fold cross-validation error (R: 0.596; MAE: 0.197; RMSE: 0.246). Scenario 3 performed the poorest of the six scenarios, which may be due to the lack of further available information from input variables that could be fully exploited (Xiao et al., 2018). Scenario 2, only containing passive microwave brightness temperature data similar to the variables used in Xu et al. (2016), had the second poorest performance. This shows that location information and topographic factors play a crucial role in snowpack distribution (Revuelto et al., 2014; Czyzowska-Wisniewski et al., 2015; Sturm et al., 2010). In this study, the retrieval method required these five basic input variables as auxiliary information in order to learn the characteristics of snow cover under different surface conditions to assist in accurately estimating snow cover properties. In contrast, in the absence of these basic input variables, the established model has no advantage in accurately predicting the characteristics of fractional snow cover under complex surface conditions. The major difference between scenarios 1, 4, 5, and 6 and scenarios 2 and 3 (Table 1) was whether or not the basic input variables (location information and topographic factors) were considered. The comparison results (i.e., scenarios 1, 4, 5, and 6 vs. scenarios 2 and 3) indicate that the effect of location information and topography need to be considered to estimate snow parameters. Moreover, when compared to scenarios 1, 4, and 5, the setting in Scenario 6, in which input variables were selected by importance, had the third poorest performance with a low R and a high MAE and RMSE. Scenarios 1, 4, and 5, generated better results; there were no significant differences in R, MAE, and RMSE for the tests on the four land cover types (Tables 4 and S4 and S5 in the Supplement). The comparison between scenarios 1, 4, and 5 indirectly indicates that the variables used in Scenario 1 may have some information redundancy that slightly weakens the efficiency of the random forest retrieval model. Although the selection methods of scenarios 4 and 5 performed well (in terms of modeling time and accuracy of predicted target), only one scenario was selected; the other may be used as an alternative in future. To this end, the variable combinations in Scenario 5 were selected for further analysis.

4.1.2 Determination of sample size

The datasets from 2014 to 2016 were used to examine the sensitivity to training sample size when the accuracy metrics used were the same as Sect. 4.1.1. As the values and variation trends of the accuracy metrics of the OOB error and 10-fold cross-validation error were almost equivalent, only the OOB error is shown in Figs. 4 and S1, S2, S3 (in the Supplement). We compared the performance of the random-forest-based models by altering the training sample size for the four land cover types.

Figure 4 The performance of random forest models when increasing the size of the training sample for shrub type.

Figure 4 illustrates that there is a slow increase in R and a slight decrease in MAE and RMSE when the training sample size increased from 0.15 % to 0.25 % on the shrub type, whilst there was a significant increase in modeling time. As the sample size increases from 0.25 % to 0.35 %, the model consistently estimates fractional snow cover accurately (higher R and lower MAE and RMSE). This finding appears to be consistent with previous studies (Colditz, 2015; Nguyen et al., 2018). An applicable and eligible sample selection scheme, which can achieve an acceptable target prediction accuracy and an adequate execution time, is essential for the implementation of a random forest model with superior predictive capability. One noticeable distinction between the three sample sizes (0.25 %–35 %) was modeling time. Interestingly, the 0.3 % training sample size had the shortest modeling time of the three sample sizes (Fig. 4); Figs. S1, S2, and S3 also exhibit similar findings on modeling time. The reasons underpinning the difference in modeling time is beyond the scope of this study and requires further research. We used the sample dataset covering 0.30 % of the study area of each class as a suitable size to randomly select training samples. We subsequently extracted the training samples for each land cover type from the 2011–2016 dataset to establish the retrieval models.

4.2 Comparison of the four retrieval methods

In this section, the independent testing datasets from 2010 were entirely used to assess the predictive performance of the random-forest-based models and the other three models (based on linear regression, ANN, and MARS). A comparison of the modeling time for the four methods (Table 5) showed that linear regression had the shortest time, with less than 1 s for the four land cover types, followed by ANN with approximately 51 s (forest), 22 s (shrub), 156 s (prairie), and 35 s (bare land). Random forest modeling times were very close to ANN modeling times for each land cover type. In contrast, MARS was the most time consuming with the longest time (approximately 6.5 h) for the prairie type and the shortest (19 min) for the shrub type. The absolute value of modeling time would vary under different computing capabilities.

Table 5 Performance of linear regression, ANN, MARS, and random forest model using test datasets from 2010 for four land cover types. FSC indicates fractional snow cover. The number outside parentheses indicates the number of pixels. The number inside brackets indicates their percentage.

Figure 5 The variation in the accuracy indexes (MAE, RMSE, and R) of the four algorithms (linear regression, ANN, MARS, and random forest) for four land covers.

Figure 6 The color density scatter plots between the estimated fractional snow cover and MODIS-derived fractional snow cover for four algorithms (linear regression, ANN, MARS, and random forest) for forest type. The accuracy metric refers to Table 5. (Note: out-of-range fractional snow cover values of linear regression, ANN, and MARS were truncated on 0 and 1). Note that all extracted records in January and February 2010 were used as the testing sample.

Table 5 and Fig. 5 present the results of the four retrieval methods for the four land cover types. The retrieval models of the shrub type predominantly have the lowest RMSE in contrast to the other three land cover types using the four methods (Table 5; cf. Figs. 5 and S4, S5, S6 in the Supplement). The random forest model had the highest R (0.916), lowest MAE (0.202) and RMSE (0.245), and no out-of-range records for the forest type (Table 5; Figs. 5 and 6). The distribution of and variation in MAE and RMSE for the four methods were similar under different land cover types with the exception of the shrub type (Table 5; Fig. 5). With the exception of ANN, the ranking of the three algorithms based on the accuracy of results (MAE and RMSE) for the shrub test data was also the same as under other land cover types (i.e., random forest > MARS > linear regression). Random forest had the greatest R value, followed by ANN, then MARS, and for most of the land cover types, the smallest R value was from the linear regression. Figure 6 illustrates that the random forest (Fig. 6d) produced a relatively small number of overestimated (∼0) and underestimated (∼1) values compared to the other three models (Fig. 6a–c). The MAE (0.315) and RMSE (0.401) of ANN were greater than MARS (MAE = 0.208, RMSE = 0.254). The number of out-of-range estimated values of ANN (36.62 %; 161 260) was also greater than MARS (2.65 %; 11 667), which may be attributed to a major underestimation of the fractional snow cover using the ANN method. The maximum and minimum values of ANN and MARS for the forest type were 0.949 (−0.52) and 2.132 (−0.122), respectively. For the other three land cover types, the number of out-of-range pixels of the four methods were almost in the same order (random forest < ANN < MARS < linear regression).

The random-forest-based models had the best performance with the highest R and lowest MAE and RMSE (Table 5). Previous studies have generally neglected the analysis and evaluation of whether the estimated value is out of range (Liang et al., 2017; Wang et al., 2017; Hao et al., 2019; Masson et al., 2018). Table 5 shows that the random forest models for the four land cover types produced reasonable fractional snow cover values ranging between 0 and 1. In comparison, the estimated fractional snow cover from the other three methods (linear regression, ANN, and MARS) was beyond this range. From the number of out-of-range records, the linear regression method generated the largest number of out-of-range fractional snow cover estimates with more than 0.85 million pixels (18.69 %). Although the number of out-of-range records of ANN (12.31 %) was less than MARS (16.39 %), both numbers exceed 0.5 million. Kuter et al. (2018) estimated fractional snow cover using the MARS and ANN techniques, which also yielded similar out-of-range values. The linear regression method had the poorest performance in estimating fractional snow cover from passive microwave data with the lowest R and the largest number of out-of-range records. These results indicate that nonlinear methods should first be used. Xiao et al. (2018) demonstrated the nonlinear relationship between passive microwave brightness and snow depth. De Lannoy et al. (2012) provided an exponential function for converting SWE to fractional snow cover. Thus, it is reasonable that a nonlinear relationship exists between fractional snow cover and passive microwave brightness temperature.

4.3 Evaluation of fractional snow cover

In this part, we analyzed the results from the training and evaluation stage for four land cover types (Table 6; Figs. 7 and 8). The independent data which were randomly selected from the data in January and February (2008, 2009, and 2017) and whose selection rule is the same as the training sample were used to further evaluate the predictive capability of random forest models in all range values. Other independent data in December (2007, 2008, 2016) were selected to examine the predictive capability of the established retrieval models in fractional snow cover to other months. To avoid the influence of the evolution of snow physical properties in the evaluation tests, we only considered December (Xiao et al., 2018).

Table 6 The performance of random forest models on training and evaluation data for four land cover types. Evaluation 1 is used to evaluate the estimation performance of the established retrieval models on fractional snow cover in January and February (Fig. 7). Evaluation 2 is used to analyze the prediction performance of the retrieval models in December (Fig. 8).

Figure 7 The color density scatter plots between the estimated fractional snow cover and MODIS-derived fractional snow cover in January and February for four land cover types (forest: a, b; shrub: c, d; prairie: e, f; bare land: g, h). The left column is the results of the training stage (a, c, e, g); the right column is the results of the evaluation stage (b, d, f, h). Statistics metrics refer to Table 6: Training and Evaluation 1. Note that 0.3 % of the independent evaluation dataset in January and February of 2008, 2009, and 2017 was randomly selected as the evaluation sample with stratified random sampling.

Figure 8 The color density scatter plots between the estimated fractional snow cover and MODIS-derived fractional snow cover in December for four land cover types: forest (a), shrub (b), prairie (c), and bare land (d). Statistics metrics refer to Table 6: Evaluation 2. Note that 0.3 % of the independent evaluation dataset in December of 2007, 2008, and 2016 was randomly selected as the evaluation sample with stratified random sampling.

We first saw the evaluation results in January and February (Evaluation 1 in Table 6). For forest type, Fig. 7a and b show that fractional snow cover around 1 is distinctly underestimated, and few are above the 1:1 line. The model for forest type had the poorest performance with the lowest R (0.636) and the highest RMSE (0.221) for the evaluation dataset (Table 6). The retrieval model on the prairie type had the best performance for the evaluation data (R: 0.752; MAE: 0.148; RMSE: 0.189). In shrub and bare land types (Fig. 7c, d, g, and h; Table 6), the retrieval models have a similar performance in the evaluation datasets (R: 0.712 and 0.719; MAE: 0.160 and 0.165; RMSE: 0.212 and 0.216, respectively). In the second evaluation experiments (Evaluation 2 in Table 6; Fig. 8), the best performance in predicted fractional snow cover is over prairie, and the relatively large underestimation can be seen over forest (MAE and RMSE). Meanwhile, we do not see striking differences between these two evaluation experiments (Evaluation 1 vs. Evaluation 2) with respect to their RMSE values. The difference in the used evaluation sample can explain the slight diversity in statistics metrics (R). Most of “true” fractional snow cover values in the training and validation datasets were distributed at two polar ends (0.0–0.3 and 0.9–1.0) in these two land cover types. When it comes to the results in the training stage and the evaluation stage, we found that the estimation performance of the retrieval model in the evaluation datasets is highly dependent on the quality of training sample which was used to establish the retrieval models. Figures 7 and 8 show that the established models in four land cover types can properly capture the characteristics of the entire range of fractional snow cover values.

Figure 9 Comparison of our estimated fractional snow cover (c; 6.25 km) with the reference MODIS fractional snow cover (b; 6.25 km) with respect to the MODIS composite binary snow cover products (a; 500 m) and a comparison example in the central Canadian area (d) on 27 February 2017 (2017058) (cf. the results in continuous value; Fig. S7 in the Supplement).

Apart from the scatter plots and statistical analysis, Fig. 9 shows the distribution pattern of snow cover from a spatial perspective on 27 February 2017, including MODIS composite binary snow cover (Fig. 9a), MODIS fractional snow cover (Fig. 9b), and the fractional snow cover estimated by the proposed algorithm (Fig. 9c). When the most rigorous pixel filtering rule at the 15×15 pixel window was applied (see Sect. 3.2), the large number of cloud-covered pixels (yellow) in Fig. 9a resulted in most areas of the MODIS fractional snow cover image (Fig. 9b) being represented by a “fill value”. Additionally, the number of intermediate values for MODIS fractional snow cover in winter would be much lower than the number of values near the two extreme values (0 and 1). In contrast, the estimated fractional snow cover from passive microwave brightness temperature data can provide almost complete coverage and continuous spatial information on snow cover (Figs. 9c and S7 in the Supplement). Figure 9 shows the comparison between our estimated fractional snow cover and the reference MODIS fractional snow cover and, more importantly, provides another perspective for snow cover identification in Sect. 4.4. Thus, Figs. 9b and 8c used the threshold 0.3 to define snow-covered and snow-free area, and this threshold was adopted through a series of experiments in Sect. 4.4. The threshold 0.5 was selected according to previous optical remote sensing studies on fractional snow cover (Tran et al., 2019; Marchane et al., 2015). This means that the pixel was identified as being snow-covered when the fractional snow cover value was greater than 0.3. From Fig. 9a–c, the spatial pattern of estimated fractional snow cover from the proposed method seems to accurately capture the distribution of snow cover from MODIS under clear-sky conditions, such as the snow-free area in most areas of North America and snow-covered areas in northern Canada. Figure 9d presents a specific example comparing these two fractional snow cover datasets and MODIS composite binary snow cover products in central Canada on 27 February 2017. We also provide another date comparison for MODIS and the estimated fractional snow cover on 10 January 2017 in the Supplement. Based on this example, we find that our estimated fractional snow cover was capable of obtaining the snow cover distribution when most of the area was covered by cloud, which was not the case for MODIS. This example also shows that the extent of snowline observed in the MODIS binary snow cover image (500 m), which was the boundary between snow-covered and snow-free areas, was well described and exhibited by the estimated fractional snow cover (6.25 km).

Thus far, we have evaluated the performance of random-forest-based models on independent datasets from 2008 to 2010 and 2017 for each land cover type. The results from the random forest (Table 5; Figs. 6d, S4–S6, and S7) show that the minimum estimates were higher than 0 and approximately 0.01. As to this overestimation, one possible error source is the inversion method used from which the final predicted outputs are obtained by averaging all results of the established sub-decision trees (see Sect. 3.4.4) (Breiman, 2001; Belgiu and Drăguţ, 2016), and this overestimation can be found in other studies applying the random forest method (Wei et al., 2019). Additionally, the attenuation and scattering of passive microwave radiation are not only caused by the snowpack. The non-snow scatterers (e.g., precipitation, cold desert, frozen ground) may be the majority error source potentially contributing to the overestimation of snow cover area as these scatters were easily misclassified as snow cover in less snow cover conditions (Grody and Basist, 1996; Dai et al., 2017). A detailed analysis on the misclassification due to non-snow scatterers is provided in Sect. 4.4. Although MODIS snow cover products are highly accurate in snow cover identification (Tran et al., 2019), the estimated results indicate that a large number of fractional snow cover values were overestimated (∼0) and underestimated (∼1). Some fractional snow cover estimates, at the individual pixel level, show a large discrete distribution near the 1:1 line (Fig. 7). These misestimates are not confined to the results of the random forest model but also appear in results of the other three methods (Figs. 6 and S4–S6 in the Supplement). Satellite sensors may provide completely different snow cover information because of different satellite overpass time. In this study, the difference in the Equator crossing time between MODIS and the passive microwave sensor was close to 6.5 and 9.5 h (see Sect. 2.1 and 2.3.1). Generally, the error caused by the differences in the satellite overpass time may easily be neglected when using multi-sensor observations for data fusion.

4.4 Validation using ground measurements

In winter with clouds and snow cover, the MCTD10A1 data still contained a large number of clouds (Fig. 9a; yellow) despite the implementation of the cloud removal and filling process for MODIS snow cover data. When applying the rigorous pixel filters (see Sect. 3.2), there were very little snow cover data for further model training and results analysis in one image (Fig. 9b). To evaluate and validate the estimated fractional snow cover in the absence of reference MODIS fractional snow cover, we conducted further analysis on snow cover detection capability. The ground snow depth measurements were utilized to investigate the accuracy of snow cover identification from two snow cover datasets: the snow cover detected by Grody's algorithm and fractional snow cover derived from random forest. We collected all available meteorological station snow depth measurements from 2008 to 2009 and 2017 (January and February) over North America, obtaining more than 900 000 pairs of records. This includes the snow depth measurements, snow cover area converted from the estimated fractional snow cover (hereafter referred to as random forest SCA), and snow cover area derived from Grody's snow cover mapping algorithm (hereafter referred to as Grody's algorithm SCA).

Figure 10 The changes in accuracy indicators (OA, precision, recall, specificity, F1 score, kappa) for snow cover detection are the result of two algorithms (a: Grody' algorithm; b: random forest) with increasing in situ snow depth value.

The sensitivity to ground-based snow depth in the snow cover detection results was tested by computing the accuracy metrics using data from 2017. Figure 10 shows that the accuracy metrics vary with increasing snow depth, whereby the metrics change significantly when snow depth exceeds 2 cmand reach a relative optimum when snow depth is equal to 2 cm. Che et al. (2008) stated that snow cover may be detected by passive microwave sensors when snow depth is greater than 2 cm. For this reason, we adopted 2 cm as the optimum depth threshold to transform ground snow depth measurements to snow-covered or snow-free information. We also conducted a series of sensitivity experiments to search for the optimum threshold for converting fractional snow cover to binary snow cover (Fig. 11). Figure 11 shows that recall and precision have opposing variation trends; the F1 score is up to the maximum value when FSC = 0.3 (FSC indicates fractional snow cover). In addition, the other two indicators (OA, kappa) also reached their maximum when the FSC values ranged between 0.3 and 0.4. As expected, 0.3 was used as the conversion threshold for fractional snow cover. Nevertheless, the conversion thresholds of snow depth and fractional snow cover need to be optimized with more data in the future.

Figure 11 The changes in accuracy indicators (OA, precision, recall, specificity, F1 score, kappa) for snow cover detection result in an increasing fractional snow cover value (FSC).

Figure 12 The accuracy indicators (OA, precision, recall, specificity, F1 score, kappa) of snow cover detection from two algorithm (Grody's algorithm; random forest).

We used a 2 cm snow depth threshold and a 0.3 fractional snow cover threshold to calculate the confusion matrix for Grody's algorithm SCA and random forest SCA against ground snow depth measurements in 2008–2009 and 2017 (Figs. 12 and S9). Figure 12 illustrates that the overall accuracy of snow cover identification had significantly improved by 20 %, from 0.71 for Grody's algorithm SCA to 0.85 for random forest SCA, indicating that the latter's results were in good agreement with ground snow cover measurements (kappa = 0.65). For the random forest SCA, precision (0.84) was lower than the recall (0.86), which means that snow cover area was more likely to be overestimated (CE = 0.16) than underestimated (OE = 0.14) with respect to in situ measurements. In contrast, for Grody's algorithm SCA, precision (0.87) was larger than the recall (0.51). By utilizing the proposed method, the OE of snow cover identification had been reduced by 71 % compared to the OE for Grody's algorithm SCA. The snow cover identification accuracy for the four land cover types is illustrated in Fig. S9 (in the Supplement) by radar charts.

Figure 13 The mixed snow cover detection map for different condition combinations of random forest SCA and Grody's algorithm SCA on 27 February 2017 (2017058). (a) SC_Grody=0 and FSC≤0.3; (b) SC_Grody=0 and FSC>0.3; (c) SC_Grody=1 and FSC≤0.3; and (d) SC_Grody=1 and FSC>0.3. SC_Grody=0 denotes snow-free (precipitation, cold desert, and frozen ground) determined by Grody's algorithm SCA; otherwise, it is snow-covered. FSC≤0.3 denotes snow-free cover detected by random forest SCA; otherwise, it is snow-covered.

We subsequently explored the influence of non-snow scatterers in estimating fractional snow cover. The CE of Grody's algorithm (CE = 0.13) was lower than that of the random forest SCA (CE = 0.16). Figure 12 shows that the overall snow-free identification capability of Grody's algorithm SCA (specificity = 0.91) was significantly superior to the random forest SCA (specificity = 0.81), which was also apparent for the four land cover types (Fig. S9). This may possibly be due to Grody's algorithm filtering out non-snow scatter signature (precipitation, cold desert, and frozen ground) (Grody and Basist, 1996). We counted the number of records in which a pixel had been detected as snow-free by the station and Grody's algorithm but was, however, considered snow-covered by random forest SCA. The records which were misclassified as snow-covered by random forest SCA, although they are non-snow scatter components (precipitation, cold desert, and frozen ground), account for 72.6 % of total misclassification records (CE = 0.16), of which 64.3 % come from precipitation, 7.0 % from cold desert, and 0.9 % from frozen ground. For forest, shrub, prairie, and bare land types, this misclassification proportion because of the non-snow scatters were 76.0 %, 92.3 %, 68.6 %, and 65.4 %, respectively. These results demonstrate that the non-snow scatterer is the major source of snow cover misclassification for random forest FSC results (Grody and Basist, 1996). Therefore, it is necessary to first distinguish the scattering signature of snow cover from other non-snow scattering signatures when using passive microwave data to identify snow cover. Similar preprocessing has been applied to snow depth estimation to minimize its uncertainties (Xiao et al., 2018; Wang et al., 2019; Tedesco and Jeyaratnam, 2016). An example is given to illustrate the inconsistency and consistency results of snow cover mapping between the random forest SCA and Grody's algorithm SCA (Fig. 13). Figure 13 shows that most areas of North America have consistent snow cover mapping results: Code A (cyan) and B (green) regions. Referring to the fractional snow cover image (Fig. S7b in the Supplement), the inconsistencies in monitoring snow cover extent always occur in mid-latitude regions with the value of fractional snow cover around 0.5. According to our experimental results analysis, the mid-value of fractional snow cover regions is one source of the omission error for Grody's algorithm SCA.

5 Discussions

5.1 Sensitivity to training sample size and quality

The size and quality of training samples may contribute to a large error at the individual pixel level (Dobreva and Klein, 2011; Kuter et al., 2018; Nguyen et al., 2018; Belgiu and Drăgu, 2016). Previous studies have investigated the sensitivity to sample size and sample quality (Nguyen et al., 2018; Colditz, 2015; Lyons et al., 2018). While some studies indicate that a larger training sample size improves the accuracy of estimates, we found that a training sample dataset covering about 0.3 % of the total study area was sufficient to achieve high accuracy in the estimation of fractional snow cover. When compared to previous sensitivity tests on sample size (Nguyen et al., 2018), the major difference was taking modeling time as an index in this study.

The estimation results of the random forest model (for the training, testing, and evaluation datasets; Sect. 4.2 and 4.3) showed that, generally, the prediction performance of the random forest model was closely related to the quality of the training sample. In this study's datasets, a greater number of records were located near the extreme values of the fractional snow cover (0 and 1). Thus, it is reasonable to use stratified random sampling (Dobreva and Klein, 2011) but not the proportional distribution of target values suggested by previous studies (Nguyen et al., 2018; Millard and Richardson, 2015). Even in these cases, the overestimation and underestimation often occur in the results of training datasets (Fig. 7a, c, e, and g) and evaluation datasets (Fig. 7b, d, f, and h), respectively. This is mainly because the established fractional snow cover retrieval model when using the training sample with relatively low diversity of fractional snow cover values does not learn the snow cover distribution characteristics of the various surface conditions well. Therefore, it is necessary for future studies to increase the number of samples by extending the study period to the snow accumulation and snow ablation stages (Xiao et al., 2018) where there is much more shallow snow and “patchy” snow cover. Another option is using data from multi-source sensors to generate reference snow cover data (e.g., Sentinel-1 Synthetic Aperture Radar data). By doing this, the proportion of fractional snow cover values in the training sample may be distributed as evenly as possible (Colditz, 2015; Jin et al., 2014; Lyons et al., 2018).

5.2 Effects of vegetation

Snow cover detection can be partially or completely obscured (or intercepted) by dense vegetation canopies. This introduces major uncertainties in the accurate detection of snow cover (Che et al., 2016; Hall et al., 2001; Parajka et al., 2012). Forest cover is an influential factor that cannot be ignored in optical and microwave remote sensing studies (Metsämäki et al., 2005; Cohen et al., 2015). It is evident that fractional snow cover retrievals typically are the least accurate under the forest type in comparison to other land cover types (Table 5; Figs. 6 and 7). There are two reasons that may be attributed to this error; one is the accuracy of the reference “true” fractional snow cover data in a forested area (Riggs and Hall, 2016), and the other is the microwave radiation attenuation caused by forests (Che et al., 2016).

Previous studies have reported that the lower accuracy of MODIS snow cover products was found in forest-covered areas and complex terrain (Hall and Riggs, 2007; Tran et al., 2019; Coll and Li, 2018). Several studies have validated and evaluated the accuracy of MODIS snow cover products, particularly in forested areas (Parajka et al., 2012; Zhang et al., 2019; Arsenault et al., 2014; Kostadinov and Lookingbill, 2015). In term of the NDSI threshold in forested areas (Sect. 3.2), we used 0.4 as a conservative threshold. According to previous studies, our operation (merely using NDSI as the criterion) in forest-covered areas would produce greater commission errors compared to using the Normalized Difference Vegetation Index (NDVI) as auxiliary information (Hall and Riggs, 2007). The retrieval results indicated that the NDSI threshold in forested areas needs optimization using numerous datasets (Riggs et al., 2017; Xin et al., 2012). In addition to the influence of forests on MODIS data, forests also hamper the upwelling microwave radiation emitted from the ground. Snow cover in forested areas may be divided into under and above forest canopy snow cover (Xin et al., 2012). This apparently distinguishes the interference effects of evergreen forests and deciduous forests on snow cover (Gascoin et al., 2019; Romanov and Tarpley, 2007). Additionally, there are major differences in the observation means for optical and passive microwave sensors in forested areas. The capacity for optical sensors to observe above forest canopy snow cover is mainly dependent on the vegetation canopy density (Kuter et al., 2018), while microwave sensor may obtain information on snow cover under the vegetation canopy (under forest canopy snow cover) (Che et al., 2016; Cohen et al., 2015). Overall, the combination of these two effects may produce low estimation accuracy for fractional snow cover.

6 Conclusions

Many previous studies have focused on estimating fractional snow cover utilizing optical remote sensing imagery which suffers from cloud contamination during data acquisition. In contrast, microwave sensors offer the attractive advantage of working under all weather conditions. In this study, we attempted to develop an algorithm for applying the enhanced-resolution passive microwave brightness temperature data (6.25 km) to fractional snow cover estimation during January and February of 2010 to 2017. Using the reference fractional snow cover stem from MODIS snow cover products as the “true” observation, we established fractional snow cover retrieval models for four land cover types (forest, shrub, prairie, and bare land) inputting 12 variables selected by 24 sensitivity experiments. The proposed algorithm accounted for a series of influential factors including topography, land cover, and location information. Compared to the other three methods (linear regression, ANN, and MARS), the random-forest-based algorithm had the best performance with high accuracy (highest R and lowest MAE and RMSE) and no out-of-range retrievals. The results of the evaluation using the reference fractional snow cover data in 2017 showed that our proposed algorithm had good retrieval performance in estimating fractional snow cover, with RMSEs ranging from 0.167 to 0.198. Moreover, in situ snow depth measurements were used to validate the accuracy of the proposed fractional snow cover estimation algorithm in snow cover mapping. The snow cover detection capability of the random-forest-based method was superior (OA = 0.85, kappa = 0.65) to that of Grody's algorithm, with overall accuracy increasing by 20 % (from 0.71 to 0.85) and omission error reducing by 71 % (from 0.48 to 0.14) when the fractional snow cover threshold was 0.3. Although the random-forest-based models achieved an acceptable accuracy, the fractional snow cover was more likely to be overestimated (CE = 0.16) than underestimated (OE = 0.14). In addition, the effect of the non-snow scatterer was evaluated on fractional snow cover predictions by means of the good prediction of Grody's algorithm on snow-free class; the results indicated that more than 70 % of CE was caused by misclassifying the non-snow scatterer (precipitation, cold desert, frozen ground) as snow cover. These models established that using several data sources in January and February had better applicability in dry snow conditions, while estimation results could be less accurate in wet snow conditions.

Numerous studies have investigated the relationship between common snowpack physical properties (e.g., snow depth and water equivalent) and passive microwave brightness temperature at different frequencies and polarizations (Chang et al., 1987; Dietz et al., 2011; Kim et al., 2019; Xiao et al., 2018). Unlike many previous studies, this study innovatively used passive microwave data to directly estimate fractional snow cover. The results showed that it is possible to directly obtain an estimated fractional snow cover with high accuracy from high-spatial-resolution passive microwave data (6.25 km) under all weather conditions. Further detailed study on the use of high-spatial-resolution passive microwave data for fractional snow cover estimation presents itself as an interesting research direction for the development of the studies on fractional snow cover estimation. Furthermore, to reduce some of the limitations (e.g., forest effects) (Cohen et al., 2015) and deficiencies (overestimation and underestimation) identified in this study, future studies should pay greater attention to the prediction of the fractional snow cover using passive microwave data. To the best of our knowledge, this study may represent the first attempt to establish a relationship between microwave brightness temperature and the reference “true” fractional snow cover using machine learning methods. However, it also contains significant limitations in understanding the physics that relates fractional snow cover to the signature of passive microwave brightness temperature (Cohen et al., 2015; Che et al., 2016). Future studies need to use physical snowpack models and radiation transfer theory to explore the physical mechanistic relationships between microwave brightness temperature and fractional snow cover (Pan et al., 2014).