Evaluation of VIIRS and MODIS Snow Cover Fraction in High-Mountain Asia Using Landsat 8 OLI

We present the first application of the Snow Covered Area and Grain size model (SCAG) to the Visible Infrared imaging Radiometer Suite (VIIRS) and assess these retrievals with finer‐resolution fractional snow cover maps from Landsat 8 Operational Land Imager (OLI). Because Landsat 8 OLI avoids saturation issues common to Landsat 1–7 in the visible wavelengths, we re-assess the accuracy of the SCAG fractional snow cover maps from Moderate Resolution Imaging Spectroradiometer (MODIS) that were previously evaluated using data from earlier Landsat sensors. Use of the fractional snow cover maps from Landsat 8 OLI shows a negative bias of −0.5% for MODSCAG and −1.3% for VIIRSCAG, whereas previous MODSCAG evaluations found a bias of −7.6% in the Himalaya. We find similar root mean squared error (RMSE) values of 0.133 and 0.125 for MODIS and VIIRS, respectively. The Recall statistic (probability of detection) for cells with more than 15% snow cover in this challenging steep topography was found to be 0.90 for both MODSCAG and VIIRSCAG, significantly higher than previous evaluations based on Landsat 5 Thematic Mapper (TM) and 7 Enhanced Thematic Mapper Plus (ETM+). In addition, daily retrievals from MODIS and VIIRS are consistent across gradients of elevation, slope, and aspect. Different native resolutions of the gridded products at 1 km and 500 m for VIIRS and MODIS, respectively, result in snow cover maps showing a slightly different distribution of values with VIIRS having more mixed pixels and MODIS having 7% more pure snow pixels. Despite the resolution differences, the snow maps from both sensors produce similar total snow-covered areas and snow-line elevations in this region, with R 2 values of 0.98 and 0.88, respectively. We find that the SCAG algorithm performs consistently across various spatial resolutions and that fractional snow cover maps from the VIIRS instruments aboard Suomi NPP, JPPS–1, and JPPS–2 can be a suitable replacement as MODIS sensors reach their ends of life.


INTRODUCTION
Earth's cryosphere is rapidly diminishing in extent and volume (Bormann et al., 2018). Seasonal duration of the snow cover has become shorter (Dettinger and Cayan, 1995;Derksen and Brown, 2012), and earlier onset of snowmelt leads to earlier streamflow in snow-dominated basins (Stewart et al., 2005). Global climate models show that more countries will experience water stress in the future, particularly those that rely on melt of snow and ice for their water resources (Arnell, 1999). Observational datasets of glaciers show mass loss on a global scale (Zemp et al., 2015). In the mid-latitudes especially, glaciers ablate faster when the winter snow cover disappears earlier . Interannual variability and trends in snow cover result from variability in snowfall and the rate of snowmelt, which is driven by the energy balance at the snow surface and is especially sensitive to snow albedo. A comprehensive view of changes and variability in snow extent and albedo can be observed directly only with satellite remote sensing.
Spaceborne sensors observing in the visible and near-infrared parts of the spectrum have long been popular instruments for detection and monitoring of snow cover (Dozier, 1989;Hall et al., 2002;Dietz et al., 2012). For clean snow, the spectral reflectance can be modeled from its grain size and the wavelength-dependent complex index of refraction of ice (Wiscombe and Warren, 1980;Dozier and Painter, 2004). For both ice and water, absorption through the solar spectrum varies by eight orders of magnitude, from highly transparent in the visible to nearly opaque in the shortwave-infrared (Warren and Brandt, 2008). Thus, snow reflects nearly all radiation in the visible spectrum but absorbs strongly in the longer wavelengths (Figure 1), making it distinguishable from most other natural substances. When dust and industrial pollutants are present in snow, its reflectance drops in the visible and part of the near-infrared spectrum (400 -876 nm), where 63% of the at-surface solar irradiance occurs (Painter et al., 2012b;Bair et al., 2019).
The Moderate-resolution Imaging Spectroradiometer (MODIS, Justice et al., 1998) and the Visible Infrared Imaging Spectrometer Suite (VIIRS, Justice et al., 2013) collect data at spatial resolutions of 250 m to 1 km depending on the instrument band. The gridded Level 2 surface reflectance products are provided at 500 m for MODIS and 1 km for VIIRS. At these moderate spatial resolutions, with a daily return interval and global coverage, these sensors have been instrumental in capturing heterogeneous snow cover at the global scale both in complex terrain and at regional to watershed scales . The MODIS and VIIRS surface reflectance data are also well suited to examine snow cover surface properties in addition to snow cover, such as grain size and light-absorbing impurities (Painter et al., 2012b) that collectively drive the variability in snow albedo.
The snow mapping technique that is applied most commonly to data from multispectral satellites and used worldwide to map snow for products like MOD10A1 (Hall and Riggs, 2007) and the VIIRS Snow Cover product (Justice et al., 2013;Key et al., 2013) uses a computationally simple algorithm developed in the 1980s, when computers were more primitive. This approach, the Normalized Difference Snow Index (NDSI), was developed using band difference ratios (Dozier, 1989) and exploits the fact that snow is bright in the visible and dark in the shortwave-infrared wavelengths. The NDSI algorithm is effective in identifying snow in homogeneous areas where the presence of mixed pixels is low and vegetation is sparse, or at the finest spatial resolutions where most pixels are either entirely snow-covered or bare. However, for snow in rugged terrain, mixed pixels predominate at the MODIS resolution (Selkowitz et al., 2014); and as such, the performance of NDSI is degraded. While the MOD10A1 algorithm produces fractional snow cover (Salomonson and Appel, 2004), the value is from a regression of MODIS NDSI against binary Landsat snow cover maps based on NDSI. Specifically, the 30 m pixels from Landsat are classified as entirely snow-covered or snow-free before aggregation to 500 m for the regression analysis as a target variable. The binary classification of the Landsat data results in a systematic underestimate for snowcovered pixels with less than 50% snow cover and an overestimate for pixels with more than 50% snow cover for MODIS data (Figure 2, in Rittger et al., 2013).
To improve snow mapping for these mixed pixels, the unique spectra of snow, vegetation, and soil/rock combined with multispectral measurements led to the development of spectral un-mixing algorithms that calculate the snow cover fraction (f SCA ) of each pixel (Nolin et al., 1993;Rosenthal and Dozier, 1996;Vikhamar and Solberg, 2003;Sirguey et al., 2009;Bair et al., 2020). Spectral unmixing solves simultaneous equations to produce not only f SCA , but also snow surface properties such as grain size and broadband albedo (Painter et al., 2003;Painter et al., 2009). Such algorithms can be applied to reflectance data from any multispectral sensor that includes spectral bands from the visible through shortwave-infrared, and improvements in computational efficiency now enable subpixel snow retrievals from moderate-resolution imagery at the global scale.
While the previous spectral un-mixing algorithms take advantage of all available bands, the algorithm selected by the National Oceanic and Atmospheric Administration (NOAA) for fractional snow cover retrieval was developed for geostationary satellites and uses only visible reflectance (Romanov et al., 2003). NASA's standard snow cover product from VIIRS relies on the NDSI based method.
The MODIS Snow-Covered Area and Grain-Size algorithm (Painter et al., 2009) produced operationally for the Western United States, the Canadian Rockies and Alaska, High-mountain Asia, the North Pole, and South America demonstrates clear benefits over the standard global operational MODIS Snow Cover Products in mixed pixels and during accumulation and melt . MODSCAG data have been used for estimating sub-canopy snow cover (Raleigh et al., 2013;Rittger et al., 2020), light absorbing impurities on snow (Painter et al., 2012b), snow albedo (Bair et al., 2019), persistent ice cover (Painter et al., 2012a), impacts of wildfire on snowmelt (Micheletty et al., 2014), and snow water equivalent (SWE) (Guan et al., 2013;Bair et al., 2016;Rittger et al., 2016). It has also been used for evaluating climate models (Wrzesien et al., 2015;Minder et al., 2016); assimilating data into regional-scale climate models (Oaida et al., 2019); and validating global to regional scale models to better understand regional dust and carbon impacts on snowmelt (Sarangi et al., 2019;Sarangi et al., 2020).
In this study, we use high-resolution fractional snow cover maps, derived from the Landsat 8 Operational Land Imager (OLI) using spectral mixture analysis, to 1) re-evaluate MODSCAG fractional snow cover maps and 2) provide the first evaluation of VIIRSCAG maps over the challenging terrain in High Mountain Asia. The evaluation was conducted using 30 Landsat 8 OLI scenes, spatially spanning two MODIS/VIIRS tiles in the sinusoidal grid in the Himalaya. For the evaluation, we use binary metrics precision, recall, and F-score and fractional metrics mean difference, median difference, and RMSE. Landsat 8 OLI images in 2013 were selected late in the melt and early in the accumulation period, previously shown to be the most challenging times for accurate snow detection. In addition to the evaluation with Landsat 8 OLI, we directly compare snow cover area and snowline elevation from MODSCAG and VIIRSCAG for a more temporally complete dataset over the years of 2012-2015.

DATA AND METHODS
The initial validation of MODSCAG and its comparison to widely used snow cover maps from MOD10A1  used Landsat 5 Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper Plus (ETM+) images from Colorado, California, and the Himalaya to understand errors in retrievals from various snow regimes. Rittger et al. (2013) found that MODSCAG RMSE was highest and the probability of detection, also known as recall, was lowest in the Himalaya, which suggested snow identification was more difficult in that region, possibly because of steeper terrain relative to other validation sites in Colorado and California. In addition, this previous work found the highest errors in all regions early in the accumulation season, important for recreation such as skiing, and late in the melt season, when river forecasts are often in error (Dozier, 2011).
For this study, we considered the spatial extent of the Himalayan region to include parts of the Amu Darya, Indus, and Ganges River basins shown in Figure 2. We searched the USGS database of Landsat 8 OLI imagery in 2013 for all images with cloud cover less than 10% in april to July (melt season) and September-November (early accumulation season). We then visually inspected all images and removed those with any cloud cover. The result was a set of 30 images from 17 unique Landsat Worldwide Reference System 2 (WRS2) path/rows covering diverse regions of the Himalaya. This list of images along with the scene-specific validation statistics listed in the previous section is detailed in Supplementary Table S1 of the Supplementary Materials. We selected MODIS and VIIRS images within ± 3 days with the smallest satellite view zenith angles corresponding to nadir viewing Landsat 8 OLI images.
The terrain of this study area provides a challenging environment for both MODIS and VIIRS, possibly more challenging for VIIRS because of the coarser spatial resolution relative to MODIS. The climatic characteristics of this domain vary from west to east with westerlies driving much of the FIGURE 1 | Spectral sampling differences between OLI, MODIS, and VIIRS shown as sensor bandpasses plotted with snow spectra for fine and coarse snow with and without light absorbing particles.  (Cao et al., 2013). With near-global daily coverage, moderate resolution retrievals at 11 spectral bands (Figure 1), and a 13:30 overpass time, Suomi NPP VIIRS provides similar information to that from MODIS for snow cover retrievals. While the daily Aqua overpass aligns temporally with Suomi NPP overpass, the mid-morning MODIS retrievals from Terra are generally preferred for snow retrieval due to cloud-cover considerations. The MODIS instrument aboard Aqua suffered from loss of 75% of its detectors for Band 6 shortly after launch (Gladkova et al., 2012). MODIS Band 6 observes radiances in the spectral range of 1.628-1.652 µm, which is a key bandwidth for common snow retrieval algorithms as shown in Figure 1. Robust discrimination between cloud and snow cover can be challenging (Dietz et al., 2012;Stillinger et al., 2019;White et al., 2020); and as such, we have visually selected clear-sky scenes for this study.
Pixel distortion created by the combination of MODIS scan geometry and the earth's curvature causes edge-of-scan pixels in MODIS surface reflectance data to be collected from an area nearly 10 times the size of those viewed from nadir (Dozier et al., 2008). Therefore, off-nadir observations image different (larger) areas on the ground than nadir views especially in topographically complex terrain where the view azimuth offsets the direction of pixel elongation . To address this shortcoming, the VIIRS instrument uses deaggregation of sub-pixels at wide scan angles (>32°, Cao et al., 2013) to sample at a more constant spatial resolution in the cross-track direction. This technique reduces edge-of-scan distortion in VIIRS surface reflectance data (Schueler et al., 2013) and, in turn, distortion in representation of perceived snow cover area at the edge of the swath. Because future JPSS missions plan to carry VIIRS optical instrumentation and MODIS's end-of-life is approaching, the future of continued global fractional snow cover observations will hinge on VIIRSbased retrievals.
For this study, the MODIS inputs used were MOD09GA L2G Surface Reflectance 500 m SIN grid from Collection 5 that are provided by the Land Processes DAAC (lpdaac.usgs.gov). We used Collection 5 instead of the more recent Collection 6 MODIS products to allow direct comparison with previous evaluations of MODSCAG . Degradation of visible and near-infrared bands in Collection 5 has been shown to impact long-term trends in geophysical variables, yet the impact is small in an absolute sense (Table 1, in Lyapustin et al., 2014) and below the stated accuracy of MODIS products. The errors in this paper for Collection 5 can be considered a high estimate relative to those if we had used Collection 6. The VIIRS inputs used were the early release products DSRF1KD L2GD Surface Reflectance 1 km SIN grid and DGA1KD L2GD Geolocation Angles 1 km SIN grid from Collection 3110, provided by the Level-1 and Atmosphere Archive and Distribution System (LAADS). While the VIIRS moderate resolution data are being resolved to 750 m in the native grid projection, the Level 2 gridded products are distributed at 1 km spatial resolution on the MODIS Sinusoidal Tile Grid for consistency with the 500 m MODIS archive. Consequently, the spatial resolution of our snow cover products derived from VIIRS are at 1 km. These surface reflectance products are likely consistent with those used to develop the snow retrieval algorithms for the Standard VIIRS Snow Products, (Key et al., 2013), although these inaugural studies do not explicitly state the data products used. However, since our analysis, Collection 3110 has been replaced by the newer VNP09GA Version 1 Surface Reflectance products. A comparison of the DSRF1KD and VNP09GA surface reflectance products shows minor differences in spectral reflectance over snow (Supplementary Figure S1).

Snow Covered Area and Grain Size
Fractional snow cover and grain size maps are derived from MODIS and VIIRS surface reflectance inputs using the Snow Cover and Grain Size (SCAG) algorithm (Painter et al., 2003;Painter et al., 2009). The algorithm disaggregates pixelaveraged surface reflectance into a group of individual land surface components using a combination of reflectance spectra from pure spectral end member libraries. Commonly known as spectral mixture analysis, the SCAG algorithm is based on a set of simultaneous linear equations that are solved for these individual components (Roberts et al., 1998).
F i is the fraction of endmember i; R λ,i is the hemispherical directional reflectance factor (HDRF, Schaepman-Strub et al., 2006) of endmember i at wavelength λ; and ε λ is a residual error. We assume a linear combination of reflectances using this method, which is appropriate above timberline where the surface is planar. Below timberline, there are nonlinear interactions with the canopy; however, we include canopylevel vegetation endmembers in the SCAG spectral library that represent multiple scattering and allow the linear assumption. The spectral library contains 109 snow (of varying grain size), 3 ice, 20 vegetation, and 17 rock spectra that were measured in situ (Painter et al., 2003;Painter et al., 2009) and allows for changes in snow spectra with viewing angle. The snow spectral HDRF come from radiative transfer modeling, constrained by snow grain size (10-1100 μm by 10 μm) and solar zenith angle by 15 degree increments (Painter et al., 2003;Painter et al., 2009). The SCAG model is run iteratively; and with each set of F i and R λ,i , we estimate the RMSE as the difference between the modeled and observed pixel-averaged reflectance. SCAG chooses the model with the lowest RMSE that provides estimates of the snow cover fraction, vegetation fraction, and bare soil or rock fraction. Fractional snow cover and other endmembers are normalized by the spectral fraction of photometric shade to account for topographic effects on irradiance. These model runs are constrained to endmember fractions in the range of −0.01-1.01 with an overall RMSE <2.5%, with no three spectrally consecutive residuals exceeding 2.5%. For pixels that do not meet these constraints, a second set of model runs is completed using looser constraints. SCAG selects the model with the fewest endmembers and the smallest RMSE. The selected endmember for snow (where it exists) has an associated grain size in the spectral library, which allows for snow grain size mapping. Whereas this grain size can be converted to a clean snow albedo or combined with impurities to estimate the true snow albedo (Wiscombe and Warren, 1980;Painter et al., 2009), this step is outside the scope of the work presented here. Historically, albedo from snow grain sizes derived with SCAG have been combined with dust impacts on albedo estimated from a separate algorithm (Painter et al., 2012b), however, recent work shows that these variables can be solved for simultaneously (Bair et al., 2020).
The SCAG products used in this study are the fractional snow cover area (f SCA ) maps. To account for the spectral differences between MODIS and VIIRS (Figure 1), the high resolution spectral libraries  were resampled to the VIIRS spectral bands using the relative spectral response data available from the Center for Satellite Applications and Research JPSS website.

Landsat TM, ETM+, and OLI
We use surface reflectance from USGS Landsat 8 OLI in this study. Landsat 5 TM and Landsat 7 ETM+ have been used to evaluate snow cover products in the past, including those from MODIS (Salomonson and Appel, 2004;Dozier et al., 2008;Painter et al., 2009;Bormann et al., 2012;Rittger et al., 2013;Metsämäki et al., 2015;Yang et al., 2015). TM and ETM+ sensors were designed for observing fine radiometric detail in dark ocean water and when compressed into an 8-bit data stream, result in saturation of the sensor over snow in bands 1, 2, and 3 (Dozier, 1989;Duguay and Ledrew, 1992). The more snow cover, the brighter the surface at visible wavelengths, and the more likely the sensor will become saturated. When saturation of the sensor occurs, the observed surface reflectance values are lower than the true reflectance. A complicating factor is the use of cubic convolution instead of nearest neighbor processing the USGS uses to resample Landsat data that can merge values from saturated pixels and unsaturated pixels making it more difficult to identify affected pixels.
For spectral mixture analysis, the use of the saturated values would introduce erroneously low reflectance values leading to an underestimation of f SCA . Therefore, when only one or two of the visible bands are saturated, spectral mixture analysis of TM and ETM+ data uses the remainder of the bands to produce the f SCA model output. When all three visible bands are saturated there is insufficient information to apply spectral mixture analysis. Previous work assumed that pixels saturated in all three visible bands were fully snow cover (Painter et al., 2009;Rittger et al., 2013), however, it is possible that they were not.
For this reason, the Landsat 5 TM and Landsat 7 ETM+ are imperfect sources of validation data leading to a systematic positive bias in f SCA from Landsat. The previous evaluation of MODSCAG using these Landsat f SCA estimates as "truth" found a mean bias of −5.2% over four regions and −7.6% in the Himalayan region . In this paper, we assess the number of pixels with saturation in previous comparisons to MODIS snow maps. NDSI and regressionbased approaches (Hall et al., 2002;Salomonson and Appel, 2004;2006) that assume accurate Landsat snow cover maps for validation and calibration do not appear to account for this issue, and the standard NASA global-fractional snow products still use regression equations based on saturated sensors.
In contrast to the 8-bit data stream from TM and ETM+, the optical sensor aboard Landsat 8 OLI, uses a 12-bit data stream that allows it to see the darkest ocean to the brightest snow without saturation, allowing the spectral mixture analysis to use all bands for even the brightest pixels. The higher accuracy of Landsat 8 data provide the opportunity to update the evaluation of MODSCAG fractional snow cover maps and contrast it with the first evaluation of VIIRSCAG performance.

Evaluation Metrics
We use a set of binary metrics previously used to evaluate MODIS snow cover maps  to facilitate comparison to the previous studies. Assuming Landsat 8 snow cover maps to be the "truth", we classify pixels into four possible outcomes: true positive (TP), true negative (TN), false positive (FP), and false negative (FN). We calculate the four outcomes for both MODIS and VIIRS snow maps and consider pixels <15% to be snow-free, again for consistency with previous studies. Common measures of a classification algorithm's performance are calculated using the four outcomes: We abstain from using the accuracy calculation (Eq. 4) because of the overwhelming number of TN pixels included when analyzing full Landsat scenes, which contain large snow-free areas. We find the F score to be the single most useful metric as it balances Precision and Recall, penalizing both missing snow and falsely identified bare areas as snow, while excluding extensive snow-free areas (Wilks, 2006). We use mean difference, median difference, and RMSE to evaluate performance of the MODIS and VIIRS fractional snow cover maps. We exclude TN occurrences from the RMSE calculation because large snow-free areas, especially in the summer, can result in very small errors . Previous studies (Hall et al., 2002;Appel, 2004, 2006;Painter et al., 2009) did include those in the error statistics, therefore yielding lower RMSE values when compared to this study. The RMSE is defined as: where f C SCA is the fractional snow cover from MODIS and VIIRS products and f T SCA is the fractional snow cover from the Landsat 8 OLI validation data.
We use mean absolute error (MAE) to quantify the discrepancy between the snow fraction for MODIS f M SCA , and that of VIIRS f V SCA : We use Landsat 8 OLI data employing the SCAG algorithm to compare with the coarser maps from both MODIS and VIIRS. After mapping fractional snow cover at 30 m spatial resolution, we coarsen the data by identifying the subset of Landsat pixels within each MODIS or VIIRS pixel (∼256 values) and then averaging those values. Similar to previous MODSCAG validation , we find that coarsening MODIS from 500 m to 2 km significantly improves both the binary and fractional statistics, indicating a MODIS geolocation error between 1 and 2 pixels as found in other studies (Wolfe et al., 2002). Similarly, we find the same reduction in errors when coarsening the VIIRS data from 1-4 km. We isolate the errors from the SCAG algorithm by performing the analysis at these coarser spatial resolutions.

RESULTS
In this section we briefly discuss the impact of saturation on 30 m snow maps from from Landsat 5 TM and Landsat 7 ETM+. We  then compare the binary and fractional statistics for snow maps in the melt and early accumulation periods from MODIS and VIIRS with that of Landsat 8 OLI and discuss the improved snow detection capability relative to Landsat 5 TM and Landsat 7 ETM+. Rittger et al. (2013) used 25 Himalayan images for validation of MODSCAG with three images from September, fifteen from October, and seven from November. The September through October dates occurred after the equinox when the Sun is lower in the sky and saturation is less likely because of lower incoming shortwave radiation. Figure 3 shows the percent of snow-covered pixels saturated in three band combinations for the 25 images. 28% of the pixels were saturated in all visible bands (1, 2, & 3), and were assumed to be fully snow covered as previously described above in Landsat TM and ETM+. Individually, band 1 saturates most frequently followed by band 3 and then band 2. Overall, 2% of the snow-covered pixels are saturated in band 1 only (not 2 & 3), and 4% of the pixels are saturated in bands 1 & 3 (not 2). This new study uses Landsat 8 OLI which does not saturate, and includes May (n 9), June (n 5), September (n 5), and October (n 11) of 2013.

Evaluation Metrics
The first column in Figure 4 shows the false color Landsat 8 OLI image at 30 m resolution for five dates in the Pamir-Alay mountain system, which runs through Tajikistan, Kyrgyzstan, and Uzbekistan. The Landsat WRS-2 path/row is 153/033, i.e. the Western most region in this study (Figure 2). This region is relatively cloud-free compared to regions to the East and provides the most cloud free images in this study. The second column in Figure 4 shows the fractional snow cover maps for the five image dates. These 30 m snow maps show that most pixels are nearly fully snow covered (white) while there are fewer partially covered snow pixels (blue). When coarsened to 500 m shown in the third column in Figure 4, the fractional snow cover increases as fully snow-covered pixels are aggregated with snow-free pixels, especially near the snowline. The fourth and fifth column in Figure 4, show the fractional snow cover maps at native resolution from MODIS and VIIRS, respectively. Visually, images appear similar from all three sensors indicating agreement.
The binary statistics used for comparison of the Landsat 8 OLI snow maps with the snow maps from MODIS and VIIRS are defined in Eqs. 2, 3, 5. Table 1 shows the summary of binary statistics which include the mean, minimum, maximum, and standard deviation for the 30 image comparisons as well as  Supplementary Table S1.
The maximum snow cover in the Hindu-Kush Himalayan region occurs in February to March (Figure 3, Armstrong et al., 2019), while minimum snow and ice cover usually occurs in October (Painter et al., 2012a). Therefore, our study covers a period during melt (April-June) and during the minimum and beginning of accumulation (September-November). Figure 5 shows the binary statistics within these time periods and are labeled as Melt Season and Accumulation Season.
Precision, which indicates the number of true positives relative to all positives was 0.871 and 0.855 for MODIS and VIIRS, respectively. For MODIS, precision was lower than previously reported; but previous work included midwinter periods where there is much more snow and high values for precision are easier to achieve. Figures 5A,B show the precision during the Melt Season and Accumulation Season, respectively. During Melt Season, a high precision is maintained, while during Accumulation Season there is more scatter with values as low as ∼0.6. Right before accumulation, at the minimum snow cover, precision is a bit higher.
Recall (i.e., the probability of detection) for MODIS and VIIRS was 0.900 and 0.889, respectively. This is considerably higher than previous estimates for MODIS, and is plotted for the Melt Season and Accumulation Season in Figures 5C,D, respectively. Recall during melt is very high for both MODIS and VIIRS except for one outlier on May 31st with a recall of 0.728, which is possibly due to unwanted cloud interference given the otherwise high accuracy across other dates. Similar to precision, recall is lower during early accumulation.
The F-score, which balances precision and recall, was 0.890 and 0.867 for MODIS and VIIRS, respectively. This value was nearly the same in the previous study for MODIS. As shown in Figures 5E,F, the F-score is higher in the Melt Season than the Accumulation Season as also seen in previous work. Compared to previous work, precision is lower but recall is higher, resulting in a similar overall F-score as assessed by binary statistics. Table 2 shows the summary (mean, minimum, maximum, and standard deviation) of fractional statistics described in previously in Evaluation Metrics. These fractional statistics on each date for both MODIS and VIIRS are included in Supplementary Table S1. The fractional metrics for the Melt Season period are shown in Figure 6A and for the Accumulation Season period in Figure 6B. The mean and median difference during the Melt Season period are nearly all less than 0.01 when averaged for all 30 images ( Table 2) for both MODIS and VIIRS f SCA . Figures 6A,B show the mean, whereas Figures 6C,D show the median over time. Median errors are smaller than mean errors, indicating a few larger errors as opposed to a more normally distributed set of errors. This also indicates that the accumulation season has more large errors than the melt season, which is also seen in the RMSE errors plotted in Figure 6F and compared to Figure 6E. The RMSE, like all the other statistics, was lower during the melt season and higher during the minimum and the start of accumulation.

DISCUSSION
We compare the distribution of fractional snow cover values from MODIS and VIIRS to each other and examine differences across elevation, slope, and aspect gradients. We present a comparison of two common metrics-total snow-covered area and regional snowline elevations-directly between MODIS and VIIRS to demonstrate continuity between the instruments.

Distribution of Fractional Snow Cover
A direct comparison of snow-covered area from the MODIS and VIIRS f SCA products (Figure 4, columns 4 and 5, respectively) on the evaluation dates (Supplementary Table S1) returns a MAE (Eq. 7) of 10.5% with an mean difference of +3.2%. When we expand the evaluation to the full set of coincident scenes for the period 2012-2015 for the two tiles shown in Figure 2 and limit to clear sky dates (<2% cloud), we find the MAE in snow cover area remains relatively consistent at 10.8%, and VIIRS tends to underestimate snow cover area by an average of −4.3%, indicating that VIIRS generally reports slightly less snow cover area than MODIS. The switch from positive VIIRS bias to negative VIIRS bias when the evaluation is expanded from 30 dates to three complete years is a reflection of the change in distribution of the months and seasons represented in the two samples. We consider the underestimation from VIIRS relative to MODIS of −4.3% to be more representative of overall performance. Note that these values are a measure of agreement between MODIS and VIIRS snow cover products rather than an evaluation with Landsat 8 OLI maps (Table 1; Figures 5,6), hence these values differ from those presented earlier.
A comparison of the frequency distribution of f SCA values between the two sensors in Figure 7 shows distributions and the overall median f SCA values to be very similar, 70% for MODSCAG and 71% for VIIRSCAG. These statistics and frequency distributions are consistent when evaluated for each tile individually, with median f SCA values of 67% and 70% for MODIS and 68% and 69% for VIIRS for the h23v05 and h24v05 tiles, respectively. Further analysis of the frequency distribution of f SCA values reveals that, overall, MODIS products find ∼7% purer snow-covered pixels than VIIRS with MODSCAG having more pixels covered (88%-100%), visible as taller red bars in Figure 7. The VIIRS products return more snow covered pixels (50%-88%), visible as taller blue bars in Figure 7.
The expected frequency distribution shift is likely a reflection of the differences in spatial resolution of the gridded reflectance data used as inputs to the snow retrieval algorithm; there are likely more mixed pixels at 1 km VIIRS resolution than at the 500 m MODIS resolution. The two MODIS tiles ( Figure 2) are exceptionally clear on September 18, 2013 (<0.8% clouds). Figures 8A,B show the MODIS and VIIRS f SCA, respectively, along with spatial f SCA differences between the two snow cover products in Figure 8C. The comparison between Figures 8A,B shows good agreement, with a higher abundance of high f SCA values in the MODIS data leading to negative differences in Figure 8C, which are spatially correlated with the 88%-100% f SCA pixels in the MODIS data (see Figure 7). Figure 8D shows the vegetation type (Arino et al., 2010) while Figures 8E,F show the elevation and slope.
These data are used for the topographic analysis shown in Figures 9, 10.

Performance Similarity Across Topography
The f SCA retrievals for both sensors show similar performance across elevation, slope, and aspect gradients in tile h23v05 ( Figures 9A-C), where we also compared Landsat 8 OLI in Results. The performance between the sensors is nearly the same for tile h24v05 (not shown), although the distribution of snow cover fraction with elevation is modulated by a more variable  (Figure 2), difference in f SCA (C). For the analysis shown in Figures 9 and 10: Globcover land cover class (D), elevation (E), and terrain slope (F). The inset shows a zoom region over a 111 × 78 km region as marked by the red viewport in panel (A). The satellite view zenith for MODIS was twice that of VIIRS in these data. hypsometric curve along with monsoonal climates that occur in the eastern portion of the Himalaya. Figure 10 shows f SCA compared across two MODIS tiles for September 18, 2013 segmented by vegetation types ( Figure 8D). We see larger f SCA values (up to 23%) retrieved from VIIRS in vegetated areas. However, in this region, there are very few pixels classified as "open forest" or "closed forest" relative to other land surfaces resulting in the relatively similar comparison shown in Figure 9. Therefore, the impact on total snow cover is relatively small. An analysis of the ancillary data shows that over the vegetated regions (open forest, closed forest, and shrub), sensor view zenith for MODIS is 44.6°and 45.9°, whereas for VIIRS they are 23.1°and 25.5°, respectively, for the two tiles. The larger view zenith for MODIS relative to VIIRS explains the higher estimate of snow cover from VIIRS. The analysis of differences related to view zenith is outside the scope of this paper and of interest for further research but is described in Evaluation Metrics of Dozier et al. (2008) and also in Rittger et al. (2020).
With respect to consistency between the MODIS and VIIRS fractional snow cover products, we compare snow cover area and regional snowline elevation across tile h23v05 at the common resolution of 1 km for 2013 to 2015. We observe a strong correlation between the snow cover area retrievals (R 2 0.98) and regional snowline elevation estimation (R 2 0.88) as shown in Figure 11. The strong agreement between MODIS and VIIRS  Figure 9B inset (VIIRS fSCA) shows more fractional snow cover than Figure 9A inset (MODIS fSCA) in open forest visible in Figure 9D supporting this conclusion.
FIGURE 11 | Comparison of MODIS and VIIRS-derived total snow cover area (A) and regional snowline elevation (B) for tile h23v05. The panels show data for 2012-2015 period on clear-sky days only i.e. < 2% cloud cover. All clear-sky days are plotted as black squares with the red crosses plotted over the observations from April-September.

CONCLUSION
An evaluation of new moderate resolution fractional snow cover maps from VIIRS retrievals (VIIRSCAG) in the Himalaya show performance statistics comparable to the more widely evaluated and available MODSCAG products. MODSCAG fractional snow cover maps re-evaluated using snow maps from Landsat 8 OLI do not show the negative bias found in a previous study. The reason for the improvement in MODSCAG performance revealed in this analysis is due to the improved Landsat 8 OLI validation data, more specifically improved sensor technology that avoids saturation in the visible bands and systematic overestimation bias in the validation maps from the assumption of 100% snow cover when bands 1, 2, and 3 were saturated in Landsat TM and ETM+. MODSCAG snow retrievals were expected to outperform those from VIIRSCAG in the steep terrain of the Himalaya, largely due to the finer spatial resolution of MODIS data. However, despite the coarser native resolution of the gridded products (1 km vs 500 m), VIIRSCAG evaluations produced slightly lower overall RMSE than for MODSCAG, perhaps because of improved representation of pixels toward the edge of the scan. The correlation and similarity in snow elevation and snow-covered area between MODIS and VIIRS observations over a three-year period, along with similar performance when compared to Landsat 8 OLI, show the ability of spectral mixture analysis to perform consistently over a range of spatial resolutions, making it a good candidate for multisatellite fusion efforts. In addition, the VIIRSCAG snow products appear to be suitable for extending the long-term SCAG timeseries (now at 20 years) well beyond the life of MODIS. Extending this record is important to performing snow cover change evaluation to monitor regional climate change. Future MODSCAG processing will benefit from MODIS Collection 6 surface reflectance products, which may allow processing of Aqua data in addition to Terra.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: (Rittger, 2021