Divergence of apparent and intrinsic snow albedo over a season at a sub-alpine site with implications for remote sensing

Intrinsic albedo is the bihemispherical reflectance independent of effects of topography or surface roughness. Conversely, the apparent albedo is the reflected radiation divided by the incident and may be affected by topography or roughness. For snow, the surface is often rough, and these two optical quantities have different uses: intrinsic albedo is used in scattering equations whereas apparent albedo should be used in energy balance models. Complementing numerous studies devoted to surface roughness and its effect on snow reflectance, this work analyzes a time series of intrinsic and apparent snow albedos over a season at a sub-alpine site using an automated terrestrial laser scanner to map the snow surface topography. An updated albedo model accounts for shade, and in situ albedo measurements from a field spectrometer are compared to those from a spaceborne multispectral sensor. A spectral unmixing approach using a shade endmember (to address the common problem of unknown surface topography) produces grain size and impurity solutions; the modeled shade fraction is compared to the intrinsic and apparent albedo difference. As expected and consistent with other studies, the results show that intrinsic albedo is consistently greater than apparent albedo. Both albedos decrease rapidly as ablation hollows form during melt, combining effects of impurities on the surface and increasing roughness. Intrinsic broadband albedos average 0.056 greater than apparent albedos, with the difference being 0.052 in the near infrared or 0.022 if the average (planar) topography is known and corrected. Field measurements of spectral surface reflectance confirm that multispectral sensors see the apparent albedo but lack the spectral resolution to distinguish between darkening from ablation hollows versus low concentrations of impurities. In contrast, measurements from the field spectrometer have sufficient resolution to discern darkening from the two sources. Based on these results, conclusions are as follows: (1) impurity estimates from multispectral sensors are only reliable for relatively dirty snow with high snow fraction; (2) a shade endmember must be used in spectral mixture models, even for in situ spectroscopic measurements; and (3) snow albedo models should produce apparent albedos by accounting for the shade fraction. The conclusion re-iterates that albedo is the most practical snow reflectance quantity for remote sensing.

1 Introduction

Snow albedo plays an important role in Earth's climate and hydrology. For example, a small (0.015 to 0.030) decrease in snow albedo over the Northern Hemisphere is twice as effective as a doubling of CO₂ at raising global air temperature (Hansen and Nazarenko, 2004). Likewise, during the COVID-19 lockdowns, a cleaner snowpack, presumably from a reduction in anthropogenic emissions, prevented 6.6 km³ of snow/ice from melting in the Indus River basin (Bair et al., 2021a), more water than is stored in the largest reservoir in California. Yet, snow albedo is difficult to measure (Bair et al., 2018), especially in the mountains where lighting conditions vary dramatically. To understand Earth's climate and the effect humans have on it, an understanding of how snow surface topography affects snow albedo is imperative. The concepts of intrinsic and apparent albedos form the basis of this study. Intrinsic albedo is the bihemispherical reflectance (Nicodemus et al., 1977; Schaepman-Strub et al., 2006) of a substance independent of effects of roughness or topography. Apparent albedo is the ratio of the reflected divided by the incident radiation and may incorporate artifacts caused by roughness or topography. Here we use the term albedo to refer to a broadband albedo, covering the solar spectrum. Albedos covering a narrower spectral range are denoted with additional descriptors such as near-infrared albedo. Since the snow surface is rarely smooth, distinction between apparent and intrinsic albedo is an important consideration that is often ignored. For example, MODIS measurements of snow albedo that comprise the National Solar Radiation Database have been found to be positively biased because they fail to account for surface roughness (Gueymard et al., 2019). Both albedos should be studied, as apparent and intrinsic albedos have different uses. An apparent albedo should be used when modeling energy budgets (Bair et al., 2016), as it dictates how much shortwave radiation is absorbed by the surface. Intrinsic albedos are needed to understand changes in snow properties that affect albedo, such as changes in grain size and darkening from light-absorbing particles like soot or dust (Clarke and Noone, 1985; Jones, 1913; Warren, 2019).

Most snow albedo models follow approaches developed 4 decades ago, based on radiative transfer (Warren, 1982). These models provide intrinsic albedos controlled by illumination angle, water equivalent when snow is shallow, and grain-scale snow properties, which have included grain size (Wiscombe and Warren, 1980), grain shape (Libois et al., 2013), snow structure (Kaempfer et al., 2007), direct and indirect effects of light-absorbing particles (Picard et al., 2020; Skiles and Painter, 2019), and vertical heterogeneity (Zhou et al., 2003). Other efforts have focused on rapid calculation (Bair et al., 2019; Flanner et al., 2021; Gardner and Sharp, 2010) and inversion from remotely sensed imagery (Bair et al., 2021b; Nolin, 2010; Painter et al., 2012a). Weiser et al. (2016) present a correction for albedometers over snow where the underlying terrain is unknown, based on modeled or measured irradiance from nearby well-leveled radiometers, but not accounting for surface roughness. A shade endmember has been introduced to account for lighting differences across surfaces (Adams et al., 1986), thereby enabling the use of an apparent albedo for quantitative spectroscopy. These shade endmembers have proven successful when applied to snow cover mapping (Bair et al., 2021b; Nolin et al., 1993; Painter et al., 2003; Rosenthal and Dozier, 1996). Yet, the widely used albedo models cited above do not account for varying illumination within the field of view, meaning their results can be positively biased.

Features that affect snow roughness include suncups (ablation hollows), penitentes, and wind-formed features like ripples, sastrugi, and dunes (Filhol and Sturm, 2015). Because of their topographic variation in solar exposure, all of these roughness features can significantly affect apparent albedo. Matthes (1934) described “suncups” as having “a honeycombed appearance, the surface being pitted with deep cell-like hollows”. However, Rhodes et al. (1987) use the term “ablation hollows” to describe these features as they are not always caused by solar radiation. Instead Rhodes et al. (1987) find that the presence of impurities on the snow surface governs the formation of ablation hollows, growing in direct sunlight for relatively clean snow and decaying in dirty snow (Lliboutry, 1964). This hypothesis was confirmed with a field experiment where an ash-covered snowfield on Mount Olympus from the Mount Saint Helens eruption was cleared. After 2 weeks, the ash-free area had developed larger ablation hollows than the rest of the ash-covered snowfield (Rhodes et al., 1987). Observations of penitentes go back to Darwin (1845, Ch. XV). Penitentes are columns of snow that point at the sun and are thought to be sublimation features (Betterton, 2001). Penitentes can be much larger than ablation hollows, with measured heights over 2 m (Lhermitte et al., 2014). Ripples, sastrugi, and dunes are formed by wind erosion whose orientation varies with the direction of the prevailing winds (Filhol and Sturm, 2015; Seligman, 1936). Warren et al. (1998) report that sastrugi can reduce albedo by altering the angle of incidence for direct solar radiation and by trapping photons through multiple reflections.

Several studies have attempted to model the reflectance of roughness features with simple shapes (Carroll, 1982; Leroux and Fily, 1998; Zhuravleva and Kokhanovsky, 2011), with more recent studies employing ray tracing of three-dimensional surface models (Larue et al., 2020; Manninen et al., 2021). A few studies have focused on the surface roughness and the implications for remote sensing by incorporating multiple viewing geometries (Corbett and Su, 2015; Kuchiki et al., 2011; Lyapustin et al., 2010; Nolin and Payne, 2007) or by measuring spatial variability within a satellite sensor pixel (Wright et al., 2014). These approaches are well-suited toward expansive high-latitude snowpacks but ill-suited towards dynamic midlatitude snowpacks with mixed pixels where the snow cover can change between satellite overpasses. The consensus in the literature is that roughness features can lower the snow albedo by up to 0.40, but decreases of a few percent are more common. To our knowledge, none of these studies have tracked the snow surface topography throughout a snow season, nor have they examined the effects of snow surface topography on spectral mixture analysis.

2 Approach

2.1 Radiometric measurements

Albedos were measured (Fig. 1) at CUES – Cold Regions Research and Engineering Laboratory and University of California, Santa Barbara Energy Site – on Mammoth Mountain, CA, USA (Bair et al., 2015). To eliminate darkening from the ground, shadowing from vegetation, and effects from high zenith angles, only clear days with a deep, optically thick snowpack were examined. Radiometer measurements were taken at the satellite overpass time (Sect. 2.3). Uplooking and downlooking Eppley precision spectral pyranometers (PSPs) with both clear (285–2800 nm) and near-infrared (700–2800 nm) domes were located on both the fixed and adjustable arms, providing redundant measurements of the incoming irradiance in both wavelength regions, and providing measurements of reflected radiation from both the fixed and adjustable arms. The adjustable arm keeps its downlooking radiometers about 1 m above the snow surface, whereas the fixed arm is mounted 8 m above the ground, so its distance from the snow surface depends on the snow depth. In measuring the reflected radiation, two artifacts must be minimized. If the downlooking radiometer is too far above the snow, the field of view is too large, so other, darker objects like the tower itself and trees, will cause the snow albedo to be too low. Conversely, if the radiometer and its arm are too close to the snow, they will cast a shadow that will also cause the albedo to be too low. By experiment, we found that the combination of these two artifacts is minimized when the radiometer is ∼ 1 m above the snow, so as the snow depth changes, we maintain the adjustable arm's height at about that distance.

Figure 1 Fixed and adjustable albedo arms at the CRREL UCSB Energy Site (CUES) in the summer (a) and winter (b).

The ratio of diffuse to direct irradiance was computed using a Delta-T SPN1 Sunshine pyranometer mounted on the fixed arm, which integrates over a slightly different spectral band (400–2700 nm) than the PSP clear. Because of the different response and biases (Habte et al., 2015; Wilcox and Myers, 2008) arising from issues such as thermal offsets (Haeffelin et al., 2001), only the diffuse ratio (used in the terrain correction described in Sect. 2.2) from the SPN1 was used. The irradiance measured by each PSP was split into direct and diffuse components using this ratio. Calculations using SMARTS v2.9.8 (Gueymard, 2019) provide an estimate of the spectral distribution of irradiance not subject to instrument error. We use the SMARTS simulations to adjust the measurements of the diffuse fraction from the SPN1 (400–2700 nm) to account for the diffuse fraction in the irradiance measurements from the PSPs with clear and near-infrared domes. The accuracy of an atmospheric radiation model depends on the accuracy of the estimates of the atmospheric properties, principally aerosols and water vapor. Errors in field radiometer measurements stem from calibration inaccuracies and siting of the instrument. The comparison between SMARTS and the measurements yields R²≥0.99 for both the PSP and the SPN1 (Fig. 2 and Table 1), suggesting sufficient relative accuracy to make both instruments suitable for albedo measurement. However, the reflected radiation is measured by downlooking PSPs – there is no downlooking SPN1 – so we used the same type of radiometers (PSPs) to measure the irradiance and reflected solar radiation.

Figure 2 Measured vs. modeled irradiance at CUES for three broadband sensors: an Eppley precision spectral pyranometer (PSP) mounted on an adjustable albedometer arm kept ∼ 1 m above the snow surface (PSP clear, adjustable), a PSP mounted ∼ 8 m above bare ground (PSP clear, fixed), and a Delta-K SPN1 Sunshine pyranometer also mounted ∼ 8 m above bare ground (SPN1 global, fixed). The differences between the instruments, particularly at high radiation values, are likely caused by different thermal responses (Haeffelin et al., 2001).

Table 1 Radiometer measurement differences shown in Fig. 2.

Reflected radiation was measured using the downlooking PSPs, in both broadband (285–2800 nm) and near-infrared (700–2800 nm) wavelengths. We mounted one pair of PSPs on the adjustable computer-controlled and self-leveling arm, kept ∼ 1 m above the snow surface to prevent non-snow objects from being seen, and the other pair on the fixed arm 8 m above the bare ground. To illustrate the effect of non-snow objects within the downlooking radiometers' fields of view, Fig. 3 shows a comparison.

Figure 3 Reflected radiation, from downlooking radiometers, and snow depth measured at CUES.

When the snowpack is deep and continuous spatially, the downlooking radiometers on the adjustable boom have greater values than those on the fixed arm (Fig. 3, 10–17 May 2021). This condition occurs because darker non-snow objects are within the radiometers' fields of view on the fixed arm. Contrast this to the snow-free condition at the end of May where reflected radiation is the same for the radiometers on both the fixed arm and adjustable arm. In patchy snow, the opposite occurs; on 19–20 May 2021, the radiation measured by the nIR PSP on the fixed arm exceeds that of the clear PSP. This condition occurs because the radiometers on the fixed boom view additional emerging vegetation with a higher nIR albedo than snow. Thus, to prevent non-snow objects from contaminating the snow albedo measurements, only the downlooking radiometers on the adjustable arm were used to measure reflected radiation.

Although a radiometer views a hemisphere, the downlooking field of view is restricted to about ∼ 150^∘ due to manufacturing constraints (Wu et al., 2018). Sailor et al. (2006) showed that the size of a radiometer's field of view that accounts for 95 % of the reflected radiation is ∼8.7h, where h is the height of the radiometer above the surface. The radiometer's height above the snow surface of h∼1.0 m translates to a footprint diameter d of 8.7 m. In comparison, the downlooking radiometers on the fixed arm 8 m above bare ground would see a footprint larger than 40 m over snow with 1 m depth. To our knowledge, CUES is the only site where snow albedo is measured using an adjustable albedo arm. Given such a large footprint, an examination of published images of tower arms at other sites where snow albedo is measured (Elder et al., 2009; Landry et al., 2014; Lejeune et al., 2019; Lhermitte et al., 2014) shows non-snow objects within the downlooking radiometer's field of view at every site.

2.2 Surface topography and corrections

A Riegl VZ-400 laser scanner automatically scanned the snow surface every hour during the 2021 water year. Point clouds were converted to surfaces as follows. Noise was removed using a filter (Rusu et al., 2008), and additional days with blowing snow were manually removed because the moving particles obscure the snow surface (Bair et al., 2012). The adjustable albedometer arm was removed from the point clouds using a morphological filter (Pingel et al., 2013). Point clouds were converted to surfaces with 1 cm spatial resolution using bilinear interpolation. A radial mask was applied to the surface to simulate the footprint seen by the downlooking PSP. Slope and aspect were computed for a plane fit to the surface. The rough surface combines with the local illumination angle to affect the apparent snow albedo.

Four broadband albedos were computed. An uncorrected apparent albedo is computed as

$$\begin{array}{}\text{(1)}& {\mathit{\alpha }}_{\mathrm{uncorrected}}={\displaystyle \frac{D_{\uparrow }}{I_{\downarrow }}},\end{array}$$

where D_↑ is the reflected radiation measured by the downlooking PSP, and I_↓ is the irradiance measured by the uplooking PSP. An albedo with a plane fit to the surface built from the point cloud is computed as

$$\begin{array}{}\text{(2)}& {\mathit{\alpha }}_{\mathrm{planar}}={\displaystyle \frac{D_{\uparrow }}{c{B}_{\downarrow }+D_{\downarrow }}},\end{array}$$

where $c=\cos {\mathit{\theta }}_{\mathrm{S}}/\cos {\mathit{\theta }}_{\mathrm{0}}$ is a correction factor of a sloped to a level surface. θ_S is the illumination angle for the plane, θ₀ is the solar zenith angle for a level surface, B_↓ is the direct irradiance, and D_↓ is the diffuse irradiance. This planar correction has been applied in previous work (Bair et al., 2018; Painter et al., 2012b). Because the ratio c is in the denominator of Eq. (2), α_planar>α_uncorrected when cos θ_S<cos θ₀ and equal when the angles are equal, less otherwise.

An albedo with a spatial correction to account for the rough surface is computed by considering the effects for a generic point on the rough surface and then averaging those effects over the downlooking radiometer's field of view, i.e., a circle with 8.7 m diameter. Every point on the surface has slope S and aspect A, and ϕ₀ is the solar azimuth. The cosine of the illumination angle at each point is

$$\begin{array}{}\text{(3)}& \cos {\mathit{\theta }}_{\mathrm{S}}=\max \left[\mathrm{0},\cos {\mathit{\theta }}_{\mathrm{0}}\cos S+\sin {\mathit{\theta }}_{\mathrm{0}}\sin S\cos \left({\mathit{\varphi }}_{\mathrm{0}}-A\right)\right].\end{array}$$

The use of the max function sets the value of cos θ_S to zero on self-shaded slopes, when otherwise the cosine would be negative. In addition to the slope affecting the magnitude of the irradiance, local horizons formed by neighboring points, in or out of the same ablation hollow, affect the illumination in two ways: (1) a neighboring high point might shade a slope that would otherwise be illuminated, and (2) the set of horizons in all directions partly obstructs the overlying hemisphere. We define the view factor V_Ω as the fraction of the hemisphere that is open to the sky; a completely unobstructed surface has a view factor V_Ω=1. Dozier (2022) describes methods to rapidly compute the horizons and the view factor.

Considering the albedo of a rough snow surface involves multiple reflections. Over a range of wavelengths, the spectral distribution changes with each reflection. Therefore, the initial approach to model this effect uses monochromatic radiation, with ρ to indicate a spectral albedo, omitting a wavelength identifier unless necessary. Setting F_dif as the fraction of the spectral irradiance that is diffuse and setting the value of the initial irradiance on a horizontal surface to I, the “spatial” spectral radiation that initially escapes into the overlying hemisphere without being re-reflected is

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& I_{\mathrm{esc}}^{\left(\mathrm{0}\right)}=\\ & {\mathrm{IV}}_{\mathrm{\Omega }}\left[{\displaystyle \frac{\cos {\mathit{\theta }}_{\mathrm{S}}}{\cos {\mathit{\theta }}_{\mathrm{0}}}}\left(\mathrm{1}-F_{\mathrm{dif}}\right){\mathit{\rho }}_{\mathrm{intrinsic}}^{\left(\text{direct}\right)}+F_{\mathrm{dif}}{\mathit{\rho }}_{\mathrm{intrinsic}}^{\left(\text{diffuse}\right)}+\left(\mathrm{1}-V_{\mathrm{\Omega }}\right){\left({\mathit{\rho }}_{\mathrm{intrinsic}}^{\left(\text{diffuse}\right)}\right)}^{\mathrm{2}}\right]\\ & \mathrm{\text{directly reflected}}\mathrm{\text{diffusely reflected}}\end{array},\end{array}$$

with ρ_intrinsic as the intrinsic spectral albedo on a level, smooth surface unaffected by topography; the superscripts designate the albedo to direct vs. diffuse irradiance. The right-hand term inside the brackets accounts for reflected radiation within a point's field of view impinging on the point. The direct and diffuse spectral albedos of snow differ slightly (Wiscombe and Warren, 1980); the major difference in the broadband values lies in the different spectral distributions of the direct and diffuse irradiance. Generally, ${\mathit{\alpha }}_{\mathrm{intrinsic}}^{\text{(diffuse)}}$ will be larger because the diffuse irradiance more heavily concentrates in the wavelengths where snow is brightest.

Not all the initially reflected radiation escapes into the overlying hemisphere. Instead, some of it re-reflects and eventually escapes or is trapped (Warren et al., 1998) by the roughness. The re-reflected radiation that does not escape is subject to possible internal reflection, its initial value being

$$\begin{array}{}\text{(5)}& I_{\mathrm{internal}}^{\left(\mathrm{0}\right)}=I_{\mathrm{esc}}^{\left(\mathrm{0}\right)}\left({\displaystyle \frac{\mathrm{1}-V_{\mathrm{\Omega }}}{V_{\mathrm{\Omega }}}}\right).\end{array}$$

To account for multiple reflections, at each reflection the value of the incident radiation is multiplied by the fraction (1−V_Ω) that accounts for the reflection remaining in the ablation hollow, the fraction V_Ω that escapes, and the spectral albedo. The albedo of the re-reflected radiation, ${\mathit{\alpha }}_{\mathrm{intrinsic}}^{\left(\text{RR}\right)}$, is biased toward the wavelengths where snow is brightest. An order-of-scattering approach to the multiple reflections lets some reflected radiation escape at each iteration n, and some remains available for re-reflection:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathrm{escaped}{I}_{\mathrm{esc}}^{\left(n\right)}& =I_{\mathrm{internal}}^{\left(n-\mathrm{1}\right)}{\mathit{\rho }}_{\mathrm{intrinsic}}^{\left(\text{diffuse}\right)}{V}_{\mathrm{\Omega }}\\ \mathrm{remaining}{I}_{\mathrm{internal}}^{\left(n\right)}& =I_{\mathrm{internal}}^{\left(n-\mathrm{1}\right)}{\mathit{\rho }}_{\mathrm{intrinsic}}^{\left(\text{diffuse}\right)}\left(\mathrm{1}-V_{\mathrm{\Omega }}\right).\end{array}\end{array}$$

This series converges in a half dozen iterations because $I_{\mathrm{internal}}^{\left(n\right)}$ declines in proportion to ${\left(\mathrm{1}-V_{\mathrm{\Omega }}\right)}^n$. The spatial spectral albedo ${\mathit{\rho }}_{\mathrm{spatial}}=\sum I_{\mathrm{esc}}/I$.

To adapt Eqs. (4) through (6) to compare modeled and measured albedo integrated over a range of wavelengths – for example the broadband and near-infrared albedos described in Sect. 2.1 – ρ_intrinsic cannot simply be replaced with α_intrinsic, because wavelength-integrated albedo depends on the convolution of the spectral albedo with spectral distribution of the irradiance. Including the spectral identifier λ, the wavelength integrated albedo is

$$\begin{array}{}\text{(7)}& \mathit{\alpha }={\displaystyle \frac{\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}\mathit{\rho }\left(\mathit{\lambda }\right)I\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}I\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},\end{array}$$

where ρ(λ) varies with wavelength, so α×I has a different spectral distribution than I itself. That distribution is weighted toward the wavelengths where ρ(λ) is larger, so each reflection causes α to increase even though ρ(λ) does not change. To address this problem, we derive an empirical function to estimate intrinsic broadband and near-infrared albedos at step n. In Eqs. (4) through (6), ρ_intrinsic is replaced with ${\mathit{\alpha }}_{\mathrm{intrinsic}}^{\left(n\right)}=f\left({\mathit{\alpha }}_{\mathrm{intrinsic}}^{\left(\mathrm{0}\right)},n\right)$. α_spatial is modeled at every point in each day's topographic grid. For each day, the mean of those values, $\stackrel{\mathrm{‾}}{{\mathit{\alpha }}_{\mathrm{spatial}}}$, over the field of view of the downlooking radiometer is equivalent to the measured α_uncorrected, so comparing the model to the measurement enables solving for the intrinsic wavelength-integrated snow albedo ${\mathit{\alpha }}_{\mathrm{intrinsic}}^{\left(\mathrm{0}\right)}={\mathit{\alpha }}_{\mathrm{intrinsic}}$.

To create $f\left({\mathit{\alpha }}_{\mathrm{intrinsic}}^{\left(\mathrm{0}\right)},n\right)$, we generated solar irradiance spectra using SMARTS (Gueymard, 2019) over observed solar zenith angles, 23 to 63^∘. We modeled spectral snow albedo (Warren, 1982) over the range of zenith angles, snow grain effective radii from 50 to 1000 µm, and mass concentrations of dust from 10⁻⁸ to 10⁻³ (i.e., 10 ng g⁻¹ to 1 g kg⁻¹), assuming an effective dust radius of 3 µm, comparable to measured values, and dust optical properties from measurements by Skiles et al. (2017) from the San Juan Mountains. This simulation thus covered spectral albedo ranges of clean to dirty snow with fine to coarse grains. The SMARTS calculations also enabled transformation of the diffuse fraction measured by the SPN1 to the wavelength ranges of the broadband and near-infrared PSP radiometers. Equation (4), without the V_Ω term, was applied and spectral albedos were multiplied by the spectral irradiance. Defining I as spectral radiation and E as wavelength-integrated radiation, initial values are

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{I}_{\mathrm{reflected}}^{\left(\mathrm{0}\right)}\left(\mathit{\lambda }\right)& =I_{\downarrow }\left(\mathit{\lambda }\right)\left\{{\mathit{\rho }}_{\mathrm{direct}}\left(\mathit{\lambda }\right)\left[\mathrm{1}-F_{\mathrm{dif}}\left(\mathit{\lambda }\right)\right]\right.\\ & \left.+{\mathit{\rho }}_{\mathrm{diffuse}}\left(\mathit{\lambda }\right)F_{\mathrm{dif}}\left(\mathit{\lambda }\right)\right\}\\ E_{\mathrm{reflected}}^{\left(\mathrm{0}\right)}& =\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}{I}_{\mathrm{reflected}}^{\left(\mathrm{0}\right)}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\\ {\mathit{\alpha }}^{\left(\mathrm{0}\right)}& =E_{\mathrm{reflected}}^{\left(\mathrm{0}\right)}/\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}{I}_{\downarrow }\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }.\end{array}\end{array}$$

then at iteration, the value of α increases in the following way (Fig. 4). Note that V_Ω is omitted from these iterations, because the interest lies in the change in wavelength-integrated albedo, not in the escaping radiation at each reflection. Moreover, all the radiation in the subsequent reflections is diffuse,

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{I}_{\mathrm{reflected}}^{\left(n\right)}\left(\mathit{\lambda }\right)& =I_{\mathrm{reflected}}^{\left(n-\mathrm{1}\right)}\left(\mathit{\lambda }\right){\mathit{\rho }}_{\mathrm{diffuse}}\left(\mathit{\lambda }\right)\\ E_{\mathrm{reflected}}^{\left(n\right)}& =\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}{I}_{\mathrm{reflected}}^{\left(n\right)}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\\ {\mathit{\alpha }}^{\left(n\right)}& =E_{\mathrm{reflected}}^{\left(n\right)}/E_{\mathrm{reflected}}^{\left(n-\mathrm{1}\right)}.\end{array}\end{array}$$

The assumption of a Lambertian surface versus the use of directional quantities differs in the snow literature. In this study, a Lambertian assumption is used, justified with the use of nadir-looking instruments with measurements taken midday and with the lack of directional knowledge of the re-reflected radiation. Further, as surface roughness increases, so does backscattering (Manninen et al., 2021), thereby counteracting some of the forward scattering in snow. Finally, ablation hollows, the largest surface roughness features observed, have no preferred orientation, unlike sastrugi or penitentes. These factors reduce the importance of angular effects (Painter and Dozier, 2004; Warren et al., 1998). Further, a goal of this study is to compare in situ with remotely sensed snow measurements. At the remote sensing scale, the average or sub-pixel-scale snow surface topography is usually unknown, thus the directional factors cannot be accurately computed. Although the snow-free topography may be known, the snow surface above can differ markedly, especially at fine (e.g., meter) scales.

Figure 4 Increase in broadband albedo caused by internal reflections within an ablation hollow. The x axis shows the initial albedo from Eq. (8) covering a range of grain sizes and concentration of light-absorbing particles, and the y axis shows the number of reflections from 1 to 10. The intensity shows the resulting increase in albedo from Eq. (9).

2.3 Remotely sensed measurements

Bottom-of-atmosphere (surface, Level 2A) reflectance estimates from the Sentinel-2A/B (S2) multispectral instrument were obtained. Nine bands (bands 2–7, 8a, and 11–12) were used with a spatial resolution of 20 m. To convert the narrow band surface reflectance estimates to broadband albedo, coefficients for snow-free and snow-covered surfaces, derived from radiative transfer simulations were used (Table 2 in Li et al., 2018). This surface reflectance product was processed using the Snow Property Inversion from Remote Sensing model (SPIReS, Bair et al., 2021b) to obtain fractional snow-covered area and surface properties. Broadband albedo uncertainty from S2 (0.036) was estimated based on maximum differences between acquisitions for a bare-ground target pixel, consisting of no trees, bare soil, and small shrubs. This uncertainty is close to a validation effort of S2 over dark and bright soils that showed band-wise errors up to 0.040 (Gascon et al., 2017).

The target pixel on Mammoth Mountain for comparison to the snow measured at CUES was selected because it is near CUES (2.2 km away), is at a similar elevation (CUES at 2916 m vs. target at 3041 m), has a slope of zero across the 20 m pixel, and was nearly 100 % snow-covered for 6 months, from mid-November through mid-May. It would have been preferable to select a pixel immediately adjacent to CUES, but none met those criteria. Thus, it is assumed that snow conditions and thus albedo were similar at the two sites, at least within the uncertainty of the remotely sensed and in situ broadband measurements. The mean local solar time for overpass from Sentinel-2 is 10:30, leading to times at CUES of 18:39 to 18:47 UTC. Thus, the corresponding in situ albedo measurements described in Sect. 2.2 were taken within that window of time.

2.4 Shade endmember simulations

Intrinsic snow albedo was modeled using a two-stream radiative transfer approximation coupled with Mie scattering as described in Sect. 2.2. Of note is that dust is assumed to be the predominant pollutant, based on chemical analyses from CUES (Sterle et al., 2013). Other endmembers used were an empirical snow-free background (for the remotely sensed solutions) and an ideal shade endmember with an albedo of zero across all bands (Adams et al., 1986).

2.5 In situ spectroscopy

A Spectra Vista HR-1024i was used with a Spectralon panel with 0.99 albedo over the 250–2500 nm wavelength range for irradiance measurement. The lens used has a 4^∘ field of view and was held about 1.5 m above the snow surface, leading to a footprint of about 5 cm. Measurements were made on days with clear skies, and the spectrometer was held plumb rather than slope parallel. Noise was smoothed using an 11-point sliding window fit with a local regression using a first-degree polynomial.

3 Results and discussion

An example of ablation hollows mapped by the laser scanner is shown in Fig. 5a, b. In situ albedos from CUES from the water year 2021 are shown in Fig. 6: uncorrected α_uncorrected, planar-corrected α_planar, and intrinsic α_intrinsic are based on the spatial calculations.

Figure 5 Snow with ablation hollows on 12 May 2021 at 10:45:00 PST. (a) Corresponding apparent albedo seen by the radiometer. (b) The uncorrected albedo is 0.54 (mean of what is shown). The albedo with a planar correction is 0.55, and the intrinsic albedo based on the spatial analysis is 0.61.

Figure 6 In situ albedos on Mammoth Mountain in the water year 2021. Shown are the uncorrected, planar-corrected, and intrinsic albedos for broadband (a) and near-infrared (b) wavelengths. Planar correction involved fitting a plane to the snow surface and using the solar illumination angle on that plane compared to that on a flat surface. Intrinsic albedos are derived from analyzing the view factors and illumination angles on the rough surface and using Eqs. (4) through (9) to solve for α_intrinsic. The difference between the intrinsic and planar albedos is shown in (c).

In situ and remotely sensed albedos on Mammoth Mountain from the water year 2021 are shown in Fig. 7. An unadjusted (i.e., not adjusted for shade or trees) fractional snow-covered area (f_sca), estimated with SPIReS (Bair et al., 2021b), from a nearby target pixel is also shown. The high f_sca confirms that mixed (snow and non-snow) pixel effects are minimal. An estimate of the broadband pixel albedo measured by Sentinel 2A/B (S2) is also shown, as described in Sect. 2.3. Finally, the surface roughness (in degrees, divided by 30 for scale) is plotted, also described in Sect. 2.2.

Figure 7 In situ and remotely sensed snow on Mammoth Mountain, water year 2021. Shown are uncorrected albedos measured at CUES, with the error bars (0.020) based on stated values from the manufacturer. The unadjusted (i.e., not adjusted for shade or trees) fractional snow-covered area (f_sca) from the Snow Property Inversion from Remote Sensing (SPIReS) model is shown. An estimate of the broadband pixel albedo measured by Sentinel 2A/B (S2) is shown. The error bar height (0.036) is the maximum difference in the bare-ground (no snow) reflectance.

In Fig. 6, the intrinsic albedo is usually greater than the uncorrected or planar-corrected albedo, agreeing with previous work over more limited timespans (e.g., Larue et al., 2020; Lhermitte et al., 2014; Manninen et al., 2021). The largest planar corrections appear in winter, when the planar sloped surface facing away from the sun receives the lowest irradiance relative to a flat surface. The spatial corrections are more nuanced because they involve the solar geometry and the roughness of the surface. As the days get longer in the spring, the solar zenith angle is smaller, but the rougher surface causes more variability in the view factor and illumination on each slope.

Warren et al. (1998) posited two mechanisms for albedo reduction caused by surface roughness: reduction of effective illumination angle and photon trapping. The difference α_uncorrected−α_intrinsic characterizes the combined contribution. In this study covering 110 d of the water year 2021 snow season, the differences amounted to −0.056 in the broadband albedo and −0.052 in the near infrared. Larue et al. (2020) estimate a decrease in spectral albedo at 1000 nm of −0.02 to −0.03 for low-SSA (specific surface area, i.e., large grain size) snow, but in the snow studied here with extensive ablation hollows, the magnitudes are greater. The difference $\left[I_{\mathrm{internal}}^{\left(\mathrm{0}\right)}-\left(\sum I_{\mathrm{esc}}-I_{\mathrm{esc}}^{\left(\mathrm{0}\right)}\right)\right]$ characterizes photon trapping, which accounts for a mean of 4 % of the lost broadband radiation and 5 % of the loss in the near infrared. In the late spring when the snow surface was quite rough, these losses exceeded 20 %. These result follows from Warren et al. (1998), who state that intermediate snow albedos will be most impacted by photon trapping.

The intrinsic albedo is generally greater than the planar-corrected albedo, showing that the planar correction that has been performed in previous research (Bair et al., 2018; Painter et al., 2012b) accounts for surface slope but not for roughness. But the planar correction is useful as the difference between the planar-corrected and the intrinsic albedo quantifies the impact of sub-slope surface roughness at this location. This difference implies that in areas where the average surface topography is accurately quantified (e.g., over 0.5–1.0 km pixels), a terrain-corrected (adjusted to level) surface reflectance can be used in a spectral mixture model in with a shade endmember to decrease uncertainty in impurity estimates. However, for sensors with finer resolution (e.g., ≤ 30 m), caution is advised with terrain corrections. If ground control points are not available, as in the case of many remote parts of the world, vertical errors in high-resolution elevation products approach the pixel size (Gottwald et al., 2017; Rodríguez et al., 2006; Shean et al., 2016). These errors are compounded when computing gradients (i.e., slope and aspect) needed for terrain corrections. These errors are especially noticeable for sharp features such as ridgelines. Thus, a shade endmember without any terrain correction may produce the most accurate results for these locations.

Narrow-to-broadband albedo conversions confirm that the apparent albedo is being seen from space. As surface roughness increases to its maximum during melt, albedo falls rapidly. This period coincides with the time of year when snow becomes dirtiest on the surface, as the albedo is no longer being refreshed with new snowfall. Thus, the darkening effects of surface roughness occur simultaneously with the build-up of impurities (Betterton, 2001; Rhodes et al., 1987), which presents a challenge for remote sensing. However, because impurities only affect visible through near-infrared snow albedo, and snow grain size only affects albedo in the nIR/SWIR, while shadowing affects the entire broadband spectrum, an instrument with sufficient spectral resolution and accuracy should be able to discriminate between the causes of darkening.

To test this hypothesis, SPIReS was run on S2 imagery with dirty snow endmembers and with a clean snow assumption. The resulting grain size and impurity concentration estimates were then used in the updated broadband snow albedo that now accounts for shade. Because the pixel is close to fully snow covered, this estimated albedo should be comparable to the narrow-to-broadband conversions shown in Fig. 7. The uncorrected albedo measured at CUES from Fig. 7 is plotted along with these two model runs (Fig. 8). With overlapping error bars for each scene, the resulting albedos are indistinguishable within measured error (Bair et al., 2021b). In the clean-snow run, the dust endmember is swapped for the shade endmember (Table 2). In situ spectroscopic measurements (also in Table 2) provide some validation but also illustrate the wide spatial variability of the snow surface just across the CUES study area.

Figure 8 Broadband snow albedo solutions from SPIReS compared to the uncorrected albedo measured at CUES (same as in Fig. 7). In the first set of SPIReS solutions, dirty snow endmembers are used, while in the other the snow is assumed clean. Both sets use a shade endmember. Error bars are ±0.025 (Bair et al., 2021b).

Table 2 Model solutions from SPIReS using measurements from Mammoth Mountain taken on 11 May 2021, the last two points with error bars shown in Fig. 8. The instruments are Sentinel 2B MSI (S2) and the Spectra Vista HR 1024i field spectrometer (SVC). One of the SPIReS runs used a clean-snow assumption to illustrate the difficulty in separating shade from dust endmembers (with low concentrations) with a multispectral instrument. The fractional snow-covered area (f_sca) and shade (f_shade) as well as the grain radius and dust concentration are unknowns that are solved for.

Importantly, the spectroscopic measurements show that, when used in a model, there is a consistent ability to discriminate between darkening caused by impurities and by shade. For example, despite the high spatial variability, neither the shade endmember nor the dust concentration is zero in any of the solutions. An example of dirty snow with a shaded solution is shown in Fig. 9.

Figure 9 Example of measured and modeled reflectance from field spectroscopy measurements from 12 May 2021 (Table 2). The model estimates (with an RMSE = 0.006) are f_shade = 0.31, grain radius = 454 µm, dust concentration = 77 ppm, α_apparent=0.45 (measured/modeled), and α_intrinsic=0.66 (modeled without shade).

Because the snow surface is rarely flat or level, shade needs to be accounted for, even when using measurements taken from a field spectrometer. Thus, shade needs to be included in snow albedo models, which often use lookup tables for rapid processing. Figure 10 shows the results of radiative transfer simulations to illustrate the effect of shade on the difference between intrinsic and apparent albedo.

Figure 10 Difference between intrinsic and apparent albedo versus shade fraction. The gray area represents the range of radiative transfer solutions using different combinations of grain sizes, solar zenith angles, and impurities.

There is a positive relationship: as f_shade increases, the difference between intrinsic and apparent albedo increases, but the scatter also increases. A simple adjustment is not possible; instead the lookup tables and albedo model presented in Bair et al. (2019) have been updated to include f_shade. The new albedo model estimates an apparent albedo as

$$\begin{array}{}\text{(10)}& {\mathit{\alpha }}_{\mathrm{apparent}}=f\left(r_{\mathrm{g}},\mathit{\mu },Z,{\mathrm{LAP}}_{\mathrm{name}},\mathit{\delta },\mathrm{1}-f_{\mathrm{shade}}\right),\end{array}$$

where α_apparent is the apparent albedo over three wavelength ranges (broadband, near-infrared, and visible), r_g is the grain radius in micrometers, μ is the cosine of the solar zenith angle, Z is the surface elevation in kilometers, LAP_name is the type of light-absorbing particles (dust or soot), and δ is the LAP concentration. Other properties such as an assumed midlatitude winter atmosphere are unchanged from Bair et al. (2019).

4 Conclusions

A time series of intrinsic and apparent snow albedos over a season at a sub-alpine site were presented. In situ albedo measurements were compared to those from a spaceborne multispectral sensor. The multispectral measurements and those from a field spectrometer were used in a spectral mixture model. As expected and consistent with other studies, the results show that intrinsic albedo is consistently greater than apparent albedo. Both albedos decrease rapidly as ablation hollows form during melt, combining effects of build-up of impurities on the surface and increasing roughness.

There are several conclusions with implications for remote sensing, but also in situ measurement of snow albedo. For multispectral sensors, darkening effects from snow surface roughness are significant and can easily be confused with those from impurities. In contrast, measurements from a field spectrometer have sufficient spectral resolution and accuracy to distinguish between the two effects. A spectral mixture model run on spectra obtained at a study site confirms significant darkening at the snow surface, simultaneously occurring from roughness and impurities, with wide variation spatially. In turn, a spectral mixture model was used with Sentinel 2A/B multispectral imagery assuming a clean snowpack and a dirty snowpack. Both model runs were able to match measured snow albedo with plausible solutions, but the clean snow model used the shade endmember in place of the dust endmember.

The 0.056 difference between intrinsic and apparent albedo is equivalent to the decrease in broadband albedo caused by 63 ppm dust for typical snow in spring. If the surface topography is known to the point where a plane can be fit, the difference between the intrinsic and planar-corrected albedo (mean of 0.022) could be used instead, equivalent to darkening by around 22 ppm dust. Thus, to improve uncertainty in impurity estimates, a terrain correction used in conjunction with a shade endmember in a spectral mixture model can be used for moderate resolution sensors (e.g., 0.4–1 km), but caution is advised for terrain corrections at finer resolutions (≤ 30 m) owing to elevation model errors. Generally, impurity estimates from multispectral sensors are only distinguishable from surface roughness effects for relatively dirty snow. Likewise, for a multispectral sensor, mixed pixels can be spectrally inseparable from pixels containing only dirty snow. Thus, only pixels with high snow fraction should be used for impurity estimates from a multispectral sensor (Bair et al., 2021b; Painter et al., 2012a). These conclusions were also reached by Warren (2013), but for black carbon on the snow surface in the Arctic.

This study emphasizes the difficulties in modeling lighting conditions on the snow surface. Because of these difficulties, a recommendation is to always use a shade endmember in unmixing models, even for in situ spectroscopic measurements. Likewise, snow albedo models should produce apparent albedos by accounting for the shade fraction. To this end, lookup tables and code have been revised to account for shade. The apparent albedo produced should be used in energy balance models where intrinsic albedos have been previously used.

In this study, albedos were used rather than directional reflectance quantities. The justifications are the use of nadir-looking instruments with measurements taken midday; that as surface roughness increases, so does backscattering, thereby counteracting the forward scattering in snow; and that ablation hollows, the largest surface roughness features observed, have no preferred orientation, unlike sastrugi or penitentes. These factors reduce the importance of angular effects. But the most compelling justification is that for snow, the average or sub-pixel-scale snow surface topography is usually unknown, so the directional factors cannot be accurately computed.

Future work could focus on testing these findings in other snow climates with different surface roughness features, mainly formed by wind (Filhol and Sturm, 2015). The findings about discrimination between darkening from surface roughness and impurities as well as detection limits for impurities from multispectral sensors require further testing. For example, results from dirtier snowpacks should be examined, although the size of the ablation hollows will be reduced (Lliboutry, 1964; Rhodes et al., 1987). These findings highlight the need for hyperspectral measurements of snow from aerial and spaceborne sensors. The NASA Earth Observing-1 Hyperion was promising in this regard, but lack of coverage, repeat passes, or a surface reflectance product limited utility. The upcoming NASA Surface Biology and Geology (SBG) and ESA Copernicus Hyperspectral Imaging Mission for the Environment (CHIME) spaceborne spectrometers may offer chances to test these findings using spectroscopic measurements from space.