How do tradeoffs in satellite spatial and temporal resolution impact snow water equivalent reconstruction?
Given the tradeoffs between spatial and temporal resolution, questions about resolution optimality are fundamental to the study of global snow. Answers to these questions will inform future scientific priorities and mission specifications. Heterogeneity of mountain snowpacks drives a need for daily snow cover mapping at the slope scale (≤30 m) that is unmet for a variety of scientific users, ranging from hydrologists to the military to wildlife biologists. But finer spatial resolution usually requires coarser temporal or spectral resolution. Thus, no single sensor can meet all these needs. Recently, constellations of satellites and fusion techniques have made noteworthy progress. The efficacy of two such recent advances is examined: (1) a fused MODIS–Landsat product with daily 30 m spatial resolution and (2) a harmonized Landsat 8 and Sentinel 2A and B (HLS) product with 3–4 d temporal and 30 m spatial resolution. State-of-the-art spectral unmixing techniques are applied to surface reflectance products from 1 and 2 to create snow cover and albedo maps. Then an energy balance model was run to reconstruct snow water equivalent (SWE). For validation, lidar-based Airborne Snow Observatory SWE estimates were used. Results show that reconstructed SWE forced with 30 m resolution snow cover has lower bias, a measure of basin-wide accuracy, than the baseline case using MODIS (463 m cell size) but greater mean absolute error, a measure of per-pixel accuracy. However, the differences in errors may be within uncertainties from scaling artifacts, e.g., basin boundary delineation. Other explanations are (1) the importance of daily acquisitions and (2) the limitations of downscaled forcings for reconstruction. Conclusions are as follows: (1) spectrally unmixed snow cover and snow albedo from MODIS continue to provide accurate forcings for snow models and (2) finer spatial and temporal resolution through sensor design, fusion techniques, and satellite constellations are the future for Earth observations, but existing moderate-resolution sensors still offer value.
Mountain snowpacks are challenging for remote sensing because they change rapidly. Moderate-resolution sensors such as MODIS and VIIRS image Earth daily but at resolutions (463–750 m) that cannot resolve slope-scale features of interest to a variety of scientific users ranging from hydrologists (Blöschl, 1999) to the military (Vuyovich et al., 2018) to wildlife biologists (Conner et al., 2018). Finer-resolution multispectral sensors such as Landsat 8 and 9 provide spatial resolutions of 30 m but at 16 d revisits, during which time the snow cover can change considerably. Because of cloud cover, useable optical imagery with such infrequent revisits can be months apart. Recognizing that no single satellite/instrument can provide fine spatial and temporal resolution, constellations of satellites with coordinated overpass times have emerged. Two examples are the Sentinel 2 A and B and Landsat 8 and 9 pairs. For optical bands, the Sentinels image Earth every 5 d at 20 m, and Landsat 8 and 9 image Earth every 8 d at 30 m. The harmonized Landsat 8 and 9 and Sentinel 2 (HLS) product (Claverie et al., 2018) improves the average revisit time to 3–4 d at 30 m spatial resolution.
Effects of snow cover estimates at finer resolution have been examined in a few studies, showing a wide range of improvements in errors. In comparing snow cover depletion curves from Landsat Multispectral Scanner (MSS; 80 m pixels; 16 d repeat) and Advanced Very High Resolution Radiometer (AVHRR) imagery (1100 m pixels; daily repeat), Baumgartner et al. (1987) found that AVHRR tended to overestimate snow cover where it was patchier (lower elevations) and underestimate snow cover where it was more widespread (higher elevations) relative to MSS. They concluded that AVHRR imagery could be used to fill in temporal gaps in depletion curves generated from Landsat MSS. Luce et al. (1998) compared a spatially explicit snow water equivalent (SWE) model at 30 m with single and two point models for a small basin in Idaho. The 30 m model showed significantly lower errors than the single and two point models. Cline et al. (1998) examined the effect of upscaling the spatial resolution of a DEM and snow cover in an energy balance SWE model at a range of resolutions: 30, 90, 250, and 500 m. Positive biases in the coarser-resolution estimates arising solely from basin delineation artifacts were reported; thus the authors advise using vector basin outlines (as was done in Sect. 2.5). When these artifacts were corrected, the SWE volumes at 90 m were overestimates, while those at coarser resolutions were underestimates. Blöschl (1999) examined scaling issues in snow hydrology and showed that pixel sizes of a few meters are needed to accurately capture basin-scale SWE. Turpin et al. (2000) examined snow cover maps derived from AVHRR and Landsat TM (Thematic Mapper; 30 m resolution; 16 d repeat) and report discrepancies, also finding that AVHRR failed to resolve patchy snow compared to TM. Durand et al. (2008) were the first to create a fused MODIS and Landsat product. For the coarse-resolution product they used binary snow cover from MODIS (Hall et al., 2002). For the fine-resolution product they used Landsat 7 ETM+ surface reflectance in a spectral unmixing algorithm (Painter et al., 2003). The authors then used a linear program, constrained to match the ETM+ fractional snow-covered area (fsca) imagery while also matching the daily changes in fsca observed by MODIS. Applying their linear program to the upper Rio Grande, the authors found differences in fsca between the ETM+ fsca, the MODIS fsca, and the fused product. When run through a snow reconstruction model, these differences equated to a 51 % reduction in mean absolute error (MAE) and a 49 % reduction in bias for SWE using the fused snow cover versus the ETM+ snow cover. Using the same reconstruction model, Molotch and Margulis (2008) report a 23 % MAE in SWE using ETM+ snow cover versus a 50 % MAE in SWE using MODIS snow cover and 89 % MAE using AVHRR snow cover. Rittger et al. (2013) examined spectrally unmixed snow cover from ETM+ (similar to the approach in Painter et al., 2009) and several approaches for mapping snow cover from MODIS, including spectral mixture analysis (Painter et al., 2009). They found that ETM+ mapped consistently more patchy snow cover than the MODIS approaches, suggesting fewer false negatives and thus a greater recall statistic. Winstral et al. (2014) examined scale in a snow energy balance model at a range of spatial resolutions and found that 100 m spatial resolution is needed to accurately simulate snowmelt. Contrary to Cline et al. (1998), Schlögl et al. (2016) report that SWE increases with DEM resolution in two alpine basins. Similarly, Baba et al. (2019) used an energy balance model with a DEM at 8–1000 m and report good agreement with fine-resolution snow cover maps up to 250 m but a loss in agreement at coarser resolution likely due to excessive smoothing of topographic effects. Rittger et al. (2021) used a random forest to fuse spectrally unmixed snow cover from MODIS with Landsat 5 and 7 ETM+. The authors' comparisons show sharper snowlines (transition from no snow to fully snow covered) in the Landsat and fused imagery compared to MODIS, again indicating that Landsat may have greater recall than MODIS in this difficult to validate region. Bouamri et al. (2021) examined differences between snowmelt models with and without solar radiation represented in the Atlas Mountains of Morocco. Although the models with solar radiation better simulated the snow cover used for validation, aggregating the simulated snow cover from 100 to 500 m suppressed those improvements. In summary, many studies have compared coarse- and fine-resolution snow cover, but only three studies to our knowledge (Cline et al., 1998; Durand et al., 2008; Molotch and Margulis, 2008) have examined the impact of resolution on SWE reconstruction, all finding significant improvements from finer spatial resolution. Since those studies, considerable advances have been made in SWE reconstruction techniques (Bair et al., 2016; Rittger et al., 2016) as well as snow cover (Stillinger et al., 2023) and albedo mapping (Bair et al., 2019), hence the justification for revisiting the effects of spatial and temporal resolution.
Three daily snow cover estimates were used to force a SWE reconstruction model at two spatial resolutions: a baseline at 463 m using MODIS with the Snow Property Inversion from Remote Sensing (SPIReS; Sect. 2.1; Bair et al., 2021) and two 30 m estimates, one from the Harmonized Landsat Sentinel (HLS) surface reflectance product (also using SPIReS; Sect. 2.2) and the other from Snow Covered Area and Grain Size (SCAG)–Fusion (Sect. 2.3). The period covered is 1 January 2018 to 31 December 2020, limited by the intersection of the availability of the SCAG–Fusion and the HLS. The domain is the Tuolumne River basin above the Hetch Hetchy Reservoir in the Sierra Nevada, USA, because of the availability of Airborne Snow Observatory (Painter et al., 2016) estimates of SWE for validation. This approach rests on the hypothesis that (1) fsca and (2) snow albedo are the two most important variables in SWE reconstruction (see Sect. 2.4). The importance of fsca in SWE reconstructions can be traced to several studies (e.g., Durand et al., 2008; Molotch and Margulis, 2008). The importance of snow albedo in SWE reconstructions is shown in Bair et al. (2019).
The baseline case uses SPIReS to map snow cover from MODIS at 463 m daily resolution, although the effective pixel size can be up to 5× as large for off-nadir acquisitions (Wolfe et al., 1998; Dozier et al., 2008). The MOD09GA daily surface reflectances (Vermote and Wolfe, 2015) are unmixed into fsca and properties used to model albedo (grain size and dust concentration). These estimates are then run through a series of filters including persistence filters for clouds and time-based smoothing and interpolation.
One daily 30 m snow cover product used also comes from the SPIReS approach applied to the HLS. It is the first snow mapping application, to our knowledge, of HLS data. Thus, we describe the workflow in more detail than SPIReS–MODIS. The Tuolumne River basin above Hetch Hetchy Reservoir straddles two Sentinel tiles, so HLS version 1.4 multi-band hierarchical data format (HDF) files from both tiles were downloaded (https://hls.gsfc.nasa.gov/, last access: 23 November 2022). For calendar years 2018–2020, four combinations of products were downloaded: two tiles (11SKB,11SKC) and two products – S30 (Harmonized Sentinel-2 MSI) and L30 (Harmonized Landsat-8 OLI). We attempted to download the newer HLS version 2.0 from NASA Earthdata Search, but as of this writing, the S30 product for those tiles only extends back to 23 September 2020 because of daily limits on the number of Sentinel-2A and B scenes that can be downloaded by NASA from ESA for reprocessing. Seven bands covering visible through shortwave–infrared wavelengths were used from each sensor: 1–4, 8A, and 11–12 for S30 and 1–7 for L30. Mean local solar geometry was obtained from the accompanying header files. The multi-band images were stacked, mosaicked, and cropped to the basin to form a 311 (depending on the year) 4D data structure, with dimensions of rows, columns, bands, and time. Each Landsat has a 16 d revisit, thereby providing imagery at 8 d intervals for each tile, and each Sentinel has a 10 d revisit, providing imagery at 5 d intervals. Thus, combined revisits ranged from around 1 to 10 d with a mean of 3.5 d (Fig. 1). Theoretical revisit times estimated before Landsat 9 was launched (Li and Roy, 2017) for a three-satellite constellation at this basin's latitude shows a mean of 3.8 d, with a minimum of less than 1 d and a maximum of 7.0 d. There are only four revisit times greater than 7 d shown in Fig. 1; all other observations lie within the theoretical revisit times.
Red–green–blue band imagery for each day was examined visually. Days with clouds or incomplete spatial coverage over the watershed (many images have large areas with no data) were discarded. After filtering, 156 or about half of the days were kept. The minimum, median, and maximum time spacings between acquisitions after filtering were 1, 5, and 40 d. The SPIReS spectral unmixing approach was then applied to these filtered surface reflectances as described for Landsat 8 OLI in Bair et al. (2021), yielding the variables fsca, snow grain size, and dust concentration. A per-pixel spline interpolation was applied to each of the variables in the time dimension to make them continuous, covering all days from 2018–2020.
A second daily 30 m snow cover product used was a MODIS–Landsat fusion, created using two random forests for classification and regression based on previous work (Rittger et al., 2021) but retrained using Landsat 8 OLI data. Standard cloud masks (Foga et al., 2017) were used to select the 100 most cloud-free level 2 surface reflectance images (USGS, 2021) for dates spanning March 2013 to March 2021 (Fig. 2). Winter months had fewer training data than summer months. Six scenes were manually removed after visual inspection of red–green–blue imagery. Because the initial filtering did not remove all clouds, a second cloud filtering step with superpixels and Gabor filtering was used (Stillinger, 2019). This second step removed all the clouds but removed some snow cover as well. These filtered Landsat 8 surface reflectances along with MODIS MOD09GA surface reflectances were unmixed into fsca and snow surface properties that affect albedo (Painter et al., 2009, 2012). This SCAG spectral unmixing differs from the SPIReS approach; it finds the best fit from an endmember library for the snow-free parts of the pixel, whereas SPIReS uses an empirical snow-free endmember. There are other differences in the treatment of light-absorbing impurities, filtering, and time–space smoothing (Rittger et al., 2020). For more details and a recent comparison between SPIReS, SCAG, and other snow mapping algorithms see Stillinger et al. (2023). Estimates of fsca and snow surface properties that affect albedo, i.e., grain size and visible albedo degradation, were used as training data. Physiographic variables, including solar illumination and land classification were used as predictors. The two-step approach consists of an initial model that classifies pixels into three cases: (1) 0 %, (2) 100 %, or (3) 1 %–99 % fsca. For case 3, a second regression random forest was used to estimate fsca on the 1 %–99 % interval. This two-step classification–regression approach was found to be less biased at predicting 100 % snow-covered pixels than using a single-step random forest predicting 0 %–100 % fsca.
2.4 Parallel energy balance
At an hourly time step, the Parallel Energy Balance model ParBal (Bair et al., 2016, 2018; Rittger et al., 2016) downscales state and flux variables solving for the surface snow energy balance. The computed melt is multiplied by the fsca and summed backward from end-of-melt to peak SWE for each pixel to estimate SWE on the ground throughout the melt season. Evaluations of ParBal (Bair et al., 2016, 2018) forced with snow cover from the MODSCAG approach (Painter et al., 2009; Rittger et al., 2020) show a mean absolute error (MAE) of 22 %–26 % using SWE from Airborne Snow Observatory (ASO) for validation. There are two significant changes to ParBal here. (1) U and V wind component forcings use the hourly MERRA-2 (GMAO, 2015) data instead of N/GLDAS (Rodell et al., 2004; Xia et al., 2012). Using U and V components with global forcings allows for terrain-based wind downscaling using curvature and slope (Liston et al., 2007), whereas GLDAS only provides wind speed. (2) A new estimate of SWE on the ground, called hybrid SWE, leverages GLDAS SWE (the GLDAS NOAH 3 h 0.25∘ v2.1 model was used) and captures the accumulation phase. Previously, ParBal estimates were limited to the ablation phase only. The concept is to identify GLDAS pixels with similar snow cover duration as the fine-scale fsca pixels, find the peak SWE day from those GLDAS pixels, and then scale the GLDAS estimates by the ParBal SWE estimate on that peak day. This process is repeated for every fine-scale pixel. The GLDAS SWE, SWEGLDAS, is extracted for the domain, in this case a bounding box covering the Tuolumne River basin above Hetch Hetchy Reservoir. Pixels with the same snow cover duration are identified by the logical vector t as
where the asterisk denotes the selected pixels and fsca is from the fine-scaled product, i.e., Sect. 2.1–2.2. The Δt1 and Δt2 indicate different time periods (days). Because there can be multiple pixels with matching snow cover duration, the daily mean of is taken. The maximum value of that daily mean and its index imax are computed. A scaling coefficient c is calculated as
where SWEParBal,imax is the value of the reconstructed SWE from ParBal (Eq. 2 in Bair et al., 2016) at the time indexed by imax and the overbar denotes an average. The following case can arise:
For example, this case can occur when ParBal models all the mass loss via sublimation. For this case, SWEGLDAS is used when fsca >0. Otherwise, the hybrid SWE prior to the peak is set at day i as
This scaling can cause unrealistic daily increases and decreases in SWE; thus a smoothing spline is applied. This hybrid SWE has yet to be evaluated throughout the accumulation season, but comparisons with the reconstructed SWE during the ASO acquisitions show negligible differences, indicating at least the imax estimate is occurring roughly at the right time of year since all of the ASO flights examined here took place during the ablation season (with the exception of 13 April 2020; Sect. 3). Figure 3 shows this hybrid GLDAS and reconstructed SWE for an example pixel in water year (WY) 2019. ParBal was run with each of the snow cover forcings, holding all other inputs constant.
2.5 Airborne Snow Observatory
ASO 50 m SWE estimates for the Tuolumne River basin above the Hetch Hetchy Reservoir, which is the most sampled basin by ASO, were downloaded from the National Snow and Ice Data Center for 2018 and 2019 and from ASO Inc. for 2020 (Table 1). The number of acquisitions per year ranged from two (2018) to four (2019) with a total of nine. The accuracy of ASO measurements at the basin scale cannot be estimated directly from data, since there is no better method for validation, but since 2021, ASO has provided basin-wide uncertainty estimates in their reports available on their website (https://www.airbornesnowobservatories.com, last access: 23 November 2022), mostly based on uncertainty in modeled density, with a small uncertainty in depth. The reported mean basin-wide uncertainty in SWE for ASO flights for the entire Tuolumne River basin for 2021 and 2022 is ±4 %, so we assume similar errors in 2018–2020 and use that uncertainty estimate.
The ASO images were resampled from a cell size of 50 to 2000 m (4 times the MODIS resolution) and 120 m (4 times the Landsat resolution), using a mean-preserving technique with a weighted resampling covering the image (mapresize; MathWorks, 2022). The ASO images were kept in their native UTM 11N projection. The upscaled cell sizes account for geolocational and sensor-to-sensor uncertainty of 1–2 pixels for MODIS and Landsat–Sentinel-2 (Tan et al., 2006; Storey et al., 2016). The ASO dates in Table 1 were extracted for each of the three SWE reconstructions. Then, the matched baseline SPIReS–MODIS images were upscaled from 463 to 2000 m and reprojected from a sinusoidal projection to UTM 11N. The matched MODIS–Landsat fusion and SPIReS–HLS images were upscaled from 30 to 120 m but kept in their native UTM 11N projection. Vectors of the Tuolumne River basin above Hetch Hetchy Reservoir were obtained from ASO Inc. These vectors were then converted into coarse-resolution masks of the basin. Waterbodies and other areas with no data either in the ASO images or in the SWE reconstructions were removed from the masks. Areas outside of the masks were then set to null values. These common masks were applied to the upscaled SWE reconstruction and ASO images such that areas outside the masks were set to null values. The upscaled ASO images were compared with the upscaled SWE reconstruction images. The following error statistics were computed for a given date: bias as a measure of basin-wide error, relative bias normalized by ASO mean SWE, mean absolute error (MAE) as an unweighted measure of per-pixel error, and relative MAE normalized by ASO mean SWE:
where N the total number of pixels and j is an individual pixel. Mean values of the four error statistics were also averaged by year. MAE is used instead of root mean squared error because it evenly weights errors, which is preferred when comparing modeled values (Willmott and Matsuura, 2005), i.e., ParBal to ASO, neither of which directly measure SWE.
2.7 Snow albedo errors
Errors in snow albedo directly impact the accuracy of reconstructed SWE (Bair et al., 2019). However, for the spectral unmixing approaches used here, the albedo errors are low, evaluated using terrain-corrected measurements from Mammoth Mountain (e.g., Bair et al., 2022), only 23 km from Mount Lyell, the highest point in the Tuolumne River basin. For example, from water years 2017–2019, the root mean squared error (RMSE) for MODIS–SPIReS, calculated using the best value for a 3×3 neighborhood around the validation site, is 2.3 % with no bias (Table 2). These albedo errors are similar to the accuracy of the hemispherical directional reflectance factor (HDRF) surface reflectance products, evaluated over dark targets (Vermote et al., 2016; Bair et al., 2022). These improvements in remotely sensed snow albedo over previous assessments, showing RMSE values of 4.6 % to 4.8 % with 0.7 %–1.3 % bias for MODIS (Bair et al., 2019, 2021), come from improved cloud snow discrimination filters and adjustments to thresholds such as the minimum grain size for dirty snow (Sect. III-J of Bair et al., 2021).
We are not dismissing errors in albedo, as these remotely sensed snow albedo errors can lead to 5 %–11 % MAE for reconstructed SWE (Bair et al., 2019), but without independent measurements of spatially distributed albedo, we lack validation data for further error evaluation.
Basin-wide mean SWE errors are shown in Fig. 5a–c and in Table 3. The lower mean MAE, by 10 %–14 % (bold in Table 3) from SPIReS–MODIS, is perhaps the most intriguing result, contradicting the result of previous studies (Sect. 1) which find that finer-spatial-resolution estimates of snow cover reduce SWE errors. The reduction of 4 %–5 % relative bias from the two 30 m snow cover forcings compared to MODIS agrees with the previous findings, although the magnitudes of the reductions are smaller than in previous studies (e.g., Durand et al., 2008). To test if the lower MAEs from MODIS are resolution artifacts, the SPIReS–HLS and SCAG–Fusion products were also upscaled to 2000 m cell sizes instead of 120 m. For mean values over all water years for these upscaled comparisons (Table A1), SPIReS–MODIS still had the lowest relative MAE, but the SCAG–Fusion relative bias dropped to 3 %, while the SPIReS–HLS relative bias increased to 9 %, equal to SPIReS–MODIS. These results suggest that the evaluations are sensitive to the upscaled pixel size, meaning that the differences in errors across the three SWE reconstructions may be within uncertainty bounds introduced by upscaling artifacts such as basin delineation. For example, in the UTM 11N projection, the shapefile obtained from ASO Inc. for the Tuolumne River basin above Hetch Hetchy Reservoir has an area of 1175 km2; a raster of the basin at 120 m has an area of 1153 km2 (−1.8 %), while a raster at 2000 m has an area of 1132 km2 (−3.6 %). Even when using vector basin outlines, as suggested by Cline et al. (1998), these artifacts are inherent in the discretization of geospatial data and cannot be eliminated.
Another explanation for the poorer MAE performance from SPIReS–MODIS is that some spatial variation in topography is lost with the coarser resolution. To test this hypothesis, a semivariogram of the terrain slope is examined, as in Baba et al. (2019). The semivariogram shows a flattening around 500 m, indicating that variation in topography, which can manifest in topographically driven variables such as direct solar illumination, is poorly captured at MODIS and coarser spatial scales. This semivariogram analysis confirms the above hypothesis. Further, downscaling coarse-scale reanalysis products (Winstral et al., 2014), e.g., the downwelling radiation from Clouds and the Earth's Radiant Energy System (Rutan et al., 2015) at 1∘ spatial resolution, has inherent limitations, often due to clouds (Lapo et al., 2017). Important to note is that ParBal does not use precipitation as a forcing and thus does not suffer from well-known biases and downscaling issues (Raleigh et al., 2015; Pflug et al., 2021).
Alternatively, the lower MAE may indicate the importance of daily imaging from MODIS compared to the HLS snow cover, which had median gaps of 5 d between revisits after filtering for clouds. In contrast, the SCAG–Fusion used daily MODIS snow cover in the prediction and training steps indicating it suffers from errors not related to revisit time.
An example of the SWE modeled by ASO on 4 May 2019 and the three reconstructions is shown in Fig. 6. The spatial distribution of the SWE from ASO matches well with all the reconstructions. Differences between the reconstructions can be seen around Mount Lyell, at the southernmost part of the basin. The ASO SWE shows high variability here, ranging from a few hundred millimeters of SWE to over 2000 mm, while the reconstructions model consistently higher amounts of SWE. The overestimates here are likely related to false-positive classifications for snow. Especially late in the summer, when melt rates are high, these false positives can lead to substantial overestimates of SWE during reconstruction (Slater et al., 2013). A close examination of the mostly snow-free areas in gray shows that only the SPIReS–HLS reconstructions replicate the small patches of thin snow in this area, likely because the SPIReS–HLS snow cover was not smoothed to the same degree as SPIReS–MODIS or SCAG–Fusion, which both use heavy smoothing to reduce noise and smearing from MODIS.
Errors are further examined by date (Fig. 7 and Table A2). Except for 13 April 2020, the bias across all the products is between −20 % and 20 % (Fig. 7a). Figure 8 shows a snow pillow (weighing gauge, California Department of Water Resources station code DAN, elevation 2987 m) and that the ASO flight on 13 April 2020 is the only flight in this study that occurred prior to peak SWE. Overestimates of SWE prior to its peak are a limitation of SWE reconstruction. The hybrid SWE method (Sect. 2.4) extends SWE estimates throughout the year, but the high biases found on this date are not surprising because snowmelt occurred prior to the flight and snow accumulation occurred after the flight. Note the missing data from DAN after 3 May 2018 in Fig. 8, but the CUES snow pillow, which is nearby and at a similar elevation (2940 m), shows clear ablation during May 2018.
Examination of the per-pixel MAE (Fig. 7b) shows that the SPIReS–HLS product has the most consistent values, with the two approaches that used MODIS data (SPIReS–MODIS and SCAG–Fusion) showing more variability, perhaps again due to the smoothing needed for the relatively noisy MODIS data or the fact that SCAG–Fusion was trained using more data outside the test period (March 2013 to December 2017 and January 2021 to March 2022) than within it (January 2018 to January 2021).
Stillinger et al. (2023) show that errors in snow cover mapping depend on canopy cover, which has to do with how much areal snow is viewable at the pixel scale by a sensor, which affects the accuracy of the SWE reconstructions (Bair et al., 2016). Thus, we examine errors in the SWE reconstructions, binned by canopy cover fraction, for each snow cover forcing. The bin centered at 5 % (range: 0 % to 9.9 %) canopy cover (containing 46 %–60 % of pixels in the basin; Table A3) shows (Fig. 9ab) relatively unbiased errors with MAE values close to the means (Table 3), but SWE biases become positive with increasing canopy cover for SPIReS–MODIS yet negative for SCAG–Fusion and for SPIReS–HLS (except for the highest canopy fractions which contain only 5 % of the basin's pixel; Table A3).
The bias and MAE with increasing canopy cover for SPIReS–HLS and SCAG–Fusion SWE reconstructions are similar to errors in fsca from Landsat 8 (Fig. 4c; Stillinger et al., 2023). These fsca biases have similar shapes to the SWE biases indicating these fsca errors cause the SWE errors. Conversely, SPIReS–MODIS shows unbiased fsca with increasing canopy cover (Fig. 5d; Stillinger et al., 2023), indicating some other source of error in the SPIReS–MODIS SWE reconstructions.
In summary, the answer to the question posed by the title of this study is that basin-wide SWE is marginally more accurate with finer spatial resolution. Specifically, the bias – arguably the most important error statistic for water resource management – was 4 %–5 % lower using the finer-resolution snow cover forcings. However, the results are mixed relative to previous studies. For example, Durand et al. (2008) and Molotch and Margulis (2008) report both lower MAE and bias with a 30 m Landsat ETM+ snow cover forcing compared to snow cover from MODIS and AVHRR. The explanation for why some previous studies showed more significant improvements going from moderate- to high-resolution forcings may be the snow mapping algorithms used. An accurate technique for dealing with mixed pixels is particularly important for moderate-resolution sensors since for midlatitude mountains most pixels are mixed at 500 m (Selkowitz et al., 2014). In Durand et al. (2008) and Molotch and Margulis (2008), the finer-resolution Landsat ETM+ snow cover used a spectral unmixing technique (Painter et al., 2003), but the MODIS snow cover was based on the normalized snow difference technique, which only uses two bands versus all available for spectral unmixing and is shown to have higher MAE and bias (Stillinger et al., 2023). In Cline et al. (1998), the only other study to specifically examine spatial scale with SWE reconstruction, a spectral mixture technique was used on 30 m Landsat ETM+ to produce snow cover estimates (Rosenthal and Dozier, 1996). In that study, the coarsened results produced basin-wide SWE above and below the control simulation used as validation, suggesting that coarsening components of the energy balance did not show a clear trend in error. The snow cover used in that study is shown to have low bias and other measures of error from [0–1] fsca (Rosenthal and Dozier, 1996), thus reducing errors from mixed pixels. Increased spatial and temporal resolution through sensor design, fusion techniques, and satellite constellations is the future of Earth observations, but this study shows how a moderate-resolution sensor such as MODIS still offers value for snow mapping and modeling.
Optimal-resolution questions are fundamental to the global study of snow and will inform future scientific priorities and mission specifications. Increasing spatial and temporal resolution mark remote sensing achievements with the implicit assumption that finer resolution provides greater accuracy. To test this assumption for snow hydrology, an energy balance SWE reconstruction model was run at two different spatial resolutions using three different snow cover forcings. Contrary to previous work, the baseline case using SPIReS–MODIS, a daily 463 m product, showed a lower MAE – a measure of per-pixel accuracy – compared to SCAG–Fusion and SPIReS–HLS, both with 30 m spatial resolution. SPIReS–HLS showed the lowest bias; however, the differences in the errors between all three products may be within the uncertainty caused by scaling artifacts such as basin boundary delineation. The improved bias with increasing spatial resolution, arguably the most important measure for water management, is a promising result; however, the increased MAE with finer spatial resolution suggests that the daily acquisitions from MODIS with finer temporal resolution provide additional accuracy and/or that there are downscaling limitations with relatively coarse reanalysis data, e.g., 105 m (1∘) downscaled to 30 m. Improvements such as the inclusion of Landsat 9 and version 2.0 of the HLS data may improve some of the errors. Future satellite missions that leverage existing and planned constellations such as Landsat Next will improve revisit times, as gaps between observations are still an issue for the HLS data. In summary, conclusions are as follows: (1) spectrally unmixed snow cover and snow albedo from MODIS continue to provide accurate forcings for snow models and (2) increased spatial and temporal resolution through sensor design, fusion techniques, and satellite constellations are the future of Earth observations, but existing moderate-resolution sensors still offer value.
All data are in accessible repositories. SPIReS–MODIS: the snow cover is part of a daily Western US product covering WY 2001–2021 (https://doi.org/10.21424/R4H05T, Bair and Stillinger, 2022). The corresponding reconstructed SWE is in a Dryad Repository (https://doi.org/10.25349/D9TK7H, Bair, 2023a). SPIReS–HLS: the snow cover and reconstructions are in a Dryad repository (https://doi.org/10.25349/D9PW47, Bair, 2023b). SCAG–Fusion: the snow cover and reconstructions are in a Dryad repository (https://doi.org/10.25349/D9PW47, Bair, 2023b).
According to CRediT taxonomy: EHB – all 14 contributor roles; JD – conceptualization, writing (review and editing); KR – conceptualization, data curation, formal analysis, funding acquisition, methodology, resources, writing (review and editing); TS – investigation, writing (review and editing); WK – data curation, formal analysis, investigation, methodology, software, supervision; RED – resources, funding acquisition.
At least one of the (co-)authors is a member of the editorial board of The Cryosphere. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We thank Marie Dumont for editing as well as Simon Gascoin and one anonymous referee for their insightful reviews.
This research has been supported by the National Aeronautics and Space Administration (grant nos. 80NSSC21K0997, 80NSSC20K1722, 80NSSC20K1349, 80NSSC18K1489, 80NSSC21K0620, 80NSSC18K0427, 80NSSC20K1721, 80NSSC22K0703, and 80NSSC22K0929), and other support is from the Broad Agency Announcement Program and the Cold Regions Research and Engineering Laboratory (ERDC-CRREL) under contract no. W913E520C0019.
This paper was edited by Marie Dumont and reviewed by Simon Gascoin and one anonymous referee.
Baba, M. W., Gascoin, S., Kinnard, C., Marchane, A., and Hanich, L.: Effect of digital elevation model resolution on the simulation of the snow cover evolution in the High Atlas, Water Resour. Res., 55, 5360–5378, https://doi.org/10.1029/2018WR023789, 2019.
Bair, E. H.: SPIReS-MODIS-ParBal Snow Water Equivalent Reconstruction: Western USA, water years 2001–2021, Dryad [data set], https://doi.org/10.25349/D9TK7H, 2023a.
Bair, E. H.: Snow cover and snow water equivalent for “How do tradeoffs in satellite spatial and temporal resolution impact snow water equivalent reconstruction?”, Dryad [data set], https://doi.org/10.25349/D9PW47, 2023b.
Bair, E. H.: ParBal, Zenodo [code], https://doi.org/10.5281/zenodo.8106305, 2023c.
Bair, E. H.: SPIRES, Zenodo [code], https://doi.org/10.5281/zenodo.8106303, 2023d.
Bair, E. H. and Stillinger, T.: SPIReS: Western USA snow cover and snow surface properties, water years 2001–2021, CUES [data set], https://doi.org/10.21424/R4H05T, 2022.
Bair, E. H., Dozier, J., Davis, R. E., Colee, M. T., and Claffey, K. J.: CUES – A study site for measuring snowpack energy balance in the Sierra Nevada, Front. Earth Sci., 3, 58, https://doi.org/10.3389/feart.2015.00058, 2015.
Bair, E. H., Rittger, K., Davis, R. E., Painter, T. H., and Dozier, J.: Validating reconstruction of snow water equivalent in California's Sierra Nevada using measurements from the NASA Airborne Snow Observatory, Water Resour. Res., 52, 8437–8460, https://doi.org/10.1002/2016WR018704, 2016.
Bair, E. H., Abreu Calfa, A., Rittger, K., and Dozier, J.: Using machine learning for real-time estimates of snow water equivalent in the watersheds of Afghanistan, The Cryosphere, 12, 1579–1594, https://doi.org/10.5194/tc-12-1579-2018, 2018.
Bair, E. H., Rittger, K., Skiles, S. M., and Dozier, J.: An examination of snow albedo estimates from MODIS and their impact on snow water equivalent reconstruction, Water Resour. Res., 55, 7826–7842, https://doi.org/10.1029/2019wr024810, 2019.
Bair, E. H., Stillinger, T., and Dozier, J.: Snow Property Inversion from Remote Sensing (SPIReS): A generalized multispectral unmixing approach with examples from MODIS and Landsat 8 OLI, IEEE T. Geosci. Remote, 59, 7270–7284, https://doi.org/10.1109/TGRS.2020.3040328, 2021.
Bair, E. H., Dozier, J., Stern, C., LeWinter, A., Rittger, K., Savagian, A., Stillinger, T., and Davis, R. E.: Divergence of apparent and intrinsic snow albedo over a season at a sub-alpine site with implications for remote sensing, The Cryosphere, 16, 1765–1778, https://doi.org/10.5194/tc-16-1765-2022, 2022.
Baumgartner, M. F., Seidel, K., and Martinec, J.: Toward snowmelt runoff forecast based on multisensor remote-sensing informnation, IEEE T. Geosci. Remote, GE-25, 746–750, https://doi.org/10.1109/TGRS.1987.289744, 1987.
Blöschl, G.: Scaling issues in snow hydrology, Hydrol. Process., 13, 2149–2175, https://doi.org/10.1002/(SICI)1099-1085(199910)13:14/15<2149::AID-HYP847>3.0.CO;2-8, 1999.
Bouamri, H., Kinnard, C., Boudhar, A., Gascoin, S., Hanich, L., and Chehbouni, A.: MODIS does not capture the spatial heterogeneity of snow cover induced by solar radiation, Front. Earth Sci., 9, 640250, https://doi.org/10.3389/feart.2021.640250, 2021.
Claverie, M., Ju, J., Masek, J. G., Dungan, J. L., Vermote, E. F., Roger, J. C., Skakun, S. V., and Justice, C.: The harmonized Landsat and Sentinel-2 surface reflectance data set, Remote Sens. Environ., 219, 145–161, https://doi.org/10.1016/j.rse.2018.09.002, 2018.
Cline, D., Elder, K., and Bales, R.: Scale effects in a distributed snow water equivalence and snowmelt model for mountain basins, Hydrol. Process., 12, 1527–1536, https://doi.org/10.1002/(SICI)1099-1085(199808/09)12:10/11<1527::AID-HYP678>3.0.CO;2-E, 1998.
Conner, M. M., Stephenson, T. R., German, D. W., Monteith, K. L., Few, A. P., and Bair, E. H.: Survival analysis: Informing recovery of Sierra Nevada bighorn sheep, J. Wildlife Manage., 82, 1442–1458, https://doi.org/10.1002/jwmg.21490, 2018.
Dozier, J., Painter, T. H., Rittger, K., and Frew, J. E.: Time-space continuity of daily maps of fractional snow cover and albedo from MODIS, Adv. Water Resour., 31, 1515–1526, https://doi.org/10.1016/j.advwatres.2008.08.011, 2008.
Durand, M., Molotch, N. P., and Margulis, S. A.: Merging complementary remote sensing datasets in the context of snow water equivalent reconstruction, Remote Sens. Environ., 112, 1212–1225, https://doi.org/10.1016/j.rse.2007.08.010, 2008.
Foga, S., Scaramuzza, P. L., Guo, S., Zhu, Z., Dilley, R. D., Beckmann, T., Schmidt, G. L., Dwyer, J. L., Joseph Hughes, M., and Laue, B.: Cloud detection algorithm comparison and validation for operational Landsat data products, Remote Sens. Environ., 194, 379–390, https://doi.org/10.1016/j.rse.2017.03.026, 2017.
Global Modeling and Assimilation Office (GMAO): tavg1_2d_flx_Nx: MERRA-2 1-Hourly, Time-Averaged, Single-Level, Assimilation, Surface Flux Diagnostics version 5.12.4, GES DISC, https://doi.org/10.5067/7MCPBJ41Y0K6, 2015.
Hall, D. K., Riggs, G. A., Salomonson, V. V., DiGirolamo, N. E., and Bayr, K. J.: MODIS snow-cover products, Remote Sens. Environ., 83, 181–194, https://doi.org/10.1016/S0034-4257(02)00095-0, 2002.
Lapo, K. E., Hinkelman, L. M., Sumargo, E., Hughes, M., and Lundquist, J. D.: A critical evaluation of modeled solar irradiance over California for hydrologic and land surface modeling, J. Geophys. Res.-Atmos., 122, 299–317, https://doi.org/10.1002/2016JD025527, 2017.
Li, J. and Roy, D. P.: A global analysis of Sentinel-2A, Sentinel-2B and Landsat-8 data revisit intervals and implications for terrestrial monitoring, Remote Sensing, 9, 902, https://doi.org/10.3390/rs9090902, 2017.
Liston, G. E., Haehnel, R. B., Sturm, M., Hiemstra, C. A., Berezovskaya, S., and Tabler, R. D.: Instruments and methods: Simulating complex snow distributions in windy environments using SnowTran-3D, J. Glaciol., 53, 241–256, https://doi.org/10.3189/172756507782202865, 2007.
Luce, C. H., Tarboton, D. G., and Cooley, K. R.: The influence of the spatial distribution of snow on basin-averaged snowmelt, Hydrol. Process., 12, 1671–1683, https://doi.org/10.1002/(SICI)1099-1085(199808/09)12:10/11<1671::AID-HYP688>3.0.CO;2-N, 1998.
MathWorks: MATLAB Mapping Toolbox 5.4: User's Guide, The MathWorks, Natick, MA, 814 pp., 2022.
Molotch, N. P. and Margulis, S. A.: Estimating the distribution of snow water equivalent using remotely sensed snow cover data and a spatially distributed snowmelt model: A multi-resolution, multi-sensor comparison, Adv. Water Resour., 31, 1503–1514, https://doi.org/10.1016/j.advwatres.2008.07.017, 2008.
Painter, T. H., Dozier, J., Roberts, D. A., Davis, R. E., and Green, R. O.: Retrieval of subpixel snow-covered area and grain size from imaging spectrometer data, Remote Sens. Environ., 85, 64–77, https://doi.org/10.1016/S0034-4257(02)00187-6, 2003.
Painter, T. H., Rittger, K., McKenzie, C., Slaughter, P., Davis, R. E., and Dozier, J.: Retrieval of subpixel snow-covered area, grain size, and albedo from MODIS, Remote Sens. Environ., 113, 868–879, https://doi.org/10.1016/j.rse.2009.01.001, 2009.
Painter, T. H., Bryant, A. C., and Skiles, S. M.: Radiative forcing by light absorbing impurities in snow from MODIS surface reflectance data, Geophys. Res. Lett., 39, L17502, https://doi.org/10.1029/2012GL052457, 2012.
Painter, T. H., Berisford, D. F., Boardman, J. W., Bormann, K. J., Deems, J. S., Gehrke, F., Hedrick, A., Joyce, M., Laidlaw, R., Marks, D., Mattmann, C., McGurk, B., Ramirez, P., Richardson, M., Skiles, S. M., Seidel, F. C., and Winstral, A.: The Airborne Snow Observatory: Fusion of scanning lidar, imaging spectrometer, and physically-based modeling for mapping snow water equivalent and snow albedo, Remote Sens. Environ., 184, 139–152, https://doi.org/10.1016/j.rse.2016.06.018, 2016.
Pflug, J. M., Hughes, M., and Lundquist, J. D.: Downscaling snow deposition using historic snow depth patterns: diagnosing limitations from snowfall biases, winter snow losses, and interannual snow pattern repeatability, Water Resour. Res., 57, e2021WR029999, https://doi.org/10.1029/2021WR029999, 2021.
Raleigh, M. S., Lundquist, J. D., and Clark, M. P.: Exploring the impact of forcing error characteristics on physically based snow simulations within a global sensitivity analysis framework, Hydrol. Earth Syst. Sci., 19, 3153–3179, https://doi.org/10.5194/hess-19-3153-2015, 2015.
Rittger, K., Painter, T. H., and Dozier, J.: Assessment of methods for mapping snow cover from MODIS, Adv. Water Resour., 51, 367–380, https://doi.org/10.1016/j.advwatres.2012.03.002, 2013.
Rittger, K., Bair, E. H., Kahl, A., and Dozier, J.: Spatial estimates of snow water equivalent from reconstruction, Adv. Water Resour., 94, 345–363, https://doi.org/10.1016/j.advwatres.2016.05.015, 2016.
Rittger, K., Raleigh, M. S., Dozier, J., Hill, A. F., Lutz, J. A., and Painter, T. H.: Canopy adjustment and improved cloud detection for remotely sensed snow cover mapping, Water Resour. Res., 56, e2019WR024914, https://doi.org/10.1029/2019WR024914, 2020.
Rittger, K., Krock, M., Kleiber, W., Bair, E. H., Brodzik, M. J., Stephenson, T. R., Rajagopalan, B., Bormann, K. J., and Painter, T. H.: Multi-sensor fusion using random forests for daily fractional snow cover at 30 m, Remote Sens. Environ., 264, 112608, https://doi.org/10.1016/j.rse.2021.112608, 2021.
Rodell, M., Houser, P. R., Jambor, U., Gottschalck, J., Mitchell, K., Meng, C. J., Arsenault, K., Cosgrove, B., Radakovich, J., Bosilovich, M., Entin, J. K., Walker, J. P., Lohmann, D., and Toll, D.: The Global Land Data Assimilation System, B. Am. Meteorol. Soc., 85, 381–394, https://doi.org/10.1175/BAMS-85-3-381, 2004.
Rosenthal, W. and Dozier, J.: Automated mapping of montane snow cover at subpixel resolution from the Landsat Thematic Mapper, Water Resour. Res., 32, 115–130, https://doi.org/10.1029/95WR02718, 1996.
Rutan, D. A., Kato, S., Doelling, D. R., Rose, F. G., Nguyen, L. T., Caldwell, T. E., and Loeb, N. G.: CERES synoptic product: Methodology and validation of surface radiant flux, J. Atmos. Ocean. Tech., 32, 1121–1143, https://doi.org/10.1175/JTECH-D-14-00165.1, 2015.
Schlögl, S., Marty, C., Bavay, M., and Lehning, M.: Sensitivity of Alpine3D modeled snow cover to modifications in DEM resolution, station coverage and meteorological input quantities, Environ. Modell. Softw., 83, 387–396, https://doi.org/10.1016/j.envsoft.2016.02.017, 2016.
Selkowitz, D. J., Forster, R. R., and Caldwell, M. K.: Prevalence of pure versus mixed snow cover pixels across spatial resolutions in alpine environments, Remote Sens., 6, 12478–12508, https://doi.org/10.3390/rs61212478, 2014.
Slater, A. G., Barrett, A. P., Clark, M. P., Lundquist, J. D., and Raleigh, M. S.: Uncertainty in seasonal snow reconstruction: Relative impacts of model forcing and image availability, Adv. Water Resour., 55, 165–177, https://doi.org/10.1016/j.advwatres.2012.07.006, 2013.
Stillinger, T., Rittger, K., Raleigh, M. S., Michell, A., Davis, R. E., and Bair, E. H.: Landsat, MODIS, and VIIRS snow cover mapping algorithm performance as validated by airborne lidar datasets, The Cryosphere, 17, 567–590, https://doi.org/10.5194/tc-17-567-2023, 2023.
Stillinger, T. C.: Observing Snow from Space: Snow/Cloud Discrimination and Opportunities in Water Supply Forecasting, University of California, Santa Barbara, United States – California, 152 pp., 2019.
Storey, J., Roy, D. P., Masek, J., Gascon, F., Dwyer, J., and Choate, M.: A note on the temporary misregistration of Landsat-8 Operational Land Imager (OLI) and Sentinel-2 Multi Spectral Instrument (MSI) imagery, Remote Sens. Environ., 186, 121–122, https://doi.org/10.1016/j.rse.2016.08.025, 2016.
Tan, B., Woodcock, C. E., Hu, J., Zhang, P., Ozdogan, M., Huang, D., Yang, W., Knyazikhin, Y., and Myneni, R. B.: The impact of gridding artifacts on the local spatial properties of MODIS data: Implications for validation, compositing, and band-to-band registration across resolutions, Remote Sens. Environ., 105, 98–114, https://doi.org/10.1016/j.rse.2006.06.008, 2006.
Turpin, O. C., Caves, R. G., Ferguson, R. I., and Johansson, B.: Verification of simulated snow cover in an Arctic basin using satellite-derived snow-cover maps, Ann. Glaciol., 31, 391–396, https://doi.org/10.3189/172756400781819932, 2000.
USGS: Collection 2 Landsat 8-9 OLI (Operational Land Imager) and TIRS (Thermal Infrared Sensor) Level-2 Science Product, USGS EROS, https://doi.org/10.5066/P9OGBGM6, 2021.
Vermote, E. and Wolfe, R. E., MOD09GA MODIS/Terra Surface Reflectance Daily L2G Global 1kmand 500m SIN Grid V006, NASA LP DAAC, https://doi.org/10.5067/MODIS/MOD09.006, 2015.
Vermote, E., Justice, C., Claverie, M., and Franch, B.: Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product, Remote Sens. Environ., 185, 46–56, https://doi.org/10.1016/j.rse.2016.04.008, 2016.
Vuyovich, C. M., Deeb, E. J., Polashenski, C., Courville, Z., Hiemstra, C. A., Wagner, A. M., Eylander, J. B., and Davis, R. E.: Snow Strategic Science Plan, ERDC/CRREL, Hanover, NH, 91, https://doi.org/10.21079/11681/29554, 2018.
Willmott, C. J. and Matsuura, K.: Advantages of the mean absolute error (MAE) over the root mean square error (RMSE) in assessing average model performance, Clim. Res., 30, 79–82, https://doi.org/10.3354/cr030079, 2005.
Winstral, A., Marks, D., and Gurney, R.: Assessing the sensitivities of a distributed snow model to forcing data resolution, J. Hydrometeorol., 15, 1366–1383, https://doi.org/10.1175/jhm-d-13-0169.1, 2014.
Wolfe, R. E., Roy, D. P., and Vermote, E.: MODIS land data storage, gridding, and compositing methodology: Level 2 grid, IEEE T. Geosci. Remote, 36, 1324–1338, https://doi.org/10.1109/36.701082, 1998.
Xia, Y., Mitchell, K., Ek, M., Sheffield, J., Cosgrove, B., Wood, E., Luo, L., Alonge, C., Wei, H., Meng, J., Livneh, B., Lettenmaier, D., Koren, V., Duan, Q., Mo, K., Fan, Y., and Mocko, D.: Continental-scale water and energy flux analysis and validation for the North American Land Data Assimilation System project phase 2 (NLDAS-2): 1. Intercomparison and application of model products, J. Geophys. Res.-Atmos., 117, D03109, https://doi.org/10.1029/2011JD016048, 2012.