Monitoring of snow surface near-infrared bidirectional reflectance factors with added light-absorbing particles

Broadband snow albedo can range from 0.3 to 0.9 depending on microphysical properties and light-absorbing particle (LAP) concentrations. Beyond the widely observed direct and visibly apparent effect of darkening snow, it is still unclear how LAPs influence snow albedo feedbacks. To investigate LAPs' indirect effect on snow albedo feedbacks, we developed and calibrated the Near-Infrared Emitting and Reflectance-Monitoring Dome (NERD) and monitored bidirectional reflectance factors (BRFs) hourly after depositing dust and black carbon (BC) particles onto experimental snow surfaces. After comparing snow infrared BRFs to snow specific surface areas (SSAs), we found that both measured and modeled snow infrared BRFs are correlated with snow SSA. These results, however, demonstrate a considerable uncertainty of ±10 m² kg⁻¹ in the determination of snow SSA from our BRF measurements. The nondestructive technique for snow SSA retrieval that we present here can be further developed for science applications that require rapid in situ snow SSA measurements. After adding large amounts of dust and BC to snow, we found more rapid decreasing of snow BRFs and SSAs in snow with added LAPs compared to natural (clean) snow but only during clear-sky conditions. These results suggest that deposition of LAPs onto snow can accelerate snow metamorphism via a net positive snow grain-size feedback.

1 Introduction

Snow cover modulates the Earth's surface energy budget by reflecting a large portion of the incident solar radiative energy flux. As snow melts, more absorptive surfaces are uncovered resulting in a positive feedback mechanism known as the snow albedo feedback (SAF) (Qu and Hall, 2007; Hall, 2004). Light-absorbing particles (LAPs) within the snowpack, including elemental (black) carbon (BC), brown carbon, dust, and microbes, directly decrease snow albedo which can initiate the SAF and accelerate snowmelt (Bond et al., 2013; Qian et al., 2015; Skiles et al., 2018). Hadley and Kirchstetter (2012) experimentally verified that the reduction of snow albedo due to BC is enhanced for larger snow grains, implying a possible positive “grain-size” feedback induced by impurities in the snow. This positive feedback can also be inferred from the spectral snow albedo models presented by Warren and Wiscombe (1980) and from the Snow, Ice, and Aerosol Radiation (SNICAR) model (Flanner et al., 2007, 2009). Generally, two-layer models that represent snowpack as a collection of spheres can reproduce spectral hemispherical reflectances that compare well with observations (Grenfell et al., 1994).

Snow grain size is often quantified by its optically equivalent sphere effective radius, r_e, which is related to specific surface area, SSA, such that $\mathrm{SSA}=\mathrm{3}/({\mathit{\rho }}_{\mathrm{ice}}{r}_{\mathrm{e}}$), where ρ_ice is the density of pure ice (917 kg m⁻³ at 0 ^∘C). Snow SSA is defined as the total ice–air interfacial surface area S to ice mass m ratio such that

$$\begin{array}{}\text{(1)}& \mathrm{SSA}=S/m={\displaystyle \frac{S}{{\mathit{\rho }}_{\mathrm{ice}}V}}\end{array}$$

is expressed in terms of its total ice mass or volume V (Hagenmuller et al., 2016; Gallet et al., 2014). In snowpack with large temperature gradients, the diffusion of vapor causes snow SSA to decrease during the natural process of snow metamorphism (Flanner and Zender, 2006; Wang and Baker, 2014). In the absence of a temperature gradient, an isothermal snowpack with snow grains with low radii of curvature undergoes coarsening in a process driven by the Kelvin effect. Ebner et al. (2015) show that measurements of snow SSA evolution in isothermal snow agree with the isothermal snow metamorphism modeling framework developed by Legagneux et al. (2004) and Legagneux and Domine (2005). These laboratory studies express snow SSA in isothermal metamorphism as a function of time t as follows:

$$\begin{array}{}\text{(2)}& \mathrm{SSA}\left(t\right)={\mathrm{SSA}}_{\mathrm{0}}{\left({\displaystyle \frac{\mathit{\tau }}{\mathit{\tau }+t}}\right)}^{\mathrm{1}/n}\end{array}$$

for initial snow SSA₀ at t=0, growth rate parameter τ, and growth exponent n.

Snow SSA strongly affects absorption of infrared radiation. This relationship is evident from measurements of infrared reflectance that are highly correlated with snow SSA for various snow types (Domine et al., 2006). Among others, Gallet et al. (2009, 2014) use this correlation to determine dry snow SSA and wet snow SSA, respectively, using 1.31 µm directional hemispherical reflectance measurements (1.55 µm for measurements of snow SSA >60 m² kg⁻¹). Arnaud et al. (2011) present the Profiler Of Snow SSA Using SWIR reflectance Measurement (POSSUM), which applies the theoretical formulations from Kokhanovsky and Zege (2004) and Picard et al. (2009) to derive snow SSA from hemispherically averaged bidirectional reflectance measurements. Other studies establish techniques to obtain snow SSA using methane gas absorption (Legagneux et al., 2002) and X-ray micro-computed tomography (X-CT) in cold rooms (Pinzer and Schneebeli, 2009). Matzl and Schneebeli (2006) also derive snow SSA using infrared photography. Other techniques that are nondestructive enable the rapid retrieval of snow effective radius r_e from field measurements. Gergely et al. (2014), for example, demonstrate a technique to quickly determine the snow optically equivalent diameter from 0.95 µm bi-hemispherical reflectance measurements. Painter et al. (2007) infer snow r_e from spectral hemispherical directional reflectance factor measurements using a contact probe and a spectrometer.

While previous studies monitor snow metamorphism in clean snow, in situ experimentation of how LAPs affect the snow grain-size feedback diurnally in natural environments is challenging. These challenges can be partly attributed to the limited availability of inexpensive snow SSA measurement devices that can operate quickly in the field, and to the numerous constraints on the reproduction of naturally occurring conditions for which LAPs in snow strongly influence the climate system. While the POSSUM is a suitable instrument for this study, here, our first objective is to demonstrate the utility of a new instrument we use to quickly obtain approximate surface snow SSA (with LAPs). Our second objective is to investigate the effects of added LAPs on the snow grain-size feedback.

In the following sections, we first describe a numerical model that simulates multiple scattering in three dimensions to aid the design and calibration of the Near-Infrared Emitting and Reflectance-Monitoring Dome (NERD). Next, we describe the NERD apparatus, which enables multiple 1.30 and 1.55 µm bidirectional reflectance factor (BRF) measurements in just minutes while minimally altering the snow structure. To calibrate the NERD with respect to snow SSA, we then compare snow BRFs with X-CT-derived SSAs and develop an empirical relationship between measured SSAs and snow 1.30 µm BRFs. Finally, we discuss results from our LAPs in snow experiments in which we monitored accelerated snow metamorphism. Overall, this study demonstrates the NERD measurement technique and conditions for which the snow grain-size feedback can be enhanced by the presence of LAPs.

2 Instrumentation and methods

2.1 Monte Carlo modeling of multiple scattering in snowpack

The Monte Carlo method is applied to numerically simulate three-dimensional (3-D) light scattering within an idealized snowpack. Gaussian distributions (with 0.085 and 0.130 µm full widths at half-maximum) of photon wavelengths (centered around 1.30 and 1.55 µm) were selected to model light emission by narrow-band infrared light-emitting diodes (LEDs). One by one, photons are initiated downward into the snow medium, as demonstrated by Kaempfer et al. (2007), and propagated in optical depth space. Extinction, absorption, and polarization are accounted for following the scattering approach for geometrical optics described by Malinka (2014). Accordingly, random numbers (RN_i) in the interval (0, 1) are generated to determine the photon optical path lengths, l_i such that

$$\begin{array}{}\text{(3)}& l_i=\ln (\mathrm{1}/{\mathrm{RN}}_i)\end{array}$$

before the first scattering event (i=0) and again after each scattering event i. Absorption (and termination) of a photon can occur during a scattering event if an additionally generated random number (between 0 and 1) is greater than the particles' single scattering albedo. When a photon is scattered, its new direction cosines are determined from an optimized “rejection method” using the particles' scattering matrices (Ramella-Roman et al., 2005).

To calculate snow BRFs, 1 000 000 photons are propagated and traced through the modeled snowpack until they are absorbed or exit the medium. The snowpack is modeled as two-phase (air and ice) media containing a regular arrangement of suspended homogenous ice particles. Aspherical particle single scattering properties, including the mass absorption cross sections, asymmetry parameters, single scattering albedos, projected areas, volumes, and scattering matrices were calculated by Yang et al. (2013) for randomly oriented droxtals and solid hexagonal columns. For spheres, Mie theory is applied but with the analytical Henyey–Greenstein phase function (van de Hulst, 1968) to improve computational efficiency. The subset of shape habits (smooth droxtals and hexagonal columns) from the large dataset provided by Yang et al. (2013) was selected because these shape habits are purely convex solid ice. Because they are solid convex bodies, their SSAs can be computed from the their projected areas and volumes (Vouk, 1948). Azimuthal mean BRFs are calculated according to the reflectance definitions presented by Dumont et al. (2010), Hudson et al. (2006), and Nicodemus et al. (1977). Accordingly, photon exit angles are grouped into 30 exit zenith angle (θ_r) bins at 3^∘ resolution. Azimuthal (ϕ) mean BRFs are calculated by zenith angle, θ_r, from the total incident photon flux, Φ_i (at a given θ_i) from

$$\begin{array}{}\text{(4)}& R({\mathit{\theta }}_{\mathrm{i}};{\mathit{\theta }}_{\mathrm{r}})=(\mathrm{2}{\mathrm{\Phi }}_{\mathrm{i}}\sin {\mathit{\theta }}_{\mathrm{r}}\cos {\mathit{\theta }}_{\mathrm{r}}{)}^{-\mathrm{1}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathrm{d}{\mathrm{\Phi }}_{\mathrm{r}}\mathrm{d}\mathit{\varphi },\end{array}$$

where dΦ_r represents the reflected photon flux through discrete solid angle bins. In the denominator, the cos θ_r factor satisfies Lambert's cosine law while sin θ_r accounts for the zenith angular dependence of the azimuthally integrated projected solid angle. Finally, the factor 2 is necessary to normalize the resulting weighting function w(θ_r)=sin θ_rcos θ_r as

$$\begin{array}{}\text{(5)}& \underset{\mathrm{0}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}\sin {\mathit{\theta }}_{\mathrm{r}}\cos {\mathit{\theta }}_{\mathrm{r}}\mathrm{d}\mathit{\theta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

Statistical uncertainty was determined by computing BRFs from simulations of Lambertian surfaces and tested using different numbers of photons. Azimuthal averaging reduces the BRFs' dimensionalities so that fewer photons are needed to mitigate Monte Carlo noise. Equation (4) was initially applied to Monte Carlo simulations of 75 000 photons reflected by various ideal Lambertian surfaces. At 3^∘ (zenith angular) resolution, 30 and 60^∘ BRFs of Lambertian surfaces were simulated accurately to within ±2 %. This margin of uncertainty was determined by computing rms differences between calculated and specified Lambertian reflectances ranging from 0 to 1. Across this range, rms differences at 30 and 60^∘ were generally less than 0.01. In subsequent test cases, simulating snowpack BRFs with up to 10 000 000 photons did not significantly change results when compared with simulations of 250 000 photons. Ultimately, it was determined that simulations with 1 000 000 photons were appropriate for the remainder of this study.

2.2 Near-infrared bidirectional reflectance factor measurements

The NERD is designed to measure 1.30 and 1.55 µm BRFs. These wavelengths are selected for snow SSA retrieval due to the strong dependence of snow albedo on grain size in these parts of the spectrum. Snow spectral albedo is simulated here using the SNICAR model to demonstrate its sensitivity (Fig. 1). While snow spectral albedo is sensitive to the snow effective radius r_e (and thus snow SSA), it is not sensitive to small black carbon concentrations at these wavelengths. Snow SSA can therefore be retrieved using 1.30 µm and 1.55 µm reflectance measurements for snow with small black carbon concentrations.

Figure 1 Black sky spectral snow albedo under nadir illumination. Snow albedo is simulated by the Snow, Ice, and Aerosol Radiation (SNICAR) model (Flanner et al., 2007). Dashed-dotted curves represent clean snow of medium-high SSA (60 m² kg⁻¹, gray) and medium-low SSA (20 m² kg⁻¹, black) to show the dependence of snow albedo on snow SSA. Dotted curves represent contaminated snow with uncoated black carbon (BC) particulate concentrations of 100 ng g⁻¹.

The design principle is similar to the DUal Frequency Integrating Sphere for Snow SSA measurements (DUFISSS) (Gallet et al., 2009). The NERD also uses 1.30 (1.31 in DUFISSS) and 1.55 µm emitters to illuminate the snow surface from nadir (15^∘ off nadir for 1.55 µm in NERD). The main distinction between the DUFISSS and the NERD is the type of reflectance measured. Gallet et al. (2009) use an integrating sphere to measure hemispherical reflectance. In the NERD, however, photodiodes are directed toward the illuminated surface in a black dome to measure BRFs, as demonstrated by Arnaud et al. (2011). The interior of the dome is painted with a flat black paint to increase absorptivity and to minimize internal reflections between the dome and snow surface. To test the near-infrared absorptivity of the black pigment, the NERD was used to measure BRFs of a painted black surface. All measured BRFs of the black surface were less than 0.03.

Figure 2 Near-Infrared Emitting and Reflectance-Monitoring Dome (NERD) schematic (a) and photographs (b). The photograph on the bottom right shows the (radial) placement of photodiodes and LEDs (toward the center) within the dome.

Four infrared light-emitting diodes (LEDs) are mounted into the 17 cm in diameter black styrene half-sphere. Two LEDs with peak emission wavelengths of 1.30 µm are mounted at nadir and 10^∘ relative to zenith and two LEDs with peak emission wavelengths of 1.55 µm are mounted at 15^∘ off nadir (see Fig. 2). The 1.30 µm LEDs have spectral line half-widths of 0.085 µm and half-intensity beam angles of 10^∘, while the 1.55 µm LEDs have half-maximum bandwidths of 0.130 µm and 20^∘ beam angles. These high powered, narrow beam infrared LEDs were selected to illuminate a small oval (estimated major axes of 1.5 cm at 1.30 µm and 3.0 cm at 1.55 µm) of the experimental surface to maximize the reflected radiance signal. The reflected radiance signal is measured using four indium gallium arsenide photodiodes mounted in two different azimuthal planes (0 and 90^∘ relative to the illumination): two each at 30 and 60^∘ relative to zenith. Photodiodes highly sensitive to light with wavelengths ranging from 0.80 to 1.75 µm and relatively large active areas (1 mm) were selected to maximize sensitivity.

Because the orientation of LEDs and photodiodes are fixed, reflectance factors of surfaces with negligible subsurface scattering can be obtained after calibration using two diffuse reflectance targets (Gallet et al., 2009; Gergely et al., 2014; Dumont et al., 2010). These Lambertian targets reflect incident light according to Lambert's cosine law and appear equally bright at all viewing angles. By comparing the measured voltage signal from the experimental (snow) surface to that measured from the reflectance targets, two BRFs at both 30 and 60^∘ viewing angles are obtained for each light source. While subsurface scattering of visible light in snow is pervasive (Smith et al., 2018), the light penetration in snow near 1.30 and 1.55 µm is at most a couple centimeters due to the strong absorption features in the near-infrared (Kaempfer et al., 2007; Grenfell et al., 1994; Brandt and Warren, 1993). Subsurface scattering is therefore assumed to be minimal and fully contained within each photodiode's field of view.

Table 1 NERD Lambertian reflectance measurements. Tabulated values represent median BRFs (R(θ_i;θ_r)) calculated for n samples of measurements from Lambertian reflectance targets with nominal reflectances of ρ_L (rms difference in parentheses). Linear regressions for each wavelength (λ) are calculated from $n_{\mathrm{1}}+n_{\mathrm{2}}=N$ samples.

To validate NERD reflectance measurement accuracy, precision, and responsiveness, measured BRFs of reflectance standards are recorded after calibration. BRFs (R) for each LED and photodiode viewing zenith angle (θ_i;θ_r) combination are measured in temperatures ranging from −20 to +2 ^∘C. In general, NERD BRFs of the Lambertian reflectance standards are accurate to within ±2 %. We quantify instrument precision (2 %) by computing root mean squared (rms) differences from repeated measurements (see Table 1). Linear regressions quantify the linear response (A) over the reflectance range of 0.41 to 0.95. Response uncertainty ranges from −2 % to +3 % and from +1 % to +3 % at 1.30 and 1.55 µm, respectively. These test results indicate the NERD's ability to obtain BRFs on smooth reflectance standards with a measurement uncertainty of 1 %–2 %.

2.3 Snow specific surface area measurements

Surface snow (just the top few centimeters) samples were collected in nature over the span of 3 years (winters 2015–2017) and transported in coolers to the nearby U.S. Army's Cold Regions Research Engineering Laboratory (CRREL) in Hanover, New Hampshire. Depth hoar samples, however, were instead grown inside the CRREL at −20 ^∘C using a forced temperature gradient. Snow samples are classified based on X-CT results according to Fierz et al. (2009) (Table 2). Snow SSA was derived from X-CT analysis and contact spectroscopy as described below.

Table 2 Snow sample physical classification (Fierz et al., 2009). Snow densities, porosities, and specific surface areas (SSAs) are calculated from X-ray micro-computed tomography.

To determine snow SSA, X-CT was conducted on a class of six snow samples according to Lieb-Lappen et al. (2017). X-ray (40–45 kV, 177–200  µA) transmission through cylindrical snow samples was measured at rotation steps of 0.3–0.4^∘. To limit scan times to 15 min, exposure time was set to 340 ms at a cubic voxel resolution of 14.9 µm. Processing software enables SSA calculations from 3-D morphology results (Pinzer and Schneebeli, 2009).

In some cases, snow samples were scanned several hours or days after snow BRFs were measured. To correct for natural isothermal snow SSA metamorphism while samples were being stored, Eq. (2) was applied with time t set equal to the total time elapsed between NERD measurements and X-CT scan times. The unitless growth rate parameters τ=721.2 or 14 400 and growth exponents n=2.15 or 0.32, respectively, were inferred from Ebner et al. (2015) depending on the snow sample type. Applying this correction yielded a SSA decrease between 3 and 5 m² kg⁻¹.

Snow SSA was also inferred from optical grain-size measurements using contact spectroscopy (Painter et al., 2007). Snow reflectance spectra were collected using an ASD FieldSpec4 spectroradiometer and a high-intensity contact probe with reference to a Spectralon white reference panel. Snow effective radius r_e is determined from the normalized area of the absorption feature centered at approximately 1 µm using a lookup table (Nolin and Dozier, 2000). These measurements were conducted inside the CRREL in a cold lab only for the depth hoar (DH_2016) and rounded grain (RG_2015) samples.

2.4 Light-absorbing particles in snow experiments

Snow BRFs and SSAs were measured throughout the day in the following dust and BC in snow experiments. Sand particles and hydrophobic BC were sifted multiple times with a salt shaker (with holes of roughly 1 mm in diameter) to filter out larger particles. The filtered LAPs were then deposited onto experimental snow plots in an open field in Hanover, New Hampshire, on 10 and 16 February 2017 shortly after fresh snowfall events. For each experiment, 1 m² plots of snow were designated as natural (control) or contaminated (experimental). Snow BRFs and SSAs were obtained using the NERD and from X-CT analyses, respectively. For each set of NERD measurements, 30 and 60^∘ BRFs were both recorded four times. BRFs were measured over two different locations within the experimental plot using two photodiodes at each viewing angle (30 and 60^∘).

2.4.1 Cloudy sky with diffuse ambient lighting (10 February 2017)

Early on 10 February, experimental plots were loaded with BC until visible darkening was apparent. Snow BRFs were measured shortly after 00:00, 03:00, 06:00 eastern standard time (EST) (during the night) and then periodically throughout the day. Because these plots were well shaded by tall trees, these measurements were used to monitor snow metamorphism without the influence of direct solar illumination. Furthermore, mostly cloudy conditions on 10 February obscured direct solar radiation so that ambient lighting was diffuse.

2.4.2 Clear sky with direct solar irradiance (17 February 2017)

On 17  February, just a pinch (<1 g) of BC and 30 g of sand were deposited (quasi-)uniformly over separate 1 m² plots. These experimental deposition fluxes were selected to mimic the most extreme LAPs' loading event observed by Skiles and Painter (2017) in the San Juan Mountains in Colorado. As in the previous experiment, snow BRFs were measured periodically throughout the day; however, all snow plots were in full view of the clear sky to maximize incident direct solar irradiance.

3 Results and discussion

3.1 Near-infrared reflectance calculations for regular homogenous ice particles

To validate the Monte Carlo model for snow applications, we simulated 1.30 and 1.55 µm narrow band black-sky albedo for snow effective radius r_e ranging from 36 to 327 µm (SSA =80 to 10 m² kg⁻¹) and compared results with the SNICAR model (Fig. 3). As expected, results from Monte Carlo simulations yield slightly higher hemispheric reflectances for droxtals (for all SSAs) and solid hexagonal columns (for SSA >40 m² kg⁻¹) than for spheres and from the SNICAR model (Dang et al., 2016). For all particle sizes, Monte Carlo simulations using spherical ice particles generated nearly identical narrow-band albedo values compared to those from Flanner et al. (2007). These results are also similar to the numerical modeling results from Picard et al. (2009) and the analytical solutions from Kokhanovsky and Zege (2004), though for different particle shapes.

Figure 3 Modeled 1.30 µm nadir (a) and 1.55 µm 15^∘ directional hemispherical reflectance (b) for various snow SSA. Solid line segments connect albedo calculations from Monte Carlo simulations of light scattering in snow media comprised of droxtals (stars) and solid hexagonal columns (diamonds). Circles connected by dashed or dotted line segments connect snow albedo calculations modeled as spherical ice particles from Monte Carlo modeling (black) and from the Snow, Ice, and Aerosol Radiation (SNICAR) online model (gray).

To guide our choice of a snow BRF to SSA calibration function, we calculated BRFs from Eq. (4) for various particle SSAs ranging from 10 to 80 m² kg⁻¹ (Fig. 4, line graphs). Generally, we found exponential (linear) relationships between 1.30 (1.55) µm BRFs and particle SSAs. Modeled 1.30 µm BRFs are slightly higher at 30^∘ than at 60^∘ for particle SSA >30 m² kg⁻¹. For a given particle size, a different shape habit can yield changes in BRFs of as much as 0.1. Monte Carlo modeling results yield the highest reflectances for droxtals and the lowest reflectances for spheres.

Figure 4 NERD SSA calibration. All panels contain NERD-measured and Monte Carlo-modeled (black) BRFs (1.30 µm, a, b; 1.55 µm, c, d; 30^∘ viewing, a, c; 60^∘ viewing, b, d scattered against snow SSA. NERD BRFs are scattered against X-CT-derived snow SSA. Line segments connect Monte Carlo modeled BRFs of snow media comprised of spheres (filled circles, dashed lines), droxtals (stars, solid lines), and solid hexagonal columns (diamonds, solid lines). Snow sample key codes, symbols, and colors conform with the physical snow classification standards defined by Fierz et al. (2009) (Table 2). Vertical error bars on NERD BRFs represent standard deviations calculated from multiple azimuthal samples. Horizontal error bars on X-CT-derived SSA, where present, represent standard deviations from multiple scans on similar snow samples.

3.2 Relating snow specific surface area to near-infrared bidirectional reflectance factors

To calibrate the NERD for snow SSA retrieval, we compared X-CT-derived snow SSA with NERD snow BRF measurements (Fig. 4). In general, measured snow BRFs are directly related to snow SSAs. At 1.30 µm, BRFs range from just under 0.2 (for low SSAs) to as high as 0.7 (for high SSAs) and are slightly higher at 60^∘ than at 30^∘. We observed 1.55 µm snow BRFs close to 0 (for low SSAs) and as high as 0.2 (for high SSAs). We recorded the highest 1.55 µm snow BRFs at 60^∘ for fresh snow (needles). At 1.30 µm, measured 30^∘ snow BRFs for varying snow SSA fall within the envelope of Monte Carlo modeled BRFs for all three shape habits (Fig. 4, top left). Modeling results are in closest agreement with measurements at 30^∘ viewing for 1.30 µm. At 1.55 µm, measured BRFs are larger than predicted from modeling across all SSA.

Our measurements and modeling results show a considerable spread of BRFs across snow samples and shape habits, respectively, for both 1.30 and 1.55 µm and for both 30 and 60^∘ viewing angles. The spread in measurements, in particular, indicates a considerable uncertainty in the ability to retrieve snow SSA from measured BRFs. While the 1.30 µm 30^∘ viewing zenith angle BRF combination most closely agrees with modeled BRFs, the 60^∘ viewing zenith angle can provide a second estimate of snow SSA. Reporting two snow SSA values using both view angles can provide an estimate of the variability in SSA retrieval resulting from the anisotropy of the snow bidirectional reflectance distribution function in the near-infrared (Dumont et al., 2010).

Our finding of the exponential relationships between snow SSAs and 1.30 µm BRFs is consistent with results from previous studies (Picard et al., 2009; Gallet et al., 2009). Gallet et al. (2009) also identify a linear relationship between 1.55  µm reflectance and snow SSA and use the longer wavelength in their DUFISSS to obtain measurements of high snow SSA (>60 m² kg⁻¹). In this study, however, nearly all snow samples are lower than this threshold.

Measured snow BRFs at 1.55 µm are higher than both hemispherical reflectance measurements by Gallet et al. (2009) and Monte Carlo calculations. This discrepancy might be explained by the NERD LEDs, which have full widths at half-maximum (FWHM) of 0.130 µm and emit non-negligible light at shorter wavelengths. Supplementary SNICAR and albedo calculations using a broadened FWHM of 0.260 µm support this hypothesis in which values closer to the measured 1.55 µm BRFs are reproduced (Fig. 3b).

In light of these empirical and numerical results, we propose the following general exponential form relating 1.30 µm snow BRFs to SSA such that

$$\begin{array}{}\text{(6)}& \mathrm{SSA}=\mathit{\alpha }\exp \left(R_{\mathrm{1.30}}\right)+\mathit{\beta }\end{array}$$

for 1.30 µm snow BRF R_1.30. Using least squares regression analysis, we estimate “best fit” values for α and β for both 30 and 60^∘ viewing zenith angles (see Fig. 5).

Figure 5 Snow 1.3 µm BRF to SSA exponential regressions. Markers with standard error bars represent SSA measured with X-CT^∗ and scattered against nadir 1.30 µm 30^∘ (a) and 60^∘ (b) BRFs measured by the NERD. Solid curves represent exponential regression functions, and dashed curves represent standard errors of the regressions such that at 30^∘ viewing, $\mathit{\alpha }=\mathrm{88.7}(\pm \mathrm{9.50})$ m² kg⁻¹ and $\mathit{\beta }=-\mathrm{103}$ m² kg⁻¹, and at 60^∘ viewing, $\mathit{\alpha }=\mathrm{91.7}(\pm \mathrm{10.13})$ m² kg⁻¹ and $\mathit{\beta }=-\mathrm{113}$ m² kg⁻¹. ^∗ Hollow triangles (blue, depth hoar; pink, rounded grains) represent snow SSA measurements derived from contact spectroscopy.

Ideally, an empirically derived calibration function would include SSA measurements from multiple methods to mitigate uncertainties associated with collection methods needed for X-CT analysis. Such collection methods can easily change the snow microphysical characteristics and lead to biases in the X-CT-derived SSAs. As a preliminary validation of Eq. (6), we compare snow SSA results to SSAs derived from snow effective radius r_e measurements conducted using contact spectroscopy in Fig. 5. Encouragingly, two out of three measurements fall within the bounds of the standard error of the regression. Contact spectroscopy snow SSA values are consistently higher than those calculated from X-CT analysis and are also higher than those determined by Eq. (6). This comparison provides a preliminary assessment of the NERD snow BRF to SSA calibration function.

Hereafter, we apply Eq. (6) in the following LAPs in snow experiments to estimate hourly snow SSAs from measured snow BRFs. Because the remainder of this study is concerned with relatively large changes in SSAs, approximate SSA retrieval using the NERD is useful to quantitatively assess how added LAPs affect snow grain-size feedback and snow metamorphism.

3.3 Light-absorbing particles' effect on snow metamorphism

First, to monitor snow metamorphism without solar heating, during the early morning (night) hours on 10 February, we deposited BC onto an experimental plot after the previous day's snowfall. Surface temperatures ranged from −14 to −9 ^∘C. We observed low to moderate wind speeds from the early morning hours through the afternoon with partly to mostly cloudy conditions during the day. The 1.30 µm snow BRF measurements conducted shortly after BC application indicate a direct darkening of 0.03 at 30^∘ viewing and 0.05 at 60^∘ viewing. While these decreased BRF values are considerable, this experiment contained heavily contaminated snow with a high BC concentration unlikely to occur in nature. Despite the high BC concentration, BRFs measured at 1.30 (1.55) µm remained within 0.5 and 0.6 (0.1 and 0.2) throughout the day in both contaminated and natural snow (see Fig. 6). X-CT analysis showed small differences in morning (49 m² kg⁻¹) and afternoon (48 m² kg⁻¹) snow SSA. Our results indicate that heavy BC loading had little to no effect on snow metamorphism during the 16 h experiment without direct solar irradiance.

Figure 6 The 10 February control results (mostly cloudy). All panels contain NERD measurement data collected from 00:00 through to 17:00 eastern standard time (EST) on 10 February, 2017. Filled circles connected by solid lines represent measurements on natural snow. Filled squares connected by dashed lines represent measurements on snow heavily contaminated by hydrophobic BC. In panels (a) and (b), blue (green) curves represent 1.30 µm (1.55 µm) BRFs at 30^∘ viewing (a) and 60^∘ viewing (b). Error bars represent standard errors calculated from sample averages from as many as eight locations within each 1 m² plot. In panels (c) and (d), blue curves represent NERD-calibrated SSAs from 30^∘ (c) and 60^∘ (d) viewing BRFs. Red curves represent SSAs derived from X-CT scans. Error bars represent NERD calibration uncertainty computed from regression analysis (standard errors of the gradients).

Figure 7 The 17 February experimental results (clear sky). All panels contain NERD measurement data collected from 07:00 through to 17:00 EST on 17 February 2017. Filled circles connected by solid lines represent measurements on natural snow. Filled squares connected by dashed lines represent measurements on snow lightly contaminated (<1 g m⁻²) by hydrophobic BC. Triangles connected by dotted lines represent measurements on snow contaminated by dust (30 g m⁻²). In panels (a) and (b), blue (green) curves represent 1.30 µm (1.55 µm) BRFs at 30^∘ viewing (a) and 60^∘ viewing (b). Error bars represent standard errors calculated from sample averages from as many as eight locations within each 1 m² plot. In panels (c) and (d), blue curves represent NERD-calibrated SSAs from 30^∘ (c) and 60^∘ (d) viewing BRFs. Red curves represent SSAs derived from X-CT scans. Error bars represent NERD calibration uncertainty computed from regression analysis (standard errors of the gradients).

Second, to monitor snow metamorphism occurring after forced large BC and dust deposition events under direct solar illumination, on 17 February, we set up a similar experiment in full view of the sun. Surface temperatures ranged from −4 to +2 ^∘C. We observed minimal wind speeds and cloud cover resulting in calm, clear-sky conditions. BRF measurements conducted shortly after LAPs application indicate minimal direct darkening at 1.30 µm and possibly a small brightening effect by dust at 1.55 µm. In natural snow, 1.30 µm BRFs remained close to 0.5 throughout the day, with the lowest values (0.49) recorded in the afternoon (13:00 EST) and the highest values (0.55) recorded in the morning (08:00 EST) and evening (17:00 EST). The 1.55 µm BRFs remained just above 0.1. In the dust loaded plot, snow 1.30 (1.55) µm BRFs decreased rapidly from above 0.5 (0.1) before 10:00 to below 0.3 (0.05) by 13:00 EST. We found less extreme metamorphism in the lightly contaminated snow with added BC, as BRF measurements decreased from above 0.5 (0.1) to below 0.45 (0.1), 1.30 µm snow BRFs slightly increased thereafter (from 13:00 to 17:00 EST) in both natural and contaminated snow (see Fig. 7). Snow SSA also decreased throughout the day. From X-CT analysis, we found morning snow SSA to be about 50 m² kg⁻¹, which thereafter decreased to 41, 23, and 18 m² kg⁻¹ in natural, BC loaded, and dust loaded snow, respectively. NERD-derived snow SSA appears to be biased low in the afternoon dust loaded plot. This bias might be an indication of the presence of liquid water that was also visible to the naked eye. X-CT scans performed on this snow sample are representative of refrozen snow and do not conform to the isothermal snow SSA correction (Eq. 2) applied to snow samples scanned several hours after collection. In BC loaded plots, we observed a large spatial heterogeneity in measurements, indicating that small BC deposition has a powerful localized effect on snow metamorphism.

These results suggest that extreme LAPs deposition can accelerate snow metamorphism. The primary cause of this accelerated process is enhanced absorption of solar radiation by LAPs. Surprisingly, added BC had little to no effect on snow metamorphism during cloudy conditions. In the clear-sky experiments, LAPs enhanced solar absorption at the surface which warmed the snowpack. As the snow surface began to melt, near-infrared reflectance decreased rapidly. Rapidly decreasing near-infrared reflectance is indicative of either the accumulation of liquid water from melting snow or decreasing snow SSA.

Accelerated snow metamorphism by dust loading is consistent with the findings of Skiles and Painter (2017). The indirect effect of LAPs on snow is also demonstrated by Hadley and Kirchstetter (2012) in which the albedo reduction due to the presence of BC in snow is amplified in snow of lower SSA. This enhancement of snow albedo reduction is another source of instability in the snowpack that increases the strength of the snow grain-size feedback. Typical BC deposition events are very small, however, so it is difficult to reproduce natural BC concentrations when adding any BC to a 1 m² plot.

4 Conclusions

LAPs in snow accelerated snow metamorphism and enhanced positive grain-size feedback during cloud free, calm weather conditions when surface air temperatures were near 0 ^∘C. To observe this effect, we engineered an instrument (i.e., the NERD) that measures 1.30 and 1.55 µm snow BRFs without destroying snow samples. We evaluated the NERD's accuracy, precision, and responsiveness by testing with idealized Lambertian surfaces before obtaining snow BRFs. Notwithstanding the limitations associated with retrieving precise snow SSA from BRFs, we proposed an analytical calibration function relating snow SSA to 1.30 µm BRFs. Our results lead to the conclusion that the NERD can provide estimates of snow SSA to within ±10 m² kg⁻¹.

The NERD will serve to further study the effects of LAPs on snow metamorphism and to explore the spatial heterogeneity of snow SSA. Because the NERD can also operate quickly, NERD measurements can complement satellite borne observations during narrow sampling windows. To fulfill these pursuits, however, a more comprehensive snow SSA measurement uncertainty assessment is needed to fully validate the NERD technique. Validation would also include additional independent measurement methods on snow samples with a larger snow SSA span from a variety of environmental conditions and further experimentation into the small-scale effects on snow BRF measurements. Further investigation of the natural variability of snow near the surface and the related quantitative uncertainties regarding the NERD snow SSA retrieval method will be the subject of a follow-up study. Such a study will also include snow of higher SSA to determine the utility of 1.55 µm snow BRFs in measuring fresh snow of extremely high SSA as demonstrated by Gallet et al. (2009). This would help expand the utility of the NERD measurement technique for future Arctic and Antarctic campaigns, building on the polar studies conducted by Legagneux et al. (2002) and Libois et al. (2015).

Appendix A: NERD photodiode current amplifiers

To detect reflected radiance signals, photodiodes are reverse biased to induce currents linearly related to the amount of light incident on their active regions. Because these light signals are reflected from the experimental surface, the currents induced by the photodiodes are very small (nano- to microampere). To measure the small currents, the photodiodes are connected to transimpedance amplifiers. The transimpedance amplifier circuits convert and amplify the small photodiode currents into measurable voltage signals.

Two NERDs are engineered with different photodiode current amplifications. Photodiode current amplification is determined by the feedback resistance in the transimpedance amplifier circuits. Active low pass filters are applied between the amplifier and the analog-to-digital converter (ADC) to reduce noise. This filter is designed to have a time constant of less than 0.5 s to achieve balance between adequate noise reduction and speed.