Snow albedo is a dominant control on snowmelt in many parts of the world. An empirical albedo decay equation, developed over 60 years ago, is still used in snowmelt models. Several empirical snow albedo models developed since show wide spread in results. Remotely sensed snow albedos have been used in a few studies, but validations are scarce because of the difficulty in making accurate in situ measurements. Reconstruction of snow water equivalent (SWE), where the snowpack is built in reverse, is especially sensitive to albedo. We present two new contributions: (1) an updated albedo model where grain size and light absorbing particle content are solved for simultaneously and (2) multiyear comparisons of remotely sensed and in situ albedo measurements from three high‐altitude sites in the western United States. Our remotely sensed albedos show 4 to 6% RMSE and negligible bias. In comparison, empirical albedo decay models, which require extensive in situ measurements, show RMSE values of 7 to 17% with biases of −6 to −14%. We examine the sensitivity of SWE reconstructions to albedo error at two sites. With no simulated error in albedo, reconstructed SWE had MAE values of 7 to 13% and 5–6% bias. The accuracy actually improved with some simulated added error, likely because of a fundamental bias in the reconstruction approach. Conversely, the best age‐based decay model showed an 18–20% MAE and bias in reconstructed SWE. We conclude that remotely sensed albedos where available are superior to age‐based approaches in all aspects except simplicity.

Even modern snow models (e.g., Anderson, ; Guan et al., ; Hedrick et al., ; Lehning et al., ) still use age‐based albedo decay equations to estimate snow albedo. This approach can be traced back to a publication over 60 years old (U.S. Army Corps of Engineers, ). Although simple to use, these empirical models require precise knowledge of when the last snowfall ended and do not account for topographic variability in the processes that change albedo, so they cannot be used in many parts of the world, and they have been shown to have RMSE values of 7 to 17% (Bair et al., , Bair, Davis, et al., ). Because snow albedo changes spatially and temporally depending on local variability in the energy balance (Molotch et al., ) and because albedo degradation caused by small, highly light‐absorbing particles (LAP) occurs independently of local processes in the snowpack (Skiles et al., ), we hypothesize that remotely sensed measurements of snow albedo, when and where available, are inherently more accurate than models based on aging.

Since 1956, radiative transfer approaches have been used to estimate the spectral albedo of both clean and dirty snow (reviewed by Warren, ). Advances in measurement of the absorption coefficients of ice (Picard et al., ; Warren & Brandt, ), dust (Haywood et al., ), and black carbon (Bond & Bergstrom, ) in the snowpack have led to more a more robust understanding of the spectral albedo of snow (Dang et al., ; Flanner & Zender, ). In addition to airborne dust and soot, snow algae (Painter et al., ), char left over from fires (Gleason & Nolin, ), and plant litter (Melloh et al., ) reduce snow albedo, particularly in the visible wavelengths.

Because running a full radiative transfer model to estimate snow albedo over large regions is impractical, statistical fits to the radiative transfer solutions have been used to estimate the broadband snow albedo *α* (Dang et al., ; Dozier et al., ; Gardner & Sharp, ; Marshall & Warren, ) which we define as*α*_{λ} is the spectral albedo of snow; *S*_{λ} is the solar spectral irradiance at the snow surface, both at wavelength *λ*; and *x*_{1} = 0.28 μm and *x*_{2} = 4.00 μm.

The spectral albedo *α*_{λ} depends only on the snow properties and not on the illumination itself, but the broadband albedo is a convolution of the spectral albedo and the spectrum of the incident radiation. The value of *α*_{λ} depends on optically equivalent grain radius *r*, optical properties and concentration of LAP, and illumination angle, while the magnitude, spectral distribution, and direct versus diffuse portions of *S*_{λ} depend on the path length through the atmosphere and cloud and atmospheric properties. While the snow community universally recognizes that illumination angle *θ*, *r*, and LAP affect *α*, the effects of the solar spectrum and atmospheric model on *α* are often ignored. Even when using the same snow reflectance model, significant differences in *α* result from using different atmospheric models and parameters (Figure and Table ). Therefore, any statistical fit to a snow albedo model is only useful if the atmospheric conditions where the model is applied are similar to those for which the model was developed. In Figure , we show *α* versus *r* for five models using an absolute air mass of 1.5, equivalent to the insolation‐weighted cosine of the sunlit portion of Earth (Cronin, ):

This value is conveniently close to the angle (49.50°) at which the direct and diffuse broadband snow albedos intersect (Wiscombe & Warren, ).

Remotely sensed observations from satellite or aerial platforms offer an ability to estimate snow albedo over large regions, but at least with sensors such as MODIS or VIIRS on polar‐orbiting satellites, retrievals provide at best a daily or twice daily snapshot. Multispectral sensors on geostationary satellites can image Earth at 1‐km spatial and 15‐min temporal resolution, but coverage is not global, and until the latest‐generation sensors on GOES‐16/17, Himawari‐8/9, and Fengyun‐4 were launched, the spectral and radiometric characteristics of the bands were inadequate for spectroscopic determination of the surface properties of snow. Thus, snow albedo models are still needed for simulations and to interpolate between measurement periods.

Because snow is so bright in the visible wavelengths, a small relative change in its albedo leads to a much larger relative change in its absorption, 1 − *α*, also known as *co‐albedo*. For example, a 5% decrease in *α*, say from 85% to 80%, represents about a 6% relative decrease (5%/85%). However, the relative change in absorption is 33% (5%/15%). These increases in absorption significantly affect the melt rate of the snowpack.

Reconstruction of snow water equivalent (SWE), a technique where the snowpack is built up in reverse from melt‐out to peak (Martinec & Rango, ), is particularly sensitive to albedo because the snowpack is built backward during the ablation phase, which is controlled by net radiation *R* in many parts of the world (Marks & Dozier, ; Mazurkiewicz et al., ; Oerlemans, ). The amount of downwelling shortwave radiation *S*_{↓} absorbed by the snowpack is a large part of net radiation *R*. For a ripe snowpack, this energy balance can be expressed as

*M*_{p} is the potential energy available for melt, *H* is the sensible heat flux, and *L* is the latent heat flux. The latent and sensible terms tend to cancel each other (Marks & Dozier, ) so we are left with net radiation *R* controlling melt, which is expressed as

*I*_{↓} is incoming longwave radiation and *I*_{↑} is emitted longwave radiation. Incoming and outgoing longwave do not cancel each other and are important terms in the energy balance during melt. Because most snow exists near its melting point (273.15 K), its longwave losses are usually greater than downward emission from the atmosphere; therefore, the magnitude of *I*_{↓} is usually smaller than that of *I*_{↑}, which causes the snowpack to lose energy in the absence of solar input (i.e., at night) and when *H*+*L* ≈ 0. Thus, in this simplified energy balance example, one can see how absorbed solar radiation *S*_{↓}(1 − *α*) is an important factor in snowpack melt.

In addition to *α*, SWE reconstruction is also sensitive to the fractional snow‐covered area *f*_{SCA}, as potential melt is distributed spatially as

*M* is an estimate of the snow or ice melt across a pixel. Although there are improvements to be made in remote sensing of *f*_{SCA}, considerable effort has gone into this problem, with demonstrated success. For example, the MODIS Snow Covered Area and Grain Size (MODSCAG) model (Painter et al., ), which we have used for *f*_{SCA} in our reconstructed SWE (Bair et al., ; Rittger et al., ), shows an RMSE of 0.09 and an F‐score (combining precision and recall) of 0.95 (Rittger et al., ).

For our remotely sensed *α* estimates, we start with the MODSCAG grain size (Painter et al., ), an estimate of *r* using spectral unmixing. Using a library of end‐members for clean snow reflectance, soil/rock, vegetation, and photometric shade, the MODSCAG algorithm uses a least squares fit to the MODIS BRDF retrievals from the MOD09GA surface reflectance product. Pixels are assumed to contain linear mixtures of the end‐members. The MODSCAG solutions contain fractional estimates of each of the end‐members and *r*.

MODSCAG has been shown to produce accurate *f*_{SCA} estimates (Rittger et al., ), but its *r* validation has been limited. Using only 12 manual measurements with a hand lens as validation, Painter et al. () show an RMSE of 64 μm and a bias of 30 μm.

For estimates of the LAP content and its effect on *α*, we use MODDRFS, the MODIS Dust and Radiative Forcing on Snow algorithm (Painter et al., ). Specifically, we use ΔVIS, the degradation of visible albedo by dust:

*α*_{vis,clean} and *α*_{vis,dirty} are the clean and dirty snow albedos in the visible part of the solar spectrum. Validation of the MODDRFS product has been limited. After bias correction, Painter et al. () report an RMSE of 33 W/m^{2} for the radiative forcing product (ΔVIS × *S*_{↓}), which translates to about a 6.4% RMSE in the visible albedo, validated with six years (2005–2011) of measurements at two study plots in the San Juan Mountains of Colorado. Their analysis used images with MODIS sensor zenith angles within 30° of nadir, which leads to less skew and other viewing geometry problems (Tan et al., ; Xiaoxiong et al., ). Our previous work (Bair et al., , Bair, Calfa, et al., ; Rittger et al., ) used these products for SWE reconstructions.

MODIS is a push‐broom scanner instrument with a 2,300‐km swath, so the sensor view zenith angle changes day‐to‐day for a given pixel, causing substantial skew and geolocational errors of up to one pixel (463 m) when the swath data are projected onto the sinusoidal grid (Tan et al., ; Xiaoxiong et al., ). Moreover, the pixel size (the ground instantaneous field of view) increases by a factor of 10 from nadir to the edge of the swath (Dozier et al., ). In addition, the reflective and thermal infrared bands on MODIS do not see through clouds, so the clouds must be identified and masked.

We filter, interpolate, and smooth the grain size and ΔVIS retrievals simultaneously with *f*_{SCA} processing (Dozier et al., ). The raw MODSCAG/MODDRFS retrievals are assembled into time‐space cubes. Pixels with gaps due to clouds, striping, and other errors are flagged. A persistence filter is applied to each pixel such that snow must be identified at a minimum number of days (e.g., 15 days) over a given time period (e.g., 90 days); otherwise, that pixel is set to *f*_{SCA} = 0 and all snow cover properties are set to null for those days. This persistence filter helps eliminate false positive snow identifications, which are caused by clouds and other bright surfaces. These constitute the *filtered* images. The flagged pixels in the filtered images are then interpolated using a weighted spline temporal interpolation, with the inverse of the sensor zenith angle for the weights (Dozier et al., ). We refer to these images as *interpolated*. The last step is a three‐dimensional (*x*, *y*, time) Gaussian smoothing for additional noise reduction, which mainly results from the day‐to‐day differences in the location and size of pixels caused by variation in the MODIS view angle. We refer to these images as *smoothed*.

The sites (Figure ) used were the Cold Regions Research and Engineering Laboratory and University of California Santa Barbara Energy site (CUES) on Mammoth Mountain, CA, USA (Bair et al., ; Bair, Davis, et al., ), and the Swamp Angel Study Plot (SASP) and Senator Beck Study Plot (SBSP) in the Senator Beck Basin, CO, USA (Landry et al., ; Skiles, ). These high‐altitude energy balance sites were selected for their long records of high‐quality energy balance measurements, in particular terrain‐corrected snow albedo (section ), partitioning of the illumination and albedo into broadband and near‐infrared components, and separating diffuse from direct irradiance. The combination of broadband and near‐infrared albedo enables separation of the effects of grain size and LAP on snow albedo.

The CUES and SASP sites are comparable in that both are subalpine, but CUES is located in a predominately maritime snow climate, with some years showing intermountain character (Bair, ), while SASP is located in a predominately continental snow climate (Landry et al., ). Of note is that the high elevation at CUES (2,940 m), compared to other maritime sites, makes midwinter rain infrequent. The SBSP alpine site is at high altitude (3,714 m), which makes it excellent for studying the snowpack energy balance but also makes access difficult so fewer snowpit measurements (i.e., SWE for the purposes of this paper) are available there.

For this study, vegetation neighboring the CUES instruments can shadow the snow beneath the downlooking radiometers, so the data must be cleaned to eliminate incorrect values (Bair, Davis, et al., ). At SASP the canopy is considerably more open, and SBSP is devoid of trees, so shadows are less of a problem at these sites, but the albedo measurements show other artifacts caused by, for example, nonsnow objects in the radiometers' field of view. However, the MODIS Terra overpass occurs at about 10:30 am, which is close to solar noon and when our albedo inversions (section ) are done, so problems with shadows are minimized.

At all three sites, the snow surface is neither level nor flat, meaning that illumination angles at the snow surface are different than those at the uplooking radiometers that measure incoming solar radiation. To correct for this problem, a plane is fit to the snow, and its local solar zenith angle ? is calculated, which can then be used to estimate a correction factor *c* to match the illumination conditions of the snow surface (Bair et al., ; Bair, Davis, et al., ; Painter et al., ):

*θ* is the local solar zenith angle and *θ*_{0} is the solar zenith angle for a level surface. Then the correction factor is applied to the downwelling direct radiation *B*_{↓}, while the downwelling diffuse radiation *D*_{↓} is not corrected because it does not depend on illumination:

*D*_{↑} is the reflected radiation from the snow surface measured by a downlooking radiometer. At CUES, *B*_{↓} and *D*_{↓} are measured automatically by a Delta‐T SPN1 Sunshine Pyranometer (Bair, Davis, et al., ), while at SASP and SBSP, diffuse and direct radiation are partitioned by applying potential irradiance values from SBDART to the measured broadband global irradiance. The snow surface correction for the near‐IR albedos does not include the diffuse irradiance term, as there is little atmospheric scattering in the near‐IR:

*α*_{nir} is the near‐IR albedo and *S*_{↓,nir} and *D*_{↑,nir} are the incoming and reflected near‐IR irradiances.

At CUES, hourly scans from a terrestrial laser scanner are used to map the snow surface, while at SASP and SBSP, manual snow depth measurements are recorded at multiple stakes with known locations. Both techniques are used to fit planes to the snow surface, from which the local ? is computed. While both approaches have pitfalls, we expect the snow surface correction at CUES to be more precise because measurements of the snow surface are taken every hour, whereas at SASP and SBSP, the manual measurements are taken less frequently.

Additional measurements at more sites would be useful for validating our approach. However, for example the data from Col de Porte in the French Alps (Lejeune et al., ) explicitly do not include the terrain correction, and the data set assembled to test snow models (Ménard et al., ) apparently does not implement this procedure needed to prevent artifacts in albedo measurements.

The snow properties themselves and the illumination angle control the spectral albedo of snow, *α*_{λ} in equation , for direct and diffuse irradiance. The broadband snow albedo, however, depends on the spectral distribution of the irradiance and the partitioning between its direct and diffuse components. Therefore, the broadband albedo depends on the spectral albedo and on the choice of the solar spectrum and atmospheric radiative transfer model used to estimate *S*_{λ} in equation .

We usually use the Simple Model for Atmospheric Radiative Transfer (SMARTS; Gueymard, , ) to model the incoming direct and diffuse solar spectral irradiances. To approximate conditions in the snow‐covered mountains of the western United States, among the myriad of model parameters, we used a midlatitude winter atmospheric profile, 3‐km surface elevation, coarse‐grained snow surface, and rural aerosol model (Shettle & Fenn, ). For comparison, we ran the same analysis with a subarctic summer atmospheric profile at sea level. As Figure shows, the atmospheric model and parameter choices substantially affect the broadband albedo; therefore, we chose parameters to best simulate the domain which we are studying.

For the snow reflectance, we compute the Mie () scattering parameters usually with the complex angular momentum approximation. Nussenzveig and Wiscombe () point out that the complex angular momentum approximation is valid for situations with a dimensionless Mie size parameter *β* = 2*πr*/*λ* > 10, where radius *r* and wavelength *λ* have the same units. The approximation is accurate enough and much faster than other Mie calculation methods, and it implicitly averages over size; hence, the complex angular momentum method is appropriate for all snow grain sizes. We resort to MatScat (Schäfer et al., ), a MATLAB adaptation of the Bohren and Huffman () FORTRAN code, only for calculations about dust or soot particles in the snowpack. Values of the complex refractive index of ice are from Warren and Brandt () with some visible and UV corrections from Picard et al. (). If the snow is dirty, Mie parameters for the dust or soot are averaged with those of the snow, with dust optical properties from Skiles et al. () and Flanner et al. (). The single scattering albedo and asymmetry factor estimated from the Mie equations are then fed into a two stream radiative transfer model, and the delta‐Eddington approximation (Joseph et al., ; Wiscombe & Warren, ) is used to solve for the spectral reflectance of snow. Once the spectral reflectance of the snow has been estimated, the broadband albedo is estimated using equation . Estimating the effect of other likely contaminants such as snow algae or forest litter requires a slightly different approach, because the optical properties of the absorbing materials are not calculated with Mie theory. Instead, the snow grain size can be recovered from reflectance in wavelengths beyond the visible spectrum (Dozier, ; Nolin & Dozier, ). Then reflectance values in the visible wavelengths indicate the degree to which forest litter (Melloh et al., ) or algae (Painter et al., ) degrade the albedo. Empirically, that degradation has been related to concentration of the contaminating materials.

Running the radiative transfer equations over an extended terrain, as would be required for large‐scale remote sensing work, is impractical. Thus, we have developed statistical fits for the curves shown in Figure . Following the methodology of Gardner and Sharp (), for a given solar zenith angle we fit the albedo versus grain size relationship for clean snow with a three‐parameter power law:

The coefficients *a*,*b*,and *d* depend on the solar zenith angle. The relationship follows a rational equation, whereby *μ*_{0} = cos*θ*_{0} and *a*,*b*,and *d* each follow a rational equation:

Solving for the values depends on the starting guesses, because a least squares or maximum likelihood estimate identifies a local minimum. Therefore, we run the fitting equations 100 times each with random starting values and select only those where the *R*^{2} value exceeds 0.999. Then if the maximum *R*^{2} value coincides with the minimum *RMSE*, we select those coefficients. Otherwise, or if *max*(*R*^{2}) < 0.999, we consider the combinations where *R*^{2} values are above the 95th percentile and *RMSE* values are less than the 5th percentile, we sort the *R*^{2} in descending order and the *RMSE* in ascending order, and we choose the combination with the minimum sum of ranks.

In matrix form, we can express this relationship as

The symbol ⊘ denotes Hadamard (element‐by‐element) division. Table shows example values for the ** P** and

In normal applications, we would apply these equations over a satellite image and coincident topography, which provide estimates of r and *μ*_{0} at every grid cell, expressed as matrices ** R** and

The symbols ⊗ and ⊘ denote Hadamard multiplication and division, and the exponents in *R*^{B} infer raising to the power element‐by‐element. If the computing language employed enables vector and matrix operations directly, without loops, equations – calculate the clean snow albedo at every grid cell. In MATLAB for example, computing the albedo over a grid needs only those four equations converted to four lines of code.

Figure shows a surface plot of the functional relationship *α* = f(*r*, cos*θ*_{0}) for the midlatitude winter atmosphere at 3‐km elevation. Compared to the modeled reflectance, the statistical fit for values in the range 0.07 ≤ *μ*_{0} ≤ 1, that is, *θ*_{0} < 86°, is almost exact, with values of *R*^{2} > 0.9999, *RMSE* = 2.11 × 10^{−4}, and a bias of 6.6 × 10^{−5}. For the Subarctic Summer, the corresponding values are *R*^{2} > 0.9999, *RMSE* = 1.96 × 10^{−4}, bias 1.29 × 10^{−4}. At solar positions so close to the horizon, the atmospheric models are less accurate. They must incorporate refraction, and the diffuse irradiance is much greater than the direct over much of the spectrum at those angles. Therefore, for solar zenith angles greater than 85°, we recommend using the *a*, *b*, and *d* values at *μ*_{0} = 0.09.

We are working on adding clouds and impurities to our statistical albedo model, using an approach similar to Gardner and Sharp (), but they are not required for this study for two reasons: (1) since MODIS is an optical instrument that cannot see through clouds, we can only validate albedo measurements on clear days and (2) we currently use a separate step for estimating the broadband albedo degradation due to LAP (section ). For the most exact albedo estimates for our SWE reconstructions, clouds must be accounted for; however under cloudy skies, longwave radiation dominates the energy balance, so we do not expect the inclusion of clouds in the albedo model to substantially affect the results.

For a slower but more flexible approach, in the Appendix we have included a lookup table function and associated data. As input, the function takes the parameters shown in Figure (solar zenith and snow radius), but adds local illumination angle, elevation, and concentration of dust or soot.

The approach above can used iteratively (Coleman & Li, ; Móre, ) to solve for the grain size and impurity concentration, given observed *α* and *α*_{nir}. These two broadband albedos are required because grain size affects albedo mostly in the near‐IR, while LAP affect the visible part of the spectrum. To interpret the albedos observed at the CUES site, the spectral range of the radiometers is used (provided in Bair, Davis, et al., ) along with *θ*_{0} in the SMARTS model to obtain the incoming spectral direct and diffuse solar irradiance. The size and concentration of the LAP together affect the spectral albedo, with smaller particles being more effective. However, one can almost perfectly mimic the spectra with different pairs of values. Therefore, we simply assume a typical size value for dust (3 μm) and solve for concentration. While this assumption does not yield the actual LAP concentration, it does provide the right spectral albedo. Once a solution is reached, the model is currently set up to produce an estimate of the LAP concentration by mass fraction rather than ΔVIS. To obtain ΔVIS, estimates of the clean and dirty snow albedos over the visible range are made and applied to equation .

Sometimes however, the solved‐for grain sizes are unrealistically large, and several plausible explanations exist. For example, the snowpack could be wet or refrozen on the surface. Wet snow forms clusters of grains (Colbeck, ), which behave optically as large grains. In laboratory measurements of reflectance, O'Brien and Munis () observed that near‐IR reflectance decreased as the snow melted, but the values did not rebound when the snow refroze. These interpretations of large grain sizes could also be caused by dark objects, like trees, in the radiometers' field of view. Further, fitting a plane is a simplification of the snow surface geometry, but necessary given that our radiative transfer model is not yet adapted to rough surfaces. Several studies show that macroscopic roughness of the snowpack such as sun cups, sastrugi, and nieve penitentes reduce the snow surface albedo (Lhermitte et al., ; Warren et al., ; Zhuravleva & Kokhanovsky, ). Warren et al. () suggest that these rough features reduce snow albedo for two reasons: (1) decrease of effective zenith angle and (2) photon trapping. Especially at higher zenith angles, these rough features reduce the effective zenith angle as they can cause more of the snow to be oriented toward the sun. Photon trapping is an increase in the probability of absorption as the number of scattering events increases in a rough versus smooth surface.

Because of different plausible causes of reduction in observed snow albedo, we set the *r* upper limit at 1,500 μm, a somewhat arbitrary value well above what has been measured using different techniques (Aoki et al., ; Matzl & Schneebeli, ; Nolin & Dozier, ; Painter et al., ; Skiles & Painter, ). Values above this limit were discarded. Our reasoning behind this high limit is that the MODSCAG/MODDRFS‐derived albedos are likely to be impacted by some of the same artifacts as the in situ observations, which could cause unrealistically large grain size estimates. Thus, by allowing these large grain sizes, we allow for some of these albedo darkening mechanisms to be accounted for, but suggest that *r* values above 1,500 μm are not reliable. Additionally, there were 44 measurements (February through December 2007) at SASP and 45 at SBSP that had zeros values for the *f*_{SCA} and *r* due to a processing error at the NASA JPL server.

Because our hypothesis is that our remotely sensed albedos are more accurate than age‐based albedo decay models, we applied the Biosphere‐Atmosphere Transfer Scheme (BATS; Dickinson et al., ), which was the most accurate of the age‐based models in a previous study (Bair et al., ) with an RMSE of 9% and a bias of −6%. As with all age‐based decay models, BATS requires a precise estimate of when the last snowfall ended. Such a value is not available in most parts of the world, as in situ measurements are needed. Manual measurements are preferable, for example to distinguish the end of a storm from snow redistributed by wind, which cannot be discerned using automated snow depth measurements. As the Mammoth Mountain Ski Patrol takes manual snow depth and weight measurements each day that it snows, a reasonably good estimate of when it stopped snowing prior to the albedo measurement can be made at CUES. Still, these manual weather measurements are only consistently made once a day. Thus, even with these careful daily measurements, there is still uncertainty in determining precisely when it stopped snowing. As in our previous study, we used a 10‐mm threshold for new snow to reset to the maximum broadband albedo value of 0.86. BATS also requires the local solar zenith angle *θ* and the snow surface temperature as inputs. The local solar zenith was computed from the lidar snow surface scans (section ). Snow surface temperature is now measured at CUES but was not available over most of the study period, so we used the lesser of the measured air temperature or 273.15 K as in our previous study (Bair et al., ).

Figure shows remotely sensed versus in situ modeled grain radius and ΔVIS at CUES (Figures a and d), SASP (Figures b and e), and SBSP (Figures c and f) for each of the filtering, interpolation, and smoothing steps described in section .

Table shows error statistics are for each of the three sites and mean values. Note that these errors represent the best match to the in situ measurements from a 9‐pixel neighborhood centered on each of the sites. This best‐of‐9‐pixel neighborhood approach is often used with MODIS measurements (Bair, Calfa, et al., ; Rittger et al., ) because of geolocational uncertainty (section ) and spatial variability of the snow surface.

Figures d–f clearly show that CUES is less dusty than SASP or SBSP with no ΔVIS values larger than 25% and many in situ observations indicating a value of zero. Also, the relative errors of ΔVIS are larger than those of grain size across all three sites. The *interpolated* product performs better than the *filtered* or *smoothed* product for every site and for both variables, except for one case (CUES ΔVIS). Overall, there is a slightly negative bias in grain size and ΔVIS, but the magnitude is small, averaging −34 μm in *r* and around −1% for ΔVIS.

A useful way to put these errors in perspective is to examine their impact on *α*. Using a 500‐μm grain radius *r* and cos*θ* = 2/3, the mean RMSE values of 118 to 179 μm across all three sites in grain sizes are comparable to a 1.2 to 1.7% difference in *α*. The impact of the 3.6 to 6.2% mean RMSE in ΔVIS can be found by multiplying ΔVIS by *f*_{VIS}, the fraction of total irradiance in the visible spectrum. Based on our results from SMARTS, we find *f*_{VIS} ≈ 0.63, relatively insensitive to illumination conditions for clear skies**.** Note that this estimate is slightly greater than our previous (e.g., Bair et al., ) and less accurate estimate of *f*_{VIS}, which was 1/2. Thus, the impact on the broadband albedo of the RMSE error in ΔVIS is 2.3% to 3.9%. Together the combined impact of the grain size and ΔVIS RMSE on broadband albedo is 3.5% to 5.6%, which is lower than the 7% to 17% RMSE found in age‐based decay models. The mean bias in grain size (−23 to −51 μm) and −2.4 to 1.2% are equivalent to about a ±1% bias at the most.

For comparison, we use *r* and ΔVIS from the *interpolated* images which we found most accurate in section to produce *α* estimates for comparison with the age‐based BATS model. Figures and a–6g and Table show that the RMSE for BATS is 2 times greater (4.8% versus 9.6%) and the bias is more than 4 times greater in magnitude (1.3% versus −5.4%) than for our remotely sensed values. The remotely sensed albedos show increasing error at the lower values of *α* at the end of the melt season (Figures a–g), which correspond to larger grain size values and/or increased contamination. The BATS *α* errors did not show a trend. We find this to be clear evidence that even with the best available in situ measurements, the best age‐based albedo model was considerably less accurate and more biased than our remotely sensed albedo estimates.

To examine the sensitivity of reconstructed SWE to albedo, we ran our Parallel Energy Balance model (ParBal; Bair et al., ) to reconstruct SWE at CUES and SASP using measured forcings from the suite of radiometers and sensors for wind, temperature, and humidity (Bair, Davis, et al., ; Skiles, ). We used the in situ measured albedo and melt out date, but also introduced two sources of error: random error and biased error. For our error values, we used error results from the previous section, that is, Table .

We did not use our remotely sensed albedo estimates or estimates from BATS directly. BATS could not be applied at SASP because there was no record of when the last snowfall ended. And applying our remotely sensed albedos would have required interpolation over cloudy days with no way to account for geolocational uncertainty and spatial variability of the snow surface across the MODIS pixel. Rather, we perturbed the in situ measured albedo by the error in Table for a tractable comparison which could be used to examine the effect of any RMSE and bias values in albedo on reconstructed SWE. For the random error, we assumed normally distributed error with a mean of zero and a standard deviation of the RMSE shown in Table . For the bias error, we added the bias error shown in Table .

For comparison, we include the measured SWE at both sites. For CUES, we used estimated SWE values directly in front of the platform without shade. We did not simply use the values measured from the snow pillow because the pillow is located under shade from trees and gets additional wind‐deposited snow. Thus, the pillow stores more SWE and melts out later than the open area by the radiometers. At first, to estimate the SWE closer to the radiometers, we tried to directly estimate the density of the snow above the pillow using the acoustic depth sensor and the pillow SWE, but the acoustically measured depths were too noisy most of the time for reliable estimates. Instead we used a linear regression between the acoustic depth sensor at the platform and the sensor directly above the pillow. Then, we adjusted the SWE using the slope coefficient, but converted the intercept from centimeter of depth to centimeter of SWE using the median snow density measured above the pillow. Note the switch from RMSE to MAE for SWE, as we use RMSE to examine errors for inputs into the ParBal, but MAE to examine errors in the reconstructed SWE (Bair, Davis, et al., ). The reasoning is that we apply greater penalties for greater errors in the ParBal inputs using RMSE, but not in the SWE output with MAE. Also note that CUES lacked a snow pillow from the late 1990s until WY 2013, so we start the reconstructions then.

At SASP and SBSP, there are no snow pillows, so SWE values were computed using manual snowpits. This method yields accurate values of total SWE, but measurements are limited to times when pits can be dug. Because of the difficulty in reaching SBSP, that site has fewer pits than SASP; thus, we only used the SASP site for the reconstructed SWE comparison. Note that both CUES and SASP have well‐measured peak SWE values and snow disappearance dates, the two most important validation points for the reconstruction method.

Figures a–e show a time series of the measured and reconstructed SWE values for CUES, and Figure shows the same for SASP. Table summarizes the errors for both sites. The salient observations are (1) the reconstructed SWE values without added noise and bias (i.e., measured albedo) agree well with the measured SWE at each site; (2) adding the remotely sensed noise and bias to the measured albedo actually improves the accuracy of the reconstructed SWE; (3) adding the BATS noise and bias led to large overestimates of peak SWE; and (4) although not shown here, the added bias had a much larger effect on the reconstructed SWE than the added noise. We note that finding (2) suggests that our additive noise/bias serves as a proxy for some other process in our reconstructed SWE model that is reducing SWE. Figures a and b show this with the red curves (remotely sensed noise and bias) being slightly lower than the blue curves (no noise, no bias). We might expect this because the SWE reconstruction technique will overestimate peak SWE when there is accumulation during the ablation season. The reason is that the melt, which is summed in reverse to estimate SWE, is constrained to be ≥0. Thus, while going backward from melt out to the peak, the method cannot decrease SWE for any given point in time, which is effectively what happens just before a significant snow storm during the melt season. Hence, reconstruction provides an upper envelope for SWE values when snow accumulation occurs during the melt season.

Finding (3) illustrates the sensitivity of the reconstructed SWE to errors in albedo; albedo noise with a 9.6% standard deviation and a −5.4% bias led to a 12–15% increase in peak SWE bias (5% to 20% for CUES and 6% to 18% for SASP) compared to the no‐noise, no‐bias reconstruction case.

We improved a broadband albedo model by accounting for atmospheric conditions and incoming solar spectra. Using results from this model shows the importance of choosing the correct atmospheric conditions. This is especially important when using a statistical fit for the broadband albedo. Additionally, different atmospheric models (SBDART and SMARTS) produce substantially different results, even when given identical parameters.

We then applied this improved albedo model to in situ albedo measurements at three high‐altitude sites in the western United States to validate remotely sensed retrievals from MODIS. To our knowledge, this is the first comprehensive validation of remotely sensed albedos over such a long time period, 1,632 total days. We would have liked to have used more sites in a broader range of snow climates, but we are not aware of any other sites in the world where broadband snow albedo, near‐infrared snow albedo, and local terrain under the radiometers have been measured for a comparable period of time.

We tested grain size and LAP content retrievals separately. This study also differs from previous studies in that different levels of filtering, interpolation, and smoothing were validated. Examination of these techniques is relevant as they must be used to create continuous and accurate time series of MODIS snow cover retrievals.

Our remotely sensed broadband albedo estimates compare favorably with the in situ estimates, with RMSE values of 4–6% and negligible bias. Errors in the LAP content were relatively higher than those in grain size. For comparison, the most accurate age‐based albedo, based on a previous study, showed a 10% RMSE and −5% bias.

We examined the impact of these errors on reconstructed SWE at two of the sites. Using the measured albedo as a baseline, we added noise and bias to simulate the errors of the remotely sensed albedos and age‐based albedos. As expected the bias error in albedo substantially impacted the reconstructed SWE. The baseline showed a 7–13% MAE when compared with measured peak SWE. The accuracy of the reconstruction actually improved when simulated remotely sensed noise and bias were added, which had the effect of raising the albedo. We suggest that this result highlights a shortcoming with reconstructed SWE, that it overestimates peak SWE when there are snow accumulation events during the ablation season. In comparison, the age‐based decay model showed large MAE and bias (18–20%) when used to reconstruct the measured peak SWE.

We find that our remotely sensed albedos are superior to age‐based models in all aspects except simplicity. Perhaps even more advantageous than the demonstrated improvements in accuracy, is the ability for remotely sensed albedos to be used globally, in places with austere surface infrastructure such as much of High Mountain Asia, as long as obscuration by clouds is not too frequent. We will focus future improvements in areas that challenge optical sensors, such as under clouds and in forests.

At the suggestion of A. W. Nolin, we have implemented a lookup table, in lieu of equations –, to calculate broadband albedo as a function of grain size (optically equivalent radius), solar zenith angle, local solar illumination angle (to account for slopes) elevation, and concentration of a light‐absorbing particle, in this case dust of size *r* = 3 μm from the Colorado Plateau, whose optical properties were estimated by Skiles et al. (), or soot of size *r* = 0.5 μm, with optical properties from (Bond & Bergstrom, ). We implement the lookup function in five dimensions using the MATLAB griddedInterpolant function with Modified Akima cubic Hermite interpolation (Akima, ), introduced in MATLAB R2017b (Moler, ). The MATLAB function AlbedoLookup.m and the lookup table albedoLUT.mat are available at 10.5281/zenodo.3228428. Computing albedo values with equations – is about 15 times faster than with the lookup table, but the lookup table provides more flexibility. It handles dust or soot concentrations, accounts for variation in altitude that affects the spectral distribution of the incoming solar radiation, and includes both solar zenith and local illumination angles and therefore accounts for variability that results from path length through the atmosphere and forward scattering at the snow surface. The posted version of the lookup table uses the SMARTS 2.9.5 midlatitude winter atmosphere, but the tables could be created to compute values for other atmospheric profiles.

The in situ measurements at CUES have been published previously (Bair, Davis, et al., ). The in situ measurements for SASP and SBSP are based off published data from an earlier study (Landry et al., ), but the updated data sets used in this study are now available at

This work was supported by NASA awards 80NSSC18K1489, 80NSSC18K0427, and NNX10AO97G; NASA JPL award 1591172; NOAA award NA18OAR4590380; and University of California award LFR‐18‐54831. We thank the Center for Snow and Avalanche Studies for the maintenance of Senator Beck Basin Study Area and their data availability. We appreciate reviews from Joseph Cook and Anne Nolin.