Measurement of snowpack density, grain size, and black carbon concentration using time-domain diffuse optics

Abstract

Diffuse optical spectroscopy (DOS) techniques characterize scattering media by examining their optical response to laser illumination. Time-domain DOS methods involve illuminating the medium with a laser pulse and using a fast photodetector to measure the time-dependent intensity of light that exits the medium after multiple scattering events. While DOS research traditionally focused on characterizing biological tissues, we demonstrate that time-domain diffuse optical measurements can also be used to characterize snow. We introduce a model that predicts the time-dependent reflectance of a dry snowpack as a function of its density, grain size, and black carbon content. We develop an algorithm that retrieves these properties from measurements at two wavelengths. To validate our approach, we assembled a two-wavelength lidar system to measure the time-dependent reflectance of snow samples with varying properties. Rather than measuring direct surface returns, our system captures photons that enter and exit the snow at different points, separated by a small distance (4–10 cm). We observe clear, linear correlations between our retrievals of density and black carbon concentration, and ground truth. For black carbon concentration the correlation is nearly one-to-one. We also find that our method is capable of distinguishing between small and large grain sizes.

1. Introduction

Snow is composed of transparent ice grains that absorb light very weakly at visible wavelengths (Warren, 2019). Because of this, photons that enter a snowpack will typically scatter many times off of a large number of ice grains before they either exit the medium or get absorbed. The study of snow optics has historically focused on the interaction of snow with sunlight, as understanding this interaction is essential to understanding snow cover's contribution to the Earth's climate (Henderson and others, 2018) and for forecasting snow melt (Painter and others, 2010), among other things. A key goal of snow optics has been to predict spectral albedo as a function of snowpack properties – such as grain size, which determines the probability that a photon will be absorbed at each scattering event (Wiscombe and Warren, 1980a), and the concentration of light absorbing particles such as dust, black carbon, or algae that are mixed into the snow (Wiscombe and Warren, 1980b; Skiles and others, 2018).

The development of accurate spectral albedo models has, in turn, led to the development of optical sensing methods that retrieve grain size (Nolin and Dozier, 2000; Gallet and others, 2009) and LAP concentrations (Zege and others, 2011; Painter and others, 2012) from spectral albedo measurements. These methods, while useful, have limitations. Snowpack albedo is largely independent of important properties such as snow density (Wiscombe and Warren, 1980a). Furthermore, spectral albedo measurements usually require passive illumination by sunlight, and as such cannot be used to retrieve snow properties at night and for several months of the year in polar regions. Albedo models developed for solar illumination assume steady-state illumination that is collimated, diffuse, or a mixture of the two. As such, they cannot fully model lidar waveform measurements, which consist of the time-dependent optical response of snowpack to focused, pulsed illumination.

Over the past few decades, in parallel with advances in snow optics, the biomedical optics community has developed a suite of techniques for characterizing biological tissue, which, like snow, is also a highly scattering medium. Collectively, these methods are referred to under the umbrella term of diffuse optical spectroscopy (DOS) (Durduran and others, 2010), which refers to the fact that the propagation of photons within the scattering medium is modeled using the diffusion approximation to the radiative transfer equation (Welch and van Gemert, 1995), and to the fact that multi-wavelength illumination is frequently used (although this is not required). In DOS techniques the tissue is probed with a focused laser source that can be time-modulated, frequency-modulated, or continuous-wave. Measurements of the tissue's optical response are then used to estimate its optical properties, such as the tissue's absorption coefficient or effective scattering coefficient. These optical properties, in turn, can be related to clinically useful properties of the tissue such as blood oxygenation (Sevick and others, 1991), organelle size (Li and others, 2008), and the concentrations of water, lipids, and collagen (Quarto and others, 2014). DOS has also been applied in non-clinical settings for the inspection of produce (Nicolaï and others, 2014), and for characterizing porous materials such as wood (Bargigia and others, 2013) and pharmaceutical tablets (Johansson and others, 2002).

Because snow is also a highly scattering medium, many of the results from diffuse optical spectroscopy can be adapted to the characterization of snowpack properties. Despite this, the adoption of diffuse optics concepts in the snow sensing community has been limited. Várnai and Cahalan (2007) proposed that the spatial spread of diffused laser light could be used to determine snow and sea ice thickness. Smith and others (2018) noted that the multiple scattering of green laser light within a snowpack should result in biases in lidar altimetry measurements. They used a combination of diffusion theory and Monte-carlo modeling to assess the dependence of this multiple scattering bias on grain size, black carbon concentration, and the choice of surface height retrieval algorithm. Smith and others (2023) used the model of Smith and others (2018) to develop an algorithm that infers snow grain size from full waveform lidar measurements collected by the Airborne Topographic Mapper (ATM). Fair and others (2024) use lidar retrievals of grain size to predict biases in snow surface heights retrieved using green (532 nm) lidar beams on IceSat-2 and ATM. As far as we are aware, prior to this work, the only direct application of DOS techniques to retrieve bulk snowpack properties was made by Allgaier and Smith (2022). In their work, the snow was illuminated with continuous-wave laser sources at two different wavelengths, and a smartphone camera was used to take images of the spatially resolved, steady-state intensity of light that exited the snowpack after diffusing within the snow. From these smartphone images, along with an independent in situ measurement of the snow's density, the authors were able to retrieve the absorption and effective scattering coefficients of the snowpack, as well as an estimate of the concentration of black carbon within it. In Ackermann and others (2006), and in a separate work by Allgaier and others (2022), time-domain diffuse optical measurements were used to estimate the absorption and scattering coefficients of glacier ice, which is optically similar to snow. A theoretical analysis of diffuse optical spectroscopy applied to glacier ice is provided in Allgaier and Smith (2021). Studinger and others (2024) use multiple scattering returns in green lidar measurements to infer the scattering length within sea ice.

In this work we introduce a new method for characterizing the bulk properties of snow that is based on time-domain diffuse optical measurements. Our instrument is effectively a photon-counting lidar system that consists of two pulsed lasers with different wavelengths (one red, one near-infrared), and a single-photon avalanche diode (SPAD) receiver. Rather than measuring surface returns, which might be used for altimetry, we measure photons that enter the snowpack at a single point on the surface and exit at a second surface point that is displaced from the point of entry by a small distance (4–10 cm). Through a series of proof-of-principle experiments, we show that our method is capable of retrieving the density (through the ice volume fraction), grain size, and the concentration of light absorbing particles of a dry snowpack, in a non-invasive way.

2. Methods

2.1. Diffusion model

The propagation of a laser pulse inside a scattering medium is described by the time-dependent radiative transfer equation (Welch and van Gemert, 1995), which models the flow of radiance (W m⁻² sr⁻¹) within a medium as a function of space and time. The scattering medium is described by a scattering coefficient μ _s (m⁻¹), a scattering phase function, an absorption coefficient μ _a (m⁻¹), and the speed of light within the medium $c_\ast$ (m s⁻¹).

Under the diffusion approximation to the radiative transfer equation, photons are modeled as particles that ‘diffuse’ through a scattering medium via random walks. This approximation accurately describes situations for which the distance scales considered are much larger than the mean free path of photons within the medium (=(μ _a + μ _s)⁻¹), and photons are typically scattered many times before they are absorbed (μ _s ≫ μ _a) (Welch and van Gemert, 1995). The photon diffusion equation can be written as follows:

(1)$${1\over c_\ast } {\partial\over \partial t} \phi( {\bf r},\; t) - D \nabla^2\phi( {\bf r},\; t) + \mu_a \phi( {\bf r},\; t) = S( {\bf r},\; t) .$$

A derivation of the photon diffusion equation can be found in Haskell and others (1994). Unlike the time-dependent radiative transfer equation, which models the time-evolution of a five-dimensional radiance field, the photon diffusion equation models the lower dimensional quantity of photon fluence ϕ(r,  t) (W m⁻²), which is the integral of radiance over all directions. The variable S represents an isotropic source term, and the diffusion constant D is defined as

(2)$$D = {1\over 3\left[\mu_a + ( 1 - g) \mu_s \right]}.$$

Here g is the asymmetry factor of the scattering phase function, which can take values between −1 and 1 depending on whether the medium is primarily backward scattering (g < 0), isotropically scattering (g = 0), or forward scattering (g > 0).

Crucially, the photon diffusion equation permits analytical solutions when the geometry of the scattering medium is sufficiently simple. We consider the scenario depicted in Fig. 1a. Here, the medium is assumed to be semi-infinite and homogeneous. The medium's surface is illuminated by a pulsed, pencil-beam source at time t = 0, and a detector observes the time-dependent intensity of light that exits the medium at a second point that is displaced from the point of illumination by a distance s. Kienle and Patterson (1997) showed that, in this scenario, Eq. (1) can be accurately solved by imposing an extrapolated boundary condition, which requires that photon fluence goes to zero along a planar boundary that lies just above the medium's surface. The source term S is approximated by a point source buried one transport mean free path beneath the surface, and the equation is then solved via the method of images. This yields the following expression for photon fluence inside the medium:

(3)$$\eqalign{\phi( s,\; z,\; t) & = {c_\ast \over ( 4\pi Dc_\ast t) ^{3/2}}\exp{\big(-\mu_a c_\ast t\big)} \cr & \quad\times \big\{\exp{\big[-{( z-z_0) ^2 + s^2\over 4 D c_\ast t}\big]} \cr & - \exp{\big[-{( z + z_0 + 2z_b) ^2 + s^2\over 4 D c_\ast t}\big]} \big\}.}$$

Figure 1. (a) Illustration of the measurement geometry employed in this work. A point on the snow surface is illuminated by a laser pulse at time t = 0. A detector observes the time-dependent intensity of light that exits the snow from a second point at distance s from the laser spot. (b) Comparison of time-dependent intensity predicted by our model (Eq. (5)), to photon time-of-flight measurements generated using a Monte-carlo simulation of a scattering medium with the same properties. Both curves are normalized to their respective peaks.

Here z denotes the distance from the surface, going down; z ₀ = [μ _a + (1 − g)μ _s]⁻¹ is the depth of the buried point source; and z _b denotes the height of the extrapolated boundary. Haskell and others (1994) proposed a value of z _b = (1 + R _eff)/(1 − R _eff)2D, where R _eff is the fraction of photons that are internally reflected at the interface between the scattering medium and the external (non-scattering) medium due to a refractive index mismatch. Because our ultimate goal is to model the optical response of a snowpack, and because the snow-air boundary of a typical snowpack is not a dielectric interface at optical wavelengths, we assume for this work that R _eff = 0 and hence, z _b = 2D.

From Eq. (3), we compute the time dependent radiosity (Wm⁻²) that exits the surface at position s using Fick's Law (Kienle and Patterson, 1997):

(4)$$J( s,\; t) = -D\nabla\phi( s,\; z,\; t) \cdot( -\hat{z}) \vert _{z = 0} .$$

The reflected flux R measured by a detector that observes the medium's surface from a distance can then be described using the following expression:

(5)$$\eqalign{R( s,\; t) & = {\alpha c_\ast \over 3( 2\pi) ^{3/2}}{z_0^2\over ( 2Dc_\ast t) ^{5/2}}\exp\left(-\mu_ac_\ast t -{s^2 + z_0^2\over 4Dc_\ast t}\right)\cr & \quad\times \left[1 + {7\over 3}\exp\left(-{10z_0^2\over 9Dc_\ast t}\right)\right],\; }$$

where α is a constant that encapsulates instrumental parameters such as transmitted laser power, detection efficiency, and the detector's etendue. We note that we have made liberal use of the substitutions D = z ₀/3 and z _b = 2z ₀/3. In deriving Eq. (5), we also assumed that the surface could be accurately described as a Lambertian emitter, which means that the radiance emitted by the surface is independent of the emission angle. Previous work has relaxed this assumption (Kienle and Patterson, 1997). We found that doing so produced nearly identical results when describing a nadir-pointing detector, but added significant complexity to the model. For this reason, we elected to use Eq. (5).

In Fig. 1b we compare the time-dependent intensity predicted by Eq. (5) to simulated photon time-of-flight measurements generated using a Monte-carlo simulation (Henley, 2020). The modeled results match the simulation very closely. In general, models derived from the diffusion approximation to the radiative transfer equation accurately describe the measurements of photons that arrive at later times ($c_\ast t\gg z_0$) and larger distances from the laser spot (s ≫ z ₀), as these photons have scattered many times before exiting the medium.

2.2. Snow scattering model

Our measurement model, defined in Eq. (5), is expressed in terms of three phenomenological parameters – the absorption coefficient μ _a, the effective scattering coefficient μ _s′ = (1 − g)μ _s, and the effective speed of light in the medium $c_\ast$. We use a scattering model derived from the geometric-optics scattering model of Kokhanovsky and Zege (2004) to define μ _a, μ _s′, and $c_\ast$ in terms of three physically meaningful snowpack parameters – $v_\ast$, the fraction of the snowpack volume that is occupied by ice; $r_\ast$ (m), the grain radius; and C _bc (kg kg⁻¹), the mass mixing ratio of black carbon in the snowpack. We note that, for a dry snowpack, the ice volume fraction $v_\ast$ is readily converted to bulk snowpack density $\rho _\ast$ (kg m⁻³) via the expression $\rho _\ast = v_\ast \rho _{ice} + ( 1-v_\ast ) \rho _{air}\approx v_\ast \rho _{ice}$, where ρ _ice and ρ _air are the intrinsic densities of ice and air, respectively. We do not consider wet snow in this work.

2.2.1. Clean snowpack

For a dry snowpack that contains optically insignificant concentrations of light absorbing particles, the scattering and absorption coefficients can be written entirely as functions of $v_\ast$ and $r_\ast$. The absorption and effective scattering coefficients are computed as follows:

(6)$$\mu_a = B\Gamma v_\ast $$

(7)$$\mu'_s = {3\over 2}( 1-g) {v_\ast \over r_\ast } ,\; $$

The grain radius $r_\ast$ can be interpreted as the characteristic size of the ice grains. As in Kokhanovsky and Zege (2004), $r_\ast$ is defined as the radius of the spherical ice grain that would have the same surface-area-to-volume ratio as the ice-air matrix that comprises the true snowpack. Explicitly $r_\ast = 3{\langle V\rangle }/{\langle \Sigma \rangle }$, where 〈V〉 is the mean ice grain volume and 〈Σ〉 is the mean ice grain surface area.

The absorption enhancement parameter B and scattering asymmetry factor g are determined by grain shape (Libois and others, 2013). A recent study by Robledano and others (2023) suggests that these parameters cluster around B = 1.7 and g = 0.825 for most real snow samples, so we use those values here. These values are approximately valid for visible and near-infrared wavelengths (400 nm to 14000 nm) (Robledano and others, 2023). Notably, the values determined by Robledano and others (2023) closely match theoretical predictions for a two-phase random mixture of ice and air, in which grains have random and irregular shapes rather than idealized shapes such as spheres or hexagonal plates (Malinka, 2014).

In Fig. 2a we visualize the range of values for μ _s′ obtained across a domain of feasible grain sizes and ice volume fractions. Figure 2(b) shows the dependence of μ _a on ice volume fraction and wavelength. Unlike μ _s′, the absorption coefficient depends strongly on wavelength, and varies by more than an order of magnitude between the red (λ = 640 nm) and near infrared (λ = 905 nm) wavelengths used in this study.

Figure 2. (a) Effective scattering coefficient μ _s′ (m⁻¹) as a function of ice volume fraction $v_\ast$ (unitless) and grain radius $r_\ast$ (mm). (b) Absorption coefficient μ _a (m⁻¹) of clean snow as a function of ice volume fraction $v_\ast$ (unitless) and wavelength λ (nm).

2.2.2. Effective speed of light in snow

The last parameter to calculate is the effective speed of light within the snowpack $c_\ast$. In many problems that involve light propagation in a scattering medium, light's speed is treated as a constant that can be computed beforehand if the medium's index of refraction is known. This approach does not work for snow, which is a heterogeneous mixture of two materials – ice and air – that have markedly different refractive indices and that may be mixed at any ratio.

The geometric scattering model proposed by Kokhanovsky and Zege (2004), from which we defined μ _a and μ _s′ (Eqs. (6), (7)), implicitly defines the distance that a photon travels through an ice grain as the distance (i.e. l _ice) of the chord that connects the points at which a photon enters (point 1 in Fig. 3) and exits (point 3) the grain. This effective ‘transportation distance’ differs from the true distance (i.e. l _ice′) traveled through the grain if the photon is internally reflected (e.g. at point 2) before it exits the grain. The absorption enhancement parameter B is approximately equal to the ratio between the true and effective transportation distances, averaged over all possible internal paths (i.e. l _ice′ ≈ Bl _ice) (Libois and others, 2019).

Figure 3. A comparison of the true path (red) traveled by a photon through an ice grain, including internal reflections, to the effective transportation path (black, dashed) of length l _ice that is implicitly assumed by our absorption and effective scattering coefficient models.

An effective light speed model that is compatible with our definitions of μ _a and μ _s′ must describe the average speed at which light advances along this effective transportation path, which is equivalent to the true photon path in the air phase, but shorter than the true photon path in the ice phase by a factor of B. If a photon travels along an effective transportation path of length L, on average that path will pass through $( 1-v_\ast ) L$ of air and $v_\ast L$ of ice. The time T require to traverse this path is

(8)$$T = {( 1-v_\ast ) L\over c_0} + {n_{ice} B v_\ast L\over c_0} ,\; $$

where n _ice is the real component of the refractive index of ice (Warren and Brandt, 2008) and c ₀ is the speed of light in air (where it's assumed that n _air = 1). The travel time within the ice phase has been increased by a factor B to account for the difference between the true and effective transportation path lengths. Dividing L by T leaves us with

(9)$$c_\ast = {c_0\over 1 + ( n_{ice}B-1) v_\ast } .$$

We stress that the effective light speed defined here is lower than the mean speed of light computed with respect to the lengths of true photon paths through snow, which due to internal reflections can include jagged paths within grains. An expression for this true mean light speed was computed by Libois and others (2019), and is equal to the effective light speed of Eq. (9), multiplied by a factor of $[ 1 + ( B-1) v_\ast ]$.

2.2.3. Effect of light absorbing impurities

Ice is an exceptionally weak absorber of light at visible wavelengths (Warren, 2019). As such, the absorption of visible light within a snowpack can be enhanced significantly – even dominated – by the presence of trace concentrations of more absorptive substances. This has the important effect of reducing snowpack albedo, which increases radiative forcing on the snow surface and subsequently enhances snow melt and metamorphism and can also influence the local climate (Skiles and others, 2018). For our purposes, the presence of small concentrations of LAPs can increase the absorption coefficient of a snowpack considerably – thus rendering Eq. (6), our model for clean snowpack absorption, insufficient. Globally, radiative forcing from LAPs is dominated by black carbon, mineral dust, organic or ‘brown’ carbon, and snow algae (Skiles and others, 2018). Here we assume that absorption by LAPs is dominated by black carbon, but note that our model could be extended to include other types of particles by modifying the LAP absorption spectrum used here.

According to Flanner and others (2012), between 32–73% of the black carbon in global surface snow is embedded within ice grains (or ‘internally mixed’), with the remainder being external to those grains (‘externally mixed’) in the air phase. The elongated paths followed by photons within ice increases the probability that photons will interact with internally mixed black carbon particles. As a consequence, internally mixed black carbon has an outsized impact on snow absorption and albedo, relative to externally mixed black carbon (Flanner and others, 2012). Models for snow's absorption coefficient that consider the mixing state of black carbon have been proposed (Liou and others, 2014; Dombrovsky and Kokhanovsky, 2020), however these models typically require idealized grain shapes such as spheres – which do not accurately represent real snow – and assign highly non-linear dependencies on black carbon concentration that are grounded in electromagnetic theory (Dombrovsky and Kokhanovsky, 2020) or stochastic simulations (Liou and others, 2014).

Here we propose a simple geometric optics model for the additional absorption due to black carbon that can be computed from the bulk density of black carbon particles embedded inside ice grains $\rho _{bc}^{in}$ (kg m⁻³), the bulk density of black carbon particles external to the grains $\rho _{bc}^{out}$ (kg m⁻³), and the wavelength-dependent mass absorption efficiency MAE_bc (m² kg⁻¹) of the black carbon particles (Grenfell and others, 2011). Under this model, the presence of black carbon in a snowpack alters its properties primarily by adding an extra term to the absorption coefficient, i.e. $\mu _a^{snow} = \mu _a^{ice} + \mu _a^{bc}$. Our proposed model is written as follows:

(10)$$\matrix{\mu_a^{bc} & \hskip-2.3pc = {\rm MAE}_{bc}\left[( 1-v_\ast ) \rho_{bc}^{out} + Bv_\ast \rho_{bc}^{in}\right]\cr & = {\rm MAE}_{bc}\rho_{ice}v_\ast \left[( 1-v_\ast ) C_{bc}^{out} + Bv_\ast C_{bc}^{in}\right]. }$$

Here absorption by internally mixed impurities is multiplied by the absorption enhancement factor B to account for the elongation of photon paths within ice grains. In the lower part of Eq. (10), we replace $\rho _{bc}^{in}$ and $\rho _{bc}^{out}$ with the products of the intrinsic density of ice (ρ _ice = 916.5 kg m⁻³), the snowpack's ice volume fraction $v_\ast$, and the mass mixing ratios $C_{bc}^{in}$ and $C_{bc}^{out}$ of internally and externally mixed black carbon, respectively.

To simplify our model further, we assume that the black carbon is evenly mixed, i.e. $C_{bc}^{in} = C_{bc}^{out}$. We then combine Eqs. (6) and (10) to obtain a complete expression for the snowpack absorption coefficient:

(11)$$\mu_a = B\Gamma v_\ast + {\rm MAE}_{bc}\rho_{ice}C_{bc}v_\ast \left[1 + ( B-1) v_\ast \right].$$

Following the example of Doherty and others (2014), we model the wavelength dependence of MAE_bc using a power law spectrum:

(12)$${\rm MAE}_{bc}( \lambda) = {\rm MAE}_{bc}( \lambda_{ref}) \left(\lambda_{ref}/\lambda\right)^{\rm \AA},\; $$

that has an Ångstrom coefficient Å =1.1 and is referenced to MAE_bc(λ _ref) = 6500 m² kg⁻¹, where λ _ref = 600 nm.

Figure 4 illustrates that, for a fixed C _bc, the fraction of absorption attributable to black carbon in snow depends strongly on the wavelength of light that interacts with the snowpack. We plot the ratio of absorption due to black carbon (using Eq. (10)) to the total absorption (from Eq. (11)) for a selection of wavelengths that range from 400 nm (blue) to 1000 nm (near infrared). In computing these ratios, we assume an ice volume fraction of $v_\ast = 0.3$. For blue light, our model suggests that absorption is entirely dominated by just 1 part per billion by weight (ppbw) of black carbon, which is comparable to mass mixing ratios found in Greenlandic snow (Warren, 2019). In contrast, at 1000 nm, absorption from black carbon only eclipses ice absorption for mass mixing ratios above 7500 ppbw – a very high level of soot that would cause the snow to appear visibly grey. This decreased sensitivity at longer wavelengths is not caused by the decreased absorption of black carbon at these wavelengths, but rather by the increased absorption efficiency of ice.

Figure 4. The ratio of light absorption due to black carbon to total absorption by the snowpack for a range of wavelengths. Ratio computation assumes ice volume fraction $v_\ast = 0.3$.

2.2.4. Effect of snow properties on time-domain response

Having obtained expressions that relate μ _a, μ _s′, and $c_\ast$ to the grain size, ice volume fraction, and black carbon concentration of a dry snowpack, we can now develop an understanding of how changes to $v_\ast$, $r_\ast$, and C _bc affect the snowpack's time-domain optical response. Upon inspection of Eq. (5), we see that the shape of a snowpack's transient response is primarily controlled by $\mu _ac_\ast$, which determines the rate of decay of the signal's tail; $2Dc_\ast$, which can be interpreted as the rate at which a Gaussian cloud of diffusing photons expands over time, and which controls the position of the signal's peak; and $z_0^2$, which influences the shape of the response at the earliest arrival times, but in practice has little effect when s ≫ z ₀.

The exponential decay rate, $\mu _ac_\ast$, depends on $v_\ast$ and C _bc. On the other hand, $2Dc_\ast$ and $z_0^2$ primarily depend on the medium's scattering coefficient, which in turn depends on the ratio $v_\ast /r_\ast$. These effects are visualized in Figs. 5a–c, where we plot the predicted transient response curves for snowpack with varying $v_\ast$, $r_\ast$, and C _bc. For these curves, the snow is probed with red (640 nm) light, and the position of the detector's focus spot is fixed at s = 8  cm.

Figure 5. Plots of predicted time-dependent flux measured by a detector observing snow following illumination by a laser pulse (λ = 640 nm). Curves were produced using Eq. (5) and μ _a, μ _s′, and $c_\ast$ were computed from $v_\ast$, $r_\ast$, and C _bc using Eqs. (11), (7), and (9), respectively. (a) Ice volume fraction $v_\ast$ is varied. $r_\ast = 100\, \mu {\rm m}$, C _bc = 0  ppbw, and s = 8  cm. (b) Grain radius $r_\ast$ is varied. $v_\ast = 0.3$, C _bc = 0  ppbw, and s = 8 cm. (c) Impurity concentration C _bc is varied. $v_\ast = 0.3$, $r_\ast = 100\, \mu {\rm m}$, and s = 8  cm. (d) Detector focus position s is varied. $v_\ast$ = 0.3, $r_\ast = 100\, \mu {\rm m}$, and C _bc = 0 ppbw.

In Fig. 5a we see that for a clean snowpack (C _bc = 0), as $v_\ast$ is increased while $r_\ast$ is held constant, the slope of the signals’ exponential tail becomes more steep as light is absorbed by the medium more quickly. The arrival time of the signal peak is also pushed back because tighter packing of the ice grains reduces the distance between photon scattering events, thus reducing the rate at which light diffuses within the medium.

In Fig. 5b, $r_\ast$ is varied while $v_\ast$ is held constant and again C _bc = 0. As grain size increases, grains must be spaced further apart to maintain the same density, thus increasing the rate of diffusion within the medium. As such, for a fixed snow density and source-detector separation, the peak of the diffusion signal will arrive earlier, and will be more intense, when the grains are large.

In Fig. 5c, $v_\ast$ and $r_\ast$ are held fixed while C _bc is varied. Black carbon content only influences the absorption coefficient of the snowpack, and so increasing C _bc steepens the exponential decay rate $\mu _ac_\ast$. At the probing wavelength of 640 nm (where ice is a relatively weak absorber) the influence of C _bc is quite dramatic when compared to the comparable influence of ice volume fraction on the exponential decay rate, shown in Fig. 5a.

In Fig. 5d, $v_\ast$, $r_\ast$, and C _bc are held fixed and the detector focus position s is varied. As s increases, the signal peak arrives later and becomes more faint. However, as time passes, all signals converge as light spreads within the medium and the distribution of emitted photons becomes nearly uniform across the observed region.

2.3. Algorithm

We fit functions of the same form as Eq. (5) to two photon time-of-flight histograms – each measured using a different laser wavelength. We implement a grid search algorithm to find the fit parameters that minimize a negative log-likelihood loss function that properly accounts for photon count statistics. The parameters of the fitted curves are then used to compute the snowpack properties $v_\ast$, $r_\ast$, and C _bc. A visualization of our retrieval algorithm is shown in Fig. 6.

Figure 6. (Left) A diffusion model (lower right) is fit to photon time-of-flight histograms measured at two wavelengths. Fit parameters α′, β, γ, and δ are determined using a grid search algorithm that minimizes a negative Poisson log-likelihood function. (Top right) Snowpack properties $v_\ast$, $r_\ast$, and C _bc are computed directly from parameters β ₁, β ₂, γ ₁, and γ ₂ by evaluating analytical expressions (Eqs. (17), (18), and (19)).

2.3.1. Fit parameterization

We re-parameterize Eq. (5) in terms of the fitting parameters $\alpha ' = ( {\alpha c_\ast }) /( {3( 2\pi ) ^{3/2}})$, $\beta = \mu _a c_\ast$, $\gamma = 2Dc_\ast$, and $\delta = z_0^2$. This allows for the model to be expressed in simplified form:

(13)$$R( s,\; t) = \alpha'{\delta\over ( \gamma t) ^{5/2}}\exp\left(-\beta t -{s^2 + \delta\over 2\gamma t}\right)\left[1 + {7\over 3}\exp\left(-{20\delta\over 9\gamma t}\right)\right].$$

Although Eq. (13) appears to have four degrees of freedom, only the exponential decay rate parameter β and the spatial spread rate γ are used to estimate snowpack properties in practice. Interpreting the scaling constant α′ requires precise calibration of the instrument and measurement geometry which is challenging in practice and which we did not attempt. Equation (13) also depends only weakly on the squared source depth δ, to the point that accurate estimates of δ are almost never obtained.

The relative contributions of $v_\ast$ and C _bc to the decay rate parameter β vary significantly as a function of wavelength. However, for visible and near infrared light ($\lambda \lesssim 1100$ nm), the spatial spread rate γ depends primarily on μ _s′ and $c_\ast$, which are largely independent of the measurement wavelength. Altogether, this means that n + 1 independent parameters can be retrieved from measurements taken at n wavelengths. If absorption due to light absorbing particles is known to be insignificant (i.e. $\mu _a^{bc} \ll \mu _a^{ice}$), then $v_\ast$ and $r_\ast$ can be retrieved from measurements at a single wavelength. Otherwise, retrieving $v_\ast$, $r_\ast$, and C _bc requires measurements at two or more wavelengths.

2.3.2. Maximum likelihood estimate of model parameters

The number of counts in a histogram timing bin centered at t _i is assumed to be a Poisson random variable with a rate parameter x _i, defined as

(14)$$x_i = R( s,\; t_i ;\; {\bf \Theta}) + \eta = \alpha'r_i + \eta ,\; $$

where R is the flux predicted by Eq. (13) at position s, time t _i, and for parameters Θ = {α′,  β,  γ,  δ}. The rate of background counts produced by ambient light, detector dark counts, and detector afterpulsing (Zappa and others, 2007) is denoted by η, and is assumed to be constant with respect to time. The variable r _i denotes the normalized predicted flux, for which α′ = 1.

The probability of observing a vector of time-binned photon counts y given a vector of predicted count rates x is

(15)$$P( {\bf y}\vert {\bf x}) = \prod^N_\Delta {x_i^{y_i}e^{-x_i}\over y_i!}$$

where N denotes the total number of timing bins in the histogram and Δ is the starting bin for the curve fit. We seek to find the parameters Θ = {α′,  β,  γ,  δ} and η that minimize the negative log-likelihood

(16)$$\matrix{{\cal L}( {\bf \Theta},\; \eta \vert {\bf y}) & \hskip-2pc= -\ln{P\left({\bf y}\vert {\bf x}\left({\bf \Theta},\; \eta\right)\right)} \cr & = \displaystyle\sum^N_\Delta x_i - y_i\ln{x_i} + \ln{y_i!} . }$$

We minimize Eq. (16) using a grid search. To reduce the dimensionality of the search, we first estimate η by computing the mean number of counts in a designated set of noise bins that reliably contains effectively zero non-background counts. For any combination of β, γ, and δ, the scaling term α′ can then be estimated using the expression $\bar {\alpha }' = \sum ^N_\Delta ( y_i - \eta ) /\sum ^N_\Delta r_i$.

We can thus define a three-dimensional search area that contains all feasible values of β, γ, and δ. The feasible range for δ is tightly constrained to (3γ/2c ₀)² < δ < (3n _iceγ/2c ₀)², which allows for a coarse fit to be obtained using an effectively two-dimensional search. We uniformly sample the search area to create a grid of candidate fit parameters, and compute the negative log-likelihood for each set of parameters on the grid.

We perform a sequence of nested searches – we first obtain a coarse fit, then define a small search range around the fitted parameters and repeat the search using a smaller grid cell size. This procedure is iterated until a fit with the desired precision is obtained. Our fitting algorithm was implemented in MATLAB on a Lenovo Thinkpad T590 laptop with 16GB of RAM. Run time per fit was typically 56 seconds for 640 nm histograms, and 19 seconds for 905 nm histograms (which had fewer timing bins). Curve fits obtained using our algorithm are shown in Fig. 6. We estimated the uncertainty in the retrieved values of β, γ, and δ by computing the inverse of the Hessian of the loss function at the estimated minimum, and then taking the diagonal terms. These terms approximate the variances in parameter fits when the loss function is approximately Gaussian near the minimum (Bevington and Robinson, 1992).

2.3.3. Computing $v_\ast$, $r_\ast$, and C _bc

When measurements are obtained at two wavelengths, λ ₁ and λ ₂, the ice volume fraction $v_\ast$ and black carbon mixing ratio C _bc can be extracted from the decay parameters β ₁ and β ₂. Each term β _i can be expressed as a function of $v_\ast$ and C _bc by taking the product of Eqs. (9) and (11). This results in a set of two equations which can be solved, first, for $v_\ast$:

(17)$$v_\ast = {b_2\beta_1 - b_1\beta_2\over c_0\left(a_1b_2 - a_2b_1\right)- d_1b_2\beta_1 + d_2b_1\beta_2} .$$

For notational simplicity we have made the substitutions a _i = BΓ_i, b _i = ρ _iceMAE_bc(λ _i), and d _i = n _ice(λ _i)B − 1. As before, the term c ₀ refers to the speed of light in air.

Once $v_\ast$ has been obtained, C _bc can be computed as follows:

(18)$$\matrix{C_{bc} & \!\!\!\!\!\!\!\!\!\!\!\!= {1\over c_0b_1( 1 + f v_\ast ) }\left[\left({1\over v_\ast } + d_1\right)\beta_1 - c_0a_1 \right]\cr &\; = {1\over c_0b_2( 1 + f v_\ast ) }\left[\left({1\over v_\ast } + d_2\right)\beta_2 - c_0a_2 \right]. }$$

Here we have subsituted f = B − 1. After $v_\ast$ and C _bc have been computed, the grain radius $r_\ast$ can be computed from the spatial spread parameter γ _i at either wavelength:

(19)$$r_\ast = e_i \left[{2c_0\over 3\gamma_iv_\ast ( 1 + d_iv_\ast ) } -a_i - b_iC_{bc}( 1 + f v_\ast ) \right]^{-1} .$$

Here a _i, b _i, d _i and f are defined as they were previously, and e _i = 3(1 − g)/2. Because $r_\ast$ can be computed using either γ ₁ or γ ₂, we evaluate Eq. (19) at both wavelengths, and then take the uncertainty-weighted average of the two values obtained in this way to arrive at our final estimate for $r_\ast$. Uncertainties in $v_\ast$, $r_\ast$, and C _bc are obtained via error propagation from uncertainties in β ₁, β ₂, γ ₁, and γ ₂.

If it is known that absorption by light absorbing particles is small compared to absorption by ice grains, then $v_\ast$ and $r_\ast$ can be computed from the fit parameters extracted from single-wavelength measurements. First $v_\ast$ can be computed from the exponential decay rate β, as follows:

(20)$$v_\ast = {\beta\over ac_0 - \beta d} ,\; $$

and then $r_\ast$ can be obtained from the spatial spread rate γ, and our estimate of $v_\ast$:

(21)$$r_\ast = e\left[{2c_0\over 3\gamma v_\ast ( 1 + dv_\ast ) } -a \right]^{-1} .$$

Here a, b, d, and e retain their meanings from Eqs. (17) and (19).

2.3.4. Evaluation using simulated measurements

We validated our algorithm using a GPU-accelerated Monte-carlo photon transport simulation (Henley, 2020), which was adapted from a simulator originally developed for tissue imaging studies (Satat, 2019). We modeled the propagation of photons within a semi-infinite, homogeneous scattering medium. The medium's properties μ _a, μ _s′, and $c_\ast$ were computed from $v_\ast$, $r_\ast$, and C _bc using Eqs. (7), (9), and (11). To simulate pencil-beam illumination, photons were launched at the origin ([x,  y,  z] = 0) at time t = 0 and at normal incidence to the snow surface. Photons scattered randomly in the medium until they were absorbed, exited the medium, or satisfied an outlier termination criterion such as maximum number of scattering events. For more details, we refer the reader to Chapter 4 of Welch and van Gemert (1995).

We simulated photon time-of-flight histograms at 640 nm and 905 nm measurement wavelengths for two snow samples with different properties. We binned photons by the transverse position s (bin width 1 cm) and time t (bin width 16 ps) that photons exited the snow surface. In the first simulation, for 640 nm measurements, photons detected at s = 8.0 ± 0.5 cm were used for curve fitting, whereas for 905 nm measurements photons detected at s = 5.0 ± 0.5  cm were used. In the second simulation, for 640 nm measurements, photons detected at s = 10.0 ± 0.5 cm were used for curve fitting, whereas for 905 nm measurements photons detected at s = 7.0 ± 0.5 cm were used. Once a histogram of signal photons was created, a random number of background counts was added to each timing bin by sampling from a Poisson distribution with rate parameter η that was chosen to be consistent with the uniform background count levels observed in experimental measurements.

Our results are shown in Fig. 7. For the first simulation, the true snowpack properties were $v_\ast = 0.465$, $r_\ast = 240\, \mu {\rm m}$, and C _bc = 50  ppbw. Our method retrieved values of $v_\ast$ = 0.46±0.02, $r_\ast = 242\pm 9\, \mu {\rm m}$, and C _bc = 48 ± 3 ppbw. For the second simulation, the true snowpack properties were $v_\ast$ = 0.162, $r_\ast = 85\, \mu {\rm m}$, and C _bc = 0  ppbw. Our method retrieved values of $v_\ast = 0.164\pm 0.004$, $r_\ast = 86\pm 2\, \mu {\rm m}$, and C _bc = 0 ± 3 ppbw. These results suggest that our algorithm should produce accurate estimates under the idealized conditions prescribed here, and if our snow scattering model is correct.

Figure 7. We used a monte-carlo photon transport simulator to validate our retrieval algorithm. Measurements were simulated for 640 nm (left) and 905 nm (center) light for two simulated snow samples. True and estimated snowpack parameters for each sample are shown on the right.

3. Materials

3.1. Apparatus

3.1.1. Lidar system

We assembled a simple lidar system to measure the time-domain optical response of a variety of snow samples. We perform time-correlated single-photon counting, in which a histogram of photon times of flight is built up over time by repeatedly illuminating the snow surface with a pulsed laser. Photographs of our experimental setup are shown on the left in Fig. 8. Our lidar system used a single-pixel SPAD detector (Microphoton Devices PDM series) with a timing jitter of ~50 ps (FWHM), and two pulsed diode laser sources – a red laser with a wavelength of 640  nm (Picoquant LDH-P-C-640B), and a near-infrared laser with a wavelength of 905 nm (Picoquant LDH-P-C-905). Each laser was operated at a pulse repetition frequency of 2.5 MHz and had a quoted pulsewidth of <90 ps. The 640 nm laser was operated at a time-averaged power of 80 μW, and the 905 nm laser was operated at a time-averaged power of 55  μW. A Picoquant Hydraharp 400 was used to synchronize the arrival times of detected photons with the laser repetition rate. The overall instrument response function (IRF) of the system was measured to be 128 and 160 ps (FWHM) for 640 and 905 nm measurements, respectively.

Figure 8. (Left) Photos depicting our experimental setup. Snow was held in a cooler placed on the floor, and illuminated using pulsed diode lasers at two wavelengths (640 nm, 905 nm). A SPAD (enclosed in pink insulating foam) measured the time-of-flight of photons that exited the snow surface at distance s from the laser spot. (Center) Schematic of experimental setup. (Right) Time-of-flight histograms measured using our system. For this test the snow sample consisted of natural snow that had been aged for nine months at −10 $^\circ$C.

3.1.2. Measurement procedure

Experiments were conducted in a cold room at −1 $^\circ$C. The room's lights were switched off and windows blacked out to reduce interference from ambient background light. A large folding mirror was used to direct the lidar beam and detector field of view (FOV) towards a cooler filled with snow that was placed on the floor. A schematic that illustrates our system's optical design is shown in the center of Fig. 8. Because only a single laser diode could be operated at any one time, 640 nm measurements were collected first. During these measurements, a red bandpass filter (Edmund Optics TECHSPEC 650 nm/50 nm) was placed in front of the detector to suppress interference from ambient background light. Following these measurements the 640 nm laser head and bandpass filter were removed and replaced with the 905 nm laser head and a near-infrared bandpass filter (Thorlabs FL905-10). We then collected a second set of measurements.

The beam from either laser head could be scanned by hand using a set of steering mirrors. A lens was placed in front of the detector to focus its FOV to a small spot (<1 cm FWHM) on the snow surface. To find this focus spot, the laser beam would be steered to the point on the surface at which detector counts were maximized. Once the focus spot was found, a laser pointer (distinct from the pulsed diode lasers) was steered to mark the position of the focus spot. The pulsed beam could then be steered to a point on the snow surface that was displaced from the focus spot by a small distance s that was measured using a ruler. The focus-marker beam would then be switched off. When the 905  nm laser was in use, a phosphorescent laser viewing card was used to find the position of the laser spot on the snow surface.

We note that even when the laser and focus spots were separated by several centimeters, interference from direct returns off of the snow surface remained significant due to phenomena such as lens flare. Although we could not suppress this interference entirely, we were able to mitigate it by placing a long lens tube in front of our detector that functioned as a baffle. When possible, we would further reduce interference by placing a larger tube (5 cm diameter) constructed from blackout material (Thorlabs BKF12) on the snow surface, surrounding the spot observed by the detector. Together, these two baffles blocked most light paths that scattered off of the snow surface while transmitting paths that traveled beneath the surface.

For each snow sample, and for each laser wavelength, we collected measurements at multiple source-detector separations s. Each measurement consisted of a histogram of photon arrival times with 16 ps timing bins that spanned a 250 ns timing window. Examples of histograms measured by our lidar system are shown on the right in Fig. 8. The first measurement would always be collected at s = 0 cm to measure the time-of-arrival of photons that scattered directly off of the snow surface. The peak of this direct return would serve as a reference time for all subsequent measurements. Direct surface returns were always measured with a 60 second integration time, with a neutral density filter placed in front of the detector to prevent saturation, and with a wooden ruler placed on the snow surface at the position of the laser spot to prevent bias due to subsurface scattering. Following this, histograms would be collected for one or more non-zero source-detector separations. We used an integration time of 10 minutes for each histogram collected with 640 nm light, and 30 minutes for each histogram collected with 905 nm light. A longer integration time was required at 905 nm because our SPAD detector was less sensitive at this wavelength, the output power of our laser was lower, and the snow itself was less reflective. Ice grains from each sample were inspected before and after each set of measurements to ensure that snow properties had not changed significantly due to metamorphism.

Before proceeding, we want to stress that our lidar system was assembled strictly for the proof-of-principle demonstrations documented in this paper. It was not optimized for ease of use or light collection efficiency. Although the integration times reported here are quite long, we expect that a cleverly engineered system might collect equivalent data with integration times that are far shorter – perhaps by several orders of magnitude. Integration time could be reduced significantly, for instance, by using a multipixel Silicon Photomultiplier (SiPM) in place of the single-pixel SPAD used here, and by using lasers with higher power and higher repetition rates. The use of laser sources and SPADs designed for a consumer electronics environment (King and others, 2023), rather than the optical bench equipment used here, would also allow for a system that was portable, rugged, and affordable. Altogether, this suggests that the development of a field-deployable system is a feasible goal – one which we hope to pursue in future work.

3.2. Samples

We performed two sets of experiments. In the first, samples had relatively low LAP concentrations but grain size and density varied significantly. In the second, the samples had varying amounts of black carbon mixed into them, but density and grain size was relatively constant.

All snow used in our experiment originated as natural snow harvested on Dartmouth College campus and was subsequently modified in various ways. When not being used for experiments, snow samples were stored in lidded coolers in a −10 $^\circ$C cold room.

3.2.1. Clean snow samples

We performed five sets of measurements on samples with varying density and grain size but relatively low LAP content. The snow used in the first set of measurements was harvested after a snowfall in March 2022 and then kept in a −10 $^\circ$C cold room for nine months. By the time measurements were taken, the snow had become more dense and the grains had metamorphosed into medium size rounded grains and rounding faceted particles (Fierz and others, 2009). The next three data collections were performed on a single snow sample that was modified between measurements. The sample was harvested 30 minutes after snow had ceased falling, immediately outside our laboratory at Dartmouth College. It was then stored overnight at −10 $^\circ$C. Measurements were collected the next morning on the unmodified sample, which had a very low density and consisted of precipitation particles (Fierz and others, 2009), with many stellar dendrites. A second set of measurements was collected after the snow had been compacted with a shovel – thus increasing its density and potentially reducing grain size by fragmenting grains. The third set of measurements was collected after the snow was aged for three weeks at −10 $^\circ$C and then for one day at 0 $^\circ$C. This aging produced a clear change in grain shape, to small rounded grains and decomposing precipitation particles, and a small increase in grain size and density. For our final set of measurements we harvested snow that had been sitting outside for weeks, where it had experienced several melt and re-freeze events. This snow had very high density and coarse grains.

At the time of data collection, all samples were held in coolers with approximate internal dimensions of 50 cm × 25 cm × 30 cm and that had matte white internal walls. Snow would fill the cooler to varying degrees, but was typically at least 20 cm deep, relative to the cooler bottom.

3.2.2. Soot addition experiments

For the second set of experiments, we filled five Styrofoam coolers (dimensions 17.5 cm × 23.5 cm × 24.0 cm) with freshly fallen snow. We then mixed small amounts of Sigma-Aldrich Fullerene Soot (PN: 572497) into the samples, such that the five respective samples had 0, 1, 2, 3, and 4 baseline units of soot. To add soot to the snow in a controlled fashion, we created a soot-water suspension with a known concentration of soot, and then applied controlled volumes of the suspension to each snow sample with a spray bottle. The soot was mixed evenly into the snow using an ice scraper.

After performing a first set of measurements on the sooty samples we found that the added soot had a weaker effect on the snowpack absorption coefficients than had been expected. Following this finding, we approximately doubled the added soot concentration in all samples and repeated the measurements.

3.3. Ground truth measurements

Ground truth ice volume fraction was measured by extracting a small core with a depth of ~5 cm and diameter of ~6 cm from the snow surface using a cylindrical polypropylene jar. We measured the volume of snow in the core. The snow was then allowed to melt, and we measured the volume of the meltwater. Ice volume fraction was computed from the snow and meltwater volumes using conservation of mass.

Ground truth grain size was measured by imaging a small, snow-filled test tube (1.4 cm internal diameter) with a SkyScan 1172 microCT scanner (40 kV, 250 μA source, 17 μm resolution). Bruker NRecon software was used to reconstruct a 3D image of the sample. Following guidance from Hagenmuller and others (2016), the image was then blurred with a Gaussian kernel (radius 1 pixel), binarized with Otsu's method, and morphologically ‘opened’ (radius 1 pixel). The surface area to volume ratio (SA/V) of the imaged sample was then computed two times using Bruker's CTAN software, following marching squares (2D analysis) and marching cubes (3D analysis) surface reconstructions. We computed the grain radius from each SA/V ratio independently, and then used the average of these two values as ground truth.

Ground truth estimates of black carbon concentration were obtained using a single particle soot photometer (SP2; Droplet Measurement Technologies), in a manner similar to that reported in Lazarcik and others (2017). Each snow sample was melted and ultrasonicated for at least 15 minutes prior to analysis. The liquid snow samples were aerosolized using an ultrasonic nebulizer (CETAC U5000AT), which removes moisture from the liquid stream before passing aerosols such as black carbon onto the SP2. The SP2 estimates black carbon particle mass via measurements of laser-induced incandescence. This system was calibrated using a series of fullerene soot standards. To avoid saturating the SP2, snow samples that were expected to have particularly high black carbon concentrations were diluted with MilliQ water by a factor of 6.

4. Results

4.1. Clean snow experiments

4.1.1. Individual snow sample

To provide insight into our data collection and fitting procedures, we first present a detailed review of all measurements collected for a single snow sample. This sample, which is described in greater detail in the Materials section, consisted of natural snow that had been aged for nine months in a $-10^\circ$C cold room.

The raw, time-of-flight histogram data collected for this sample, as well as our curve fits to those measurements, are shown in Fig. 9. Measurements were taken at four different source-detector separations for each wavelength: s = 4, 6, 8 and 10 cm at 640 nm, and s = 4, 5, 6 and 7 cm at 905 nm. As a rule of thumb, we would start each fit at a timing bin that corresponded to the peak of the diffusion signal. This was done to avoid fitting to the earliest arriving photons, which are poorly described by our diffusion model. Histograms were collected at multiple s values because it was not known a priori what range of s values would yield good diffusion curve fits. If s and μ _s′ were both small, then photons in the signal peak would be poorly described by our diffusion model because they would exit the snowpack after too few scattering events. On the other hand, if s and μ _s′ or μ _a were too large, the diffusion signal would be faint relative to background interference, and the fit would be poor.

Figure 9. Raw measurements collected for a single snow sample. Time-of-flight histograms were measured at two wavelengths (640 nm, 905 nm) and for four source-detector separations per wavelength. A diffusion model was fit to each histogram. The pair of curves with the best goodness of fit was used to compute snowpack properties.

In Fig. 10, we show how the retrieved snow properties varied with respect to our choices of source-detector separation at each wavelength. In general, estimates of $v_\ast$, $r_\ast$, and C _bc did not vary substantially if good curve fits were obtained at both wavelengths, but diverged from the typical value when one or both of the curve fits were poor. As an example, it is evident in Fig. 10c that C _bc estimates are biased high at s ₆₄₀ = 4 cm, but are otherwise relatively insensitive to changes in s at either wavelength. A relatively small amount of variance was observed in snow property estimates even when good fits were obtained. The sources of this variance are unconfirmed, but could be explained by instrumental phenomena such as photon pile-up distortion (Coates, 1968); the varied influence of unmodelled phenomena such as surface roughness, finite cooler size, or interference from direct surface returns; the metamorphism of snow between data collections; or true spatial variation in the snow's properties.

Figure 10. Dependence of $v_\ast$, $r_\ast$, and C _bc estimates on choice of source-detector separation s for each measurement wavelength. Estimates that correspond to the pair of curve fits with the lowest reduced deviance (McCullagh, 2019) are outlined in red.

To arrive at a single estimate for $v_\ast$, $r_\ast$, and C _bc, we chose the curve fit at each wavelength with the lowest reduced deviance (McCullagh, 2019). Deviance is a goodness of fit metric that is appropriate for data that follows Poisson statistics, and that is asymptotically equivalent to χ ² goodness of fit when the number of counts in all histogram bins is high. For the data collection described here, the best fits corresponded to s = 10 cm at 640 nm and s = 7 cm at 905 nm. From the parameters of these two fits we estimated that $v_\ast = 0.361\pm 0.004$, $r_\ast = 379\pm 4\mu {\rm m}$, and C _bc = 91 ± 1 ppbw. As described previously, the reported uncertainties correspond to statistical uncertainties in the curve fit parameters, propagated through Eqs. (17), (18), and (19). They do not account for potential inaccuracies in the diffusion or scattering models. The ground truth measurements of $v_\ast$, $r_\ast$, and C _bc were 0.465, 242.5 μm, and 30.7 ppbw, respectively.

4.1.2. Full results

We now present a summary of all results obtained for the clean snow samples. The properties of the snow samples used in these tests varied widely, from light, fine-grained fresh powder to dense, coarse-grained snow that had experienced several melt and re-freeze events. In Fig. 11, we show a scatter plot of the densities and grain sizes estimated using our method, as well as ground truth values.

Figure 11. Summary of the estimated and ground truth ice volume fractions $v_\ast$ and grain sizes $r_\ast$ of five clean snow samples. Error bars indicate one standard deviation. Estimates and ground truth values are matched by color.

Our estimates of $v_\ast$, $r_\ast$, and C _bc are plotted with respect to ground truth in Fig. 12. In Fig. 12a, we see a clear positive and nearly linear relationship between the ice volume fraction estimated using our technique, and ground truth, although estimates appear to be biased towards lower densities. The trends for $r_\ast$ and C _bc are less clear, although our method appears to be capable of distinguishing between small and large grain sizes, and low and moderate impurity concentrations. To the extent that trends can be observed, there appears to be an approximately 1:1 relationship between $r_\ast$ estimates and ground truth, whereas C _bc appears to be over-estimated by a factor of ~2.5. All statistical uncertainties in $v_\ast$ and $r_\ast$ estimates are 1–2% of the estimated value. All statistical uncertainties in C _bc estimates are 1–2  ppbw. Notably, these uncertainties are comparable to the statistical uncertainties reported by Allgaier and others (2022) in their estimates of black carbon concentrations in glacier ice.

Figure 12. Ground truth versus estimated values of ice volume fraction $v_\ast$, grain size $r_\ast$, and black carbon mass mixing ratio C _bc for five clean snow samples. Blue marks indicate estimates obtained using two measurement wavelengths (640 nm, 905 nm). Error bars indicate one standard deviation. Red marks indicate estimates of $v_\ast$ and $r_\ast$ obtained from 905 nm measurements only. Uncertainties were not computed for 905-only estimates.

In addition to dual-wavelength estimates of $v_\ast$, $r_\ast$, and C _bc, we also show estimates of $v_\ast$ and $r_\ast$ computed using only 905 nm measurements. The single-wavelength results match the dual-wavelength results very closely. Ice volume fraction estimates are slightly higher, which is consistent with excess absorption due to unmodeled LAPs. Single-wavelength grain size estimates are alternately higher or lower than corresponding dual-wavelength estimates.

Considering the very small statistical uncertainties in our results, we expect that the biases seen here are most likely attributable to model mismatch. In particular, the excess black carbon content predicted by our method is plausibly explained by the presence of other kinds of light absorbing impurities such as dust. The samples used in this test were collected outdoors and were handled with shovels, ice scrapers, and various other equipment that may have been coated with dust or dirt. Further investigation is needed to understand the bias in estimates of $v_\ast$, which appear to be underestimated by a factor of approximately 3/4. One possible explanation is that the measured signal was influenced by unmodeled reflections from the white side-walls of the cooler. A deeper analysis would be required to confirm this. However, one would expect such reflections to reduce the observed decay rate by scattering photons back into the probed region instead of allowing them to escape. There is a notable outlier among the $r_\ast$ estimates that is approximately 50% higher than its ground truth measurement (estimate: 379.0 ± 4.1 μm, truth: 242.5 μm). The origins of this outlier are unclear. By inspection of Eq. (7), we see that the estimated $r_\ast$ value would be reduced by 50% if the modeled scattering asymmetry factor was increased from 0.825 to 0.883. It is thus possible that the outlier snow sample – which consisted of medium size rounded grains and rounding faceted particles that had been aged for nine months in a −10 $^\circ$C cold room – had a higher scattering asymmetry factor than the others. However it is not clear why this would be so.

4.2. Soot addition experiments

Here we present the results of the soot addition experiments described in the Materials section, where the snow samples contained varying concentrations of black carbon. For these tests, the source-detector separation was held fixed at s = 8 cm for 640 nm measurements. For 905 nm measurements, a value of s = 5  cm was typically used, although this was occasionally reduced to 4 cm if the measured signal would otherwise be too faint to yield a good curve fit.

The primary goal of these experiments was to assess the accuracy and sensitivity of the estimates of black carbon mass mixing ratio produced by our method. To this end, in Fig. 13 we show a plot of the C _bc values estimated with our method versus ground truth estimates obtained using an SP2. Blue data points correspond to the first set of measurements, for which the soot concentrations were relatively low, and red data points correspond to a second set of measurements that was collected after the added black carbon concentration in each snow sample had been approximately doubled.

Figure 13. Retrieved values of black carbon mass mixing ratio C _bc plotted with respect to ground truth estimates obtained by a single-particle soot photometer (SP2). Increasing quantities of Sigma-Aldrich fullerene soot were added to five snow samples. Results from the first set of measurements are shown in blue. Soot concentrations were then approximately doubled for all samples and a second set of measurements was taken. Results from the second set of measurements are shown in red. Error bars indicate one standard deviation. SP2 errors are typically too small to be visible.

Upon inspection we see a clear correlation between the estimated and ground truth C _bc values. The correlation is approximately linear and nearly one-to-one. Two outlier data points (with ground truth C _bc of 58, 59 ppbw) lie off of the one-to-one line. We expect that the outliers are the result of an error in the ground truth estimates. It is possible that our mixing process did not uniformly distribute the black carbon content throughout the snow and that the region sampled for SP2 analysis was unusually clean.

The range of estimated C _bc values indicates that our technique is sensitive to concentrations above 100 ppbw. Notably, this is significantly more sensitive than estimates derived from remote, multi-spectral albedo measurements, which are unreliable for black carbon concentrations below 1000 ppbw (Zege and others, 2011; Warren, 2013). Methods that infer black carbon concentration from in situ, hyperspectral albedo are sensitive to black carbon concentrations above 50 ppbw (Dumont and others, 2017), which is comparable to the sensitivity achieved here. The sensitivity of our method could be improved, perhaps significantly, by using a blue or green laser source in place of the red laser used here. As shown in Fig. 4, the influence of black carbon on the absorption coefficient of snow is much stronger at these wavelengths.

For good measure, we also show the estimates of ice volume fraction, grain radius, and black carbon concentration obtained for all samples. Estimates are plotted in Fig. 14 as a function of total units of soot-water solution that were applied to each sample using a spray bottle. Although ground truth $v_\ast$ and $r_\ast$ were not collected, results in Fig. 14a and b indicate that the density and grain size of the snow samples was relatively consistent, but did have some variance. This variance may have been caused by differences in how each sample was mixed, or from interaction with the liquid water in the soot-water suspension. Regardless, we see in Fig. 14c that estimated C _bc increases approximately linearly and nearly monotonically as a function of the amount of added soot, with no clear dependence on density or grain size.

Figure 14. Estimates of $v_\ast$, $r_\ast$, and C _bc obtained during soot addition experiments. Estimates are plotted as a function of the number of soot units mixed into the snow. One soot unit corresponds to a fixed volume of soot-water suspension that is applied to the snow sample with a spray bottle and then mixed into the snow volume. The amount of soot per unit in the second run (red) was approximately double the soot per unit in the first run (blue). Error bars indicate one standard deviation.

5. Discussion

In this work we have introduced a new method for measuring the density, grain size, and black carbon content of a dry snowpack using non-invasive, time-domain diffuse optical measurements. We have presented a model for the time-domain optical response of a snowpack that was adapted from the biomedical optics literature (Kienle and Patterson, 1997; Haskell and others, 1994). Our model was obtained by solving the photon diffusion equation – an approximation to the radiative transfer equation that accurately describes the propagation of light in highly scattering media (Welch and van Gemert, 1995). We used a geometric scattering model to relate the parameters of our photon diffusion model to dry snowpack density, grain size, and black carbon concentration. Our scattering model was derived from a well-known snow-optics model (Kokhanovsky and Zege, 2004), but extended to account for the effective speed of light within snow, as well as the mixing state of black carbon. We then developed an algorithm to retrieve the snowpack properties from time-domain optical measurements collected at two wavelengths.

We were able to validate our method in a series of proof-of-principle experiments in which we measured the properties of real snow samples using a photon-counting lidar system. The results of these experiments are encouraging. We see a clear, nearly linear correlation between the snowpack densities estimated by our method, and ground truth. When the LAPs in the snow were known to be black carbon particles, we observed a nearly one-to-one correlation between the black carbon mass mixing ratios estimated using our method, and those measured using a single-particle soot photometer (Schwarz and others, 2006). A nearly one-to-one correlation was also found between the grain sizes measured by our method and those determined from micro-CT images – although this correlation was not as strong. Our goal in this work was to obtain proof-of-principle results. More experiments are required to comprehensively assess our method's accuracy, biases, and failure modes.

Although our results are encouraging, we believe that the primary contribution of our work is not necessarily the exact method that we have proposed, but rather that we have been able to clearly demonstrate that time-domain diffuse optics is an appropriate sensing modality for measuring snowpack properties. Previous works have used ray-tracing simulations to explore the feasibility of inferring snow properties from time-domain diffuse optical signals (Libois and others, 2019), and to predict the relationship between snow properties and lidar altimetry biases (Smith and others, 2018). However, as far as we are aware, our work is the first to provide clear experimental evidence that the optical response of a snowpack that has been illuminated by a laser pulse can be accurately described using a photon diffusion model; and also that this response is measurably influenced by changes to important snowpack properties like grain size, density, and impurity content.

Our method could be improved in several ways. More sophisticated models that incorporate features such as liquid water content in the snow, finite snow depth, or surface roughness might be developed to enable retrieval of snow properties in a broader set of circumstances. Our measurement procedure could also be improved. In our experiments, processing a single snow sample required between 40 minutes to several hours of data collection time. This could be dramatically reduced by improving our instrument design to incorporate multi-pixel SPAD detectors, higher power lasers, or by simply placing the laser and detector closer to the snow surface. The integration times used in our experiments were also conservative – further analysis could determine the minimum number of photons required to accurately retrieve snow properties. Finally, using our method in the field would require the development of a rugged and portable instrument. The dramatic decrease in the cost and size of pulsed lasers and photon counting detectors in recent years makes this possible. All components required for such a device can be found in the current model of the iPhone Pro (King and others, 2023).

Although non-invasive, optical methods for measuring snow grain size (Nolin and Dozier, 2000; Gallet and others, 2009) and LAP concentrations (Dumont and others, 2017; Allgaier and Smith, 2022) have been demonstrated previously, our work provides what is, to our knowledge, the first experimental demonstration of snow density estimation from non-invasive optical measurements. Previous works used ray-tracing simulations to demonstrate non-invasive porosity ($= 1-v_\ast$) measurements in arbitrary porous media (Libois and others, 2019), or estimated snow density using invasive optical transmission measurements (Gergely and others, 2010). Our method could potentially provide field measurements of snow-water-equivalent (SWE) – the product of snow depth and density – or surface density, which might in turn prove useful in hydrological or ecological studies, or for validating remote sensing techniques (Kinar and Pomeroy, 2015). Diffuse optical methods may also enable more sensitive field measurements of LAPs, particularly if shorter wavelength (e.g. blue or green) illumination is used. Field measurement of black carbon concentration from snow's hyperspectral albedo has been demonstrated (Dumont and others, 2017). Our method infers black carbon concentration from a decay rate parameter that is proportional to snow's absorption coefficient, which is far more sensitive to impurity content than albedo. The trace concentrations of impurities found in remote snowpacks reduce snow's albedo by at most a few percent (Warren, 2013) whereas, in theory, the absorption coefficient of snow at green wavelengths should be doubled by just a few ppbw of black carbon (see Fig. 4). We note that the spatially-resolved diffuse optical technique demonstrated by Allgaier and Smith (2022) also infers LAP concentrations via the absorption coefficient, and so should be able to achieve comparable sensitivities.

Although the instrument used in this study was assembled from the same components that make up a typical photon counting lidar system, our measurements were effectively in situ because our lidar was always placed within a meter of the snow's surface. In the future, we hope to develop true remote lidar sensing techniques that are grounded in time-domain diffuse optics models. Such methods would enable important capabilities such as the remote mapping of SWE or impurity concentrations. However, the leap from in situ to remote measurements poses new challenges that include dramatically lower photon counts, wider beam footprints, and confocal measurement geometries. Further analysis is required to determine which snow properties can be feasibly retrieved under these conditions. An alternative direction for future work is the development of more advanced algorithms for processing in situ measurements that leverage decades of advances in diffuse optical spectroscopy research (Konugolu Venkata Sekar and others, 2019). In particular, the adaption of diffuse optical tomography methods (Okawa and Hoshi, 2023) to snow would enable non-invasive retrieval of snow stratigraphy, or even full 3D mapping of snowpack properties within a probed region. Observations of snow stratigraphy are often made to assess the structural stability of the snowpack to predict avalanche risk, as well as the history of snow deposition and metamorphism (Nienow and Campbell, 2011).

6. Conclusions

We have developed a new technique to estimate the density, grain size, and black carbon concentration of a dry snowpack using time-resolved measurements of laser light that has scattered multiple times beneath the snow surface. Our method was inspired by diffuse optical spectroscopy techniques that were originally developed for biomedical applications. We validated our method in a series of proof-of-principle experiments. Our results revealed strong, nearly linear correlations between our estimates of snow density and black carbon concentration, and independent ground truth measurements. Additionally, our method successfully distinguished between small and large grain sizes. Our results indicate that non-invasive optical measurement of snow density, grain size, and black carbon concentration is possible. However, further refinement of our instrument design is needed for practical field use. More broadly, our work provides the first clear experimental evidence that time-dependent scattering of laser light by snow is well described by a diffuse optical model. This could pave the way for future algorithms that retrieve snow properties from remote lidar measurements as well as more advanced in situ techniques, such as methods that infer snow stratigraphy.