Canopy Adjustment and Improved Cloud Detection for Remotely Sensed Snow Cover Mapping

Maps of snow cover serve as early indicators for hydrologic forecasts and as inputs to hydrologic models that inform water management strategies. Advances in snow cover mapping have led to increasing accuracy, but unsatisfactory treatment of vegetation's interference when mapping snow has led to maps that have limited utility for water forecasting. Vegetation affects snow mapping because ground surfaces not visible to the satellite produce uncertainty as to whether the ground is snow covered. At nadir, the forest canopy obscures the satellite view below the canopy. At oblique viewing angles, the forest floor is obscured by both the canopy and the projection of tree profiles onto the forest floor. We present a canopy correction method based on Moderate Resolution Imaging Spectroradiometer satellite imagery validated with field observations that mitigates geometric and georegistration issues associated with changing satellite acquisition angles in forested areas. The largest effect from a variable viewing zenith angle on the viewable gap fraction in forested areas occurs in moderately forested areas with 30–40% tree canopy coverage. Cloud cover frequently causes errors in snow identification, with some clouds identified as snow and some snow identified as cloud. A snow‐cloud identification method utilizes a time series of fractional vegetation and rock land‐surface data to flag snow‐cloud identification errors and improve snow‐map accuracy reducing bias by 20% over previous methods. Together, these contributions to snow‐mapping techniques could advance hydrologic forecasting in forested, snow‐dominated basins that comprise an estimated one fifth of Northern Hemisphere snow‐covered areas.


Introduction
Water stored in seasonal snow serves as an important water resource to downstream areas across much of the Northern Hemisphere during seasons when rain-fed supply does not meet the demand (Mankin et al., 2015). Tracking changes to snow-cover patterns is ever more critical for water-resource management in the current era of hydrologic nonstationarity amidst a warming world (Bormann et al., 2018;Milly et al., ©2019. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2019WR024914
Special Section: Advances in remote sensing, measurement, and simulation of seasonal snow Key Points: • The canopy correction methodology improves remotely sensed snow mapping accuracy in forested environments, reducing bias by 20% • A novel cloud detection algorithm that uses fractional vegetation and rock coverage can reduce snow-cloud identification errors • Improved performance over previous snow mapping approaches enable melt out date identification within days of actual snow disappearance Supporting Information: • Supporting Information S1 • Data Set S1 • Figure S1 • Table S1 • Table S2 • Table S3 • Table S4 • Table S5 • Table S6 • Table S7 • Table S8 • Table S9 • Table S10 • Table S11 • Table S12 • Table S13 • Table S14 • Table S15 • Table S16 2008). Historical norms of snow-cover extent (Mudryk et al., 2015), accumulation (Kunkel et al., 2016), and melt patterns (Musselman et al., 2015;Skiles et al., 2015) no longer reliably indicate snowmelt magnitude or timing. In recent decades, snow monitoring has been revolutionized by new satellite capabilities, especially imagery from the National Aeronautics and Space Administration's Moderate Resolution Imaging Spectroradiometer (MODIS) enabling the MOD10A1 snow products at daily, global coverage (Hall et al., 2002) and the MODIS Snow Covered-Area and Grain Size (MODSCAG) products that provide fractional cover and snow grain size over selected areas (Painter et al., 2009). Snow-cover maps serve as both early indicators for hydrologic forecasts and as inputs to hydrologic models that inform efficient water management strategies. Advances in snow-cover mapping have led to increasing accuracy of forecasts, especially during hydrologically important accumulation and melt periods .
Even with improved snow-mapping techniques, clouds and forests present persistent challenges to snow mapping. Cloud cover can create errors in spectral snow-detection algorithms because of the similar reflectance of cloud and snow at wavelengths utilized for snow identification. Clouds are often misclassified as snow leading to an overestimation of snow-covered area on cloudy days, and snow is often misclassified as cloud leading to an underestimation. Discriminating snow from clouds and developing cloud-free products has seen much attention in the literature (Allen et al., 1990;Dozier, 1989;Gafurov & Bárdossy, 2009;Gao et al., 2010;Hall et al., 2010;Parajka & Blöschl, 2008). Forested areas challenge remotely sensed snow-mapping techniques because snow maps only represent the viewable snow-covered area fraction (f SCA ) of a satellite pixel. Optical sensors only view snow that is visible beneath leafless deciduous stands (such as aspen), in clearings and forest gaps between coniferous trees (i.e., pine, spruce, and fir) and through thin foliage or needles, or snow that has been intercepted by the forest canopy. Often, viewable f SCA maps do not represent the actual snow cover on the land surface in forested regions for two reasons: (1) The forest floor is obscured directly by the canopy or tall understory shrubs, thus shielding snow on the ground from nadir satellite view, and (2) for whiskbroom sensors like MODIS with wide swath widths, the oblique perspective introduced by off-nadir satellite angles decreases the viewable gap fraction (VGF) and stretches the pixel as the viewing zenith angle (VZA) increases (Liu et al., 2004(Liu et al., , 2008Xin et al., 2012) because vertical features like trees become projected, obscuring the sensor view at oblique angles ( Figure 1).
Binary snow mapping (Hall et al., 1995) assessments indicate 10% of snow cannot be detected on a single date in the Sierra Nevada under forest with 98% accuracy, but only for pixels covered with 60% snow or more. Treatment of vegetation has been incorporated into this snow-mapping algorithm based on findings in Klein et al. (1998), who used the normalized difference vegetation index to determine where trees were impacting retrievals. Other fractional snow estimates for MODIS that account for forest exist (Metsamaki et al., 2012) but do not explicitly deal with off-nadir viewing. Past studies have attempted vegetation corrections utilizing multiple satellite products with mixed results (Dressler et al., 2006;Durand et al., 2008;Raleigh et al., 2013) and interference from vegetation degrades hydrologic prediction and simulated streamflow in forested, snowy regions (Yatheendradas et al., 2012). With nearly one fifth of the Northern Hemisphere seasonally snow-covered area estimated to intersect with boreal forests (Rutter et al., 2009) and with high proportions of catchments serving anthropogenic use being snow covered (Klein et al., 1998), poor snow maps over forested areas could lead to false interpretation of the extent of seasonal snow resources. More reliable water forecasts in snow-fed, forested basins require improved remotely sensed snow maps that better account for clouds and vegetation.
To address snow-mapping inaccuracies in forested regions caused by reductions in viewable gaps imposed by off-nadir VZA (Figure 1), we present a method to adjust viewable f SCA from MODIS to account for vegetation's effect on viewing geometry. The VZA is the view angle to the sensor from the ground, measured from zenith. Because of Earth curvature, the maximum VZA for MODIS, 65°, is larger than the maximum nadir view angle, 55°, from the satellite. The VGF-VZA relationship necessitates including the area within the stretched pixels that changes daily as VZA and viewing azimuth angle vary. Our method improves correction for vegetation because (1) it does not require explicit, higher-resolution data about the canopy itself and (2) it utilizes products derived from a single satellite (in this case MODIS), making it wholly self-consistent across data sets. The fractional areas of vegetation, soil, and snow calculated via the canopy-correction process presented here can also be utilized to detect misidentification of clouds as snow. This new technique for improved cloud detection utilizes a MODIS time series that flags abrupt changes to surface reflectance that indicate the presence of clouds.

Evolution of Satellite-Derived Snow-Cover Products and Their Use in Hydrologic Models
Snow-cover mapping techniques have evolved to improve temporal and spatial resolutions using a variety of sensors. Initial binary snow-mapping techniques classifying a pixel as "snow" or "no snow" (Dozier, 1989;Hall et al., 1995) gave way to more sophisticated techniques allowing the determination of pixelscale fractional snow-covered area (f SCA ; Salomonson & Appel, 2004 based on the Normalized Difference Snow Index (NDSI), which leverages the difference of snow's reflectance in the visible and shortwave-infrared (SWIR) wavelengths. Snow mapping further improved with the development of physically based spectral mixture analysis algorithms that utilize seven bands of data (Painter et al., 2009;Sirguey et al., 2009). These algorithms provide more accurate estimates of viewable f SCA than those based on NDSI, especially during hydrologically important periods  of most interest to water managers.
The use of satellite snow maps to aid hydrologic prediction of snowmelt runoff has been explored since the 1970s (Rango et al., 1977;Rango & Martinec, 1979). Applications of snow-cover maps include their use for validating modeled snow cover and snow-depletion curves (Shamir & Georgakakos, 2006;Wrzesien et al., 2015) and reconstructing peak snow water equivalent (SWE) distribution when combined with spatially distributed snow-melt models Bair et al., 2018;Durand et al., 2008;Rittger et al., 2016). Using snow maps as inputs to hydrologic models has yielded mixed results. Direct insertion approaches have shown improvements in estimating SWE (Rodell & Houser, 2004) and have reduced error and improved performance for streamflow forecasting (McGuire et al., 2006). More sophisticated assimilation approaches (e.g., ensemble Kalman filter, particle filter) have shown only slight improvements in accuracy for estimating SWE during melt (Andreadis & Lettenmaier, 2006) or modeled streamflow (Clark et al., 2006;Liu et al., 2013). Assimilation of MODIS f SCA into the U.S. National Weather Service SNOW-17 model resulted in a degradation in simulated streamflow due to the lack of MODIS f SCA canopy adjustment (Yatheendradas et al., 2012). The diversity of results achieved by integrating MODIS snow maps into existing models, in part, highlights the inadequate treatment of MODIS f SCA products in river basins with varied vegetation patterns and densities.

Previous
Approaches to Address Outstanding Challenges 2.2.1. Viewable and Forest Floor f SCA Assumptions are often made that the snow cover under the canopy is the same as in the gaps, yielding the canopy adjusted snow cover to be However, the assumption that viewable snow cover mimics hidden snow is problematic because vegetation affects snow processes in myriad ways. Vegetation intercepts falling snow that reduces accumulation and increases sublimation (Essery et al., 2003;Pomeroy et al., 1993;Winkler et al., 2005), changes the energy balance in relation to radiation sources for melt (Niu & Yang, 2004), and modifies turbulent fluxes at the snow surface that affect sublimation (Varhola et al., 2010). These dynamics vary at small scales over time and space (Clark et al., 2011) and are often climate specific . Furthermore, the effect of tree cover on snow inherently depends on the nature of the canopy profile-in forests with frequent fire regimes (i.e., ponderosa pine or mixed-conifer), foliage may be arranged well above the snow surface, and in forests with infrequent fire regimes (i.e., spruce/fir forests), foliage may extend down to or even below the level of settled snow (Fites-Kaufman et al., 2006). All of these factors make snow cover obscured by a forest canopy difficult to predict using rule-based algorithms or physically based models (Clark et al., 2011). Because snow accumulation and ablation patterns in forested areas often differ from open regions, the snow cover from openings seen by optical remote sensors may not accurately represent the surfaces obscured by vegetation. In this paper we assume equation (1) holds, and we adjust undercanopy snow surfaces from a temporally and spatially variable VGF, but we acknowledge this is an important area for further research.

VZA's Effects on VGF
The relationship between VZA and VGF introduces georegistration inaccuracies, daily changes to the area within stretched pixels acquired off-nadir, and diminishing VGF with increasing VZA. The first issue is that georegistration error of nadir MODIS pixels is on the order of 50 to 150 m and can increase nonlinearly with increasing VZA (Wolfe et al., 2002). This leads to a degradation of the georegistration in conversion to the sinusoidal projection used for data distribution and further degradation by transformation to the user's selected projection. Without correcting for the georegistration error, modeled pixels may not match the actual location of the MODIS pixels. We largely side-step georegistration errors in the vegetation correction because, in estimating VGF, we utilize a vegetation fraction calculated from the same algorithm, as opposed to integrating a separate satellite data set. Second, the instantaneous field of view of MODIS is about 500 m at nadir, but as the MODIS VZA increases to 65°at the edge of the swath, pixels grow by nearly a factor of 10, 2 times in the along-track direction and 4.8 times in the cross-track . Consequently, our estimates of canopy properties need to include the area within the stretched pixels that change daily.
The final and more complex obstacle relates to accommodating the decreasing VGF with increasing VZA (Liu et al., 2004; Figure 1). Previous geometric approaches to address VGF-VZA relationships include geometric optical (GO) models that parameterize tree crowns as cones (Li & Strahler, 1985, 1986 or ellipses (Li & Strahler, 1992). GO models capture the basic shape of the relationship between VGF and VZA, but hemispherical photos show they underestimate the VGF because the models do not account for "withincanopy gaps," solar radiation that reaches the ground surface despite being filtered and scattered through thin foliage (Liu et al., 2004). Efforts to capture the VGF by using a hybrid GO radiative transfer model (Li et al., 1995) that accounts for within-canopy gaps have shown improvements over the simpler GO model. However, in addition to readily available data like terrain slope and aspect, VZA, and satellite viewing azimuth angle, GO and GO radiative transfer models require detailed knowledge of the tree canopy (e.g., vertical and horizontal crown radius, tree density, and foliage area volume density). Ground-based data and lidar have been used to estimate relevant forest properties (Varhola & Coops, 2013); but lidar is not typically available for larger regions, remote locations, and the global scale that MODIS covers.
A more direct way of approximating the VGF used to adjust snow cover for modeling SWE combines the fine resolution (30 m) of the National Land Cover Dataset (NLCD; Homer et al., 2015) to estimate f VEG with the temporal variation of moderate-resolution imagery to estimate viewable f SCA . The assumption is that VGF = 1 − f VEG so that equation (1) becomes An alternate approach to account for pixel elongation and to reduce the effect of VGF, when interpolating daily fractional snow cover, is to assign less weight to observations at larger VZA and to rely more heavily on the observations from nadir Dozier & Frew, 2009). Still, the off-nadir views have some weight, so snow cover is likely to be underestimated when off-nadir views are cloud-free and nadir views are cloudy. For example, previous research using viewable MODSCAG f SCA data and a static VGF derived from NLCD showed improved accuracy over viewable f SCA but underestimated it compared to ground observations . We refer to this as the "static method" in following sections.

MODSCAG
The MODSCAG algorithm (Painter et al., 2009) uses linear spectral mixture analysis to determine subpixel coverage by snow, vegetation, and soil. MODSCAG models each pixel as the linear combination of three physical endmembers: (1) snow (2) vegetation, and (3) soil in two endmember mixtures to match apparent surface reflectance data from MOD09GA (Collection 5 in this paper). A shade fraction is calculated as the additive complement to 1 to the pixel's sum of physical endmember fractions to account for changes in atsurface irradiance expressed in the apparent surface reflectance of MOD09GA. Both snow and vegetation or soil endmembers are scaled by the shade fraction conserving the retrieved ratio such that the sum of physical endmember fractions equals 1.0. The output includes f VEG and f SOIL (fractional soil) maps in addition to those of f SCA and grain size. An example of this output is shown in Figures 2b-2d as compared to the MODIS false color image ( Figure 2a). We refer to "viewable MODSCAG f SCA " simply as "MODSCAG f SCA " hereafter.

VZA-Dependent Canopy Adjustment
We developed an observation-based approach that utilizes the relationship of MODSCAG f SCA and MODSCAG f VEG to estimate VGF. Using MODSCAG f VEG , VGF is then estimated as VGF = 1 − f VEG , and canopy-adjusted f SCA is calculated by equation (1).
This approach improves over previous uses of equation (1) because of the internal consistency we attain by deriving f VEG and viewable f SCA from a single sensor. By using simultaneous observations from MODSCAG to estimate f VEG and viewable f SCA , we also impose a physical constraint (total fraction of surfaces = 1), whereas this is not guaranteed if combining f VEG and viewable f SCA products derived across multiple sensors (e.g., adjusting MODIS f SCA with Landsat NLCD f VEG via the static method).
We explore the relationship of viewable MODSCAG f SCA and f VEG with VZA to test how VGF affects retrievals. To evaluate the relationship of f SCA with VZA, an independent estimate of VGF at nadir is needed. "Percent tree canopy" (PTC) from the NLCD (Homer et al., 2015) shows the best agreement to nadir f VEG from MODSCAG as compared to other available products and was selected as the independent estimate of VGF.
We utilize viewable MODSCAG f SCA maps for 2001-2012 for the Sierra Nevada together with VZA from the MODIS atmospherically corrected reflectance product, MOD09GA (Collection 5), to examine VZA's effect on f SCA . This range of years includes relatively wet years such as 2006 and 2011 and relatively dry years such as 2007 and 2012. The data are analyzed in increments of 10% tree canopy cover using NLCD's PTC data. We match each pixel's f SCA and f VEG with its VZA for the 12-year period and calculate the mean f SCA in VZA increments of 5°. However, f VEG from MODSCAG includes both tree cover and other shorter vegetation (e.g., low shrubs and herbaceous plants) that would not obscure the viewing angle. To distinguish visible shrubs and herbs from tall shrubs and tree cover, we use tree height data from the LANDFIRE 30-m data set (Rollins, 2009) to filter f VEG to MODIS pixels that have mean canopy heights greater than 2.5 m.

Complementary Cloud Detection
This analysis of MODSCAG creates a new procedure for detecting clouds utilizing f VEG and f SOIL that flags snow-cloud confusion errors. Discriminating clouds and snow has been a topic of research since the inception of optical remote sensing (see section 1). In a binary classification of pixels, clouds are brighter in the SWIR; so these wavelengths can be used to differentiate clouds from snow (Dozier, 1989). If viewable f SCA is above a fractional pixel threshold, typically 0.9, the retrieved grain size can also be used to distinguish snow from clouds . For pixels that have snow-cover fractions less than this threshold, cloud discrimination remains a consistent problem, especially for ice clouds (Allen et al., 1990).
Here we present an additional cloud-detection approach that can be used together with existing algorithms to improve snow-cloud discrimination. This new approach uses the daily sequence of f VEG and f SOIL to detect clouds that would otherwise remain undetected. To filter clouds that are not previously eliminated in forested regions (vegetation height > 2.5 m) through the SWIR or grain-size methods described above, we calculate the 70th percentile value of f VEG for each pixel over the period of interpolation and remove days with less than 5% of the 70th percentile f VEG value. These f VEG thresholds were determined based on analysis that showed they reject the most errors; refining these thresholds is an area for further research. Understanding a pixel's most likely vegetation fraction gives us insight into cloud detection described in section 4. After the cloud detection step, we apply the adjustment as in equation (2).

Canopy-Adjusted, STC MODSCAG f SCA
We then apply an interpolation with smoothing splines in the time dimension that weighs nadir observations more than off-nadir observations as in Dozier et al. (2008) to the canopy-adjusted f SCA . In the final step, we apply three-dimensional smoothing with a [3 3 3] box convolution kernel, thereby creating a canopy adjusted spatially and temporally complete (STC)-MODSCAG f SCA that we hereafter refer to as "STC-MODSCAG f SCA ." The spatial smoothing essentially enables fractional pixel values for the accumulation and melt season in forested areas. Nearby snow-covered alpine and meadow pixels included in the box convolution where f SOIL is retrieved instead of f VEG also contribute in this way.

Field Observations for Canopy-Adjusted Snow Map Validation
Grids of temperature measurements at 5-cm depth below the soil surface in four plots ( (Lutz et al., 2012), is analyzed here for the first time. All hourly ground temperature data and derived daily f SCA time series are included with this paper (see Supporting Information S1).
The temperature sensors detect snow, after the snowpack is established, by taking advantage of the insulating effect of the snow blanket from ambient temperatures. When snow is present, the near-surface ground temperature remains nearly constant, usually near 0°C, with little diurnal variability (Lundquist & Lott, 2008;Tyler et al., 2008). This approach has been successfully used in several studies for mapping snow duration and snow cover in forests and mountains (Dickerson-Lange et al., 2015;Ford et al., 2013;Schmid et al., 2012;Selkowitz et al., 2014).
Accuracy of snow-presence detection with the sensors is challenged in the early season when tracking of diurnal temperature variation can be confused by ephemeral snowfall events (i.e., episodic snow cover that persists for 1-3 days and then melts away). Prior studies have shown agreement to within ±1 day in both snow duration and the timing of snow disappearance derived from ground-temperature sensors versus time-lapse cameras when the pixel location of temperature sensors is known in the camera view. Likewise, f SCA depletion derived from temperature sensors has achieved high-temporal correlation (R 2 = 0.98) with f SCA derived from time-lapse imagery and agreed to within 4% of f SCA from 15-m maps from Advanced Spaceborne Thermal Emission and Reflection Radiometer . Point-scale comparisons between temperature sensors and automatic snow sensors are more difficult to interpret due to (1) differences in measurement support (i.e., 1.7-cm iButton vs. 100-to 300-cm sensing area for acoustic snow depth sensors or snow pillows) and (2) ambiguities in automated observations of snow disappearance (e. g., acoustic rangers may measure grass height rather than snow depth once snow cover becomes thin). However, temperature sensors installed near automatic snow measurements show reasonable correspondence in observing snow disappearance, typically within 2 or 3 days (Lundquist & Lott, 2008;Merio et al., 2018). When snow depth declines to less than~30 cm, the insulating effect of the snowpack diminishes (Taras et al., 2002), and damped diurnal variations emerge in the temperature of the snow-soil interface. Our methodology for identifying snow presence from ground temperature accounts for the reduced insulation in thinning snow by permitting diurnal variations of 1°C (daily range), as Raleigh et al. (2013) advise. We apply the same technique to the temperature data to derive daily snow presence at each temperature sensor and then estimate daily fractional snow cover based on the number of sensors reporting snow presence in a network.
3.6. Evaluation Metrics 3.6.1. Binary Metrics To quantitatively assess errors between observed f SCA (i.e., inferred from temperature sensor networks) and STC-MODSCAG f SCA , we use a set of binary metrics that rely on the common classification of four possible outcomes in identifying a MODIS pixel as containing snow or not: true positive (TP), true negative (TN), false positive (FP), and false negative (FN). These outcomes are used to calculate the binary statistics Precision, Recall, and F score in equations (3)-(5) typically used to asses snow cover accuracy (Masson et al., 2018;Painter et al., 2009 ;Raleigh et al., 2013 ;Rittger et al., 2013). These three statistics do not rely on TNs, when both STC-MODSCAG f SCA and ground-observed f SCA are zero to avoid inflating statistical performance.
Publicly available viewable MODSCAG f SCA products mask out values less than 0.15 to eliminate FPs; hence, previous MODSCAG-based work does not include these low f SCA values (Painter et al., 2009;Raleigh et al., 2013;Rittger et al., 2013). In a recent study, Masson et al. (2018) assessed products from another spectral mixture algorithm for snow cover, MODiMLAb, that retain values for f SCA less than 0.15 by using an NDSI threshold greater than −0.2 to mask out FPs instead of using the <0.15 f SCA threshold. In our analysis, we calculate the full range of values, 0 to 1, in a similar manner to Masson et al. (2018) shown in Table 1.
By including values of 0.01 to 0.14 in the analysis, we are essentially making it harder to score well in Precision because differentiating snow in that range from spectrally bright soils can be difficult resulting in more FPs. Simultaneously, Recall would be easier to score well because there should be fewer FN as the f SCA is not set to zero. These differences balance each other in F scores, but we considered evaluating the full range as more justified.

Fractional Metrics
To quantitatively assess the fractional errors between observed temperature sensor f SCA and STC-MODSCAG f SCA , we calculate mean and median difference and root mean squared error (RMSE) shown in equation (6) where N is the number of observations.
Like Rittger et al. (2013), we compare only observations where either MODIS f SCA or ground-observed f SCA are greater than zero (i.e., excluding TNs like with the binary analysis); but here we include values between 0.01 and 0.14 (Masson et al., 2018). Excluding these TNs results in fewer observations to compare in the fractional analysis than the binary analysis as indicated in lower fractional counts than binary counts in Table 1.

Relationship Between Viewable MODSCAG f SCA and VZA
The advantage of using the f VEG approach with MODSCAG instead of the static adjustment is that we can apply a per-pixel adjustment with a concurrent daily estimate of f VEG to account for VGF using equation (1). Figure 4 shows the viewable MODSCAG f SCA (also interpolated and smoothed), STC-MODSCAG f SCA , and f VEG for 1 March and 1 July 2012 in the Tuolumne and Merced River basins. On 1 March, nearly half the basin's snow cover needs to be adjusted. In contrast, on 1 July, after much of the snow has melted in the subalpine forests, the alpine snow cover does not require an adjustment. Also, visible on 1 July is an increase in f VEG relative to 1 March most likely from the spring greening of the understory in forested areas. An apparent increase in f VEG due to VZA is a less plausible explanation because VZA was 10°smaller for 1 July than 1 March.
While Figure 4 shows an example of our snow and vegetation maps, Figures 5a and 5b demonstrate the relationship of viewable f SCA and f VEG with VZA for the entire region covered in Figure 2. Land-surface changes (e.g., f SCA and f VEG ) due to PTC are shown for the spectrum of VZAs in Figure 5c. The analysis used data across the Sierra Nevada ( Figure 2) from 2001 to 2012. Changes to viewable MODSCAG f SCA and f VEG caused by VZA are smallest for dense forest canopies (≥75 PTC), with only an 8% change in land-cover fraction (f SCA or f VEG ) resulting from 65°VZA as compared to nadir (Figures 5a and 5b). This is expected in dense forest canopies where most of the forest floor is obscured regardless of the satellite's position overhead, so there is little possible change in land-cover fractions with increasing VZA (Figures 1 and 5c). The greatest effect of VZA on land-surface fraction occurs between 30% and 40% PTC (Figure 5c). At these moderate forest densities, there is a linear response to f SCA and f VEG from VZA (Figures 5a and 5b). In contrast, in sparse canopies, the relationship between land-cover fraction and VZA is nonlinear. For PTC of 5%, there is little change (<0.02) to f SCA and f VEG land-cover fraction at VZAs less than 40°but a 0.2 change for VZAs between 40°and 65°.
For f SCA (Figure 5a), the y intercepts are smaller than the VGF calculated with NLCD; and for f VEG (Figure 5b), the y intercepts are larger. The overestimate of f VEG may be the result of MODSCAG's three-endmember approach, where the first endmember is snow, the second is vegetation, and the third is shade, and a no-soil endmember is used if the residual from spectral mixing is small enough. Occasional snow intercepted by the tree canopy may also explain the f SCA underestimate.
In a parallel analysis for MOD10A1 f SCA , the effect of VZA was twice as large as for viewable MODSCAG f SCA ( Figure S1) possibly because the NDSI method is not as sensitive to detecting snow in forests because of the limited use of spectral information. This analysis sheds light on the snow-research community's current priorities (e.g., National Aeronautics and Space Administration's SnowEx) to identify limitations of satellite observations and algorithms in snowy vegetated areas.

MODSCAG Endmembers for Cloud Detection
The improved cloud detection method is demonstrated by an example showing a daily sequence of adjusted f SCA and viewable f SCA , f VEG , and f SOIL (Figure 6). In Figure 6a, the area shaded in gray identifies that the images from 13 through 18 May 2011 are correctly classified as cloudy by either quality flags in MOD09GA or grain size filtering. Figure 6b shows the sequence of images centered on the pixel in Figure 6a, with noticeable clouds on 19 and 23 May, dates that were not flagged as cloudy in MOD09GA or grain size filtering. Figure 6c shows the reflectance of the pixel for each day, with clearly elevated near-infrared (0.86 μm) and SWIR (≥1.24 μm) reflectance on 19 and 23 May. Those cloudy days, which were not correctly filtered by quality flags or grain size, instead produce incorrect fractional land-cover maps (Figure 6a).
While viewable f SCA appears reasonable on 19 and 23 May given viewable f SCA in the adjacent days (Figure 6a), it is unlikely the solution is achieved for the right reasons. Instead of previously retrieved f VEG , we see higher values for f SOIL on 19 and 23 May as compared to the correctly identified days. Given the higher SWIR reflectance from the presence of clouds, the MODSCAG algorithm is compelled to address a physically unrealistic subspace of the spectral vector space (for assumption of reflectance only from the land surface) and in turn mistakenly results in producing f SOIL . Previous algorithms  do not address these un-mixing errors. To catch these errors, an additional flagging technique is implemented that identifies anomalous f VEG data points (refer to section 3.3) over the period of interpolation.
On 19 and 23 May, adjusted f SCA using equation (1) turns out to be the same as viewable f SCA because f VEG is erroneously assigned to zero due to cloud cover. Figure 6 shows the VZA on the top axis. These acquisitions are closer to nadir than some of their neighbors and, if not identified with this new technique, are weighted heavily in interpolating across time, resulting in a false decrease in canopy-adjusted f SCA . The data need to be discarded because the spectral mixture analysis is again forced to explore a nonphysical vector subspace with respect to land surface reflectance. This new algorithm presented here flags days like 19 and 23 May and excludes them from the temporal interpolation.

Comparison of Adjusted f SCA With Ground Observations 4.3.1. Seasonal Comparison
We compare the observed f SCA from the temperature sensors in four sites to the STC-MODSCAG f SCA in water years 2010, 2011, and 2012, starting on 1 December until melt out (Figure 7). Given the geolocation accuracies described in section 2.2.2, we use the mean f SCA of all MODIS pixels that contain at least one temperature sensor as shown in Figure 3.
As described above, the accuracy of temperature sensor data requires a sustained snowpack. For 2010 and 2011, the snowpack was established by 1 December, whereas in 2012, the snowpack was not established until late January. Viewable MODSCAG f SCA retrievals have best accuracy above a threshold of 0.15 , and by the time f SCA drops to this threshold (gray band in Figures 7 and 8), the largest SWE volume has already melted. We show this threshold so that readers can see how well STC-MODSCAG f SCA performs in periods where inputs from MODSCAG f SCA are uncertain and previously unused.
In 2010, the observed f SCA is approximately 1.0 for the entire season until melt onset for all sites. STC-MODSCAG f SCA slightly underestimates snow cover for TUM and DAN, the least forested sites. Raleigh et al. (2013) describe rock outcrops at TUM that remain snow free during the season, which can explain this difference since the temperature sensors were placed in the soils of meadows and forests, not over the outcrops. Similar to TUM, DAN also has rock outcrops. At the OCR site, because of the relatively dense forest cover and absence of rock outcrops, this apparent underestimate does not occur. For OCR in 2010, STC-MODSCAG f SCA matches the observed f SCA until the onset of melt in May, when it overestimates f SCA until a late-season storm fully covers the pixels again prior to a second onset of melt. Melt out begins mid-May for OCR, early-May for TUM, and early-June for DAN and is discussed in more detail in section 4.3.2.
In 2011, we have a fourth study site, the Yosemite Dynamics Forest Plot (FDP) with the densest forest cover. For TUM, DAN, and OCR, observed f SCA is again approximately 1.0 for the entire season until melt with similar small underestimates from STC-MODSCAG f SCA relative to observed f SCA for TUM and DAN but not OCR. For FDP, observed f SCA shows three melt events in mid-winter during early December, mid-January, and early February not seen in the STC-MODSCAG f SCA . In this dense canopy, STC-MODSCAG f SCA underestimates in late March and shows complete melt out well before the observed snow cover disappears. STC-MODSCAG f SCA does capture the late season storm in May but not the magnitude of the storm indicating a possible need for reducing or removing smoothing (Bair et al., 2019). In 2011 for OCR, STC-MODSCAG f SCA has similar performance as in 2010, but the observed f SCA does not have an early melt signal like 2010, so the satellite observations better match the ground data. In 2012, the snowpack was not well established until late January. Variability between early December and late January is difficult to interpret because the errors could be either from the temperature sensor or the STC-MODSCAG f SCA . After the snowpack is well established, performance is similar to the two previous years. Notably, two melt-season storms are captured by both the temperature sensors and MODIS, which are explored more in the following section.

Melt Out
The most hydrologically important part of the season to water managers is when the snow melts and disappears. Many forecasts rely on snowmelt depletion curves that describe the relationship between snow-covered area and SWE. The subplots in Figure 8 show the mean f SCA values defined in the previous section during the period from full snow cover to no snow (melt out).
TUM has the least PTC coverage. In 2010, the mean of the STC-MODSCAG f SCA drops to zero on 7 June, matching the threshold (0.15) melt out date. Actual full melt out occurs approximately 10 days later. In 2011, MODSCAG f SCA performs better than in 2010, with the individual STC-MODSCAG f SCA pixel values enveloping (not shown) the true day of melt out with the mean corresponding well to the actual melt-out date. In 2012, the threshold melt-out date is the same, but snow cover again persists for a week with values below 0.15.
At DAN, the STC-MODSCAG f SCA from individual pixels envelopes the true observed f SCA melt-out date all 3 years. The STC-MODSCAG f SCA are slightly elevated just prior to melt most likely from the influence of shrubs and grass in forests on the f VEG signal boosting STC-MODSCAG f SCA but possibly also from data gaps that are interpolated.
At  DAN and TUM in relation to the 25 May 2012 storm; however, the surrounding observations (in time) of zero have an unfavorable effect on the interpolated f SCA . Given the short duration of these events, there is likely some noise in the ground data as well. In summary, accounting for the sensitivity limitation of 0.15, STC-MODSCAG f SCA can estimate the day of melt (±3 days) except for the most vegetated plot, FDP, where it melts out 2 weeks early.

Binary Statistics
In Table 1, we summarize the binary statistics for all dates as well as seasonally segmenting statistics for the winter (December, January, and February; DJF), spring (March, April, and May; MAM), summer (June, July, and August; JJA), and fall (September, October, and November; SON).
Over the available years and sites, Precision, Recall, and F score range from 0.87 to 0.95, 0.86 to 0.97, and 0.86 to 0.96, respectively. Precision is highest for DAN followed by TUM, OCR, and FDP. Recall is highest for OCR, closely followed by DAN and TUM, then FDP. F score, which balances Precision and Recall, is highest for DAN followed by TUM, OCR, and FDP. F score is highest for DAN followed by TUM, OCR, and FDP, indicating the two sites with the highest levels of PTC (OCR and FDP) were the hardest to map correctly.

Fractional Statistics
RMSE is lowest for the least forested site (TUM) and highest for the most forested site (FDP) and ranges from 0.17 to 0.36. Seasonally, RMSE is lowest in the spring followed by summer, winter, and fall for all sites except at FDP, which also has the highest errors in the fall, but in contrast to the other three sites has relatively high errors in the spring as well.
The mean and median differences are calculated as the STC-MODSCAG f SCA minus the ground-observed f SCA . Mean differences for all dates range from −0.05 to 0.05, and median differences range from −0.07 to 0. Seasonally, these two statistics show that STC-MODSCAG f SCA overestimates in the fall; while the winter, spring, and summer are mixed with more underestimates than overestimates. The early melt out for FDP seen in Figure 8 accounts for large negative mean and median differences of −0.34 and −0.36 in the spring.
Mean and median differences for some sites differ in sign, indicating occasional larger overestimates but regular small underestimates; for example, the small underestimates at TUM and DAN attributed to rock outcrops unsampled in the ground data. Large overestimates as seen in January 2012 (Figure 7) could be from still undetected clouds or, alternatively, errors in the early season ground data before an established snowpack.
At OCR, overestimates and underestimates balance errors resulting in very small mean and median differences except in the fall where the not-yet-established snowpack in 2012 drives up the differences. Similarly, TUM and DAN have small mean and median differences throughout the season except in the fall.

Improvements in Snow Mapping Accuracy
Previous validation efforts that use gap-filled time series from MODIS and evaluate the remote sensing observations under various canopy densities are scarce. The commonly used MOD10C1 product has been 10.1029/2019WR024914 gap filled and evaluated by assimilation into a land-surface model but not compared to ground observations of f SCA (Hall et al., 2010). Alternative gap-filling approaches also exist (Parajka et al., 2010). These studies did not adjust f VEG dynamically with VZA, whereas this study accounted for these effects (Figure 1). Previous studies also did not use a physically based snow mapping model like MODSCAG with the advantage of concurrent f VEG and f SOIL estimates that can be used for better cloud filtering. We evaluated the adjusted f SCA using the same metrics as previous studies but for the full range of f SCA values from 0 to 1. Based on our melt-out comparisons shown in section 4.3.2 and the statistical comparison, STC-MODSCAG f SCA is valid below 0.15.
A recent assessment of existing methodologies for snow cover fraction from Masson et al. (2018)  also under tougher-to-predict cloudy sky conditions. To best understand improvements gained by the methods presented in this paper, we can directly compare to previous work on adjusting MODSCAG snow cover for vegetation using NLCD . In Table 1, "All" is most comparable to tables 2 and 3 in Raleigh et al. (2013). Finally, previous comparisons of STC-MODSCAG f SCA to airborne-derived snow cover shows similar accuracy even in extreme drought years such as 2015 .
With regard to Precision, values reported in Table 1 are 0.08, 0.01, 0.12, and 0.09 lower than values reported by Raleigh et al. (2013). This can be explained by Raleigh et al. (2013) setting all MODIS values of 0.15 to zero, decreasing the number of FPs and increasing the Precision (equation (3)). In contrast, values of Recall in Table 1 were 0.11, 0.00, 0.13, and 0.24 larger than in Raleigh et al. (2013) for precisely the same reason: We include the values below 0.15 for STC-MODSCAG f SCA boosting the number of TPs while simultaneously decreasing FNs resulting in a higher Recall (equation (4)). The F score that balances Precision and Recall is equal to or greater (0.02, 0.00, 0.01, and 0.11) in this study than in Raleigh et al. (2013), showing the clear improvement even while considering the full range of possible snow fractions from 0 to 1.
RMSE in this study is 0.08, 0.03, 0.01, and 0.19 less than in the previous study at TUM, DAN, OCR, and FDP, respectively, showing the most significant improvements at the least and most forested sites. Most notably, the mean bias has been reduced to 0.02, 0.03, 0.05, and −0.05 from previous values of −0.22, −0.09, −0.09, and −0.37 at TUM, DAN, OCR, and FDP, respectively, for a 20% overall reduction in bias. In addition to the binary statistics, these fractional statistics also show the clear improvement using the new techniques presented in this paper.
Although not obvious from these statistics, there is improved ability to estimate f SCA at OCR when f SCA is less than the PTC of 66%. Scatterplots (Figure 9) presented over the same time frames as in Raleigh et al. (2013) show the method presented in this paper can provide more accurate results everywhere but show a considerable improvement at moderately forested sites like OCR. In contrast, the static method did not perform well when f SCA was less than the complement of the f VEG (i.e., to the left of the black lines in Figure 9). This issue is only visible in FDP with the largest PTC of our sites. At FDP, canopy cover is both taller (maximum 67 m) and more variable than at the other sites. Lower values of f SCA shown here are generally in later months and near melt out; the improved cloud detection appears to help considerably. For example, in Raleigh et al. (2013) in 2010 and 2011, OCR melted out 2 weeks too early; but our method captures the threshold melt out to within 3 days and actual melt out to within 5 days.
In summary, STC-MODSCAG performs favorably compared to the previous ground-based study with a higher F score and a lower bias and RMSE on average over most of the season, likely from the improved vegetation adjustment and more accurate melt-out dates due to improved cloud detection. The assessment provides errors segmented by season and forest density that can be used for data assimilation in future studies and applications utilizing the STC-MODSCAG product.

Limitations of the Approach and Future Improvements
Several outstanding areas of research not addressed in our approach may produce snow-mapping errors. As discussed in section 2.2.1, we assume that viewable f SCA in open areas can be applied to undercanopy forest floors. This increases snow-mapping error where and when the assumption is not true. Previous studies have shown that snow-cover differences between open areas and forests vary regionally . Developing a more informed approach to estimating f SCA on the forest floor remains an outstanding challenge that either requires additional remote-sensing observations, for example, lidar that can delineate the vertical profile of canopy cover (Kostadinov et al., 2019;Painter et al., 2016), or more reliance on an integrative framework that incorporates models with existing remotely sensed data to discriminate differences between open and subcanopy f SCA .
Snow intercepted by the forest canopy can also increase snow-map errors since viewable snow in the canopy is not differentiated from viewable snow on the ground. The issue is made more complicated by the differential time that snow remains in the canopy in warmer or cooler temperatures. We expect that the STC-MODSCAG f SCA will not be extremely sensitive to this issue. If there is snow in the canopy, there is likely to be snow on the ground, but we acknowledge viewable f SCA will be enhanced and f VEG will be reduced during periods with snow in the canopy. Discriminating snow in the canopy from snow on the ground is a possible avenue for future development and could be useful for evaluating model representations of snowforest interactions. Validation of our product is provided via four Sierra Nevada sites over 3 years, including only 1 year in the most heavily forested site. This limits our ability to extrapolate the accuracy of the product in heavily treed environments and over regions other than the Sierra Nevada that experience different patterns in snow processes and forest effects on snow cover. The temperature sensor network data we use provide ideal temporal information and sufficient spatial sampling for validation but lack the spatially intensive snow-cover mapping data afforded by other remote-sensing data sources such as lidar or very high-resolution optical products which are only available for cloud-free days. Possible future validation approaches include comparing the STC-MODSCAG f SCA to f SCA derived from lidar data sets in areas with varying forest cover or from time-lapse camera networks (e.g., SnowEx).
Finally, there is a sustained need for improved cloud masks. The approach presented here detects cloudy conditions not found in the original MOD09GA cloud mask, but geostationary or commercial satellites with sufficient spatial resolution and repeat cycles hold promise in helping to further develop and evaluate more advanced cloud masking algorithms.

Conclusions
We test an adjustment for MODIS snow-cover maps to account for VZA dependencies in forest cover and to better detect clouds. For STC-MODSCAG f SCA , estimates of snow-covered area are consistent with ground observations except for the densest canopy (i.e., forest cover exceeding 80%). The method of observing snow cover on the ground using temperature sensors is of great utility for evaluating optical satellite data and interpolation algorithms. However, more ground data are necessary to better understand the threshold at which optical data become unreliable for various tree species and densities. Rock outcrops should be quantitatively assessed during the nonsnow-covered period, for example, using high spatial-resolution imagery to improve the accuracy of the ground observations. Large-scale efforts such as SnowEx in 2017 that focused on identifying the forest density threshold for reliable snow remote sensing provide well-timed complimentary initiatives and field data upon which to further validate our canopy-adjustment approach. MODSCAG is a physically based model producing vegetation and soil endmembers that allows us to constrain cloud identification in ways not possible with conventional index-based algorithms like NDSI. Together, these contributions to snow-mapping techniques have the potential to advance natural resource management for quantifying the effects of forest disturbance on canopy structure due to insect-and pathogen-related infestation (Baker et al., 2017) and fire (Micheletty et al., 2014) as well as for hydrologic forecasting in the snowy forested basins that comprise one fifth of Northern Hemisphere snow-covered areas.

Data
Data for all study sites can be found in the supporting information as Excel (.xlsx) files. The files include temperature sensor data and derived f SCA as well as corresponding STC-MODSCAG f SCA . Viewable MODSCAG f SCA was obtained from JPL (https://snow.jpl.nasa.gov/portal/data/). STC-MODSCAG data sets used in this study can be found on the CWEST snow server at CU Boulder (ftp://snowserver.colorado.edu/pub/ fromRittger/201901_STC-MODSCAG_WRR_data). If the paper is accepted, these data will be submitted to an appropriate repository such as the NSIDC DAAC or a FAIR-aligned repository.