# Measurement of Snowpack Density, Grain Size, and Black Carbon Concentration
Using Time-domain Diffuse Optics

Connor HENLEY  1,2 Joseph L. HOLLMANN  2 Colin R. MEYER  3 Ramesh RASKAR1
<EMAIL_ADDRESS>1MIT Media Lab, Massachusetts Institute of Technology, 75
Amherst Street, Cambridge, MA, USA
2The Charles Stark Draper Laboratory, Inc., 555 Technology Square, Cambridge,
MA, USA
3Thayer School of Engineering, Dartmouth College, 15 Thayer Drive, Hanover,
NH, USA
Correspondence: Connor Henley

###### Abstract

Diffuse optical spectroscopy (DOS) techniques aim to 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 snowpack
as a function of its density, grain size, and black carbon content, and we
develop an algorithm that retrieves these properties from measurements at two
wavelengths. To validate our approach, we use a two-wavelength lidar system
and 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-10cm). We find strong, linear correlations between our retrievals
of density and black carbon concentration, and ground truth measurements.
Although the correlation is not as strong, we also find that our method is
capable of distinguishing between small and large grains.

##  1 Introduction

Snow is composed of transparent ice grains that absorb light very weakly at
visible and near-infrared 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
intrinsic 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, which precludes their use at night and for
several months of the year in polar regions. Albedo models also assume steady-
state illumination that is either collimated or diffuse. As such, they cannot
fully model the optical response of snowpack to the focused, time-modulated
illumination employed by lidar systems.

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).

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 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. Fair and others, (2022) experimentally confirmed
the predictions of Smith and others, (2018) by comparing biases in the snow
surface heights retrieved by IceSat-2 and by the 532 nm and 1064 nm beams on
the Airborne Topographic Mapper. 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
when it is rich in air bubbles.

In this work we introduce what is, to our knowledge, the first 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 snowpack
density (through the ice volume fraction), grain size, and the concentration
of light absorbing particles, in a non-invasive way. As far as we are aware,
ours is the first method to estimate snowpack _density_ using non-invasive
optical reflectance measurements. Although our system uses the same components
as a photon-counting lidar, our measurements are effectively _in situ_ , as
the instrument is always placed within a meter of the snow surface. However,
by demonstrating that important bulk snowpack properties can be retrieved from
time-domain measurements of multiply scattered photons, we hope that our work
will motivate the future development of techniques that retrieve snowpack
properties from remote lidar measurements. We also believe that time-domain
diffuse optical measurements, in general, represent a new frontier for
studying the optical properties of snow.

##  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_
($\text{W}\cdot\text{m}^{\text{-2}}\cdot\text{s}r^{\text{-1}}$) within a
medium as a function of space and time. The scattering medium is described by
a scattering coefficient $\mu_{s}$ ($\text{m}^{\text{-1}}$), a scattering
phase function, an absorption coefficient $\mu_{a}$ ($\text{m}^{\text{-1}}$),
and the speed of light within the medium $c_{*}$ (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 (= $(\mu_{a}+\mu_{s})^{-1}$), and photons are
typically scattered many times before they are absorbed ($\mu_{s}\gg\mu_{a}$)
(Welch and van Gemert,, 1995). The photon diffusion equation can be written as
follows:

$\frac{1}{c_{*}}\frac{\partial}{\partial
t}\phi(\mathbf{r},t)-D\nabla^{2}\phi(\mathbf{r},t)+\mu_{a}\phi(\mathbf{r},t)=S(\mathbf{r},t).$
(1)

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_
$\phi(\mathbf{r},t)$ ($\text{W}\cdot\text{m}^{\text{-2}}$), which is the
integral of radiance over all directions. Here
$D=3\left[\mu_{a}+(1-g)\mu_{s}\right]^{-1}$ is the diffusion constant and $g$
is the anisotropy 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$) or forward scattering ($g>0$). The variable $S$
represents an isotropic source term.

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 at
distance $s$ from the laser spot. (b) Comparison of time-dependent intensity
predicted by our model (Eq. 4), 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.

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. 1(a). 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 small distance $s$. Kienle and Patterson,
(1997) showed that, in this scenario, Eq. 1 can be accurately solved by
imposing an an _extrapolated boundary condition_ , which requires that photon
fluence goes to zero along a planar boundary that lies just above the medium’s
surface. This yields the following expression for photon fluence inside the
medium:

$\phi(s,z,t)=\frac{c_{*}}{(4\pi
Dc_{*}t)^{3/2}}\exp{\left(-\mu_{a}ct\right)}\left\\{\exp{\left[-\frac{(z-z_{0})^{2}+s^{2}}{4Dc_{*}t}\right]}-\exp{\left[-\frac{(z+z_{0}+2z_{b})^{2}+s^{2}}{4Dc_{*}t}\right]}\right\\}.$
(2)

Here $z$ denotes the distance from the surface, going down;
$z_{0}=[\mu_{a}+(1-g)\mu_{s}]^{-1}$ is the transport mean free path of photons
within the medium; and $z_{b}$ denotes the height of the extrapolated
boundary. Haskell and others, (1994) proposed a value of
$z_{b}=\frac{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 time-dependent 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. 2, we compute the time dependent radiosity
($\text{W}\cdot\text{m}^{\text{-2}}$) that exits the surface at position $s$
using Fick’s Law (Kienle and Patterson,, 1997)

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

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:

$R(s,t)=\frac{\alpha
c_{*}}{3(2\pi)^{3/2}}\frac{z_{0}^{2}}{(2Dc_{*}t)^{5/2}}\exp\left(-\mu_{a}c_{*}t-\frac{s^{2}+z_{0}^{2}}{4Dc_{*}t}\right)\left[1+\frac{7}{3}\exp\left(-\frac{10z_{0}^{2}}{9Dc_{*}t}\right)\right],$
(4)

where $\alpha$ 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_{0}/3$ and
$z_{b}=2z_{0}/3$. In deriving Eq. 4, 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. 4.

In Fig. 1(b) we compare the time-dependent intensity predicted by Eq. 4 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_{*}t\gg z_{0}$) and larger distances from the
laser spot ($s\gg z_{0}$), as these photons have scattered many times before
exiting the medium.

###  2.2 Snow Scattering Model

Our measurement model, defined in Eq. 4, is expressed in terms of three
phenomenological parameters—the absorption coefficient $\mu_{a}$, the
effective scattering coefficient $\mu_{s}^{\prime}=(1-g)\mu_{s}$, and the
speed of light in the medium $c_{*}$. As is often done when modeling snow
albedo (Wiscombe and Warren,, 1980a), we assume the snowpack can be modeled as
a monodispersion of spherical ice grains. This assumption allows us to use a
Mie scattering model to define $\mu_{a}$, $\mu_{s}^{\prime}$, and $c_{*}$ in
terms of three physically meaningful snowpack parameters—$v_{*}$, the fraction
of the snowpack volume that is occupied by ice; $r_{*}$ (m), the grain radius;
and $C_{bc}$ (kg/kg), the mass mixing ratio of black carbon in the snowpack.
We note that the ice volume fraction $v_{*}$ is readily converted to snowpack
density $\rho_{*}$ ($kg/m^{3}$) via the expression
$\rho_{*}=v_{*}\rho_{ice}+(1-v_{*})\rho_{air}\approx v_{*}\rho_{ice}$, where
$\rho_{ice}$ and $\rho_{air}$ are the bulk densities of ice and air,
respectively.

####  2.2.1 Clean Snowpack

For a snowpack that contains optically insignificant concentrations of light
absorbing particles, the scattering and absorption coefficients can be written
entirely as functions of $v_{*}$ and $r_{*}$. We use Mie codes (Sumlin and
others,, 2018) to compute the scattering efficiency $Q_{s}$, the absorption
efficiency $Q_{a}$, and the scattering asymmetry factor $g$ as a function of
grain radius, wavelength, and the complex refractive index at the chosen
wavelength (Warren and Brandt,, 2008). The absorption and effective scattering
coefficients of the snowpack are then computed as follows:

$\mu_{a}=Q_{a}\pi r_{*}^{2}N_{*}=\frac{3Q_{a}v_{*}}{4r_{*}}$ (5)

$\mu^{\prime}_{s}=Q_{s}(1-g)\pi
r_{*}^{2}N_{*}=\frac{3Q_{s}(1-g)v_{*}}{4r_{*}},$ (6)

where $N_{*}$ is the number density of ice grains in the snowpack. In theory,
the values of $Q_{a}$, $Q_{s}$, and $g$ display a highly non-linear and
oscillatory dependence on the grain size $r_{*}$. This oscillatory “ringing”
is an artifact of the assumption that all grains have precisely the same
radius—a condition that is never satisfied by real snow. When the oscillations
are smoothed out, $Q_{s}$, $g$, and the _ratio_ $Q_{a}/r_{*}$ become nearly
independent of grain size at visible and near-infrared wavelengths, so long as
the grains are not unreasonably large. We thus treat these values as
wavelength-dependent constants, computed by taking their average value across
a range of feasible grain sizes at the wavelength of interest. Doing so
dramatically simplifies our model—$\mu_{a}$ becomes independent of grain size
and depends linearly on $v_{*}$, and $\mu_{s}^{\prime}$ depends only on the
ratio $v_{*}/r_{*}$.

Figure 2: (a) Effective scattering coefficient $\mu_{s}^{\prime}$ (1/m) as a
function of ice volume fraction $v_{*}$ (unitless) and grain radius $r_{*}$
(mm). Values computed for 640 nm light, such that $Q_{s}=2.0053$ and
$g=0.8926$ (Sumlin and others,, 2018). (b) Absorption coefficient $\mu_{a}$
(1/m) of clean snow as a function of ice volume fraction $v_{*}$ (unitless)
and wavelength $\lambda$ (nm).

In Fig. 2(a) we visualize the range of values for $\mu_{s}^{\prime}$ obtained
across a domain of feasible grain sizes and ice volume fractions. The values
shown here were computed for a wavelength of 640 nm—however, because $Q_{s}$
and $g$ are only weakly dependent on wavelength, the $\mu_{s}^{\prime}$ values
shown in Fig. 2(a) are comparable to those that would be obtained across a
large fraction of the visible and near-infrared spectrum. Figure 2(b) shows
the dependence of $\mu_{a}$, obtained using Eq. 5, on ice volume fraction and
wavelength. Unlike $\mu_{s}^{\prime}$, the absorption coefficient depends
strongly on wavelength, and varies by more than an order of magnitude between
the red ($\lambda=$ 640 nm) and near infrared ($\lambda=$ 905 nm) wavelengths
used in this study.

The last parameter that we must calculate is the speed of light within the
snowpack $c_{*}$. 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 heterogenous mixture of two materials—ice and
air—that have markedly different refractive indices and that may be mixed at
any ratio.

We follow the lead of Smith and others, (2018) and use the _mean_ speed of
light within the snowpack, which can be derived from the expected fraction of
a photon’s path that passes through ice, and through air. This value is
necessarily dependent upon the ice volume fraction, and can be written as

$c_{*}=\frac{c_{0}}{1+(n_{ice}-1)v_{*}},$ (7)

where $n_{ice}$ is the real component of the index of refraction of ice and
$c_{0}$ is the speed of light in air (where it’s assumed that $n=1$).

Having obtained expressions that relate $\mu_{a}$, $\mu_{s}^{\prime}$, and
$c_{*}$ to the grain size and ice volume fraction of a clean snowpack, we can
now develop an understanding of how changes to $v_{*}$ and $r_{*}$ affect the
snowpack’s time-domain optical response. Upon inspection of Eq. 4, we see that
the shape of a snowpack’s transient response is primarily controlled by
$\mu_{a}c_{*}$, which determines the rate of decay of the signal’s tail;
$2Dc_{*}$, 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\gg
z_{0}$.

The exponential decay rate, $\mu_{a}c_{*}$, only depends on $v_{*}$. On the
other hand, $2Dc_{*}$ and $z_{0}^{2}$ primarily depend on the medium’s
scattering coefficient, which in turn depends on the ratio $v_{*}/r_{*}$.
These effects are visualized in Figs 3(a) and (b), where we plot the predicted
transient response curves for snowpack with varying $v_{*}$ and $r_{*}$. 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 3: Plots of predicted time-dependent flux measured by a detector
observing snow following illumination by a laser pulse ($\lambda=$ 640 nm).
Curves were produced using Eq. 4 and $\mu_{a}$, $\mu_{s}^{\prime}$, and
$c_{*}$ were computed from $v_{*}$, $r_{*}$, and $C_{bc}$ using Eqs. 9, 6, and
7, respectively. (a) Ice volume fraction $v_{*}$ is varied. $r_{*}$ = 100
$\upmu$m, $C_{bc}$ = 0 ppbw, and $s$ = 8 cm. (b) Grain radius $r_{*}$ is
varied. $v_{*}$ = 0.3, $C_{bc}$= 0 ppbw, and $s$ = 8 cm. (c) Impurity
concentration $C_{bc}$ is varied. $v_{*}$ = 0.3, $r_{*}$ = 100 $\upmu$m, and
$s$ = 8 cm. (d) Detector focus position $s$ is varied. $v_{*}$ = 0.3, $r_{*}$
= 100 $\upmu$m, and $C_{bc}$ = 0 ppbw.

In Fig. $\ref{fig:var}(a)$, we see that as $v_{*}$ is increased while $r_{*}$
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. 3(b), $r_{*}$ is varied
while $v_{*}$ is held constant. 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. 3(d), $v_{*}$ and $r_{*}$ 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.2.2 The Effect of Light Absorbing Impurities

Ice is an exceptionally weak absorber of light, particularly at visible
wavelengths (Warren,, 2019). As such, the absorption of 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. 5, 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.

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}$. This additional term can be
computed from the dispersed (rather than bulk) density $\rho_{bc}$ and the
wavelength-dependent mass absorption efficiency $\textrm{MAE}_{bc}$
($\text{m}^{\text{2}}/\text{kg}$) of the black carbon particles (Grenfell and
others,, 2011), as follows:

$\mu_{a}^{bc}=\textrm{MAE}_{bc}\rho_{bc}=\textrm{MAE}_{bc}\rho_{ice}C_{bc}v_{*}.$
(8)

On the right side of Eq 8, we replace $\rho_{bc}$ with the product of the bulk
density of ice ($\rho_{ice}=$ 916.5 kg/$\text{m}^{\text{3}}$), the mass mixing
ratio of black carbon in the snow $C_{bc}$, and the snowpack’s ice volume
fraction $v_{*}$. We then combine Eqs. 5 and 8 to obtain a complete expression
for the snowpack absorption coefficient:

$\mu_{a}=\left(\frac{3Q_{a}}{4r_{*}}+\textrm{MAE}_{bc}\rho_{ice}C_{bc}\right)v_{*}.$
(9)

Following the example of Doherty and others, (2014), we model the wavelength
dependence of $\textrm{MAE}_{bc}$ using a power law spectrum:

$\textrm{MAE}_{bc}(\lambda)=\textrm{MAE}_{bc}(\lambda_{ref})\left(\lambda_{ref}/\lambda\right)^{\textrm{{\AA}}},$
(10)

that has an Ångstrom coefficient Å $=1.1$ and is referenced to
$\textrm{MAE}_{bc}(\lambda_{ref})=6500$ $\text{m}^{\text{2}}/\text{kg}$, where
$\lambda_{ref}=600$ nm.

Figure 4: The ratio of light absorption due to black carbon (Eq. 8) to total
absorption by the snowpack (Eq. 9) for a range of wavelengths.

Fig. 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. 8) to the total absorption (from Eq. 9) for a
selection of wavelengths that ranges from 400 nm (blue) to 1000 nm (near
infrared). For blue light, our model suggests that absorption is entirely
dominated by just 1 ng/g 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 ng/g—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.

In Fig. 3(c), we show how changes to $C_{bc}$ affect that snowpack’s time-
domain optical response. We plot transient response curves for snow samples
with a variety of black carbon concentrations. These curves assume a laser
wavelength of 640 nm and a fixed grain radius, ice volume fraction, and focus
spot position. Black carbon content only influences the absorption coefficient
of the snowpack, and the effect of this is to steepen the exponential decay
rate $\mu_{a}c_{*}$. At the probing wavelength of 640 nm this effect is quite
dramatic when compared to the comparable influence of ice volume fraction on
the exponential decay rate, which is shown in in Fig. 3(a).

####  2.2.3 Non-spherical Grains

Although the assumption of spherical ice grains has been used successfully for
many applications (Wiscombe and Warren,, 1980a; Nolin and Dozier,, 2000), it
is not strictly correct as real snow often contains strikingly non-spherical
grain shapes (Bentley and Humphreys,, 1931). Scattering models for non-
spherical ice grains have been developed (Kokhanovsky and Zege,, 2004; Flanner
and others,, 2021). Non-spherical grain models introduce two additional
degrees of freedom, the scattering asymmetry factor $g$ and an “absorption
enhancement parameter” $B$ that accounts for internal reflections within the
ice grains (Kokhanovsky and Zege,, 2004).

We elected to use a spherical grain model because we assumed that these
additional degrees of freedom would prevent direct retrieval of $v_{*}$,
$r_{*}$, or $C_{bc}$. Recently, (Robledano and others,, 2023) argued that
while $g$ and $B$ vary significantly among idealized grain shapes such as
spheres, cubes and fractals, $g$ and $B$ for natural snow samples do not vary
very much, and cluster around the predicted values for a random two-phase
medium. Such findings suggest that it should be possible to both accurately
account for grain asphericity and unambiguously retrieve $v_{*}$, $r_{*}$, or
$C_{bc}$. In future work, we will apply the findings of (Robledano and
others,, 2023) to our inversions and determine the effect of asphericity on
the results.

###  2.3 Algorithm

We fit functions of the same form as Eq. 4 to two photon time-of-flight
histograms—each measured using a different laser wavelength. We use a grid
search to find the fit parameters that minimize a negative Poisson log-
likehilood function that properly accounts for photon count statistics. The
parameters of the fitted curves are then used to compute the snowpack
properties $v_{*}$, $r_{*}$, and $C_{bc}$. A visualization of our retrieval
algorithm is shown in Fig. 5.

Figure 5: (Left) A diffusion model (lower right) is fit to photon time-of-
flight histograms measured at two wavelengths. Fit parameters
$\alpha^{\prime}$, $\beta$, $\gamma$, and $\delta$ are determined using a grid
search algorithm that minimizes a negative Poisson log-likelihood function.
(Top right) Snowpack properties $v_{*}$, $r_{*}$, and $C_{bc}$ are computed
directly from parameters $\beta_{1}$, $\beta_{2}$, $\gamma_{1}$, and
$\gamma_{2}$ by evaluating analytical expressions (Eqs. 15, 16, and 17)

####  2.3.1 Fit Parameterization

We re-parameterize Eq. 4 in terms of the fitting parameters
$\alpha^{\prime}=\dfrac{\alpha c_{*}}{3(2\pi)^{3/2}}$, $\beta=\mu_{a}c_{*}$,
$\gamma=2Dc_{*}$, and $\delta=z_{0}^{2}$. This allows for the model to be
expressed in simplified form:

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

Although Eq. 11 appears to have four degrees of freedom, only the exponential
decay rate parameter $\beta$ and the spatial spread rate $\gamma$ are used to
estimate snowpack properties in practice. Interpreting the scaling constant
$\alpha^{\prime}$ requires precise calibration of the instrument and
measurement geometry which is challenging in practice and which we did not
attempt. Eq. 11 also depends only weakly on the squared source depth $\delta$,
to the point that accurate estimates of $\delta$ are almost never obtained.

The spatial spread rate $\gamma$ depends primarily on $\mu_{s}^{\prime}$ and
$c_{*}$, which are largely independent of the measurement wavelength for
visible and near-infrared light. 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_{*}$ and $r_{*}$ can be retrieved
from measurements at a single wavelength. Otherwise, retrieving $v_{*}$,
$r_{*}$, and $C_{bc}$ requires measurements at least two 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

$x_{i}=R(s,t_{i};\mathbf{\Theta})+b=\alpha^{\prime}r_{i}+b,$ (12)

where $R$ is the flux predicted by Eq. 11 at position $s$, time $t_{i}$, and
for parameters $\mathbf{\Theta}=\\{\alpha^{\prime},\beta,\gamma,\delta\\}$.
The rate of background counts produced by ambient light, detector dark counts,
and detector afterpulsing (Zappa and others,, 2007) is denoted by $b$, and is
assumed to be constant with respect to time. The variable $r_{i}$ denotes the
normalized predicted flux, for which $\alpha^{\prime}=1$.

The probability of observing a vector of time-binned measurements $\mathbf{y}$
given a vector of predicted count rates $\mathbf{x}$ is

$P(\mathbf{y}|\mathbf{x})=\prod^{N}_{\Delta}\frac{x_{i}^{y_{i}}e^{-x_{i}}}{y_{i}!}$
(13)

where $N$ denotes the total number of timing bins in the histogram and
$\Delta$ is the starting bin for the curve fit. We seek to find the parameters
$\alpha^{\prime}$, $\beta$, $\gamma$, $\delta$, and $b$ that minimize the
negative Poisson likelihood:

$\mathcal{L}(\mathbf{x},\mathbf{\Theta},b|\mathbf{y})=-\ln{P(\mathbf{y}|\mathbf{x})}=\sum^{N}_{\Delta}x_{i}-y_{i}\ln{x_{i}}+\ln{y_{i}!}.$
(14)

We find the set of parameters that minimizes the negative log likelihood using
a grid search. To reduce the dimensionality of the search, we first estimate
$b$ 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 $\beta$, $\gamma$, and $\delta$, the scaling term
$\alpha^{\prime}$ can then be estimated using the expression
$\bar{\alpha}^{\prime}=\sum^{N}_{\Delta}\left(y_{i}-b\right)/\sum^{N}_{\Delta}r_{i}$.

We can thus define a three-dimensional search area that contains all feasible
values of $\beta$, $\gamma$, and $\delta$. The feasible range for $\delta$ is
tightly constrained to
$\left(\dfrac{3\gamma}{2c_{0}}\right)^{2}<\delta<\left(\dfrac{3n_{ice}\gamma}{2c_{0}}\right)^{2}$,
which allows for a coarse fit to be obtained using an effectively two-
dimensional search.

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. 5. We estimated the uncertainty in the retrieved
values of $\beta$, $\gamma$, and $\delta$ by computing the inverse of the
Hessian of the loss function at the estimated minimum, and then taking the
diagonal terms.

####  2.3.3 Computing $v_{*}$, $r_{*}$, and $C_{bc}$

When measurements are obtained at two wavelengths, $\lambda_{1}$ and
$\lambda_{2}$, the ice volume fraction $v_{*}$ and black carbon mixing ratio
$C_{bc}$ can be extracted from the decay parameters $\beta_{1}$ and
$\beta_{2}$. Each term $\beta_{i}$ can be expressed as a function of $v_{*}$
and $C_{bc}$ by taking the product of Eqs. 7 and 9. This results in a set of
two equations which can be solved, first, for $v_{*}$:

$v_{*}=\frac{b_{2}\beta_{1}-b_{1}\beta_{2}}{c_{0}\left(a_{1}b_{2}-a_{2}b_{1}\right)-d_{1}b_{2}\beta_{1}+d_{2}b_{1}\beta_{2}}.$
(15)

For notational simplicity we have made the substitutions
$a_{i}=\frac{3}{4}\left[\frac{Q_{a}}{r_{*}}\right]_{i}$,
$b_{i}=\rho_{ice}\textrm{MAE}_{bc}(\lambda_{i})$, and
$d_{i}=n_{ice}(\lambda_{i})-1$. As was justified previously, we treat the
ratio $\left[\frac{Q_{a}}{r_{*}}\right]$ as a constant with respect to $r_{*}$
that depends only on wavelength. As before, the term $c_{0}$ refers to the
speed of light in air.

Once $v_{*}$ has been obtained, $C_{bc}$ can be computed as follows:

$C_{bc}=\frac{1}{c_{0}b_{1}}\left[\left(\frac{1}{v_{*}}+d_{1}\right)\beta_{1}-c_{0}a_{1}\right]=\frac{1}{c_{0}b_{2}}\left[\left(\frac{1}{v_{*}}+d_{2}\right)\beta_{2}-c_{0}a_{2}\right].$
(16)

After $v_{*}$ and $C_{bc}$ have been computed, the grain radius $r_{*}$ can be
computed from the spatial spread parameter $\gamma_{i}$ at either wavelength:

$r_{*}=e_{i}\left[\frac{2c_{0}}{3\gamma_{i}v_{*}(1+d_{i}v_{*})}-a_{i}-b_{i}C_{bc}\right]^{-1}.$
(17)

Here $a_{i}$, $b_{i}$, and $d_{i}$ are defined as they were previously, and
$e_{i}=3Q_{s,i}(1-g_{i})/4$. Because $r_{*}$ can be computed using either
$\gamma_{1}$ or $\gamma_{2}$, we evaluate Eq. 17 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_{*}$. Uncertainties in $v_{*}$,
$r_{*}$, and $C_{bc}$ are obtained via error propagation from uncertainties in
$\beta_{1}$, $\beta_{2}$, $\gamma_{1}$, and $\gamma_{2}$.

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

$v_{*}=\frac{\beta}{ac_{0}-\beta d},$ (18)

and then $r_{*}$ can be obtained from the spatial spread rate $\gamma$, and
our estimate of $v_{*}$:

$r_{*}=e\left[\frac{2c_{0}}{3\gamma v_{*}(1+dv_{*})}-a\right]^{-1}.$ (19)

Here $a$, $b$, $d$, and $e$ retain their meanings from Eqs. 15 and 17.

####  2.3.4 Evaluation Using Simulated Measurements

Figure 6: 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.

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 $\mu_{a}$, $\mu_{s}^{\prime}$, and $c_{*}$ were
computed from $v_{*}$, $r_{*}$, and $C_{bc}$ using Eqs. 6, 7, and 9. To
simulate pencil-beam illumination, photons were launched at the origin
($[x,y,z]=\mathbf{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. For 640 nm measurements, photons
detected at $s=10.0\pm 0.5$ cm were used for curve fitting, whereas for 905 nm
measurements photons detected at at $s=7.0\pm 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 $b$ that was chosen to be consistent with the uniform background
count levels observed in experimental measurements.

Our results are shown in Fig. 6. For the first simulation, the true snowpack
properties were $v_{*}$ = 0.465, $r_{*}$ = 240 $\upmu$m, and $C_{bc}$ = 50
ppbw. Our method retrieved values of $v_{*}$ = 0.445$\pm$0.011, $r_{*}$ =
235.2$\pm$5.5 $\upmu$m, and $C_{bc}$ = 52.0$\pm$3.0 ppbw. For the second
simulation, the true snowpack properties were $v_{*}$ = 0.162, $r_{*}$ = 85
$\upmu$m, and $C_{bc}$ = 0 ppbw. Our method retrieved values of $v_{*}$ =
0.162$\pm$0.002, $r_{*}$ = 86.0$\pm$0.9 $\upmu$m, and $C_{bc}$ = 1.3$\pm$2.1
ppbw. These results suggest that our algorithm should produce accurate
estimates of snow properties if our snow scattering model is correct.

##  3 Materials

###  3.1 Apparatus

####  3.1.1 Lidar System

Figure 7: (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 single-photon avalanche diode measured the
time-of-flight of photons that exited the snow surface at distance $s$ from
the laser spot. (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$ ∘C.

We assembled a simple lidar system to measure the time-domain optical response
of a variety of snow samples. Photographs of our experimental setup are shown
on the left in Fig 7. Our lidar system used a single-pixel SPAD detector
(Microphoton Devices PDM series) with a timing jitter of $\sim$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 $\upmu$W, and the 905 nm laser was
operated at a time-averaged power of 55 $\upmu$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.

####  3.1.2 Measurement Procedure

Experiments were conducted in a cold room at $-1$ ∘C. 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. 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 650nm/50nm) 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 steering
mirror. 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.

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. 7. 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.

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$ ∘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$ ∘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. Measurements were first collected immediately
after the snow had fallen, when the snow 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 it’s density but leaving grain size
and shape relatively unchanged. The third set of measurements was collected
after the snow was aged for three weeks at $-10$ ∘C and then for one day at
$0$ ∘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$\times$25 cm$\times$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$\times$23.5 cm$\times$24.0 cm) with freshly fallen snow.
We then mixed small amounts of Sigma-Aldrich Fullerene Soot 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
(depth $\sim$5 cm) from the snow surface. 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
$\upmu$A source, 17 $\upmu$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 Fresh Snow Sample

To provide insight into our data collction 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$ ∘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. 8. 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 $\mu_{s}^{\prime}$ 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
$\mu_{s}^{\prime}$ or $\mu_{a}$ were too large, the diffusion signal would be
faint relative to background interference, and the fit would be poor.

Figure 8: 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. 9, 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_{*}$, $r_{*}$, and $C_{bc}$ were relatively stable so long as
good curve fits were obtained at both wavelengths, but diverged from the
stable value when one or both of the curve fits were poor. As an example, it
is evident in Fig. 9(c) that $C_{bc}$ estimates are biased high at $s$ = 4 cm,
but are otherwise relatively insensitive to changes in $s$ at either
wavelength.

To arrive at a single estimate for $v_{*}$, $r_{*}$, 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 $\chi^{2}$ 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_{*}$ = 0.370$\pm$0.003,
$r_{*}$ = 181.9$\pm$1.4 $\upmu$m, and $C_{bc}$ = 84.9$\pm$1.2 ppbw. The ground
truth measurements of these properties were 0.465, 242.5 $\upmu$m, and 30.7
ppbw, respectively.

Figure 9: Dependence of $v_{*}$, $r_{*}$, 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.

####  4.1.2 Full Results

We now present a summary of all results obtained for the clean snow samples
described in the Materials section. 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. 10, we show a scatter plot of the densities and grain sizes estimated
using our method, as well as ground truth values.

Figure 10: Summary of the estimated and ground truth ice volume fractions
$v_{*}$ and grain sizes $r_{*}$ of five clean snow samples. Error bars
indicate one standard deviation. Figure 11: Ground truth versus estimated
values of ice volume fraction $v_{*}$, grain size $r_{*}$, 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_{*}$
and $r_{*}$ obtained from 905 nm measurements only. Uncertainties were not
computed for 905-only estimates.

Our estimates of $v_{*}$, $r_{*}$, and $C_{bc}$ are plotted with respect to
ground truth in Fig. 11. In Fig. 11(a), 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_{*}$ 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, $r_{*}$ appears to be underestimated by a factor of
$\sim$2, whereas $C_{bc}$ appears to be over-estimated by a factor of
$\sim$2.5.

In addition to dual-wavelength estimates of $v_{*}$, $r_{*}$, and $C_{bc}$, we
also show estimates of $v_{*}$ and $r_{*}$ computed using only 905 nm
measurements. Notably, 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. The smaller than expected grain size
estimates are partly explained by the smaller than expected ice volume
fraction estimates—in our method, $v_{*}$ is estimated first and used as an
input to Eq. 17 or Eq. 19 to compute $r_{*}$. Even accounting for this, grain
size estimates are still lower than expected. The remaining bias is plausibly
explained by the fact that our model assumes a scattering asymmetry factor $g$
that is appropriate for spherical grains, but that is higher than than the
theoretical $g$ value for the most common non-spherical grain shapes including
spheroids, hexagonal plates, and fractals (Libois and others,, 2013). By
inspection of Eq. 6, we see that when $g$ is assumed to be higher than the
true value, predictions of $r_{*}$ will be lower than the true value. Further
investigation is needed to understand the biases in estimates of $v_{*}$. They
are likely not explained by aspherical grain shapes, as the inclusion of a
larger “absorption enhancement factor” that appropriately models common non-
spherical grain shapes would have increased, not decreased, estimates of
$v_{*}$. The effect of “close packing” of grains has been reported to reduce
the absorption enhancement factor of snow (Libois and others,, 2014). However,
it is not clear that the snow samples used in this study were close-packed to
an extraordinary degree.

###  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.

Figure 12: 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.

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. 12 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.

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.

Figure 13: Estimates of $v_{*}$, $r_{*}$, and $C_{bc}$ obtained for “double
soot” samples. 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. Error bars indicate one standard deviation.

For good measure, we also show the estimates of ice volume fraction, grain
radius, and black carbon concentration obtained for the five “double-soot”
samples described previously. Estimates are plotted in Fig. 13 as a function
of total units of soot-water solution that were applied to each sample using a
spray bottle. Although ground truth $v_{*}$ and $r_{*}$ were not collected,
results in Fig. 13(a) 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. 13(c) 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.

##  5 Discussion

In this work we have introduced a new method for measuring the density, grain
size, and black carbon content of a 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. 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. We used a Mie scattering
model to relate the parameters of our photon diffusion model to snowpack
density, grain size, and black carbon concentration, and we developed an
algorithm to retrieve these 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 mixing ratios estimated using our method,
and those measured using a single-particle soot photometer. A 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.
It is clear from our measurements that the optical response of a snowpack that
has been illuminated by a laser pulse can be accurately described using a
photon diffusion model. Furthermore, it is clear 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 many ways. The biases observed in our results
suggest that the use of more sophisticated scattering optics models should be
explored. Such models might adopt known features of current snowpack albedo
models, such as non-spherical grain shapes, but could also consider new
physics, such as the effect that internal reflections have on the mean speed
of light within the snow. 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 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 snow-water equivalent or impurity concentrations. However
the remote sensing environment introduces new challenges that include
dramatically lower photon counts, wider beam footprints, and confocal
measurement geometries. 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. In
particular, the adaption of _diffuse optical tomography_ methods to snow would
enable non-invasive retrieval of snow stratigraphy, or even full 3D mapping of
snowpack properties within a probed region. Finally, time-domain optical
measurements represent a new opportunity to develop a deeper understanding of
the scattering optics of snow. Our work has already shown that, unlike
spectral albedo, time-domain diffuse optical signals depend strongly on
snowpack density. It’s possible that more sophisticated measurements, such as
the time-resolved measurement of light penetration or bidirectional
reflectance, might be sensitive to other snowpack properties—for example
scattering asymmetry, or the absorption enhancement factor—that have been
challenging to measure by other means.

##  6 Acknowledgements

We thank Brent Minchew and Joanna Millstein for the helpful discussions that
inspired this project and for support with initial theoretical work. Andrii
Murdza provided general assistance with micro-CT measurements and cold room
experiments. Anna Valentine prepared the snow samples used in our soot
addition experiments. Jacob Chalif ran our SP2 measurements and helped us
interpret the results. Connor Henley was supported by a Draper Scholarship and
by a grant from the Office of Naval Research (N00014-21-C-1040). Any opinions,
findings and conclusions or recommendations expressed in this material are
those of the author(s) and do not necessarily reflect the views of the Office
of Naval Research. Colin R. Meyer was supported by the Heising-Simons
Foundation (#2020-1911), the Army Research Office (78811EG), and the National
Science Foundation (2024132).

## References

  * Ackermann and others, (2006) Ackermann M and others (2006) Optical properties of deep glacial ice at the south pole. _Journal of Geophysical Research: Atmospheres_ , 111(D13) (https://doi.org/10.1029/2005JD006687)
  * Allgaier and Smith, (2022) Allgaier M and Smith BJ (2022) Smartphone-based measurements of the optical properties of snow. _Appl. Opt._ , 61(15), 4429–4436 (10.1364/AO.457976)
  * Allgaier and others, (2022) Allgaier M, Cooper MG, Carlson AE, Cooley SW, Ryan JC and Smith BJ (2022) Direct measurement of optical properties of glacier ice using a photon-counting diffuse lidar. _Journal of Glaciology_ , 68(272), 1210–1220 (10.1017/jog.2022.34)
  * Bentley and Humphreys, (1931) Bentley W and Humphreys W (1931) _Snow Crystals_. McGraw-Hill, New York, NY
  * Doherty and others, (2014) Doherty SJ, Dang C, Hegg DA, Zhang R and Warren SG (2014) Black carbon and other light-absorbing particles in snow of central north america. _Journal of Geophysical Research: Atmospheres_ , 119(22), 12,807–12,831 (https://doi.org/10.1002/2014JD022350)
  * Durduran and others, (2010) Durduran T, Choe R, Baker W and Yodh A (2010) Diffuse optics for tissue monitoring and tomography. _Reports on progress in physics_ , 73(7), 076701
  * Fair and others, (2022) Fair Z, Flanner M, Vuyovich C, Adam MS and others (2022) Quantifying volumetric scattering bias in icesat-2 and operation icebridge altimetry over snow-covered surfaces. _Authorea Preprints_
  * Fierz and others, (2009) Fierz C, Armstrong R, Durand Y, Etchevers P, Greene E, McClung D, Nishimura K, Satyawali P and Sokratov S (2009) The international classification for seasonal snow on the ground, the international classification for seasonal snow on the ground, ihp-vii technical documents in hydrology no 83, iacs contribution no 1
  * Flanner and others, (2021) Flanner MG, Arnheim JB, Cook JM, Dang C, He C, Huang X, Singh D, Skiles SM, Whicker CA and Zender CS (2021) Snicar-adv3: a community tool for modeling spectral snow albedo. _Geoscientific Model Development_ , 14(12), 7673–7704 (10.5194/gmd-14-7673-2021)
  * Gallet and others, (2009) Gallet JC, Domine F, Zender CS and Picard G (2009) Measurement of the specific surface area of snow using infrared reflectance in an integrating sphere at 1310 and 1550 nm. _The Cryosphere_ , 3(2), 167–182 (10.5194/tc-3-167-2009)
  * Grenfell and others, (2011) Grenfell TC, Doherty SJ, Clarke AD and Warren SG (2011) Light absorption from particulate impurities in snow and ice determined by spectrophotometric analysis of filters. _Appl. Opt._ , 50(14), 2037–2048 (10.1364/AO.50.002037)
  * Hagenmuller and others, (2016) Hagenmuller P, Matzl M, Chambon G and Schneebeli M (2016) Sensitivity of snow density and specific surface area measured by microtomography to different image processing algorithms. _The Cryosphere_ , 10(3), 1039–1054 (10.5194/tc-10-1039-2016)
  * Haskell and others, (1994) Haskell RC, Svaasand LO, Tsay TT, Feng TC, McAdams MS and Tromberg BJ (1994) Boundary conditions for the diffusion equation in radiative transfer. _J. Opt. Soc. Am. A_ , 11(10), 2727–2741 (10.1364/JOSAA.11.002727)
  * Henderson and others, (2018) Henderson GR, Peings Y, Furtado JC and Kushner PJ (2018) Snow–atmosphere coupling in the northern hemisphere. _Nature Climate Change_ , 8(11), 954–963 (10.1038/s41558-018-0295-6)
  * Henley, (2020) Henley C (2020) Snow lidar monte-carlo simulator. Available at https://github.com/co24401/SnowLiDARMonteCarlo
  * Kienle and Patterson, (1997) Kienle A and Patterson MS (1997) Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium. _J. Opt. Soc. Am. A_ , 14(1), 246–254 (10.1364/JOSAA.14.000246)
  * King and others, (2023) King F, Kelly R and Fletcher CG (2023) New opportunities for low-cost lidar-derived snow depth estimates from a consumer drone-mounted smartphone. _Cold Regions Science and Technology_ , 207, 103757, ISSN 0165-232X (https://doi.org/10.1016/j.coldregions.2022.103757)
  * Kokhanovsky and Zege, (2004) Kokhanovsky AA and Zege EP (2004) Scattering optics of snow. _Appl. Opt._ , 43(7), 1589–1602 (10.1364/AO.43.001589)
  * Lazarcik and others, (2017) Lazarcik J, Dibb JE, Adolph AC, Amante JM, Wake CP, Scheuer E, Mineau MM and Albert MR (2017) Major fraction of black carbon is flushed from the melting new hampshire snowpack nearly as quickly as soluble impurities. _Journal of Geophysical Research: Atmospheres_ , 122(1), 537–553 (https://doi.org/10.1002/2016JD025351)
  * Li and others, (2008) Li C, Grobmyer SR, Massol N, Liang X, Zhang Q, Chen L, Fajardo LL and Jiang H (2008) Noninvasive in vivo tomographic optical imaging of cellular morphology in the breast: Possible convergence of microscopic pathology and macroscopic radiology. _Medical Physics_ , 35(6Part1), 2493–2501
  * Libois and others, (2013) Libois Q, Picard G, France JL, Arnaud L, Dumont M, Carmagnola CM and King MD (2013) Influence of grain shape on light penetration in snow. _The Cryosphere_ , 7(6), 1803–1818 (10.5194/tc-7-1803-2013)
  * Libois and others, (2014) Libois Q, Picard G, Dumont M, Arnaud L, Sergent C, Pougatch E, Sudul M and Vial D (2014) Experimental determination of the absorption enhancement parameter of snow. _Journal of Glaciology_ , 60(222), 714–724 (10.3189/2014JoG14J015)
  * McCullagh, (2019) McCullagh P (2019) _Generalized linear models_. Routledge
  * Nolin and Dozier, (2000) Nolin A and Dozier J (2000) A hyperspectral method for remotely sensing the grain size of snow. _Remote sensing of Environment_ , 74(2), 207–216
  * Painter and others, (2010) Painter TH, Deems JS, Belnap J, Hamlet AF, Landry CC and Udall B (2010) Response of colorado river runoff to dust radiative forcing in snow. _Proceedings of the National Academy of Sciences_ , 107(40), 17125–17130 (10.1073/pnas.0913139107)
  * Painter and others, (2012) Painter TH, Bryant AC and Skiles SM (2012) Radiative forcing by light absorbing impurities in snow from modis surface reflectance data. _Geophysical Research Letters_ , 39(17) (https://doi.org/10.1029/2012GL052457)
  * Quarto and others, (2014) Quarto G, Spinelli L, Pifferi A, Torricelli A, Cubeddu R, Abbate F, Balestreri N, Menna S, Cassano E and Taroni P (2014) Estimate of tissue composition in malignant and benign breast lesions by time-domain optical mammography. _Biomed. Opt. Express_ , 5(10), 3684–3698 (10.1364/BOE.5.003684)
  * Robledano and others, (2023) Robledano A, Picard G, Dumont M, Flìn F, Arnaud L and Libois Q (2023) Unraveling the optical shape of snow. _Nature Communications_ , 14(1), 3955 (10.1038/s41467-023-39671-3)
  * Satat, (2019) Satat G (2019) _All photons imaging : time-resolved computational imaging through scattering for vehicles and medical applications with probabilistic and data-driven algorithms_. Ph.D. thesis, Massachusetts Institute of Technology, Cambridge, MA
  * Sevick and others, (1991) Sevick E, Chance B, Leigh J, Nioka S and Maris M (1991) Quantitation of time- and frequency-resolved optical spectra for the determination of tissue oxygenation. _Analytical Biochemistry_ , 195(2), 330–351, ISSN 0003-2697 (https://doi.org/10.1016/0003-2697(91)90339-U)
  * Skiles and others, (2018) Skiles S, Flanner M, Cook J, Dumont M and Painter T (2018) Radiative forcing by light-absorbing particles in snow. _Nature Clim. Change_ , 8(11), 964–971
  * Smith and others, (2018) Smith BE, Gardner A, Schneider A and Flanner M (2018) Modeling biases in laser-altimetry measurements caused by scattering of green light in snow. _Remote Sensing of Environment_ , 215, 398–410 (https://doi.org/10.1016/j.rse.2018.06.012)
  * Sumlin and others, (2018) Sumlin BJ, Heinson WR and Chakrabarty RK (2018) Retrieving the aerosol complex refractive index using pymiescatt: A mie computational package with visualization capabilities. _Journal of Quantitative Spectroscopy and Radiative Transfer_ , 205, 127–134 (https://doi.org/10.1016/j.jqsrt.2017.10.012)
  * Várnai and Cahalan, (2007) Várnai T and Cahalan RF (2007) Potential for airborne offbeam lidar measurements of snow and sea ice thickness. _Journal of Geophysical Research: Oceans_ , 112(C12) (https://doi.org/10.1029/2007JC004091)
  * Warren, (2019) Warren SG (2019) Optical properties of ice and snow. _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 377(2146), 20180161 (10.1098/rsta.2018.0161)
  * Warren and Brandt, (2008) Warren SG and Brandt RE (2008) Optical constants of ice from the ultraviolet to the microwave: A revised compilation. _Journal of Geophysical Research_ , 113(D14)
  * Welch and van Gemert, (1995) Welch AJ and van Gemert MJ (1995) _Optical-thermal response of laser-irradiated tissue_. Springer, New York
  * Wiscombe and Warren, (1980a) Wiscombe WJ and Warren SG (1980a) A model for the spectral albedo of snow. i: Pure snow. _Journal of Atmospheric Sciences_ , 37, 2712–2733
  * Wiscombe and Warren, (1980b) Wiscombe WJ and Warren SG (1980b) A model for the spectral albedo of snow. ii: Snow containing atmospheric aerosols. _Journal of Atmospheric Sciences_ , 37(12), 2735–2745
  * Zappa and others, (2007) Zappa F, Tisa S, Tosi A and Cova S (2007) Principles and features of single-photon avalanche diode arrays. _Sensors and Actuators A: Physical_ , 140(1), 103–112, ISSN 0924-4247 (https://doi.org/10.1016/j.sna.2007.06.021)
  * Zege and others, (2011) Zege E, Katsev I, Malinka A, Prikhach A, Heygster G and Wiebe H (2011) Algorithm for retrieval of the effective snow grain size and pollution amount from satellite measurements. _Remote Sensing of Environment_ , 115(10), 2674–2685, ISSN 0034-4257 (https://doi.org/10.1016/j.rse.2011.06.001)