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.

1 Introduction

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.

2 Approach

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).

2.1 SPIReS–MODIS

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.

2.2 SPIReS–HLS

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 $\mathrm{1119}\times \mathrm{1297}\times \mathrm{7}\times$ 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.

Figure 1 Revisit times covering the Tuolumne River basin above Hetch Hetchy Reservoir for the HLS S30 and L30 combined.

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.

2.3 SCAG–Fusion

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.

Figure 2 Cumulative number of images over 8 years acquired from Landsat 8. Lighter to darker lines indicate increasing cloud coverage from 10 %–50 %. Numbers after the cloud cover percentage in the legend correspond to the total number of images. Images with 12 % or less cloud cover were selected.

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, SWE_GLDAS, 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

$$\begin{array}{}\text{(1)}& \begin{array}{rl}t& =\left({\mathrm{SWE}}_{\mathrm{GLDAS},\mathrm{\Delta }{t}_{\mathrm{1}}}>\mathrm{0},{\mathrm{fsca}}_{\mathrm{\Delta }{t}_{\mathrm{1}}}>\mathrm{0}\right)\\ & \mathrm{\&}\left({\mathrm{SWE}}_{\mathrm{GLDAS},\mathrm{\Delta }{t}_{\mathrm{2}}}=\mathrm{0},{\mathrm{fsca}}_{\mathrm{\Delta }{t}_{\mathrm{2}}}=\mathrm{0}\right)\\ {\mathrm{SWE}}_{\mathrm{GLDAS}}^{\ast }& ={\mathrm{SWE}}_{\mathrm{GLDAS},t},\end{array}\end{array}$$

where the asterisk denotes the selected pixels and fsca is from the fine-scaled product, i.e., Sect. 2.1–2.2. The Δt₁ and Δt₂ indicate different time periods (days). Because there can be multiple pixels with matching snow cover duration, the daily mean of ${\mathrm{SWE}}_{\mathrm{GLDAS}}^{\ast }$ is taken. The maximum value of that daily mean and its index i_max are computed. A scaling coefficient c is calculated as

$$\begin{array}{}\text{(2)}& c={\mathrm{SWE}}_{\mathrm{ParBal},i\max }/max\left(\stackrel{\mathrm{‾}}{{\mathrm{SWE}}_{\mathrm{GLDAS}}^{\ast }}\right),\end{array}$$

where SWE_ParBal,imax is the value of the reconstructed SWE from ParBal (Eq. 2 in Bair et al., 2016) at the time indexed by i_max and the overbar denotes an average. The following case can arise:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}c& =\mathrm{0},{\mathrm{SWE}}_{\mathrm{ParBal},i_{\max }}=\mathrm{0},\left({\mathrm{SWE}}_{\mathrm{GLDAS}}>\mathrm{0}\right)\\ & \mathrm{\&}\left(\mathrm{fsca}>\mathrm{0}\right).\end{array}\end{array}$$

For example, this case can occur when ParBal models all the mass loss via sublimation. For this case, SWE_GLDAS is used when fsca >0. Otherwise, the hybrid SWE prior to the peak is set at day i as

$$\begin{array}{}\text{(4)}& {\mathrm{SWE}}_{\mathrm{hybrid},i}=c\times {\mathrm{SWE}}_{\mathrm{GLDAS},i},i\le i\max .\end{array}$$

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 i_max 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.

Figure 3 Hybrid SWE estimates for the accumulation season combining reconstruction and GLDAS for an example pixel using the SPIReS–HLS snow cover.

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.

Table 1 Tuolumne River basin above Hetch Hetchy Reservoir SWE estimates for 2018–2020 for the Airborne Snow Observatory.

2.6 Analysis

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:

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathrm{bias}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{j=\mathrm{1}}^N{\mathrm{SWE}}_{\mathrm{ParBal},j}-{\mathrm{SWE}}_{\mathrm{ASO},j},\text{(6)}& {\displaystyle }\mathrm{relative}\mathrm{bias}={\displaystyle \frac{\frac{\mathrm{1}}{N}{\sum }_{j=\mathrm{1}}^N{\mathrm{SWE}}_{\mathrm{ParBal},j}-{\mathrm{SWE}}_{\mathrm{ASO},j}}{\frac{\mathrm{1}}{N}{\sum }_{j=\mathrm{1}}^N{\mathrm{SWE}}_{\mathrm{ASO},j}}},\text{(7)}& {\displaystyle }\mathrm{mean}\mathrm{absolute}\mathrm{error}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{j=\mathrm{1}}^N\left|{\mathrm{SWE}}_{\mathrm{ParBal},j}-{\mathrm{SWE}}_{\mathrm{ASO},j}\right|,\text{(8)}& {\displaystyle }\begin{array}{rl}\mathrm{relative}& \mathrm{mean}\mathrm{absolute}\mathrm{error}=\\ & {\displaystyle \frac{\frac{\mathrm{1}}{N}\sum _{j=\mathrm{1}}^N\left|{\mathrm{SWE}}_{\mathrm{ParBal},j}-{\mathrm{SWE}}_{\mathrm{ASO},j}\right|}{\frac{\mathrm{1}}{N}\sum _{j=\mathrm{1}}^N{\mathrm{SWE}}_{\mathrm{ASO},j}}},\end{array}\end{array}$$

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).

Table 2 Snow albedo from SPIReS–MODIS validated with terrain-corrected snow albedo from the CUES site on Mammoth Mountain (Bair et al., 2015) taken using an adjustable arm to keep the radiometers 1 m above the snow surface (Bair et al., 2022). A best of 3×3 pixel neighborhood was used to account for geolocational uncertainty.

Figure 4 Semivariance of terrain slope of the Tuolumne River basin above Hetch Hetchy. The slope of the semivariance (not the terrain slope itself) shows a flattening around 500 m or about the MODIS pixel size.

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.

Figure 5 Basin-wide SWE values by date for the three SWE reconstructions (a–c) compared to Airborne Snow Observatory estimates. Assumed uncertainty in the ASO measurements is ±4 % (Sect. 2.5) and is shaded in gray.

Table 3 Error statistics by year for the three SWE reconstructions. Mean values for all years are shown in bold. More detailed errors by date are given in the Appendix.

3 Results and discussion

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 km²; a raster of the basin at 120 m has an area of 1153 km² (−1.8 %), while a raster at 2000 m has an area of 1132 km² (−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.

Figure 6 SWE in the Tuolumne River basin above Hetch Hetchy Reservoir for 4 May 2019 modeled by the Airborne Snow Observatory (a) along with reconstructions from SPIReS–MODIS (b), SPIReS–HLS (c), and SCAG–Fusion (d).

Figure 7 SWE errors by date for the three SWE reconstructions. The relative bias is shown in (a) and the relative MAE in (b).

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).

Figure 8 Snow pillow DAN in the Tuolumne River basin showing daily SWE. The X markers show the dates of ASO flights. Circled is the 13 April 2020 ASO flight, which is the only flight that occurred prior to peak SWE. The pillow was not reported from 3 May to 21 November 2018, but a nearby snow pillow shows consistent ablation in May 2018.

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).

Figure 9 SWE errors for all dates for the three SWE reconstructions binned by canopy cover percent. Labeled are the bin centers. The relative bias is shown in (a) and the relative MAE in (b).

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.

4 Conclusions

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., 10⁵ 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.

Appendix A

Table A1 Error statistics by date for the three SWE reconstructions but with all pixels upscaled to 2000 m. The SPIRES–MODIS rows are identical to those in Table 2 and are shown for comparison. Mean values for all years are shown in bold.

Table A2 Error statistics by date for the three SWE reconstructions.

Table A3 Error statistics by canopy cover, for all dates, for the three SWE reconstructions.