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 m2 kg−1 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.


Introduction
Common light absorbing impurities (LAI) in snow include elemental (black) carbon (BC), brown carbon, and dust particulate matter, all of which play an important role in the climate system (Bond et al., 2013;Qian et al., 2015). In particular, these LAI 15 at the snow surface can reduce snow albedo and enhance snow melt via the snow albedo feedback (Qu and Hall, 2007;Skiles and Painter, 2017). Hadley and Kirchstetter (2012) experimentally verify this effect and confirm what Warren and Wiscombe (1980) and the widely used Snow, Ice, and Aerosol Radiation (SNICAR) model (Flanner et al., 2007(Flanner et al., , 2009 predict in that snow albedo reduction by BC is enhanced for larger snow effective radii. Numerical snow albedo studies, including those from Picard et al. (2009), typically model snowpack as a semi-infinite medium of suspended 20 spherical ice particles. These models are highly accurate for spectral snow albedo calculations when the snow effective radius is a tunable parameter. Spherical snow grain size is related to specific surface area (SSA) by its effective radius, such that SSA = 3/(ρ ice Re). Expressions of sphere effective radii, r ef f and Re, defined by ice particle surface area S versus ice particle projected area A, respectively, are equivalent for convex bodies (Vouk, 1948). Snow SSA is defined as the total ice-air interfacial surface area S to ice mass m ratio, expressed in terms of its total ice volume V such that, where ρ ice is the density of pure ice (917 kg/m 3 at 0 • C) (Hagenmuller et al., 2016;Gallet et al., 2014). Snow SSA strongly affects absorption of infrared radiation. This relationship is evident from measurements of infrared reflectance that are highly 5 correlated with snow SSA for various snow types (Domine et al., 2006). Among others, Gallet et al. (2009) and Gallet et al. (2014) exploit this correlation in the accurate determination of 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 2 kg −1 ). Other studies establish techniques to accurately obtain snow SSA using methane gas absorption (Legagneux et al., 2002) and X-ray microcomputed tomography (X-CT) in cold rooms (Pinzer and Schneebeli, 2009). Matzl and Schneebeli (2006) also derive snow 10 SSA using infrared photography. Other techniques that are nondestructive enable the rapid retrieval of snow optically equivalent grain size from field measurements. Gergely et al. (2014), for example, demonstrate an accurate technique to quickly determine the snow optically equivalent diameter from 0.95 µm bi-hemispherical reflectance measurements. Painter et al. (2007) infer snow optical grain radius (r ef f ) from spectral hemispherical directional reflectance factor measurements using a contact probe and a spectrometer. 15 In snowpacks with high 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 isothermal snow, highly faceted snow grains with relatively high SSA and low radii of curvature undergo 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 studies express snow 20 SSA in isothermal metamorphism as a function of time t as follows, for initial snow SSA 0 at t = 0 and adjustable parameters τ and n. Domine et al. (2009), however, observe increasing snow SSA due to the fragmentation of surface snow grains mobilized by wind.
Currently, it is still unclear how solar heating of snow by LAI affects possible negative feedbacks relating to temperature 25 gradient metamorphism. The combined net effect of positive and negative feedback mechanisms within snowpacks is unknown and difficult to study in nature because measurement techniques easily disturb the natural snow structure. The question remains whether or not enhanced solar heating from LAI at the surface could slow metamorphism by weakening the temperature gradient and associated vertical flux of vapor. In light of these current unknowns, the purpose of this study is twofold: One, to investigate the effects of added LAI on snow albedo feedback; and two, to demonstrate the utility of a new instrument we use 30 to obtain approximate snow SSA. We hypothesize that if we add dust and BC to snow surfaces, then we will induce measurable snow albedo positive feedback.
To test this hypothesis, we first describe the design principle and calibration of the Near-Infrared Emitting and Reflectance Monitoring Dome (NERD), an instrument that is placed gently onto the snow surface to obtain snow SSA. The NERD enables multiple 1.30 and 1.55 µm bidirectional reflectance measurements in just minutes with minimal alteration of the snow structure.
To calibrate with respect to snow SSA, we compare snow BRFs with X-CT derived SSA and find an exponential relationship between SSA and snow 1.30 µm BRFs. These relationships are also explored using a three dimensional Monte Carlo photon 5 transport model. We then present results from our LAI in snow experiments. We discuss these results and their implications.
Overall, this study demonstrates conditions for which snow metamorphism can be enhanced by the presence of LAI.
2 Instrumentation and Methods

The Near-Infrared Emitting and Reflectance-Monitoring Dome (NERD)
The NERD is designed to measure 1.30 and 1.55 µm BRFs. These wavelengths are selected for snow SSA retrieval due to the 10 strong dependence of snow albedo on snow optical grain size (i.e. sphere equivalent radius). Snow spectral albedo is simulated here using the SNICAR model to demonstrate this sensitivity ( Fig. 1). While snow spectral albedo is sensitive to snow optical grain size (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.
The design principle is similar to the DUal Frequency Integrating Sphere for Snow SSA measurements (DUFISSS) (Gallet 15 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 degrees 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. The interior of the dome is painted with a flat black paint to increase absorptivity and minimize internal reflections between the dome and snow surface. Four infrared 20 light emitting diodes (LEDs) are mounted into a 17 cm diameter black styrene half-sphere. Two LEDs with peak emission wavelengths of 1.30 µm are mounted at nadir and ten degrees relative to zenith and two LEDs with peak emission wavelengths of 1.55 µm are mounted at 15 degrees off nadir (see Fig. 2). 1.30 µm LEDs have spectral line half widths of 85 nm and half intensity beam angles of ten degrees, while 1.55 µm LEDs have half-maximum bandwidths of 130 nm and 20 degree beam angles. These high powered, narrow beam infrared LEDs are selected to illuminate a small oval (estimated major axes of 1.5 25 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 InGaAs photodiodes mounted in two different azimuthal planes (0 and 90 degrees relative to the illumination); two each at 30 and 60 degrees relative to zenith. Photodiodes highly sensitive to light ranging from 800 to 1750 nm and relatively large active areas (1 mm) are selected to maximize sensitivity.
Because the orientation of LEDs and photodiodes are fixed, reflectance factors of surfaces with negligible subsurface scat-30 tering can be obtained after calibration using two diffuse reflectance targets in a manner similar to that used by Gallet et al. (2009), Gergely et al. (2014), and 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 degree viewing angles are obtained for each light source. While subsurface scattering of visible light in snow is pervasive, 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. Subsurface scattering is assumed to be minimal and fully contained within each photodiode's field of view. Therefore, this procedure enables simultaneous measurements of multiple snow BRFs at 1.30 and 1.55 µm. 5 To validate NERD reflectance measurement accuracy, precision, and responsiveness, measured BRFs of reflectance standards are recorded after calibration. Ten BRFs (R) for each LED / 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 10 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 %.

Snow samples
Surface snow (just the top few centimeters) samples were collected in nature over the span of three years (winters 2015-2017) 15 and transported in coolers to the nearby US 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).

X-ray micro-computed tomography (X-CT)
To determine snow SSA, X-CT was conducted on a class of six snow samples according to Lieb-Lappen et al. (2017). X-ray 20 (40-45 kV, 177-200 micro-Amps) transmission through cylindrical snow samples was measured at rotation steps of 0.3-0.4 degrees. To limit scan times to 15 minutes, exposure time was set to 340 ms at a cubic voxel resolution of 14.9 µm. Processing software enables SSA calculations from three dimensional 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 t set equal to the total time 25 elapsed between NERD measurements and X-CT scan times and with τ and n inferred from Ebner et al. (2015).

Contact spectroscopy
Snow SSA was also inferred from optical grain size measurements using contact spectroscopy (Painter et al., 2007). Snow reflectance spectra are collected using an ASD FieldSpec4 and high-intensity contact probe with reference to a Spectralon white reference panel. The effective radius is determined from the normalized area of the absorption feature centered at approximately 30 1 nm using a look up table (Nolin and Dozier, 2000). These measurements were conducted on depth hoar created in a cold lab (DH_2016) and on rounded grains (RG_2015).

Monte Carlo Modeling of Bidirectional Reflectance Factors
The Monte Carlo method is applied in this study to numerically simulate light emission by the NERD and the resulting three dimensional light scattering within modeled snow packs. Arrays of photons with wavelengths generated at random using 5 Gaussian distributions are used to mimic the 85 and 130 nm full width at half-maximum spectral emission characteristics of the narrow-band LEDs mounted in the NERD. These LEDs are modeled as photon emitters according to their orientation in the dome. Photons are initiated downward into the snow medium, as demonstrated by Kaempfer et al. (2007), and propagated in optical depth space. Photon particle interactions are determined using random number generators. Photons can either be absorbed or scattered with the probability determined by the particle single scattering albedo. Photons are terminated upon 10 absorption and followed if scattered. When a photon is scattered, its new direction cosines are determined by the particle scattering phase function. 1,000,000 photons per simulation are propagated and followed through the snow medium until they are absorbed or exit the medium. The snow packs are modeled as homogenous matrices of suspended particles with input data containing the particle mass absorption cross section, asymmetry parameter, single scattering albedo, projected area, volume, and scattering matrix. These scattering properties are calculated by Yang et al. (2013) for randomly oriented ice particle shape 15 habits that include droxtals, solid hexagonal columns, and spheres. For spheres, we apply the Henyey-Greenstein phase function where θ is the scattering angle and g is the relevant asymmetry parameter. We select these subset of shape habits from the larger dataset provided by Yang et al. (2013) (1977). Accordingly, photon exit angles are grouped into 30 exit zenith angle (θ r ) bins at three degree resolution. Azimuthal mean BRFs are calculated by zenith angle θ r from the total incident photon flux Φ i by where Φ r represents the azimuthally integrated photon flux through each θ r bin. 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 two is necessary to normalize the resulting weighting function w(θ r ) = sin θ r cos θ r , as π/2 0 sin θ r cos θ r dθ r = 1 2 .
the BRFs' dimensionality, so that fewer photons are needed to mitigate Monte Carlo noise. Equation (5)  Early on February 10, experimental plots were loaded with BC until visible darkening was apparent. Snow BRFs were measured shortly after 00:00, 03:00, 06:00, and then periodically throughout the day. Because these plots were well shaded by tall trees, 20 these measurements were used to monitor snow metamorphism without the influence of direct solar illumination. Furthermore, mostly cloudy conditions on February 10 diffused incoming radiation so that ambient lighting was nearly isotropic.

February 17 experiment (clear sky / direct solar heating)
On February 17, just a pinch of BC and 30g of sand were deposited on separate experimental plots. These surface fluxes were selected to mimic extreme LAI deposition events observed by Skiles and Painter (2017). As in the previous experiment, snow 25 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.

Monte Carlo modeling
To validate the Monte Carlo model for snow applications, 1.30 µm narrow band black-sky albedo was calculated and compared to the SNICAR model for snow Re ranging from 36 to 327 µm (SSA = 80 to 10 m 2 kg −1 ) (see Fig. 3

, left). As expected, Monte
Carlo results from snow modeled as spherical ice particles were consistent with narrow-band albedo calculations from Flanner 5 et al. (2007). These results show slightly higher hemispheric reflectances for droxtals (for all SSA) and solid hexagonal columns (for SSA > 40 m 2 kg −1 ) than those calculated from equal SSA spheres and from the SNICAR model. We hypothesize that the variations in these albedo calculations, and also BRFs, across particle shapes are inversely related to the particles' asymmetry parameters.
To inform on our choice of a snow SSA calibration function, Monte Carlo simulated snow BRFs are calculated for various 10 (homogeneous) particle SSA ranging from 10 to 80 m 2 kg −1 and plotted against measurement data in Fig. 4 BRFs are larger than predicted from modeling across all SSA.

NERD Snow SSA Calibration
To calibrate the NERD for snow SSA retrieval, we compare X-CT derived snow SSA with NERD snow BRF measurements (see Fig. 4). In general, snow BRFs are directly related to snow SSA. At 1.30 µm BRFs range from just under 0.2 (for low SSA) to as high as 0.7 (for high SSA) and are slightly higher at 60 degrees than at 30 degrees. We observe 1.55 µm snow BRFs reflectance and snow SSA and use the longer wavelength in their DUFISSS to obtain measurements of high snow SSA (> include snow of higher SSA to determine the utility of 1.55 µm snow BRFs in measuring fresh snow of extremely high SSA particularly common in the extremely cold Arctic and Antarctic environments and observed by Legagneux et al. (2002) and Libois et al. (2015).
Snow BRFs at 1.55 µm observed by the NERD are higher than both hemispherical reflectance measurements by Gallet et al. 5 (2009) and those predicted from Monte Carlo modeling. NERD LEDs have full width at half maximums of 130 nm and emit non-negligible light at wavelengths much shorter, toward the near-infrared. We hypothesize that higher than expected measured 1.55 µm BRFs are caused by reflected light at shorter wavelengths. Additional SNICAR modeling results and Monte Carlo simulations support this hypothesis (see Fig. 3, right).
In light of these empirical and numerical results, we propose the following general exponential form relating 1.30 µm snow 10 BRFs to SSA, such that for 1.30 µm snow BRF R 1.30 . Using least squares regression analysis, we compute parameters α and β for both 30 and 60 degree viewing zenith angles (see Fig. 5).
Ideally, an empirically derived calibration function would include SSA measurements from multiple methods to mitigate un- 15 certainties 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 SSA. As a preliminary validation of eq. 6, we compare snow SSA results to SSA derived from snow optical effective radii 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. As expected, contact spectroscopy snow SSA values are consistently higher than those calculated from X-CT analysis and therefore eq. 6.

20
Because contact spectroscopy measurements typically yield higher SSA values than those derived from other optical methods, these comparisons, though preliminary, offer some initial validation of the NERD snow SSA calibration.
Hereafter, we apply eq. 6 in the following LAI in snow experiments to estimate hourly snow SSA from measured snow BRFs. Because the remainder of this study is concerned with relatively large changes in SSA, approximate SSA retrieval using the NERD are sufficient to quantitatively assess snow metamorphism in the presence of LAI.

LAI enhanced snow metamorphism
First, to monitor snow metamorphism without solar heating, during the early morning (night) hours on February 10, we deposit BC onto an experimental plot after the previous day's snow fall. 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. 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 30 and natural snow (see Fig. 6). X-CT analysis showed small differences in morning (49 m 2 kg −1 ) and afternoon (48 m 2 kg −1 ) snow SSA. Our results from this experiment indicate that heavy BC loading had little to no effect on snow metamorphism without direct solar irradiance.
Second, to monitor snow metamorphism occurring after forced large BC and dust deposition events under direct solar illumination, on February 17, 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. In natural snow, 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). 1.55 µm BRFs remained just above 0.1. In the 5 dust loaded plot, snow 1.30 (1.55) µm BRFs decreased rapidly from above 0.5 (0.1) before 10:00am to to below 0.3 (0.05) by 1:00pm 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 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 2 kg −1 , which thereafter decreased to 41, 23, and 18 m 2 kg −1 in natural, BC loaded, and 10 dust loaded snow, respectively. NERD derived 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. 15 These results suggest that realistic LAI deposition can accelerate snow metamorphism. The primary cause of this accelerated process is enhanced solar radiation by LAI. Surprisingly, added BC had little to no effect on snow metamorphism during cloudy conditions. One possible explanation of this surprising result is that adding BC to snow only initiated melting during clear sky conditions. In the clear sky experiments, LAI 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 20 indicative of either the accumulation of liquid water from melting snow or decreasing snow surface SSA.
Accelerated snow metamorphism by dust loading is consistent with the findings of Skiles and Painter (2017). The indirect effect of LAI on snow is also demonstrated by Hadley and Kirchstetter (2012), where 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 snow pack that increases the strength of snow internal albedo feedback. Typical BC deposition events are very small, 25 however, so it is difficult to reproduce natural BC concentrations when adding any BC to a one square meter plot.

Conclusions
Taken together, these results indicate that LAI deposition can accelerate snow metamorphism and enhance positive snow albedo feedback, especially during cloud free, calm weather conditions when surface air temperatures are near 0 • C. To obtain quick, repeatable measurements of snow SSA without destroying samples, we engineered an instrument (i.e., the NERD) 30 that measures 1.30 and 1.55 µm BRFs of snow. We evaluated NERD 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 2 kg −1 .
The NERD will serve to further study the effects of LAI on snow metamorphism and to explore the spatial heterogeneity of snow SSA. Because it 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 validation is 5 needed. Additional independent measurement methods that include snow samples with a larger snow SSA span from a variety of environmental conditions and further experimentation into the small scale effects on NERD snow BRF measurements are needed to fully justify the NERD as an accurate snow SSA measurement technique. Further investigation into the microphysical limitations and quantitative uncertainties associated with the precise retrieval of snow SSA from near-infrared BRF measurements is the subject of a follow-on study.
10 Code and data availability. Plot data referenced in this manuscript and associated Python scripts used to generate figures are made available via the University of Michigan's Deep Blue data repository (Schneider and Flanner, 2018).

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 its active region. Because these light signals are reflected from the experimental surface, the currents induced 15 by the photodiodes are very small (nano-to micro-Amps). 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 20 and the analog-to-digital converter (ADC) to reduce noise. This filter is designed to have a time constant of less than 0.5 seconds to achieve balance between adequate noise reduction and speed.
Competing interests. Mark Flanner is currently an editor for The Cryosphere. We are not aware of any other competing interests associated with the publication of this manuscript.
facilities at CRREL. We also thank Ross Lieblappen for sharing his micro-computed tomography expertise through providing a thorough tutorial for running and analyzing snow scans.
The authors are thankful for reviews from three referees (including Marie Dumont). These detailed reviews were helpful in reorganizing the manuscript and improving the presentation of the key results.