Case study of spatial and temporal variability of snow cover, grain size, albedo and radiative forcing in the Sierra Nevada and Rocky Mountain snowpack derived from imaging spectroscopy

Abstract

Quantifying the spatial distribution and temporal change in mountain snow cover, microphysical and optical properties is important to improve our understanding of the local energy balance and the related snowmelt and hydrological processes. In this paper, we analyze changes of snow cover, optical-equivalent snow grain size (radius), snow albedo and radiative forcing by light-absorbing impurities in snow and ice (LAISI) with respect to terrain elevation and aspect at multiple dates during the snowmelt period. These snow properties are derived from the NASA/JPL Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) data from 2009 in California's Sierra Nevada and from 2011 in Colorado's Rocky Mountains, USA.

Our results show a linearly decreasing snow cover during the ablation period in May and June in the Rocky Mountains and a snowfall-driven change in snow cover in the Sierra Nevada between February and May. At the same time, the snow grain size is increasing primarily at higher elevations and north-facing slopes from 200 microns to 800 microns on average. We find that intense snowmelt renders the mean grain size almost invariant with respect to elevation and aspect. Our results confirm the inverse relationship between snow albedo and grain size, as well as between snow albedo and radiative forcing by LAISI. At both study sites, the mean snow albedo value decreases from approximately 0.7 to 0.5 during the ablation period. The mean snow grain size increased from approximately 150 to 650 microns. The mean radiative forcing increases from 20 W m${}^{-\mathrm{2}}$ up to 200 W m${}^{-\mathrm{2}}$ during the ablation period. The variability of snow albedo and grain size decreases in general with the progression of the ablation period. The spatial variability of the snow albedo and grain size decreases through the melt season while the spatial variability of radiative forcing remains constant.

Introduction

Snow is important to the climate and hydrology of mountain and arctic regions, as well as to water resources management. Snow influences the Earth's energy balance through its strong albedo feedback. Changes to the snow cover extent, snow albedo and snow water equivalent (SWE) therefore have direct impacts on the environment, water availability and economics. In the western United States, seasonal snowmelt provides more than 70 % of water supply. The Sierra Nevada in California and the Rocky Mountains in Colorado are examples of two mountain ranges where reservoirs store snowmelt for use in late spring and throughout the summer and early fall when precipitation is less frequent. As such, these regions are of great interest to snow research and water management .

Snowmelt in midlatitude mountains is largely dominated by the absorbed solar radiation . The net shortwave radiation is a function of the solar zenith, cloud cover, aerosol attenuation, surface albedo as well as slope and aspect of the terrain. The solar irradiance at the surface is highly variable in mountainous areas mainly due to the complex terrain, cloud cover and vegetation above the snow surface . Besides snowmelt, snowfall impacts snowpack properties very rapidly and, thus, it is important to capture temporal changes in addition to spatial patterns. For example, showed that spatially distributed snowmelt models benefit significantly from incorporating snow albedo estimates derived from the NASA/JPL Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) data instead of using empirical basin-wide snow albedo assumptions.

Snow microstructure is a indicator of physical processes occurring in a snowpack. Snow grain size is generally used as a proxy for the microstructure due to the lack of a direct measurement technique. Throughout this paper, we refer to the optical-equivalent snow grain radius ($r$), or simply snow grain size for brevity. It can be estimated with handheld spectroscopy instruments as well as airborne imaging spectroscopy . The former can be used to resolve the vertical distribution of snow grain sizes within a snowpack, and the latter is very useful to resolve the spatial and temporal distribution of not only snow grain size but also albedo and radiative effects on scales useful for snow hydrology. It should be noted that $r$ should not be directly compared with the snow grain size derived using other measurement techniques, such as the traditional visually estimated snow grain size or the effective radius of a grain size distribution corresponding to a specific surface area (SSA). and suggest a conversion of $r$ $=$ $\mathrm{3}/({\mathit{\rho }}_{\text{ice}}$ SSA), where ${\mathit{\rho }}_{\text{ice}}$ is the density of ice (917 kg m${}^{-\mathrm{2}}$). An optical-equivalent snow grain radius of 100 and 650 $\mathrm{\mu }$m corresponds therefore to an SSA of 33 and 5 m${}^{\mathrm{2}}$ kg${}^{-\mathrm{1}}$, respectively.

Snow albedo is a key parameter to the Earth's energy balance. The snow albedo feedback mechanism has a strong impact on climate on all temporal and spatial scales. In particular, snow albedo is a dominant driver of snowmelt, ecohydrological processes and water supplies. Snow albedo is mainly changed by snowfall, metamorphism, surface roughness and light-absorbing impurities in snow and ice (LAISI), such as dust and black carbon (BC). While fresh snow with small grains has a high albedo value, aged snow has a lower albedo value in conjunction with snow grain size growth and change in optical path length. LAISI, such as black carbon, dust and tree litter, have the most profound impacts on snow albedo .

Geographic extent of the AVIRIS image data. The study area is shown hatched. It is based on a mosaic of subsets, such that a common area exists for all flight dates. The elevation model is hill shaded in the background with an illumination corresponding to the mean solar position for all dates.

[Figure omitted. See PDF]

This “darkening” of snow increases the absorption of solar irradiance, which reduces the cold content of the snowpack or enhances snowmelt if its temperature is at 0 ${}^{\circ }$C. The related consequences of LAISI to snow hydrological and climatological processes became more evident by the work of , , , and others. Novel remote sensing algorithms were developed by to retrieve the radiative forcing by LAISI at the snow surface, hereafter referred to as simply radiative forcing.

The objective of this research is to use imaging spectroscopy to identify the spatiotemporal distribution of quantitative snow properties at study sites in the Sierra Nevada and Rocky Mountains. AVIRIS data and the study sites are described below in Sect. , followed by our retrieval algorithms in Sect. . Results are presented in Sect.  with the most interesting findings discussed in Sect.  about the relationship between elevation and grain size for melting snow.

California and Colorado study areas and data description

Terrain

We use the United States Geological Survey National Elevation Dataset (NED) with 1 arc second spatial resolution (roughly 30 m) as provided by the AVIRIS data product. A hill shaded version is used as background to Fig. . The Sierra Nevada study region covers a larger distribution of elevations with a vertical range of 2400 m while the Colorado study area has a vertical range of only 1200 m (Fig. ). The Colorado study area reaches higher elevations with a maximum at 4111 m, while the maximum elevation in the California study area is only 3533 m. The mean elevation is over 1000 m higher in the Colorado study area (Table ). The Kaweah and Kings river basins in Sierra Nevada study area drain to the west as seen in the distribution of aspects with a peak around 270${}^{\circ }$ in Fig. . The Uncompahgre River basin drains to the north – hence the lower percentage of 180${}^{\circ }$ slopes. A larger proportion of the Uncompahgre slopes face east and west, draining toward the center of the basin and explaining the larger percentage of the basin with aspects of 90 and 270${}^{\circ }$.

Data parameters of the study areas derived as a mosaic from the subsets of flight lines (see Table  and Fig. ).

Sierra Nevada, CA, 2009
Dates	27 Feb, 5 Mar, 4 Apr, 7 May, 8 May
Pixel size $($m$)$	14.5
Image area $($km${}^{\mathrm{2}}$, pixels$)$	252.03067, 1.2 $\times$ 10${}^{\mathrm{6}}$
Elevation min./mean/max. $($m$)$	1140/2575/3533
Rocky Mountains, CO, 2011
Dates	13 May, 16 May, 9 Jun, 15 Jun, 23 Jun
Pixel size $($m$)$	13.8
Image area $($km${}^{\mathrm{2}}$, pixels$)$	54.721982, 2.9 $\times$ 10${}^{\mathrm{5}}$
Elevation min./mean/max. $($m$)$	2950/3605/4111

Snow hydrology

In the Sierra Nevada of California, the study area covers part of the Kings and Kaweah river basins. In the Rocky Mountains of Colorado, the study area covers the Senator Beck drainage, which is part of the Uncompahgre watershed . The Kings and Kaweah have maritime snowpacks (typically warmer and deeper) while the Uncompahgre has a continental mountain snowpack (typically colder and more shallow) , allowing us to test our algorithms in different snow regimes.

We use snow pillow data of SWE at high and low elevation as a proxy for snow accumulation or snowmelt as shown in Fig. . We choose the Big Meadow (36.72${}^{\circ }$ N, 118.84${}^{\circ }$ W; 2300 m, Kings) and the Farewell Gap (36.41${}^{\circ }$ N, 118.58${}^{\circ }$ W; 2900 m, Kaweah) stations to represent the California study area. For Colorado, we choose the Idarado (37.93${}^{\circ }$ N, 107.68${}^{\circ }$ W; 3000 m, Uncompahgre) and the Red Mountain Pass (37.90${}^{\circ }$ N, 107.72${}^{\circ }$ W; 3400 m, Uncompahgre) SNOTEL stations. Based on those SWE data, 2009 was an average year in the Sierra Nevada and 2011 was an above average year in the Rocky Mountains. Peak SWE and length of the ablation period were fairly similar at the two high-elevation stations in each area, but separated by a month. At Farewell Gap, SWE peaked at approximately 80 cm on 17 April 2009 and snow melt out occurred on 19 May 2009. At Red Mountain Pass, SWE peaked at approximately 85 cm on 25 May 2011 and melt out occurred on 20 June 2011.

Histograms of the terrain aspect (top panels) and elevation (bottom panels) for the two study sites, Kings/Kaweah Basin, Sierra Nevada, CA (left panels), and Uncompahgre Basin, Rocky Mountains, CO (right panels). The values are normalized by the size of the study area given in Table . The peak occurrences of terrain aspect give an indication to the main drainage direction of a river basin.

[Figure omitted. See PDF]

In addition to regional and interannual variability, there was some distinct differences in melt patterns. In Kings/Kaweah multiple snowfall events refreshed the snow surface over the period of flight, delaying the main snowmelt pulse until just prior to snow melt out. In the Uncompahgre melt initiated just after peak SWE and steadily increased over the ablation period, with the exception of a minor snowfall event at the end of May. While the measurement sites may not be wholly representative of the two study areas, they are helpful in interpreting some of the spatial and temporal variability observed and discussed in Sect. .

Snow accumulation (blue) and snowmelt (red) indicated with the differential snow water equivalent ($\mathrm{\Delta }$SWE) with respect to time. The AVIRIS flight dates are marked with black vertical lines and a symbol indicating the ER-2 or Twin Otter airborne platform. The upper plot shows snow sensor pillow data from two stations representing a higher and a lower location in the Sierra Nevada, CA, study area in 2009. The lower plot shows the same information, but from stations close to the Rocky Mountains, CO, study area in 2011.

[Figure omitted. See PDF]

Dust concentrations in snow

The Senator Beck Basin is a research catchment established to monitor the hydrological impacts of regularly observed episodic events of dust deposition on snow in the Upper Colorado River basin . Dust is particularly effective at altering the snow energy balance in this region because the majority of mass is deposited in the spring coinciding with snowmelt onset (80 % of dust deposition events occur March–May) and, additionally, dust is inefficiently scavenged by melt water and layers coalesce at the surface as snowmelt, compounding albedo decay. Over the 10-year record (2005–2015) dust concentrations toward the end of the snow season have ranged from 0.8 to 4.8 mg g${}^{-\mathrm{1}}$ in the top 30 cm of the snow column and exhibit a bimodal distribution, with about 1 mg g${}^{-\mathrm{1}}$ in lesser dust years and over 4 mg g${}^{-\mathrm{1}}$ in extreme dust years . There were nine observed dust events in 2011 after 1 March and the end of year dust concentration was 1.3 mg g${}^{-\mathrm{1}}$, making it a lesser dust year. Black carbon has only been monitored across the 2013 ablation period, and concentrations were so low ($<$ 1 to 20 ng g${}^{-\mathrm{1}}$) that it was concluded that its contribution to radiative forcing in this region is negligible in the presence of such heavy dust loading .

Unfortunately, similar dust deposition data do not exist in the Sierra Nevada, but given the predominant wind patterns and location of dust emission sources in the western US, dust deposition is expected to be relatively low. There was a study conducted in the eastern Sierra Nevada in 2009 by which found surface black carbon concentrations of $<$ 1 to 75 ng g${}^{-\mathrm{1}}$ and dust concentrations of 1 to 44 $\mathrm{\mu }$g g${}^{-\mathrm{1}}$ over the ablation period. Despite the lower concentrations sampled in the Sierra Nevada relative to the Rocky Mountains, this study also concluded that dust dominated the radiative forcing.

Imaging spectroscopy data

Results provided in this paper were derived from georeferenced at-sensor radiances acquired by AVIRIS flown on NASA's ER-2 high-altitude aircraft and on a Twin Otter aircraft. Table  shows specific information for each flight line including aircraft altitude and heading, terrain elevation, pixel size, solar geometry and columnar water vapor retrieved by using the Atmospheric/Topographic Correction for Airborne Imagery (ATCOR-4) software . The pixel size is directly related to the terrain height and the aircraft altitude, which varied significantly between acquisition dates. We therefore resampled the AVIRIS data to the common coarsest resolutions of 14.5 and 13.8 m for California and Colorado, respectively. For each acquisition date, the images are mosaicked and cropped to their common area (252 km${}^{\mathrm{2}}$ in California and 54.7 km${}^{\mathrm{2}}$ in Colorado). This enables intercomparisons between acquisition dates at the pixel level. Figure  shows the outline of each flight line overlaid on the terrain for each study area (intersection of all flight dates) as a hatched line.

Remote sensing retrievals of quantitative snow properties

The quantitative snow properties shown and analyzed in this paper are extracted from imaging spectrometer data using retrievals based on a legacy of algorithms, including but not limited to , , , and . Our retrieval algorithms are described in details in and summarized in this section for reference.

Summary of the individual AVIRIS scenes used in this paper. The scene number corresponds to the AVIRIS “run ID”. The aircraft altitude, heading and the terrain elevation are averages per scene. The solar geometry is valid at the center of a scene. Water vapor is retrieved using ATCOR-4 (see Sect. ).

AVIRIS		Aircraft		Mean	Pixel	Solar		Water
Flight	Scene	Altitude	Heading	elevation	size	Zenith	Azimuth	vapor
date	nr.	$($km$)$	$($deg.$)$	$($km$)$	$($m$)$	$($deg.$)$	$($deg.$)$	$($cm$)$
Sierra Nevada, CA, 2009
27 Feb	6	19.799	354	2.547	14.1	48.3	205.7	0.32
	8	19.836	185	2.925	13.9	56.1	205.4	0.23
5 Mar	5	19.806	181	2.056	14.4	45.7	155.3	0.22
	6	19.765	000	2.554	14.0	44.7	158.7	0.13
	7	19.773	181	2.879	13.7	44.1	162.7	0.09
4 Apr	5	19.799	182	2.082	14.4	41.0	131.6	0.16
	6	19.674	359	2.552	13.9	39.5	134.4	0.10
	7	19.902	181	2.911	13.8	38.5	137.3	0.07
	8	19.796	000	1.865	14.5	37.4	139.7	0.24
7 May	6	14.262	348	3.020	9.0	23.4	144.0	0.30
	7	14.274	193	2.650	9.3	22.4	148.6	0.43
	8	14.273	349	2.277	9.3	21.5	154.4	0.62
8 May	5	14.141	343	3.178	8.9	27.6	129.0	0.36
	6	14.154	194	2.642	9.2	26.4	132.2	0.55
	7	14.154	345	2.383	9.2	25.0	136.4	0.73
	8	14.144	195	2.035	9.6	23.9	140.7	0.89
	9	14.157	346	2.715	9.2	22.6	146.1	0.58
Rocky Mountains, CO, 2011
13 May	5	6.496	352	3.560	2.1	19.6	173.0	0.25
	6	6.536	354	3.526	2.1	19.5	182.9	0.25
	7	6.533	354	3.683	2.1	19.9	192.6	0.21
16 May	7	19.866	180	3.237	13.5	20.5	153.5	0.22
9 Jun	12	20.170	179	3.260	13.7	15.0	179.1	0.33
15 Jun	12	20.216	180	3.196	13.8	20.8	129.7	0.51
23 Jun	5	20.139	001	3.271	13.7	38.3	102.0	0.54

Directional reflectance

Reflected and scattered solar radiation are observed by the airborne spectrometer and transformed into calibrated radiance for each wavelength band ${\mathit{\lambda }}_i$, based on band-wise calibration. The retrieval of surface properties requires that we convert the measured radiance first to reflectance at the sensor level $R_{\text{sensor}}^{\text{obs}}$ and then to surface directional reflectance by compensation of atmospheric scattering and absorption. We use the ATCOR-4 software to model the reflectance of the atmosphere at the sensor level $R_{\text{atm}}^{\text{mdl}}$, the total transmittance $T^{\updownarrow ,\text{mdl}}$ and the spherical albedo at the surface level $S_{\text{sfc}}^{\text{mdl}}$ to derive the spectral hemispherical–directional reflectance factor at the surface level from $R_{\text{sensor}}^{\text{obs}}$ : $${\text{HDRF}}_{\text{sfc}}^{\text{obs}}={\displaystyle \frac{R_{\text{sensor}}^{\text{obs}}-R_{\text{atm}}^{\text{mdl}}}{T^{\updownarrow ,\text{mdl}}+S_{\text{sfc}}^{\text{mdl}}\left(R_{\text{sensor}}^{\text{obs}},-,R_{\text{atm}}^{\text{mdl}}\right)}}.$$

Note that the rugged terrain in both study areas (see Sect. ) requires us to account for the local solar and viewing angles at each terrain facet. Local zenith angles are derived by $${\mathit{\theta }}_{[\odot ,v]}^{\prime }=\arccos \left(\cos ,{\mathit{\theta }}_{[\odot ,v]},\cdot ,\cos ,\text{slp},+,\sin ,{\mathit{\theta }}_{[\odot ,v]},\cdot ,\sin ,\text{slp},\cdot ,\cos ,{\mathit{\phi }}^{\prime }\right),$$ where ${\mathit{\phi }}^{\prime }$ $=$ $\mathit{\phi }$ $-$ asp, and slp and asp denote the terrain's slope and aspect, respectively. All angles are defined in degrees. Although the following geometries are locally defined with respect to the terrain, we omit its notation for readability. Note that we mask pixels affected by terrain shadows in our retrievals, except for the snow cover determination.

Snow cover

A byproduct derived from HDRF${}_{\text{sfc}}^{\text{obs}}$ is the snow cover map. We apply the normalized difference snow index (NDSI): $$\text{NDSI}\equiv {\displaystyle \frac{{\text{HDRF}}_{\text{sfc}}^{\text{obs}}\left({\mathit{\lambda }}_{\text{VIS}}\right)-{\text{HDRF}}_{\text{sfc}}^{\text{obs}}\left({\mathit{\lambda }}_{\text{SWIR}}\right)}{{\text{HDRF}}_{\text{sfc}}^{\text{obs}}\left({\mathit{\lambda }}_{\text{VIS}}\right)+{\text{HDRF}}_{\text{sfc}}^{\text{obs}}\left({\mathit{\lambda }}_{\text{SWIR}}\right)}},$$ where ${\mathit{\lambda }}_{\text{VIS}}$ $=$ 0.668 $\mathrm{\mu }$m and ${\mathit{\lambda }}_{\text{SWIR}}$ $=$ 1.502 $\mathrm{\mu }$m . To mask partially covered snow pixels, we empirically determined NDSI thresholds of 0.90 for the Sierra Nevada and 0.93 for the Rocky Mountains by visual interpretation with the corresponding true-color images. The high spatial resolution of the AVIRIS imagery (see Table ) allows the use of this simple technique without rejecting too many pixels. The multiplication of the number of remaining snow pixels by the pixel size yields the approximate (low-biased) snow-covered area. We sum the snow-covered area and normalize it with the study area to report the approximate percent snow cover.

Snow grain size

The optical-equivalent snow grain size, reported here as radius $r$, is an important snow property and prerequisite for the snow albedo retrieval. Based on , we derive the snow grain size by searching for the best fit between the observed HDRF${}_{\text{snow}}^{\text{obs}}$ and modeled HDRF${}_{\text{snow}}^{\text{mdl}}$ snow spectra with the minimum sum of absolute differences: $$\begin{array}{cc}r& :\left[\text{min},\sum _{i=1.17\mathrm{\mu }\text{m}}^{1.27\mathrm{\mu }\text{m}},\left|{\text{HDRF}}_{\text{snow}}^{\text{mdl}},\left({\mathit{\theta }}_{\odot },,,{\mathit{\theta }}_v,,,\mathit{\phi },;,r,;,{\mathit{\lambda }}_i\right)\right.\right.\\ & & \left..\left..-,{\text{HDRF}}_{\text{snow}}^{\text{obs}},\left({\mathit{\theta }}_{\odot },,,{\mathit{\theta }}_v,,,\mathit{\phi },;,r,;,{\mathit{\lambda }}_i\right)\right|\right].\end{array}$$ We assume to observe a single layer with homogeneous spherical snow grains. Using the ice absorption feature centered around 1.27 $\mathrm{\mu }$m allows us to probe the uppermost layer in the snowpack with a depth in the order of 1 cm . HDRF${}_{\text{snow}}^{\text{mdl}}$ is modeled using the radiative transfer code DISORT . The two spectra are matched at the corresponding solar and viewing geometry for each pixel in the terrain. The fit is performed at AVIRIS band numbers 86 to 99, which correspond to the wavelength interval $\mathit{\lambda }$ $=$ [1.17, 1.27 $\mathrm{\mu }$m]. The use of these bands reduces the bias by liquid water on light absorption by frozen water. This simple yet powerful approach allows more robust retrievals of the optical snow grain size under snowmelt conditions . Nevertheless, strong snowmelt with percolation of meltwater is assumed to introduce yet unknown bias to the snow grain size retrieval. More detailed studies of optical properties of melting snow would be required, such as .

Snow albedo

We derive the spectral snow albedo ${\mathit{\alpha }}_{\text{snow}}^{\text{obs}}$ from HDRF${}_{\text{snow}}^{\text{obs}}$ by applying the anisotropy factor, which is given by the ratio between the modeled spectral albedo and the spectral HDRF${}_{\text{snow}}^{\text{mdl}}$: $${\mathit{\alpha }}_{\text{snow}}^{\text{obs}}(r,\mathit{\lambda })={\displaystyle \frac{{\text{HDRF}}_{\text{snow}}^{\text{obs}}\left({\mathit{\theta }}_{\odot },,,{\mathit{\theta }}_v,,,\mathit{\phi },;,r,;,\mathit{\lambda }\right){\mathit{\alpha }}_{\text{snow}}^{\text{mdl}}(r,\mathit{\lambda })}{{\text{HDRF}}_{\text{snow}}^{\text{mdl}}\left({\mathit{\theta }}_{\odot },,,{\mathit{\theta }}_v,,,\mathit{\phi },;,r,;,\mathit{\lambda }\right)}}.$$ The integration of Eq. () normalized by the topographically modulated spectral total (direct $+$ diffuse) irradiance at the surface level $E_{\text{sfc}}^{\text{mdl}}(\mathit{\lambda }$, ${\mathit{\theta }}_{\odot })$ derives the broadband (spectrally integrated) albedo as a function of the snow grain size $r$: $${\mathit{\alpha }}_{\text{snow}}^{\text{obs}}\left(r\right)={\displaystyle \frac{\underset{\mathit{\lambda }=0.36\mathrm{\mu }\text{m}}{\overset{2.5\mathrm{\mu }\text{m}}{\int }}{E}_{\text{sfc}}^{\text{mdl}}\left(\mathit{\lambda },,,{\mathit{\theta }}_{\odot }\right){\mathit{\alpha }}_{\text{snow}}^{\text{obs}}(r,\mathit{\lambda })\mathrm{\Delta }\mathit{\lambda }}{\underset{\mathit{\lambda }=0.36\mathrm{\mu }\text{m}}{\overset{2.5\mathrm{\mu }\text{m}}{\int }}{E}_{\text{sfc}}^{\text{mdl}}\left(\mathit{\lambda },,,{\mathit{\theta }}_{\odot }\right)\mathrm{\Delta }\mathit{\lambda }}}.$$ $E_{\text{sfc}}^{\text{mdl}}(\mathit{\lambda }$, ${\mathit{\theta }}_{\odot })$ is provided with the direct and the diffuse component by the ATCOR-4 software . Note that albedo values must range between 0 and 1, while HDRF can range up to 1.4 or even more.

Snow radiative forcing by light-absorbing impurities

We calculate the radiative forcing by LAISI (RF${}_{\text{snow}}$) based on the difference between the observed spectral albedo and a modeled spectral clean snow albedo of the same grain radius . This reduction of spectral snow albedo by LAISI is multiplied by the direct and diffuse instantaneous spectral solar irradiance entering the snow surface and integrated over the given spectrum, such that the additional radiative flux into the snow pack is given by $$\begin{array}{cc}{\text{RF}}_{\text{snow}}& =\underset{\mathit{\lambda }=0.36\mathrm{\mu }\text{m}}{\overset{1.08\mathrm{\mu }\text{m}}{\int }}{E}_{\text{sfc}}^{\text{mdl}}\left(\mathit{\lambda },,,{\mathit{\theta }}_{\odot }\right)\left[{\mathit{\alpha }}_{\text{snow}}^{\text{mdl}},(,r,,,\mathit{\lambda },)\right.\\ & & \left..-,{\mathit{\alpha }}_{\text{snow}}^{\text{obs}},(,r,,,\mathit{\lambda },),c,\left(r,,,{\mathit{\lambda }}_c\right)\right]\mathrm{\Delta }\mathit{\lambda },\end{array}$$ where $c(r$, ${\mathit{\lambda }}_c)$ $=$ ${\mathit{\alpha }}_{\text{snow}}^{\text{mdl}}(r$, ${\mathit{\lambda }}_c)/{\mathit{\alpha }}_{\text{snow}}^{\text{obs}}(r$, ${\mathit{\lambda }}_c)$ denotes a scalar to adjust the modeled spectra to the measurements for a potential bias by terrain modeling imperfections. $c$ is derived at ${\mathit{\lambda }}_c$ $=$ 1.08 $\mathrm{\mu }$m where the influence of dust is relatively small .

Results

Figures  to  show maps of snow grain size, albedo and radiative forcing in the Sierra Nevada (upper panel) and Rocky Mountains (lower panel) derived using algorithms described in Sect.  and . The three dates represent an early, middle and late observation during the melt period (see Fig. ). These results are further analyzed in Sect.  regarding the temporal changes and in Sects.  to  regarding the spatial distribution and correlation with respect to elevation and aspect of the terrain.

Maps of snow grain size retrieved using Eq. () for three dates during the melt period in the Sierra Nevada (upper panel), and the Rocky Mountains (lower panel). The upper and lower map represent 252 and 55 km${}^{\mathrm{2}}$, respectively (Table ).

[Figure omitted. See PDF]

Temporal changes of regional optical snow properties

In general, the snow line is significantly lower in the Sierra Nevada study area as compared to the Rocky Mountains because of the maritime snowpack. However, interannual variability also modulates the snow-covered area. We calculate the snow-covered area of the study area and normalize it by the basin area to serve as a proxy for accumulation and melt (Sect. ). Figure  shows that the sequence of images in Kings/Kaweah captures the transition from an accumulating to melting snowpack, while it is mainly the melt period in the Uncompahgre (refer also to Fig. ). We can observe a difference between the two areas in the rate of snow cover change with respect to time, which is an indicator for the intensity of snowmelt. This is likely related to both the more frequent snowfall in the Kings/Kaweah in April and difference in solar irradiance between the March–April and May–June time periods.

Same as Fig.  but for snow albedo.

[Figure omitted. See PDF]

Same as Fig.  but for snow radiative forcing.

[Figure omitted. See PDF]

Time series of snow cover, albedo, grain size and impurity radiative forcing in the Sierra Nevada, CA, USA, in 2009 (upper panel) and in the Rocky Mountains, CO, USA, in 2011 (lower panel). The snow cover is normalized by the size of the study area and the other values represent the average value per date surrounded by the corresponding standard deviation.

[Figure omitted. See PDF]

A snowfall event between the first two airborne observations in the Sierra Nevada (Fig. ) decreased the mean snow grain size from approximately 330 to 180 $\mathrm{\mu }$m (Fig. ). The grain size varied little between the second and third observation, after which it increased by the fourth and fifth observation to 650 $\mathrm{\mu }$m within a month of time. We know from the SWE measurements at the surface sites that the transition to a melting snowpack likely occurred after the second airborne observation (Fig. ). See Sect.  for a discussion on the inverse relationship between snow grain size and elevation under intense snowmelt conditions.

Patterns in mean broadband albedo are similar to those observed for grain size. The Sierra Nevada study site remained constant at around 0.7 due to late season snowfall until April and then decreased to 0.55 by the last two acquisitions in May, when the snowpack was melting. In the Rocky Mountain study area the mean snow albedo decreased from 0.75 in May to 0.5 in late June. In both areas the largest decreases were observed in May, when both sites were likely transitioning to isothermal snowpack conditions (at 0 ${}^{\circ }$C) and melting.

The time series of radiative forcing are very similar for the two study areas, except that the values in the Rocky Mountains are roughly 20 W m${}^{-\mathrm{2}}$ larger. This is consistent with the high dust concentrations in the Rocky Mountains in combination with the increased solar irradiance later in the season.

Spatial distribution of snow grain size

Snow grain sizes at top of the snow pack are affected by snowfall, metamorphism, wind effects and melt processes. We would expect to find smaller snow grains at higher elevation and north-facing slopes. With time, we also expect snow surface grain size to grow at varying rates but decrease abruptly with the addition of snowfall. The analysis of our retrievals using Eq. () is consistent with this common theory. The two-dimensional histograms in Figs.  and  show the relation of grain size with elevation and with terrain aspect.

For both study areas we find, for the early ablation period, that the snow grain size at high elevations is less than 200 $\mathrm{\mu }$m and at low elevations is 500 to 800 $\mathrm{\mu }$m. The snow grain size is inversely correlated with elevation in the May data of the Rocky Mountains. Snow grain size gradually increases over time at higher elevation as shown in the left panels of Figs.  and . Late in the ablation period the snow grain size is nearly constant with elevation, ranging from 500 to 700 $\mathrm{\mu }$m and becomes positively, but less sensitively, correlated with elevation. Potential reasons for this remarkable change over time are discussed in Sect. . The snow grain size distribution with respect to the terrain orientation (aspect) is shown in the right panels of Figs.  and . As expected, we retrieve small snow grain sizes, on the order of about 100 $\mathrm{\mu }$m on north-facing slopes, versus 400 $\mathrm{\mu }$m or more on southeast slopes. On 5 March 2009, just after a snow event, the snow grain sizes are significantly smaller and have relatively low spatial variability, especially at higher elevations, especially at higher elevations. More pixels at lower elevations were covered with snow (see Fig. ).

Two-dimensional histogram of snow grain size with respect to terrain elevation (left panels) and aspect (right panels) in the Sierra Nevada in spring 2009.

[Figure omitted. See PDF]

Same as Fig.  but for the Rocky Mountains in spring 2011. Note: a rare retrieval artifact is found on 15 June where Eq. () has no unique solution within a snow grain size range of 540 to 570  $\mathrm{\mu }$m, causing a gap in the corresponding plots.

[Figure omitted. See PDF]

Snow albedo along a north–south (left panels) and west–east transect (right panels) through the middle of the observed Sierra Nevada area. The values are smoothed with a boxcar average of 30 pixels. The snow-free surface is shown in black.

[Figure omitted. See PDF]

Spatial distribution of snow albedo

Snow albedo is correlated with elevation and terrain aspect early in the ablation period at both study areas (see Fig.  through Fig. ). Our retrievals in the early ablation period show albedo values of 0.6 to 0.8 at higher elevations and 0.5 to 0.6 at lower elevations for both study areas. A rather large spatial variance is observed in February, March and April in the Sierra Nevada as shown in Fig. () with albedo varying from 0.5 to almost 0.9 within less than 1 km. In each case the higher albedo values are found on north- and steep west-facing slopes. Later in the ablation period, albedo values decreased to approximately 0.5 and 0.4 at high and low elevations, respectively. In general the albedo increase with elevation is more distinct at the Sierra Nevada study area than in the Rocky Mountain study area. A strong dependence of snow albedo with terrain slope is found in February and March in the Sierra Nevada. North-facing slopes show values ranging between 0.8 and 0.9, whereas south-facing slopes are clearly darker with approximately 0.6 albedo values (Figs.  and ). We found similar, although less distinct, dependencies in the Rocky Mountains (Figs.  and ). This is likely because northern slopes experience larger relative changes in irradiance during spring due to a decrease in solar zenith angle and the related reduction of terrain shadows. The Rocky Mountain study area observations were started in May, so the terrain effects are less prominent than those of the Sierra Nevada study area, which were started in February. During the late ablation period the differences in snow albedo between north- and south-facing slopes was less than 0.1.

In both time series the spatial variance in the retrieved albedo corresponds with the spatial variance in retrieved snow grain sizes (see Fig. ). For example, the retrievals on 5 March 2009, just after the snowfall event, show an increased variance in both albedo and grain size (see Fig. ).

Spatial distribution of radiative forcing

We find the radiative forcing increased during ablation from almost 0 to spatial mean values of 170 W m${}^{-\mathrm{2}}$ in the Sierra Nevada and 190 W m${}^{-\mathrm{2}}$ in the Rocky Mountains (Figs. , , and ). The standard deviations are in the order of 60 W m${}^{-\mathrm{2}}$ in both study areas across the full period. The Colorado data were acquired later in spring than the California data and thus with more solar irradiance. A direct quantitative comparison of radiative forcing between the two data sets is therefore difficult. The higher elevations in Colorado contribute to a later onset of snowmelt as compared to the Sierra Nevada, which leads to a later exposure of LAISI at the surface of the snowpack. We should bear in mind that each year is different in terms of dust deposition and exposure at the surface . Nevertheless a qualitative comparison with changes in SWE (Fig. ), as an indicator for snow accumulation and snowmelt, shows a relationship of longer snowmelt periods with an increase in radiative forcing. This could be dominated by the SWE magnitude with a longer time period over which impurities can be deposited and then accumulated at the surface. It could also be due to higher impurity content in the Rocky Mountains or a combination of both. In 2011 only three of the observed dust events occurred between the first and last flights in the Rocky Mountains, but all of the dust deposited over the full season would have surfaced as snow started to melt, likely contributing to the higher radiative forcings in this area.

Same as Fig.  but for the Rocky Mountains.

[Figure omitted. See PDF]

Snow albedo with respect to terrain elevation (left panels) and terrain aspect (right panels) in the Sierra Nevada in spring 2009.

[Figure omitted. See PDF]

Same as Fig.  but for the Rocky Mountains in spring 2011.

[Figure omitted. See PDF]

Radiative forcing of snow impurities with respect to terrain elevation (left panels) and aspect (right panels) in the Sierra Nevada in spring 2009.

[Figure omitted. See PDF]

Same as Fig.  but for the Rocky Mountains in spring 2011.

[Figure omitted. See PDF]

In California, we found almost no dependence of radiative forcing on elevation (Fig. ). In contrast, the radiative forcing retrieved in Colorado is clearly correlated with respect to elevation with a distinct trend of larger radiative forcing at lower elevations (Fig. ). Column 2 of Figs.  and  shows the radiative forcing with respect to the terrain aspect. In the Northern Hemisphere, south-facing slopes experience stronger irradiance and decreased snow albedo values (Sect. ). The combination of both of these factors, in addition to LAISI deposition, led to increased radiative forcing. For example, in the Sierra Nevada, the radiative forcing is close to 0 until May at northern slopes and up to about 150 W m${}^{-\mathrm{2}}$ at south-southeastern aspects. In May, we retrieve up to 150 W m${}^{-\mathrm{2}}$ at northern slopes and about 200 W m${}^{-\mathrm{2}}$ at southern-oriented terrain slopes (Fig. ). The analysis of the Rocky Mountain study area shows a similar trend over time, but with radiative forcing values 50 to 100 W m${}^{-\mathrm{2}}$ larger than in the Sierra Nevada (Fig. ). We also find that the dependence of terrain aspect angle with radiative forcing decreases until the summer solstice due to the decreasing solar zenith angle and the corresponding irradiance increase on steep slopes in mountainous terrain. Our results show radiative forcing differences between north- and south-facing slopes being smoothed out over time. An exception is found in late June in Colorado with significantly larger radiative forcing on southeast slopes of up to 300 W m${}^{-\mathrm{2}}$. This is presumably due to the earlier melt out of south to southwest slopes and may also be due to dust transport effects.

Uncertainty and error analysis

reported baseline retrieval uncertainties based on the AVIRIS sensor performance (noise-equivalent changes in radiance or NE$\mathrm{\Delta }$L) of $\pm$0.001, $\pm$3.6 $\mathrm{\mu }$m (at the 1.03 $\mathrm{\mu }$m ice absorption feature) and $\pm$0.2 W m${}^{-\mathrm{2}}$ for albedo, grain size and radiative forcing, respectively. When compared to two in situ measurements, found an albedo retrieval error of less than 0.005 and a radiative forcing error of about 2 $\pm$ 5 W m${}^{-\mathrm{2}}$. The grain size retrieval error for the smaller grains is assumed to be on the order of $\pm$25 $\mathrm{\mu }$m .

As mentioned in Sect. , investigations are needed to further quantify the grain size retrieval accuracy under melting snow conditions and concentrations of dust and BC. However, showed in a controlled field experiment that grain size retrievals leveraging the 1.03 $\mathrm{\mu }$m ice absorption feature were invariant from dry snow to melting snow conditions. Likewise, grain size retrievals from AVIRIS data across larger regions in the Sierra Nevada remained largely stable under change from dry to melting conditions . The sensitivity of grain size retrievals to the presence of dust and BC has been investigated . While the 1.03 $\mathrm{\mu }$m ice absorption feature clearly distorts when dust concentrations increase markedly , it appears to be only at concentrations far greater than those observed in normal field conditions where aeolian transport dominates . Further investigation is certainly needed. However, an additional line of evidence suggests that grain size does begin to plateau or rebound to smaller values after initial coarsening, as discussed in the next section.

While inherent snow albedo is controlled by grain size and impurity content over the course of a day, the solar zenith angle, cloud cover and topography also act to influence albedo and net solar radiation . For our acquisitions, we account for local illumination and viewing geometries in the albedo retrievals (Eqs.  and ). Additional error sources are related to terrain shadows and the digital elevation model (Sect. ). Up to 10 % of the snow pixels in the February data are in shadow due to the large solar zenith angle. Snow properties can be retrieved in shadows using the diffuse irradiance. However, we suspect that large uncertainties are introduced due to adjacency effects, which are difficult to model in snow-covered mountains. We therefore remove retrievals in shadow areas from our analysis, except for the snow cover calculation. Pixels very close to a ridge can be assigned an incorrect aspect value, which causes retrieval outliers. We found only very few snow pixels affected by this artifact and they do not show up in the histogram plots (Figs. , and –) due to a threshold of 15 occurrences. Changes in slope and aspect of the surface due to snow accumulation are not taken into account in this study and further analysis would be required to quantify associated errors.

To estimate the data quality and retrieval robustness, we performed an analysis of the variance of the retrieved quantities on timescales shorter than the natural variability. For example, the last two Sierra Nevada data sets are a single day apart and the snow products were retrieved with almost identical values. Between 7 and 8 May the snow-covered area has decreased by 2.8 km${}^{\mathrm{2}}$ or about 1 % of the mean. With a sample size of more than 280 000 snow pixels, the mean, percent difference (in brackets) and $\pm$1 standard deviation of the differences in albedo, grain size and radiative forcing are roughly $-$0.01 (0.5 %) $\pm$ 0.04, $+$34 (1.3 %) $\pm$ 112 $\mathrm{\mu }$m and $+$5 (0.8 %) $\pm$ 26 W m${}^{-\mathrm{2}}$ (see also Fig. ), respectively. These values are in the same order of magnitude as the expected retrieval accuracy . Nevertheless, the quantities and their signs show plausible values for changes in snow properties from one day to another at the end of the ablation period.

Inverse relationship of grain size with elevation of melting snow

As snow grain metamorphism is primarily driven by temperature and vapor pressure gradients in dry snow and then markedly reinforced by liquid water in wet snow, we expect increased snow grain sizes at lower elevations because of higher temperatures and a greater positive snowpack energy balance. This is in agreement with our observations shown in Figs.  and  of the initial ablation period in both study areas. However, we find the snow grain sizes at higher elevations increase significantly late in the ablation period, while they remain mostly unchanged at lower elevations. In the Sierra Nevada this is apparent in the May acquisitions relative to the April acquisition. In the Rocky Mountains this is apparent in the June acquisitions. Later in the ablation period, with increased solar irradiance, the snow grain size distribution becomes independent of the terrain aspect.

We suggest that there are two processes that account for these observations. First, changes in grain size are slowest when the snowpack is both cold and when it is fully isothermal, and most rapid during the transition between these two phases . When the grain growth is slow at lower elevations, but rapid at higher elevations, this indicates that lower elevations have become isothermal while the upper elevations are still undergoing that transition. Grain growth can be rapid in cold snowpacks, but this tends to occur near the ground where the temperature gradient between the ground and snow is greatest, and not at the surface. When liquid water is present, snow grains tend to form clusters but optically the absorbing path length seems to change very little. This same slowing of grain growth rates after the isothermal transition was observed in daily vertical profiles of optical snow grain size over the 2013 ablation period in the Senator Beck Basin .

Second, when there is an actual inversion of grain sizes with elevation, we suggest this is accounted for by rapid melt induced by high impurity content at the surface. Although the change in grain size is on the order of the uncertainty in grain size retrieval, it is not the absolute values that are important but rather the trend. A decrease in surface grain size in the presence of high impurity content was also observed in the daily field observations in 2013 by , when the dust layers coalesced at the surface and skies were clear both grain sizes and density in the surface layers decreased and surface roughness increased, consistent with observations in other regions . and as well as the results in this paper suggest that the intensification of snowmelt by dust in the visible wavelengths results in destructive metamorphism near the snow–atmosphere interface and snowmelt infiltration with refreezing and snow grain coarsening at depths of approximately 2 to 10 cm. The observed snow grains in the near-surface layer can therefore be smaller under intense snowmelt. The results could also suggest that sublimation could form frost flower-like crystals on the surface of a refreezing surface layer, resulting in an increase of SSA and a decrease in snow grain size .

In addition, our analyses confirm the inverse correlation between the snow albedo and grain size distribution . Smaller snow grains have smaller absorbing path lengths and smaller single-scattering asymmetry parameters (more side-scattering), which combined lead to higher snow albedo. We find in our results a similar inverse correlation between snow albedo and radiative forcing. The quantitative relationships between radiative forcing and the spatial and temporal parameters can vary with respect to snow albedo due to the differences in the local instantaneous solar irradiance.

As mentioned in the previous section, in situ albedo measurements suggest that the plateau in grain size with impurities is not just a result of the spectroscopy retrieval. In , it was found that near infrared (NIR)/shortwave infrared (SWIR) albedo measured at the detailed energy balance and radiation towers of the Senator Beck Basin would decrease as visible albedo dropped due to impurities, until the visible albedo decrease reached approximately 17 %. Once this threshold is reached, NIR/SWIR albedo would then start to increase again, consistent with a drop in grain size. Moreover, our observations in the field have shown that smaller grains are left when liquid water from intense melt percolates to deeper layers that are not detectable with the NIR/SWIR part of the spectrum. One hypothesis is that dust reflectance at higher concentrations may be raising the NIR/SWIR spectral albedo (Painter, personal observation). However, the concentrations necessary to drive visible albedos down 17–20 % are unlikely to be sufficient to cause a raise in NIR/SWIR albedo. Finally, it is possible that sublimation processes could be contributing or dominating the drop in the near-surface grain size, as discussed in . Given the negative feedback and the potential impact on our understanding of cryosphere energy balance, further study should investigate the controls on grain growth in the presence of impurities, with stratification according to impurity type, meteorological conditions and solar geometry.

Conclusions

Snow cover across the world's mountainous regions is important for both regional climate and hydrology. Snow albedo, itself controlled by impurity content and grain size, is particularly relevant for understanding melt initiation and melt rates. Ongoing consistent measurements of snow albedo only occur at six snow energy balance sites across the western USA, with even sparser or nonexistent measurements globally . The analysis of the snow remote sensing products presented here, derived from a series of imaging spectroscopy data, provides new insights into the spatial distribution of snow albedo and snow optical properties at five time steps across the accumulation/ablation transition in the Sierra Nevada and across the ablation period in the Rocky Mountains. The focus of most of the publications on snow property retrievals from imaging spectrometer data (see Sect. ) has been on algorithm development and technology demonstration. Here we take the next step, utilizing the proven data and algorithms to assess spatial and temporal patterns. While some of our results are intuitive, such as smaller grain sizes on north-facing slopes and the inverse relationship between grain size and albedo, other were counterintuitive, including faster grain growth rates at higher elevations and higher grain sizes at high elevations relative to lower elevations late in the ablation period.

The retrieval of impurity radiative forcing is especially relevant as efforts are now in place to reduce short-lived climate pollutants, such as black carbon, brown carbon, methane, tropospheric ozone and some hydrofluorocarbons by the Climate and Clean Air Coalition managed by the United Nations Environmental Programme . It will become critical to monitor these short-lived climate pollutants with respect to their spatial distribution on the local to global scale, as well as their change in time. Additional studies of the spatial and temporal distribution of impurity content in snow and ice of large mountain ranges, boreal and arctic regions are needed to enable improved understanding of snow in the framework of complex climate–terrain feedbacks. The application of our methods and algorithms to airborne and upcoming spaceborne imaging spectrometer data, such as the Hyperspectral Infrared Imager (HyspIRI), are potential tools to address this important question.

Felix C. Seidel wrote and applied the retrieval codes to derive the snow products, which are based on original algorithms by Thomas H. Painter. Karl Rittger prepared the AVIRIS and snow pillow data. S. McKenzie Skiles contributed content related to snow hydrology. Noah P. Molotch was the investigator for the Southern Sierra Snow Survey flights. Felix C. Seidel prepared the manuscript with contributions from all co-authors.

Acknowledgements

This research was performed at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Thanks to the AVIRIS and ER-2 teams as well as to Peter Kirchner for useful discussions. Edited by: F. Dominé